added all the functions to compute netmets metric ie. objtonetwork conversion,kd…

… tree generation, comparison of network and kdtree

added all the functions to compute netmets metric ie. objtonetwork conversion,kd…
… tree generation, comparison of network and kdtree
pranathivemuri
1 parent 018bda33
Showing 1 changed file with 279 additions and 32 deletions Show diff stats
stim/visualization/network.h
 #ifndef STIM_NETWORK_H
 #define STIM_NETWORK_H
  
+#include <list>
+#include <stdlib.h>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+#include <string.h>
+#include <math.h>
 #include <stim/math/vector.h>
 #include <stim/visualization/obj.h>
-#include <list>
+#include <stim/visualization/fiber.h>
 #include <ANN/ANN.h>
-#include "fiber.h"
-#include <sstream>
+#include <boost/tuple/tuple.hpp>
  
+#define multfactor 2.51
  
 namespace stim{
 /** This is the a class that interfaces with gl_spider in order to store the currently
@@ -77,10 +84,7 @@ class network{
 		setEndpoints(int id1, int id2)
 		{
 			Nodes[1] = id1; Nodes[2] = id2;
-		}
-
-	
-			
+		}		
  
 	};	
  
@@ -480,40 +484,74 @@ class network{
  
 	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
 	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
-	///adds an edge from two std::vectors with positions and radii.
-	void
-	addEdge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
+	///returns the forest as a std::string. For testing only.
+	std::string
+	edges_to_str()
 	{
-		edge *a = new edge(pos,radii);
-		E.push_back(a);
+		std::stringstream ss;
+		for(unsigned int i = 0; i < E.size(); i++)
+		{
+			ss << i << ": " << E[i]->str() << std::endl;
+		}
+		return(ss.str());
 	}
-	void
-	addNode(stim::vec<T> nodes)
+	// total number of points on all edges!=fibers in the network
+	int
+	totalEdges()
 	{
-			V.push_back(nodes);
+		int totalEdgeVertices=0;int N=0;
+		for (int i=0; i < sizeE(); i ++)
+		{// FOR N points on the fiber N-1 edges are possible
+			N = E[i]->n_pts();
+			totalEdgeVertices = totalEdgeVertices + N- 1;
+		}
+		return totalEdgeVertices;
 	}
-
-	///adds an edge from two std::lists with positions and radii.
-	/// NOT NECESSARY.
-	void
-	addEdge(std::list< stim::vec<T> > pos, std::list<T> radii)
+	// sum of all the fiber lengths
+	float
+	lengthOfNetwork()
 	{
-		edge a = edge(pos,radii);
-		E.push_back(a);
+		float networkLength=0;float N=0;
+		for (int i=0; i < sizeE(); i ++)
+		{
+			N = E[i]->length();
+			networkLength = networkLength + N;
+		}
+		return networkLength;
 	}
-
-	///returns the forest as a std::string. For testing only.
+	/// get index of a node on a fiber 
+	// by matching the node on fiber to already set vertices (both strings)
+	// used in obj file conversion
+	int 
+	getIndexes(std::string* input, std::string searched, int sizeV) 
+	{
+		int result = 0;
+		for (int i = 0; i < sizeV; i++)
+		{
+			if (input[i] == searched) 
+			{
+				result = i + 1;
+			}
+		}
+		return result;
+    }
+	// strObj returns a string of fiber indices corresponding to vectors of positions in the fiber including intermediate nodes
 	std::string
-	edges_to_str()
+	strObj(std::string* strArray, int sizeV)
 	{
 		std::stringstream ss;
-		for(unsigned int i = 0; i < E.size(); i++)
+		std::stringstream oss;
+		for(unsigned int i = 0; i < N; i++)
 		{
-			ss << i << ": " << E[i]->str() << std::endl;
+			ss.str(std::string());
+			for(unsigned int d = 0; d < 3; d++)
+			{
+				ss<<c[i][d];
+			}
+			oss<<getIndexes(strArray, ss.str(), sizeV)<<" ";
 		}
-		return(ss.str());
+		return oss.str();
 	}
-
 	///@param i: index of the required fiber.
 	///Returns an std::vector with the centerline of the ith fiber in the edgelist.
 	std::vector< stim::vec<T> >
@@ -541,8 +579,158 @@ class network{
 		ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]";
 		return ss.str();
 	}
-	
+	//load an obj file into a network class
+	stim::network<T>
+	objToNetwork(stim::obj<T> objInput)
+	{
+		stim::network<T> nwc;
+		// function to convert an obj file loaded using stim/visualization/obj.h
+		// to a 3D network class using network.h methods.
+		std::vector< stim::vec<T> > fiberPositions; //initialize positions on the fiber to be added to network
+		objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file
+		// to verify if the nodes are already pushed into node list
+		std::vector<bool> validList;
+		validList.assign(objInput.numV(), false);
+		// go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes
+		// using addNode function for adding nodes
+		// and the line as an edge(belongs to fiber class) using addEdge function
+		std::vector<unsigned> vertexIndices(objInput.numV());
+		std::vector< stim::vec<T> > allVerticesList = objInput.getAllV(vertexIndices);
+		for (unsigned i =1; i< objInput.numL() + 1; i++) 
+			{
+			    // edges/fibers could be added to the network class
+				std::vector< stim::vec<T> > fiberPositions;
+
+			    objInput.getLine(i, fiberPositions);
+			    // finding size to allocate radii
+			    int numPointsOnFiber = fiberPositions.size();
+			    // to extend it to a 3D postion if it is a 1D/2D vertex in space
+				std::vector< stim::vec<T> > fiberPositions1(numPointsOnFiber);
+			   // 2D case append and assign the last dimension to zero
+			   if (numDim == 2)
+			   {//  2D case append and assign the last dimension to zero repeating it for all the points on fiber
+				   for (int j = 0; j < numPointsOnFiber; j ++)
+				   {
+					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
+					   fiberPositions1[j][numDim -1 ] = fiberPositions[j][1];
+					   fiberPositions1[j][numDim] = 0;
+				   }
+			   }
+			   // 1D case append and assign last two dimensions to zero
+			   else if (numDim == 1)
+			   {
+				   for (int j = 0; j < numPointsOnFiber; j ++)
+				   {
+					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
+					   fiberPositions1[j][numDim -1 ] = 0;
+					   fiberPositions1[j][numDim] = 0;
+				   }
+			   }
+			   else
+			   {
+				   fiberPositions1 = fiberPositions;
+			   }
+			   std::vector<T> radii(numPointsOnFiber); // allocating radii to zero
+		   		// then add edge 
+			 //   edge* a = new edge(fiberPositions1,radii); 
+				//E.push_back(a);
+			    nwc.addEdge(fiberPositions1, radii);
+				// add nodes
+				stim::vec<T> v0(3);stim::vec<T> vN(3);
+				int endId = numPointsOnFiber -1;
+				v0[0] = fiberPositions1[0][0];v0[1] = fiberPositions1[0][1];v0[2] = fiberPositions1[0][2];
+				vN[0] = fiberPositions1[endId][0];vN[1] = fiberPositions1[endId][1];vN[2] = fiberPositions1[endId][2];
+				// VISITED INDEXES OF the nodes are set to true
+				if(!validList[objInput.getIndex(allVerticesList, v0)])
+				{
+					//V.push_back(node(v0));
+					nwc.addNode(v0);
+					validList[objInput.getIndex(allVerticesList, v0)] = true;
+				}
+				if(!validList[objInput.getIndex(allVerticesList, vN)])
+				{
+					//V.push_back(node(vN));
+					nwc.addNode(vN);
+					validList[objInput.getIndex(allVerticesList, vN)] = true;
+				}
+			}
+		return nwc;
+		}
+	// convert ground and test case to network ,kd trees
+	boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<T>, stim::network<T> > 
+	LoadNetworks(std::string gold_filename, std::string test_filename)
+	{
+		ANNkd_tree* kdTree1;ANNkd_tree* kdTree2;
+		using namespace stim;
+		network network1;network network2;
+		stim::obj<T> objFile1(gold_filename);
+		network1 = objToNetwork(objFile1);
+		kdTree1 = generateKdTreeFromNetwork(network1);
+		std::cout<<"Gold Standard:network and kdtree generated"<<std::endl;
+		stim::obj<T> objFile2(test_filename);
+		network2 = objToNetwork(objFile2);
+		kdTree2 = generateKdTreeFromNetwork(network2);
+		std::cout<<"Test Network:network and kdtree generated"<<std::endl;
+		return boost::make_tuple(kdTree1, kdTree2, network1, network2);
+	}
+	//distance between two points
+	double dist(std::vector<T> p0, std::vector<T> p1)
+	{
  
+		double sum = 0;
+		for(unsigned int d = 0; d < 3; d++)
+			sum += p0[d] * p1[d];
+		return sqrt(sum);
+	}
+	// sum of elements in a vector
+	double sum(stim::vec<T> metricList)
+	{
+		float sumMetricList = 0;
+		for (int count=0; count<metricList.size(); count++)
+		{ sumMetricList += metricList[count];}
+		return sumMetricList;
+	}
+	//generate a kd tree from network
+	ANNkd_tree* 
+	generateKdTreeFromNetwork(stim::network<T> network)			//kd-tree stores all points in the fiber for fast searching
+	{       
+		ANNkd_tree* kdt;
+		double **c;						//centerline (array of double pointers)
+		float* r;						   // array of fiber radii
+		unsigned int n_data = network.sizeV(); //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
+			{c[i] = (double*) malloc(sizeof(double) * 3);}	
+		r = (float*) malloc(sizeof(float) * n_data);			//allocate space for the radii
+		stim::vec<T> node;
+		for(int i = 0; i < n_data; i++)
+			{
+				node = network.V[i].getPosition();
+				//convert_to_double
+				for(unsigned int d = 0; d < 3; d++){		//for each dimension
+					c[i][d] = double(node[d]);				//copy the coordinate
+				}
+			r[i] = r[i];						//copy the radius
+			}
+		ANNpointArray pts = (ANNpointArray)c;           //create an array of data points
+		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
+		return kdt;
+	}
+
+	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
+	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
+	///adds an edge from two std::vectors with positions and radii.
+	void
+	addEdge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
+	{
+		edge *a = new edge(pos,radii);
+		E.push_back(a);
+	}
+	void
+	addNode(stim::vec<T> nodes)
+	{
+		V.push_back(nodes);
+	}
 	///exports the graph.
 	void
 	to_csv()
@@ -575,6 +763,65 @@ class network{
 		} 
 		ofs.close();
 	}
+	// gaussian function
+	float gaussianFunction(float x, float std=25)
+	{
+		float normalization = multfactor * std;
+		float evaluate = 1.0 - (exp(-x/(2*std*std)));
+		return evaluate;
+	}
+	boost::tuple< float, float > 
+	compareSkeletons(boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<int>, stim::network<int> > networkKdtree)
+	{
+		float gFPR, gFNR;
+		gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>());
+		gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>());
+		return boost::make_tuple(gFPR, gFNR);
+	}
+		// gaussian distance of points on network to Kdtree
+	float 
+	CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network<T> networkTruth)
+	{
+		std::vector< stim::vec<T> > fiberPoints;
+		float netmetsMetric=0;
+		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; int numQueryPoints = networkTruth.totalEdges();
+		float totalNetworkLength = networkTruth.lengthOfNetwork();
+		stim::vec<float> fiberMetric(networkTruth.sizeE());
+		//for each fiber
+		for (unsigned i=0; i < networkTruth.sizeE(); i ++)
+		{
+			std::vector<T> p1(3); std::vector<T> p2(3);//temporary variables to store point positions
+			fiberPoints = networkTruth.E[i]->centerline();
+			N = networkTruth.E[i]->n_pts();
+			stim::vec<float> segmentMetric(N-1);
+			// for each segment in the fiber
+			for (unsigned j = 0; j < N - 1; j++)
+			{
+				ANNpoint queryPt1; queryPt1 = annAllocPt(3);
+				ANNpoint queryPt2; queryPt2 = annAllocPt(3);
+				//for each dimension of the point
+				for(unsigned int d = 0; d < 3; d++)
+				{ 
+					queryPt1[d] = double(fiberPoints[j][d]);
+					p1[d] = double(fiberPoints[j][d]);
+					queryPt2[d] = double(fiberPoints[j + 1][d]);
+					p2[d] = double(fiberPoints[j + 1][d]);// for each dimension
+				}
+				kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound
+				kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound
+				float dist1 = gaussianFunction(dists1[0]);float dist2 = gaussianFunction(dists2[0]);
+				segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ;
+			}
+			fiberMetric[i] = sum(segmentMetric);
+		}
+		netmetsMetric =  sum(fiberMetric)/totalNetworkLength ;
+		return netmetsMetric;	
+	}
 	void removeCharsFromString(std::string &str, char* charsToRemove ) {
 	   for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
 		  str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
@@ -588,10 +835,10 @@ class network{
 		ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app);
 		int num = V.size();
 		string *strArray = new string[num];
-		for(int i = 0; i < V.size(); i++)
+		for(int i = 0; i < allVerticesList.size(); i++)
 		{
 			std::string str;
-			str = V[i].str();
+			str = allVerticesList[i].str();
 			removeCharsFromString(str, "[],");
 			ofs << "v " << str << "\n";
 			removeCharsFromString(str," ");