cylinder can build OBJ objects, simplified class hierarchy

David Mayerich
1 parent df171bd2
Showing 6 changed files with 1300 additions and 352 deletions Show diff stats
stim/biomodels/centerline.h
stim/biomodels/centerline_dep.h
stim/math/circle.h
stim/math/circle_dep.h
stim/visualization/cylinder.h
stim/visualization/cylinder_dep.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){
@@ -216,123 +143,84 @@ public:
 			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;
 			//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,
+		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
-		std::vector<stim::vec<T> > fiberPositions = centerline();
-		std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
+
+		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
@@ -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
@@ -3,56 +3,163 @@
 #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
-
-		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;
+class cylinder : public centerline<T> {
+protected:
+	//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::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)
+	//std::vector< T > L;								//length of the cylinder at each position (pre-integration)
+
+
+	//using stim::centerline<T>::c;
+	//using stim::centerline<T>::N;
+	using stim::centerline<T>::size;
+	using stim::centerline<T>::at;
+	using stim::centerline<T>::findIdx;
-		///default init	
-		void
-		init()
-		{
-
+	//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
 		}
-
-		///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;
+		U[size() - 1] = c.U;							//for the last point, duplicate the final frenet frame vector
+	}
+
+	///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
+	}
+
+public:
+	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();
+	}
+
+	cylinder(centerline& c, T r = 0) : 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 = 0) {
+		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, size_t 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 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;
+	}
-		///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 +353,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 +378,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 +418,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,7 +600,7 @@ class cylinder
 			return cylinder<T>(result);
-		}
+		}*/
 };
+#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