From 8761b649ae2c4f4d0921fff6bff6a1a4aa465508 Mon Sep 17 00:00:00 2001
From: pgovyadi <pgovyadi@gmail.com>
Date: Wed, 6 Apr 2016 16:04:12 -0500
Subject: [PATCH] moved a lot of the tracing functionality into the spider class. From the outside just attach spider and run trace()

---
 stim/biomodels/network.h              | 597 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 stim/gl/gl_spider.h                   | 141 ++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------
 stim/visualization/aaboundingbox.h    |  72 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 stim/visualization/gl_aaboundingbox.h |  60 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 395 insertions(+), 475 deletions(-)
 create mode 100644 stim/visualization/aaboundingbox.h
 create mode 100644 stim/visualization/gl_aaboundingbox.h

diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h
index f604361..29f4fe5 100644
--- a/stim/biomodels/network.h
+++ b/stim/biomodels/network.h
@@ -1,459 +1,312 @@
 #ifndef STIM_NETWORK_H
 #define STIM_NETWORK_H
 
+#include <stdlib.h>
+#include <assert.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/visualization/cylinder.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){
-
-			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
-
+	///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 cylinder<T>
+	{
+		public:
+		unsigned v[2];		//unique id's designating the starting and ending
+		// default constructor
+		 edge() : cylinder<T>(){v[1] = -1; v[0] = -1;}
+		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
+
+		///@param p is an array of positions in space
+		edge(std::vector< stim::vec<T> > p) : cylinder<T>(p){}
+
+		/// Copy constructor creates an edge from a fiber
+		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
+
+		/// Resamples an edge by calling the fiber resampling function
+		edge resample(T spacing){
+			edge e(cylinder<T>::resample(spacing));	//call the fiber->edge constructor
+			e.v[0] = v[0];					//copy the vertex data
+			e.v[1] = v[1];
+
+			return e;						//return the new edge
 		}
 
+		/// 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<<"("<<cylinder<T>::size()<<")\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();
-
-
-
-		}
 	};
 
-
-//---------------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++;
-	}
+	std::vector<edge> E;       //list of edges
+	std::vector<vertex> V;	    //list of vertices.
 
 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;
-		}
+	///Returns the number of edges in the network.
+	unsigned int edges(){
+		return E.size();
+	}
 
-		return ss.str();
+	///Returns the number of nodes in the network.
+	unsigned int vertices(){
+		return V.size();
 	}
 
-	/// Load a network from an OBJ object
+	stim::cylinder<T> get_cylinder(unsigned f){
+		return E[f];					//return the specified edge (casting it to a fiber)
+	}
 
-	/// @param object is the object file to be used as the basis for the network
-	void load( stim::obj<T> object){
+	//load a network from an OBJ file
+	void load_obj(std::string filename){
 
-		//get the number of vertices in the object
-		unsigned int nV = object.numV();
+		stim::obj<T> O;			//create an OBJ object
+		O.load(filename);		//load the OBJ file as an object
 
-		//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 new_edge = c;							//create an edge from the given centerline
+			unsigned int I = new_edge.size();			//calculate the number of points on the 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 new_vertex = new_edge[0];				//create a new vertex, assign it a position
+				new_vertex.e[0].push_back(E.size());			//add the current edge as outgoing
+				new_edge.v[0] = V.size();						//add the new vertex to the edge
+				V.push_back(new_vertex);						//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
+				new_edge.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 new_vertex = new_edge[I-1];							//create a new vertex, assign it a position
+				new_vertex.e[1].push_back(E.size());						//add the current edge as incoming
+				new_edge.v[1] = V.size();
+				V.push_back(new_vertex);									//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
+				new_edge.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(new_edge);								//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();
+	/// Output the network as a string
+	std::string str(){
 
-			vi++;										//increment the array index
+		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 the resulting array
-		return result;
-	}
+		ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
+		for(unsigned e = 0; e < E.size(); e++){
+			ss<<"\t"<<e<<E[e].str()<<std::endl;
+		}
 
-	std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
-		return f->geometry();
+		return ss.str();
 	}
 
-	/// Generate an OBJ file from the network
-
-	stim::obj<T> obj(){
-
-		//create an OBJ object
-		stim::obj<T> object;
-
-		//name the nodes
-		set_names();
+	/// This function resamples all fibers in a network given a desired minimum spacing
+	/// @param spacing is the minimum distance between two points on the network
+	stim::network<T> resample(T spacing){
+		stim::network<T> n;					//create a new network that will be an exact copy, with resampled fibers
+		n.V = V;							//copy all vertices
 
-		//retrieve a list of terminal node positions
-		std::vector< stim::vec<T> > node_pos = get_node_positions();
+		n.E.resize(edges());				//allocate space for the edge list
 
-		//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++){
+		//copy all fibers, resampling them in the process
+		for(unsigned e = 0; e < edges(); e++){				//for each edge in the edge list
+			n.E[e] = E[e].resample(spacing);				//resample the edge and copy it to the new network
+		}
 
-			//get an array of fiber points
-			std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
+		return n;							              //return the resampled network
+	}
 
-			//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());
+	/// Calculate the total number of points on all edges.
+	unsigned total_points(){
+		unsigned n = 0;
+		for(unsigned e = 0; e < E.size(); e++)
+			n += E[e].size();
+		return n;
+	}
 
-			//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;
+	// gaussian function
+	float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
 
-			//add the intermediate vertex indices to the line array
-			for(unsigned int i = 0; i < fiber_subset.size(); i++){
-				line[1 + i] = nP + i;
-			}
+    // stim 3d vector to annpoint of 3 dimensions
+	void stim2ann(ANNpoint &a, stim::vec<T> b){
+		a[0] = b[0];
+		a[1] = b[1];
+		a[2] = b[2];
+	}
 
-			//add the line list to the object class
-			object.addLine(line);
+	/// Calculate the average magnitude across the entire network.
+	/// @param m is the magnitude value to use. The default is 0 (usually radius).
+	T average(unsigned m = 0){
 
+		T M, L;												//allocate space for the total magnitude and length
+		M = L = 0;											//initialize both the initial magnitude and length to zero
+		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+			M += E[e].integrate(m);							//get the integrated magnitude
+			L += E[e].length();								//get the edge length
 		}
 
-		return object;
+		return M / L;					//return the average magnitude
 	}
 
-	/// 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){
+	/// This function compares two networks and returns the percentage of the current network that is missing from A.
 
-		//determine the number of fibers in the network
-		unsigned int nF = F.size();
+	/// @param A is the network to compare to - the field is generated for A
+	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
+	stim::network<T> compare(stim::network<T> A, float sigma){
 
-		//allocate the necessary space to store the fiber information
-		n0.resize(nF);
-		n1.resize(nF);
-		length.resize(nF);
-		radius.resize(nF);
+		stim::network<T> R;								//generate a network storing the result of the comparison
+		R = (*this);									//initialize the result with the current network
 
-		//assign names (identifiers) to the network components
-		set_names();
+		//generate a KD-tree for network A
+		float metric = 0.0;                               // initialize metric to be returned after comparing the networks
+		ANNkd_tree* kdt;                                 // initialize a pointer to a kd tree
+		double **c;						                 // centerline (array of double pointers) - points on kdtree must be double
+		unsigned int n_data = A.total_points();          // 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 of 3 dimensions
+			c[i] = (double*) malloc(sizeof(double) * 3);
 
-		//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;
+		unsigned t = 0;
+		for(unsigned e = 0; e < A.E.size(); e++){					//for each edge in the network
+			for(unsigned p = 0; p < A.E[e].size(); p++){			//for each point in the edge
+				for(unsigned d = 0; d < 3; d++){				//for each coordinate
 
-			r = f->radius(l);		//get the length and radius of the capillary (calculated at the same time)
+					c[t][d] = A.E[e][p][d];
+				}
+				t++;
+			}
+		}
 
-			radius[i] = r;			//store the radius in the output array
-			length[i] = l;			//store the length in the output array
+		//compare each point in the current network to the field produced by A
+		ANNpointArray pts = (ANNpointArray)c;           // create an array of data points of type double
+		kdt = new ANNkd_tree(pts, n_data, 3);			// build a KD tree using the annpointarray
+		double eps = 0; // error bound
+		ANNdistArray dists = new ANNdist[1];     // near neighbor distances
+		ANNidxArray nnIdx = new ANNidx[1];				// near neighbor indices // allocate near neigh indices
+
+		stim::vec<T> p0, p1;
+		float m0, m1;
+		float M = 0;											//stores the total metric value
+		float l;												//stores the segment length
+		float L = 0;											//stores the total network length
+		ANNpoint queryPt = annAllocPt(3);
+		for(unsigned e = 0; e < R.E.size(); e++){					//for each edge in A
+			R.E[e].add_mag(0);									//add a new magnitude for the metric
+
+			for(unsigned p = 0; p < R.E[e].size(); p++){			//for each point in the edge
+
+				p1 = R.E[e][p];									//get the next point in the edge
+				stim2ann(queryPt, p1);
+				kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);	//find the distance between A and the current network
+				m1 = 1.0f - gaussianFunction(dists[0], sigma);		//calculate the metric value based on the distance
+				R.E[e].set_mag(m1, p, 1);						//set the error for the second point in the segment
 
-			i++;					//increment the array index
+			}
 		}
 
-
+		return R;		//return the resulting network
 	}
 
-};
-
-};	//end namespace stim
+	/// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
+	unsigned nmags(){
+		return E[0].nmags();
+	}
 
 
+};		//end stim::network class
+};		//end stim namespace
 #endif
diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h
index 8fe05bd..665e5ea 100644
--- a/stim/gl/gl_spider.h
+++ b/stim/gl/gl_spider.h
@@ -96,7 +96,7 @@ class gl_spider : public virtual gl_texture<T>
 		std::vector< stim::vec<float> > cD;	//Direction of line currently being traced.
 		std::vector< stim::vec<float> > cM;	//Magnitude of line currently being traced.
 
-		stim::glObj<float> sk;			//object to store the skeleton.
+//		stim::glObj<float> sk;			//object to store the skeleton.
 		stim::glnetwork<float> nt;		//object for storing the network.
 
 		stim::vec<float> rev;			//reverse vector;
@@ -1157,7 +1157,7 @@ class gl_spider : public virtual gl_texture<T>
                    		{
                            		float x, y, z, u, v, w, m;
                            		myfile >> x >> y >> z >> u >> v >> w >> m;
-					setSeed(x, y , z);
+					setSeed(x, y, z);
 					setSeedVec(u, v, w);
 					setSeedMag(m);
                    		}
@@ -1171,14 +1171,28 @@ class gl_spider : public virtual gl_texture<T>
 		void
 		saveNetwork(std::string name)
 		{
+			stim::glObj<float> sk;
+			for(int i = 0; i < nt.sizeE(); i++)
+			{
+				std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
+                 		std::vector<stim::vec< float > > ce = nt.getEdgeCenterLine(i);
+				sk.Begin(stim::OBJ_LINE);
+				for(int j = 0; j < ce.size(); j++)
+				{
+					sk.TexCoord(cm[j][0]);
+					sk.Vertex(ce[j][0], ce[j][1], ce[j][2]);
+				}
+				sk.End();
+			}	
 			sk.save(name);
 		}
 
+		///Depreciated, but might be reused later()
 		///returns a COPY of the entire stim::glObj object.
 		stim::glObj<float>
 		getNetwork()
 		{
-			return sk;
+//			return sk;
 		}
 
 		///returns a COPY of the entire stim::glnetwork object.
@@ -1326,114 +1340,35 @@ class gl_spider : public virtual gl_texture<T>
 		void
 		trace(int min_cost)
 		{	
-			Bind();
-			rev = stim::vec<float>(0.0,0.0,1.0);
+//			rev = stim::vec<float>(0.0,0.0,1.0);
 			bool sEmpty = true;
 			float lastmag = 16.0;;
-			while(!seeds.empty())
+			stim::vec<float> curSeed; 
+			stim::vec<float> curSeedVec;			
+			float curSeedMag;
+			while(!Empty())
 			{
 				//clear the currently traced line and start a new one.
 				cL.clear();
 				cM.clear();
-				sk.Begin(stim::OBJ_LINE);
-				stim::vec<float> curSeed = seeds.top();
-//				std::cout << "The current seeds is " << curSeed << std::endl;
-				stim::vec<float> curSeedVec = seedsvecs.top();
-				float curSeedMag = seedsmags.top();
+				cD.clear();
+				curSeed = seeds.top();
+				curSeedVec = seedsvecs.top();
+				curSeedMag = seedsmags.top();
 				seeds.pop();
 				seedsvecs.pop();
 				seedsmags.pop();
 //				std::cout << "The current seed Vector is " << curSeedVec << std::endl;
 				setPosition(curSeed);
 				setDirection(curSeedVec);
-				cL.push_back(curSeed);
-				cM.push_back(curSeedMag);
-				sk.createFromSelf(GL_SELECT);
-				traceLine(min_cost);
-
-				sk.rev();
-		//		std::cout << "reversed" << std::endl;
-				std::reverse(cL.begin(), cL.end());
-				std::reverse(cM.begin(), cM.end());
-				setPosition(curSeed);
-				setDirection(-rev);
-				setMagnitude(16.0);
-				sk.createFromSelf(GL_SELECT);
-				traceLine(min_cost);
-				sk.End();
+				setMagnitude(curSeedMag);
+//				cL.push_back(curSeed);
+//				cM.push_back(curSeedMag);
+//				cD.push_back(curSeedMag);
+				pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost);
 			}
-			Unbind();
 		}
 
-		///@param min_cost the cost value used for tracing
-		///traces the seedpoint passed to completion in one directions.
-		void
-		traceLine(int min_cost)
-		{
-			stim::vec<float> pos;
-			stim::vec<float> mag;
-			int h;
-			bool started = false;
-			bool running = true;
-			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
-			while(running)
-			{
-				int cost = Step();
-				if (cost > min_cost){
-					running = false;
-					break;
-				} else {
-					//Have we found an edge?
-					pos = getPosition();
-					if(pos[0] > size[0] || pos[1] > size[1]
-					 || pos[2] > size[2] || pos[0] < 0
-					 || pos[1] < 0 || pos[2] < 0)
-					{
-//					       std::cout << "Found Edge" << std::endl;
-						running = false;
-						break;
-					}
-					//If this is the first step in the trace,
-					// save the direction
-					//(to be used later to trace the fiber in the opposite direction)
-					if(started == false){
-						rev = -getDirection();
-						started = true;
-					}
-//					std::cout << i << p << std::endl;
-					m = getMagnitude();
-					//Has the template size gotten unreasonable?
-					if(m[0] > 75 || m[0] < 1){
-//						std::cout << "Magnitude Limit" << std::endl;
-						running = false;
-						break;
-					}
-					else
-					{
-						h = selectObject(pos, getDirection(), m[0]);
-						//Have we hit something previously traced?
-						if(h != -1){
-							std::cout << "I hit a line" << h << std::endl;
-							running = false;
-							break;
-						}
-						else {          
-							cL.push_back(stim::vec<float>(p[0], p[1],p[2]));//
-							sk.TexCoord(m[0]);
-							sk.Vertex(p[0], p[1], p[2]);
-							Bind(btexbufferID, bfboID, 27);
-							CHECK_OPENGL_ERROR
-							branchDetection();
-							CHECK_OPENGL_ERROR
-							Unbind();
-							CHECK_OPENGL_ERROR
-						}
-				 	}
-                 		}
-         		}
-		}	
-
-
 		int
 		selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag) 
 		{
@@ -1597,7 +1532,7 @@ class gl_spider : public virtual gl_texture<T>
 			stim::vec<float> sdir = getDirection();	
 
 			Bind();
-			sk.Begin(stim::OBJ_LINE);
+//			sk.Begin(stim::OBJ_LINE);
 
 
 //			sk.createFromSelf(GL_SELECT);
@@ -1618,7 +1553,7 @@ class gl_spider : public virtual gl_texture<T>
 				if (cost > min_cost){
 					std::cout << "Cost Limit" << std::endl;
 					running = false;
-					sk.End();
+//					sk.End();
 		branchDetection2();
 					pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
 					addToNetwork(a, spos, smag, sdir);
@@ -1634,7 +1569,7 @@ class gl_spider : public virtual gl_texture<T>
 						std::cout << "Edge Limit" << std::endl;
 //					       std::cout << "Found Edge" << std::endl;
 						running = false;
-						sk.End();
+//						sk.End();
 		branchDetection2();
 						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
 						addToNetwork(a, spos, smag, sdir);
@@ -1654,7 +1589,7 @@ class gl_spider : public virtual gl_texture<T>
 					if(mag[0] > 75 || mag[0] < 1){
 						std::cout << "Magnitude Limit" << std::endl;
 						running = false;
-						sk.End();
+//						sk.End();
 		branchDetection2();
 						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
 						addToNetwork(a, spos, smag, sdir);
@@ -1668,7 +1603,7 @@ class gl_spider : public virtual gl_texture<T>
 						if(h != -1){
 							std::cout << "Hit Limit" << std::endl;
 							running = false;
-							sk.End();
+//							sk.End();
 		branchDetection2();
 							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
 							addToNetwork(a, spos, smag, sdir);
@@ -1680,8 +1615,8 @@ class gl_spider : public virtual gl_texture<T>
 							cM.push_back(stim::vec<float>(m[0], m[0]));
 //							cM.push_back(m[0]);
 
-							sk.TexCoord(m[0]);
-							sk.Vertex(p[0], p[1], p[2]);
+//							sk.TexCoord(m[0]);
+//							sk.Vertex(p[0], p[1], p[2]);
 							Bind(btexbufferID, bfboID, 27);
 							CHECK_OPENGL_ERROR
 //							branchDetection();
diff --git a/stim/visualization/aaboundingbox.h b/stim/visualization/aaboundingbox.h
new file mode 100644
index 0000000..a0f54cd
--- /dev/null
+++ b/stim/visualization/aaboundingbox.h
@@ -0,0 +1,72 @@
+#ifndef STIM_AABB
+#define STIM_AABB
+
+namespace stim{
+
+
+/// This class describes a structure for an axis-aligned bounding box
+template< typename T >
+class aaboundingbox{
+
+public:
+	bool set;				//has the bounding box been set to include any points?
+	stim::vec<T> A;			//minimum point in the bounding box
+	stim::vec<T> B;			//maximum point in the bounding box
+
+	aaboundingbox(){					//constructor generates an empty bounding box
+		set = false;
+	}
+
+
+	/// Test if a point is inside of the bounding box and returns true if it is.
+
+	/// @param p is the point to be tested
+	bool test(stim::vec<T> p){
+
+		for(unsigned d = 0; d < p.size(); p++){		//for each dimension
+			if(p[d] < A[d]) return false;			//if the point is less than the minimum bound, return false
+			if(p[d] > B[d]) return false;			//if the point is greater than the max bound, return false
+		}
+		return true;
+	}
+
+	/// Expand the bounding box to include the specified point.
+
+	/// @param p is the point to be included
+	void expand(stim::vec<T> p){
+
+		if(!set){							//if the bounding box is empty, fill it with the current point
+			A = B = p;
+			set = true;
+		}
+
+		for(unsigned d = 0; d < p.size(); d++){		//for each dimension
+			if(p[d] < A[d]) A[d] = p[d];			//expand the bounding box as necessary
+			if(p[d] > B[d]) B[d] = p[d];
+		}
+	}
+
+	/// Return the center point of the bounding box as a stim::vec
+	stim::vec<T> center(){
+		return (B + A) * 0.5;
+	}
+
+	/// Return the size of the bounding box as a stim::vec
+	stim::vec<T> size(){
+		return (B - A);
+	}
+
+	/// Generate a string for the bounding box
+	std::string str(){
+		std::stringstream ss;
+		ss<<A.str()<<"----->"<<B.str();
+		return ss.str();
+	}
+
+
+};		//end stim::aabb
+
+
+};		//end namespace stim
+
+#endif
\ No newline at end of file
diff --git a/stim/visualization/gl_aaboundingbox.h b/stim/visualization/gl_aaboundingbox.h
new file mode 100644
index 0000000..de81ae8
--- /dev/null
+++ b/stim/visualization/gl_aaboundingbox.h
@@ -0,0 +1,60 @@
+#ifndef STIM_GL_AABB
+#define STIM_GL_AABB
+
+#include <stim/visualization/aaboundingbox.h>
+#include <GL/gl.h>
+
+namespace stim{
+
+template <typename T>
+class gl_aaboundingbox : public aaboundingbox<T>{
+
+public:
+
+	//default constructor
+	gl_aaboundingbox() : stim::aaboundingbox<T>(){}
+
+	//constructor takes an AABB
+	gl_aaboundingbox(stim::aaboundingbox<T> b) : stim::aaboundingbox<T>(b){}
+
+
+	/// Specifies vertices of the bounding box using CW winding. Use GL_LINE_LOOP for wireframe or GL_QUADS for a solid.
+	void glWire(){
+
+		//front plane (in A[2])
+		glBegin(GL_LINE_LOOP);
+			glVertex3f(A[0], A[1], A[2]);
+			glVertex3f(A[0], B[1], A[2]);
+			glVertex3f(B[0], B[1], A[2]);
+			glVertex3f(B[0], A[1], A[2]);
+		glEnd();
+
+		//back plane (in B[2])
+		glBegin(GL_LINE_LOOP);
+			glVertex3f(B[0], B[1], B[2]);
+			glVertex3f(A[0], B[1], B[2]);
+			glVertex3f(A[0], A[1], B[2]);			
+			glVertex3f(B[0], A[1], B[2]);
+		glEnd();
+
+		//fill out the rest of the lines to connect the two faces
+		glBegin(GL_LINES);
+			glVertex3f(A[0], B[1], A[2]);
+			glVertex3f(A[0], B[1], B[2]);
+			glVertex3f(B[0], B[1], B[2]);
+			glVertex3f(B[0], B[1], A[2]);
+			glVertex3f(B[0], A[1], A[2]);
+			glVertex3f(B[0], A[1], B[2]);
+			glVertex3f(A[0], A[1], B[2]);
+			glVertex3f(A[0], A[1], A[2]);
+		glEnd();
+
+	}
+
+
+};		//end stim::gl_aabb
+
+
+};		//end namespace stim
+
+#endif
\ No newline at end of file
--
libgit2 0.21.4