diff --git a/stim/biomodels/fiber.h b/stim/biomodels/fiber.h index b1361ac..8e92183 100644 --- a/stim/biomodels/fiber.h +++ b/stim/biomodels/fiber.h @@ -444,7 +444,7 @@ public: return ss.str(); } /// Returns the number of centerline points in the fiber - unsigned int n_pts(){ + unsigned int size(){ return N; } @@ -461,7 +461,7 @@ public: return get_vec(N-1); } ////resample a fiber in the network - stim::fiber resample(T spacing=25.0) + stim::fiber resample(T spacing) { @@ -473,7 +473,7 @@ public: std::vector > fiberPositions = centerline(); std::vector > newPointList; // initialize list of new resampled points on the fiber // for each point on the centerline (skip if it is the last point on centerline) - unsigned int N = fiberPositions.size(); // number of points on the fiber + //unsigned int N = fiberPositions.size(); // number of points on the fiber for(unsigned int f=0; f< N-1; f++) { diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 7f08e40..f41d898 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -44,7 +44,7 @@ class network{ /// Resamples an edge by calling the fiber resampling function edge resample(T spacing){ - edge e(fiber::resample()); //call the fiber->edge constructor + edge e(fiber::resample(spacing)); //call the fiber->edge constructor e.v[0] = v[0]; //copy the vertex data e.v[1] = v[1]; @@ -238,10 +238,13 @@ class network{ // 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(); + unsigned total_points(){ + unsigned n = 0; + for(unsigned e = 0; e < E.size(); e++) + n += E[e].size(); + //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 @@ -264,40 +267,84 @@ class network{ return sqrt(sum); } + void stim2ann(ANNpoint &a, stim::vec b){ + a[0] = b[0]; + a[1] = b[1]; + a[2] = b[2]; + } + /// This function compares two networks and returns a metric - float compare(stim::network networkTruth, float sigma){ + float compare(stim::network A, 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 + c[i] = (double*) malloc(sizeof(double) * 3); + //std::vector > pointsVector = points_afterResampling(); //get vertices in the network after resampling + unsigned t = 0; + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + for(unsigned p = 0; p < E[e].size(); p++){ //for each point in the edge + for(unsigned d = 0; d < 3; d++){ //for each coordinate + + c[t][d] = E[e][p][d]; } + t++; } + } + 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 + ANNdistArray dists = new ANNdist[1]; // near neighbor distances + ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + + stim::vec p0, p1; + float m0, m1; + float M = 0; //stores the total metric value + float l; //stores the segment length + float L = 0; //stores the total network length + ANNpoint queryPt = annAllocPt(3); + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A + M = 0; //total metric value for the fiber + p0 = A.E[e][0]; //get the first point in the edge + stim2ann(queryPt, p0); + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network + m0 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance + + for(unsigned p = 1; p < A.E[e].size(); p++){ //for each point in the edge + + p1 = A.E[e][p]; //get the next point in the edge + stim2ann(queryPt, p1); + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network + m1 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance + + //average the metric and scale + l = (p1 - p0).len(); //calculate the length of the segment + L += l; //add the length of this segment to the total network length + M += (m0 + m1)/2 * l; //scale the metric based on segment length + + //swap points + p0 = p1; + m0 = m1; + } + } + + return M/L; + + + /*int N; // number of points on the fiber in the second network + float totalNetworkLength = A.lengthOfNetwork(); + stim::vec fiberMetric(A.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 + for (unsigned int i=0; i < A.edges(); i ++) // loop through all the edges/fibers in the network { - FIBER = networkTruth.get_fiber(i); // Get the fiber corresponding to the index i + FIBER = A.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 + N = FIBER.size(); // 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++) @@ -317,15 +364,12 @@ class network{ 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; + return metric;*/ } }; //end stim::network class }; //end stim namespace -- libgit2 0.21.4