merged the two stim versions to work for volume-spider, and ONLY volume spider. …

…circle.h, cylinder.h, centerline.h and some other have been changed. Other compilation errors on linux also fixxed in spharmonic.h and gl_spharmonic.h

merged the two stim versions to work for volume-spider, and ONLY volume spider. …
…circle.h, cylinder.h, centerline.h and some other have been changed. Other compilation errors on linux also fixxed in spharmonic.h and gl_spharmonic.h
Pavel Govyadinov
1 parent a8ad9192
Showing 11 changed files with 551 additions and 596 deletions Show diff stats
stim/biomodels/centerline.h
stim/cuda/filter.cuh
stim/cuda/ivote/local_max.cuh
stim/gl/gl_spider.h
stim/image/bmp.h
stim/image/image.h
stim/math/circle.h
stim/math/random.h
stim/math/spharmonics.h
stim/visualization/cylinder.h
stim/visualization/gl_spharmonics.h
@@ -3,6 +3,7 @@
 #include <vector>
 #include <stim/math/vec3.h>
+//#include <ANN/ANN.h>
 namespace stim{
@@ -11,134 +12,195 @@ namespace stim{
  *	class to describe an interconnected (often biological) network.
  */
 template<typename T>
-class centerline : public std::vector< stim::vec3<T> >{
+class centerline{
 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;
+	}
-	std::vector<T> L;										//stores the integrated length along the fiber (used for parameterization)
-
-	///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)).norm();	//the last direction vector is oriented towards the last line segment
-
-		//all other direction vectors are the average direction of the two joined line segments
-		vec3<T> a = at(i) - at(i - 1);
-		vec3<T> b = at(i + 1) - at(i);
-		vec3<T> ab = a.norm() + b.norm();
-		return ab.norm();
-	}
-	//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++;
-		}
+	/// 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
+//		kdt = NULL;
+		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);
-		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
+		///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
 		}
+//		if(kd)
+//			gen_kdtree();							//generate the kd tree for the new fiber
 	}
-	void init() {
-		if (size() == 0) return;								//return if there aren't any points
-		update_L();
+	/// 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;
+		}
+		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;
+//	}
+
 	/// 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){
-		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) {
-		for (size_t i = 1; i < L.size(); i++) {				//for each point in the centerline
-			if (L[i] > l) return i - 1;						//if we have passed the desired length value, return i-1
-		}
-		return L.size() - 1;
-		/*size_t i = L.size() / 2;
-		size_t max = L.size() - 1;
-		size_t min = 0;
-		while (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;*/
-	}
+		stim::vec3<T> r;
+		r.resize(3);
+		r[0] = c[i][0];
+		r[1] = c[i][1];
+		r[2] = c[i][2];
-	///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);
+		return r;
 	}
 public:
-	using std::vector< stim::vec3<T> >::at;
-	using std::vector< stim::vec3<T> >::size;
-
-	centerline() : std::vector< stim::vec3<T> >() {
+	centerline(){
 		init();
 	}
-	centerline(size_t n) : std::vector< stim::vec3<T> >(n){
-		init();
+
+	/// Copy constructor
+	centerline(const stim::centerline<T> &obj){
+		copy(obj);
 	}
-	//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);
+	//temp constructor for graph visualization
+	centerline(int n)
+	{
+		init(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) return at(0);			//return the first element
-		if (pvalue >= 1.0) return back();			//return the last element
+	/// 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++){
-		T l = pvalue*L[L.size() - 1];
-		int idx = findIdx(l);
-		return p(l, idx);
+			//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();
 	}
-	///Update centerline internal parameters (currently the L vector)
-	void update() {
-		init();
+	/// 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];
+			//set the radius
+		}
+
+		//generate a kd tree
+		//if(kd)
+		//	gen_kdtree();
+	}
+
+	/// Assignment operation
+	centerline& operator=(const centerline &rhs){
+		if(this == &rhs) return *this;			//test for and handle self-assignment
+		copy(rhs);
+		return *this;
 	}
-	///Return the length of the entire centerline
-	T length() {
-		return L.back();
+
+
+	/// 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;
 	}
 	/// 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
-		size_t N = size();
+		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){
@@ -154,84 +216,123 @@ public:
 			fl.resize(2);								//set the array size to 2
-			fl[0] = stim::centerline<T>(N1);			//set the size of each fiber
-			fl[1] = stim::centerline<T>(N2);
+			fl[0].init(N1);								//set the size of each fiber
+			fl[1].init(N2);
 			//copy both halves of the fiber
-			unsigned int i;
+			unsigned int i, d;
 			//first half
-			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
+			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++)
-				fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
-			fl[1].init();							//initialize the length vector
+			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
 	}
+	/// 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;
-		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;
+
+		//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();
 	}
-
-	/// Back method returns the last point in the fiber
-	stim::vec3<T> back(){
-		return std::vector< stim::vec3<T> >::back();
+	/// Returns the number of centerline points in the fiber
+	unsigned int size(){
+		return N;
 	}
-		////resample a fiber in the network
-	stim::centerline<T> resample(T spacing)
-	{	
-		//std::cout<<"fiber::resample()"<<std::endl;
-		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
+	/// Bracket operator returns the element at index i
-		size_t N = size();
+	/// @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);
+	}
-		centerline<T> new_c; // initialize list of new resampled points on the fiber
+	/// 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)
+		//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;
+		{
-			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);
+			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
-				new_c.push_back(at(f));
-		}
-		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;
+				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;
 	}
 };
@@ -180,6 +180,7 @@ namespace stim
 		#endif
+//		stim::cuda::gpu_local_max<float>(centers, res, conn, DIM_X, DIM_Y);
 		stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
 		cudaDeviceSynchronize();
@@ -10,24 +10,24 @@ namespace stim{
 		// this kernel calculates the local maximum for finding the cell centers
 		template<typename T>
-		__global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, size_t x, size_t y){
+		__global__ void cuda_local_max(T* gpuCenters, T* gpuVote, T final_t, int conn, int x, int y){
 			// calculate the 2D coordinates for this current thread.
-			size_t xi = blockIdx.x * blockDim.x + threadIdx.x;
-			size_t yi = blockIdx.y * blockDim.y + threadIdx.y;
+			int xi = blockIdx.x * blockDim.x + threadIdx.x;
+			int yi = blockIdx.y * blockDim.y + threadIdx.y;
 			if(xi >= x || yi >= y)
 				return;
 			// convert 2D coordinates to 1D
-			size_t i = yi * x + xi;
+			int i = yi * x + xi;			
 			gpuCenters[i] = 0;		//initialize the value at this location to zero
 			T val = gpuVote[i];
 			//compare to the threshold
-			//if(val < final_t) return;
+			if(val < final_t) return;
 			//define an array to store indices with same vote value
 			/*int * IdxEq;
@@ -56,25 +56,24 @@ namespace stim{
 					return;
 				}
 			}	*/	
-			//gpuCenters[i] = 1;
-			gpuCenters[i] = gpuVote[i];
+			gpuCenters[i] = 1;
 		}
 		template<typename T>
-		void gpu_local_max(T* gpuCenters, T* gpuVote,  unsigned int conn, size_t x, size_t y){
+		void gpu_local_max(T* gpuCenters, T* gpuVote, T final_t, unsigned int conn, unsigned int x, unsigned int y){
 			unsigned int max_threads = stim::maxThreadsPerBlock();
 			/*dim3 threads(max_threads, 1);
 			dim3 blocks(x/threads.x + (x %threads.x == 0 ? 0:1) , y);*/
-			dim3 threads((unsigned int)sqrt(max_threads), (unsigned int)sqrt(max_threads) );
-			dim3 blocks((unsigned int)x/threads.x + 1, (unsigned int)y/threads.y + 1);
+			dim3 threads( sqrt(max_threads), sqrt(max_threads) );
+			dim3 blocks(x/threads.x + 1, y/threads.y + 1);
 			//call the kernel to find the local maximum.
-			cuda_local_max <<< blocks, threads >>>(gpuCenters, gpuVote, conn, x, y);
+			cuda_local_max <<< blocks, threads >>>(gpuCenters, gpuVote, final_t, conn, x, y);
 		}
 		template<typename T>
-		void cpu_local_max(T* cpuCenters, T* cpuVote, unsigned int conn, unsigned int x, unsigned int y){
+		void cpu_local_max(T* cpuCenters, T* cpuVote, T final_t, unsigned int conn, unsigned int x, unsigned int y){
 			//calculate the number of bytes in the array
 			unsigned int bytes = x * y * sizeof(T);
@@ -91,7 +90,7 @@ namespace stim{
 			HANDLE_ERROR(cudaMemcpy(gpuVote, cpuVote, bytes, cudaMemcpyHostToDevice));
 			//call the GPU version of the local max function
-			gpu_local_max<T>(gpuCenters, gpuVote, conn, x, y);
+			gpu_local_max<T>(gpuCenters, gpuVote, final_t, conn, x, y);
 			//copy the cell centers data to the CPU
 			cudaMemcpy(cpuCenters, gpuCenters, bytes, cudaMemcpyDeviceToHost) ;
@@ -2,33 +2,47 @@
 #define STIM_GL_SPIDER_H
 //#include <GL/glew.h>
+///basic GL and Cuda libs
 #include <GL/glut.h>
 #include <cuda.h>
 #include <cuda_gl_interop.h>
 #include <cudaGL.h>
-#include <math.h>
+
+///stim .*h for gl tracking and visualization
 #include <stim/gl/gl_texture.h>
-#include <stim/visualization/camera.h>
 #include <stim/gl/error.h>
+#include <stim/visualization/camera.h>
+#include <stim/math/matrix_sq.h>
+#include <stim/cuda/spider_cost.cuh>
+#include <stim/visualization/glObj.h>
+#include <stim/cuda/cudatools/glbind.h>
+
+///stim *.h for CUDA
+#include <stim/cuda/cuda_texture.cuh>
+#include <stim/cuda/cudatools.h>
+
+///stim *.h for branch detection
+#include <stim/visualization/cylinder.h>
+#include <stim/cuda/branch_detection.cuh>
+#include "../../../volume-spider/glnetwork.h"
+
+/// *.h for math
+#include <math.h>
 #include <stim/math/vector.h>
 #include <stim/math/vec3.h>
 #include <stim/math/rect.h>
-#include <stim/math/matrix.h>
 #include <stim/math/constants.h>
-#include <stim/cuda/spider_cost.cuh>
-#include <stim/cuda/cudatools/glbind.h>
 #include <stim/cuda/arraymath.cuh>
-#include <stim/cuda/cuda_texture.cuh>
-#include <stim/cuda/cudatools.h>
-#include <stim/cuda/ivote.cuh>
-#include <stim/visualization/glObj.h>
+#include <stim/math/random.h>
+
+
+///c++ libs for everything
 #include <vector>
 #include <stack>
-#include <stim/cuda/branch_detection.cuh>
-#include "../../../volume-spider/glnetwork.h"
-#include <stim/visualization/cylinder.h>
 #include <iostream>
 #include <fstream>
+
+///Timing and debugging
 #ifdef TIMING
 	#include <ctime>
 	#include <cstdio>	
@@ -39,6 +53,7 @@
 	#include <ctime>
 #endif
+//	#include <stim/cuda/testKernel.cuh>
 #ifdef DEBUG
 	#include <stim/cuda/testKernel.cuh>
 #endif
@@ -61,6 +76,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			double cost_time;// = 0;
 			double network_time;// = 0;
 			double hit_time;// = 0;
+			int iter_t;
 		#endif
 		stim::vec3<float> p;  	//vector designating the position of the spider.
@@ -72,7 +88,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		std::vector<float> mV;	//A list of all the size vectors.
 		std::vector<float> lV;	//A list of all the size vectors.
-		stim::matrix<float, 4> cT;			//current Transformation matrix (tissue)->(texture)
+		stim::matrix_sq<float, 4> cT;			//current Transformation matrix (tissue)->(texture)
 		GLuint texID;						//OpenGL ID for the texture to be traced
 		stim::vec3<float> S;					//Size of a voxel in the volume.
 		stim::vec3<float> R;					//Dimensions of the volume.
@@ -339,31 +355,6 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			#endif
 		}
-
-		float uniformRandom()
-		{
-			return (  (float)(rand()))/(  (float)(RAND_MAX));							///generates a random number between 0 and 1 using the uniform distribution.
-		}
-
-		float normalRandom()
-		{
-			float u1 = uniformRandom();
-			float u2 = uniformRandom();
-			return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));							///generate a random number using the normal distribution between 0 and 1.
-		}
-
-		stim::vec3<float> uniformRandVector()
-		{
-			stim::vec3<float> r(uniformRandom(), uniformRandom(), 1.0);						///generate a random vector using the uniform distribution between 0 and 1.
-			return r;
-		}
-
-		stim::vec3<float> normalRandVector()
-		{
-			stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);						///generate a random vector using the normal distribution between 0 and 1.
-			return r;
-		}
-
 //--------------------------------------------------------------------------//
 //---------------------TEMPLATE CREATION METHODS----------------------------//
@@ -400,8 +391,8 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			glNewList(dList, GL_COMPILE); 			///create a display list of all the direction templates.
 			for(int i = 0; i < numSamples; i++)		///for each sample
 			{
-				z = uniformRandom()*range + Z[1];	///generate a z coordinate
-				theta = uniformRandom()*stim::TAU;	///generate a theta coordinate
+				z = stim::Random<float>::uniformRandom()*range + Z[1];	///generate a z coordinate
+				theta = stim::Random<float>::uniformRandom()*stim::TAU;	///generate a theta coordinate
 				phi = acos(z);				///generate a phi from the z.
 				stim::vec3<float> sph(1, theta, phi);	///combine into a vector in spherical coordinates.
 				stim::vec3<float> cart = sph.sph2cart();///convert to cartesian.	
@@ -460,7 +451,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
                         UpdateBuffer(0.0, 0.0+0*n_pixels);
 			for(int i = 1; i < numSamplesPos; i++)		///for the number of position samples
 			{
-				stim::vec3<float> temp = uniformRandVector();	///generate a random point on a plane.
+				stim::vec3<float> temp = stim::Random<float>::uniformRandVector();	///generate a random point on a plane.
 				temp[0] = temp[0]*delta;
 				temp[1] = temp[1]*2*stim::PI;
@@ -953,6 +944,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			cost_time = 0;
 			network_time = 0;
 			hit_time = 0;
+			iter_t = 0;
 #endif
 #ifdef DEBUG
 			iter = 0;
@@ -960,7 +952,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			iter_dir = 0;
 			iter_siz = 0;
 #endif
-			stepsize = 6.0;
+			stepsize = 3.0;
 			n_pixels = 16.0;
 			srand(100);	
@@ -1274,12 +1266,19 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		{
 			std::vector <double> ret;
 			ret.resize(7);
+		//	ret[0] = branch_time/(float)iter;
 			ret[0] = branch_time;
+		//	ret[1] = direction_time/(float)iter;
 			ret[1] = direction_time;
+		//	ret[2] = position_time/(float)iter;
 			ret[2] = position_time;
+		//	ret[3] = size_time/(float)iter;
 			ret[3] = size_time;
+		//	ret[4] = cost_time/(float)iter;
 			ret[4] = cost_time;
+		//	ret[5] = network_time/(float)iter;
 			ret[5] = network_time;
+		//	ret[6] = hit_time/(float)iter;
 			ret[6] = hit_time;
 			return ret;
@@ -1311,6 +1310,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 					setSeed(x, y, z);
 					setSeedVec(u, v, w);
 					setSeedMag(m);
+					std::cout << getLastSeed() << getLastSeedVec() << std::cout;
                    		}
                    	myfile.close();
 	          	} else {
@@ -1320,7 +1320,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		///Saves the network to a file.
 		void
-		saveNetwork(std::string name)
+		saveNetwork(std::string name, int xoffset = 0, int yoffset = 0, int zoffset = 0)
 		{
 			stim::glObj<float> sk1;
 			for(int i = 0; i < nt.sizeE(); i++)
@@ -1331,7 +1331,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 				for(int j = 0; j < ce.size(); j++)
 				{
 					sk1.TexCoord(cm[j]);
-					sk1.Vertex(ce[j][0], ce[j][1], ce[j][2]);
+					sk1.Vertex(ce[j][0]+xoffset, ce[j][1]+yoffset, ce[j][2]+zoffset);
 				}
 				sk1.End();
 			}	
@@ -1408,6 +1408,11 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			#ifdef DEBUG
 			std::cerr << std::endl;
 			#endif
+			
+			#ifdef TIMING
+			iter_t++;
+			#endif
+
 			return current_cost;
 		}
@@ -1422,36 +1427,6 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 //--------------------------------------------------------------------------//
 //-----------------------------EXPERIMENTAL METHODS-------------------------//
 //--------------------------------------------------------------------------//
-
-		void
-		MonteCarloDirectionVectors(int nSamples, float solidAngle = stim::TAU)
-		{
-//			float PHI[2];//, Z[2];//, range;
-//			PHI[0] = asin(solidAngle/2);
-//			PHI[1] = asin(0);
-			
-//			Z[0] = cos(PHI[0]);
-//			Z[1] = cos(PHI[1]);
-			
-//			range = Z[0] - Z[1];
-
-			std::vector<stim::vec3<float> > vecsUni;
-			for(int i = 0; i < numSamplesPos; i++)
-			{
-				stim::vec3<float> a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0);
-				a[0] = a[0]-0.4;
-				a[1] = a[1]-0.4;
-				vecsUni.push_back(a);
-			}
-
-			stringstream name;
-			for(int i = 0; i < numSamplesPos; i++)
-				name << vecsUni[i].str() << std::endl;
-			
-			std::ofstream outFile;
-			outFile.open("New_Pos_Vectors.txt");
-			outFile << name.str().c_str();
-		}
 /*
 		void
 		DrawCylinder()
@@ -1570,9 +1545,6 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		{
 		//Define the varibles and turn on Selection Mode
-			#ifdef TIMING
-				gpuStartTimer();
-			#endif
 			GLuint selectBuf[2048];			///size of the selection buffer in bytes.
 			GLint hits;				///hit fibers
@@ -1602,6 +1574,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 				glPushMatrix();
 				glLoadIdentity();
+			#ifdef TIMING
+				gpuStartTimer();
+			#endif
 				CHECK_OPENGL_ERROR
 				gluLookAt(ps[0], ps[1], ps[2],
 					 ds[0], ds[1], ds[2],
@@ -1849,6 +1824,65 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 				std::cout << "I broke out!" << std::endl;
 				#endif
 		}
+//--------------------------------------------------------------------------//
+//--------------------------------DEPRECIATED-------------------------------//
+//--------------------------------------------------------------------------//
+/*
+		float uniformRandom()
+		{
+			return (  (float)(rand()))/(  (float)(RAND_MAX));							///generates a random number between 0 and 1 using the uniform distribution.
+		}
+
+		float normalRandom()
+		{
+			float u1 = uniformRandom();
+			float u2 = uniformRandom();
+			return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));							///generate a random number using the normal distribution between 0 and 1.
+		}
+
+		stim::vec3<float> uniformRandVector()
+		{
+			stim::vec3<float> r(uniformRandom(), uniformRandom(), 1.0);						///generate a random vector using the uniform distribution between 0 and 1.
+			return r;
+		}
+
+		stim::vec3<float> normalRandVector()
+		{
+			stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);						///generate a random vector using the normal distribution between 0 and 1.
+			return r;
+		}
+
+		void
+		MonteCarloDirectionVectors(int nSamples, float solidAngle = stim::TAU)
+		{
+//			float PHI[2];//, Z[2];//, range;
+//			PHI[0] = asin(solidAngle/2);
+//			PHI[1] = asin(0);
+			
+//			Z[0] = cos(PHI[0]);
+//			Z[1] = cos(PHI[1]);
+			
+//			range = Z[0] - Z[1];
+
+			std::vector<stim::vec3<float> > vecsUni;
+			for(int i = 0; i < numSamplesPos; i++)
+			{
+				stim::vec3<float> a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0);
+				a[0] = a[0]-0.4;
+				a[1] = a[1]-0.4;
+				vecsUni.push_back(a);
+			}
+
+			stringstream name;
+			for(int i = 0; i < numSamplesPos; i++)
+				name << vecsUni[i].str() << std::endl;
+			
+			std::ofstream outFile;
+			outFile.open("New_Pos_Vectors.txt");
+			outFile << name.str().c_str();
+		}
+
+*/
 };
 }
 #endif
@@ -187,7 +187,7 @@ namespace stim {
 		}
 		bool open(std::string filename) {								//open the bitmap file and read the header data
-			file.open(filename, std::ifstream::binary);
+			file.open(filename.c_str(), std::ifstream::binary);
 			if (!file) {
 				std::cout << "stim::bmp ERROR: error opening file: " << filename.c_str() << std::endl;
 				return false;
@@ -246,7 +246,7 @@ namespace stim {
 		info_header.bcSize = sizeof(tagBITMAPCOREHEADER);
 		info_header.bcPlanes = 1;
-		std::ofstream outfile(filename, std::ios::binary);										//open the output file for binary writing
+		std::ofstream outfile(filename.c_str(), std::ios::binary);										//open the output file for binary writing
 		outfile.write((char*)&file_header, sizeof(tagBITMAPFILEHEADER));						//write the file header
 		outfile.write((char*)&info_header, sizeof(tagBITMAPCOREHEADER));						//write the information header
@@ -259,4 +259,4 @@ namespace stim {
 		free(pad);
 		return true;
 	}
-}
 \ No newline at end of file
+}
@@ -606,7 +606,7 @@ public:
 	//crop regions given by an array of 1D index values
 	std::vector< image<T> > crop_idx(size_t w, size_t h, std::vector<size_t> idx) {
-		std::vector<image<T>> result(idx.size());										//create an array of image files to return
+		std::vector<image<T> > result(idx.size());										//create an array of image files to return
 		for (size_t i = 0; i < idx.size(); i++) {										//for each specified index point
 			size_t y = idx[i] / X();													//calculate the y coordinate from the 1D index (center of ROI)
 			size_t x = idx[i] - y * X();												//calculate the x coordinate (center of ROI)
@@ -18,14 +18,13 @@ class circle : plane&lt;T&gt;
 private:
-	//stim::vec3<T> Y;
-	T R;								//radius of the circle
+	stim::vec3<T> Y;
-	/*CUDA_CALLABLE void
+	CUDA_CALLABLE void
 	init()
 	{
 		Y = U.cross(N).norm();
-	}*/
+	}
 public:
 	using stim::plane<T>::n;
@@ -35,8 +34,6 @@ public:
 	using stim::plane<T>::rotate;
 	using stim::plane<T>::setU;
-	using stim::plane<T>::init;
-
 	///base constructor
 	///@param th value of the angle of the starting point from 0 to 360.
 	CUDA_CALLABLE
@@ -45,28 +42,26 @@ public:
 		init();
 	}
-	///create a circle given a size and position in Z space.
+	///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 radius, T z_pos = (T)0) : plane<T>(z_pos)
+	circle(T size, T z_pos = (T)0) : plane<T>()
 	{
-		//center(stim::vec3<T>(0, 0, z_pos));
-		//scale(size);
-		//init();
-		R = radius;
+		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>(n, c)
+	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
 	{
-		//center(c);
-		//normal(n);
-		//init();
-		R = (T)1;
+		center(c);
+		normal(n);
+		init();
 	}
 	/*///create a rectangle from a center point, normal, and size
@@ -89,18 +84,14 @@ 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 r, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
+	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
 	{
-		P = c;
-		N = n;
+		init();
 		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
@@ -108,111 +99,86 @@ public:
 	CUDA_CALLABLE
 	void scale(T factor)
 	{
-		//U *= factor;
-		//Y *= factor;
-		R *= factor;
-	}
-
-	///set the radius of circle to a certain value 
-	///@param value: the new radius of the circle
-	CUDA_CALLABLE
-	void set_R(T value) 
-	{
-		R = value;
+		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);
+	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){
-		P = p;
+		this->P = p;
 	}
 	///boolean comparison
 	bool
 	operator==(const circle<T> & rhs)
 	{
-		if(P == rhs.P && U == rhs.U)
+		if(P == rhs.P && U == rhs.U && Y == rhs.Y)
 			return true;
 		else
 			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 r, T theta)
+	CUDA_CALLABLE stim::vec3<T> operator()(T a, T 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);
+		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> > 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
+	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
-	///the only difference between points and glpoints is that the first point appears twice in the returning lists
-	std::vector< stim::vec3<T> > glpoints(unsigned n) {
-		std::vector< stim::vec3<T> > result(n + 1);
-		float dt = stim::TAU / n;
-		for (unsigned i = 0; i < n; i++)
-			result[i] = p(i * dt);
-		result[n] = p(0);									//close 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 << "r = "<<R<<"  (P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ")";
+		ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str();
 		return ss.str();
 	}
@@ -17,7 +17,7 @@ protected:
 		if(seed <= 0)
 			srand(time(NULL));
 		else if(seed > 0)
-			srand(time(seed));
+			srand(seed);
 		else
 			std::cout << "Error: Unknown value: in STIM::RANDOM"  << std::endl;
 	}
@@ -181,7 +181,7 @@ namespace stim {
 		/// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
 		/// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1.
 		/// @param 
-		/*void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>())
+		void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>())
 		{
 			std::vector<double> weights;		///the weight at each point on the surface of the sphere.
 												//		weights.resize(n);
@@ -223,7 +223,7 @@ namespace stim {
 				}
 			}
 			mcEnd();
-		}*/
+		}
 		std::string str() {
@@ -313,6 +313,52 @@ namespace stim {
 		T operator()(T theta, T phi) {
 			return p(theta, phi);			
 		}
+
+		//overload arithmetic operations
+
+		spharmonics<T> operator*(T rhs) const {
+
+			spharmonics<T> result(C.size());	//create a new spherical harmonics object
+
+			for (size_t c = 0; c < C.size(); c++)	//for each coefficient
+
+				result.C[c] = C[c] * rhs;	//calculate the factor and store the result in the new spharmonics object
+
+			return result;
+
+		}
+
+
+
+		spharmonics<T> operator+(spharmonics<T> rhs) {
+
+			size_t low = std::min(C.size(), rhs.C.size());		//store the number of coefficients in the lowest object
+			size_t high = std::max(C.size(), rhs.C.size());		//store the number of coefficients in the result
+			bool rhs_lowest = false;				//true if rhs has the lowest number of coefficients
+			if (rhs.C.size() < C.size()) rhs_lowest = true;		//if rhs has a lower number of coefficients, set the flag
+
+
+
+			spharmonics<T> result(high);								//create a new object
+
+			size_t c;
+			for (c = 0; c < low; c++)		//perform the first batch of additions
+				result.C[c] = C[c] + rhs.C[c];	//perform the addition
+
+			for (c = low; c < high; c++) {
+				if (rhs_lowest)
+					result.C[c] = C[c];
+				else
+					result.C[c] = rhs.C[c];
+			}
+			return result;
+		}
+
+
+
+		spharmonics<T> operator-(spharmonics<T> rhs) {
+			return (*this) + (rhs * (T)(-1));
+		}
 		/// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi]
 		void get_func(T* data, size_t X, size_t Y) {
 			T dt = stim::TAU / (T)X;			//calculate the step size in each direction
@@ -351,4 +397,4 @@ namespace stim {
 }
-#endif
 \ No newline at end of file
+#endif
@@ -3,268 +3,56 @@
 #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> {
-protected:
-	
-	using stim::centerline<T>::d;
+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> > U;					//stores the array of U vectors defining the Frenet frame
-	std::vector< T > R;				//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
+		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>::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
-		R.resize(size());
-	
-		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
-	}
-
-public:
-
-	using stim::centerline<T>::size;
-	using stim::centerline<T>::at;
-	using stim::centerline<T>::L;
-	using stim::centerline<T>::length;
-
-	cylinder() : centerline<T>(){}
-
-	cylinder(centerline<T>c) : centerline<T>(c) {
-		init();
-	}
-
-	//cylinder(centerline<T>c, T r) : centerline(c) {
-	//	init();
-	//	//add_mag(r);
-	//}
-
-	//initialize a cylinder with a list of points and magnitude values
-	//cylinder(centerline<T>c, std::vector<T> r) : centerline(c){
-	//	init();
-	//	//add_mag(r);
-	//}
-
-	//copy the original radius
-	void copy_r(std::vector<T> radius) {
-		for (unsigned i = 0; i < radius.size(); i++)
-			R[i] = radius[i];
-	}
-
-	///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 r(T l, int idx) {
-		T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
-		T v2 = (R[idx] + (R[idx + 1] - R[idx])*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 rl(T pvalue) {
-		if (pvalue <= 0.0) return R[0];
-		if (pvalue >= 1.0) return R[size() - 1];
-
-		T l = pvalue*L[L.size() - 1];
-		int idx = findIdx(l);
-		return r(l, idx);
-	}
-
-	///	Returns the magnitude at the given index
-	///	@param i is the index of the desired point
-	/// @param r is the index of the magnitude value
-	T r(unsigned i) {
-		return R[i];
-	}
-
-	///adds a magnitude to each point in the cylinder
-	/*void add_mag(V 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
-			R[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
-			R[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_r(size_t i, T r) {
-		R[i] = r;
-	}
-
-	/*size_t nmags() {
-		if (M.size() == 0) return 0;
-		else return R[0].size();
-	}*/
-
-	///Returns a circle representing the cylinder cross section at point i
-	stim::circle<T> circ(size_t i) {
-		return stim::circle<T>(at(i), R[i], 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) {
-		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 = " << R[i] << "   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() {
-		T sum = 0;							//initialize the integral to zero
-		if (L.size() == 1)
-			return sum;
-		else {					
-			T m0, m1;						//allocate space for both magnitudes in a single segment
-			m0 = R[0];					//initialize the first point and magnitude to the first point in the cylinder
-			T len = L[1];					//allocate space for the segment length
+		using stim::centerline<T>::c;
+		using stim::centerline<T>::N;
+	
+		///default init	
+		void
+		init()
+		{
-			for (unsigned p = 1; p < size(); p++) {				//for every consecutive point in the cylinder
-				m1 = R[p];
-				if (p > 1) len = (L[p] - L[p - 1]);		//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.R[i] = r(t);								//retrieve the interpolated magnitude from the current cylinder and store it in the new one
-			//}
-		}
-		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].R.resize(idx);
-		for (; i < idx + 1; i++) {
-			//for(unsigned d = 0; d < 3; d++)
-			//C[0][i][d] = LL[0].c[i][d];
-			C[0].R[i] = R[i];
-			//C[0].R[i].resize(1);
-		}
-		if (C.size() == 2) {
-			C[1] = LL[1];
-			i--;
-			//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].R[i - idx] = R[i];
-				//C[1].M[i - idx].resize(1);
-			}
 		}
-		return C;
-	}
-
+		///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;
+			}
+*/		}
-		/*
-		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM)
+		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
 		void
-		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
+		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());
@@ -286,7 +74,7 @@ public:
 			}
 			stim::vec3<T> dr = (inP[1] - inP[0]).norm();
-			s = stim::circle<T>(inP[0], inR[0][0], dr, stim::vec3<T>(1,0,0));
+			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;
@@ -300,7 +88,7 @@ public:
 				norms[i] = s.N;
-				s.scale(inR[i][0]/inR[i-1][0]);
+				s.scale(inM[i][0]/inM[i-1][0]);
 				e[i] = s;
 				Us[i] = e[i].U;
 			}
@@ -312,7 +100,7 @@ public:
 			norms[j] = s.N;
-			s.scale(inR[j][0]/inR[j-1][0]);
+			s.scale(inM[j][0]/inM[j-1][0]);
 			e[j] = s;
 			Us[j] = e[j].U;
 		}
@@ -439,7 +227,7 @@ public:
 			stim::vec<T> zero(0.0,0.0);
 			temp.resize(inM.size(), zero);
 			for(int i = 0; i < inM.size(); i++)
-				temp[i][0] = inR[i];
+				temp[i][0] = inM[i];
 			init(inP, temp);
 		}
@@ -458,18 +246,13 @@ public:
 			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).
@@ -483,7 +266,22 @@ public:
 			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.
@@ -523,7 +321,10 @@ public:
 			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.
@@ -705,10 +506,10 @@ public:
 			return cylinder<T>(result);
-		}*/
+		}
 };
 }
-#endif
 \ No newline at end of file
+#endif
@@ -51,6 +51,13 @@ namespace stim {
 			magnitude = true;
 		}
+		gl_spharmonics(stim::spharmonics<T> disp, stim::spharmonics<T> color, size_t slices) :
+			 gl_spharmonics<T>(slices)
+		{
+			Sc = color;
+			Sd = disp;
+		}
+
 		~gl_spharmonics() {
 			if (dlist) glDeleteLists(dlist, 1);		//delete the display list when the object is destroyed
 		}
@@ -166,13 +173,13 @@ namespace stim {
 		}
 		/// Project a set of samples onto the basis
-		void project(std::vector<vec3<float>>& vlist, size_t nc) {
+		void project(std::vector<vec3<float> >& vlist, size_t nc) {
 			Sd.project(vlist, nc);
 			Sc = Sd;
 		}
 		/// Calculate a density function from a list of points in spherical coordinates
-		void pdf(std::vector<stim::vec3<T>>& vlist, size_t nc) {
+		void pdf(std::vector<stim::vec3<T> >& vlist, size_t nc) {
 			Sd.pdf(vlist, nc);
 			Sc = Sd;
 		}
@@ -196,4 +203,4 @@ namespace stim {
-#endif
 \ No newline at end of file
+#endif