codebase / stimlib

#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H

#include <stim/math/vector.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