Merge branch 'pranathi_stimlib' of https://git.stim.ee.uh.edu/codebase/stimlib into pranathi_stimlib

pranathivemuri
2 parents ad2de689 7f27eafa
Showing 4 changed files with 636 additions and 1306 deletions Show diff stats
stim/visualization/fiber.h → stim/biomodels/fiber.h
stim/biomodels/network.h
stim/biomodels/network_dep.h
stim/visualization/network.h
@@ -13,6 +13,7 @@ namespace stim{
 template<typename T>
 class fiber{
+protected:
 	unsigned int N;					//number of points in the fiber
 	double **c;						//centerline (array of double pointers)
@@ -70,11 +71,16 @@ class fiber{
 		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
 	}
+
 	double dist(double* p0, double* p1){
 		double sum = 0;
-		for(unsigned int d = 0; d < 3; d++)
-			sum += p0[d] * p1[d];
+		float v;
+		for(unsigned int d = 0; d < 3; d++){
+			v = p1[d] - p0[d];
+			sum +=v * v;
+
+		}
 		return sqrt(sum);
 	}
@@ -92,6 +98,19 @@ class fiber{
 		return *idx;
 	}
+	/// Returns a stim::vec representing the point at index i
+
+	/// @param i is an index of the desired centerline point
+	stim::vec<T> get_vec(unsigned i){
+		stim::vec<T> r;
+		r.resize(3);
+		r[0] = c[i][0];
+		r[1] = c[i][1];
+		r[2] = c[i][2];
+
+		return r;
+	}
+
@@ -114,6 +133,25 @@ public:
 		init(n);
 	}
+	/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
+	fiber(std::vector< stim::vec<T> > p){
+		init(p.size());		//initialize the fiber
+
+		//for each point, set the centerline position and radius
+		for(unsigned int i = 0; i < N; i++){
+
+			//set the centerline position
+			for(unsigned int d = 0; d < 3; d++)
+				c[i][d] = (double) p[i][d];
+
+			//set the radius
+			r[i] = 0;
+		}
+
+		//generate a kd tree
+		gen_kdtree();
+	}
+
 	/// constructor takes a list of points and radii
 	fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){
 		init(pos.size());		//initialize the fiber
@@ -440,6 +478,19 @@ public:
 	unsigned int n_pts(){
 		return N;
 	}
+
+	/// Bracket operator returns the element at index i
+
+	/// @param i is the index of the element to be returned as a stim::vec
+	stim::vec<T> operator[](unsigned i){
+		return get_vec(i);		
+	}
+
+	/// Back method returns the last point in the fiber
+	stim::vec<T> back(){
+		return get_vec(N-1);
+	}
+
 };
 #ifndef STIM_NETWORK_H
 #define STIM_NETWORK_H
+#include <list>
+#include <stdlib.h>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+#include <string.h>
+#include <math.h>
 #include <stim/math/vector.h>
 #include <stim/visualization/obj.h>
-#include <list>
+#include <stim/biomodels/fiber.h>
 #include <ANN/ANN.h>
+#include <boost/tuple/tuple.hpp>
+
 namespace stim{
+/** This is the a class that interfaces with gl_spider in order to store the currently
+ *   segmented network. The following data is stored and can be extracted:
+ *   1)Network geometry and centerline.
+ *   2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
+*/
-/** 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){
+	///Each edge is a fiber with two nodes.
+	///Each node is an in index to the endpoint of the fiber in the nodes array.
+	class edge : public fiber<T>
+	{
+		public:
+		unsigned v[2];		//unique id's designating the starting and ending
-			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
-
-		}
+		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
+		///@param p is a position in space
+		edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
+			
+		/// Output the edge information as a string
 		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;
-			}
-
+			ss<<"("<<fiber<T>::N<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
 			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++){
+	};	
+	
+	///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
+	class vertex : public stim::vec<T>
+	{
+		public:
+			//std::vector<unsigned int> edges;		//indices of edges connected to this node.
+			std::vector<unsigned int> e[2];			//indices of edges going out (e[0]) and coming in (e[1])
+			//stim::vec<T> p;							//position of this node in physical space.
+
+			//constructor takes a stim::vec
+			vertex(stim::vec<T> p) : stim::vec<T>(p){}
+
+			/// Output the vertex information as a string
+			std::string	str(){
+				std::stringstream ss;
+				ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
+
+				if(e[0].size() > 0){
+					ss<<"\t> ";
+					for(unsigned int o = 0; o < e[0].size(); o++)
+						ss<<e[0][o]<<" ";
+				}
+				if(e[1].size() > 0){
+					ss<<"\t< ";
+					for(unsigned int i = 0; i < e[1].size(); i++)
+						ss<<e[1][i]<<" ";
+				}
-				if(f != out.begin())
-					ss<<", ";
-				ss<<(*f)->n[1]->id;
+				return ss.str();
 			}
-
-
-			return ss.str();
-
-
-
-		}
+			
 	};
+	private:
-//---------------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++;
+	std::vector<edge> E;       //list of edges
+	std::vector<vertex> V;	    //list of vertices.
+	
+	public:
-		i = 0;
-		for(fiber_i fi = F.begin(); fi != F.end(); fi++)
-			fi->id = i++;
+	///Returns the number of edges in the network.
+	unsigned int edges(){
+		return E.size();
 	}
-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();
+	///Returns the number of nodes in the network.
+	unsigned int vertices(){
+		return V.size();
 	}
-	/// Load a network from an OBJ object
+	//load a network from an OBJ file
+	void load_obj(std::string filename){
-	/// @param object is the object file to be used as the basis for the network
-	void load( stim::obj<T> object){
+		stim::obj<T> O;			//create an OBJ object
+		O.load(filename);		//load the OBJ file as an object
-		//get the number of vertices in the object
-		unsigned int nV = object.numV();
+		//grab each line in the OBJ object - these will be fibers
+		std::cout<<O.str()<<std::endl;
-		//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());
+		std::vector<unsigned> id2vert;	//this list stores the OBJ vertex ID associated with each network vertex
-		unsigned int nL = object.numL();			//get the number of lines in the OBJ
+		unsigned i[2];						//temporary, IDs associated with the first and last points in an OBJ line
-		//for each line in the OBJ structure
-		for(unsigned int li = 0; li < nL; li++){
+		//for each line in the OBJ object
+		for(unsigned int l = 1; l <= O.numL(); l++){
-			F.push_back(fiber());				//push a new fiber onto the fiber list
+			std::vector< stim::vec<T> > c;				//allocate an array of points for the vessel centerline
+			O.getLine(l, c);							//get the fiber centerline
-			fiber_i f = --(F.end());			//get an iterator to the new fiber
+			edge e(c);									//create an edge from the given centerline
-			//----------Handle the terminating nodes for the fiber
+			//get the first and last vertex IDs for the line
+			std::vector< unsigned > id;					//create an array to store the centerline point IDs
+			O.getLinei(l, id);							//get the list of point IDs for the line			
+			i[0] = id.front();							//get the OBJ ID for the first element of the line
+			i[1] = id.back();							//get the OBJ ID for the last element of the line
-			//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;
+			std::vector<unsigned>::iterator it;			//create an iterator for searching the id2vert array
+			unsigned it_idx;							//create an integer for the id2vert entry index
-			//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
+			//find out if the nodes for this fiber have already been created
+			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
+			it_idx = std::distance(id2vert.begin(), it);
+			if(it == id2vert.end()){							//if i[0] hasn't already been used
+				vertex v = e[0];								//create a new vertex, assign it a position
+				v.e[0].push_back(E.size());						//add the current edge as outgoing
+				e.v[0] = V.size();								//add the new vertex to the edge
+				V.push_back(v);									//add the new vertex to the vertex list
+				id2vert.push_back(i[0]);						//add the ID to the ID->vertex conversion list
 			}
-			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
+			else{												//if the vertex already exists
+				V[it_idx].e[0].push_back(E.size());					//add the current edge as outgoing
+				e.v[0] = it_idx;
 			}
-			//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);
+			it = find(id2vert.begin(), id2vert.end(), i[1]);	//look for the second ID
+			it_idx = std::distance(id2vert.begin(), it);
+			if(it == id2vert.end()){							//if i[1] hasn't already been used
+				vertex v = e.back();							//create a new vertex, assign it a position
+				v.e[1].push_back(E.size());						//add the current edge as incoming
+				e.v[1] = V.size();
+				V.push_back(v);									//add the new vertex to the vertex list
+				id2vert.push_back(i[1]);						//add the ID to the ID->vertex conversion list
 			}
-			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);
+			else{												//if the vertex already exists
+				V[it_idx].e[1].push_back(E.size());					//add the current edge as incoming
+				e.v[1] = it_idx;
 			}
-			//-------------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);
+			E.push_back(e);										//push the edge to the list
-			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);
+	/// Output the network as a string
+	std::string str(){
+		std::stringstream ss;
+		ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
+		for(unsigned int v = 0; v < V.size(); v++){
+			ss<<"\t"<<v<<V[v].str()<<std::endl;
 		}
-		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
+		ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
+		for(unsigned e = 0; e < E.size(); e++){
+			ss<<"\t"<<e<<E[e].str()<<std::endl;
 		}
-
+		return ss.str();
 	}
-
-};
-
-};	//end namespace stim
-
-
+};		//end stim::network class
+};		//end stim namespace
 #endif
+#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
-#ifndef STIM_NETWORK_H
-#define STIM_NETWORK_H
-
-#include <list>
-#include <stdlib.h>
-#include <sstream>
-#include <fstream>
-#include <algorithm>
-#include <string.h>
-#include <math.h>
-#include <stim/math/vector.h>
-#include <stim/visualization/obj.h>
-#include <stim/visualization/fiber.h>
-#include <ANN/ANN.h>
-#include <boost/tuple/tuple.hpp>
-
-
-namespace stim{
-/** This is the a class that interfaces with gl_spider in order to store the currently
- *   segmented network. The following data is stored and can be extracted:
- *   1)Network geometry and centerline.
- *   2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
-*/
-
-
-template<typename T>
-class network{
-
-	///Each edge is a fiber with two nodes.
-	///Each node is an in index to the endpoint of the fiber in the nodes array.
-	class edge : public fiber<T>
-	{
-		public:
-		int Nodes[2];		//unique id's designating the starting and ending
-
-		///default constructor
-		edge() : fiber<T>()
-		{
-			Nodes[1] = -1; Nodes[2] = -1;
-		}
-
-		///sized costructor with two known nodes.
-		///@param startId: int storing idx assigned to Nodes[1].
-		///@param endId: int storing idx assigned to Nodes[2].
-		///@param n: int for the number of points in the fiber.
-		edge(int startId, int endId, int n)
-		 :fiber<T>(n)
-		{
-			Nodes[1] = startId; Nodes[2] = endId;
-		}
-		
-		///constructor from a std::vector of stim::vecs of positions and radii.
-		///@param pos: Vector of stim vecs storing the positions of the fiber.
-		///@param mag: Vector of stim vecs storing the radii of the fiber.
-		edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii)
-			 : fiber<T>(pos, radii)
-		{
-			Nodes[1] = -1; Nodes[2] = -1;
-		}
-		
-		///constructor from an std::vector of stim::vecs of positions and radii as well as a known starting and ending node.
-		///@param pos: Vector of stim vecs storing the positions of the fiber.
-		///@param mag: Vector of stim vecs storing the radii of the fiber.
-		///@param id1: int storing idx assigned to Nodes[1].
-		///@param id2: int storing idx assigned to Nodes[2].
-		edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii, int id1, int id2)
-			 : fiber<T>(pos, radii)
-		{
-			Nodes[1] = id1; Nodes[2] = id2;
-		}
-
-
-		edge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
-			 : fiber<T>(pos, radii)
-		{
-			Nodes[1] = -1; Nodes[2] = -1;
-		}
-
-		///sets the endpoints to the two int values.
-		///@param int id1: index of Nodes[1].
-		///@param int id2: index of Nodes[2].
-		void
-		setEndpoints(int id1, int id2)
-		{
-			Nodes[1] = id1; Nodes[2] = id2;
-		}		
-
-	};	
-	
-	///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
-	class node
-	{
-		private:
-			std::vector<int> edges;		//indices of edges connected to this node.
-			stim::vec<T> p;			//position of this node in physical space.
-		public:
-			//no default constructor
-
-			///@param pos: stim vec with the x, y, z position of the edge.
-			///stim::vec constructure with a position but no attached edges.
-			node(stim::vec<T> pos)
-			{
-				p = pos;
-			}
-
-			///@param pos: stim vec with the x, y, z position of the edge.
-			///@param i: int i storing the index of an attached edge.
-			//stim::vec constructor with a position and an attached edge.
-			node(stim::vec<T> pos, int i)
-			{
-				p = pos;
-				edges.push_back(i);
-			}
-
-			///@param x: x coordinate of the node..
-			///@param y: y coordinate of the node..
-			///@param z: z coordinate of the node..
-			//float value constructor. 
-			node(T x, T y, T z)
-			{
-				p = stim::vec<T>(x,y,z);
-			}
-			///@param x: x coordinate of the node..
-			///@param y: y coordinate of the node..
-			///@param z: z coordinate of the node..
-			///@param i: int i storing the index of an attached edge.
-			//float value constructor with an attached edge.
-			node(T x, T y, T z, int i)
-			{
-				p = stim::vec<T>(x,y,z);
-				edges.push_back(i);
-			}
-
-			///@param id: int index of the fiber to attach to this node.
-			///adds the fiber to the rest of the fibers connected to this node.
-			void
-			addEdge(int id)
-			{
-			 	edges.push_back(id);
-			}
-	
-			///returns the position of the node.
-			stim::vec<T>
-			getPosition()
-			{	
-				return p;
-			}
-		
-			///@param id: int index of the fiber to detach to this node.
-			///removes the edge from the list of the edges attached to this node.
-			void
-			removeEdge(int id)
-			{
-				for(int i = 0; i < edges.size(); i++)
-				{
-					if(edges[i] == id)
-						edges.erase(edges.begin()+i);
-				}
-			}
-			
-			///returns and std::string with the position of this node.
-			std::string
-			str()
-			{
-				return p.str();
-			}
-
-			///returns and std::string with the list of edges attached to this node.
-			std::string
-			edges_to_str()
-			{
-				std::ostringstream ss;
-				std::cout<<"here"<<std::endl;
-//				ss << "[";
-				for(int i = 0; i < edges.size()-1; i++)
-				{
-					std::cout<<"here"<<i<<std::endl;
-					ss << edges[i] << ";";
-					std::cout<<edges.size()-1<<std::endl;
-				}
-				ss << edges[edges.size()-1];
-//				ss << "]";
-				return ss.str();
-			}
-			
-	};
-
-	public:
-
-	std::vector<edge*> E;       //list of pointers to edges.
-	std::vector<node> V;	    //list of nodes.
-	std::vector< stim::vec<T> > allVerticesList; // all nodes before sampling
-	std::vector<stim::vec<T> > allVerticesAfterSampling ; //all vertices after sampling
-	
-	///Returns the number of edges in the network.
-	unsigned int
-	sizeE()
-	{
-		return E.size();
-	}
-
-	///Returns the number of nodes in the network.
-	unsigned int
-	sizeV()
-	{
-		return V.size();
-	}
-	unsigned int
-	sizeAfterSamplingV()
-	{
-		return allVerticesAfterSampling.size();
-	}
-/*	//adds an edge from two std::vectors with positions and radii.
-	void
-	addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int endId)
-	{
-		
-		edge a(pos,radii, endId);
-		E.push_back(a);
-	} */
-
-	///A complicated method that adds an edge to the network.
-	///Most of this functionality will be moved into fiber.
-	void
-	addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int startId, int endId)
-	{
-		//
-		if(startId == -1 && endId == -1)
-		{
-			//case where the edge is neither starting nor ending in a fiber.
-			//i. first fiber.
-			
-			//Add two nodes to the nodes vector
-			V.push_back(node(pos[pos.size()-1]));
-			V.push_back(node(pos[0]));
-
-			//the edge will be connected to the nodes
-			edge *a = new edge(pos,radii,(V.size()-2), (V.size()-1));
-	
-			//add fiber to fiber list.
-			E.push_back(a);
-
-			//The last two nodes added to V will "own" the last edge added to E.
-			V[V.size()-1].addEdge(E.size()-1);
-			V[V.size()-2].addEdge(E.size()-1);
-			
-		}
-	
-		else if(startId != -1 && endId == -1)
-		{
-			//case where the edge is starting with a fiber, but not ending in one.
-			
-			//split the fiber behind you into two.
-			unsigned int point = E[startId]->nearest_idx(pos[0]);
-
-			//split the hit fiber at point two parts temp[0], temp[1]
-			std::vector < stim::fiber <T> > temp = E[startId]->split(point);
-			if(temp.size() > 1)
-			{
-				//add the nearest point in the behind fiber into the hitting fiber.
-				pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
-				stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
-				radii.insert(radii.begin(), v);
-
-				//reset the fiber at the endId to be a shorter fiber(temp[0]).
-				std::vector<stim::vec<T> > ce = temp[0].centerline();
-				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
-
-				//remake the edge, such that the starting point of this edge
-				//is the same the split point, but the end edge is the old end.
-				V.push_back(node(ce[ce.size()-1]));
-				int tempNodeId = E[startId]->Nodes[1];
-				E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
-				V[V.size()-1].addEdge(startId);
-				
-
-				//add the other part of the fiber (temp[1])
-				ce = temp[1].centerline();
-				cm = temp[1].centerlinemag();
-				E.push_back(new edge(ce, cm,tempNodeId ,(V.size()-1)));
-				V[V.size()-1].addEdge(E.size()-1);
-
-				V[tempNodeId].removeEdge(startId);
-				V[tempNodeId].addEdge(E.size()-1);
-//				V[V.size()-1].removeEdge(startId);
-
-				//make the new hitting fiber..
-				V.push_back(node(pos[pos.size()-1]));
-				edge *a = new edge(pos, radii, (V.size()-2), (V.size()-1));
-				E.push_back(a);
-				V[V.size()-1].addEdge(E.size()-1);
-				V[V.size()-2].addEdge(E.size()-1);
-			} else {
-				stim::vec<T> pz = E[startId]->nearest(pos[0]);
-				if((V[E[startId]->Nodes[1]].getPosition() - pz).len() <
-					(V[E[startId]->Nodes[2]].getPosition() - pz).len())
-				{
-					unsigned int point = E[startId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
-					stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[startId]->Nodes[1], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[startId]->Nodes[1]].addEdge(E.size()-1);
-
-				} 
-				else
-				{
-					unsigned int point = E[startId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
-					stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[startId]->Nodes[2], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[startId]->Nodes[2]].addEdge(E.size()-1);
-				}
-			}
-		}	
-
-		//case where the edge is starting at a seedpoint but ends in a fiber.
-		if(startId == -1 && endId != -1 && endId < sizeE())
-		{	
-			//split the hit fiber into two.
-			unsigned int point = E[endId]->nearest_idx(pos[pos.size() -1]);
-			
-			//split the hit fiber at point into two parts temp[0], temp[1]
-			std::vector < stim::fiber <T> > temp
-				 = E[endId]->split(point);
-			if(temp.size() > 1)
-			{
-				//add the nearest point in the hit fiber into the hitting fiber.
-				pos.push_back(E[endId]->nearest(pos[pos.size()-1]));
-	//			pos.insert(pos.end(), E[endId].nearest(pos[pos.size()-1]));
-				stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
-				radii.push_back(v);
-
-				//split the hit fiber at point into two parts temp[0], temp[1]
-				std::vector < stim::fiber <T> > temp
-					 = E[endId]->split(point);
-				
-				//reset the fiber at endId to be a shorter fiber (temp[0]).
-				std::vector<stim::vec<T> > ce = temp[0].centerline();
-				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
-
-				//remake the edge, such that the ending point of this edge
-				//is the same as before, but the starting edge is new.
-				V.push_back(node(ce[ce.size()-1])); //split point.
-				int tempNodeId = E[endId]->Nodes[2];
-				E[endId] = new edge(ce, cm, E[endId]->Nodes[1], (V.size()-1));
-				V[V.size()-1].addEdge(endId);
-
-				//add that other part of the fiber (temp[1])
-				//such that it begins with the middle node, and ends with 
-				//the last node of the previous fiber.
-				ce = temp[1].centerline();
-				cm = temp[1].centerlinemag();
-				E.push_back(new edge(ce, cm, (V.size()-1), tempNodeId));
-				V[V.size()-1].addEdge(E.size()-1);
-	//			V[V.size()-1].removeEdge(endId);
-				
-				
-		
-				//make the new hitting fiber..
-				V.push_back(pos[0]);
-				edge *a = new edge(pos,radii,(V.size()-1), (V.size()-2));
-				E.push_back(a);
-				V[V.size()-1].addEdge(E.size()-1);
-				V[V.size()-2].addEdge(E.size()-1);
-
-				//in the end we have added two new nodes and two new edges.
-			}
-			else {
-				stim::vec<T> pz = E[endId]->nearest(pos[0]);
-				if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
-					(V[E[endId]->Nodes[2]].getPosition() - pz).len())
-				{
-					unsigned int point = E[endId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
-					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[endId]->Nodes[1]].addEdge(E.size()-1);
-
-				} 
-				else
-				{
-					unsigned int point = E[endId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
-					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
-				}
-			}
-		}
-
-		if(startId != -1 && endId != -1 && endId < sizeE())
-		{
-			//case where the edge is starting with a fiber, and ends in one.
-			
-			//split the fiber behind you into two.
-			unsigned int point = E[startId]->nearest_idx(pos[0]);
-//			std::cout << "in merge to both" << std::endl;
-
-			//split the hit fiber at point two parts temp[0], temp[1]
-			std::vector < stim::fiber <T> > temp = E[startId]->split(point);
-			if(temp.size() > 1)
-			{
-				//extend the previous fiber (guaranteed to be added last) by one
-				//and add the
-				pos = E[E.size()-1]->centerline();
-				radii = E[E.size()-1]->centerlinemag();
-				pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
-				stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
-				radii.insert(radii.begin(), v);
-				V.erase(V.end());
-				V.push_back(node(pos[0]));
-
-
-				//something weird is going on here.
-				E[E.size()-1] = new edge(pos, radii, (V.size()-2), (V.size()-1));
-				V[V.size()-1].addEdge(E.size()-1);
-
-				//reset the fiber at the endId to be a shorter fiber(temp[0]).
-				std::vector<stim::vec<T> > ce = temp[0].centerline();
-				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
-
-//				std::cout << ce.size() << std::endl;
-				//remake the edge, such that the starting point of this edge
-				//is the same as before, but the end edge is new.
-				int tempNodeId = E[startId]->Nodes[1];
-				E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
-				V[V.size()-1].addEdge(startId);
-
-				//add the other part of the fiber (temp[1])
-				ce = temp[1].centerline();
-				cm = temp[1].centerlinemag();
-				E.push_back(new edge(ce, cm,tempNodeId, (V.size()-1)));
-				V[V.size()-1].addEdge(E.size()-1);
-				V[tempNodeId].removeEdge(startId);
-				V[tempNodeId].addEdge(E.size()-1);
-//				V[V.size()-1].removeEdge(startId);
-			}
-			else {
-				stim::vec<T> pz = E[endId]->nearest(pos[0]);
-				if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
-					(V[E[endId]->Nodes[2]].getPosition() - pz).len())
-				{
-					unsigned int point = E[endId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
-					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[endId]->Nodes[1]].addEdge(E.size()-1);
-
-				} 
-				else
-				{
-					unsigned int point = E[endId]->nearest_idx(pos[0]);
-					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
-					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
-					radii.insert(radii.begin(), v);
-					V.push_back(node(pos[pos.size()-1]));
-					edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
-					E.push_back(a);
-					
-					V[V.size()-1].addEdge(E.size()-1);
-					V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
-				}
-			}
-		}
-		
-	}
-
-	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
-	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
-	///returns the forest as a std::string. For testing only.
-	std::string
-	edges_to_str()
-	{
-		std::stringstream ss;
-		for(unsigned int i = 0; i < E.size(); i++)
-		{
-			ss << i << ": " << E[i]->str() << std::endl;
-		}
-		return(ss.str());
-	}
-	// total number of points on all edges!=fibers in the network
-	int
-	totalEdges()
-	{
-		int totalEdgeVertices=0;int N=0;
-		for (unsigned int i=0; i < sizeE(); i ++)
-		{// FOR N points on the fiber N-1 edges are possible
-			N = E[i]->n_pts();
-			totalEdgeVertices = totalEdgeVertices + N- 1;
-		}
-		return totalEdgeVertices;
-	}
-	// sum of all the fiber lengths
-	float
-	lengthOfNetwork()
-	{
-		float networkLength=0;float N=0;
-		for (unsigned int i=0; i < sizeE(); i ++)
-		{
-			N = E[i]->length();
-			networkLength = networkLength + N;
-		}
-		return networkLength;
-	}
-	///@param i: index of the required fiber.
-	///Returns an std::vector with the centerline of the ith fiber in the edgelist.
-	std::vector< stim::vec<T> >
-	getEdgeCenterLine(int i)
-	{
-		std::vector < stim::vec<T> > a;
-		return E[i]->centerline();
-	}
-
-	///@param i: index of the required fiber.
-	///Returns an std::vector with the centerline radii of the ith fiber in the edgelist.
-	std::vector< stim::vec<T> >
-	getEdgeCenterLineMag(int i)
-	{
-		std::vector < stim::vec<T> > a;
-		return E[i]->centerlinemag();
-	}
-	
-	///@param i: index of the required fibers nodes..
-	///Returns an std::string with the start and end points of this node..
-	std::string
-	nodes_to_str(int i)
-	{
-		std::stringstream ss;
-		ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]";
-		return ss.str();
-	}
-	//load an obj file into a network class
-	stim::network<T>
-	objToNetwork(stim::obj<T> objInput)
-	{
-		stim::network<T> nwc;
-		//network network2;
-		// function to convert an obj file loaded using stim/visualization/obj.h
-		// to a 3D network class using network.h methods.
-		std::vector< stim::vec<T> > fiberPositions; //initialize positions on the fiber to be added to network
-		objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file
-		// to verify if the nodes are already pushed into node list
-		std::vector<bool> validList;
-		validList.assign(objInput.numV(), false);
-		// go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes
-		// using addNode function for adding nodes
-		// and the line as an edge(belongs to fiber class) using addEdge function
-		std::vector<unsigned> vertexIndices(objInput.numV());
-		std::vector< stim::vec<T> > vertices = objInput.getAllV(vertexIndices);
-		nwc.addVertices(vertices);
-		for (unsigned i =1; i< objInput.numL() + 1; i++) 
-			{
-			    // edges/fibers could be added to the network class
-				std::vector< stim::vec<T> > fiberPositions;
-
-			    objInput.getLine(i, fiberPositions);
-			    // finding size to allocate radii
-			    int numPointsOnFiber = fiberPositions.size();
-			    // to extend it to a 3D postion if it is a 1D/2D vertex in space
-				std::vector< stim::vec<T> > fiberPositions1(numPointsOnFiber);
-			   // 2D case append and assign the last dimension to zero
-			   if (numDim == 2)
-			   {//  2D case append and assign the last dimension to zero repeating it for all the points on fiber
-				   for (int j = 0; j < numPointsOnFiber; j ++)
-				   {
-					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
-					   fiberPositions1[j][numDim -1 ] = fiberPositions[j][1];
-					   fiberPositions1[j][numDim] = 0;
-				   }
-			   }
-			   // 1D case append and assign last two dimensions to zero
-			   else if (numDim == 1)
-			   {
-				   for (int j = 0; j < numPointsOnFiber; j ++)
-				   {
-					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
-					   fiberPositions1[j][numDim -1 ] = 0;
-					   fiberPositions1[j][numDim] = 0;
-				   }
-			   }
-			   else
-			   {
-				   fiberPositions1 = fiberPositions;
-			   }
-		   		// then add edge 
-			    //edge* a = new edge(fiberPositions1,radii); 
-				//std::cout<<"here"<<std::endl;
-			    //E.push_back(a);
-			    std::vector<stim::vec<T> > newPointList;
-			    newPointList  = Resample(fiberPositions1);
-				int numPointsOnnewFiber = newPointList.size();
-				nwc.addVerticesAfterSamplimg(newPointList);
-				std::vector<T> radii(numPointsOnnewFiber); // allocating radii to zero
-				nwc.addEdge(newPointList, radii);
-				// add nodes
-				stim::vec<T> v0(3);stim::vec<T> vN(3);
-				int endId = numPointsOnnewFiber -1;
-				v0[0] = newPointList[0][0];v0[1] = newPointList[0][1];v0[2] = newPointList[0][2];
-				vN[0] = newPointList[endId][0];vN[1] = newPointList[endId][1];vN[2] = newPointList[endId][2];
-				// VISITED INDEXES OF the nodes are set to true
-				if(!validList[objInput.getIndex(vertices, v0)])
-				{
-					//V.push_back(node(v0));
-					nwc.addNode(v0);
-					validList[objInput.getIndex(vertices, v0)] = true;
-				}
-				if(!validList[objInput.getIndex(vertices, vN)])
-				{
-					//V.push_back(node(vN));
-					nwc.addNode(vN);
-					validList[objInput.getIndex(vertices, vN)] = true;
-				}
-			}
-		return nwc;
-		}
-	// convert ground and test case to network ,kd trees
-	boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<T>, stim::network<T> > 
-	LoadNetworks(std::string gold_filename, std::string test_filename)
-	{
-		ANNkd_tree* kdTree1;ANNkd_tree* kdTree2;
-		using namespace stim;
-		network network1;network network2;
-		stim::obj<T> objFile1(gold_filename);
-		network1 = objToNetwork(objFile1);
-		kdTree1 = generateKdTreeFromNetwork(network1);
-		std::cout<<"Gold Standard:network and kdtree generated"<<std::endl;
-		stim::obj<T> objFile2(test_filename);
-		network2 = objToNetwork(objFile2);
-		kdTree2 = generateKdTreeFromNetwork(network2);
-		std::cout<<"Test Network:network and kdtree generated"<<std::endl;
-		return boost::make_tuple(kdTree1, kdTree2, network1, network2);
-		std::cout<<"out of this loop"<<std::endl;
-	}
-	//distance between two points
-	double dist(std::vector<T> p0, std::vector<T> p1)
-	{
-
-		double sum = 0;
-		for(unsigned int d = 0; d < 3; d++)
-			sum += p0[d] * p1[d];
-		return sqrt(sum);
-	}
-	// sum of elements in a vector
-	double sum(stim::vec<T> metricList)
-	{
-		float sumMetricList = 0;
-		for (unsigned int count=0; count<metricList.size(); count++)
-		{ sumMetricList += metricList[count];}
-		return sumMetricList;
-	}
-	//generate a kd tree from network
-	ANNkd_tree* 
-	generateKdTreeFromNetwork(stim::network<T> network)			//kd-tree stores all points in the fiber for fast searching
-	{       
-		std::cout<<"kd trees generated"<<std::endl;
-		ANNkd_tree* kdt;
-		double **c;						//centerline (array of double pointers)
-		float* r;						   // array of fiber radii
-		unsigned int n_data = network.sizeAfterSamplingV(); //set the number of points
-		c = (double**) malloc(sizeof(double*) * n_data);  //allocate the array pointer
-		for(unsigned int i = 0; i < n_data; i++)		 //allocate space for each point
-			{c[i] = (double*) malloc(sizeof(double) * 3);}	
-		r = (float*) malloc(sizeof(float) * n_data);			//allocate space for the radii
-		//stim::vec<T> node;
-		for(unsigned int i = 0; i < n_data; i++)
-			{
-				//node = network.V[i].getPosition();
-				//convert_to_double
-				for(unsigned int d = 0; d < 3; d++){		//for each dimension
-					c[i][d] = double(network.allVerticesAfterSampling[i][d]);				//copy the coordinate
-				}
-			r[i] = r[i];						//copy the radius
-			}
-		ANNpointArray pts = (ANNpointArray)c;           //create an array of data points
-		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
-		return kdt;
-	}
-
-	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
-	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
-	///adds an edge from two std::vectors with positions and radii.
-	void
-	addEdge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
-	{
-		edge *a = new edge(pos,radii);
-		E.push_back(a);
-	}
-	void
-	addNode(stim::vec<T> nodes)
-	{
-		V.push_back(nodes);
-	}
-	void
-	addVertices(std::vector< stim::vec<T> > vertices)
-	{
-		for (unsigned int i=0; i < vertices.size(); i ++)
-		{
-			allVerticesList.push_back(vertices[i]);
-		}
-	}
-	void
-	addVerticesAfterSamplimg(std::vector< stim::vec<T> > vertices)
-	{
-		for (unsigned int i=0; i < vertices.size(); i ++)
-		{
-			allVerticesAfterSampling.push_back(vertices[i]);
-		}
-	}
-	// gaussian function
-	float gaussianFunction(float x, float std=25)
-	{
-		float evaluate = 1.0f - ((exp(-x/(2*std*std))));
-		return evaluate;
-	}
-	// compares point on a network to a kd tree for two skeletons
-	boost::tuple< float, float > 
-	compareSkeletons(boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<float>, stim::network<float> > networkKdtree)
-	{
-		float gFPR, gFNR;
-		gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>());
-		gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>());
-		return boost::make_tuple(gFPR, gFNR);
-	}
-	// gaussian distance of points on network to Kdtree
-	float 
-	CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network<T> networkTruth)
-	{
-		std::vector< stim::vec<T> > fiberPoints;
-		float netmetsMetric=0;
-		double eps = 0; // error bound
-		ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances
-		ANNdistArray dists2; dists2 = new ANNdist[1];// near neighbor distances
-		ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
-		ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
-		int N; int numQueryPoints = networkTruth.totalEdges();
-		float totalNetworkLength = networkTruth.lengthOfNetwork();
-		stim::vec<float> fiberMetric(networkTruth.sizeE());
-		//for each fiber
-		for (unsigned int i=0; i < networkTruth.sizeE(); i ++)
-		{
-			std::vector<T> p1(3); std::vector<T> p2(3);//temporary variables to store point positions
-			fiberPoints = networkTruth.E[i]->centerline();
-			N = networkTruth.E[i]->n_pts();
-			stim::vec<float> segmentMetric(N-1);
-			// for each segment in the fiber
-			for (unsigned int j = 0; j < N - 1; j++)
-			{
-				ANNpoint queryPt1; queryPt1 = annAllocPt(3);
-				ANNpoint queryPt2; queryPt2 = annAllocPt(3);
-				//for each dimension of the point
-				for(unsigned int d = 0; d < 3; d++)
-				{ 
-					queryPt1[d] = double(fiberPoints[j][d]);
-					p1[d] = double(fiberPoints[j][d]);
-					queryPt2[d] = double(fiberPoints[j + 1][d]);
-					p2[d] = double(fiberPoints[j + 1][d]);// for each dimension
-				}
-				kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound
-				kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound
-				std::cout<<float(dists1[0])<<"dist1-----"<<float(dists2[0])<<"dist2---"<<std::endl;
-				float dist1 = gaussianFunction(float(dists1[0]));float dist2 = gaussianFunction(float(dists2[0]));
-				segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ;
-			}
-			fiberMetric[i] = sum(segmentMetric)/;
-		}
-		netmetsMetric =  sum(fiberMetric)/totalNetworkLength ;
-		return netmetsMetric;	
-	}
-    //resample a fiber in the network
-	std::vector<stim::vec<T> > Resample(std::vector< stim::vec<T> > fiberPositions, float spacing=25.0)
-	{
-		std::vector<T> p1(3), p2(3), v(3);
-		stim::vec<T> p(3);
-		std::vector<stim::vec<T> > newPointList;
-		for(unsigned int f=0; f<fiberPositions.size()-1; f++)
-		{
-			T lengthSegment = dist(p1,p2);
-			for(int d=0; d<3;d++)
-			{
-				p1[d] = T(fiberPositions[f][d]);
-			    p2[d] = T(fiberPositions[f + 1][d]);
-			}// for each dimension
-			if( lengthSegment >= spacing )
-				{	for (int dim=0; dim<3;dim++) //find the direction of travel
-						{v[dim] = p2[dim] - p1[dim];}
-					//set the step size to the voxel size
-					T step;
-					for(step=0.0; step<lengthSegment; step+=spacing)
-					{
-						for(int i=0; i<3;i++)
-							{p[i] = p1[i] + v[i]*(step/lengthSegment);}
-						newPointList.push_back(p);
-					}
-				}
-			else 
-				newPointList = fiberPositions;
-			}
-		return newPointList;
-	}
-	// modulus of a vector
-	T
-	Length(std::vector<T> v)
-	{
-	 	T sum=0;
-        for (int i=0; i<v.size(); i++)
-			sum += v[i] * v[i];
-	    return sum;
-	}
-	// function to remove characters from a string
-	void removeCharsFromString(std::string &str, char* charsToRemove ) {
-	   for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
-		  str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
-	   }
-	}
-	///exports the network graph to obj
-	void
-	to_obj()
-	{
-		std::ofstream ofs;
-		ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app);
-		int num = allVerticesList.size();
-		std::string *strArray = new std::string[num];
-		for(unsigned int i = 0; i < allVerticesList.size(); i++)
-		{
-			std::string str;
-			str = allVerticesList[i].str();
-			removeCharsFromString(str, "[],");
-			ofs << "v " << str << "\n";
-			removeCharsFromString(str," ");
-			strArray[i] = str;
-		}
-		for(unsigned int i = 0; i < E.size(); i++)
-		{
-			 std::string str;
-			 str = E[i]->strObj(strArray, num);
-			 //removeCharsFromString(str,"0");
-             ofs << "l " << str << "\n";
-		} 
-		ofs.close();
-	}	
-		///exports the graph.
-	void
-	to_csv()
-	{
-		std::ofstream ofs;
-		ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app);
-		std::cout<<"here"<<std::endl;
-		for(int i = 0; i < V.size(); i++)
-		{
-			std::cout<<"here"<<i<<std::endl;
-			ofs << V[i].edges_to_str() << "\n";
-		}
-		ofs.close();
-	}
-	///exports the graph.
-	void
-	to_gdf()
-	{
-		std::ofstream ofs;
-		ofs.open("Graph.gdf", std::ofstream::out | std::ofstream::app);
-		ofs << "nodedef>name VARCHAR\n";
-		for(int i = 0; i < V.size(); i++)
-		{
-			ofs << i << "\n";
-		}
-		ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n";
-		for(int i = 0; i < E.size(); i++)
-		{
-			ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <<E[i]->n_pts()
-			 << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n";
-		} 
-		ofs.close();
-	}
-
-};
-};
-#endif