add many function in order to work with netmets

Jiaming Guo
1 parent 57914be5
Showing 10 changed files with 1849 additions and 397 deletions Show diff stats
stim/biomodels/centerline.h
stim/biomodels/centerline_dep.h
stim/biomodels/network.h
stim/math/circle.h
stim/math/circle_dep.h
stim/math/filters/conv2.h
stim/visualization/aabbn.h
stim/visualization/cylinder.h
stim/visualization/cylinder_dep.h
stim/visualization/gl_network.h
@@ -3,7 +3,6 @@
 #include <vector>
 #include <stim/math/vec3.h>
-//#include <ANN/ANN.h>
 namespace stim{
@@ -12,195 +11,123 @@ namespace stim{
  *	class to describe an interconnected (often biological) network.
  */
 template<typename T>
-class centerline{
+class centerline : public std::vector< stim::vec3<T> >{
 protected:
-	unsigned int N;					//number of points in the fiber
-	double **c;						//centerline (array of double pointers)
-//	ANNkd_tree* kdt;				//kd-tree stores all points in the fiber for fast searching
-	/// Initialize an empty fiber
-	void init()
-	{
-		N=0;
-		c=NULL;
-//		kdt = NULL;
-	}
-
-	/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
-	void init(unsigned int n)
-	{
+	std::vector<T> L;										//stores the integrated length along the fiber (used for parameterization)
-		N = n;												//set the number of points
-//		kdt = NULL;
-		c = (double**) malloc(sizeof(double*) * N);			//allocate the array pointer
+	///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
+	vec3<T> d(size_t i) {
+		if (size() <= 1) return vec3<T>(0, 0, 0);						//if there is insufficient information to calculate the direction, return a null vector
+		if (size() == 2) return (at(1) - at(0)).norm();					//if there are only two points, the direction vector at both is the direction of the line segment
+		if (i == 0) return (at(1) - at(0)).norm();						//the first direction vector is oriented towards the first line segment
+		if (i == size() - 1) return at(size() - 1) - at(size() - 2);	//the last direction vector is oriented towards the last line segment
-		for(unsigned int i = 0; i < N; i++)					//allocate space for each point
-			c[i] = (double*) malloc(sizeof(double) * 3);
+		//all other direction vectors are the average direction of the two joined line segments
+		return ((at(i) - at(i - 1)).norm() + (at(i + 1) - at(i)).norm()).norm();
 	}
-
-	/// Copies an existing fiber to the current fiber
-
-	/// @param cpy stores the new copy of the fiber
-	void copy( const stim::centerline<T>& cpy, bool kd = 0){
-
-		///allocate space for the new fiber
-		init(cpy.N);
-
-		///copy the points
-		for(unsigned int i = 0; i < N; i++){
-			for(unsigned int d = 0; d < 3; d++)		//for each dimension
-				c[i][d] = cpy.c[i][d];				//copy the coordinate
+	//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
+	void update_L(size_t start = 0) {
+		L.resize(size());									//allocate space for the L array
+		if (start == 0) {
+			L[0] = 0;											//initialize the length value for the first point to zero (0)
+			start++;
 		}
-//		if(kd)
-//			gen_kdtree();							//generate the kd tree for the new fiber
-	}
-	/// generate a KD tree for points on fiber
-//	void gen_kdtree()
-//	{
-//		int n_data = N; //create an array of data points
-//		ANNpointArray pts = (ANNpointArray)c;			//cast the centerline list to an ANNpointArray
-//		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
-//	}
-
-	/// find distance between two points
-	double dist(double* p0, double* p1){
-
-		double sum = 0; // initialize variables
-		float v;
-		for(unsigned int d = 0; d < 3; d++)
-		{
-			v = p1[d] - p0[d];
-			sum +=v * v;
+		stim::vec3<T> d;
+		for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
+			d = at(i) - at(i - 1);
+			L[i] = L[i - 1] + d.len();				//calculate the running length total
 		}
-		return sqrt(sum);
 	}
-	/// This function retreives the index for the fiber point closest to q
-
-	/// @param q is a reference point used to find the closest point on the fiber center line
-//	unsigned int ann( stim::vec<double> q ){
-//		ANNidxArray idx = new ANNidx[1];			//variable used to hold the nearest point
-//		ANNdistArray sq_dist = new ANNdist[1];		//variable used to hold the squared distance to the nearest point
-//		kdt->annkSearch(q.data(), 1, idx, sq_dist);	//search the KD tree for the nearest neighbor
-//		return *idx;
-//	}
+	void init() {
+		if (size() == 0) return;								//return if there aren't any points
+		update_L();
+	}
 	/// 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::vec3<T> r;
-		r.resize(3);
-		r[0] = c[i][0];
-		r[1] = c[i][1];
-		r[2] = c[i][2];
+		return std::vector< stim::vec3<T> >::at(i);
+	}
+
+	///finds the index of the point closest to the length l on the lower bound.
+	///binary search.
+	size_t findIdx(T l) {
+		size_t i = L.size() / 2;
+		size_t max = L.size() - 1;
+		size_t min = 0;
+		while (i > 0 && i < L.size() - 1){
+			if (l < L[i]) {
+				max = i;
+				i = min + (max - min) / 2;
+			}
+			else if (L[i] <= l && L[i + 1] >= l) {
+				break;
+			}
+			else {
+				min = i;
+				i = min + (max - min) / 2;
+			}
+		}
+		return i;
+	}
-		return r;
+	///Returns a position vector at the given length into the fiber (based on the pvalue).
+	///Interpolates the radius along the line.
+	///@param l: the location of the in the cylinder.
+	///@param idx: integer location of the point closest to l but prior to it.
+	stim::vec3<T> p(T l, int idx) {
+		T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
+		stim::vec3<T> v1 = at(idx);
+		stim::vec3<T> v2 = at(idx + 1);
+		return(v1 + (v2 - v1)*rat);
 	}
 public:
-	centerline(){
-		init();
-	}
+	using std::vector< stim::vec3<T> >::at;
+	using std::vector< stim::vec3<T> >::size;
-	/// Copy constructor
-	centerline(const stim::centerline<T> &obj){
-		copy(obj);
+	centerline() : std::vector< stim::vec3<T> >() {
+		init();
 	}
-
-	//temp constructor for graph visualization
-	centerline(int n)
-	{
-		init(n);
+	centerline(size_t n) : std::vector< stim::vec3<T> >(n){
+		init();
 	}
-	/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
-	centerline(std::vector< stim::vec<T> > p, bool kd = 0){
-		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
-		}
-		//generate a kd tree
-//		if(kd)
-//			gen_kdtree();
+	//overload the push_back function to update the length vector
+	void push_back(stim::vec3<T> p) {
+		std::vector< stim::vec3<T> >::push_back(p);
+		update_L(size() - 1);
 	}
-	/// constructor takes a list of points
-	centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){
-		init(pos.size());		//initialize the fiber
+	///Returns a position vector at the given p-value (p value ranges from 0 to 1).
+	///interpolates the position along the line.
+	///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
+	stim::vec3<T> p(T pvalue) {
+		if (pvalue <= 0.0) return at(0);			//return the first element
+		if (pvalue >= 1.0) return back();			//return the last element
-		//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) pos[i][d];
-			//set the radius
-		}
-
-		//generate a kd tree
-		//if(kd)
-		//	gen_kdtree();
+		T l = pvalue*L[L.size() - 1];
+		int idx = findIdx(l);
+		return p(l, idx);
 	}
-	/// Assignment operation
-	centerline& operator=(const centerline &rhs){
-		if(this == &rhs) return *this;			//test for and handle self-assignment
-		copy(rhs);
-		return *this;
-	}
-
-
-	/// Return the point on the fiber closest to q
-	/// @param q is the query point used to locate the nearest point on the fiber centerline
-//	stim::vec<T> nearest(stim::vec<T> q){
-//
-//		stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]);
-//
-//		unsigned int idx = ann(temp);		//determine the index of the nearest neighbor
-//
-//		return stim::vec<T>((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]);	//return the nearest centerline point
-//	}
-
-	/// Return the point index on the fiber closest to q
-	/// @param q is the query point used to locate the nearest point on the fiber centerline
-//	unsigned int nearest_idx(stim::vec<T> q){
-//
-//		stim::vec<double> temp((double) q[0], (double) q[1], (double) q[2]);
-//
-//		unsigned int idx = ann(temp);		//determine the index of the nearest neighbor
-//
-//		return idx;	//return the nearest centerline point index
-//	}
-
-	/// Returns the fiber centerline as an array of stim::vec points
-	std::vector< stim::vec<T> > get_centerline(){
-
-		//create an array of stim vectors
-		std::vector< stim::vec3<T> > pts(N);
-
-		//cast each point to a stim::vec, keeping only the position information
-		for(unsigned int i = 0; i < N; i++)
-			pts[i] = stim::vec3<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
-
-		//return the centerline array
-		return pts;
+	///Return the length of the entire centerline
+	T length() {
+		return L.back();
 	}
 	/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
 	std::vector< stim::centerline<T> > split(unsigned int idx){
-		std::vector< stim::centerline<T> > fl;		//create an array to store up to two fibers
+		std::vector< stim::centerline<T> > fl;				//create an array to store up to two fibers
+		size_t N = size();
 		//if the index is an end point, only the existing fiber is returned
 		if(idx == 0 || idx == N-1){
@@ -211,128 +138,89 @@ public:
 		//if the index is not an end point
 		else{
-			unsigned int N1 = idx + 1;					//calculate the size of both fibers
+			unsigned int N1 = idx;					//calculate the size of both fibers
 			unsigned int N2 = N - idx;
 			fl.resize(2);								//set the array size to 2
-			fl[0].init(N1);								//set the size of each fiber
-			fl[1].init(N2);
+			fl[0] = stim::centerline<T>(N1);			//set the size of each fiber
+			fl[1] = stim::centerline<T>(N2);
 			//copy both halves of the fiber
-			unsigned int i, d;
+			unsigned int i;
 			//first half
-			for(i = 0; i < N1; i++){					//for each centerline point
-				for(d = 0; d < 3; d++)
-					fl[0].c[i][d] = c[i][d];			//copy each coordinate
-			}
+			for(i = 0; i < N1; i++)					//for each centerline point
+				fl[0][i] = std::vector< stim::vec3<T> >::at(i);
+			fl[0].init();							//initialize the length vector
 			//second half
-			for(i = 0; i < N2; i++){
-				for(d = 0; d < 3; d++)
-					fl[1].c[i][d] = c[idx + i][d];
-
-			}
+			for(i = 0; i < N2; i++)
+				fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
+			fl[1].init();							//initialize the length vector
 		}
 		return fl;										//return the array
 	}
-	/// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f
-
-	/// @param f is the fiber that will be connected to the current fiber
-/*	std::vector< stim::centerline<T> > connect( stim::centerline<T> &f, double dist){
-
-		double min_dist;
-		unsigned int idx0, idx1;
-
-		//go through each point in the query fiber, looking for the indices for the closest points
-		for(unsigned int i = 0; i < f.n_pts(); i++){
-			//Run through all points and find the index with the closest point, then partition the fiber and return two fibers.
-
-		}
-
-
-
-	}
-*/
 	/// Outputs the fiber as a string
 	std::string str(){
 		std::stringstream ss;
-
-		//create an iterator for the point list
-		//typename std::list< point<T> >::iterator i;
-		for(unsigned int i = 0; i < N; i++){
-			ss<<"  [  ";
-			for(unsigned int d = 0; d < 3; d++){
-				ss<<c[i][d]<<"  ";
-			}
-		}
+		size_t N = std::vector< stim::vec3<T> >::size();
+		ss << "---------[" << N << "]---------" << std::endl;
+		for (size_t i = 0; i < N; i++)
+			ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
+		ss << "--------------------" << std::endl;
 		return ss.str();
 	}
-	/// Returns the number of centerline points in the fiber
-	unsigned int size(){
-		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);
+	stim::vec3<T> back(){
+		return std::vector< stim::vec3<T> >::back();
 	}
+
 		////resample a fiber in the network
 	stim::centerline<T> resample(T spacing)
 	{
 		std::cout<<"fiber::resample()"<<std::endl;
-		std::vector<T> v(3);    //v-direction vector of the segment
-		stim::vec<T> p(3);      //- intermediate point to be added
-		stim::vec<T> p1(3);   // p1 - starting point of an segment on the fiber,
-		stim::vec<T> p2(3);   // p2 - ending point,
-		double sum=0;  //distance summation
-		std::vector<stim::vec<T> > fiberPositions = centerline();
-		std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
+		stim::vec3<T> v;    //v-direction vector of the segment
+		stim::vec3<T> p;      //- intermediate point to be added
+		stim::vec3<T> p1;   // p1 - starting point of an segment on the fiber,
+		stim::vec3<T> p2;   // p2 - ending point,
+		//double sum=0;  //distance summation
+
+		size_t N = size();
+
+		centerline<T> new_c; // initialize list of new resampled points on the fiber
 		// for each point on the centerline (skip if it is the last point on centerline)
-		//unsigned int N = fiberPositions.size(); // number of points on the fiber
 		for(unsigned int f=0; f< N-1; f++)
-		{
+		{			
+			p1 = at(f); 
+			p2 = at(f+1);
+			v = p2 - p1;
-			p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
-			for(unsigned int d = 0; d < 3; d++){
-				sum +=v[d] * v[d];}              //length of segment-distance between starting and ending point
-
-			T lengthSegment = sqrt(sum);  //find Length of the segment as distance between the starting and ending points of the segment
-
-			if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
-				{
-					// repeat resampling until accumulated stepsize is equsl to length of the segment
-					for(T step=0.0; step<lengthSegment; step+=spacing)
-					{
-						// calculate the resampled point by travelling step size in the direction of normalized gradient vector
-						for(unsigned int i=0; i<3;i++)
-							{
-								p[i] = p1[i] + v[i]*(step/lengthSegment);
-							} //for each dimension
-						// add this resampled points to the new fiber list
-						newPointList.push_back(p);
-					}
+			T lengthSegment = v.len();			//find Length of the segment as distance between the starting and ending points of the segment
+
+			if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample
+				
+				// repeat resampling until accumulated stepsize is equsl to length of the segment
+				for(T step=0.0; step<lengthSegment; step+=spacing){
+					// calculate the resampled point by travelling step size in the direction of normalized gradient vector
+					p = p1 + v * (step / lengthSegment);
+					
+					// add this resampled points to the new fiber list
+					new_c.push_back(p);
 				}
-			else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
-				newPointList.push_back(fiberPositions[f+1]);
 			}
-		newPointList.push_back(fiberPositions[N-1]);   //add the last point on the fiber to the new fiber list
-		centerline newFiber(newPointList);
-		return newFiber;
+			else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
+				new_c.push_back(at(f+1));
+		}
+		new_c.push_back(at(N-1));   //add the last point on the fiber to the new fiber list
+		//centerline newFiber(newPointList);
+		return new_c;
 	}
 };
+#ifndef STIM_CENTERLINE_H
+#define STIM_CENTERLINE_H
+
+#include <vector>
+#include <stim/math/vec3.h>
+
+namespace stim{
+
+/**	This class stores information about a single fiber represented as a set of geometric points
+ *	between two branch or end points. This class is used as a fundamental component of the stim::network
+ *	class to describe an interconnected (often biological) network.
+ */
+template<typename T>
+class centerline{
+
+protected:
+	unsigned int N;					//number of points in the fiber
+	double **c;						//centerline (array of double pointers)
+
+	/// Initialize an empty fiber
+	void init(){
+		N=0;
+		c=NULL;
+	}
+
+	/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
+	void init(unsigned int n){
+
+		N = n;												//set the number of points
+		c = (double**) malloc(sizeof(double*) * N);			//allocate the array pointer
+
+		for(unsigned int i = 0; i < N; i++)					//allocate space for each point
+			c[i] = (double*) malloc(sizeof(double) * 3);
+	}
+
+	/// Copies an existing fiber to the current fiber
+
+	/// @param cpy stores the new copy of the fiber
+	void copy( const stim::centerline<T>& cpy, bool kd = 0){
+
+		///allocate space for the new fiber
+		init(cpy.N);
+
+		///copy the points
+		for(unsigned int i = 0; i < N; i++){
+			for(unsigned int d = 0; d < 3; d++)		//for each dimension
+				c[i][d] = cpy.c[i][d];				//copy the coordinate
+		}
+	}
+
+	/// find distance between two points
+	double dist(double* p0, double* p1){
+
+		double sum = 0; // initialize variables
+		float v;
+		for(unsigned int d = 0; d < 3; d++)
+		{
+			v = p1[d] - p0[d];
+			sum +=v * v;
+		}
+		return sqrt(sum);
+	}
+
+	/// 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::vec3<T> r;
+		r.resize(3);
+		r[0] = c[i][0];
+		r[1] = c[i][1];
+		r[2] = c[i][2];
+
+		return r;
+	}
+
+
+public:
+
+	centerline(){
+		init();
+	}
+
+	/// Copy constructor
+	centerline(const stim::centerline<T> &obj){
+		copy(obj);
+	}
+
+	//initialize a centerline with n points
+	centerline(int n){
+		init(n);
+	}
+
+	/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
+	centerline(std::vector< stim::vec<T> > p, bool kd = 0){
+		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];
+		}
+	}
+
+	/// constructor takes a list of points
+	centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){
+		init(pos.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) pos[i][d];
+		}
+	}
+
+	/// Assignment operation
+	centerline& operator=(const centerline &rhs){
+		if(this == &rhs) return *this;			//test for and handle self-assignment
+		copy(rhs);
+		return *this;
+	}
+
+
+	/// Returns the fiber centerline as an array of stim::vec points
+	std::vector< stim::vec<T> > get_centerline(){
+
+		//create an array of stim vectors
+		std::vector< stim::vec3<T> > pts(N);
+
+		//cast each point to a stim::vec, keeping only the position information
+		for(unsigned int i = 0; i < N; i++)
+			pts[i] = stim::vec3<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
+
+		//return the centerline array
+		return pts;
+	}
+
+	/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
+	std::vector< stim::centerline<T> > split(unsigned int idx){
+
+		std::vector< stim::centerline<T> > fl;		//create an array to store up to two fibers
+
+		//if the index is an end point, only the existing fiber is returned
+		if(idx == 0 || idx == N-1){
+			fl.resize(1);							//set the size of the fiber to 1
+			fl[0] = *this;							//copy the current fiber
+		}
+
+		//if the index is not an end point
+		else{
+
+			unsigned int N1 = idx + 1;					//calculate the size of both fibers
+			unsigned int N2 = N - idx;
+
+			fl.resize(2);								//set the array size to 2
+
+			fl[0].init(N1);								//set the size of each fiber
+			fl[1].init(N2);
+
+			//copy both halves of the fiber
+			unsigned int i, d;
+
+			//first half
+			for(i = 0; i < N1; i++){					//for each centerline point
+				for(d = 0; d < 3; d++)
+					fl[0].c[i][d] = c[i][d];			//copy each coordinate
+			}
+
+			//second half
+			for(i = 0; i < N2; i++){
+				for(d = 0; d < 3; d++)
+					fl[1].c[i][d] = c[idx + i][d];
+
+			}
+		}
+
+		return fl;										//return the array
+
+	}
+
+	/// Outputs the fiber as a string
+	std::string str(){
+		std::stringstream ss;
+
+		//create an iterator for the point list
+		//typename std::list< point<T> >::iterator i;
+		for(unsigned int i = 0; i < N; i++){
+			ss<<"  [  ";
+			for(unsigned int d = 0; d < 3; d++){
+				ss<<c[i][d]<<"  ";
+			}
+		}
+
+		return ss.str();
+	}
+	/// Returns the number of centerline points in the fiber
+	unsigned int size(){
+		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);
+	}
+		////resample a fiber in the network
+	stim::centerline<T> resample(T spacing)
+	{
+		std::cout<<"fiber::resample()"<<std::endl;
+
+		std::vector<T> v(3);    //v-direction vector of the segment
+		stim::vec<T> p(3);      //- intermediate point to be added
+		stim::vec<T> p1(3);   // p1 - starting point of an segment on the fiber,
+		stim::vec<T> p2(3);   // p2 - ending point,
+		double sum=0;  //distance summation
+		std::vector<stim::vec<T> > fiberPositions = centerline();
+		std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
+		// for each point on the centerline (skip if it is the last point on centerline)
+		for(unsigned int f=0; f< N-1; f++)
+		{
+			
+			p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
+			for(unsigned int d = 0; d < 3; d++){
+				sum +=v[d] * v[d];}              //length of segment-distance between starting and ending point
+
+			T lengthSegment = sqrt(sum);  //find Length of the segment as distance between the starting and ending points of the segment
+
+			if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
+				{
+					// repeat resampling until accumulated stepsize is equsl to length of the segment
+					for(T step=0.0; step<lengthSegment; step+=spacing)
+					{
+						// calculate the resampled point by travelling step size in the direction of normalized gradient vector
+						for(unsigned int i=0; i<3;i++)
+							{
+								p[i] = p1[i] + v[i]*(step/lengthSegment);
+							} //for each dimension
+						// add this resampled points to the new fiber list
+						newPointList.push_back(p);
+					}
+				}
+			else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
+				newPointList.push_back(fiberPositions[f+1]);
+			}
+		newPointList.push_back(fiberPositions[N-1]);   //add the last point on the fiber to the new fiber list
+		centerline newFiber(newPointList);
+		return newFiber;
+	}
+
+};
+
+
+
+}	//end namespace stim
+
+
+
+#endif
@@ -15,6 +15,39 @@
 #include <stim/cuda/cudatools/timer.h>
+#ifdef __CUDACC__
+//device gaussian function
+__device__ float gaussianFunction(float x, float std){ return expf(-x/(2*std*std));}
+
+//compute metric in parallel
+template <typename T>
+__global__ void d_metric(T* M, size_t n, T* D, float sigma){
+	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
+	if(x >= n) return;
+	M[x] = 1.0f - gaussianFunction(D[x], sigma);
+}
+
+//find the corresponding edge index from array index
+__global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){
+	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
+	if(x >= n) return;
+	unsigned i = 0;
+	size_t N = 0;
+	for(unsigned e = 0; e < ne; e++){
+		N += E[e];
+		if(I[x] < N){
+			R[x] = i;
+			break;
+		}
+		i++;
+	}
+}
+#endif
+
+//hard-coded factor
+int threshold_fac = 10;
+float metric_fac = 0.6f;							// might be related to the distance field
+
 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:
@@ -22,7 +55,6 @@ namespace stim{
  *   2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
 */
-
 template<typename T>
 class network{
@@ -33,16 +65,15 @@ class network{
 		public:
 		unsigned v[2];		//unique id's designating the starting and ending
 		// default constructor
-		edge() : cylinder<T>()
-		{
+		edge() : cylinder<T>() {
 			v[1] = (unsigned)(-1); v[0] = (unsigned)(-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::vec3<T> > p) : cylinder<T>(p){}
+		edge(centerline p) : cylinder<T>(p){}
-		/// Copy constructor creates an edge from a fiber
+		/// Copy constructor creates an edge from a cylinder
 		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
 		/// Resamples an edge by calling the fiber resampling function
@@ -61,6 +92,19 @@ class network{
 			return ss.str();
 		}
+		std::vector<edge> split(unsigned int idx){
+		
+			std::vector< stim::cylinder<T> > C;
+			C.resize(2);
+			C =	(*this).cylinder<T>::split(idx);
+			std::vector<edge> E(C.size());
+
+			for(unsigned e = 0; e < E.size(); e++){
+				E[e] = C[e];
+			}
+			return E;
+		}
+
 	};
 	///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.
@@ -291,7 +335,7 @@ public:
 			std::vector< stim::vec<T> > c;						//allocate an array of points for the vessel centerline
 			O.getLine(l, c);							//get the fiber centerline
-			std::vector< stim::vec3<T> > c3(c.size());
+			stim::centerline<T> c3(c.size());
 			for(size_t j = 0; j < c.size(); j++)
 				c3[j] = c[j];
@@ -377,7 +421,8 @@ public:
 		return n;							              	//return the resampled network
 	}
-
+	// host gaussian function
+	__host__ float gaussianFunction(float x, float std = 25){ return exp(-x/(2*std*std));}  // std default value is 25
 	/// Calculate the total number of points on all edges.
 	unsigned total_points(){
@@ -403,9 +448,6 @@ public:
 		}
 	}
-	// gaussian function
-	float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
-
     // convert vec3 to array
 	void stim2array(float *a, stim::vec3<T> b){
 		a[0] = b[0];
@@ -413,6 +455,16 @@ public:
 		a[2] = b[2];
 	}
+	// convert vec3 to array in bunch
+	void edge2array(T* a, std::vector< typename stim::vec3<T> > b){
+		size_t n = b.size();
+		for(size_t i = 0; i < n; i++){
+			a[i * 3 + 0] = b[i][0];
+			a[i * 3 + 1] = b[i][1];
+			a[i * 3 + 2] = b[i][2];	 
+		}
+	}
+
 	/// 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){
@@ -420,7 +472,7 @@ public:
 		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
+			M += E[e].integrate();							//get the integrated magnitude
 			L += E[e].length();							//get the edge length
 		}
@@ -428,9 +480,9 @@ public:
 	}
 	/// This function compares two networks and returns the percentage of the current network that is missing from A.
-
 	/// @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
+	/// @param device is the device that user want to use
 	stim::network<T> compare(stim::network<T> A, float sigma, int device){
 		stim::network<T> R;										//generate a network storing the result of the comparison
@@ -452,61 +504,363 @@ public:
 		}
 		//generate a KD-tree for network A
-		//float metric = 0.0;                               	// initialize metric to be returned after comparing the network
 		size_t MaxTreeLevels = 3;								// max tree level
 #ifdef __CUDACC__
 		cudaSetDevice(device);
-		stim::cuda_kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
+		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
+
+		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
+
+		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
+			size_t errormag_id = R.E[e].nmags() - 1;		//get the id for the new magnitude
+			
+			size_t n = R.E[e].size();						// the number of points in current edge
+			T* queryPt = new T[n];
+			T* m1 = new T[n];
+			T* dists = new T[n];
+			size_t* nnIdx = new size_t[n];
+
+			T* d_dists;										
+			T* d_m1;										
+			cudaMalloc((void**)&d_dists, n * sizeof(T));
+			cudaMalloc((void**)&d_m1, n * sizeof(T));
+
+			edge2array(queryPt, R.E[e]);
+			kdt.search(queryPt, n, nnIdx, dists);		
+
+			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
+
+			// configuration parameters
+			size_t threads = (1024>n)?n:1024;
+			size_t blocks = n/threads + (n%threads)?1:0;
+
+			d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma);					//calculate the metric value based on the distance
+
+			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
+
+			for(unsigned p = 0; p < n; p++){
+				R.E[e].set_mag(errormag_id, p, m1[p]);
+			}
+
+			//d_set_mag<<<blocks, threads>>>(R.E[e].M, errormag_id, n, m1);
+		}
+
+#else
+		stim::kdtree<T, 3> kdt;
+		kdt.create(c, n_data, MaxTreeLevels);
+	
+		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
+			size_t errormag_id = R.E[e].nmags() - 1;
+
+			size_t n = R.E[e].size();						// the number of points in current edge
+			T* query = new T[n];
+			T* m1 = new T[n];
+			T* dists = new T[n];
+			size* nnIdx = new size_t[n];
+
+			edge2array(queryPt, R.E[e]);
+
+			kdt.cpu_search(queryPt, n, nnIdx, dists);			//find the distance between A and the current network
+
+			for(unsigned p = 0; p < R.E[e].size(); p++){
+				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);		//calculate the metric value based on the distance
+				R.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
+			}
+		}
+#endif
+		return R;		//return the resulting network
+	}
+
+	/// This function compares two networks and split the current one according to the nearest neighbor of each point in each edge
+	/// @param A is the network to split
+	/// @param B is the corresponding mapping network
+	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
+	/// @param device is the device that user want to use
+	void split(stim::network<T> A, stim::network<T> B, float sigma, int device){
+
+		T *c;						                 	
+		size_t n_data = B.total_points();          				
+		c = (T*) malloc(sizeof(T) * n_data * 3); 				
+
+		size_t NB = B.E.size();								// the number of edges in B
+		unsigned t = 0;
+		for(unsigned e = 0; e < NB; e++){					// for every edge in B			
+			for(unsigned p = 0; p < B.E[e].size(); p++){	// for every points in B.E[e]
+				for(unsigned d = 0; d < 3; d++){				
+
+					c[t * 3 + d] = B.E[e][p][d];			// convert to array
+				}
+				t++;
+			}
+		}
+		size_t MaxTreeLevels = 3;							// max tree level
+
+#ifdef __CUDACC__
+		cudaSetDevice(device);
+		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
 		//compare each point in the current network to the field produced by A
 		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
-		T *dists = new T[1];								// near neighbor distances
-		size_t *nnIdx = new size_t[1];							// near neighbor indices // allocate near neigh indices
-		stim::vec3<T> p0, p1;
-		T m1;
-		//float M = 0;									//stores the total metric value
-		//float L = 0;									//stores the total network length
+		std::vector<std::vector<unsigned>> relation;		// the relationship between GT and T corresponding to NN
+		relation.resize(A.E.size());										
+
+		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
+			A.E[e].add_mag(0);								//add a new magnitude for the metric
+			size_t errormag_id = A.E[e].nmags() - 1;
+			
+			size_t n = A.E[e].size();						// the number of edges in A
+
+			T* queryPt = new T[n];							// set of all the points in current edge
+			T* m1 = new T[n];								// array of metrics for every point in current edge
+			T* dists = new T[n];							// store the distances for every point in current edge
+			size_t* nnIdx = new size_t[n];					// store the indices for every point in current edge
+			
+			// define pointers in device
+			T* d_dists;														
+			T* d_m1;
+			size_t* d_nnIdx;
+
+			// allocate memory for defined pointers
+			cudaMalloc((void**)&d_dists, n * sizeof(T));
+			cudaMalloc((void**)&d_m1, n * sizeof(T));
+			cudaMalloc((void**)&d_nnIdx, n * sizeof(size_t));
+
+			edge2array(queryPt, A.E[e]);						// convert edge to array
+			kdt.search(queryPt, n, nnIdx, dists);				// search the tree to find the NN for every point in current edge
+
+			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
+			cudaMemcpy(d_nnIdx, nnIdx, n * sizeof(size_t), cudaMemcpyHostToDevice);				// copy Idx from host to device
+
+			// configuration parameters
+			size_t threads = (1024>n)?n:1024;													// test to see whether the number of point in current edge is more than 1024
+			size_t blocks = n/threads + (n%threads)?1:0;
+
+			d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma);								// calculate the metrics in parallel
+
+			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
+
+			for(unsigned p = 0; p < n; p++){
+				A.E[e].set_mag(errormag_id, p, m1[p]);											// set the error(radius) value to every point in current edge
+			}
+
+			relation[e].resize(n);																// resize every edge relation size
+
+			unsigned* d_relation;
+			cudaMalloc((void**)&d_relation, n * sizeof(unsigned));								// allocate memory
+
+			std::vector<size_t> edge_point_num(NB);												// %th element is the number of points that %th edge has
+			for(unsigned ee = 0; ee < NB; ee++)
+				edge_point_num[ee] = B.E[ee].size();
+
+			size_t* d_edge_point_num;
+			cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t));
+			cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice);
+
+			d_findedge<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB);			// find the edge corresponding to the array index in parallel
+
+			cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost);	//copy relationship from device to host
+		}
+#else
+		stim::kdtree<T, 3> kdt;
+		kdt.create(c, n_data, MaxTreeLevels);
+	
+		std::vector<std::vector<unsigned>> relation;		// the mapping relationship between two networks
+		relation.resize(A.E.size());										
+		for(unsigned i = 0; i < A.E.size(); i++)
+			relation[i].resize(A.E[i].size());
+
+		std::vector<size_t> edge_point_num(NB);				//%th element is the number of points that %th edge has
+		for(unsigned ee = 0; ee < NB; ee++)
+			edge_point_num[ee] = B.E[ee].size();
+
 		T* queryPt = new T[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 e = 0; e < A.E.size(); e++){			//for each edge in A
+			A.E[e].add_mag(0);								//add a new magnitude for the metric
+			size_t errormag_id = A.E[e].nmags() - 1;
+
+			size_t n = A.E[e].size();						// the number of edges in A
+
+			T* queryPt = new T[n];
+			T* m1 = new T[n];
+			T* dists = new T[n];							//store the distances
+			size_t* nnIdx = new size_t[n];					//store the indices
+			
+			edge2array(queryPt, A.E[e]);
+			kdt.search(queryPt, n, nnIdx, dists);		
+
+			for(unsigned p = 0; p < A.E[e].size(); p++){
+				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
+				A.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
+				
+				unsigned id = 0;																	//mapping edge's idx
+				size_t num = 0;																		//total number of points before #th edge
+				for(unsigned i = 0; i < NB; i++){
+					num += B.E[i].size();
+					if(nnIdx[p] < num){																//find the edge it belongs to
+						relation[e][p] = id;
+						break;
+					}
+					id++;																			//current edge won't be the one, move to next edge
+				}
+			}
+		}
+#endif
+		E = A.E;
+		V = A.V;
+
+		unsigned int id = 0;									// split value								
+		for(unsigned e = 0; e < E.size(); e++){					// for every edge
+			for(unsigned p = 0; p < E[e].size() - 1; p++){		// for every point in each edge
+				if(relation[e][p] != relation[e][p + 1]){		// find the nearest edge changing point
+					id = p + 1;									// if there is no change in NN
+					if(id < E[e].size()/threshold_fac || (E[e].size() - id) < E[e].size()/threshold_fac)				// set tolerance to 10, tentatively--should find out the mathematical threshold!
+						id = E[e].size() - 1;																			// extreme situation is not acceptable
+					else
+						break;
+				}
+				if(p == E[e].size() - 2)																// if there is no splitting index, set the id to the last point index of current edge
+					id = p + 1;
+			}
+			unsigned errormag_id = E[e].nmags() - 1;
+			T G = 0;																					// test to see whether it has its nearest neighbor
+			for(unsigned i = 0; i < E[e].size(); i++)
+				G += E[e].m(i, errormag_id);															// won't split special edges
+			if(G / E[e].size() > metric_fac)															// should based on the color map
+				id = E[e].size() - 1;																	// set split idx to outgoing direction vertex
+
+			std::vector<edge> tmpe;
+			tmpe.resize(2);
+			tmpe = E[e].split(id);
+			vertex tmpv = stim::vec3<T>(-1, -1, 0);														// store the split point as vertex
+			if(tmpe.size() == 2){
+				relation.resize(relation.size() + 1);
+				for(unsigned d = id; d < E[e].size(); d++)
+					relation[relation.size() - 1].push_back(relation[e][d]);
+				tmpe[0].v[0] = E[e].v[0];																// begining vertex of first half edge -> original begining vertex
+				tmpe[1].v[1] = E[e].v[1];																// ending vertex of second half edge -> original ending vertex
+				tmpv = E[e][id];
+				V.push_back(tmpv);
+				tmpe[0].v[1] = (unsigned)V.size() - 1;													// ending vertex of first half edge -> new vertex
+				tmpe[1].v[0] = (unsigned)V.size() - 1;													// begining vertex of second half edge -> new vertex
+				edge tmp(E[e]);
+				E[e] = tmpe[0];																			// replace original edge by first half edge
+				E.push_back(tmpe[1]);																	// push second half edge to the last
+				V[V.size() - 1].e[1].push_back(e);														// push first half edge to the incoming of new vertex
+				V[V.size() - 1].e[0].push_back((unsigned)E.size() - 1);									// push second half edge to the outgoing of new vertex
+				for(unsigned i = 0; i < V[tmp.v[1]].e[1].size(); i++)									// find the incoming edge of original ending vertex
+					if(V[tmp.v[1]].e[1][i] == e)
+						V[tmp.v[1]].e[1][i] = (unsigned)V.size() - 1;
+			}
+		}
+	}
-			for(unsigned p = 0; p < R.E[e].size(); p++){				//for each point in the edge
+	/// This function compares two splitted networks and yields a mapping relationship between them according to NN
+	/// @param B is the network that the current network is going to map to
+	/// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B
+	/// @param device is the device that user want to use
+	void mapping(stim::network<T> B, std::vector<unsigned> &C, int device){
+		stim::network<T> A;								//generate a network storing the result of the comparison
+		A = (*this);
-				p1 = R.E[e][p];							//get the next point in the edge
-				stim2array(queryPt, p1);
-				kdt.search(queryPt, 1, nnIdx, dists);			//find the distance between A and the current network
+		size_t n = A.E.size();							// the number of edges in A
+		size_t NB = B.E.size();							// the number of edges in B
+
+		C.resize(A.E.size());	
+
+		T *c;						                 	// centerline (array of double pointers) - points on kdtree must be double
+		size_t n_data = B.total_points();          		// set the number of points
+		c = (T*) malloc(sizeof(T) * n_data * 3); 				
+
+		unsigned t = 0;
+		for(unsigned e = 0; e < NB; e++){					// for each edge in the network
+			for(unsigned p = 0; p < B.E[e].size(); p++){	// for each point in the edge
+				for(unsigned d = 0; d < 3; d++){			// for each coordinate
+
+					c[t * 3 + d] = B.E[e][p][d];
+				}
+				t++;
+			}
+		}
-				m1 = 1.0f - gaussianFunction((T)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
+		//generate a KD-tree for network A
+		//float metric = 0.0;                               		// initialize metric to be returned after comparing the network
+		size_t MaxTreeLevels = 3;									// max tree level
+		
+#ifdef __CUDACC__
+		cudaSetDevice(device);
+		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
+	
+		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
+		for(unsigned e = 0; e < n; e++){					//for each edge in A
+			size_t errormag_id = A.E[e].nmags() - 1;		//get the id for the new magnitude
+			
+			//pre-judge to get rid of impossibly mapping edges
+			T M = 0;
+			for(unsigned p = 0; p < A.E[e].size(); p++)
+				M += A.E[e].m(p, errormag_id);
+			M = M / A.E[e].size();
+			if(M > metric_fac)
+				C[e] = (unsigned)-1;						//set the nearest edge of impossibly mapping edges to maximum of unsigned
+			else{
+				T* queryPt = new T[1];
+				T* dists = new T[1];
+				size_t* nnIdx = new size_t[1];
+
+				stim2array(queryPt, A.E[e][A.E[e].size()/2]);
+				kdt.search(queryPt, 1, nnIdx, dists);
+				
+				unsigned id = 0;							//mapping edge's idx
+				size_t num = 0;								//total number of points before #th edge
+				for(unsigned i = 0; i < NB; i++){
+					num += B.E[i].size();
+					if(nnIdx[0] < num){
+						C[e] = id;
+						break;
+					}
+					id++;
+				}
 			}
 		}
 #else
-		stim::cpu_kdtree<T, 3> kdt;
+		stim::kdtree<T, 3> kdt;
 		kdt.create(c, n_data, MaxTreeLevels);
 		T *dists = new T[1];								// near neighbor distances
 		size_t *nnIdx = new size_t[1];						// near neighbor indices // allocate near neigh indices
 		stim::vec3<T> p0, p1;
-		T m1;
 		T* queryPt = new T[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
+		for(unsigned e = 0; e < R.E.size(); e++){			// for each edge in A
+			T M;											// the sum of metrics of current edge
+			unsigned errormag_id = R.E[e].nmags() - 1;
+			for(unsigned p = 0; p < R.E[e].size(); p++)
+				M += A.E[e].m(p, errormag_id);
+			M = M / A.E[e].size();
+			if(M > metric_fac)								// if the sum is larger than the metric_fac, we assume that it doesn't has corresponding edge in B
+				C[e] = (unsigned)-1;
+			else{											// if it should have corresponding edge in B, then...
+				p1 = R.E[e][R.E[e].size()/2];							
 				stim2array(queryPt, p1);
-				kdt.cpu_search(queryPt, 1, nnIdx, dists);			//find the distance between A and the current network
-
-				m1 = 1.0f - gaussianFunction((T)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
+				kdt.cpu_search(queryPt, 1, nnIdx, dists);	// search the tree		
+				
+				unsigned id = 0;							//mapping edge's idx
+				size_t num = 0;								//total number of points before #th edge
+				for(unsigned i = 0; i < NB; i++){
+					num += B.E[i].size();
+					if(nnIdx[0] < num){
+						C[e] = id;
+						break;
+					}
+					id++;
+				}
 			}
 		}
 #endif
-		return R;		//return the resulting network
 	}
 	/// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
@@ -616,4 +970,4 @@ public:
 	}
 };		//end stim::network class
 };		//end stim namespace
 -#endif
+#endif
 \ No newline at end of file
@@ -18,13 +18,14 @@ class circle : plane&lt;T&gt;
 private:
-	stim::vec3<T> Y;
+	//stim::vec3<T> Y;
+	T R;								//radius of the circle
-	CUDA_CALLABLE void
+	/*CUDA_CALLABLE void
 	init()
 	{
 		Y = U.cross(N).norm();
-	}
+	}*/
 public:
 	using stim::plane<T>::n;
@@ -42,26 +43,28 @@ public:
 		init();
 	}
-	///create a rectangle given a size and position in Z space.
+	///create a circle given a size and position in Z space.
 	///@param size: size of the rectangle in ND space.
 	///@param z_pos z coordinate of the rectangle.
 	CUDA_CALLABLE
-	circle(T size, T z_pos = (T)0) : plane<T>()
+	circle(T radius, T z_pos = (T)0) : plane<T>(z_pos)
 	{
-		center(stim::vec3<T>(0,0,z_pos));
-		scale(size);
-		init();
+		//center(stim::vec3<T>(0, 0, z_pos));
+		//scale(size);
+		//init();
+		R = radius;
 	}
 	///create a rectangle from a center point, normal
 	///@param c: x,y,z location of the center.
 	///@param n: x,y,z direction of the normal.	
 	CUDA_CALLABLE
-	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
+	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>(n, c)
 	{
-		center(c);
-		normal(n);
-		init();
+		//center(c);
+		//normal(n);
+		//init();
+		R = (T)1;
 	}
 	/*///create a rectangle from a center point, normal, and size
@@ -84,14 +87,18 @@ public:
 	///@param n: x,y,z direction of the normal.
 	///@param u: x,y,z direction for the zero vector (from where the rotation starts)
 	CUDA_CALLABLE
-	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
+	circle(vec3<T> c, T r, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
 	{
-		init();
+		P = c;
+		N = n;
 		setU(u);
+		R = r;
+		//init();
+		//setU(u);
 //		U = u;
-		center(c);
-		normal(n);
-		scale(s);
+		//center(c);
+		//normal(n);
+		//scale(s);
 	}
 	///scales the circle by a certain factor
@@ -99,23 +106,23 @@ public:
 	CUDA_CALLABLE
 	void scale(T factor)
 	{
-		U *= factor;
-		Y *= factor;
+		//U *= factor;
+		//Y *= factor;
+		R *= factor;
 	}
 	///sets the normal for the cirlce
 	///@param n: x,y,z direction of the normal.
 	CUDA_CALLABLE void
-	normal(vec3<T> n)
-	{
-		rotate(n, Y);
+	normal(vec3<T> n){
+		rotate(n);
 	}
 	///sets the center of the circle.
 	///@param n: x,y,z location of the center.
 	CUDA_CALLABLE void
 	center(vec3<T> p){
-		this->P = p;
+		P = p;
 	}
 	///boolean comparison
@@ -128,57 +135,63 @@ public:
 			return false;
 	}
+	//returns the point in world space corresponding to the polar coordinate (r, theta)
+	CUDA_CALLABLE stim::vec3<T>
+		p(T r, T theta) {
+		T u = r * cos(theta);				//calculate the coordinates in the planar space defined by the circle
+		T v = r * sin(theta);
+
+		vec3<T> V = U.cross(N);				//calculate the orthogonal vector V
+		return P + U * u + V * v;			//calculate the cartesian coordinate of the specified point
+	}
+
+	//returns the point in world space corresponding to the value theta at radius R
+	CUDA_CALLABLE stim::vec3<T>
+		p(T theta) {
+		return p(R, theta);
+	}
+
+	//get the world space value given the polar coordinates (r, theta)
+
 	///get the world space value given the planar coordinates a, b in [0, 1]
-	CUDA_CALLABLE stim::vec3<T> p(T a, T b)
+	/*CUDA_CALLABLE stim::vec3<T> p(T a, T b)
 	{
 		stim::vec3<T> result;
 		vec3<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
 		result = A + this->U * a + Y * b;
 		return result;
-	}
+	}*/
 	///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
-	CUDA_CALLABLE stim::vec3<T> operator()(T a, T b)
+	CUDA_CALLABLE stim::vec3<T> operator()(T r, T theta)
 	{
-		return p(a,b);
+		return p(r, theta);
+	}
+
+	//parenthesis operator returns the world space coordinate at the edge of the circle given theta
+	CUDA_CALLABLE stim::vec3<T> operator()(T theta) {
+		return p(theta);
 	}
 	///returns a vector with the points on the initialized circle.
 	///connecting the points results in a circle.
 	///@param n: integer for the number of points representing the circle.
-	std::vector<stim::vec3<T> >
-	getPoints(int n)
-	{
-		std::vector<stim::vec3<T> > result;
-		stim::vec3<T> point;
-		T x,y;
-		float step = 360.0/(float) n;
-		for(float j = 0; j <= 360.0; j += step)
-		{
-			y = 0.5*cos(j*stim::TAU/360.0)+0.5;
-			x = 0.5*sin(j*stim::TAU/360.0)+0.5;
-			result.push_back(p(x,y));
-		}
+	std::vector< stim::vec3<T> > points(unsigned n) {
+		std::vector< stim::vec3<T> > result(n);				//initialize a vector of n points
+
+		float dt = stim::TAU / n;
+		for (unsigned i = 0; i < n; i++)
+			result[i] = p(i * dt);							//calculate a point on the edge of the circle
 		return result;
 	}	
-	///returns a vector with the points on the initialized circle.
-	///connecting the points results in a circle.
-	///@param n: integer for the number of points representing the circle.
-	CUDA_CALLABLE stim::vec3<T>
-	p(T theta)
-	{
-		T x,y;
-		y = 0.5*cos(theta*STIM_TAU/360.0)+0.5;
-		x = 0.5*sin(theta*STIM_TAU/360.0)+0.5;
-		return p(x,y);
-	}
+	
 	std::string str() const
 	{
 		std::stringstream ss;
-		ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str();
+		ss << "r = "<<R<<"  (P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ")";
 		return ss.str();
 	}
+#ifndef STIM_CIRCLE_H
+#define STIM_CIRCLE_H
+
+#include <stim/cuda/cudatools/callable.h>
+#include <stim/math/plane.h>
+#include <stim/math/vector.h>
+#include <stim/math/triangle.h>
+#include <stim/math/constants.h>
+#include <assert.h>
+#include <algorithm>
+#include <iostream>
+
+namespace stim{
+
+template <typename T>
+class circle : plane<T>
+{
+
+private:
+	
+	stim::vec3<T> Y;
+
+	CUDA_CALLABLE void
+	init()
+	{
+		Y = U.cross(N).norm();
+	}
+
+public:
+	using stim::plane<T>::n;
+	using stim::plane<T>::P;
+	using stim::plane<T>::N;
+	using stim::plane<T>::U;
+	using stim::plane<T>::rotate;
+	using stim::plane<T>::setU;
+
+	///base constructor
+	///@param th value of the angle of the starting point from 0 to 360.
+	CUDA_CALLABLE
+	circle() : plane<T>()
+	{
+		init();
+	}
+
+	///create a rectangle given a size and position in Z space.
+	///@param size: size of the rectangle in ND space.
+	///@param z_pos z coordinate of the rectangle.
+	CUDA_CALLABLE
+	circle(T size, T z_pos = (T)0) : plane<T>()
+	{
+		center(stim::vec3<T>(0,0,z_pos));
+		scale(size);
+		init();
+	}
+
+	///create a rectangle from a center point, normal
+	///@param c: x,y,z location of the center.
+	///@param n: x,y,z direction of the normal.	
+	CUDA_CALLABLE
+	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
+	{
+		center(c);
+		normal(n);
+		init();
+	}
+
+	/*///create a rectangle from a center point, normal, and size
+	///@param c: x,y,z location of the center.
+	///@param s: size of the rectangle.
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE 
+	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
+	{
+		init();
+		center(c);
+		rotate(n, U, Y);
+		scale(s);
+	}
+	*/
+
+	///create a rectangle from a center point, normal, and size
+	///@param c: x,y,z location of the center.
+	///@param s: size of the rectangle.
+	///@param n: x,y,z direction of the normal.
+	///@param u: x,y,z direction for the zero vector (from where the rotation starts)
+	CUDA_CALLABLE
+	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
+	{
+		init();
+		setU(u);
+//		U = u;
+		center(c);
+		normal(n);
+		scale(s);
+	}
+
+	///scales the circle by a certain factor
+	///@param factor: the factor by which the dimensions of the shape are scaled.
+	CUDA_CALLABLE
+	void scale(T factor)
+	{
+		U *= factor;
+		Y *= factor;
+	}
+
+	///sets the normal for the cirlce
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE void
+	normal(vec3<T> n)
+	{
+		rotate(n, Y);
+	}
+
+	///sets the center of the circle.
+	///@param n: x,y,z location of the center.
+	CUDA_CALLABLE void
+	center(vec3<T> p){
+		this->P = p;
+	}
+
+	///boolean comparison
+	bool
+	operator==(const circle<T> & rhs)
+	{
+		if(P == rhs.P && U == rhs.U && Y == rhs.Y)
+			return true;
+		else
+			return false;
+	}
+
+	///get the world space value given the planar coordinates a, b in [0, 1]
+	CUDA_CALLABLE stim::vec3<T> p(T a, T b)
+	{
+		stim::vec3<T> result;
+
+		vec3<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
+		result = A + this->U * a + Y * b;
+		return result;
+	}
+
+	///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
+	CUDA_CALLABLE stim::vec3<T> operator()(T a, T b)
+	{
+		return p(a,b);
+	}
+
+	///returns a vector with the points on the initialized circle.
+	///connecting the points results in a circle.
+	///@param n: integer for the number of points representing the circle.
+	std::vector<stim::vec3<T> >
+	getPoints(int n)
+	{
+		std::vector<stim::vec3<T> > result;
+		stim::vec3<T> point;
+		T x,y;
+		float step = 360.0/(float) n;
+		for(float j = 0; j <= 360.0; j += step)
+		{
+			y = 0.5*cos(j*stim::TAU/360.0)+0.5;
+			x = 0.5*sin(j*stim::TAU/360.0)+0.5;
+			result.push_back(p(x,y));
+		}
+		return result;
+	}	
+	
+	///returns a vector with the points on the initialized circle.
+	///connecting the points results in a circle.
+	///@param n: integer for the number of points representing the circle.
+	CUDA_CALLABLE stim::vec3<T>
+	p(T theta)
+	{
+		T x,y;
+		y = 0.5*cos(theta*STIM_TAU/360.0)+0.5;
+		x = 0.5*sin(theta*STIM_TAU/360.0)+0.5;
+		return p(x,y);
+	}
+
+	std::string str() const
+	{
+		std::stringstream ss;
+		ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str();
+		return ss.str();
+	}
+
+};
+}
+#endif
@@ -46,7 +46,7 @@ namespace stim {
 	}
 	//Performs a convolution of a 2D image using the GPU. All pointers are assumed to be to memory on the current device.
-	//@param out is a pointer to the output image
+	//@param out is a pointer to the output image, which is of size (sx - kx + 1) x (sy - ky + 1)
 	//@param in is a pointer to the input image
 	//@param sx is the size of the input image along X
 	//@param sy is the size of the input image along Y
@@ -30,12 +30,18 @@ struct aabbn{
 		high[0] = x1;
 	}
-	CUDA_CALLABLE aabbn(T x0, T y0, T x1, T y1) : aabbn(x0, x1) {
+	CUDA_CALLABLE aabbn(T x0, T y0, T x1, T y1) {
+		low[0] = x0;
+		high[0] = x1;
 		low[1] = y0;
 		high[1] = y1;
 	}
-	CUDA_CALLABLE aabbn(T x0, T y0, T z0, T x1, T y1, T z1) : aabbn(x0, y0, x1, y1) {
+	CUDA_CALLABLE aabbn(T x0, T y0, T z0, T x1, T y1, T z1) {
+		low[0] = x0;
+		high[0] = x1;
+		low[1] = y0;
+		high[1] = y1;
 		low[2] = z0;
 		high[2] = z1;
 	}
@@ -3,56 +3,255 @@
 #include <iostream>
 #include <stim/math/circle.h>
 #include <stim/biomodels/centerline.h>
+#include <stim/visualization/obj.h>
-/*
-	
-*/
 namespace stim
 {
 template<typename T>
-class cylinder
- : public centerline<T>
-{
-	private:
-		stim::circle<T> s;							//an arbitrary circle
-		std::vector<stim::circle<T> > e;			//an array of circles that store the centerline
+class cylinder : public centerline<T> {
+protected:
+
+	std::vector< stim::vec3<T> > U;					//stores the array of U vectors defining the Frenet frame
+	std::vector< std::vector<T> > M;				//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
+
+	using stim::centerline<T>::findIdx;
+	
+	//calculates the U values for each point to initialize the frenet frame
+	//	this function assumes that the centerline has already been set
+	void init() {
+		U.resize(size());								//allocate space for the frenet frame vectors
+	
+		stim::circle<T> c;								//create a circle
+		stim::vec3<T> d0, d1;
+		for (size_t i = 0; i < size() - 1; i++) {		//for each line segment in the centerline
+			c.rotate(d(i));								//rotate the circle to match that normal
+			U[i] = c.U;									//save the U vector from the circle
+		}
+		U[size() - 1] = c.U;							//for the last point, duplicate the final frenet frame vector
+	}
-		std::vector<stim::vec3<T> > norms;
-		std::vector<stim::vec<T> > Us;
-		std::vector<stim::vec<T> > mags;			//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
-		std::vector< T > L;							//length of the cylinder at each position (pre-integration)
+public:
+	using stim::centerline<T>::size;
+	using stim::centerline<T>::at;
-		using stim::centerline<T>::c;
-		using stim::centerline<T>::N;
+	cylinder() : centerline(){}
-		///default init	
-		void
-		init()
-		{
+	cylinder(centerline& c) : centerline(c) {
+		init();
+	}
+
+	cylinder(centerline& c, T r) : centerline(c) {
+		init();
+		add_mag(r);
+	}
+
+	//initialize a cylinder with a list of points and magnitude values
+	cylinder(centerline& c, std::vector<T> r) : centerline(c){
+		init();
+		add_mag(r);
+	}
+
+	///Returns magnitude i at the given length into the fiber (based on the pvalue).
+	///Interpolates the position along the line.
+	///@param l: the location of the in the cylinder.
+	///@param idx: integer location of the point closest to l but prior to it.
+	T m(T l, int idx, size_t i) {
+		T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
+		T v2 = (M[idx][i] + (M[idx + 1][i] - M[idx][i])*a);
+		
+		return v2;
+	}
+
+	///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1).
+	///interpolates the radius along the line.
+	///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
+	T m(T pvalue, unsigned i = 0) {
+		if (pvalue <= 0.0) return M[0][i];
+		if (pvalue >= 1.0) return M[size() - 1][i];
+
+		T l = pvalue*L[L.size() - 1];
+		int idx = findIdx(l);
+		return m(l, idx, i);
+	}
+
+	///	Returns the magnitude at the given index
+	///	@param i is the index of the desired point
+	/// @param m is the index of the magnitude value
+	T m(unsigned i, unsigned m) {
+		return M[i][m];
+	}
+
+	///adds a magnitude to each point in the cylinder
+	void add_mag(T val = 0) {
+		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
+		for (size_t i = 0; i < size(); i++)		//for each point
+			M[i].push_back(val);				//add this value to the magnitude vector at each point
+	}
+
+	//adds a magnitude based on a list of magnitudes for each point
+	void add_mag(std::vector<T> val) {
+		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
+		for (size_t i = 0; i < size(); i++)		//for each point
+			M[i].push_back(val[i]);				//add this value to the magnitude vector at each point
+	}
+
+	//sets the value of magnitude m at point i
+	void set_mag(size_t m, size_t i, T v) {
+		M[i][m] = v;
+	}
+
+	size_t nmags() {
+		if (M.size() == 0) return 0;
+		else return M[0].size();
+	}
+
+	///Returns a circle representing the cylinder cross section at point i
+	stim::circle<T> circ(size_t i, size_t m = 0) {
+		return stim::circle<T>(at(i), M[i][m], d(i), U[i]);
+	}
+
+	///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m
+	stim::obj<T> OBJ(size_t N, size_t m = 0) {
+		stim::obj<T> out;								//create an OBJ object
+		stim::circle<T> c0, c1;
+		std::vector< stim::vec3<T> > p0, p1;			//allocate space for the point sets representing the circles bounding each cylinder segment
+		T u0, u1, v0, v1;											//allocate variables to store running texture coordinates
+		for (size_t i = 1; i < size(); i++) {			//for each line segment in the cylinder
+			c0 = circ(i - 1);							//get the two circles bounding the line segment
+			c1 = circ(i);
+
+			p0 = c0.points(N);							//get t points for each of the end caps
+			p1 = c1.points(N);
+
+			u0 = L[i - 1] / length();						//calculate the texture coordinate (u, v) where u runs along the cylinder length
+			u1 = L[i] / length();
+				
+			for (size_t n = 1; n < N; n++) {				//for each point in the circle
+				v0 = (double)(n-1) / (double)(N - 1);			//v texture coordinate runs around the cylinder
+				v1 = (double)(n) / (double)(N - 1);
+				out.Begin(OBJ_FACE);						//start a face (quad)
+					out.TexCoord(u0, v0);
+					out.Vertex(p0[n - 1][0], p0[n - 1][1], p0[n - 1][2]);	//output the points composing a strip of quads wrapping the cylinder segment
+					out.TexCoord(u1, v0);
+					out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]);
+				
+					out.TexCoord(u0, v1);
+					out.Vertex(p1[n][0], p1[n][1], p1[n][2]);
+					out.TexCoord(u1, v1);
+					out.Vertex(p0[n][0], p0[n][1], p0[n][2]);
+				out.End();
+			}
+			v0 = (double)(N - 2) / (double)(N - 1);			//v texture coordinate runs around the cylinder
+			v1 = 1.0;
+			out.Begin(OBJ_FACE);
+				out.TexCoord(u0, v0);
+				out.Vertex(p0[N - 1][0], p0[N - 1][1], p0[N - 1][2]);	//output the points composing a strip of quads wrapping the cylinder segment
+				out.TexCoord(u1, v0);
+				out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]);
+
+				out.TexCoord(u0, v1);
+				out.Vertex(p1[0][0], p1[0][1], p1[0][2]);
+				out.TexCoord(u1, v1);
+				out.Vertex(p0[0][0], p0[0][1], p0[0][2]);
+			out.End();
+		}
+		return out;
+	}
+
+	std::string str() {
+		std::stringstream ss;
+		size_t N = std::vector< stim::vec3<T> >::size();
+		ss << "---------[" << N << "]---------" << std::endl;
+		for (size_t i = 0; i < N; i++)
+			ss << std::vector< stim::vec3<T> >::at(i) << "   r = " << M[i][0] << "   u = " << U[i] << std::endl;
+		ss << "--------------------" << std::endl;
+
+		return ss.str();
+	}
+
+	/// Integrates a magnitude value along the cylinder.
+	/// @param m is the magnitude value to be integrated (this is usually the radius)
+	T integrate(size_t m = 0) {
+
+		T sum = 0;						//initialize the integral to zero
+		T m0, m1;						//allocate space for both magnitudes in a single segment
+		m0 = M[0][m];					//initialize the first point and magnitude to the first point in the cylinder
+		T len = L[0];					//allocate space for the segment length
+
+											
+		for (unsigned p = 1; p < size(); p++) {				//for every consecutive point in the cylinder
+			m1 = M[p][m];
+			if (p > 1) len = (L[p - 1] - L[p - 2]);		//calculate the segment length using the L array
+			sum += (m0 + m1) / (T)2.0 * len;				//add the average magnitude, weighted by the segment length
+			m0 = m1;									//move to the next segment by shifting points
+		}
+		return sum;			//return the integral
+	}
+
+	/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
+	///		centerline points are guaranteed to exist in the new cylinder
+	/// @param spacing is the maximum spacing allowed between sample points
+	cylinder<T> resample(T spacing) {
+		cylinder<T> c = stim::centerline<T>::resample(spacing);			//resample the centerline and use it to create a new cylinder
+
+		size_t nm = nmags();											//get the number of magnitude values in the current cylinder
+		if (nm > 0) {													//if there are magnitude values
+			std::vector<T> magvec(nm, 0);							//create a magnitude vector for a single point
+			c.M.resize(c.size(), magvec);									//allocate space for a magnitude vector at each point of the new cylinder
+		}
+		T l, t;
+		for (size_t i = 0; i < c.size(); i++) {							//for each point in the new cylinder
+			l = c.L[i];													//get the length along the new cylinder
+			t = l / length();										//calculate the parameter value along the new cylinder
+			for (size_t mag = 0; mag < nm; mag++) {							//for each magnitude value
+				c.M[i][mag] = m(t, mag);								//retrieve the interpolated magnitude from the current cylinder and store it in the new one
+			}
 		}
-		///Calculate the physical length of the fiber at each point in the fiber.
-		void
-		calculateL()
-		{
-/*			L.resize(inP.size());
-			T temp = (T)0;
-			L[0] = 0;
-			for(size_t i = 1; i < L.size(); i++)
-			{
-				temp += (inP[i-1] - inP[i]).len();
-				L[i] = temp;
+
+		return c;
+	}
+
+	std::vector< stim::cylinder<T> > split(unsigned int idx){
+			
+		unsigned N = size();
+		std::vector< stim::centerline<T> > LL;
+		LL.resize(2);
+		LL = (*this).centerline<T>::split(idx);
+		std::vector< stim::cylinder<T> > C(LL.size());
+		unsigned i = 0;
+		C[0] = LL[0];
+		C[0].M.resize(idx);
+		for( ; i < idx; i++){
+			//for(unsigned d = 0; d < 3; d++)
+				//C[0][i][d] = LL[0].c[i][d];
+			C[0].M[i] = M[i];
+			C[0].M[i].resize(1);
+		}
+		if(C.size() == 2){
+			C[1] = LL[1];
+			C[1].M.resize(N  - idx);
+			for( ; i < N; i++){
+				//for(unsigned d = 0; d < 3; d++)
+					//C[1][i - idx][d] = LL[1].c[i - idx][d];
+				C[1].M[i - idx] = M[i];
+				C[1].M[i - idx].resize(1);
 			}
-*/		}
+		}
+
+		return C;
+	}
+
-		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
+		/*
+		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM)
 		void
-		init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
-		{
-			mags = inM;
+		init(centerline inP, std::vector< std::vector<T> > inM) {
+			M = inM;									//the magnitude vector can be copied directly
+			(*this) = inP;								//the centerline can be copied to this class directly
 			stim::vec3<float> v1;
 			stim::vec3<float> v2;
 			e.resize(inP.size());
@@ -246,13 +445,18 @@ class cylinder
 			init(inP, inM);
 		}
+		//assignment operator creates a cylinder from a centerline (default radius is zero)
+		cylinder& operator=(stim::centerline<T> c) {
+			init(c);
+		}
+
 		///Returns the number of points on the cylinder centerline
 		unsigned int size(){
 			return N;
 		}
-
+		
 		///Returns a position vector at the given p-value (p value ranges from 0 to 1).
 		///interpolates the position along the line.
 		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
@@ -266,22 +470,7 @@ class cylinder
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
 			return (p(l,idx));
-/*				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-				stim::vec3<T> v1(
-						c[idx][0],	
-						c[idx][1],	
-						c[idx][2]
-						);	
-						
-				stim::vec3<T> v2(
-						c[idx+1][0],	
-						c[idx+1][1],	
-						c[idx+1][2]
-						);
-
-//			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
-			return( v1 + (v2 - v1)*rat);
-*/		}
+		}
 		///Returns a position vector at the given length into the fiber (based on the pvalue).
 		///Interpolates the radius along the line.
@@ -321,10 +510,7 @@ class cylinder
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
 			return (r(l,idx));
-/*			T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-//			return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
-			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
-*/		}
+		}
 		///Returns a radius at the given length into the fiber (based on the pvalue).
 		///Interpolates the position along the line.
@@ -506,10 +692,10 @@ class cylinder
 			return cylinder<T>(result);
-		}
+		}*/
 };
 }
 -#endif
+#endif
 \ No newline at end of file
+#ifndef STIM_CYLINDER_H
+#define STIM_CYLINDER_H
+#include <iostream>
+#include <stim/math/circle.h>
+#include <stim/biomodels/centerline.h>
+
+
+namespace stim
+{
+template<typename T>
+class cylinder
+ : public centerline<T>
+{
+	private:
+		stim::circle<T> s;							//an arbitrary circle
+		std::vector<stim::circle<T> > e;			//an array of circles that store the centerline
+
+		std::vector<stim::vec3<T> > norms;
+		std::vector<stim::vec<T> > Us;
+		std::vector<stim::vec<T> > mags;			//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
+		std::vector< T > L;							//length of the cylinder at each position (pre-integration)
+
+
+		using stim::centerline<T>::c;
+		using stim::centerline<T>::N;
+	
+		///default init	
+		void
+		init()
+		{
+
+		}
+
+		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
+		void
+		init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
+		{
+			mags = inM;
+			stim::vec3<float> v1;
+			stim::vec3<float> v2;
+			e.resize(inP.size());
+
+			norms.resize(inP.size());
+			Us.resize(inP.size());
+
+			if(inP.size() < 2)
+				return;
+
+			//calculate each L.
+			L.resize(inP.size());						//the number of precomputed lengths will equal the number of points
+			T temp = (T)0;								//length up to that point
+			L[0] = temp;
+			for(size_t i = 1; i < L.size(); i++)
+			{
+				temp += (inP[i-1] - inP[i]).len();
+				L[i] = temp;
+			}
+
+			stim::vec3<T> dr = (inP[1] - inP[0]).norm();
+			s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec3<T>(1,0,0));
+			norms[0] = s.N;
+			e[0] = s;
+			Us[0] = e[0].U;
+			for(size_t i = 1; i < inP.size()-1; i++)
+			{
+				s.center(inP[i]);
+				v1 = (inP[i] - inP[i-1]).norm();
+				v2 = (inP[i+1] - inP[i]).norm();
+				dr = (v1+v2).norm();
+				s.normal(dr);
+
+				norms[i] = s.N;
+
+				s.scale(inM[i][0]/inM[i-1][0]);
+				e[i] = s;
+				Us[i] = e[i].U;
+			}
+			
+			int j = inP.size()-1;
+			s.center(inP[j]);
+			dr = (inP[j] - inP[j-1]).norm();
+			s.normal(dr);
+
+			norms[j] = s.N;
+
+			s.scale(inM[j][0]/inM[j-1][0]);
+			e[j] = s;
+			Us[j] = e[j].U;
+		}
+		
+		///returns the direction vector at point idx.
+		stim::vec3<T>
+		d(int idx)
+		{
+			if(idx == 0)
+			{
+				stim::vec3<T> temp(
+						c[idx+1][0]-c[idx][0],	
+						c[idx+1][1]-c[idx][1],	
+						c[idx+1][2]-c[idx][2]
+						);	
+//				return (e[idx+1].P - e[idx].P).norm();
+				return (temp.norm());
+			}
+			else if(idx == N-1)
+			{
+				stim::vec3<T> temp(
+						c[idx][0]-c[idx+1][0],	
+						c[idx][1]-c[idx+1][1],	
+						c[idx][2]-c[idx+1][2]
+						);	
+			//	return (e[idx].P - e[idx-1].P).norm();
+				return (temp.norm());
+			}
+			else
+			{
+//				return (e[idx+1].P - e[idx].P).norm();
+//				stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
+				stim::vec3<T> v1(
+						c[idx][0]-c[idx-1][0],	
+						c[idx][1]-c[idx-1][1],	
+						c[idx][2]-c[idx-1][2]
+						);	
+						
+//				stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
+				stim::vec3<T> v2(
+						c[idx+1][0]-c[idx][0],	
+						c[idx+1][1]-c[idx][1],	
+						c[idx+1][2]-c[idx][2]
+						);
+					
+				return (v1.norm()+v2.norm()).norm();			
+			} 
+	//		return e[idx].N;	
+
+		}
+
+		stim::vec3<T>
+		d(T l, int idx)
+		{
+			if(idx == 0 || idx == N-1)
+			{
+				return norms[idx];
+//				return e[idx].N;
+			}
+			else
+			{
+				
+				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+				return(	norms[idx] + (norms[idx+1] - norms[idx])*rat);
+//				return(	e[idx].N + (e[idx+1].N - e[idx].N)*rat);
+			} 
+		}
+
+
+		///finds the index of the point closest to the length l on the lower bound.
+		///binary search.
+		int
+		findIdx(T l)
+		{
+			unsigned int i = L.size()/2;
+			unsigned int max = L.size()-1;
+			unsigned int min = 0;
+			while(i > 0 && i < L.size()-1)
+			{
+//				std::cerr << "Trying " << i << std::endl;
+//				std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
+				if(l < L[i])
+				{
+					max = i;
+					i = min+(max-min)/2;
+				}
+				else if(L[i] <= l && L[i+1] >= l)
+				{
+					break;
+				}
+				else
+				{
+					min = i;
+					i = min+(max-min)/2;
+				}
+			}
+			return i;
+		}
+
+	public:
+		///default constructor
+		cylinder()
+		// : centerline<T>()
+		{
+
+		}
+
+		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
+		///@param inP:  Vector of stim vec3 composing the points of the centerline.
+		///@param inM:  Vector of stim vecs composing the radii of the centerline.
+		cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
+			: centerline<T>(inP)
+		{
+			init(inP, inM);
+		}
+
+		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
+		///@param inP:  Vector of stim vec3 composing the points of the centerline.
+		///@param inM:  Vector of stim vecs composing the radii of the centerline.
+		cylinder(std::vector<stim::vec3<T> > inP, std::vector< T > inM)
+			: centerline<T>(inP)
+		{
+			std::vector<stim::vec<T> > temp;
+			stim::vec<T> zero(0.0,0.0);
+			temp.resize(inM.size(), zero);
+			for(int i = 0; i < inM.size(); i++)
+				temp[i][0] = inM[i];
+			init(inP, temp);
+		}
+
+
+		///Constructor defines a cylinder with centerline inP and magnitudes of zero
+		///@param inP: Vector of stim vec3 composing the points of the centerline
+		cylinder(std::vector< stim::vec3<T> > inP)
+			: centerline<T>(inP)
+		{
+			std::vector< stim::vec<T> > inM;						//create an array of arbitrary magnitudes
+
+			stim::vec<T> zero;
+			zero.push_back(0);
+
+			inM.resize(inP.size(), zero);								//initialize the magnitude values to zero
+			init(inP, inM);
+		}
+
+		///Returns the number of points on the cylinder centerline
+
+		unsigned int size(){
+			return N;
+		}
+
+
+		///Returns a position vector at the given p-value (p value ranges from 0 to 1).
+		///interpolates the position along the line.
+		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
+		stim::vec3<T>
+		p(T pvalue)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0)
+			{
+				return stim::vec3<float>(-1,-1,-1);
+			}
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			return (p(l,idx));
+/*				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+				stim::vec3<T> v1(
+						c[idx][0],	
+						c[idx][1],	
+						c[idx][2]
+						);	
+						
+				stim::vec3<T> v2(
+						c[idx+1][0],	
+						c[idx+1][1],	
+						c[idx+1][2]
+						);
+
+//			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
+			return( v1 + (v2 - v1)*rat);
+*/		}
+
+		///Returns a position vector at the given length into the fiber (based on the pvalue).
+		///Interpolates the radius along the line.
+		///@param l: the location of the in the cylinder.
+		///@param idx: integer location of the point closest to l but prior to it.
+		stim::vec3<T>
+		p(T l, int idx)
+		{
+				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+				stim::vec3<T> v1(
+						c[idx][0],	
+						c[idx][1],	
+						c[idx][2]
+						);	
+						
+				stim::vec3<T> v2(
+						c[idx+1][0],	
+						c[idx+1][1],	
+						c[idx+1][2]
+						);
+//			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
+			return(	v1 + (v2-v1)*rat);
+//			return(
+//			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
+		}
+
+		///Returns a radius at the given p-value (p value ranges from 0 to 1).
+		///interpolates the radius along the line.
+		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
+		T
+		r(T pvalue)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0){
+				std::cerr<<"Error, value "<<pvalue<<" is outside of [0 1]."<<std::endl;
+				exit(1);
+			}
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			return (r(l,idx));
+/*			T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+//			return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
+			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
+*/		}
+
+		///Returns a radius at the given length into the fiber (based on the pvalue).
+		///Interpolates the position along the line.
+		///@param l: the location of the in the cylinder.
+		///@param idx: integer location of the point closest to l but prior to it.
+		T
+		r(T l, int idx)
+		{
+				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+			T v1 = (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
+			T v3 = (Us[idx].len() + (Us[idx+1].len() - Us[idx].len())*rat);
+			T v2 = (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
+//			std::cout << (float)v1 = (float) v2 << std::endl;
+			if(abs(v3 - v1) >= 10e-6)
+			{
+				std::cout << "-------------------------" << std::endl;
+				std::cout << e[idx].str() << std::endl << std::endl;
+				std::cout << Us[idx].str() << std::endl;
+				std::cout << (float)v1 - (float) v2 << std::endl;
+				std::cout << "failed" << std::endl;
+			}
+//			std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl;
+//			std::cout << v2 << std::endl;
+			return(v2);
+//			return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
+	//	(
+		}
+
+		///	Returns the magnitude at the given index
+		///	@param i is the index of the desired point
+		/// @param m is the index of the magnitude value
+		T ri(unsigned i, unsigned m = 0){
+			return mags[i][m];
+		}
+
+		/// Adds a new magnitude value to all points
+		/// @param m is the starting value for the new magnitude
+		void add_mag(T m = 0){
+			for(unsigned int p = 0; p < N; p++)
+				mags[p].push_back(m);
+		}
+
+		/// Sets a magnitude value
+		/// @param val is the new value for the magnitude
+		/// @param p is the point index for the magnitude to be set
+		/// @param m is the index for the magnitude
+		void set_mag(T val, unsigned p, unsigned m = 0){
+			mags[p][m] = val;
+		}
+
+		/// Returns the number of magnitude values at each point
+		unsigned nmags(){
+			return mags[0].size();
+		}
+
+		///returns the position of the point with a given pvalue and theta on the surface
+		///in x, y, z coordinates. Theta is in degrees from 0 to 360.
+		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
+		///@param theta: the angle to the point of a circle.
+		stim::vec3<T>
+		surf(T pvalue, T theta)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0)
+			{
+				return stim::vec3<float>(-1,-1,-1);
+			} else {
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			stim::vec3<T> ps = p(l, idx); 
+			T m = r(l, idx);
+			s = e[idx];
+			s.center(ps);
+			s.normal(d(l, idx));
+			s.scale(m/e[idx].U.len());
+			return(s.p(theta));
+			}
+		}
+
+		///returns a vector of points necessary to create a circle at every position in the fiber.
+		///@param sides: the number of sides of each circle.	
+		std::vector<std::vector<vec3<T> > >
+		getPoints(int sides)
+		{
+			std::vector<std::vector <vec3<T> > > points;
+			points.resize(N);
+			for(int i = 0; i < N; i++)
+			{
+				points[i] = e[i].getPoints(sides);
+			}
+			return points;
+		}
+
+		///returns the total length of the line at index j.
+		T
+		getl(int j)
+		{
+			return (L[j]);
+		}
+		/// Allows a point on the centerline to be accessed using bracket notation
+
+		vec3<T> operator[](unsigned int i){
+			return e[i].P;
+		}
+
+		/// Returns the total length of the cylinder centerline
+		T length(){
+			return L.back();
+		}
+
+		/// Integrates a magnitude value along the cylinder.
+		/// @param m is the magnitude value to be integrated (this is usually the radius)
+		T integrate(unsigned m = 0){
+
+			T M = 0;						//initialize the integral to zero
+			T m0, m1;						//allocate space for both magnitudes in a single segment
+
+			//vec3<T> p0, p1;					//allocate space for both points in a single segment
+
+			m0 = mags[0][m];				//initialize the first point and magnitude to the first point in the cylinder
+			//p0 = pos[0];
+
+			T len = L[0];						//allocate space for the segment length
+
+			//for every consecutive point in the cylinder
+			for(unsigned p = 1; p < N; p++){
+
+				//p1 = pos[p];							//get the position and magnitude for the next point
+				m1 = mags[p][m];
+
+				if(p > 1) len = (L[p-1] - L[p-2]);		//calculate the segment length using the L array
+
+				//add the average magnitude, weighted by the segment length
+				M += (m0 + m1)/(T)2.0 * len;
+
+				m0 = m1;								//move to the next segment by shifting points
+			}
+			return M;			//return the integral
+		}
+
+		/// Averages a magnitude value across the cylinder
+		/// @param m is the magnitude value to be averaged (this is usually the radius)
+		T average(unsigned m = 0){			
+
+			//return the average magnitude
+			return integrate(m) / L.back();
+		}
+
+		/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
+		///		centerline points are guaranteed to exist in the new cylinder
+		/// @param spacing is the maximum spacing allowed between sample points
+		cylinder<T> resample(T spacing){
+
+			std::vector< vec3<T> > result;
+
+			vec3<T> p0 = e[0].P;									//initialize p0 to the first point on the centerline
+			vec3<T> p1;
+			unsigned N = size();									//number of points in the current centerline
+
+			//for each line segment on the centerline
+			for(unsigned int i = 1; i < N; i++){
+				p1 = e[i].P;										//get the second point in the line segment
+
+				vec3<T> v = p1 - p0;								//calculate the vector between these two points
+				T d = v.len();										//calculate the distance between these two points (length of the line segment)
+
+				size_t nsteps = (size_t)std::ceil(d / spacing);		//calculate the number of steps to take along the segment to meet the spacing criteria
+				T stepsize = (T)1.0 / nsteps;						//calculate the parametric step size between new centerline points
+
+				//for each step along the line segment
+				for(unsigned s = 0; s < nsteps; s++){
+					T alpha = stepsize * s;							//calculate the fraction of the distance along the line segment covered
+					result.push_back(p0 + alpha * v);				//push the point at alpha position along the line segment
+				}
+
+				p0 = p1;											//shift the points to move to the next line segment
+			}
+
+			result.push_back(e[size() - 1].P);						//push the last point in the centerline
+
+			return cylinder<T>(result);
+
+		}
+
+		
+};
+
+}
+#endif
@@ -49,20 +49,72 @@ public:
 	void glCenterline(unsigned m = 0){
 		if(!glIsList(dlist)){					//if dlist isn't a display list, create it
-			dlist = glGenLists(1);				//generate a display list
+			dlist = glGenLists(3);				//generate a display list
 			glNewList(dlist, GL_COMPILE);		//start a new display list
 			for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+				unsigned errormag_id = E[e].nmags() - 1;
 				glBegin(GL_LINE_STRIP);
 				for(unsigned p = 0; p < E[e].size(); p++){				//for each point on that edge
 					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
-					glTexCoord1f(E[e].ri(p, m));						//set the texture coordinate based on the specified magnitude index
+					glTexCoord1f(E[e].m(p, errormag_id));						//set the texture coordinate based on the specified magnitude index
 				}
 				glEnd();
 			}
 			glEndList();						//end the display list
 		}		
 		glCallList(dlist);					// render the display list
+	}
+	void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
+		if(!glIsList(dlist1)){
+			glNewList(dlist1, GL_COMPILE);
+			for(unsigned e = 0; e < E.size(); e++){
+				if(map[e] != unsigned(-1)){
+					glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
+					glBegin(GL_LINE_STRIP);
+					for(unsigned p = 0; p < E[e].size(); p++){
+						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
+					}
+					glEnd();
+				}
+				else{
+					glColor3f(1.0, 1.0, 1.0);
+					glBegin(GL_LINE_STRIP);
+					for(unsigned p = 0; p < E[e].size(); p++){
+						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
+					}
+					glEnd();
+				}
+			}
+			glEndList();
+		}
+		glCallList(dlist1);
+	}
+
+	void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap){
+		if(!glIsList(dlist2)){
+			glNewList(dlist2, GL_COMPILE);
+			for(unsigned e = 0; e < E.size(); e++){
+				if(map[e] != unsigned(-1)){
+					glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
+					glBegin(GL_LINE_STRIP);
+					for(unsigned p = 0; p < E[e].size(); p++){
+						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
+					}
+					glEnd();
+				}
+				else{
+					glColor3f(1.0, 1.0, 1.0);
+					glBegin(GL_LINE_STRIP);
+					for(unsigned p = 0; p < E[e].size(); p++){
+						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
+					}
+					glEnd();
+				}
+			}
+			glEndList();
+		}
+		glCallList(dlist2);
 	}
 };		//end stim::gl_network class