added the network class, made relevant changes to the obj class to support it

David Mayerich
1 parent ac788020
Showing 2 changed files with 375 additions and 0 deletions Show diff stats
stim/biomodels/network.h
stim/visualization/obj.h
+#ifndef STIM_NETWORK_H
+#define STIM_NETWORK_H
+
+#include <stim/math/mathvec.h>
+#include <list>
+#include <ANN/ANN.h>
+
+namespace stim{
+
+/** This class provides an interface for dealing with biological networks.
+ *  It takes the following aspects into account:
+ *  	1) Network geometry and centerlines
+ *  	2) Network connectivity (a graph structure can be extracted)
+ *  	3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
+ */
+
+template<typename T>
+class network{
+
+	//helper classes
+	/// Stores information about a geometric point on the network centerline (including point position and radius)
+	//template<typename T>
+	class point : public stim::vec<T>{
+
+	public:
+		T r;
+
+		point() : stim::vec<T>(){}
+
+		//casting constructor
+		point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
+	};
+
+	//template<typename T>
+	class t_node;
+	class fiber;
+
+	//create typedefs for the iterators to simplify the network code
+	typedef typename std::list< fiber >::iterator fiber_i;
+	typedef typename std::list< t_node >::iterator t_node_i;
+
+	/// Stores information about a single capillary (a length of vessel between two branch or end points)
+	//template<typename T>
+	class fiber{
+
+	public:
+		std::list< point > P;		//geometric point positions
+
+		typename std::list< t_node >::iterator n[2];				//indices to terminal nodes
+		unsigned int id;
+
+	public:
+
+		/// Calculate the length of the fiber and return it.
+		T length(){
+
+			point p0, p1;
+			T l = 0;				//initialize the length to zero
+
+			//for each point
+			typename std::list< point >::iterator i;	//create a point iterator
+			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+
+				if(i == P.begin())						//if this is the first point, just store it
+					p1 = *i;
+				else{									//if this is any other point
+					p0 = p1;							//shift p1->p0
+					p1 = *i;							//set p1 to the new point
+					l += (p1 - p0).len();				//add the length of p1 - p0 to the running sum
+				}
+			}
+
+			return l;									//return the length
+		}
+
+		T radius(T& length){
+
+			point p0, p1;				//temporary variables to store point positions
+			T r0, r1;					//temporary variables to store radii at points
+			T l, r;						//temporary variable to store the length and average radius of a fiber segment
+			T length_sum = 0;			//initialize the length to zero
+			T radius_sum = 0;			//initialize the radius sum to zero
+
+			//for each point
+			typename std::list< point >::iterator i;	//create a point iterator
+			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+
+				if(i == P.begin()){						//if this is the first point, just store it
+					p1 = *i;
+					r1 = i->r;
+				}
+				else{									//if this is any other point
+					p0 = p1;							//shift p1->p0 and r1->r0
+					r0 = r1;
+					p1 = *i;							//set p1 to the new point
+					r1 = i->r;								//and r1
+
+					l = (p1 - p0).len();				//calculate the length of the p0-p1 segment
+					r = (r0 + r1) / 2;					//calculate the average radius of the segment
+
+					radius_sum += r * l;				//add the radius scaled by the length to a running sum
+					length_sum += l;					//add the length of p1 - p0 to the running sum
+				}
+			}
+
+			length = length_sum;						//store the total length
+			return radius_sum / length;					//return the average radius of the fiber
+		}
+
+		std::string str(){
+			std::stringstream ss;
+
+			//create an iterator for the point list
+			typename std::list<point>::iterator i;
+			for(i = P.begin(); i != P.end(); i++){
+				ss<<i->str()<<"  r = "<<i->r<<std::endl;
+			}
+
+			return ss.str();
+		}
+	};
+
+	/// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.
+	//template<typename T>
+	class t_node{
+
+	public:
+
+		unsigned int id;
+
+		//lists of edge indices for capillaries
+			//the "in" and "out" just indicate how the geometry is defined:
+			//		edges in the "in" list are geometrically oriented such that the terminal node is last
+			//		edges in the "out" list are geometrically oriented such that the terminal node is first
+		std::list< fiber_i > in;			//edge indices for incoming capillaries
+		std::list< fiber_i > out;		//edge indices for outgoing capillaries
+
+		std::string str(){
+
+			std::stringstream ss;
+
+			ss<<id<<": ";						//output the node ID
+
+			//output the IDs for both lists
+			typename std::list< fiber_i >::iterator f;
+
+			for(f = in.begin(); f != in.end(); f++){
+
+				if(f != in.begin())
+					ss<<", ";
+				ss<<(*f)->n[0]->id;
+			}
+
+			//if there are nodes in both lists, separate them by a comma
+			if(out.size() > 0 && in.size() > 0)
+				ss<<", ";
+
+			for(f = out.begin(); f != out.end(); f++){
+
+				if(f != out.begin())
+					ss<<", ";
+				ss<<(*f)->n[1]->id;
+			}
+
+
+			return ss.str();
+
+
+
+		}
+	};
+
+
+
+
+protected:
+
+	//list of terminal nodes
+	std::list<t_node> N;
+
+	//list of fibers
+	std::list<fiber> F;
+
+	/// Sets a unique ID for each terminal node and fiber
+	void set_names(){
+
+		unsigned int i;
+
+		i = 0;
+		for(t_node_i ti = N.begin(); ti != N.end(); ti++)
+			ti->id = i++;
+
+		i = 0;
+		for(fiber_i fi = F.begin(); fi != F.end(); fi++)
+			fi->id = i++;
+	}
+
+public:
+
+	std::string str(){
+
+
+		//assign names to elements of the network
+		set_names();
+
+		//create a stringstream for output
+		std::stringstream ss;
+
+		//output the nodes
+		ss<<"Nodes----------------------------"<<std::endl;
+		for(t_node_i i = N.begin(); i != N.end(); i++){
+			ss<<i->str()<<std::endl;
+		}
+
+		//output the fibers
+		ss<<std::endl<<"Fibers---------------------------"<<std::endl;
+
+		T length, radius;
+		//output every fiber
+		for(fiber_i f = F.begin(); f != F.end(); f++){
+
+			//calculate the length and average radius
+			radius = f->radius(length);
+
+			//output the IDs of the terminal nodes
+			ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<",  average radius = "<<radius<<std::endl;
+		}
+
+		return ss.str();
+	}
+
+	/// Load a network from an OBJ object
+	void load( stim::obj<T> object){
+
+		//get the number of vertices in the object
+		unsigned int nV = object.numV();
+
+		//allocate an array of pointers to nodes, which will be used to preserve connectivity
+			//initiate all values to T.end()
+		std::vector< t_node_i > node_hash(nV, N.end());
+
+		unsigned int nL = object.numL();			//get the number of lines in the OBJ
+
+		//for each line in the OBJ structure
+		for(unsigned int li = 0; li < nL; li++){
+
+			F.push_back(fiber());				//push a new fiber onto the fiber list
+
+			fiber_i f = --(F.end());			//get an iterator to the new fiber
+
+			//----------Handle the terminating nodes for the fiber
+
+			//get the indices of the line vertices
+			std::vector< unsigned int > Li = object.getL_Vi(li);
+			unsigned int i0 = Li.front() - 1;
+			unsigned int i1 = Li.back() - 1;
+
+			//deal with the first end point of the capillary
+			if(node_hash[i0] != N.end()){			//if the node has been used before
+				(*f).n[0] = node_hash[i0];			//assign the node to the new capillary
+				(*node_hash[i0]).out.push_back(f);	//add an out pointer to the existing node
+			}
+			else{									//otherwise
+				N.push_back(t_node());				//create a new node and add it to the node list
+				t_node_i t = --(N.end());			//get an iterator to the new node
+				node_hash[i0] = t;					//add a pointer to the new node to the hash list
+				(*f).n[0] = t;						//add a pointer to the new node to the capillary
+				(*t).out.push_back(f);				//add a pointer to the capillary to the new node
+			}
+
+			//deal with the last end point of the capillary
+			if(node_hash[i1] != N.end()){
+				(*f).n[1] = node_hash[i1];
+				(*node_hash[i1]).in.push_back(f);
+			}
+			else{
+				N.push_back(t_node());
+				t_node_i t = --(N.end());
+				node_hash[i1] = t;			//add the new node to the hash list
+				(*f).n[1] = t;
+				(*t).in.push_back(f);
+			}
+
+			//-------------Handle the geometric points for the fiber
+			std::vector< vec<T> > L = object.getL_V(li);
+			std::vector< vec<T> > R = object.getL_VT(li);
+
+			unsigned int nP = L.size();				//get the number of geometric points in the fiber
+			//for each vertex in the fiber
+			for(unsigned int pi = 0; pi < nP; pi++){
+				point p = (point)L[pi];					//move the geometric coordinates into a point structure
+				p.r = R[pi][0];							//store the radius
+				f->P.push_back(p);						//push the point onto the current fiber
+			}
+		}
+
+	}	//end load()
+
+	/// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.
+
+	/// @param n0 is a array which will contain the list of source nodes
+	/// @param n1 is a array which will contain the list of destination nodes
+	/// @param length is a array containing the lengths of fibers in the network
+	/// @param radius is a array containing the average radii of fibers in the network
+	void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){
+
+		//determine the number of fibers in the network
+		unsigned int nF = F.size();
+
+		//allocate the necessary space to store the fiber information
+		n0.resize(nF);
+		n1.resize(nF);
+		length.resize(nF);
+		radius.resize(nF);
+
+		//assign names (identifiers) to the network components
+		set_names();
+
+		//fill the arrays
+		unsigned int i = 0;
+		T l, r;
+		for(fiber_i f = F.begin(); f != F.end(); f++){
+			n0[i] = f->n[0]->id;	//get the identifiers for the first and second nodes for the current fiber
+			n1[i] = f->n[1]->id;
+
+			r = f->radius(l);		//get the length and radius of the capillary (calculated at the same time)
+
+			radius[i] = r;			//store the radius in the output array
+			length[i] = l;			//store the length in the output array
+
+			i++;					//increment the array index
+		}
+
+
+	}
+
+};
+
+};	//end namespace stim
+
+
+#endif
@@ -600,9 +600,42 @@ public:
 		}
  
 		return l;
+	}
+
+	/// Returns a vector containing the list of texture coordinates associated with each point of a line.
+
+	/// @param i is the index of the desired line
+	std::vector< stim::vec<T> > getL_VT(unsigned int i){
+
+		//get the number of points in the specified line
+		unsigned int nP = L[i].size();
+
+		//create a vector of points
+		std::vector< stim::vec<T> > l;
+
+		//set the size of the vector
+		l.resize(nP);
+
+		//copy the points from the point list to the stim vector
+		unsigned int pi;
+		for(unsigned int p = 0; p < nP; p++){
+
+			//get the index of the geometry point
+			pi = L[i][p][1] - 1;
+
+			//get the coordinates of the current point
+			stim::vec<T> newP = VT[pi];
+
+			//copy the point into the vector
+			l[p] = newP;
+		}
+
+		return l;
  
 	}
  
+
+
 };