#ifndef STIM_NETWORK_H #define STIM_NETWORK_H #include #include #include #include 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 class network{ //-------------------HELPER CLASSES----------------------- /// Stores information about a geometric point on the network centerline (including point position and radius) //template class point : public stim::vec{ public: T r; point() : stim::vec(){} //casting constructor point(stim::vec rhs) : stim::vec(rhs){} }; //template 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 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 > geometry(){ std::vector< stim::vec > 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::iterator i; for(i = begin(); i != end(); i++){ ss<str()<<" r = "<r< 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<::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 N; //list of fibers std::list 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----------------------------"<str()<radius(length); //output the IDs of the terminal nodes ss<n[0]->id<<" -- "<n[1]->id<<": length = "<