cleaned up the stim::network comparison calculation

David Mayerich
1 parent 4df9ec0e
Showing 2 changed files with 72 additions and 28 deletions Show diff stats
stim/biomodels/fiber.h
stim/biomodels/network.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<T> resample(T spacing=25.0)
+	stim::fiber<T> resample(T spacing)
 	{
  
  
@@ -473,7 +473,7 @@ public:
 		std::vector<stim::vec<T> > fiberPositions = centerline();
 		std::vector<stim::vec<T> > 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++)
 		{
  
@@ -44,7 +44,7 @@ class network{
  
 		/// Resamples an edge by calling the fiber resampling function
 		edge resample(T spacing){
-			edge e(fiber<T>::resample());	//call the fiber->edge constructor
+			edge e(fiber<T>::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<T> 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<T> networkTruth, float sigma){
+	float compare(stim::network<T> 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<stim::vec<T> > 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<stim::vec<T> > 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<float> 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<T> 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<float> fiberMetric(A.edges());    // allocate space for accumulating fiber metric as each fiber is evaluated
 	    stim::fiber<T> 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<T> > 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<float> 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