compare networks and give netmets metric function added in network class

pranathivemuri
1 parent 3800c27f
Showing 2 changed files with 141 additions and 24 deletions Show diff stats
stim/biomodels/fiber.h
stim/biomodels/network.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<T> q){
  
-		stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]);
+		stim::vec<double> 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<T> v(3);    //v-direction vector of the segment
 		stim::vec<T> p(3);      //- intermediate point to be added
-		std::vector<T> p1(3);   // p1 - starting point of an segment on the fiber,
-		std::vector<T> p2(3);   // p2 - ending point,
+		stim::vec<T> p1(3);   // p1 - starting point of an segment on the fiber,
+		stim::vec<T> p2(3);   // p2 - ending point,
 		double sum=0;  //distance summation
 		std::vector<stim::vec<T> > fiberPositions = centerline();
 		std::vector<stim::vec<T> > 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
  
 #ifndef STIM_NETWORK_H
 #define STIM_NETWORK_H
  
-#include <list>
 #include <stdlib.h>
+#include <assert.h> 
 #include <sstream>
 #include <fstream>
 #include <algorithm>
@@ -33,8 +33,7 @@ class network{
 		public:
 		unsigned v[2];		//unique id's designating the starting and ending
 		// default constructor
-		edge() : fiber<T>(){v[1] = -1; v[0] = -1;}
-
+		 edge() : fiber<T>(){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<unsigned> 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<T> 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<T> 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<stim::vec<T> > points_afterResampling(){
+		std::vector<stim::vec<T> > pointsVector; // points in the resampled network
+		stim::fiber<T> FIBER; // initialize a fiber used in looping through all edges casting into fibers in the network
+		std::vector<stim::vec<T> > 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<T> metricList){
+		float sumMetricList = 0;           // Initialize variable to calculate sum
+		for (unsigned int count=0; count<metricList.size(); count++) // for each element in the stim vector
+			{ sumMetricList += metricList[count];}                   // running sum of values
+		return sumMetricList;
+	}
+
+
+	// distance between two points
+	double dist(stim::vec<T> p0, stim::vec<T> p1){
+		double sum = 0;                             // initialize variables
+		stim::vec<T> 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<T> 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<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
+				}
+			}
+		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
+	    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
+		{
+			FIBER = networkTruth.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
+			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++)
+			{
+				stim::vec<T> p1 = fiberPoints[j];stim::vec<T> 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