STABLE: made template size a variable, changed the gl_spider and cuda code to support that.

Pavel Govyadinov
1 parent f2c883d8
Showing 5 changed files with 245 additions and 106 deletions Show diff stats
stim/cuda/cuda_texture.cuh
stim/cuda/spider_cost.cuh
stim/cuda/testKernel.cuh
stim/gl/gl_spider.h
stim/math/vector.h
@@ -134,6 +134,9 @@ namespace stim
 			HANDLE_ERROR(
 				cudaDestroyTextureObject(tObj)
 			);
+//			HANDLE_ERROR(
+//				cudaFreeArray(srcArray)
+//			);
 		}
 //-------------------------------------------------------------------------//
@@ -17,7 +17,7 @@ namespace stim{
 	stim::cuda::cuda_texture t;  //texture object.
 	float* result;
-	float* print;
+//	float* print;
 	///Initialization function, allocates the memory and passes the necessary
 	///handles from OpenGL and Cuda.
@@ -42,7 +42,9 @@ namespace stim{
 	__device__
 	float Template(int x)
 	{
-		if(x < 16/6 || x > 16*5/6 || (x > 16*2/6 && x < 16*4/6)){
+		if(x < 32/6 || x > 32*5/6 || (x > 32*2/6 && x < 32*4/6)){
+//		if(x < 2 || x > 13 || (x > 5 && x < 10)){
+//	if(x < ceilf(16/6) || x > floorf(16*5/6) || (x > floorf(16*2/6) && x < ceilf(16*4/6))){
 			return 1.0;
 		}else{
 			return 0.0;
@@ -57,14 +59,15 @@ namespace stim{
 	///@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 *result)
+	void get_diff (cudaTextureObject_t texIn, float *result, int dx, int dy)
 	{       
-		__shared__ float shared[16][8];
+//		__shared__ float shared[32][16];
+		extern __shared__ float shared[];
 		int x   = threadIdx.x + blockIdx.x * blockDim.x;
 		int y   = threadIdx.y + blockIdx.y * blockDim.y;
 		int x_t = threadIdx.x;
 		int y_t = threadIdx.y;
-//		int idx = y*16+x;
+		int idx = y_t*dx+x_t;
 		int g_idx = blockIdx.y;
 		float valIn             = tex2D<unsigned char>(texIn, x, y)/255.0;
@@ -72,7 +75,7 @@ namespace stim{
 //		print[idx]              = abs(valIn);             ///temporary
-		shared[x_t][y_t]        = abs(valIn-valTemp);
+		shared[idx]        = abs(valIn-valTemp);
 		__syncthreads();
@@ -81,7 +84,8 @@ namespace stim{
 			__syncthreads();
 			if (x_t < step)
 			{
-				shared[x_t][y_t] += shared[x_t + step][y_t];
+//				shared[x_t][y_t] += shared[x_t + step][y_t];
+				shared[idx] += shared[y_t*dx+x_t+step];
 			}
 		__syncthreads();
 		}
@@ -92,13 +96,14 @@ namespace stim{
 			__syncthreads();
 			if(y_t < step)
 			{
-				shared[x_t][y_t] += shared[x_t][y_t + step];
+//				shared[x_t][y_t] += shared[x_t][y_t + step];
+				shared[idx] += shared[(y_t+step)*dx+x_t];
 			}
 		__syncthreads();
 		}
 		__syncthreads();
 		if(x_t == 0 && y_t == 0)
-			result[g_idx] = shared[0][0];
+			result[g_idx] = shared[0];
 	//      //result[idx]           = abs(valIn);
@@ -112,7 +117,7 @@ namespace stim{
 	///				  may work with other gl texture types, but untested.
 	///@param DIM_Y, the number of samples in the template.
 	extern "C"
-	stim::vec<int> get_cost(GLint texbufferID, GLenum texType, int DIM_Y)
+	stim::vec<int> get_cost(GLint texbufferID, GLenum texType, int DIM_Y,int dx = 16, int dy = 8)
 	{
 		//Bind the Texture in GL and allow access to cuda.
@@ -132,10 +137,10 @@ namespace stim{
 		//cuda launch variables.
 		dim3 numBlocks(1, DIM_Y);
-		dim3 threadsPerBlock(16, 8);
+		dim3 threadsPerBlock(dx, dy);
-		get_diff <<< numBlocks, threadsPerBlock >>> (t.getTexture(), result);
+		get_diff <<< numBlocks, threadsPerBlock, dx*dy*sizeof(float) >>> (t.getTexture(), result, dx, dy);
 		HANDLE_ERROR(
 			cudaMemcpy(output, result, DIM_Y*sizeof(float), cudaMemcpyDeviceToHost)
@@ -155,6 +160,7 @@ namespace stim{
 		t.UnmapCudaTexture();
 		cleanUP();
 		ret[0] = idx; ret[1] = (int) output[idx];
+//		std::cout << "The cost is " << output[idx] << std::endl;
 		free(output);
 		return ret;
 	}
@@ -29,7 +29,7 @@
         __device__
         float templ(int x)                                                                      
         {
-                if(x < 16/6 || x > 16*5/6 || (x > 16*2/6 && x < 16*4/6)){
+                if(x < 32/6 || x > 32*5/6 || (x > 32*2/6 && x < 32*4/6)){
                         return 1.0;
                 }else{
                         return 0.0;
@@ -48,9 +48,9 @@
 	{       
 		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*16+x;
+//		int idx = y*64+x;
+		int idx = y*32+x;
+	//	int idx = y*16+x;
 		float valIn             = tex2D<unsigned char>(texIn, x, y);
 		float templa		= templ(x);
@@ -96,3 +96,33 @@
 		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);
+
+		cudaDeviceSynchronize();
+		stringstream name;      //for debugging
+		name << nam.c_str();
+		stim::gpu2image<float>(print, name.str(),x,y,0,255);
+	  
+		tx.UnmapCudaTexture();
+		cleanUP();
+	}
+
@@ -24,7 +24,7 @@
 #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 <stim/cuda/testKernel.cuh>
@@ -103,6 +103,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		stim::vec<float> ps;
 		stim::vec<float> ups;
 		stim::vec<float> ds;
+
+		static const float t_length = 16.0;
 //--------------------------------------------------------------------------//
@@ -117,9 +119,12 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			setMatrix();			//create the transformation matrix.
 			glCallList(dList);		//move the templates to p, d, m.
+			#ifdef TESTING
+			test(texbufferID, GL_TEXTURE_2D,2*t_length,numSamples*t_length, "Final_Cost_Direction.bmp");
+			#endif
 			int best = getCost(texbufferID,numSamples);		//find min cost.
 			stim::vec<float> next(		//find next vector.
- 			dV[best][0]*S[0]*R[0],
+			dV[best][0]*S[0]*R[0],
 			dV[best][1]*S[1]*R[1],
 			dV[best][2]*S[2]*R[2],
 			0);
@@ -139,6 +144,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			setMatrix();			//create the transformation matrix.
 			glCallList(dList+1);		//move the templates to p, d, m.
+			#ifdef TESTING
+			test(ptexbufferID, GL_TEXTURE_2D,2*t_length, numSamplesPos*t_length, "Final_Cost_Position.bmp");
+			#endif
 			int best = getCost(ptexbufferID, numSamplesPos);		//find min cost.
 //			std::cerr << best << std::endl;
 			stim::vec<float> next(		//find next position.
@@ -162,6 +170,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			setMatrix();			//create the transformation.
 			glCallList(dList+2);		//move the templates to p, d, m.
+			#ifdef TESTING
+			test(mtexbufferID, GL_TEXTURE_2D, 2*t_length, numSamplesMag*t_length, "Final_Cost_Position.bmp");
+			#endif
 			int best = getCost(mtexbufferID, numSamplesMag);		//get best cost.
 			setMagnitude(m[0]*mV[best][0]);	//adjust the magnitude.
 		}
@@ -263,11 +274,90 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			}
 		}
+
+		float uniformRandom()
+		{
+			return (  (float)(rand()))/(  (float)(RAND_MAX));
+		}
+
+		float normalRandom()
+		{
+			float u1 = uniformRandom();
+			float u2 = uniformRandom();
+			return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));
+		}
+
+		stim::vec<float> uniformRandVector()
+		{
+			stim::vec<float> r(uniformRandom(), uniformRandom(), 1.0);
+			return r;
+		}
+
+		stim::vec<float> normalRandVector()
+		{
+			stim::vec<float> r(normalRandom(), normalRandom(), 1.0);
+			return r;
+		}
+
 //--------------------------------------------------------------------------//
 //---------------------TEMPLATE CREATION METHODS----------------------------//
 //--------------------------------------------------------------------------//
+		void
+		MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
+		{
+///equal area projection.
+//			std::vector<stim::vec<float> > vecsUni;
+//			for(int i = 0; i < nSamples; i++)
+//			{
+//				stim::vec<float> a = uniformRandVector();
+//				a[1] = a[1]*2*M_PI;
+//				a[0] = a[0]*(solidAngle);
+//				a = a.cyl2cart();
+//				stim::vec<float> b(
+//					sqrt(1.-(pow(a[0],2)+(pow(a[1],2)))/4)*a[0],
+//					sqrt(1.-(pow(a[0],2)+(pow(a[1],2)))/4)*a[1],
+//					1.-(pow(a[0],2)+pow(a[1],2))/2);
+//					2.*a[0]/(1.+pow(a[0],2)+pow(a[1],2)),
+//					2.*a[1]/(1.+pow(a[0],2)+pow(a[1],2)),
+//					(-1.0+pow(a[0],2)+pow(a[1],2))/(1.0+pow(a[0],2)+pow(a[1],2)));
+//				b = b.norm();
+//				b[2] = -b[2]; //this doeswn't want to change
+//				vecsUni.push_back(b);
+//			}
+			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];
+
+
+			float z, theta, phi;
+			
+			std::vector<stim::vec<float> > vecsUni;
+			for(int i = 0; i < numSamplesPos; i++)
+			{
+				stim::vec<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();
+		}
+
+
 		///@param solidAngle, the size of the arc to sample.
 		///Method for populating the vector arrays with sampled vectors.
 		///Objects created are rectangles the with the created directions.
@@ -277,46 +367,45 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		genDirectionVectors(float solidAngle = 5/M_PI*4)
 		{
+
 			//Set up the vectors necessary for Rectangle creation.
 			vec<float> Y(1.0,0.0,0.0);		//orthogonal vec.
 			vec<float> pos(0.0,0.0,0.0);
 			vec<float> mag(1.0, 1.0, 1.0);
 			vec<float> dir(0.0, 0.0, 1.0);
-			//Set up the variable necessary for vector creation.
-			vec<float> d_s = d.cart2sph().norm();
-			vec<float> temp(0,0,0);
-			int dim = (sqrt(numSamples)-1)/2;
-			float p0  	= -M_PI;				 //phi angle in spherical coordinates.
-			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0)); //step size in Theta.
-			float dp  	= p0/(2.0*((float)dim + 1.0));		 //step size in Phi.
+			float PHI[2], Z[2], range;
+			PHI[0] = solidAngle/2;
+			PHI[1] = asin(0);
+			Z[0] = cos(PHI[0]);
+			Z[1] = cos(PHI[1]);
+			
+			range = Z[0] - Z[1];
+
+			float z, theta, phi;
 			glNewList(dList, GL_COMPILE); 		
-			//Loop over the above defined space creating distinct vectors. 
-			int idx = 0;
-			for(int i = -dim; i <= dim; i++){
-				for(int j = -dim; j <= dim; j++){
-					
-					//Create linear index
-					idx = (j+dim)+(i+dim)*((dim*2)+1);	
-					temp[0] = d_s[0]; 			//rotate vector
-					temp[1] = d_s[1]+dp*(float) i;
-					temp[2] = d_s[2]+dt*(float) j;
-					
-					temp = (temp.sph2cart()).norm();	//back to cart
-					dV.push_back(temp);
-				 	if(cos(Y.dot(temp))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
-					else{Y[0] = 1.0; Y[1] = 0.0;}
-
-                                       	hor = stim::rect<float>(mag,
-						 pos, temp,
-                       				((Y.cross(temp)).cross(temp)).norm());
-       				 	ver = stim::rect<float>(mag,
-						 pos, temp,
-                       				hor.n());
-					UpdateBuffer(0.0, 0.0+idx*8.0);
-					CHECK_OPENGL_ERROR
+			for(int i = 0; i < numSamples; i++)
+			{
+				z = uniformRandom()*range + Z[1];
+				theta = uniformRandom()*2*M_PI;
+				phi = acos(z);
+				stim::vec<float> sph(1, theta, phi);
+				stim::vec<float> cart = sph.sph2cart();
+				dV.push_back(cart);
+				if(cos(Y.dot(cart)) < 0.087)
+				{
+						Y[0] = 0.0; Y[1] = 1.0;
+				}else{
+					Y[0] = 1.0; Y[1] = 0.0;
 				}
+                                hor = stim::rect<float>(mag,
+					 pos, cart,
+                       			((Y.cross(cart)).cross(cart)).norm());
+       				ver = stim::rect<float>(mag,
+					 pos, cart,
+                       			hor.n());
+				UpdateBuffer(0.0, 0.0+i*t_length);
 			}
 			glEndList();
 		}
@@ -337,36 +426,21 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			vec<float> dir(0.0, 0.0, 1.0);
 			//Set up the variable necessary for vector creation.
-			vec<float> temp(0,0,0);
-			int dim = (sqrt(numSamplesPos)-1)/2;	//number of position vectors.
-			stim::rect<float> samplingPlane =	//plane from which we pull position samples
-				 stim::rect<float>(p, d);
-			samplingPlane.scale(mag[0]*delta, mag[0]*delta);
-			float step = 1.0/(dim);			//step size.
-
-			//Loop over the samples, keeping the original p samples in the center of the resulting texture to create a large number of position vectors.
-			int idx;
 			glNewList(dList+1, GL_COMPILE);
-			for(int i = -dim; i <= dim; i++){
-				for(int j = -dim; j <= dim; j++){
-					//Create linear index
-					idx = (j+dim)+(i+dim)*((dim*2)+1);	
-
-					temp = samplingPlane.p(
-							0.5+step*i,
-								 0.5+step*j
-										);
-					pV.push_back(temp);
-                			hor = stim::rect<float>(mag,
-						 temp, dir,
-                      				((Y.cross(d)).cross(d))
-						.norm());
-                			ver = stim::rect<float>(mag,
-						 temp, dir,
-                        			hor.n());
-					UpdateBuffer(0.0, 0.0+idx*8.0);
-				CHECK_OPENGL_ERROR
-				}
+			for(int i = 0; i < numSamplesPos; i++)
+			{
+				stim::vec<float> temp = uniformRandVector();
+				temp = temp*delta*2.0 - delta/2.0;
+				temp[2] = 0.0;
+				pV.push_back(temp);
+                		hor = stim::rect<float>(mag,
+					 temp, dir,
+                      			((Y.cross(d)).cross(d))
+					.norm());
+                		ver = stim::rect<float>(mag,
+					 temp, dir,
+                        		hor.n());
+				UpdateBuffer(0.0, 0.0+i*t_length);
 			}
 			glEndList();
 		}
@@ -382,8 +456,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			//Set up the vectors necessary for Rectangle creation.
-			vec<float> Y(1.0,0.0,0.0);		//orthogonal vec.
-			vec<float> pos(0.0,0.0,0.0);
+			vec<float> Y(1.0, 0.0, 0.0);		//orthogonal vec.
+			vec<float> pos(0.0, 0.0, 0.0);
 			vec<float> mag(1.0, 1.0, 1.0);
 			vec<float> dir(0.0, 0.0, 1.0);
@@ -408,7 +482,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
                			ver = stim::rect<float>(temp,
 					 pos, dir,
                        			hor.n());
-				UpdateBuffer(0.0, 0.0+i*8.0);
+				UpdateBuffer(0.0, 0.0+i*t_length);
 			CHECK_OPENGL_ERROR
 			}
 			glEndList();
@@ -421,7 +495,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		UpdateBuffer(float v_x, float v_y)
 		{	
-			float len = 8.0;
 			stim::vec<float>p1; 
         	        stim::vec<float>p2; 
 	                stim::vec<float>p3; 
@@ -442,19 +515,19 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					p2[1],
 					p2[2]
 					);
-				glVertex2f(v_x+len, v_y);
+				glVertex2f(v_x+t_length, v_y);
 				glTexCoord3f(
 					p3[0],
 					p3[1],
 					p3[2]
 					);
-				glVertex2f(v_x+len, v_y+len);
+				glVertex2f(v_x+t_length, v_y+t_length);
 				glTexCoord3f(
 					p4[0],
 					p4[1],
 					p4[2]
 					);
-				glVertex2f(v_x, v_y+len);
+				glVertex2f(v_x, v_y+t_length);
 			 glEnd();
 			 p1 = ver.p(1,1);
@@ -467,25 +540,25 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					p1[1],
 					p1[2]
 					);
-				glVertex2f(v_x+len, v_y);
+				glVertex2f(v_x+t_length, v_y);
 				glTexCoord3f(
 					p2[0],
 					p2[1],
 					p2[2]
 					);
-				glVertex2f(v_x+2.0*len, v_y);
+				glVertex2f(v_x+2.0*t_length, v_y);
 				glTexCoord3f(
 					p3[0],
 					p3[1],
 					p3[2]
 					);
-				glVertex2f(v_x+2.0*len, v_y+len);
+				glVertex2f(v_x+2.0*t_length, v_y+t_length);
 				glTexCoord3f(
 					p4[0],
 					p4[1],
 					p4[2]
 					);
-				glVertex2f(v_x+len, v_y+len);
+				glVertex2f(v_x+t_length, v_y+t_length);
 			glEnd(); 
 		}
@@ -522,7 +595,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			delete[] texels;
 			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
 			glBindTexture(GL_TEXTURE_2D, 0);
-			CHECK_OPENGL_ERROR
 		}
 		///@param uint width sets the width of the buffer.
@@ -604,7 +676,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///Method for controling the buffer and texture binding.
 		///Clears the buffer upon binding.
 		void
-		Bind(float len = 8.0)
+		Bind()
 		{
 			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
 			glFramebufferTexture2D(
@@ -623,8 +695,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			glLoadIdentity();
 			glMatrixMode(GL_MODELVIEW);
 			glLoadIdentity();
-			glViewport(0,0,2.0*len, numSamples*len);
-			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
+			glViewport(0,0,2.0*t_length, numSamples*t_length);
+			gluOrtho2D(0.0,2.0*t_length,0.0,numSamples*t_length);
 			glEnable(GL_TEXTURE_3D);
 			glBindTexture(GL_TEXTURE_3D, texID);
@@ -639,6 +711,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
 		{
+			
+//			std::cout << glGetIntegerv(GL_MAX_TEXTURE_SIZE) << std::endl;
+			GLint max;
+			glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
+//			std::cout << max << std::endl;
 			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
 			glFramebufferTexture2D(
 				GL_FRAMEBUFFER,
@@ -742,7 +819,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				start = std::clock();
 			#endif
 			stim::vec<int> cost = 
-				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n);
+				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n, 2*t_length, t_length);
 			cudaDeviceSynchronize();
 			#ifdef TESTING
 				duration_cuda = duration_cuda +
@@ -760,6 +837,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		initCuda()
 		{	
 			stim::cudaSetDevice();
+			MonteCarloDirectionVectors(500);
 			//GLint max;
 			//glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
 			//std::cout << max << std::endl;
@@ -790,12 +868,15 @@ 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 = vec<float>(0.0, 0.0, 0.0);
 			d = vec<float>(0.0, 0.0, 1.0);
 			m = vec<float>(1.0, 1.0);
 			S = vec<float>(1.0, 1.0, 1.0);
 			R = vec<float>(1.0, 1.0, 1.0);
+//			std::cout << samples << std::endl;
 			numSamples = samples;
+//			std::cout << numSamples << std::endl;
 			numSamplesPos = samplespos;
 			numSamplesMag = samplesmag;
 		}
@@ -863,24 +944,29 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		attachSpider(GLuint id)
 		{
+			srand(100);	
 			texID = id;
 			//GenerateFBO(16, numSamples*8);
-			GenerateFBO(16, numSamples*8, texbufferID, fboID); 
-			GenerateFBO(16, numSamplesPos*8, ptexbufferID, pfboID); 
-			GenerateFBO(16, numSamplesMag*8, mtexbufferID, mfboID); 
+			GenerateFBO(t_length*2, numSamples*t_length, texbufferID, fboID); 
+			CHECK_OPENGL_ERROR
+			GenerateFBO(t_length*2, numSamplesPos*t_length, ptexbufferID, pfboID); 
+			CHECK_OPENGL_ERROR
+			GenerateFBO(t_length*2, numSamplesMag*t_length, mtexbufferID, mfboID); 
+			CHECK_OPENGL_ERROR
 			GenerateFBO(16, 216, btexbufferID, bfboID); 
+			CHECK_OPENGL_ERROR
 //			setDims(0.6, 0.6, 1.0);
 //			setSize(1024.0, 1024.0, 1024.0);
 			setMatrix();
 			dList = glGenLists(3);
 			glListBase(dList);
-			Bind(texbufferID, fboID, numSamples);
+			Bind(texbufferID, fboID, numSamples, t_length);
 				genDirectionVectors(5*M_PI/4);
 			Unbind();
-			Bind(ptexbufferID, pfboID, numSamplesPos);
+			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
 				genPositionVectors();
 			Unbind();
-			Bind(mtexbufferID, mfboID, numSamplesMag);
+			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
 				genMagnitudeVectors();
 			Unbind();
 			Bind(btexbufferID, bfboID, 27);
@@ -1058,7 +1144,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		setSeedMag(float mag)
 		{
-			std::cout << "sMag: " << mag << std::endl;
+//			std::cout << "sMag: " << mag << std::endl;
 			seedsmags.push(mag);
 		}
@@ -1071,7 +1157,7 @@ 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;
+//			std::cout << "sPos: " << x << " " << y << " " << z << std::endl;
 			seeds.push(stim::vec<float>(x, y, z));
 		}
@@ -1082,7 +1168,7 @@ 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;
+//			std::cout << "sDir: " << x << " " << y << " " << z << std::endl;
 			seedsvecs.push(stim::vec<float>(x, y, z));
 		}
@@ -1222,17 +1308,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		int
 		Step()
 		{
-			Bind(texbufferID, fboID, numSamples);
+			Bind(texbufferID, fboID, numSamples, t_length);
 			CHECK_OPENGL_ERROR
 			#ifdef TESTING
 				start = std::clock();
 			#endif
 				findOptimalDirection();
 			Unbind();
-			Bind(ptexbufferID, pfboID, numSamplesPos);
+			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
 				findOptimalPosition();
 			Unbind();
-			Bind(mtexbufferID, mfboID, numSamplesMag);
+			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
 				findOptimalScale();
 			Unbind();
 			CHECK_OPENGL_ERROR
@@ -1548,6 +1634,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					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;
 					break;
 				} else {
@@ -1561,6 +1648,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						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;
 						break;
 					}
@@ -1578,6 +1666,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						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;
 						break;
 					}
@@ -1586,6 +1675,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						h = selectObject(p, getDirection(), m[0]);
 						//Have we hit something previously traced?
 						if(h != -1){
+//							std::cout << "I hit the fiber " << h << std::endl;
 							running = false;
 							branchDetection2();
 							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
@@ -133,6 +133,16 @@ struct vec : public std::vector&lt;T&gt;
 	}
+	
+	vec<T> cyl2cart() const
+	{
+		vec<T> cyl;
+		cyl.push_back(at(0)*std::sin(at(1)));
+		cyl.push_back(at(0)*std::cos(at(1)));
+		cyl.push_back(at(2));
+		return(cyl);
+		
+	}
 	/// Convert the vector from cartesian to spherical coordinates (x, y, z -> r, theta, phi where theta = [0, 2*pi])
 	vec<T> cart2sph() const
 	{