stim/biomodels/network.h

#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H
#include <stim/math/vector.h>
#include <stim/visualization/obj.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 >{
		using std::list< point >::begin;
		using std::list< point >::end;
		using std::list< point >::size;
	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 = begin(); i != end(); i++){		//for each point in the fiber
				if(i == 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 = begin(); i != end(); i++){		//for each point in the fiber
				if(i == 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
			//if the total length is zero, store a radius of zero
			if(length == 0)
				return 0;
			else
				return radius_sum / length;					//return the average radius of the fiber
		}
		std::vector< stim::vec<T> > geometry(){
			std::vector< stim::vec<T> > result;				//create an array to store the fiber geometry
			result.resize( size() );					//pre-allocate the array
			typename std::list< point >::iterator p;				//create a list iterator
			unsigned int pi = 0;						//create an index into the result array
			//for each geometric point on the fiber centerline
			for(p = begin(); p != end(); p++){
				result[pi] = *p;
				pi++;
			}
			return result;			//return the geometry array
		}
		std::string str(){
			std::stringstream ss;
			//create an iterator for the point list
			typename std::list<point>::iterator i;
			for(i = begin(); i != 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();
		}
	};
//---------------NETWORK CLASS-----------------------------
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
	/// @param object is the object file to be used as the basis for the network
	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->push_back(p);						//push the point onto the current fiber
			}
		}
	}	//end load()
	/// Returns an array of node positions
	std::vector< stim::vec<T> > get_node_positions(){
		std::vector< stim::vec<T> > result;				//create an array to store the result
		result.resize(N.size());						//set the array size
		t_node_i ni;									//create a terminal node iterator
		unsigned int vi = 0;							//vertex index into the result array
		//for every terminal node
		for(ni = N.begin(); ni != N.end(); ni++){
			//create a vector based on the node position
			//if the number of outgoing nodes is nonzero
			if(ni->out.size() != 0)
				result[vi] = ni->out.front()->front();
			else if(ni->in.size() != 0)
				result[vi] = ni->in.front()->back();
			vi++;										//increment the array index
		}
		//return the resulting array
		return result;
	}
	std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
		return f->geometry();
	}
	/// Generate an OBJ file from the network
	stim::obj<T> obj(){
		//create an OBJ object
		stim::obj<T> object;
		//name the nodes
		set_names();
		//retrieve a list of terminal node positions
		std::vector< stim::vec<T> > node_pos = get_node_positions();
		//add the nodes to the obj file
		object.addV(node_pos);
		//counter for vertex indices in the object class
		unsigned int nP;
		//for each fiber
		fiber_i fi;						//create a fiber iterator
		for(fi = F.begin(); fi != F.end(); fi++){
			//get an array of fiber points
			std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
			//create a subset of this array
			typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;
			typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;
			typename std::vector< stim::vec<T> > fiber_subset(start, end);
			//add this subset to the geometry object
			nP = object.addV(fiber_subset);
			//create an array to hold vertex indices for a line
			std::vector<unsigned int> line;
			line.resize(fiber_p.size());
			//add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)
			line[0] = fi->n[0]->id + 1;
			line[line.size() - 1] = fi->n[1]->id + 1;
			//add the intermediate vertex indices to the line array
			for(unsigned int i = 0; i < fiber_subset.size(); i++){
				line[1 + i] = nP + i;
			}
			//add the line list to the object class
			object.addLine(line);
		}
		return object;
	}
	/// 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