stim/biomodels/network_dep.h

/*
Copyright <2017> <David Mayerich>
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#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