added function to resample a fiber in the network, adding resampled edges to the network

pranathivemuri
1 parent 1c3c8f1a
Showing 1 changed file with 132 additions and 88 deletions Show diff stats
stim/visualization/network.h
@@ -14,7 +14,6 @@
 #include <ANN/ANN.h>
 #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
@@ -171,10 +170,13 @@ class network{
 			edges_to_str()
 			{
 				std::ostringstream ss;
+				std::cout<<"here"<<std::endl;
 //				ss << "[";
 				for(int i = 0; i < edges.size()-1; i++)
 				{
+					std::cout<<"here"<<i<<std::endl;
 					ss << edges[i] << ";";
+					std::cout<<edges.size()-1<<std::endl;
 				}
 				ss << edges[edges.size()-1];
 //				ss << "]";
@@ -187,6 +189,8 @@ class network{
 	std::vector<edge*> E;       //list of pointers to edges.
 	std::vector<node> V;	    //list of nodes.
+	std::vector< stim::vec<T> > allVerticesList; // all nodes before sampling
+	std::vector<stim::vec<T> > allVerticesAfterSampling ; //all vertices after sampling
 	///Returns the number of edges in the network.
 	unsigned int
@@ -201,7 +205,11 @@ class network{
 	{
 		return V.size();
 	}
-
+	unsigned int
+	sizeAfterSamplingV()
+	{
+		return allVerticesAfterSampling.size();
+	}
 /*	//adds an edge from two std::vectors with positions and radii.
 	void
 	addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int endId)
@@ -500,7 +508,7 @@ class network{
 	totalEdges()
 	{
 		int totalEdgeVertices=0;int N=0;
-		for (int i=0; i < sizeE(); i ++)
+		for (unsigned 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;
@@ -512,46 +520,13 @@ class network{
 	lengthOfNetwork()
 	{
 		float networkLength=0;float N=0;
-		for (int i=0; i < sizeE(); i ++)
+		for (unsigned int i=0; i < sizeE(); i ++)
 		{
 			N = E[i]->length();
 			networkLength = networkLength + N;
 		}
 		return networkLength;
 	}
-	/// 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
-	strObj(std::string* strArray, int sizeV)
-	{
-		std::stringstream ss;
-		std::stringstream oss;
-		for(unsigned int i = 0; i < N; i++)
-		{
-			ss.str(std::string());
-			for(unsigned int d = 0; d < 3; d++)
-			{
-				ss<<c[i][d];
-			}
-			oss<<getIndexes(strArray, ss.str(), sizeV)<<" ";
-		}
-		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> >
@@ -584,6 +559,7 @@ class network{
 	objToNetwork(stim::obj<T> objInput)
 	{
 		stim::network<T> nwc;
+		//network network2;
 		// 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
@@ -595,7 +571,8 @@ class network{
 		// 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);
+		std::vector< stim::vec<T> > vertices = objInput.getAllV(vertexIndices);
+		nwc.addVertices(vertices);
 		for (unsigned i =1; i< objInput.numL() + 1; i++) 
 			{
 			    // edges/fibers could be added to the network class
@@ -630,28 +607,33 @@ class network{
 			   {
 				   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);
+			    //edge* a = new edge(fiberPositions1,radii); 
+				//std::cout<<"here"<<std::endl;
+			    //E.push_back(a);
+			    std::vector<stim::vec<T> > newPointList;
+			    newPointList  = Resample(fiberPositions1);
+				int numPointsOnnewFiber = newPointList.size();
+				nwc.addVerticesAfterSamplimg(newPointList);
+				std::vector<T> radii(numPointsOnnewFiber); // allocating radii to zero
+				nwc.addEdge(newPointList, 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];
+				int endId = numPointsOnnewFiber -1;
+				v0[0] = newPointList[0][0];v0[1] = newPointList[0][1];v0[2] = newPointList[0][2];
+				vN[0] = newPointList[endId][0];vN[1] = newPointList[endId][1];vN[2] = newPointList[endId][2];
 				// VISITED INDEXES OF the nodes are set to true
-				if(!validList[objInput.getIndex(allVerticesList, v0)])
+				if(!validList[objInput.getIndex(vertices, v0)])
 				{
 					//V.push_back(node(v0));
 					nwc.addNode(v0);
-					validList[objInput.getIndex(allVerticesList, v0)] = true;
+					validList[objInput.getIndex(vertices, v0)] = true;
 				}
-				if(!validList[objInput.getIndex(allVerticesList, vN)])
+				if(!validList[objInput.getIndex(vertices, vN)])
 				{
 					//V.push_back(node(vN));
 					nwc.addNode(vN);
-					validList[objInput.getIndex(allVerticesList, vN)] = true;
+					validList[objInput.getIndex(vertices, vN)] = true;
 				}
 			}
 		return nwc;
@@ -672,6 +654,7 @@ class network{
 		kdTree2 = generateKdTreeFromNetwork(network2);
 		std::cout<<"Test Network:network and kdtree generated"<<std::endl;
 		return boost::make_tuple(kdTree1, kdTree2, network1, network2);
+		std::cout<<"out of this loop"<<std::endl;
 	}
 	//distance between two points
 	double dist(std::vector<T> p0, std::vector<T> p1)
@@ -686,7 +669,7 @@ class network{
 	double sum(stim::vec<T> metricList)
 	{
 		float sumMetricList = 0;
-		for (int count=0; count<metricList.size(); count++)
+		for (unsigned int count=0; count<metricList.size(); count++)
 		{ sumMetricList += metricList[count];}
 		return sumMetricList;
 	}
@@ -694,21 +677,22 @@ class network{
 	ANNkd_tree* 
 	generateKdTreeFromNetwork(stim::network<T> network)			//kd-tree stores all points in the fiber for fast searching
 	{       
+		std::cout<<"kd trees generated"<<std::endl;
 		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
+		unsigned int n_data = network.sizeAfterSamplingV(); //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++)
+		//stim::vec<T> node;
+		for(unsigned int i = 0; i < n_data; i++)
 			{
-				node = network.V[i].getPosition();
+				//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
+					c[i][d] = double(network.allVerticesAfterSampling[i][d]);				//copy the coordinate
 				}
 			r[i] = r[i];						//copy the radius
 			}
@@ -731,54 +715,38 @@ class network{
 	{
 		V.push_back(nodes);
 	}
-	///exports the graph.
 	void
-	to_csv()
+	addVertices(std::vector< stim::vec<T> > vertices)
 	{
-		std::ofstream ofs;
-		ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app);
-		for(int i = 0; i < V.size(); i++)
+		for (unsigned int i=0; i < vertices.size(); i ++)
 		{
-			ofs << V[i].edges_to_str() << "\n";
+			allVerticesList.push_back(vertices[i]);
 		}
-		ofs.close();
 	}
-
-	///exports the graph.
 	void
-	to_gdf()
+	addVerticesAfterSamplimg(std::vector< stim::vec<T> > vertices)
 	{
-		std::ofstream ofs;
-		ofs.open("Graph.gdf", std::ofstream::out | std::ofstream::app);
-		ofs << "nodedef>name VARCHAR\n";
-		for(int i = 0; i < V.size(); i++)
+		for (unsigned int i=0; i < vertices.size(); i ++)
 		{
-			ofs << i << "\n";
+			allVerticesAfterSampling.push_back(vertices[i]);
 		}
-		ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n";
-		for(int i = 0; i < E.size(); i++)
-		{
-			ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <<E[i]->n_pts()
-			 << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n";
-		} 
-		ofs.close();
 	}
 	// gaussian function
 	float gaussianFunction(float x, float std=25)
 	{
-		float normalization = multfactor * std;
-		float evaluate = 1.0 - (exp(-x/(2*std*std)));
+		float evaluate = 1.0f - ((exp(-x/(2*std*std))));
 		return evaluate;
 	}
+	// compares point on a network to a kd tree for two skeletons
 	boost::tuple< float, float > 
-	compareSkeletons(boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<int>, stim::network<int> > networkKdtree)
+	compareSkeletons(boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<float>, stim::network<float> > 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
+	// gaussian distance of points on network to Kdtree
 	float 
 	CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network<T> networkTruth)
 	{
@@ -793,14 +761,14 @@ class network{
 		float totalNetworkLength = networkTruth.lengthOfNetwork();
 		stim::vec<float> fiberMetric(networkTruth.sizeE());
 		//for each fiber
-		for (unsigned i=0; i < networkTruth.sizeE(); i ++)
+		for (unsigned int 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++)
+			for (unsigned int j = 0; j < N - 1; j++)
 			{
 				ANNpoint queryPt1; queryPt1 = annAllocPt(3);
 				ANNpoint queryPt2; queryPt2 = annAllocPt(3);
@@ -814,14 +782,56 @@ class network{
 				}
 				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]);
+				std::cout<<float(dists1[0])<<"dist1-----"<<float(dists2[0])<<"dist2---"<<std::endl;
+				float dist1 = gaussianFunction(float(dists1[0]));float dist2 = gaussianFunction(float(dists2[0]));
 				segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ;
 			}
-			fiberMetric[i] = sum(segmentMetric);
+			fiberMetric[i] = sum(segmentMetric)/;
 		}
 		netmetsMetric =  sum(fiberMetric)/totalNetworkLength ;
 		return netmetsMetric;	
 	}
+    //resample a fiber in the network
+	std::vector<stim::vec<T> > Resample(std::vector< stim::vec<T> > fiberPositions, float spacing=25.0)
+	{
+		std::vector<T> p1(3), p2(3), v(3);
+		stim::vec<T> p(3);
+		std::vector<stim::vec<T> > newPointList;
+		for(unsigned int f=0; f<fiberPositions.size()-1; f++)
+		{
+			T lengthSegment = dist(p1,p2);
+			for(int d=0; d<3;d++)
+			{
+				p1[d] = T(fiberPositions[f][d]);
+			    p2[d] = T(fiberPositions[f + 1][d]);
+			}// for each dimension
+			if( lengthSegment >= spacing )
+				{	for (int dim=0; dim<3;dim++) //find the direction of travel
+						{v[dim] = p2[dim] - p1[dim];}
+					//set the step size to the voxel size
+					T step;
+					for(step=0.0; step<lengthSegment; step+=spacing)
+					{
+						for(int i=0; i<3;i++)
+							{p[i] = p1[i] + v[i]*(step/lengthSegment);}
+						newPointList.push_back(p);
+					}
+				}
+			else 
+				newPointList = fiberPositions;
+			}
+		return newPointList;
+	}
+	// modulus of a vector
+	T
+	Length(std::vector<T> v)
+	{
+	 	T sum=0;
+        for (int i=0; i<v.size(); i++)
+			sum += v[i] * v[i];
+	    return sum;
+	}
+	// function to remove characters from a string
 	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() );
@@ -833,9 +843,9 @@ class network{
 	{
 		std::ofstream ofs;
 		ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app);
-		int num = V.size();
-		string *strArray = new string[num];
-		for(int i = 0; i < allVerticesList.size(); i++)
+		int num = allVerticesList.size();
+		std::string *strArray = new std::string[num];
+		for(unsigned int i = 0; i < allVerticesList.size(); i++)
 		{
 			std::string str;
 			str = allVerticesList[i].str();
@@ -848,10 +858,44 @@ class network{
 		{
 			 std::string str;
 			 str = E[i]->strObj(strArray, num);
+			 //removeCharsFromString(str,"0");
              ofs << "l " << str << "\n";
 		} 
 		ofs.close();
 	}	
+		///exports the graph.
+	void
+	to_csv()
+	{
+		std::ofstream ofs;
+		ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app);
+		std::cout<<"here"<<std::endl;
+		for(int i = 0; i < V.size(); i++)
+		{
+			std::cout<<"here"<<i<<std::endl;
+			ofs << V[i].edges_to_str() << "\n";
+		}
+		ofs.close();
+	}
+	///exports the graph.
+	void
+	to_gdf()
+	{
+		std::ofstream ofs;
+		ofs.open("Graph.gdf", std::ofstream::out | std::ofstream::app);
+		ofs << "nodedef>name VARCHAR\n";
+		for(int i = 0; i < V.size(); i++)
+		{
+			ofs << i << "\n";
+		}
+		ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n";
+		for(int i = 0; i < E.size(); i++)
+		{
+			ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <<E[i]->n_pts()
+			 << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n";
+		} 
+		ofs.close();
+	}
 };
 };