diff --git a/stim/biomodels/fiber.h b/stim/biomodels/fiber.h index 731467b..b1361ac 100644 --- a/stim/biomodels/fiber.h +++ b/stim/biomodels/fiber.h @@ -193,7 +193,6 @@ public: if(this == &rhs) return *this; //test for and handle self-assignment copy(rhs); - } /// Calculate the length of the fiber and return it. @@ -331,7 +330,7 @@ public: /// @param q is the query point used to locate the nearest point on the fiber centerline unsigned int nearest_idx(stim::vec q){ - stim::vec temp( (double) q[0], (double) q[1], (double) q[2]); + stim::vec temp((double) q[0], (double) q[1], (double) q[2]); unsigned int idx = ann(temp); //determine the index of the nearest neighbor @@ -449,6 +448,7 @@ public: return N; } + /// Bracket operator returns the element at index i /// @param i is the index of the element to be returned as a stim::vec @@ -467,8 +467,8 @@ public: std::vector v(3); //v-direction vector of the segment stim::vec p(3); //- intermediate point to be added - std::vector p1(3); // p1 - starting point of an segment on the fiber, - std::vector p2(3); // p2 - ending point, + stim::vec p1(3); // p1 - starting point of an segment on the fiber, + stim::vec p2(3); // p2 - ending point, double sum=0; //distance summation std::vector > fiberPositions = centerline(); std::vector > newPointList; // initialize list of new resampled points on the fiber @@ -476,14 +476,10 @@ public: unsigned int N = fiberPositions.size(); // number of points on the fiber for(unsigned int f=0; f< N-1; f++) { - for(unsigned int d = 0; d < 3; d++) - { - p1[d] = fiberPositions[f][d]; // starting point of the segment on the centerline - p2[d] = fiberPositions[f + 1][d]; // ending point of the segment on the centerline - v[d] = p2[d] - p1[d]; //direction vector - sum +=v[d] * v[d]; //length of segment-distance between starting and ending point - } - //newPointList.push_back(fiberPositions[f]); //always push the starting point + + p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1; + for(unsigned int d = 0; d < 3; d++){ + sum +=v[d] * v[d];} //length of segment-distance between starting and ending point T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index ec68567..7f08e40 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -1,8 +1,8 @@ #ifndef STIM_NETWORK_H #define STIM_NETWORK_H -#include #include +#include #include #include #include @@ -33,8 +33,7 @@ class network{ public: unsigned v[2]; //unique id's designating the starting and ending // default constructor - edge() : fiber(){v[1] = -1; v[0] = -1;} - + edge() : fiber(){v[1] = -1; v[0] = -1;} /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor ///@param p is an array of positions in space @@ -122,7 +121,7 @@ class network{ std::vector id2vert; //this list stores the OBJ vertex ID associated with each network vertex - unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line //for each line in the OBJ object for(unsigned int l = 1; l <= O.numL(); l++){ @@ -145,10 +144,10 @@ class network{ it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node it_idx = std::distance(id2vert.begin(), it); if(it == id2vert.end()){ //if i[0] hasn't already been used - vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position - new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing - new_edge.v[0] = V.size(); //add the new vertex to the edge - V.push_back(new_vertex); //add the new vertex to the vertex list + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = V.size(); //add the new vertex to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list } else{ //if the vertex already exists @@ -170,7 +169,7 @@ class network{ new_edge.v[1] = it_idx; } - E.push_back(new_edge); //push the edge to the list + E.push_back(new_edge); //push the edge to the list } } @@ -197,14 +196,136 @@ class network{ stim::network n; //create a new network that will be an exact copy, with resampled fibers n.V = V; //copy all vertices - n.E.resize(E.size()); //allocate space for the edge list + n.E.resize(edges()); //allocate space for the edge list //copy all fibers, resampling them in the process - for(unsigned e = 0; e < E.size(); e++){ //for each edge in the edge list + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network } - return n; //return the resampled network + return n; //return the resampled network + } + + + // this function gives sum of lengths of all the fibers in the network + float lengthOfNetwork(){ + stim::fiber FIBER; // initialize a fiber used in looping through all edges casting into fibers in the network + float networkLength=0;float N=0; // initialize variables for finding total length of all the fibers + // for each edge in the network + for (unsigned int i=0; i < E.size(); i ++) + { + FIBER = get_fiber(i); // cast each edge to fiber + N= FIBER.length(); // find length of the fiber + networkLength = networkLength + N; // running sum of fiber lengths + } + return networkLength; + } + + + // list of all the points after resampling -- function used only on a resampled network + std::vector > points_afterResampling(){ + std::vector > pointsVector; // points in the resampled network + stim::fiber FIBER; // initialize a fiber used in looping through all edges casting into fibers in the network + std::vector > pointsFiber; // resampled points on each fiber of the network + for(unsigned e = 0; e < E.size(); e++){ // for each edge in the edge list + FIBER = get_fiber(e); // Cast edge to a fiber + pointsFiber = FIBER.centerline(); // find points on the edge returns a stim vec + for (unsigned v = 0; v < FIBER.n_pts(); v++){ // iterate one point at a time from the stim::vec + pointsVector.push_back(pointsFiber[v]);} //add the points on centerline to the stim::vec points list + } + return pointsVector; + } + + + // total number of points on all fibers after resampling -- function used only on a resampled network + unsigned int total_points(){ + unsigned int n = points_afterResampling().size(); + return n; + } + + // gaussian function + float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25 + + // sum of values in a stim vector + float sum(stim::vec metricList){ + float sumMetricList = 0; // Initialize variable to calculate sum + for (unsigned int count=0; count p0, stim::vec p1){ + double sum = 0; // initialize variables + stim::vec v = p1 - p0;; // direction vector + for(unsigned int d = 0; d < 3; d++){ // for each dimension + sum += v[d] * v[d];} // running sum of modulus of direction vector + return sqrt(sum); + } + + + /// This function compares two networks and returns a metric + float compare(stim::network networkTruth, float sigma){ + float metric = 0.0; // initialize metric to be returned after comparing the networks + ANNkd_tree* kdt; // initialize a pointer to a kd tree + double **c; // centerline (array of double pointers) - points on kdtree must be double + unsigned int n_data = total_points(); // set the number of points + c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer + for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions + {c[i] = (double*) malloc(sizeof(double) * 3);} + std::vector > pointsVector = points_afterResampling(); //get vertices in the network after resampling + for(unsigned int i = 0; i < n_data; i++) // loop through all the vertices after resampling + { + for(unsigned int d = 0; d < 3; d++){ // for each dimension + c[i][d] = double(pointsVector[i][d]); // copy the coordinate and convert_to_double + } + } + ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double + kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray + double eps = 0; // error bound + ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances + ANNdistArray dists2; dists2 = new ANNdist[1]; // near neighbor distances + ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + int N; // number of points on the fiber in the second network + float totalNetworkLength = networkTruth.lengthOfNetwork(); + stim::vec fiberMetric(networkTruth.edges()); // allocate space for accumulating fiber metric as each fiber is evaluated + stim::fiber FIBER; // Initialize a fiber will be used in calculating the metric + //for each fiber in the second network, find nearest distance in the kd tree + for (unsigned int i=0; i < networkTruth.edges(); i ++) // loop through all the edges/fibers in the network + { + FIBER = networkTruth.get_fiber(i); // Get the fiber corresponding to the index i + std::vector< stim::vec > fiberPoints = FIBER.centerline(); // Get the points on the fiber + N = FIBER.n_pts(); // number of points on the fiber + stim::vec segmentMetric(N-1); // allocate space for a vec array that stores metric for each segmen in the fiber + // for each segment in the fiber + for (unsigned j = 0; j < N - 1; j++) + { + stim::vec p1 = fiberPoints[j];stim::vec p2 = fiberPoints[j+1]; // starting and ending points on the segments + ANNpoint queryPt1; queryPt1 = annAllocPt(3); // allocate memory for query points + ANNpoint queryPt2; queryPt2 = annAllocPt(3); + + //for each dimension of the points connecting segment + for (unsigned d = 0; d < 3; d++){ + queryPt1[d] = double(fiberPoints[j][d]); // starting point on segment in network whose closest distance on KD tree should be found + queryPt2[d] = double(fiberPoints[j + 1][d]); // ending point on segment in network whose closest distance on KD tree should be found + } + kdt->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // search the nearest point in KD tree to query point(point on network), find the distance + kdt->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // search the nearest point in KD tree to query point(point on network), find the distance + // find the gaussian of the distance and subtract it from 1 to calculate the error + float error1 = 1.0f - gaussianFunction(float(dists1[0]), sigma);float error2 = 1.0f - gaussianFunction(float(dists2[0]), sigma); + // average the error and scale it with distance between the points + segmentMetric[j] = ((error1 + error2) / 2) * dist(p1, p2) ; + /*segmentDists[j] = dist(p1, p2); */ + } + fiberMetric[i] = sum(segmentMetric); + /*distsList[i] = sum(segmentDists);*/ + } + /*assert (sum(distsList)==totalNetworkLength);*/ + metric = sum(fiberMetric)/totalNetworkLength; //normalize the scaled error of all the points with total length of the network + assert (0<=metric <=1); //assert metroc os always less than or equal to one and greater than or equal to zero + return metric; } }; //end stim::network class }; //end stim namespace -- libgit2 0.21.4