Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

David Mayerich
2 parents d11e2fb6 fac43319
Showing 9 changed files with 685 additions and 748 deletions Show diff stats
stim/biomodels/fiber.h → stim/biomodels/centerline.h
stim/cuda/testKernel.cuh
stim/gl/gl_spider.h
stim/gl/gl_texture.h
stim/math/circle.h
stim/math/vec3.h
stim/visualization/cylinder.h
stim/visualization/glObj.h
stim/visualization/obj.h
-#ifndef STIM_FIBER_H
-#define STIM_FIBER_H
+#ifndef STIM_CENTERLINE_H
+#define STIM_CENTERLINE_H
  
 #include <vector>
-#include <ANN/ANN.h>
+#include <stim/math/vec3.h>
+//#include <ANN/ANN.h>
  
 namespace stim{
  
@@ -11,22 +12,19 @@ namespace stim{
  *	class to describe an interconnected (often biological) network.
  */
 template<typename T>
-class fiber{
+class centerline{
  
 protected:
 	unsigned int N;					//number of points in the fiber
 	double **c;						//centerline (array of double pointers)
-
-	T* r;						// array of fiber radii
-	ANNkd_tree* kdt;			//kd-tree stores all points in the fiber for fast searching
+//	ANNkd_tree* kdt;				//kd-tree stores all points in the fiber for fast searching
  
 	/// Initialize an empty fiber
 	void init()
 	{
-		kdt = NULL;
-		c=NULL;
-		r=NULL;
 		N=0;
+		c=NULL;
+//		kdt = NULL;
 	}
  
 	/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
@@ -34,19 +32,17 @@ protected:
 	{
  
 		N = n;												//set the number of points
-		kdt = NULL;
+//		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);
-
-		r = (T*) malloc(sizeof(T) * N);			//allocate space for the radii
 	}
  
 	/// Copies an existing fiber to the current fiber
  
 	/// @param cpy stores the new copy of the fiber
-	void copy( const stim::fiber<T>& cpy ){
+	void copy( const stim::centerline<T>& cpy, bool kd = 0){
  
 		///allocate space for the new fiber
 		init(cpy.N);
@@ -55,20 +51,18 @@ protected:
 		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
-
-			r[i] = cpy.r[i];						//copy the radius
 		}
-
-		gen_kdtree();							//generate the kd tree for the new fiber
+//		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
-	}
+//	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){
@@ -79,7 +73,6 @@ protected:
 		{
 			v = p1[d] - p0[d];
 			sum +=v * v;
-
 		}
 		return sqrt(sum);
 	}
@@ -87,21 +80,18 @@ protected:
 	/// 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;
-	}
+//	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){
-		stim::vec<T> r;
+		stim::vec3<T> r;
 		r.resize(3);
 		r[0] = c[i][0];
 		r[1] = c[i][1];
@@ -113,25 +103,23 @@ protected:
  
 public:
  
-	fiber(){
+	centerline(){
 		init();
 	}
  
 	/// Copy constructor
-	fiber(const stim::fiber<T> &obj){
-
+	centerline(const stim::centerline<T> &obj){
 		copy(obj);
-
 	}
  
 	//temp constructor for graph visualization
-	fiber(int n)
+	centerline(int n)
 	{
 		init(n);
 	}
  
 	/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
-	fiber(std::vector< stim::vec<T> > p){
+	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
@@ -142,200 +130,77 @@ public:
 				c[i][d] = (double) p[i][d];
  
 			//set the radius
-			r[i] = 0;
 		}
-
 		//generate a kd tree
-		gen_kdtree();
+//		if(kd)
+//			gen_kdtree();
 	}
  
-	/// constructor takes a list of points and radii
-	fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){
+	/// 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
-			r[i] = radii[i];
 		}
  
 		//generate a kd tree
-		gen_kdtree();
-	}
-
-	/// constructor takes an array of points and radii
-	//		this function is used when the radii are represented as a stim::vec,
-	//		since this may be easier when importing OBJs
-	fiber(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii){
-
-		init(pos.size());
-
-		//for each point, set the position and radius
-		for(unsigned int i = 0; i < N; i++){
-			//at(i) = (double*)malloc(sizeof(double) * 3);
-			for(unsigned int d = 0; d < 3; d++)
-				c[i][d] = (double)  pos[i][d];
-
-			r[i] = radii[i][(unsigned int)0];
-		}
-
-		gen_kdtree();
+		//if(kd)
+		//	gen_kdtree();
 	}
  
 	/// Assignment operation
-	fiber& operator=(const fiber &rhs){
-
+	centerline& operator=(const centerline &rhs){
 		if(this == &rhs) return *this;			//test for and handle self-assignment
-
 		copy(rhs);
+		return *this;
 	}
  
-	/// Calculate the length of the fiber and return it.
-	double length(){
-
-		double* p0;
-		double *p1;
-		double l = 0;				//initialize the length to zero
-
-		//for each point
-		//typename std::list< point<T> >::iterator i;	//create a point iterator
-		for(unsigned int i = 0; i < N; i++){		//for each point in the fiber
-
-			if(i == 0)						//if this is the first point, just store it
-				p1 = c[0];
-			else{									//if this is any other point
-				p0 = p1;							//shift p1->p0
-				p1 = c[i];							//set p1 to the new point
-				l += dist(p0, p1);				//add the length of p1 - p0 to the running sum
-			}
-		}
-
-		return l;									//return the length
-	}
-
-	/// Calculates the length and average radius of the fiber
-
-	/// @param length is filled with the fiber length
-	T radius(T& length){
-
-		double* p0;				//temporary variables to store point positions
-		double* p1;
-		T r0, r1;					//temporary variables to store radii at points
-		double l;
-		T r_mean;						//temporary variable to store the length and average radius of a fiber segment
-		double length_sum = 0;			//initialize the length to zero
-		T radius_sum = 0;			//initialize the radius sum to zero
-
-		//for each point
-		//typename std::list< point<T> >::iterator i;	//create a point iterator
-		for(unsigned int i = 0; i < N; i++){		//for each point in the fiber
-
-			if(i == 0){						//if this is the first point, just store it
-				p1 = c[0];
-				r1 = r[0];
-			}
-			else{									//if this is any other point
-				p0 = p1;							//shift p1->p0 and r1->r0
-				r0 = r1;
-				p1 = c[i];							//set p1 to the new point
-				r1 = r[i];
-
-				l = dist(p0, p1);				//calculate the length of the p0-p1 segment
-				r_mean = (r0 + r1) / 2;					//calculate the average radius of the segment
-
-				radius_sum += r_mean * (T) l;				//add the radius scaled by the length to a running sum
-				length_sum += l;					//add the length of p1 - p0 to the running sum
-			}
-		}
-
-		length = length_sum;						//store the total length
-
-		//if the total length is zero, store a radius of zero
-		if(length == 0)
-			return 0;
-		else
-			return (T)(radius_sum / length);					//return the average radius of the fiber
-	}
-	T average_radius()
-	{
-		T r_sum = 0.;
-		for(unsigned int i = 0; i < N; i++)
-		{
-			r_sum = r_sum + r[i];
-		}
-		return r_sum/((T) N);
-	}
-
-	/// Calculates the average radius of the fiber
-	T radius(){
-		T length;
-		return radius(length);
-	}
-
-	/// Returns the radius at index idx.
-	T radius(int idx){
-		return r[idx];
-	}
  
 	/// 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
-	}
+//	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
-	}
+//	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> > centerline(){
+	std::vector< stim::vec<T> > get_centerline(){
  
 		//create an array of stim vectors
-		std::vector< stim::vec<T> > pts(N);
+		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::vec<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
-
-		//return the centerline array
-		return pts;
-	}
-
-	/// Returns the fiber centerline magnitudes as an array of stim::vec points
-	std::vector< stim::vec<T> > centerlinemag(){
-
-		//create an array of stim vectors
-		std::vector< stim::vec<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::vec<T>(r[i], r[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::fiber<T> > split(unsigned int idx){
+	std::vector< stim::centerline<T> > split(unsigned int idx){
  
-		std::vector< stim::fiber<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
  
 		//if the index is an end point, only the existing fiber is returned
 		if(idx == 0 || idx == N-1){
@@ -361,14 +226,13 @@ public:
 			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
-				fl[0].r[i] = r[i];						//copy the corresponding radius
 			}
  
 			//second half
 			for(i = 0; i < N2; i++){
 				for(d = 0; d < 3; d++)
 					fl[1].c[i][d] = c[idx + i][d];
-				fl[1].r[i] = r[idx + i];
+
 			}
 		}
  
@@ -379,7 +243,7 @@ public:
 	/// 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::fiber<T> > connect( stim::fiber<T> &f, double dist){
+/*	std::vector< stim::centerline<T> > connect( stim::centerline<T> &f, double dist){
  
 		double min_dist;
 		unsigned int idx0, idx1;
@@ -393,7 +257,7 @@ public:
  
  
 	}
-
+*/
 	/// Outputs the fiber as a string
 	std::string str(){
 		std::stringstream ss;
@@ -405,7 +269,6 @@ public:
 			for(unsigned int d = 0; d < 3; d++){
 				ss<<c[i][d]<<"  ";
 			}
-			ss<<"]  r = "<<r[i]<<std::endl;
 		}
  
 		return ss.str();
@@ -428,7 +291,7 @@ public:
 		return get_vec(N-1);
 	}
 		////resample a fiber in the network
-	stim::fiber<T> resample(T spacing)
+	stim::centerline<T> resample(T spacing)
 	{
 		std::cout<<"fiber::resample()"<<std::endl;
  
@@ -468,7 +331,7 @@ public:
 				newPointList.push_back(fiberPositions[f+1]);
 			}
 		newPointList.push_back(fiberPositions[N-1]);   //add the last point on the fiber to the new fiber list
-		fiber newFiber(newPointList);
+		centerline newFiber(newPointList);
 		return newFiber;
 	}
  
@@ -8,7 +8,7 @@
 #include <stim/cuda/cudatools/devices.h>
 #include <stim/cuda/cudatools/threads.h>
 #include <stim/cuda/cuda_texture.cuh>
-	stim::cuda::cuda_texture tx;  //texture object.
+	
 	float* print;
  
 	///Initialization function, allocates the memory and passes the necessary
@@ -26,16 +26,16 @@
 			cudaFree(print);         ///temporary
 	}  
  
-        __device__
-        float templ(int x)                                                                      
-        {
-                if(x < 32/6 || x > 32*5/6 || (x > 32*2/6 && x < 32*4/6)){
-                        return 1.0;
-                }else{
-                        return 0.0;
-                }
-
-        }
+	 __device__
+	float templ(int x, int max_x)
+	{
+	if(x < max_x/6 || x > max_x*5/6 || (x > max_x*2/6 && x < max_x*4/6))                                                                                                                             
+        	{
+        	        return 1.0;
+        	}else{
+        	        return 0.0;
+        	}
+	}
  
 	///Find the difference of the given set of samples and the template
 	///using cuda acceleration.
@@ -44,33 +44,24 @@
 	///@param float* result			--a pointer to the memory that stores the result.
 	__global__
 	//void get_diff (float *result)
-	void get_diff (cudaTextureObject_t texIn, float *print)
+	void get_diff (cudaTextureObject_t texIn, float *print, int dx)
 	{       
 		int x   = threadIdx.x + blockIdx.x * blockDim.x;
 		int y   = threadIdx.y + blockIdx.y * blockDim.y;
-//		int idx = y*64+x;
-		int idx = y*32+x;
+		int idx = y*dx+x;
 	//	int idx = y*16+x;
  
 		float valIn             = tex2D<unsigned char>(texIn, x, y);
-		float templa		= templ(x);
-		print[idx]              = valIn;             ///temporary
+		float templa		= templ(x, 32)*255.0;
+		print[idx]              = abs(valIn-templa);             ///temporary
 		//print[idx]              = abs(templa);             ///temporary
  
 	}
  
-
-	///External access-point to the cuda function
-	///@param GLuint texbufferID 	--GLtexture (most be contained in a framebuffer object)
-	///				  that holds the data that will be handed to cuda.
-	///@param GLenum texType	--either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D
-	///				  may work with other gl texture types, but untested.
-	///@param DIM_Y, the number of samples in the template.
-	void test(GLint texbufferID, GLenum texType, int x, int y)
+	void test(cudaTextureObject_t tObj, int x, int y, std::string nam)
 	{
  
 		//Bind the Texture in GL and allow access to cuda.
-		tx.MapCudaTexture(texbufferID, texType);
  
 		//initialize the return arrays.
  
@@ -85,44 +76,13 @@
  
  
 //		get_diff <<< blocks, threads >>> (tx.getTexture(), print);
-		get_diff <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), print);
-
-		cudaDeviceSynchronize();
-		stringstream name;      //for debugging
-		name << "FromTex.bmp";
-		stim::gpu2image<float>(print, name.str(),x,y,0,255);
-	  
-		tx.UnmapCudaTexture();
-		cleanUP();
-	}
-
-	void test(GLint texbufferID, GLenum texType, int x, int y, std::string nam)
-	{
-
-		//Bind the Texture in GL and allow access to cuda.
-		tx.MapCudaTexture(texbufferID, texType);
-
-		//initialize the return arrays.
-
-		initArray(x,y);
-		dim3 numBlocks(1, y);
-		dim3 threadsPerBlock(x, 1);
-		int max_threads = stim::maxThreadsPerBlock();
-		//dim3 threads(max_threads, 1);
-		//dim3 blocks(x / threads.x + 1, y);	
-		//dim3 numBlocks(2, 2);
-		//dim3 threadsPerBlock(8, 108);
-
-
-//		get_diff <<< blocks, threads >>> (tx.getTexture(), print);
-		get_diff <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), print);
+		get_diff <<< numBlocks, threadsPerBlock >>> (tObj, print, x);
  
 		cudaDeviceSynchronize();
 		stringstream name;      //for debugging
 		name << nam.c_str();
-		//stim::gpu2image<float>(print, name.str(),x,y,0,255);
+		stim::gpu2image<float>(print, name.str(),x,y,0,255);
  
-		tx.UnmapCudaTexture();
 		cleanUP();
 	}
  
@@ -23,14 +23,11 @@
 #include <stim/cuda/ivote.cuh>
 #include <stim/visualization/glObj.h>
 #include <vector>
+#include <stack>
 #include <stim/cuda/branch_detection.cuh>
-#include "../../../volume-spider/fiber.h"
 #include "../../../volume-spider/glnetwork.h"
 #include <stim/visualization/cylinder.h>
 #include <stim/cuda/testKernel.cuh>
-
-//#include <stim/cuda/testKernel.cuh>
-
 #include <iostream>
 #include <fstream>
 #ifdef TIMING
@@ -39,7 +36,6 @@
 #endif
  
 #ifdef TESTING
-	#include <iostream>
 	#include <cstdio>
 	#include <ctime>
 #endif
@@ -49,7 +45,7 @@ namespace stim
 {
  
 template<typename T>
-class gl_spider : public virtual gl_texture<T>
+class gl_spider // : public virtual gl_texture<T>
 {
 	//doen't use gl_texture really, just needs the GLuint id.
 	//doesn't even need the texture iD really.
@@ -65,72 +61,87 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			double hit_time;// = 0;
 		#endif
  
-		//
 		stim::vec3<float> p;  	//vector designating the position of the spider.
-		stim::vec3<float> d;	//vector designating the orientation of the spider
-						//always a unit vector.
-		stim::vec<float> m;	//magnitude of the spider vector.
-						//mag[0] = length.
-						//mag[1] = width.
+		stim::vec3<float> d;	//normalized direction of travel
+		float m;			//size of the spider in tissue space.
+
 		std::vector<stim::vec3<float> > dV;	//A list of all the direction vectors.
-		std::vector<stim::vec3<float> > pV;	//A list of all the position vectors.
-		std::vector<stim::vec3<float> > mV;	//A list of all the size vectors.
+		std::vector<stim::vec3<float> > pV;	//A list of all test positions (relative to p)
+		std::vector<float> mV;	//A list of all the size vectors.
  
-		stim::matrix<float, 4> cT;		//current Transformation matrix
-							//From tissue space to texture space.
-		GLuint texID;
-		stim::vec<float> S;			//Size of a voxel in the volume.
-		stim::vec<float> R;			//Dimensions of the volume.
+		stim::matrix<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.
  
  
 		//GL and Cuda variables
-		GLuint dList;				//displaylist ID
-		GLuint fboID;				//framebuffer ID
-		GLuint texbufferID;			//texbuffer ID, only necessary for 
-							//cuda aspect of the calculation.
-		GLuint pfboID;				//buffer object for position tracking.
-		GLuint ptexbufferID;			//texture object for position tracking.
-
-		GLuint mfboID;				//buffer object for magnitude adjustment.
-		GLuint mtexbufferID;			//texture object for magnitude adjustment.
-		GLuint bfboID;				//buffer object for position adjustment.
-		GLuint btexbufferID;			//buffer object for position adjustment.
+		GLuint dList;				//ID of the starting display lists (series of 4)
+									//dList + 0 = direction template rectangles
+									//dList + 1 = position template rectangles
+									//dList + 2 = size template rectangles
+									//dList + 3 = branch detection cylinder around the fiber
+
+		GLuint fboID;				//framebuffer ID for direction templates
+		GLuint texbufferID;			//texture ID for direction templates
+		GLuint direction_buffID;				//framebuffer ID, position templates
+		GLuint direction_texID;		//texture ID, position templates
+
+		GLuint position_buffID;				//framebuffer ID, position templates
+		GLuint position_texID;		//texture ID, position templates
+
+		GLuint radius_buffID;				//framebuffer ID, radius templates
+		GLuint radius_texID;		//texture ID, radius templates
+
+		GLuint cylinder_buffID;				//framebuffer ID, cylinder (surrounding fiber)
+		GLuint cylinder_texID;		//texture ID, cylinder
  
 		int numSamples;				//The number of templates in the buffer.
 		int numSamplesPos;
 		int numSamplesMag;
  
-		float stepsize;// = 5.0;			//Step size.
-//		float stepsize = 3.0;			//Step size.
-		int current_cost;			//variable to store the cost of the current step.
-		
+		float length;				//this will be a function of the radius
+		float stepsize;				//this will be a function of the length
+
+		int current_cost;			//variable to store the cost of the current step		
  
 		//Tracing variables.
-		std::stack< stim::vec3<float> > seeds;		//seed positions.
-		std::stack< stim::vec3<float> > seedsvecs;	//seed directions.
-		std::stack< float > 		seedsmags;	//seed magnitudes.
+		std::stack< stim::vec3<float> > seeds;		//seed positions
+		std::stack< stim::vec3<float> > seedsvecs;	//seed directions
+		std::stack< float > 		seedsmags;		//seed magnitudes
+
+		std::vector< stim::vec3<float> > cL;	//centerline up to the current point
+		std::vector< stim::vec3<float> > cD;	//directions up to the current point (debugging)
+		std::vector< float > cM;		//radius up to the current point
  
-		std::vector< stim::vec3<float> > cL;	//Positions of line currently being traced.
-		std::vector< stim::vec3<float> > cD;	//Direction of line currently being traced.
-		std::vector< stim::vec<float> > cM;	//Magnitude of line currently being traced.
+		stim::glnetwork<float> nt;				//network object holding the currently traced centerlines
+		stim::glObj<float> sk;					//OBJ file storing the network (identical to above)
  
-		stim::glnetwork<float> nt;		//object for storing the network.
+		//consider replacing with two seed points facing opposite directions
+		stim::vec<float> rev;					//reverse vector
  
-		stim::vec<float> rev;			//reverse vector;
-		stim::camera camSel;
-		stim::vec3<float> ps;
-		stim::vec3<float> ups;
-		stim::vec3<float> ds;
+		//selection mode - detecting fiber intersections
+		stim::camera camSel;					//camera for selection mode (detecting collisions)
+		stim::vec3<float> ps;					//position for the selection camera
+		stim::vec3<float> ups;					//up direction for the selection camera
+		stim::vec3<float> ds;					//direction for the selection camera
  
-		//static const float t_length = 16.0;
-		float t_length;
+		float n_pixels;							//length of the template (in pixels)
  
 		//cuda texture variables that keep track of the binding.
-		stim::cuda::cuda_texture t_dir;
-		stim::cuda::cuda_texture t_pos;
-		stim::cuda::cuda_texture t_mag;
+		stim::cuda::cuda_texture t_dir;				//cuda_texture object used as an interface between OpenGL and cuda for direction vectors.
+		stim::cuda::cuda_texture t_pos;				//cuda_texture object used as an interface between OpenGL and cuda for position vectors.
+		stim::cuda::cuda_texture t_mag;				//cuda_texture object used as an interface between OpenGL and cuda for size vectors.
  
  
+		#ifdef DEBUG
+			stringstream name;
+			int iter;
+			int iter_pos;
+			int iter_dir;
+			int iter_siz;
+		#endif
+
 //--------------------------------------------------------------------------//
 //-------------------------------PRIVATE METHODS----------------------------//
 //--------------------------------------------------------------------------//
@@ -142,31 +153,32 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		findOptimalDirection()
 		{
 			#ifdef TIMING
-				gpuStartTimer();
+				gpuStartTimer();	//Timer for profiling
 			#endif
 			setMatrix();			//create the transformation matrix.
 			glCallList(dList);		//move the templates to p, d, m.
-			glFinish();
-//			glFlush();
+			glFinish();			//flush the pipeline
 			#ifdef TIMING
-				direction_time += gpuStopTimer();
-			#endif
-			#ifdef TESTING
-//				test(texbufferID, GL_TEXTURE_2D,2*t_length,numSamples*t_length, "Final_Cost_Direction.bmp");
+				direction_time += gpuStopTimer();	//profiling
 			#endif
  
 			int best = getCost(t_dir.getTexture(), t_dir.getAuxArray() ,numSamples);		//find min cost.
-			stim::vec<float> next(		//find next vector.
+			#ifdef DEBUG
+				name.str("");
+				name << "Final_Cost_Direction_fiber_"<< iter << "_" << iter_dir << ".bmp";	
+				test(t_dir.getTexture(), n_pixels*2.0, numSamples*n_pixels, name.str());
+				iter_dir++;
+			#endif
+			stim::vec<float> next(		///calculate the next vector.
 			dV[best][0]*S[0]*R[0],
 			dV[best][1]*S[1]*R[1],
 			dV[best][2]*S[2]*R[2],
 			0);
-			next = (cT*next).norm();	//find next vector.
-			setPosition(	p[0]+next[0]*m[0]/stepsize,
-					p[1]+next[1]*m[0]/stepsize,
-					p[2]+next[2]*m[0]/stepsize);
-			setDirection(next[0], next[1], next[2]);
-							//move forward and change direction.
+			next = (cT*next).norm();	///transform the next vector into Tissue space.
+			setPosition(	p[0]+next[0]*m/stepsize,
+					p[1]+next[1]*m/stepsize,
+					p[2]+next[2]*m/stepsize);	
+			setDirection(next[0], next[1], next[2]);		//move forward and change direction.
 		}
  
 		/// Method for finding the best d (direction) for the spider.
@@ -176,27 +188,29 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		findOptimalPosition()
 		{
 			#ifdef TIMING
-				gpuStartTimer();
+				gpuStartTimer();	//timer for profiling
 			#endif
 			setMatrix();			//create the transformation matrix.
 			glCallList(dList+1);		//move the templates to p, d, m.
-			glFinish();
+			glFinish();			//flush the pipeline
 //			glFlush();
 			#ifdef TIMING
-				position_time += gpuStopTimer();
+				position_time += gpuStopTimer();		///timer for profiling
 			#endif
  
-			#ifdef TESTING
-//				test(ptexbufferID, GL_TEXTURE_2D,2*t_length, numSamplesPos*t_length, "Final_Cost_Position.bmp");
-			#endif
 			int best = getCost(t_pos.getTexture(), t_pos.getAuxArray(), numSamplesPos);		//find min cost.
-//			std::cerr << best << std::endl;
+			#ifdef DEBUG
+				name.str("");
+				name << "Final_Cost_Position_" << iter << "_" << iter_pos << ".bmp";	
+				test(t_pos.getTexture(), n_pixels*2.0, numSamplesPos*n_pixels, name.str());
+				iter_pos++;
+			#endif
 			stim::vec<float> next(		//find next position.
  				pV[best][0],
 				pV[best][1],
 				pV[best][2],
 				1);
-			next = cT*next;			//find next position.
+			next = cT*next;			//transform the next position vector into tissue space.
 			setPosition(
 					next[0]*S[0]*R[0],
 					next[1]*S[1]*R[1],
@@ -206,7 +220,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
  
 		/// Method for finding the best scale for the spider.
 		/// changes the x, y, z size of the spider to minimize the cost
-		/// function. */
+		/// function. 
 		void
 		findOptimalScale()
 		{
@@ -215,64 +229,21 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			#endif
 			setMatrix();			//create the transformation.
 			glCallList(dList+2);		//move the templates to p, d, m.
-			glFinish();
-//			glFlush();
+			glFinish();			//flush the drawing pipeline.
 			#ifdef TIMING
 				size_time += gpuStopTimer();
 			#endif
-			#ifdef TESTING
-//				test(mtexbufferID, GL_TEXTURE_2D, 2*t_length, numSamplesMag*t_length, "Final_Cost_Position.bmp");
-			#endif
 			int best = getCost(t_mag.getTexture(), t_mag.getAuxArray(), numSamplesMag);		//get best cost.
-			setMagnitude(m[0]*mV[best][0]);	//adjust the magnitude.
+			#ifdef DEBUG
+				name.str("");
+				name << "Final_Cost_Size_" << iter << "_" <<  iter_siz << ".bmp";	
+				test(t_mag.getTexture(), n_pixels*2.0, numSamplesMag*n_pixels, name.str());
+				iter_siz++;
+			#endif
+			setMagnitude(m*mV[best]);	//adjust the magnitude.
 		}
  
  
-		///subject to change.
-		///finds branches.
-		///depreciated
-		void
-		branchDetection()
-		{
-			setMatrix();
-			glCallList(dList+3);
-			std::vector< stim::vec<float> > result = find_branch(
-					btexbufferID, GL_TEXTURE_2D, 16, 216);
-			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
-			if(!result.empty())
-			{
-				for(int i = 1; i < result.size(); i++)
-				{
-					stim::vec<float> cylp(
-						0.5 * cos(2*M_PI*(result[i][1])),
-						0.5 * sin(2*M_PI*(result[i][1])),
-						result[i][0]-0.5,
-						1.0);
-					cylp = cT*cylp;
-					
-					stim::vec3<float> vec(
-						cylp[0]*S[0]*R[0],
-						cylp[1]*S[1]*R[1],
-						cylp[2]*S[2]*R[2]);
-					stim::vec3<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],
-						   -p[1] + cylp[1]*S[1]*R[1],
-						   -p[2] + cylp[2]*S[2]*R[2]);
-						seeddir = seeddir.norm();
-					float seedm = m[0];
-// Uncomment for global run 
-					if(
-					 !(vec[0] > size[0] || vec[1] > size[1]
-					 || vec[2] > size[2] || vec[0] < 0
-					 || vec[1] < 0 || vec[2] < 0))
-					{
-						setSeed(vec);
-						setSeedVec(seeddir);
-						setSeedMag(seedm);
-					}
-				}
-			}    
-			
-		}
  
  
 		///finds all the branches in the a given fiber.
@@ -281,78 +252,78 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		branchDetection2(int n = 8, int l_template = 8, int l_square = 8)
 		{
 			#ifdef TIMING
-				gpuStartTimer();
+				gpuStartTimer();		///timer for performance analysis
 			#endif
  
-			if(cL.size() < 4){}
+			if(cL.size() < 4){}		///if the size of the fiber is less then 4 we do nothing.
 			else{
-			setMatrix(1);
-			DrawLongCylinder(n, l_template, l_square);
+			setMatrix(1);			///finds the current transformation matrix
+			DrawLongCylinder(n, l_template, l_square);	///Draw the cylinder.
 			stim::cylinder<float> cyl(cL, cM);
-			std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
-			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
-			float pval;
-			if(!result.empty())
+			std::vector< stim::vec<float> > result = find_branch(cylinder_texID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);		///find all the centers in cuda
+			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);			///the borders of the texture.
+			float pval;									//pvalue associated with the points on the cylinder.
+			if(!result.empty())								///if we have any points
 			{
-				for(int i = 0; i < result.size(); i++)
+				for(int i = 0; i < result.size(); i++)					///for each point
 				{
 					int id = result[i][2];
-					if(fmod(result[i][2], id) != 0 && id != 0)
+					if(fmod(result[i][2], id) != 0 && id != 0)			///if the remainer is odd
 					{
  
 						pval = ((cyl.getl(id+1)-cyl.getl(id))*
-							(fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);
+							(fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);		///calculate pvalue
 					}
-					else if(id == 0)
+					else if(id == 0)									///if the point is on the edge
 					{
-						pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);
+						pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);			
 					}
 					else
 					{
-						pval = (cyl.getl(id)/cyl.getl(cL.size()-1));
+						pval = (cyl.getl(id)/cyl.getl(cL.size()-1));					///if the point is somewhere on the surface of the cylinder other than the edge
 					}
-					stim::vec3<float> v = cyl.surf(pval, result[i][0]);
-					stim::vec3<float> di = cyl.p(pval);
-					float rad = cyl.r(pval);
+					stim::vec3<float> v = cyl.surf(pval, result[i][0]);					///find the coordinates of the point at pval on the surface in tissue space.
+					stim::vec3<float> di = cyl.p(pval);							///find the coord of v in tissue space projected on the centerline.
+					float rad = cyl.r(pval)/2;								///find the radius at the pvalue's location
 					if(
 					 !(v[0] > size[0] || v[1] > size[1]
 					 || v[2] > size[2] || v[0] < 0
-					 || v[1] < 0 || v[2] < 0))
+					 || v[1] < 0 || v[2] < 0))								///if the v point is INSIDE the volume
 					{
-						setSeed(v);
-						setSeedVec((v-di).norm());
-						setSeedMag(rad);
+						setSeed(v);									///add a seedpoint's position.
+						setSeedVec((v-di).norm());							///add a seedpoints direction
+						setSeedMag(rad);								///add the starting radius.
 					}
 				}
 			}
 			}
 			#ifdef TIMING
-				branch_time += gpuStopTimer();
+				branch_time += gpuStopTimer();									///timer for performance.
 			#endif
 		}
  
  
 		float uniformRandom()
 		{
-			return (  (float)(rand()))/(  (float)(RAND_MAX));
+			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));
+			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);
+			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);
+			stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);						///generate a random vector using the normal distribution between 0 and 1.
 			return r;
 		}
  
@@ -370,49 +341,57 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///Stored in a display list.
 		///uses the default d vector <0,0,1>
 		void
-		genDirectionVectors(float solidAngle = M_PI/2)
+		genDirectionVectors(float solidAngle = stim::PI/2)
 		{
  
 			//Set up the vectors necessary for Rectangle creation.
 			stim::vec3<float> Y(1.0,0.0,0.0);		//orthogonal vec.
-			stim::vec3<float> pos(0.0,0.0,0.0);
-			stim::vec3<float> mag(1.0, 1.0, 1.0);
-			stim::vec3<float> dir(0.0, 0.0, 1.0);
+			stim::vec3<float> pos(0.0,0.0,0.0);		//center point of a rectangle
+			float mag = 1.0;				//size of the generated rectangle.
+			stim::vec3<float> dir(0.0, 0.0, 1.0);		//normal of the rectangle
  
-			float PHI[2], Z[2], range;
-			PHI[0] = solidAngle/2;
+			float PHI[2], Z[2], range;			
+			PHI[0] = solidAngle/2;				///Project the solid angle into spherical coordinates
 			PHI[1] = asin(0);
  
-			Z[0] = cos(PHI[0]);
+			Z[0] = cos(PHI[0]);				///Project the z into spherical coordinates
 			Z[1] = cos(PHI[1]);
  
-			range = Z[0] - Z[1];
+			range = Z[0] - Z[1];				///The range the possible values can be.
  
 			float z, theta, phi;
-			glNewList(dList, GL_COMPILE); 		
-			for(int i = 0; i < numSamples; i++)
+			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];
-				theta = uniformRandom()*2*M_PI;
-				phi = acos(z);
-				stim::vec3<float> sph(1, theta, phi);
-				stim::vec3<float> cart = sph.sph2cart();
-				dV.push_back(cart);
-				if(cos(Y.dot(cart)) < 0.087)
+				z = uniformRandom()*range + Z[1];	///generate a z coordinate
+				theta = 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.	
+				dV.push_back(cart);			///save the generated vector for further use.
+				#ifdef DEBUG
+//					std::cout << cart << std::endl;
+				#endif
+				if(cos(Y.dot(cart)) < 0.087)		///make sure that the Y is not parallel to the new vector.
 				{
 						Y[0] = 0.0; Y[1] = 1.0;
 				}else{
 					Y[0] = 1.0; Y[1] = 0.0;
 				}
-                                hor = stim::rect<float>(mag,
+
+                                hor = stim::rect<float>(mag,		///generate a rectangle with the new vectro as a normal.
 					 pos, cart,
                        			((Y.cross(cart)).cross(cart)).norm());
-       				ver = stim::rect<float>(mag,
+
+				#ifdef DEBUG
+			//		std::cout << hor.n() << std::endl;
+				#endif
+       				ver = stim::rect<float>(mag,		///generate another rectangle that's perpendicular the first but parallel to the cart vector.
 					 pos, cart,
                        			hor.n());
-				UpdateBuffer(0.0, 0.0+i*t_length);
+				UpdateBuffer(0.0, 0.0+i*n_pixels);	///Put the necessary points into the diplaylist.
 			}
-			glEndList();
+			glEndList();					///finilize the display list.
 		}
  
 		///@param float delta, How much the rectangles vary in position.
@@ -426,28 +405,45 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			//Set up the vectors necessary for Rectangle creation.
 			stim::vec3<float> Y(1.0,0.0,0.0);		//orthogonal vec.
-			stim::vec3<float> pos(0.0,0.0,0.0);
-			stim::vec3<float> mag(1.0, 1.0, 1.0);
-			stim::vec3<float> dir(0.0, 0.0, 1.0);
+			stim::vec3<float> pos(0.0,0.0,0.0);		//center point of a rectangle
+			float mag = 1.0;				///size of each rectangle
+			stim::vec3<float> dir(0.0, 0.0, 1.0);		///normal of the rectangle plane.
  
 			//Set up the variable necessary for vector creation.
-			glNewList(dList+1, GL_COMPILE);
-			for(int i = 0; i < numSamplesPos; i++)
+			glNewList(dList+1, GL_COMPILE);			///generate a new display list.
+			pV.push_back(pos);
+			hor = stim::rect<float>(mag,                    ///generate a rec     tangle with the new vector as a normal.  
+                	        pos, dir,
+	                       	((Y.cross(d)).cross(d))
+	                        .norm());
+                        ver = stim::rect<float>(mag,                    ///generate anoth     er rectangle that's perpendicular the first but parallel to the cart vector.
+        	                pos, dir,
+	                        hor.n()); 
+									///The first vector is always in the center.
+                        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();
-				temp = temp*delta*2.0 - delta/2.0;
+				stim::vec3<float> temp = uniformRandVector();	///generate a random point on a plane.
+				temp[0] = temp[0]*delta;
+				temp[1] = temp[1]*2*stim::PI;
+				
 				temp[2] = 0.0;
-				pV.push_back(temp);
-                		hor = stim::rect<float>(mag,
+				temp = temp.cyl2cart();
+				pV.push_back(temp);				///save the point for further use.
+                		hor = stim::rect<float>(mag,			///generate a rectangle with the new vector as a normal.
 					 temp, dir,
                       			((Y.cross(d)).cross(d))
 					.norm());
-                		ver = stim::rect<float>(mag,
+                		ver = stim::rect<float>(mag,			///generate another rectangle that's perpendicular the first but parallel to the cart vector.
 					 temp, dir,
                         		hor.n());
-				UpdateBuffer(0.0, 0.0+i*t_length);
+				UpdateBuffer(0.0, 0.0+i*n_pixels);		///sample the necessary points and put them into a display list.
 			}
-			glEndList();
+			glEndList();					///finilize the display list.
+			#ifdef DEBUG
+				for(int i = 0; i < numSamplesPos; i++)
+					std::cout << pV[i] << std::endl;
+			#endif
 		}
  
 		///@param float delta, How much the rectangles are allowed to expand.
@@ -462,35 +458,31 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
  
 			//Set up the vectors necessary for Rectangle creation.
 			stim::vec3<float> Y(1.0, 0.0, 0.0);		//orthogonal vec.
-			stim::vec3<float> pos(0.0, 0.0, 0.0);
-			stim::vec3<float> mag(1.0, 1.0, 1.0);
-			stim::vec3<float> dir(0.0, 0.0, 1.0);
+			stim::vec3<float> pos(0.0, 0.0, 0.0);		//center of the future rect.
+			float mag = 1.0;				///size of the rectangle
+			stim::vec3<float> dir(0.0, 0.0, 1.0);		///normal of the rectangle plane.
  
 			//Set up the variable necessary for vector creation.
-			int dim = (sqrt(numSamplesMag)-1)/2;
-			float min 	= 1.0-delta;
-			float max 	= 1.0+delta;
-			float step	= (max-min)/(numSamplesMag-1);
+			float min 	= 1.0-delta;			///smallest size
+			float max 	= 1.0+delta;			///largers size.
+			float step	= (max-min)/(numSamplesMag-1);	///the size variation from one rect to the next.
 			float factor;
-			stim::vec3<float> temp(0.0,0.0,0.0);
-
 			glNewList(dList+2, GL_COMPILE);
-			for(int i = 0; i < numSamplesMag; i++){
+			for(int i = 0; i < numSamplesMag; i++){		///for the number of position samples
 				//Create linear index
-				factor = (min+step*i)*mag[0];
-				temp = factor;
-				mV.push_back(temp);	
-				hor = stim::rect<float>(temp,
+				factor = (min+step*i)*mag;		///scaling factor
+				mV.push_back(factor);			///save the size factor for further use.
+				hor = stim::rect<float>(factor,		///generate a rectangle with the new vector as a normal.
 					 pos, dir, 
        	       				((Y.cross(d)).cross(d))
 					.norm());
-               			ver = stim::rect<float>(temp,
+               			ver = stim::rect<float>(factor,		///generate another rectangle that's perpendicular the first but parallel to the cart vector.
 					 pos, dir,
                        			hor.n());
-				UpdateBuffer(0.0, 0.0+i*t_length);
+				UpdateBuffer(0.0, 0.0+i*n_pixels);	///sample the necessary points and put them into a display list.
 			CHECK_OPENGL_ERROR
 			}
-			glEndList();
+			glEndList();					///finilize the displaylist.
 		}
  
 		///@param float v_x x-coordinate in buffer-space,
@@ -500,15 +492,15 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		UpdateBuffer(float v_x, float v_y)
 		{	
-			stim::vec3<float>p1; 
-        	        stim::vec3<float>p2; 
-	                stim::vec3<float>p3; 
-                	stim::vec3<float>p4;	
-			p1 = hor.p(1,1);
-			p2 = hor.p(1,0);
-			p3 = hor.p(0,0);
-			p4 = hor.p(0,1);
-			glBegin(GL_QUADS);
+			stim::vec3<float>p1; 				///first point.
+        	        stim::vec3<float>p2; 				///second point.
+	                stim::vec3<float>p3; 				///third point.
+                	stim::vec3<float>p4;				///fourth point.
+			p1 = hor.p(1,1);				///generate the top right point from the horizontal template.
+			p2 = hor.p(1,0);				///generate the bottom right point from the horizonatal template.
+			p3 = hor.p(0,0);				///generate the bottom left point from the horizontal template.
+			p4 = hor.p(0,1);				///generate the top left point from the horizonatal template.
+			glBegin(GL_QUADS);				///generate the Quad from the 4 points.
 				glTexCoord3f(
 					p1[0],
 					p1[1],
@@ -520,51 +512,51 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					p2[1],
 					p2[2]
 					);
-				glVertex2f(v_x+t_length, v_y);
+				glVertex2f(v_x+n_pixels, v_y);
 				glTexCoord3f(
 					p3[0],
 					p3[1],
 					p3[2]
 					);
-				glVertex2f(v_x+t_length, v_y+t_length);
+				glVertex2f(v_x+n_pixels, v_y+n_pixels);
 				glTexCoord3f(
 					p4[0],
 					p4[1],
 					p4[2]
 					);
-				glVertex2f(v_x, v_y+t_length);
-			 glEnd();
-
-			 p1 = ver.p(1,1);
-			 p2 = ver.p(1,0);
-			 p3 = ver.p(0,0);
-			 p4 = ver.p(0,1);
-		 	 glBegin(GL_QUADS);
+				glVertex2f(v_x, v_y+n_pixels);
+			 glEnd();					///finish the quad.
+
+			 p1 = ver.p(1,1);				///generate the top right point from the vertical template.
+			 p2 = ver.p(1,0);                               ///generate the bottom right point from the vertical template.
+			 p3 = ver.p(0,0);                               ///generate the bottom left point from the vertical template.
+			 p4 = ver.p(0,1);                               ///generate the top left point from the vertical template.
+		 	 glBegin(GL_QUADS);				///generate the Quad from the 4 points.
 				glTexCoord3f(
 					p1[0],
 					p1[1],
 					p1[2]
 					);
-				glVertex2f(v_x+t_length, v_y);
+				glVertex2f(v_x+n_pixels, v_y);
 				glTexCoord3f(
 					p2[0],
 					p2[1],
 					p2[2]
 					);
-				glVertex2f(v_x+2.0*t_length, v_y);
+				glVertex2f(v_x+2.0*n_pixels, v_y);
 				glTexCoord3f(
 					p3[0],
 					p3[1],
 					p3[2]
 					);
-				glVertex2f(v_x+2.0*t_length, v_y+t_length);
+				glVertex2f(v_x+2.0*n_pixels, v_y+n_pixels);
 				glTexCoord3f(
 					p4[0],
 					p4[1],
 					p4[2]
 					);
-				glVertex2f(v_x+t_length, v_y+t_length);
-			glEnd(); 
+				glVertex2f(v_x+n_pixels, v_y+n_pixels);
+			glEnd(); 					///finish the quad.
 		}
  
  
@@ -582,50 +574,24 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		GenerateFBO(unsigned int width, unsigned int height, GLuint &textureID, GLuint &framebufferID)
 		{
-			glDeleteFramebuffers(1, &framebufferID);
-			glGenFramebuffers(1, &framebufferID);
-			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
+			glDeleteFramebuffers(1, &framebufferID);				///clear the framebuffer.
+			glGenFramebuffers(1, &framebufferID);					///generate a clean buffer.
+			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);			///bind the new buffer.
 //			int numChannels = 1;
 //			unsigned char* texels = new unsigned char[width * height * numChannels];
-			glGenTextures(1, &textureID);
+			glGenTextures(1, &textureID);						///generate a texture that will attach to the buffer.
 			glBindTexture(GL_TEXTURE_2D, textureID);
  
 			//Textures repeat and use linear interpolation, luminance format.
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);		///Set up the texture to repeat at edges.
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);		
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);	///Set up the texture to use Linear interpolation
 			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
-			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,
+			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,				///Create the texture with no data.
 				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, NULL);   
 //			delete[] texels;
-			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
-			glBindTexture(GL_TEXTURE_2D, 0);
-		}
-
-		///@param uint width sets the width of the buffer.
-		///@param uint height sets the height of the buffer.
-		///Function for setting up the 2D buffer that stores the samples.
-		void
-		GenerateFBO(unsigned int width, unsigned int height)
-		{
-			glGenFramebuffers(1, &fboID);
-			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
-//			int numChannels = 1;
-//			unsigned char* texels = new unsigned char[width * height * numChannels];
-			glGenTextures(1, &texbufferID);
-			glBindTexture(GL_TEXTURE_2D, texbufferID);
-
-			//Textures repeat and use linear interpolation, luminance format.
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
-			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,
-				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, NULL);   
-//			delete[] texels;
-			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
-			glBindTexture(GL_TEXTURE_2D, 0);
-			CHECK_OPENGL_ERROR
+			glBindFramebuffer(GL_FRAMEBUFFER, 0); 					///Bind the frontbuffer
+			glBindTexture(GL_TEXTURE_2D, 0);					///Unbind the texture.
 		}
  
  
@@ -657,13 +623,12 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				//rotate to the current direction of the spider.
 				glRotatef(rot[0], rot[1], rot[2], rot[3]);
 				//scale to the magnitude of the spider.
-				glScalef(m[0],
-					 m[0],
-					 m[0]);
+				glScalef(m,
+					 m,
+					 m);
 				//get and store the current transformation matrix for later use.
 				glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
 				cT.set(curTrans);
-			//	printTransform();
  
 				CHECK_OPENGL_ERROR
 				//revert back to default gl mode.
@@ -683,31 +648,31 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		Bind()
 		{
-			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
+			glBindFramebuffer(GL_FRAMEBUFFER, direction_buffID);//set up GL buffer		
 			glFramebufferTexture2D(
 				GL_FRAMEBUFFER,
 				GL_COLOR_ATTACHMENT0,
 				GL_TEXTURE_2D,
-				texbufferID,
+				direction_texID,
 				0);
-			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
+			glBindFramebuffer(GL_FRAMEBUFFER, direction_buffID);
 			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
 			glDrawBuffers(1, DrawBuffers);
-			glBindTexture(GL_TEXTURE_2D, texbufferID);
+			glBindTexture(GL_TEXTURE_2D, direction_texID);
 			glClearColor(1,1,1,1);
 			glClear(GL_COLOR_BUFFER_BIT);
 			glMatrixMode(GL_PROJECTION);
 			glLoadIdentity();
 			glMatrixMode(GL_MODELVIEW);
 			glLoadIdentity();
-			glViewport(0,0,2.0*t_length, numSamples*t_length);
-			gluOrtho2D(0.0,2.0*t_length,0.0,numSamples*t_length);
+			glViewport(0,0,2.0*n_pixels, numSamples*n_pixels);
+			gluOrtho2D(0.0,2.0*n_pixels,0.0,numSamples*n_pixels);
 			glEnable(GL_TEXTURE_3D);
 			glBindTexture(GL_TEXTURE_3D, texID);
  
 			CHECK_OPENGL_ERROR
 		}
-		
+
 		///Method for controling the buffer and texture binding.
 		///Clears the buffer upon binding.
 		///@param GLuint &textureID, texture to be bound.
@@ -717,25 +682,25 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
 		{
  
-			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
-			glFramebufferTexture2D(
+			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);	///Bind the framebuffer.
+			glFramebufferTexture2D(					///associate it with the texture
 				GL_FRAMEBUFFER,
 				GL_COLOR_ATTACHMENT0,
 				GL_TEXTURE_2D,
 				textureID,
 				0);
-			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
-			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
-			glDrawBuffers(1, DrawBuffers);
-			glBindTexture(GL_TEXTURE_2D, textureID);
-			glMatrixMode(GL_PROJECTION);
+			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);	///Bind the framebuffer.
+			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};		///generate the drawbuffer.
+			glDrawBuffers(1, DrawBuffers);				///set the drawbuffer.
+			glBindTexture(GL_TEXTURE_2D, textureID);		///Bind the texture passed.
+			glMatrixMode(GL_PROJECTION);				///clear out the draw matrices
 			glLoadIdentity();
 			glMatrixMode(GL_MODELVIEW);
 			glLoadIdentity();
-			glViewport(0,0,2.0*len, nSamples*len);
-			gluOrtho2D(0.0,2.0*len,0.0,nSamples*len);
-			glEnable(GL_TEXTURE_3D);
-			glBindTexture(GL_TEXTURE_3D, texID);
+			glViewport(0,0,2.0*len, nSamples*len);			///set up viewport
+			gluOrtho2D(0.0,2.0*len,0.0,nSamples*len);		///Set up ortho
+			glEnable(GL_TEXTURE_3D);				
+			glBindTexture(GL_TEXTURE_3D, texID);			///bind the main texture (return to original state).
  
 			CHECK_OPENGL_ERROR
 		}
@@ -745,34 +710,16 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		Unbind()
 		{
 			//Finalize GL_buffer
-			glBindTexture(GL_TEXTURE_3D, 0);                      
+			glBindTexture(GL_TEXTURE_3D, 0);                      ///Bind the front buffer.
 			CHECK_OPENGL_ERROR
-			glBindTexture(GL_TEXTURE_2D, 0);                      
+			glBindTexture(GL_TEXTURE_2D, 0);                      ///Bind the default GL texture.
 			CHECK_OPENGL_ERROR
-			glBindFramebuffer(GL_FRAMEBUFFER, 0);
+			glBindFramebuffer(GL_FRAMEBUFFER, 0);		      ///Bind the defautl framebuffer.
 			CHECK_OPENGL_ERROR
-			glDisable(GL_TEXTURE_3D);
+			glDisable(GL_TEXTURE_3D);			      ///Turn off texturing.
 			CHECK_OPENGL_ERROR
 		}
  
-		///Makes the spider take a step.
-		///starting with the current p, d, m, find the next optimal p, d, m.
-		///Performs the branch detection on each step.
-		int
-		StepP()
-		{
-			Bind();
-			CHECK_OPENGL_ERROR
-			findOptimalDirection();
-			findOptimalPosition();
-			findOptimalScale();
-			Unbind();
-			Bind(btexbufferID, bfboID, 27);
-			branchDetection();
-			Unbind();
-			return current_cost;
-		}
-
  
  
  
@@ -785,42 +732,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///finds the minimum cost and sets the current_cost to that value.
 		/// and returns the index of the template with the minimal cost. 
 		int
-		getCost()
-		{
-			stim::vec<int> cost = 
-//				stim::cuda::get_cost(texbufferID, GL_TEXTURE_2D, numSamples);
-			cudaDeviceSynchronize();
-			current_cost = cost[1];
-			return cost[0];
-		}
-
-//		int
-//		getCost(GLuint tID, int n)
-//		{
-//			#ifdef TIMING
-//				gpuStartTimer();
-//			#endif
-//			stim::vec<int> cost = 
-//				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n, 2*t_length, t_length);
-//			#ifdef TIMING
-//				cost_time += gpuStopTimer();
-//			#endif
-//			current_cost = cost[1];
-//			return cost[0];
-//		}
-
-		int
 		getCost(cudaTextureObject_t tObj, float* result, int n)
 		{
 			#ifdef TIMING
 				gpuStartTimer();
 			#endif
 			stim::vec<int> cost = 
-				stim::cuda::get_cost(tObj, result, n, 2*t_length, t_length);
+				stim::cuda::get_cost(tObj, result, n, 2*n_pixels, n_pixels);	///call the cuda function with the appropriate texture buffer.
 			#ifdef TIMING
 				cost_time += gpuStopTimer();
 			#endif
-			current_cost = cost[1];
+			current_cost = cost[1];							///current cost.
 			return cost[0];
 		}
  
@@ -854,17 +776,20 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		gl_spider
 		(int samples = 1089, int samplespos = 441,int samplesmag = 144)
 		{
-//			std::cout << "I ran this constructor" << std::endl;
 			p = stim::vec3<float>(0.0, 0.0, 0.0);
 			d = stim::vec3<float>(0.0, 0.0, 1.0);
-			m = stim::vec<float>(1.0, 1.0);
+			m = 1.0;
 			S = stim::vec3<float>(1.0, 1.0, 1.0);
 			R = stim::vec3<float>(1.0, 1.0, 1.0);
-//			std::cout << samples << std::endl;
 			numSamples = samples;
-//			std::cout << numSamples << std::endl;
 			numSamplesPos = samplespos;
 			numSamplesMag = samplesmag;
+		#ifdef DEBUG
+			iter = 0;
+			iter_pos = 0;
+			iter_dir = 0;
+			iter_siz = 0;
+		#endif
 		}
  
 		///Position constructor: floats.
@@ -882,12 +807,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			p = stim::vec3<float>(pos_x, pos_y, pos_z);
 			d = stim::vec3<float>(dir_x, dir_y, dir_z);
-			m = stim::vec<float>(mag_x, mag_x, mag_x);
+			m = mag_x;
 			S = stim::vec3<float>(1.0,1.0,1.0);
 			R = stim::vec3<float>(1.0,1.0,1.0);
 			numSamples = numsamples;
 			numSamplesPos = numsamplespos;
 			numSamplesMag = numsamplesmag;
+		#ifdef DEBUG
+			iter = 0;
+			iter_pos = 0;
+			iter_dir = 0;
+			iter_siz = 0;
+		#endif
 		}
  
 		///Position constructor: vecs of floats.
@@ -900,12 +831,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			p = pos;
 			d = dir;
-			m = vec<float>(mag, mag, mag);
+			m = mag;
 			S = vec3<float>(1.0,1.0,1.0);
 			R = vec3<float>(1.0,1.0,1.0);
 			numSamples = samples;
 			numSamplesPos = samplesPos;
 			numSamplesMag = samplesMag;
+		#ifdef DEBUG
+			iter = 0;
+			iter_pos = 0;
+			iter_dir = 0;
+			iter_siz = 0;
+		#endif
 		}
  
 		///destructor	
@@ -913,14 +850,14 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		(void)
 		{
 			Unbind();
-			glDeleteTextures(1, &texbufferID);
-			glDeleteBuffers(1, &fboID);
-			glDeleteTextures(1, &ptexbufferID);
-			glDeleteBuffers(1, &pfboID);
-			glDeleteTextures(1, &mtexbufferID);
-			glDeleteBuffers(1, &mfboID);     
-			glDeleteTextures(1, &btexbufferID);
-			glDeleteBuffers(1, &bfboID);
+			glDeleteTextures(1, &direction_texID);
+			glDeleteBuffers(1, &direction_buffID);
+			glDeleteTextures(1, &position_texID);
+			glDeleteBuffers(1, &position_buffID);
+			glDeleteTextures(1, &radius_texID);
+			glDeleteBuffers(1, &radius_buffID);     
+			glDeleteTextures(1, &cylinder_texID);
+			glDeleteBuffers(1, &cylinder_buffID);
 		}
  
 		///@param GLuint id, texture that is going to be sampled.
@@ -939,45 +876,45 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			network_time = 0;
 			hit_time = 0;
 #endif
-			stepsize = 2.5;
-			t_length = 16.0;
+#ifdef DEBUG
+			iter = 0;
+			iter_pos = 0;
+			iter_dir = 0;
+			iter_siz = 0;
+#endif
+			stepsize = 10.0;
+			n_pixels = 16.0;
  
 			srand(100);	
 			texID = id;
-			//GenerateFBO(16, numSamples*8);
-			GenerateFBO(t_length*2, numSamples*t_length, texbufferID, fboID); 
-			std::cout << numSamples << std::endl;
+			GenerateFBO(n_pixels*2, numSamples*n_pixels, direction_texID, direction_buffID); 
 			CHECK_OPENGL_ERROR
-			GenerateFBO(t_length*2, numSamplesPos*t_length, ptexbufferID, pfboID); 
-			std::cout << numSamplesPos << std::endl;
+			GenerateFBO(n_pixels*2, numSamplesPos*n_pixels, position_texID, position_buffID); 
 			CHECK_OPENGL_ERROR
-			GenerateFBO(t_length*2, numSamplesMag*t_length, mtexbufferID, mfboID); 
-			std::cout << numSamplesMag << std::endl;
+			GenerateFBO(n_pixels*2, numSamplesMag*n_pixels, radius_texID, radius_buffID); 
 			CHECK_OPENGL_ERROR
-			GenerateFBO(16, 216, btexbufferID, bfboID); 
+			GenerateFBO(16, 216, cylinder_texID, cylinder_buffID); 
 			CHECK_OPENGL_ERROR
-			t_dir.MapCudaTexture(texbufferID, GL_TEXTURE_2D);
+			t_dir.MapCudaTexture(direction_texID, GL_TEXTURE_2D);
 			t_dir.Alloc(numSamples);
-			t_pos.MapCudaTexture(ptexbufferID, GL_TEXTURE_2D);
+			t_pos.MapCudaTexture(position_texID, GL_TEXTURE_2D);
 			t_pos.Alloc(numSamplesPos);
-			t_mag.MapCudaTexture(mtexbufferID, GL_TEXTURE_2D);
+			t_mag.MapCudaTexture(radius_texID, GL_TEXTURE_2D);
 			t_mag.Alloc(numSamplesMag);
-//			setDims(0.6, 0.6, 1.0);
-//			setSize(1024.0, 1024.0, 1024.0);
 			setMatrix();
 			dList = glGenLists(3);
 			glListBase(dList);
-			Bind(texbufferID, fboID, numSamples, t_length);
-				genDirectionVectors(5*M_PI/4);
+			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
+				genDirectionVectors(5*stim::PI/4);
 			Unbind();
-			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
-				genPositionVectors();
+			Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
+				genPositionVectors(0.2);
 			Unbind();
-			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
+			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
 				genMagnitudeVectors();
 			Unbind();
-			Bind(btexbufferID, bfboID, 27);
-			DrawCylinder();
+			Bind(cylinder_texID, cylinder_buffID, 27);
+//				DrawCylinder();
 			Unbind();
 		}
  
@@ -999,7 +936,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
  
 		///Returns the m vector.
-		stim::vec<float>
+		float	
 		getMagnitude()
 		{
 			return m;
@@ -1045,22 +982,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			d[2] = z;
 		}
  
-		///@param stim::vec<float> dir, the new d.
-		///Sets the m vector to the input vector mag.	
-		void
-		setMagnitude(stim::vec<float> mag)
-		{
-			m[0] = mag[0];
-			m[1] = mag[0];
-		}
  
 		///@param float mag, size of the sampled region.
 		///Sets the m vector to the input mag for both templates.
 		void
 		setMagnitude(float mag)
 		{
-			m[0] = mag;
-			m[1] = mag;
+			m = mag;
 		}
  
 		///@param float x, voxel size in the x direction.
@@ -1111,17 +1039,27 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		getRotation(stim::vec3<float> dir)
 		{
 			stim::vec<float> out(0.0,0.0,0.0,0.0);
-			stim::vec<float> from(0.0,0.0,1.0);
-			out[0] = acos(dir.dot(from))*180/M_PI;
-			if(out[0] < 1.0){
+			stim::vec3<float> from(0.0,0.0,1.0);
+			out[0] = acos(dir.dot(from))*180/stim::PI;
+			#ifdef DEBUG
+				std::cout << "out is " << out << std::endl;
+				std::cout << "when rotating from " << from << " to " << dir << std::endl;
+			#endif
+			if(out[0] < 0.01){
 				out[0] = 0.0;
 				out[1] = 0.0;
 				out[2] = 0.0;
 				out[3] = 1.0;
+			}
+			else if(out[0] < -180.0+0.01)
+			{
+				out[0] = 180.0;
+				out[1] = 1.0;
+				out[2] = 0.0;
+				out[3] = 0.0;
 			} else {
-				stim::vec<float> temp(0.0, 0.0, 0.0);;
-				stim::vec<float> dir1(dir[0], dir[1], dir[2]);
-				temp = (from.cross(dir1)).norm();
+				stim::vec3<float> temp(0.0, 0.0, 0.0);;
+				temp = (from.cross(dir)).norm();
 				out[1] = temp[0];
 				out[2] = temp[1];
 				out[3] = temp[2];
@@ -1152,7 +1090,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		setSeedMag(float mag)
 		{
-//			std::cout << "sMag: " << mag << std::endl;
 			seedsmags.push(mag);
 		}
  
@@ -1165,7 +1102,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		setSeed(float x, float y, float z)
 		{
-//			std::cout << "sPos: " << x << " " << y << " " << z << std::endl;
 			seeds.push(stim::vec<float>(x, y, z));
 		}
  
@@ -1176,7 +1112,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		setSeedVec(float x, float y, float z)
 		{
-//			std::cout << "sDir: " << x << " " << y << " " << z << std::endl;
 			seedsvecs.push(stim::vec<float>(x, y, z));
 		}
  
@@ -1297,7 +1232,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		saveNetwork(std::string name)
 		{
-			stim::glObj<float> sk;
+/*			stim::glObj<float> sk;
 			for(int i = 0; i < nt.sizeE(); i++)
 			{
 				std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
@@ -1310,7 +1245,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				}
 				sk.End();
 			}	
-			sk.save(name);
+*/			sk.save(name);
 		}
  
 		///Depreciated, but might be reused later()
@@ -1318,7 +1253,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		stim::glObj<float>
 		getNetwork()
 		{
-
+			return sk;
 		}
  
 		///returns a COPY of the entire stim::glnetwork object.
@@ -1333,7 +1268,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		GLuint
 		getFB()
 		{
-			return bfboID;
+			return cylinder_buffID;
 		}
  
 //--------------------------------------------------------------------------//
@@ -1357,14 +1292,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		int
 		Step()
 		{
-			Bind(texbufferID, fboID, numSamples, t_length);
+			#ifdef DEBUG
+			std::cerr << "Took a step" << std::endl;
+			#endif
+			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
 			CHECK_OPENGL_ERROR
 				findOptimalDirection();
 			Unbind();
-			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
+			Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
 				findOptimalPosition();
 			Unbind();
-			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
+			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
 				findOptimalScale();
 			Unbind();
 			CHECK_OPENGL_ERROR
@@ -1385,16 +1323,16 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 //--------------------------------------------------------------------------//
  
 		void
-		MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
+		MonteCarloDirectionVectors(int nSamples, float solidAngle = stim::TAU)
 		{
-			float PHI[2], Z[2], range;
-			PHI[0] = asin(solidAngle/2);
-			PHI[1] = asin(0);
+//			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]);
+//			Z[0] = cos(PHI[0]);
+//			Z[1] = cos(PHI[1]);
  
-			range = Z[0] - Z[1];
+//			range = Z[0] - Z[1];
  
 			std::vector<stim::vec3<float> > vecsUni;
 			for(int i = 0; i < numSamplesPos; i++)
@@ -1428,8 +1366,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			 int j = 0;
 			 	for(float i = step; i <= 360.0; i += step)
 			 	{
-					 x=r0*cos(i*2.0*M_PI/360.0);
-					 y=r0*sin(i*2.0*M_PI/360.0);
+					 x=r0*cos(i*stim::TAU/360.0);
+					 y=r0*sin(i*stim::TAU/360.0);
 					 glTexCoord3f(x,y,z0); 
 					 glVertex2f(0.0, j*6.4+6.4);
 //					 glVertex2f(0.0, j*27.0+27.0);
@@ -1451,12 +1389,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
  
 ///need to return the cylinder.
+///SOMETHING MIGHT BE GOING ON HERE IN GENERATE BUFFER.
 		void
 		DrawLongCylinder(int n = 8, int l_template = 8,int l_square = 8)
 		{
 			int cylLen = cL.size()-1;
-			GenerateFBO(n*l_square, cylLen*l_template, btexbufferID, bfboID);
-			Bind(btexbufferID, bfboID, cylLen, l_template*l_square/2.0);
+			GenerateFBO(n*l_square, cylLen*l_template, cylinder_texID, cylinder_buffID);
+			Bind(cylinder_texID, cylinder_buffID, cylLen, l_template*l_square/2.0);
 			stim::cylinder<float> cyl(cL, cM);
 			std::vector<std::vector<stim::vec3<float> > > p = cyl.getPoints(n);
 			for(int i = 0; i < p.size()-1; i++)	///number of circles
@@ -1502,14 +1441,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				seeds.pop();
 				seedsvecs.pop();
 				seedsmags.pop();
-//				std::cout << "The current seed Vector is " << curSeedVec << std::endl;
 				setPosition(curSeed);
 				setDirection(curSeedVec);
 				setMagnitude(curSeedMag);
+
+				#ifdef DEBUG
+					std::cout << "The new seed " << curSeed << curSeedVec << curSeedMag << std::endl;
+				#endif
+				
+//			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
+//			CHECK_OPENGL_ERROR
+//				findOptimalDirection();
+//			Unbind();
+//THIS IS EXPERIMENTAL
+			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
+				findOptimalScale();
+			Unbind();
+//THIS IS EXPERIMENTAL
+
 //				cL.push_back(curSeed);
 //				cM.push_back(curSeedMag);
 //				cD.push_back(curSeedMag);
-				pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost);
+				traceLine(p, m, min_cost);
 			}
 		}
  
@@ -1556,17 +1509,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					 ds[0], ds[1], ds[2],
 					 ups[0], ups[1], ups[2]);
  
-//				sk.Render();
-				nt.Render();								
+				sk.Render();
+//				nt.Render();								
  
 				CHECK_OPENGL_ERROR
  
  
-//				glLoadName((int) sk.numL());
-				glLoadName(nt.sizeE());
+				glLoadName((int) sk.numL());
+//				glLoadName(nt.sizeE());
  
-//				sk.RenderLine(cL);
-				nt.RenderLine(cL);	
+				sk.RenderLine(cL);
+//				nt.RenderLine(cL);	
  
 //				glPopName();
 				glFlush();
@@ -1615,8 +1568,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			cL.clear();
 			cM.clear();
-		}	
+		}
  
+/*
 		void
 		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos, 
 				stim::vec<float> smag, stim::vec3<float> sdir)
@@ -1664,13 +1618,26 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				network_time += nt * 1000.0;
 			#endif
 		}
+*/
+		void
+		addToNetwork(std::vector<stim::vec3<float> > L, std::vector<float > M)
+		{
+			if(L.size() > 3)
+			{
+				sk.Begin(stim::OBJ_LINE);
+				for(int i = 0; i < L.size(); i++)
+				{
+					sk.TexCoord(M[i]);
+					sk.Vertex(L[i][0], L[i][1], L[i][2]);
+				}
+				sk.End();
  
-//		void
-//		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos, 
-//				stim::vec<float> smag, stim::vec3<float> sdir)
-//		{
-//
-//		}
+				nt.addEdge(L,M);
+				#ifdef DEBUG
+					iter++;
+				#endif
+			}
+		}
  
  
 		void
@@ -1678,23 +1645,22 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			std::cout << nt.sizeE() << " edges " << std::endl;
 			std::cout << nt.sizeV() << " nodes " << std::endl;
-			
 		}
  
-		std::pair<stim::fiber<float>, int >
-		traceLine(stim::vec3<float> pos, stim::vec<float> mag, int min_cost)
+		void
+		traceLine(stim::vec3<float> pos, float mag, int min_cost)
 		{
 			//starting (seed) position and magnitude.
 			stim::vec3<float> spos = getPosition();
-			stim::vec<float> smag = getMagnitude();
+			float smag = getMagnitude();
 			stim::vec3<float> sdir = getDirection();	
  
-			Bind();
+//			Bind();
 //			sk.Begin(stim::OBJ_LINE);
  
  
-//			sk.createFromSelf(GL_SELECT);
-			nt.createFromSelf(GL_SELECT);
+			sk.createFromSelf(GL_SELECT);
+//			nt.createFromSelf(GL_SELECT);
  
 			cL.push_back(pos);
 			cM.push_back(mag);
@@ -1710,12 +1676,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				int cost = Step();
 				if (cost > min_cost){
 					running = false;
-//					sk.End();
 					branchDetection2();
-					pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
-					addToNetwork(a, spos, smag, sdir);
-//					std::cout << "the cost of " << cost << " > " << min_cost << std::endl;
-					return a;
+					addToNetwork(cL, cM);
+					#ifdef DEBUG
+					std::cerr << "the cost of " << cost << " > " << min_cost << std::endl;
+					#endif
 					break;
 				} else {
 					//Have we found the edge of the map?
@@ -1726,10 +1691,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					{
 						running = false;
 						branchDetection2();
-						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
-						addToNetwork(a, spos, smag, sdir);
-//						std::cout << "I hit and edge" << std::endl;
-						return a;
+						addToNetwork(cL, cM);
+						#ifdef DEBUG
+						std::cerr << "I hit and edge" << std::endl;
+						#endif
 						break;
 					}
 					//If this is the first step in the trace,
@@ -1741,39 +1706,40 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					}
 					//Has the template size gotten unreasonable?
 					mag = getMagnitude();
-					if(mag[0] > 75 || mag[0] < 1){
+					if(mag > 75 || mag < 1){
 						running = false;
 						branchDetection2();
-						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
-						addToNetwork(a, spos, smag, sdir);
-//						std::cout << "The templates are too big" << std::endl;
-						return a;
+						addToNetwork(cL, cM);
+						#ifdef DEBUG
+						std::cerr << "The templates are too big" << std::endl;
+						#endif
 						break;
 					}
 					else
 					{
-						h = selectObject(p, getDirection(), m[0]);
+						h = selectObject(p, getDirection(), m);
 						//Have we hit something previously traced?
 						if(h != -1){
-//							std::cout << "I hit the fiber " << h << std::endl;
+							#ifdef DEBUG
+							std::cerr << "I hit the fiber " << h << std::endl;
+							#endif
 							running = false;
 							branchDetection2();
-							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
-							addToNetwork(a, spos, smag, sdir);
-							return a;
+							addToNetwork(cL, cM);
 							break;
 						}
 						else {  
 							cL.push_back(stim::vec3<float>(p[0], p[1],p[2]));
-							cM.push_back(stim::vec<float>(m[0], m[0]));
-//							Bind(btexbufferID, bfboID, 27);
-							Unbind();
+							cM.push_back(m);
+//							Unbind();
 							CHECK_OPENGL_ERROR
-							
 						}
 				 	}
                  		}
          		}
+				#ifdef DEBUG
+				std::cout << "I broke out!" << std::endl;
+				#endif
 		}
 };
 }
@@ -64,7 +64,7 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 			case 4:
 				return GL_RGBA;
 			default:
-				std::cout<<"Error in stim::gl_texture - unable to guess texture format based on number of channels ("<<R[4]<<")"<<std::endl;
+				std::cout<<"Error in stim::gl_texture - unable to guess texture format based on number of channels" << std::endl;
 				exit(1);
 			}
 		}
@@ -236,7 +236,7 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 			if(texID == 0) generate_texture();					//generate the texture if it doesn't already exist
 			else{
 				std::cout<<"Texture has already been attached to a context."<<std::endl;
-				exit(1);
+				//exit(1);     Why do we need to quit if it's attached. print the message and continue.
 			}
 		}
  
@@ -48,9 +48,9 @@ public:
 	CUDA_CALLABLE
 	circle(T size, T z_pos = (T)0) : plane<T>()
 	{
-		init();
 		center(stim::vec3<T>(0,0,z_pos));
 		scale(size);
+		init();
 	}
  
 	///create a rectangle from a center point, normal
@@ -88,6 +88,7 @@ public:
 	{
 		init();
 		setU(u);
+//		U = u;
 		center(c);
 		normal(n);
 		scale(s);
@@ -174,6 +175,13 @@ public:
 		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
@@ -80,6 +80,15 @@ public:
 		return sph;
 	}
  
+	CUDA_CALLABLE vec3<T> cart2cyl() const{
+		vec3<T> cyl;
+		cyl.ptr[0] = sqrt(pow(ptr[0],2) + pow(ptr[1],2));
+		cyl.ptr[1] = atan(ptr[1]/ptr[0]);
+		cyl.ptr[2] = ptr[2];
+
+		return cyl;
+	}
+
 	/// Convert the vector from cartesian to spherical coordinates (r, theta, phi -> x, y, z where theta = [0, 2*pi])
 	CUDA_CALLABLE vec3<T> sph2cart() const{
 		vec3<T> cart;
@@ -90,6 +99,16 @@ public:
 		return cart;
 	}
  
+	/// Convert the vector from cylindrical to cart coordinates (r, theta, z -> x, y, z where theta = [0, 2*pi])
+	CUDA_CALLABLE vec3<T> cyl2cart() const{
+		vec3<T> cart;
+		cart.ptr[0] = ptr[0] * std::cos(ptr[1]);
+		cart.ptr[1] = ptr[0] * std::sin(ptr[1]);
+		cart.ptr[2] = ptr[2];
+
+		return cart;
+	}
+
 	/// Computes the normalized vector (where each coordinate is divided by the L2 norm)
 	CUDA_CALLABLE vec3<T> norm() const{
         vec3<T> result;
@@ -2,19 +2,27 @@
 #define STIM_CYLINDER_H
 #include <iostream>
 #include <stim/math/circle.h>
-#include <stim/math/vec3.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;
+		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;
 		std::vector< T > L;			//length of the cylinder at each position.
+
+
+		using stim::centerline<T>::c;
+		using stim::centerline<T>::N;
  
 		///default init	
 		void
@@ -23,6 +31,20 @@ class cylinder
  
 		}
  
+		///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 radii (inM)
 		void
 		init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
@@ -31,6 +53,10 @@ class cylinder
 			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;
  
@@ -46,7 +72,9 @@ class cylinder
  
 			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]);
@@ -54,16 +82,24 @@ class cylinder
 				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.
@@ -72,18 +108,42 @@ class cylinder
 		{
 			if(idx == 0)
 			{
-				return (e[idx+1].P - e[idx].P).norm();
+				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 == e.size()-1)
+			else if(idx == N-1)
 			{
-				return (e[idx].P - e[idx-1].P).norm();
+				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<float> v2 = (e[idx+1].P-e[idx].P).norm();
-				return (v1+v2).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;	
  
@@ -92,14 +152,17 @@ class cylinder
 		stim::vec3<T>
 		d(T l, int idx)
 		{
-			if(idx == 0 || idx == e.size()-1)
+			if(idx == 0 || idx == N-1)
 			{
-				return e[idx].N;
+				return norms[idx];
+//				return e[idx].N;
 			}
 			else
 			{
+				
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-				return(	e[idx].N + (e[idx+1].N - e[idx].N)*rat);
+				return(	norms[idx] + (norms[idx+1] - norms[idx])*rat);
+//				return(	e[idx].N + (e[idx+1].N - e[idx].N)*rat);
 			} 
 		}
  
@@ -137,6 +200,7 @@ class cylinder
 	public:
 		///default constructor
 		cylinder()
+		// : centerline<T>()
 		{
  
 		}
@@ -144,15 +208,33 @@ class cylinder
 		///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){
+		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){
-			std::vector< stim::vec<T> > inM;						//create an array of arbitrary magnitudes
+		cylinder(std::vector< stim::vec3<T> > inP)
+			: centerline<T>(inP)
+		{
+			std::vector< T > inM;						//create an array of arbitrary magnitudes
  
 			stim::vec<T> zero;
 			zero.push_back(0);
@@ -164,7 +246,7 @@ class cylinder
 		///Returns the number of points on the cylinder centerline
  
 		unsigned int size(){
-			return e.size();
+			return N;
 		}
  
  
@@ -180,9 +262,23 @@ class cylinder
 			}
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
-				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
-		}
+			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.
@@ -192,7 +288,19 @@ class cylinder
 		p(T l, int idx)
 		{
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
+				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])));
 		}
@@ -209,8 +317,11 @@ class cylinder
 			}
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
-			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*((l-L[idx])/(L[idx+1]- L[idx])));
-		}
+			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.
@@ -220,7 +331,23 @@ class cylinder
 		r(T l, int idx)
 		{
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
+			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
@@ -233,7 +360,7 @@ class cylinder
 		/// 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 < e.size(); p++)
+			for(unsigned int p = 0; p < N; p++)
 				mags[p].push_back(m);
 		}
  
@@ -279,8 +406,8 @@ class cylinder
 		getPoints(int sides)
 		{
 			std::vector<std::vector <vec3<T> > > points;
-			points.resize(e.size());
-			for(int i = 0; i < e.size(); i++)
+			points.resize(N);
+			for(int i = 0; i < N; i++)
 			{
 				points[i] = e[i].getPoints(sides);
 			}
@@ -319,7 +446,7 @@ class cylinder
 			T len = L[0];						//allocate space for the segment length
  
 			//for every consecutive point in the cylinder
-			for(unsigned p = 1; p < e.size(); p++){
+			for(unsigned p = 1; p < N; p++){
  
 				//p1 = pos[p];							//get the position and magnitude for the next point
 				m1 = mags[p][m];
@@ -5,6 +5,7 @@
 #include <GL/glut.h>
 #include <stim/math/vector.h>
 #include <stim/visualization/obj.h>
+#include <stim/gl/error.h>
  
  
 namespace stim
@@ -35,9 +36,9 @@ private:
 		glListBase(dList);
  
 		glMatrixMode(GL_PROJECTION);
-		glLoadIdentity;
+		glLoadIdentity();
 		glMatrixMode(GL_MODELVIEW);
-		glLoadIdentity;
+		glLoadIdentity();
  
 	}
  
@@ -129,7 +130,7 @@ public:
 	}
  
 	void
-	RenderLine(std::vector< stim::vec<T> > l)
+	RenderLine(std::vector< stim::vec3<T> > l)
 	{
 		glColor3f(0.5, 1.0, 0.5);
 		glLineWidth(3.0);
@@ -322,13 +322,6 @@ protected:
 		return OBJ_INVALID;
 	}
  
-	//private method returns the vertex indices for a line
-	std::vector<unsigned> get_l_v(unsigned l){
-
-
-
-	}
-
 public:
 	/// Constructor loads a Wavefront OBJ file
 	obj(std::string filename){