minor changes to network

Pavel Govyadinov
1 parent 6ada8448
Showing 3 changed files with 82 additions and 60 deletions Show diff stats
stim/biomodels/network.h
stim/math/vec3.h
stim/visualization/gl_aaboundingbox.h
@@ -33,7 +33,10 @@ class network{
 		public:
 		unsigned v[2];		//unique id's designating the starting and ending
 		// default constructor
-		 edge() : cylinder<T>(){v[1] = -1; v[0] = -1;}
+		edge() : cylinder<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
@@ -54,7 +57,7 @@ class network{
 		/// Output the edge information as a string
 		std::string str(){
 			std::stringstream ss;
-			ss<<"("<<cylinder<T>::size()<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
+			ss<<"("<<cylinder<T>::size()<<")\tl = "<<this.length()<<"\t"<<v[0]<<"----"<<v[1];
 			return ss.str();
 		}
  
@@ -64,9 +67,9 @@ class network{
 	class vertex : public stim::vec3<T>
 	{
 		public:
-			//std::vector<unsigned int> edges;		//indices of edges connected to this node.
-			std::vector<unsigned int> e[2];			//indices of edges going out (e[0]) and coming in (e[1])
-			//stim::vec3<T> p;						//position of this node in physical space.
+			//std::vector<unsigned int> edges;					//indices of edges connected to this node.
+			std::vector<unsigned int> e[2];						//indices of edges going out (e[0]) and coming in (e[1])
+			//stim::vec3<T> p;							//position of this node in physical space.
  
 			//constructor takes a stim::vec
 			vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
@@ -99,6 +102,18 @@ protected:
  
 public:
  
+	///default constructor
+	network()
+	{
+		
+	}
+
+	///constructor with a file to load.
+	network(std::string fileLocation)
+	{
+		load_obj(fileLocation);
+	}
+
 	///Returns the number of edges in the network.
 	unsigned int edges(){
 		return E.size();
@@ -110,71 +125,74 @@ public:
 	}
  
 	stim::cylinder<T> get_cylinder(unsigned f){
-		return E[f];					//return the specified edge (casting it to a fiber)
+		return E[f];									//return the specified edge (casting it to a fiber)
 	}
  
 	//load a network from an OBJ file
 	void load_obj(std::string filename){
  
-		stim::obj<T> O;			//create an OBJ object
-		O.load(filename);		//load the OBJ file as an object
+		stim::obj<T> O;									//create an OBJ object
+		O.load(filename);								//load the OBJ file as an object
  
-		std::vector<unsigned> id2vert;	//this list stores the OBJ vertex ID associated with each network vertex
+		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++){
  
-			std::vector< stim::vec<T> > c;				//allocate an array of points for the vessel centerline
+			std::vector< stim::vec<T> > c;						//allocate an array of points for the vessel centerline
 			O.getLine(l, c);							//get the fiber centerline
  
 			std::vector< stim::vec3<T> > c3(c.size());
 			for(size_t j = 0; j < c.size(); j++)
 				c3[j] = c[j];
  
-			edge new_edge = c3;							//create an edge from the given centerline
-			unsigned int I = new_edge.size();			//calculate the number of points on the centerline
+	//		edge new_edge = c3;		///This is dangerous.
+			edge new_edge(c3);
+					
+			//create an edge from the given centerline
+			unsigned int I = new_edge.size();					//calculate the number of points on the centerline
  
 			//get the first and last vertex IDs for the line
-			std::vector< unsigned > id;					//create an array to store the centerline point IDs
+			std::vector< unsigned > id;						//create an array to store the centerline point IDs
 			O.getLinei(l, id);							//get the list of point IDs for the line
 			i[0] = id.front();							//get the OBJ ID for the first element of the line
 			i[1] = id.back();							//get the OBJ ID for the last element of the line
  
-			std::vector<unsigned>::iterator it;			//create an iterator for searching the id2vert array
+			std::vector<unsigned>::iterator it;					//create an iterator for searching the id2vert array
 			unsigned it_idx;							//create an integer for the id2vert entry index
  
 			//find out if the nodes for this fiber have already been created
-			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
+			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
+			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
-				id2vert.push_back(i[0]);						//add the ID to the ID->vertex conversion list
+				new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
+				new_edge.v[0] = V.size();					//add the new edge 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
+			else{									//if the vertex already exists
 				V[it_idx].e[0].push_back(E.size());				//add the current edge as outgoing
 				new_edge.v[0] = it_idx;
 			}
  
-			it = find(id2vert.begin(), id2vert.end(), i[1]);	//look for the second ID
+			it = find(id2vert.begin(), id2vert.end(), i[1]);			//look for the second ID
 			it_idx = std::distance(id2vert.begin(), it);
-			if(it == id2vert.end()){							//if i[1] hasn't already been used
-				vertex new_vertex = new_edge[I-1];							//create a new vertex, assign it a position
-				new_vertex.e[1].push_back(E.size());						//add the current edge as incoming
+			if(it == id2vert.end()){						//if i[1] hasn't already been used
+				vertex new_vertex = new_edge[I-1];				//create a new vertex, assign it a position
+				new_vertex.e[1].push_back(E.size());				//add the current edge as incoming
 				new_edge.v[1] = V.size();
-				V.push_back(new_vertex);									//add the new vertex to the vertex list
-				id2vert.push_back(i[1]);						//add the ID to the ID->vertex conversion list
+				V.push_back(new_vertex);					//add the new vertex to the vertex list
+				id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
 			}
-			else{												//if the vertex already exists
+			else{									//if the vertex already exists
 				V[it_idx].e[1].push_back(E.size());				//add the current edge as incoming
 				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
  
 		}
 	}
@@ -199,17 +217,17 @@ public:
 	/// This function resamples all fibers in a network given a desired minimum spacing
 	/// @param spacing is the minimum distance between two points on the network
 	stim::network<T> resample(T spacing){
-		stim::network<T> n;					//create a new network that will be an exact copy, with resampled fibers
-		n.V = V;							//copy all vertices
+		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(edges());				//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 < 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
+		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
 	}
  
  
@@ -236,14 +254,14 @@ public:
 	/// @param m is the magnitude value to use. The default is 0 (usually radius).
 	T average(unsigned m = 0){
  
-		T M, L;												//allocate space for the total magnitude and length
-		M = L = 0;											//initialize both the initial magnitude and length to zero
-		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+		T M, L;										//allocate space for the total magnitude and length
+		M = L = 0;									//initialize both the initial magnitude and length to zero
+		for(unsigned e = 0; e < E.size(); e++){						//for each edge in the network
 			M += E[e].integrate(m);							//get the integrated magnitude
-			L += E[e].length();								//get the edge length
+			L += E[e].length();							//get the edge length
 		}
  
-		return M / L;					//return the average magnitude
+		return M / L;									//return the average magnitude
 	}
  
 	/// This function compares two networks and returns the percentage of the current network that is missing from A.
@@ -256,17 +274,17 @@ public:
 		R = (*this);									//initialize the result with the current network
  
 		//generate a KD-tree for network A
-		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 = A.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
+		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 = A.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);
  
 		unsigned t = 0;
 		for(unsigned e = 0; e < A.E.size(); e++){					//for each edge in the network
-			for(unsigned p = 0; p < A.E[e].size(); p++){			//for each point in the edge
+			for(unsigned p = 0; p < A.E[e].size(); p++){				//for each point in the edge
 				for(unsigned d = 0; d < 3; d++){				//for each coordinate
  
 					c[t][d] = A.E[e][p][d];
@@ -276,27 +294,27 @@ public:
 		}
  
 		//compare each point in the current network to the field produced by A
-		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
+		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 dists = new ANNdist[1];     // near neighbor distances
-		ANNidxArray nnIdx = new ANNidx[1];				// near neighbor indices // allocate near neigh indices
+		ANNdistArray dists = new ANNdist[1];     					// near neighbor distances
+		ANNidxArray nnIdx = new ANNidx[1];						// near neighbor indices // allocate near neigh indices
  
 		stim::vec3<T> p0, p1;
 		float m1;
-		float M = 0;											//stores the total metric value
-		float L = 0;											//stores the total network length
+		float M = 0;									//stores the total metric value
+		float L = 0;									//stores the total network length
 		ANNpoint queryPt = annAllocPt(3);
 		for(unsigned e = 0; e < R.E.size(); e++){					//for each edge in A
-			R.E[e].add_mag(0);									//add a new magnitude for the metric
+			R.E[e].add_mag(0);							//add a new magnitude for the metric
  
-			for(unsigned p = 0; p < R.E[e].size(); p++){			//for each point in the edge
+			for(unsigned p = 0; p < R.E[e].size(); p++){				//for each point in the edge
  
-				p1 = R.E[e][p];									//get the next point in the edge
+				p1 = R.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
+				kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);		//find the distance between A and the current network
 				m1 = 1.0f - gaussianFunction((float)dists[0], sigma);		//calculate the metric value based on the distance
-				R.E[e].set_mag(m1, p, 1);						//set the error for the second point in the segment
+				R.E[e].set_mag(m1, p, 1);					//set the error for the second point in the segment
  
 			}
 		}
@@ -3,6 +3,7 @@
  
  
 #include <stim/cuda/cudatools/callable.h>
+#include <cmath>
  
  
 namespace stim{
@@ -11,6 +11,9 @@ class gl_aaboundingbox : public aaboundingbox&lt;T&gt;{
  
 public:
  
+	using stim::aaboundingbox<T>::A;
+	using stim::aaboundingbox<T>::B;
+
 	//default constructor
 	gl_aaboundingbox() : stim::aaboundingbox<T>(){}
  
@@ -57,4 +60,4 @@ public:
  
 };		//end namespace stim
  
-#endif
 \ No newline at end of file
+#endif