Merge branch 'Glnetwork'

David Mayerich
2 parents 5bccf89d bf23ee36
Showing 25 changed files with 2925 additions and 517 deletions Show diff stats
stim/cuda/arraymath/array_cart2polar.cuh
stim/cuda/branch_detection.cuh
stim/cuda/branch_detection2.cuh
stim/cuda/cuda_texture.cuh
stim/cuda/cudatools/devices.h
stim/cuda/ivote/update_dir.cuh
stim/cuda/ivote/vote.cuh
stim/cuda/sharedmem.cuh
stim/cuda/spider_cost.cuh
stim/cuda/templates/conv2.cuh
stim/cuda/templates/conv2sep.cuh
stim/cuda/templates/gaussian_blur.cuh
stim/cuda/testKernel.cuh
stim/gl/gl_spider.h
stim/gl/gl_texture.h
stim/math/plane.h
stim/math/plane_old.h
stim/math/quaternion.h
stim/math/rect.h
stim/math/rect_old.h
@@ -4,7 +4,7 @@
 namespace stim{
 	namespace cuda{
 		template<typename T>
-		__global__ void cuda_cart2polar(T* a, int x, int y){
+		__global__ void cuda_cart2polar(T* a, int x, int y, float rotation){
  
  
 			// calculate the 2D coordinates for this current thread.
@@ -20,21 +20,21 @@ namespace stim{
 				float yl = a[i * 2 + 1];
 				float theta = atan2( yl, xl ) ;
 				float r = sqrt(xl * xl + yl * yl);
-				a[i * 2 + 0] = theta;
+				a[i * 2 + 0] = theta + rotation;
 				a[i * 2 + 1] = r;
  
 		}
  
  
 		template<typename T>
-		void gpu_cart2polar(T* gpuGrad, unsigned int x, unsigned int y){
+		void gpu_cart2polar(T* gpuGrad, unsigned int x, unsigned int y, float rotation = 0){
  
 			unsigned int max_threads = stim::maxThreadsPerBlock();
 			dim3 threads(max_threads, 1);
 			dim3 blocks(x/threads.x + (x %threads.x == 0 ? 0:1) , y);
  
 			//call the kernel to do the multiplication
-			cuda_cart2polar <<< blocks, threads >>>(gpuGrad, x, y);
+			cuda_cart2polar <<< blocks, threads >>>(gpuGrad, x, y, rotation);
  
 		}
  
@@ -67,4 +67,4 @@ namespace stim{
 	}
 }
  
-#endif
 \ No newline at end of file
+#endif
+#include <iostream>
+#include <fstream>
+#include <cuda_runtime.h>
+#include <stim/math/vector.h>
+//#include <math.h>
+#include <stim/visualization/colormap.h>
+#include <stim/cuda/cuda_texture.cuh>
+#include <stim/cuda/templates/gradient.cuh>
+#include <stim/cuda/templates/gaussian_blur.cuh>
+#include <stim/cuda/arraymath.cuh>
+#include <stim/cuda/ivote.cuh>
+#include <stim/cuda/testKernel.cuh>
+typedef unsigned int uint;
+typedef unsigned int uchar;
+
+stim::cuda::cuda_texture t;	
+float*		gpuTable;
+float*		gpuGrad;
+float*		gpuVote;	
+float*		gpuI;
+float*		gpuCenters;
+
+void atan_2d(float* cpuTable, unsigned int rmax)
+{
+	//initialize the width and height of the window which atan2 are computed in.
+	int xsize = 2*rmax +1;
+	int ysize = 2*rmax +1;
+	
+	// assign the center coordinates of the atan2 window to yi and xi
+	int yi = rmax;
+	int xi = rmax;
+	
+
+	for (int xt = 0; xt < xsize; xt++){
+
+		for(int yt = 0; yt < ysize; yt++){
+
+			//convert the current 2D coordinates to 1D
+			int id = yt * xsize + xt;
+			// calculate the distance between the pixel and the center of the atan2 window
+			float xd = xi - xt;
+			float yd = yi - yt;
+
+			// calculate the angle between the pixel and the center of the atan2 window and store the result.
+			float atan_2d_vote = atan2(yd, xd);
+			cpuTable[id] = atan_2d_vote;
+		}
+	}
+
+}
+
+void initCuda(unsigned int bytes_table, unsigned int bytes_ds)
+{
+	HANDLE_ERROR(
+		cudaMalloc((void**) &gpuTable, bytes_table)
+		);
+	HANDLE_ERROR(
+		cudaMalloc((void**) &gpuI, bytes_ds)
+		);
+	HANDLE_ERROR(
+		cudaMalloc((void**) &gpuGrad,  bytes_ds*2)
+		);
+	HANDLE_ERROR(
+		cudaMalloc((void**) &gpuVote,  bytes_ds)
+		);
+	HANDLE_ERROR(
+		cudaMalloc((void**) &gpuCenters, bytes_ds)
+		);
+}
+
+void cleanCuda()
+{
+	HANDLE_ERROR(
+		cudaFree(gpuTable)
+	);
+	HANDLE_ERROR(
+		cudaFree(gpuGrad)
+	);
+	HANDLE_ERROR(
+		cudaFree(gpuVote)
+	);
+	HANDLE_ERROR(
+		cudaFree(gpuCenters)
+	);
+	HANDLE_ERROR(
+		cudaFree(gpuI)
+	);
+}
+
+std::vector< stim::vec<float> >
+find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
+{
+	float 		phi	 	= 15.1*M_PI/180;
+	int		iter		= 5;
+	float 		dphi		= phi/iter;
+	float 		rmax 		= 10;
+	float		sigma		= 4;
+	unsigned int 	pixels 		= x * y;
+	unsigned int 	bytes  		= sizeof(float) * pixels;
+	unsigned int 	bytes_table	= sizeof(float) * (2*rmax + 1) * (2*rmax + 1);
+	unsigned int 	x_ds		= (x + (x % 1 == 0 ? 0:1));
+	unsigned int 	y_ds		= (y + (x % 1 == 0 ? 0:1));
+	unsigned int	bytes_ds	= sizeof(float) * x_ds * y_ds;
+	unsigned int	conn		= 5;
+	float		final_t		= 200.0;
+	float*		cpuTable	= (float*) malloc(bytes_table);
+	float*		cpuCenters	= (float*) malloc(bytes_ds);
+
+	stringstream name;
+
+
+
+
+	std::vector<stim::vec<float> >  output;
+	initCuda(bytes_table, bytes_ds); 
+
+	atan_2d(cpuTable, rmax);
+	cudaMemcpy(gpuTable, cpuTable, bytes_table, cudaMemcpyHostToDevice);
+
+
+	t.MapCudaTexture(texbufferID, texType);
+	cudaDeviceSynchronize();
+	stim::cuda::tex_gaussian_blur2<float>(
+		gpuI, sigma, x, y, t.getTexture(), t.getArray()
+		);
+	cudaDeviceSynchronize();
+
+
+	stim::cuda::gpu_gradient_2d<float>(
+		gpuGrad, gpuI, x, y
+		);
+	cudaDeviceSynchronize();
+	
+	stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y);
+	cudaDeviceSynchronize();
+
+	cudaDeviceSynchronize();
+	for (int i = 0; i < iter; i++)
+	{
+		stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
+	cudaDeviceSynchronize();
+		stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
+	cudaDeviceSynchronize();
+		phi = phi - dphi;
+	}
+	
+	cudaDeviceSynchronize();
+	stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, final_t, conn, x, y);
+	cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
+	for(int i = 0; i < pixels; i++)
+	{
+		int ix = (i % x);
+		int iy = (i / x);
+		if((cpuCenters[i] == 1) && (ix > 4) && (ix < x-4))
+		{
+
+			float x_v = (float) ix;
+			float y_v = (float) iy;
+			output.push_back(stim::vec<float>((x_v/(float)x),
+							  (y_v/(float)y), 0.0));	
+
+		}
+	}
+
+
+	t.UnmapCudaTexture();
+	cleanCuda();
+	free(cpuTable);
+	free(cpuCenters);
+	return output;
+}
+#include <stim/cuda/templates/gaussian_blur.cuh>
+#include <stim/cuda/templates/gradient.cuh>
+#include <stim/cuda/arraymath.cuh>
+#include <stim/cuda/ivote.cuh>
+
+
+
+
+
+
+
+
+
+
+void atan_2(float* cpuTable, unsigned int rmax){
+
+	//initialize the width and height of the window which atan2 are computed in.
+	int xsize = 2*rmax +1;
+	int ysize = 2*rmax +1;
+	
+	// assign the center coordinates of the atan2 window to yi and xi
+	int yi = rmax;
+	int xi = rmax;
+	
+
+	for (int xt = 0; xt < xsize; xt++){
+
+		for(int yt = 0; yt < ysize; yt++){
+
+			//convert the current 2D coordinates to 1D
+			int id = yt * xsize + xt;
+			// calculate the distance between the pixel and the center of the atan2 window
+			float xd = xi - xt;
+			float yd = yi - yt;
+
+			// calculate the angle between the pixel and the center of the atan2 window and store the result.
+			float atan_2d_vote = atan2(yd, xd);
+			cpuTable[id] = atan_2d_vote;
+		}
+	}
+
+}
+std::vector<stim::vec<float> > 
+find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
+{
+
+	float* cpuTable		= (float
+
+	unsigned int pixels = x * y;
+	unsigned int bytes = sizeof(float) * pixels;
+
+	//calculate the number of bytes in the atan2 table
+
+	unsigned int bytes_table = (2*rmax+1) * (2*rmax+1) * sizeof(float);
+
+
+
+	//allocate space on the GPU for the atan2 table
+
+	float* gpuTable;
+
+	cudaMalloc(&gpuTable, bytes_table);
+
+
+
+	cudaMemcpy(gpuTable, cpuTable, bytes_table, cudaMemcpyHostToDevice);
+
+	unsigned int sigma_ds = 1/resize;
+	unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));
+	unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
+	unsigned int bytes_ds = sizeof(float) * x_ds * y_ds;
+	
+
+	float* gpuI;
+	cudaMalloc(&gpuI, bytes_ds);
+
+	
+	float* gpuGrad;
+	cudaMalloc(&gpuGrad, bytes_ds*2);
+
+	float* gpuVote;
+	cudaMalloc(&gpuVote, bytes_ds);
+
+	// allocate space on the GPU for the detected cell centes
+
+	float* gpuCenters;
+
+	cudaMalloc(&gpuCenters, bytes_ds);		
+
+
+	stim::cuda::gpu_down_sample<float>(gpuI, gpuI0, resize, x , y);
+	cudaMemcpy(cpuResize, gpuI, bytes_ds, cudaMemcpyDeviceToHost);
+
+x = x_ds;
+	y = y_ds;
+	t = t * resize;
+	//sigma = sigma * resize;
+
+	cudaDeviceSynchronize();
+	stim::cuda::gpu_gaussian_blur2<float>(gpuI,sigma, x, y);
+	cudaDeviceSynchronize();
+	cudaMemcpy(cpuBlur, gpuI, bytes_ds, cudaMemcpyDeviceToHost);
+	cudaDeviceSynchronize();
+	
+	stim::cuda::gpu_gradient_2d<float>(gpuGrad, gpuI, x, y);
+	cudaDeviceSynchronize();
+	cudaMemcpy(cpuGradient, gpuGrad, bytes_ds*2, cudaMemcpyDeviceToHost);
+
+	stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y);
+	cudaDeviceSynchronize();
+	cudaMemcpy(cpuCart2Polar, gpuGrad, bytes_ds*2, cudaMemcpyDeviceToHost);
+	
+
+	//multiply the gradient by a constant and calculate the absolute value (to save an image)	
+
+	stim::cuda::cpu_multiply<float>(cpuCart2Polar, 40, x * y * 2);
+
+	cudaDeviceSynchronize();
+
+	stim::cuda::cpu_abs<float>(cpuCart2Polar, x * y * 2);
+
+	cudaDeviceSynchronize();
+
+		
+	for (int i =0; i<iter; i++){
+		
+		stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
+		cudaDeviceSynchronize();
+		stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
+		cudaDeviceSynchronize();
+		switch (i){
+		case 0 : cudaMemcpy(cpuVote1, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		case 1 : cudaMemcpy(cpuVote2, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		case 2 : cudaMemcpy(cpuVote3, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		case 3 : cudaMemcpy(cpuVote4, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		case 4 : cudaMemcpy(cpuVote5, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		default : cudaMemcpy(cpuVote5, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);
+			break;
+		}
+		phi = phi - dphi;
+	}
+	
+	stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, t, conn, x, y);
+	cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
+	
+}
+#ifndef STIM_CUDA_TEXTURE_H
+#define STIM_CUDA_TEXTURE_H
+
+#include <assert.h>
+#include <stim/cuda/cudatools/error.h>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <stdio.h>
+#include <GL/glew.h>
+#include <GL/glut.h>
+#include <sstream>
+#include <stim/visualization/colormap.h>
+#include <stim/cuda/cudatools/devices.h>
+#include <stim/cuda/cudatools/threads.h>
+#include <stim/math/vector.h>
+
+///A container for the texture based methods used by the spider class.
+namespace stim
+{
+	namespace cuda
+	{
+		class cuda_texture
+		{
+		public:
+			cudaArray* srcArray;
+			cudaGraphicsResource_t resource;
+			struct cudaResourceDesc resDesc;
+			struct cudaTextureDesc texDesc;
+			cudaTextureObject_t tObj;
+
+
+			///basic constructor that creates the texture with default parameters.
+			cuda_texture()
+			{
+			memset(&texDesc, 0, sizeof(texDesc));
+				texDesc.addressMode[0] 		= cudaAddressModeWrap;
+				texDesc.addressMode[1] 		= cudaAddressModeWrap;
+				texDesc.filterMode		= cudaFilterModePoint;
+				texDesc.readMode		= cudaReadModeElementType;
+				texDesc.normalizedCoords	= 0;
+			}
+
+//-------------------------------------------------------------------------//
+//-------------------------------CUDA_MAPPING------------------------------//
+//-------------------------------------------------------------------------//
+//Methods for creating the cuda texture.
+		///@param GLuint tex	 	 -- GLtexture (must be contained in a frame buffer object)
+		///				    that holds that data that will be handed to cuda.
+		///@param GLenum target		 -- either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D
+		///				    map work with other gl texture types but untested.
+		///Maps the gl texture in cuda memory, binds that data to a cuda array, and binds the cuda
+		///array to a cuda texture.
+		void
+		MapCudaTexture(GLuint tex, GLenum target)
+		{
+			HANDLE_ERROR(
+				cudaGraphicsGLRegisterImage(
+					&resource,
+					tex,
+					target,
+//					cudaGraphicsMapFlagsReadOnly
+					cudaGraphicsRegisterFlagsNone	
+				)
+			);
+
+			HANDLE_ERROR(
+				cudaGraphicsMapResources(1, &resource)
+			);
+			
+			HANDLE_ERROR(
+				cudaGraphicsSubResourceGetMappedArray(&srcArray, resource, 0, 0)
+			);
+			
+			memset(&resDesc, 0, sizeof(resDesc));
+			resDesc.resType = cudaResourceTypeArray;
+			resDesc.res.array.array = srcArray;
+			HANDLE_ERROR(
+				cudaCreateTextureObject(&tObj, &resDesc, &texDesc, NULL)
+			);
+		}
+
+		///Unmaps the gl texture,  binds that data to a cuda array, and binds the cuda
+		///array to a cuda texture.
+		void
+		UnmapCudaTexture()
+		{
+			HANDLE_ERROR(
+				cudaGraphicsUnmapResources(1, &resource)
+			);
+			HANDLE_ERROR(
+				cudaGraphicsUnregisterResource(resource)
+			);	
+			HANDLE_ERROR(
+				cudaDestroyTextureObject(tObj)
+			);
+		}
+
+//-------------------------------------------------------------------------//
+//------------------------------GET/SET METHODS----------------------------//
+//-------------------------------------------------------------------------//
+
+///Returns the bound texture object.
+		cudaTextureObject_t
+		getTexture()
+		{
+			return tObj;
+		}
+
+		cudaArray*
+		getArray()
+		{
+			return srcArray;
+		}
+	};
+}
+}
+
+
+#endif
@@ -4,7 +4,7 @@
 #include <cuda.h>
  
 namespace stim{
-
+extern "C"
 int maxThreadsPerBlock()
 {
 	int device;
@@ -14,6 +14,7 @@ int maxThreadsPerBlock()
 	return props.maxThreadsPerBlock;
 }
  
+extern "C"
 int sharedMemPerBlock()
 {
 	int device;
@@ -164,6 +164,9 @@ namespace stim{
 			//free allocated memory
 			cudaFree(gpuDir);
  
+			cudaDestroyTextureObject(texObj);
+			cudaFreeArray(cuArray);
+
 		}
  
 		template<typename T>
@@ -211,4 +214,4 @@ namespace stim{
 	}
 }
  
-#endif
 \ No newline at end of file
+#endif
@@ -124,6 +124,9 @@ namespace stim{
  
 			cuda_vote <<< blocks, threads,share_bytes >>>(gpuVote, texObj, gpuTable, phi, rmax, x , y);
  
+			cudaDestroyTextureObject(texObj);
+			cudaFreeArray(cuArray);
+
 		}
  
  
@@ -169,4 +172,4 @@ namespace stim{
 	}
 }
  
-#endif
 \ No newline at end of file
+#endif
@@ -34,9 +34,38 @@ namespace stim{
 				}
 			}
 		}
+
+		template<typename T, typename D>
+		__device__ void sharedMemcpy_tex2D(T* dest, cudaTextureObject_t src,
+										 unsigned int x, unsigned int y, unsigned int X, unsigned int Y,
+										 dim3 threadIdx, dim3 blockDim){
+
+			//calculate the number of iterations required for the copy
+			unsigned int xI, yI;
+			xI = X/blockDim.x + 1;				//number of iterations along X
+			yI = Y/blockDim.y + 1;				//number of iterations along Y
+
+			//for each iteration
+			for(unsigned int xi = 0; xi < xI; xi++){
+				for(unsigned int yi = 0; yi < yI; yi++){
+
+					//calculate the index into shared memory
+					unsigned int sx = xi * blockDim.x + threadIdx.x;
+					unsigned int sy = yi * blockDim.y + threadIdx.y;
+
+					//calculate the index into the texture
+					unsigned int tx = x + sx;
+					unsigned int ty = y + sy;
+
+					//perform the copy
+					if(sx < X && sy < Y)
+						dest[sy * X + sx] = abs(255 - tex2D<D>(src, tx, ty));
+				}
+			}
+		}
  
 	}
 }
  
  
-#endif
 \ No newline at end of file
+#endif
+#ifndef STIM_SPIDER_COST_H
+#define STIM_SPIDER_COST_H
+
+#include <assert.h>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <stdio.h>
+#include <stim/visualization/colormap.h>
+#include <sstream>
+#include <stim/math/vector.h>
+#include <stim/cuda/cudatools/devices.h>
+#include <stim/cuda/cudatools/threads.h>
+#include <stim/cuda/cuda_texture.cuh>
+namespace stim{
+	namespace cuda
+	{
+	
+	stim::cuda::cuda_texture t;  //texture object.
+	float* result;
+	float* print;
+	
+	///Initialization function, allocates the memory and passes the necessary
+	///handles from OpenGL and Cuda.
+	///@param DIM_Y			--integer controlling how much memory to allocate.
+	void initArray(int DIM_Y)
+	{
+//			cudaMalloc( (void**) &print, DIM_Y*16*sizeof(float));     ///temporary
+			cudaMalloc( (void**) &result, DIM_Y*sizeof(float));
+	}
+
+	///Deinit function that frees the memery used and releases the texture resource
+	///back to OpenGL.
+	void cleanUP()
+	{
+			cudaFree(result);
+//			cudaFree(print);         ///temporary
+	}  
+
+	///A virtual representation of a uniform template.
+	///Returns the value of the template pixel.
+	///@param int x			--location of a pixel.
+	__device__
+	float Template(int x)
+	{
+		if(x < 16/6 || x > 16*5/6 || (x > 16*2/6 && x < 16*4/6)){
+			return 1.0;
+		}else{
+			return 0.0;
+		}
+
+	}
+		
+	///Find the difference of the given set of samples and the template
+	///using cuda acceleration.
+	///@param stim::cuda::cuda_texture t	--stim texture that holds all the references
+	///					  to the data.
+	///@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)
+	{       
+		__shared__ float shared[16][8];
+		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 g_idx = blockIdx.y;
+
+		float valIn             = tex2D<unsigned char>(texIn, x, y)/255.0;
+		float valTemp           = Template(x);
+
+//		print[idx]              = abs(valIn);             ///temporary
+
+		shared[x_t][y_t]        = abs(valIn-valTemp);
+
+		__syncthreads();
+
+		for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
+		{
+			__syncthreads();
+			if (x_t < step)
+			{
+				shared[x_t][y_t] += shared[x_t + step][y_t];
+			}
+		__syncthreads();
+		}
+		__syncthreads();
+
+		for(unsigned int step = blockDim.y/2; step >= 1; step >>= 1)
+		{
+			__syncthreads();
+			if(y_t < step)
+			{
+				shared[x_t][y_t] += shared[x_t][y_t + step];
+			}
+		__syncthreads();
+		}
+		__syncthreads();
+		if(x_t == 0 && y_t == 0)
+			result[g_idx] = shared[0][0];
+
+
+	//      //result[idx]           = abs(valIn);
+	}
+
+
+	///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.
+	extern "C"
+	stim::vec<int> get_cost(GLint texbufferID, GLenum texType, int DIM_Y)
+	{
+
+		//Bind the Texture in GL and allow access to cuda.
+		t.MapCudaTexture(texbufferID, texType);
+
+		//initialize the return arrays.
+		float* output;	
+		output = (float* ) malloc(DIM_Y*sizeof(float));
+
+		stim::vec<int> ret(0, 0);
+		initArray(DIM_Y);
+		
+
+		//variables for finding the min.
+		float mini = 10000000000000000.0;
+		int idx = 0;
+	
+		//cuda launch variables.
+		dim3 numBlocks(1, DIM_Y);
+		dim3 threadsPerBlock(16, 8);
+
+
+		get_diff <<< numBlocks, threadsPerBlock >>> (t.getTexture(), result);
+
+		HANDLE_ERROR(
+			cudaMemcpy(output, result, DIM_Y*sizeof(float), cudaMemcpyDeviceToHost)
+			);
+
+		for( int i = 0; i<DIM_Y; i++){
+			if(output[i] < mini){
+				mini = output[i];
+				idx = i;                                                                   
+			}
+		}
+
+//		stringstream name;      //for debugging
+//		name << "Test.bmp";
+//		stim::gpu2image<float>(print, name.str(),16,218,0,256);
+	  
+		t.UnmapCudaTexture();
+		cleanUP();
+		ret[0] = idx; ret[1] = (int) output[idx];
+		free(output);
+		return ret;
+	}
+
+	}
+}
+
+
+#endif
@@ -102,8 +102,10 @@ namespace stim{
 			dim3 blocks(w / threads + 1, h);
  
 			//call the kernel to do the multiplication
-			cuda_conv2 <<< blocks, threads >>>(mask, copy, texObj, w, h, M);
-
+			//cuda_conv2 <<< blocks, threads >>>(img, mask, copy, w, h, M);
+			cuda_conv2 <<< blocks, threads >>>(img, mask, copy, texObj, w, h, M);
+			cudaDestroyTextureObject(texObj);
+			cudaFreeArray(cuArray);
 		}
  
 		template<typename T>
@@ -139,4 +141,4 @@ namespace stim{
 }
  
  
-#endif
 \ No newline at end of file
+#endif
@@ -30,7 +30,7 @@ namespace stim{
 			int byi = blockIdx.y;
  
 			//copy the portion of the image necessary for this block to shared memory
-			stim::cuda::sharedMemcpy_tex2D(s, in, bxi - kr, byi, 2 * kr + blockDim.x, 1, threadIdx, blockDim);
+			stim::cuda::sharedMemcpy_tex2D<float, unsigned char>(s, in, bxi - kr, byi, 2 * kr + blockDim.x, 1, threadIdx, blockDim);
  
 			//calculate the thread index
 			int ti = threadIdx.x;
@@ -88,7 +88,7 @@ namespace stim{
 			int byi = blockIdx.y * blockDim.y;
  
 			//copy the portion of the image necessary for this block to shared memory
-			stim::cuda::sharedMemcpy_tex2D(s, in, bxi, byi - kr, 1, 2 * kr + blockDim.y, threadIdx, blockDim);
+			stim::cuda::sharedMemcpy_tex2D<float, unsigned char>(s, in, bxi, byi - kr, 1, 2 * kr + blockDim.y, threadIdx, blockDim);
  
 			//calculate the thread index
 			int ti = threadIdx.y;
@@ -213,6 +213,8 @@ namespace stim{
 			//free allocated memory
 			cudaFree(cuArray);
  
+			cudaDestroyTextureObject(texObj);
+
 		}
  
 		/// Applies a Gaussian blur to a 2D image stored on the CPU
@@ -257,4 +259,4 @@ namespace stim{
 	};
 };
  
-#endif
 \ No newline at end of file
+#endif
@@ -7,7 +7,6 @@
 #include <stim/cuda/sharedmem.cuh>
 #include <stim/cuda/templates/conv2sep.cuh>		//GPU-based separable convolution algorithm
  
-#define pi	3.14159
  
 namespace stim{
 	namespace cuda{
@@ -37,12 +36,14 @@ namespace stim{
  
 			//copy the kernel to the GPU
 			T* gpuKernel0;
+			HANDLE_ERROR(cudaMalloc(&gpuKernel0, kwidth*sizeof(T)));
 			HANDLE_ERROR(cudaMemcpy(gpuKernel0, kernel0, kwidth * sizeof(T), cudaMemcpyHostToDevice));
  
 			//perform the gaussian blur as a separable convolution
 			stim::cuda::tex_conv2sep(out, x, y, texObj, cuArray, gpuKernel0, kwidth, gpuKernel0, kwidth);
  
 			HANDLE_ERROR(cudaFree(gpuKernel0));
+			free(kernel0);
  
 		}
  
@@ -58,7 +59,7 @@ namespace stim{
  
 			//copy the kernel to the GPU
 			T* gpuKernel0;
-			HANDLE_ERROR(cudaMalloc(&gpuKernel0, kwidth * sizeof(T)));
+			HANDLE_ERROR(cudaMalloc(&gpuKernel0, kwidth*sizeof(T)));
 			HANDLE_ERROR(cudaMemcpy(gpuKernel0, kernel0, kwidth * sizeof(T), cudaMemcpyHostToDevice));
  
 			//perform the gaussian blur as a separable convolution
@@ -87,4 +88,4 @@ namespace stim{
 	};
 };
  
-#endif
 \ No newline at end of file
+#endif
+#include <assert.h>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <stdio.h>
+#include <stim/visualization/colormap.h>
+#include <sstream>
+#include <stim/math/vector.h>
+#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
+	///handles from OpenGL and Cuda.
+	///@param DIM_Y			--integer controlling how much memory to allocate.
+	void initArray()
+	{
+			cudaMalloc( (void**) &print, 216*16*sizeof(float));     ///temporary
+	}
+
+	///Deinit function that frees the memery used and releases the texture resource
+	///back to OpenGL.
+	void cleanUP()
+	{
+			cudaFree(print);         ///temporary
+	}  
+
+        __device__
+        float templ(int x)                                                                      
+        {
+                if(x < 16/6 || x > 16*5/6 || (x > 16*2/6 && x < 16*4/6)){
+                        return 1.0;
+                }else{
+                        return 0.0;
+                }
+
+        }
+		
+	///Find the difference of the given set of samples and the template
+	///using cuda acceleration.
+	///@param stim::cuda::cuda_texture t	--stim texture that holds all the references
+	///					  to the data.
+	///@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)
+	{       
+		int x   = threadIdx.x + blockIdx.x * blockDim.x;
+		int y   = threadIdx.y + blockIdx.y * blockDim.y;
+		int idx = y*16+x;
+
+		float valIn             = tex2D<unsigned char>(texIn, x, y);
+		float templa		= templ(x);
+		//print[idx]              = abs(valIn);             ///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)
+	{
+
+		//Bind the Texture in GL and allow access to cuda.
+		tx.MapCudaTexture(texbufferID, texType);
+
+		//initialize the return arrays.
+
+		initArray();
+		
+		int x = 16;
+		int y = 27*8;
+		y = 8* 1089; 
+		int max_threads = stim::maxThreadsPerBlock();
+		//dim3 threads(max_threads, 1);
+		//dim3 blocks(x / threads.x + 1, y);	
+		dim3 numBlocks(1, 1089);
+		dim3 threadsPerBlock(16, 8);
+		//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 << "FromTex.bmp";
+		stim::gpu2image<float>(print, name.str(),16,1089*8,0,1.0);
+	  
+		tx.UnmapCudaTexture();
+		cleanUP();
+	}
+
@@ -13,50 +13,101 @@
 #include "stim/math/vector.h"
 #include "stim/math/rect.h"
 #include "stim/math/matrix.h"
-#include "stim/cuda/cost.h"
+#include "stim/cuda/spider_cost.cuh"
 #include <stim/cuda/cudatools/glbind.h>
-#include <stim/visualization/obj.h>
+#include <stim/cuda/arraymath.cuh>
+#include <stim/cuda/cudatools.h>
+#include <stim/cuda/ivote.cuh>
+#include <stim/visualization/glObj.h>
 #include <vector>
+#include <stim/cuda/branch_detection.cuh>
+#include "../../../volume-spider/fiber.h"
+#include "../../../volume-spider/glnetwork.h"
+//#include <stim/cuda/testKernel.cuh>
+
+//#include <stim/cuda/testKernel.cuh>
  
 #include <iostream>
 #include <fstream>
+#ifdef TESTING
+	#include <iostream>
+	#include <cstdio>
+	#include <ctime>
+#endif
  
  
-
-/* Technically since gl_spider inherits from gl_texture, we could
-	call the init with a path to an image stack, and upload
-	the images while creating the spider (calling init) */
 namespace stim
 {
  
 template<typename T>
-class gl_spider
+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.
 	private:
+
+		//
 		stim::vec<float> p;  	//vector designating the position of the spider.
 		stim::vec<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.
-		std::vector<stim::vec<float> > dV;
-		std::vector<stim::vec<float> > pV;
-		std::vector<stim::vec<float> > mV;
-		//currentTransform
-		stim::matrix<float, 4> cT;
+		std::vector<stim::vec<float> > dV;	//A list of all the direction vectors.
+		std::vector<stim::vec<float> > pV;	//A list of all the position vectors.
+		std::vector<stim::vec<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;
-		stim::vec<float> R;
-		cudaGraphicsResource_t resource;
+		stim::vec<float> S;			//Size of a voxel in the volume.
+		stim::vec<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.
+
+		int numSamples;				//The number of templates in the buffer.
+		int numSamplesPos;
+		int numSamplesMag;
+
+//		float stepsize = 4.0;			//Step size.
+		float stepsize = 3.0;			//Step size.
+		int current_cost;			//variable to store the cost of the current step.
+		
+		
+		//Tracing variables.
+		std::stack< stim::vec<float> > 	seeds;		//seed positions.
+		std::stack< stim::vec<float> > 	seedsvecs;	//seed directions.
+		std::stack< float > 		seedsmags;	//seed magnitudes.
  
-		GLuint dList;
-		GLuint fboID;
-		GLuint texbufferID;
-		int numSamples;
-		float stepsize = 3.0;
-		int current_cost;
+		std::vector< stim::vec<float> > cL;	//Positions of line currently being traced.
+		std::vector< stim::vec<float> > cD;	//Direction of line currently being traced.
+		std::vector< stim::vec<float> > cM;	//Magnitude of line currently being traced.
+
+		stim::glObj<float> sk;			//object to store the skeleton.
+		stim::glnetwork<float> nt;		//object for storing the network.
+
+		stim::vec<float> rev;			//reverse vector;
+		stim::camera camSel;
+		stim::vec<float> ps;
+		stim::vec<float> ups;
+		stim::vec<float> ds;
+		
+
+//--------------------------------------------------------------------------//
+//-------------------------------PRIVATE METHODS----------------------------//
+//--------------------------------------------------------------------------//
  
 		/// Method for finding the best scale for the spider.
 		/// changes the x, y, z size of the spider to minimize the cost
@@ -64,42 +115,43 @@ class gl_spider
 		void
 		findOptimalDirection()
 		{
-			setMatrix();
-			glCallList(dList);
-			int best = getCost();
-			stim::vec<float> next(
+			setMatrix();			//create the transformation matrix.
+			glCallList(dList);		//move the templates to p, d, m.
+			int best = getCost(texbufferID,numSamples);		//find min cost.
+			stim::vec<float> next(		//find 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();
-			//next = (cT*next);
+			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.
 		}
  
-		/// Method for finding the best d for the spider.
-		/// Not sure if necessary since the next p for the spider
+		/// Method for finding the best d (direction) for the spider.
+		/// Not sure if necessary since the next p (position) for the spider
 		/// will be at d * m.
 		void
 		findOptimalPosition()
 		{
-			setMatrix();
-			glCallList(dList+1);
-			int best = getCost();
-			stim::vec<float> next(
- 			pV[best][0],
-			pV[best][1],
-			pV[best][2],
-			1);
-			next = cT*next;	
+			setMatrix();			//create the transformation matrix.
+			glCallList(dList+1);		//move the templates to p, d, m.
+			int best = getCost(ptexbufferID, numSamplesPos);		//find min cost.
+			std::cerr << best << std::endl;
+			stim::vec<float> next(		//find next position.
+ 				pV[best][0],
+				pV[best][1],
+				pV[best][2],
+				1);
+			next = cT*next;			//find next position.
 			setPosition(
 					next[0]*S[0]*R[0],
 					next[1]*S[1]*R[1],
 					next[2]*S[2]*R[2]
-				   );
+				   );			//adjust position.
 		}
  
 		/// Method for finding the best scale for the spider.
@@ -108,33 +160,64 @@ class gl_spider
 		void
 		findOptimalScale()
 		{
-			setMatrix();
-			glCallList(dList+2);
-			int best = getCost();
-			setMagnitude(m[0]*mV[best][0]);
+			setMatrix();			//create the transformation.
+			glCallList(dList+2);		//move the templates to p, d, m.
+			int best = getCost(mtexbufferID, numSamplesMag);		//get best cost.
+			setMagnitude(m[0]*mV[best][0]);	//adjust the magnitude.
 		}
  
+
+		///subject to change.
+		///finds branches.
 		void
 		branchDetection()
 		{
-			Bind();
 			setMatrix();
 			glCallList(dList+3);
-			
-		//	int best = getCost();
+			std::vector< stim::vec<float> > result = find_branch(
+					btexbufferID, GL_TEXTURE_2D, 16, 216);
+			stim::vec<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::vec<float> vec(
+						cylp[0]*S[0]*R[0],
+						cylp[1]*S[1]*R[1],
+						cylp[2]*S[2]*R[2]);
+					stim::vec<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]/2.0;
+					float seedm = m[0];
+// Uncomment for global run 
+/*					stim::vec<float> lSeed = getLastSeed();
+					if(sqrt(pow((lSeed[0] - vec[0]),2)
+					 + pow((lSeed[1] - vec[1]),2) + 
+					 pow((lSeed[2] - vec[2]),2)) > m[0]/4.0
+					 &&  */
+					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);
+					}
+				}
+			}    
  
 		}
  
-
-		
-		void
-		Optimize()
-		{
-			/*find the optimum d and scale */ 
-		}
-
-		
-		
  
 //--------------------------------------------------------------------------//
 //---------------------TEMPLATE CREATION METHODS----------------------------//
@@ -142,14 +225,15 @@ class gl_spider
  
 		///@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.
+		///All points are sampled from a texture.
+		///Stored in a display list.
 		///uses the default d vector <0,0,1>
 		void
-		genDirectionVectors(float solidAngle = 3*M_PI/2)
+		genDirectionVectors(float solidAngle = 5/M_PI*4)
 		{
-			//ofstream file;
-			//file.open("dvectors.txt");
 			//Set up the vectors necessary for Rectangle creation.
-			vec<float> Y(1.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);
@@ -158,12 +242,12 @@ class gl_spider
 			vec<float> d_s = d.cart2sph().norm();
 			vec<float> temp(0,0,0);
 			int dim = (sqrt(numSamples)-1)/2;
-			float p0  	= -M_PI;
-			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0));
-			float dp  	= p0/(2.0*((float)dim + 1.0));
+			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.
  
-			glNewList(dList, GL_COMPILE);
-			//Loop over the space
+			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++){
@@ -192,28 +276,30 @@ class gl_spider
 			glEndList();
 		}
  
-		///@param solidAngle, the size of the arc to sample.
+		///@param float delta, How much the rectangles vary in position.
 		///Method for populating the buffer with the sampled texture.
+		///Objects created are rectangles the with the created positions.
+		///All points are sampled from a texture.
+		///Stored in a display list.
 		///uses the default vector <0,0,0>
 		void
 		genPositionVectors(float delta = 0.4)
 		{
 			//Set up the vectors necessary for Rectangle creation.
-			vec<float> Y(1.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);
  
 			//Set up the variable necessary for vector creation.
 			vec<float> temp(0,0,0);
-			int dim = (sqrt(numSamples)-1)/2;
-			stim::rect<float> samplingPlane =
+			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);
+			float step = 1.0/(dim);			//step size.
  
-			//Loop over the samples, keeping the original p sample
-			//in the center of the resulting texture.
+			//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++){
@@ -240,30 +326,32 @@ class gl_spider
 			glEndList();
 		}
  
-		///@param solidAngle, the size of the arc to sample.
+		///@param float delta, How much the rectangles are allowed to expand.
 		///Method for populating the buffer with the sampled texture.
+		///Objects created are rectangles the with the created sizes.
+		///All points are sampled from a texture.
+		///Stored in a display list.
 		///uses the default m <1,1,0>
 		void
 		genMagnitudeVectors(float delta = 0.70)
-//		genMagnitudeVectors(float delta = 0.50)
 		{
  
 			//Set up the vectors necessary for Rectangle creation.
-			vec<float> Y(1.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);
  
 			//Set up the variable necessary for vector creation.
-			int dim = (sqrt(numSamples)-1)/2;
+			int dim = (sqrt(numSamplesMag)-1)/2;
 			float min 	= 1.0-delta;
 			float max 	= 1.0+delta;
-			float step	= (max-min)/(numSamples-1);
+			float step	= (max-min)/(numSamplesMag-1);
 			float factor;
 			vec<float> temp(0.0,0.0,0.0);
  
 			glNewList(dList+2, GL_COMPILE);
-			for(int i = 0; i < numSamples; i++){
+			for(int i = 0; i < numSamplesMag; i++){
 				//Create linear index
 				factor = (min+step*i)*mag[0];
 				temp = factor;
@@ -280,10 +368,11 @@ class gl_spider
 			}
 			glEndList();
 		}
-		///@param v_x x-coordinate in buffer-space,
-		///@param v_y y-coordinate in buffer-space.
-		///Samples the texturespace and places a sample in the provided coordinates
-		///of bufferspace.
+
+		///@param float v_x x-coordinate in buffer-space,
+		///@param float v_y y-coordinate in buffer-space.
+		///Samples the texture space.
+ 		///places a sample in the provided coordinates of bufferspace.
 		void
 		UpdateBuffer(float v_x, float v_y)
 		{	
@@ -361,8 +450,37 @@ class gl_spider
 //--------------------------------GL METHODS--------------------------------//
 //--------------------------------------------------------------------------//
  
-		///@param width sets the width of the buffer.
-		///@param height sets the height of the buffer.
+		///@param uint width sets the width of the buffer.
+		///@param uint height sets the height of the buffer.
+		///@param GLuint &textureID gives the texture ID of the texture to be initialized.
+		///@param GLuint &framebufferID gives the buffer ID of the texture to be initialized.
+		///Function for setting up the 2D buffer that stores the samples.
+		///Initiates and sets parameters.
+		void
+		GenerateFBO(unsigned int width, unsigned int height, GLuint &textureID, GLuint &framebufferID)
+		{
+			glGenFramebuffers(1, &framebufferID);
+			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
+			int numChannels = 1;
+			unsigned char* texels = new unsigned char[width * height * numChannels];
+			glGenTextures(1, &textureID);
+			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_MAG_FILTER, GL_LINEAR);
+			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,
+				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, texels);   
+			delete[] texels;
+			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
+			glBindTexture(GL_TEXTURE_2D, 0);
+			CHECK_OPENGL_ERROR
+		}
+
+		///@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)
@@ -373,6 +491,8 @@ class gl_spider
 			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);
@@ -382,39 +502,150 @@ class gl_spider
 			delete[] texels;
 			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
 			glBindTexture(GL_TEXTURE_2D, 0);
+			CHECK_OPENGL_ERROR
 		}
  
  
-		///Method for using the gl manipulation to alighn templates from
+		///Method for using the gl manipulation to align templates from
 		///Template space (-0.5 0.5) to Texture space (0.0, 1.0),
 		///Based on the p of the spider in real space (arbitrary).
+		///All transformation happen in glMatrixMode(GL_TEXTURE).
 		void setMatrix()
 		{
-			float curTrans[16];
-			stim::vec<float> rot = getRotation(d);
+			float curTrans[16];			//array to store the matrix values.
+			stim::vec<float> rot = getRotation(d);	//get the rotation parameters for the current direction vector.
 			glMatrixMode(GL_TEXTURE);
 			glLoadIdentity();
-			glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
-
  
+			//Scale by the voxel size and number of slices.
+			glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
+			//translate to the current position of the spider in the texture.
 			glTranslatef(p[0],
 				     p[1],
 				     p[2]);
-
+			//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]);
-
+			//get and store the current transformation matrix for later use.
 			glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
 			cT.set(curTrans);
-//			printTransform();
+		//	printTransform();
  
 			CHECK_OPENGL_ERROR
+			//revert back to default gl mode.
 			glMatrixMode(GL_MODELVIEW);
 		}
+
+		///Method for controling the buffer and texture binding.
+		///Clears the buffer upon binding.
+		void
+		Bind()
+		{
+			float len = 8.0;
+			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
+			glFramebufferTexture2D(
+				GL_FRAMEBUFFER,
+				GL_COLOR_ATTACHMENT0,
+				GL_TEXTURE_2D,
+				texbufferID,
+				0);
+			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
+			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
+			glDrawBuffers(1, DrawBuffers);
+			glBindTexture(GL_TEXTURE_2D, texbufferID);
+			glClearColor(1,1,1,1);
+			glClear(GL_COLOR_BUFFER_BIT);
+			glMatrixMode(GL_PROJECTION);
+			glLoadIdentity();
+			glMatrixMode(GL_MODELVIEW);
+			glLoadIdentity();
+			glViewport(0,0,2.0*len, numSamples*len);
+			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
+			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.
+		///@param GLuint &framebufferID, framebuffer used for storage.
+		///@param int nSamples, number of rectanges to create. 
+		void
+		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples)
+		{
+			float len = 8.0;
+			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
+			glFramebufferTexture2D(
+				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);
+//			glClearColor(1,1,1,1);
+//			glClear(GL_COLOR_BUFFER_BIT);
+			glMatrixMode(GL_PROJECTION);
+			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);
+
+			CHECK_OPENGL_ERROR
+		}
+		
+		///Unbinds all texture resources.
+		void
+		Unbind()
+		{
+			//Finalize GL_buffer
+			glBindTexture(GL_TEXTURE_3D, 0);                      
+			CHECK_OPENGL_ERROR
+			glBindTexture(GL_TEXTURE_2D, 0);                      
+			CHECK_OPENGL_ERROR
+			glBindFramebuffer(GL_FRAMEBUFFER, 0);
+			CHECK_OPENGL_ERROR
+			glDisable(GL_TEXTURE_3D);
+			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
+			#ifdef TESTING
+				start = std::clock();
+			#endif
+			findOptimalDirection();
+			findOptimalPosition();
+			findOptimalScale();
+			Unbind();
+			Bind(btexbufferID, bfboID, 27);
+			branchDetection();
+			Unbind();
+
+			#ifdef TESTING
+				duration_sampling = duration_sampling +
+					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
+				num_sampling = num_sampling + 1.0; 
+			#endif
+			return current_cost;
+		}
+
  
  
  
@@ -422,95 +653,150 @@ class gl_spider
 //--------------------------------CUDA METHODS------------------------------//
 //--------------------------------------------------------------------------//
  
-		/// Method for registering the texture with Cuda for shared
-		///	access.
-		void
-		createResource()
-		{
-			HANDLE_ERROR(
-				cudaGraphicsGLRegisterImage(
-					 &resource,
-				 	texbufferID,
-				 	GL_TEXTURE_2D,
-				 	//CU_GRAPHICS_REGISTER_FLAGS_NONE)
-					cudaGraphicsMapFlagsReadOnly)
-			);
-		} 
-		
-		///Method for freeing the texture from Cuda for gl access.
-		void
-		destroyResource()
+
+		///Entry-point into the cuda code for calculating the cost of a given samples array (in texture form)
+		///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()
 		{
-			HANDLE_ERROR(
-				cudaGraphicsUnregisterResource(resource)
-			);		
+			#ifdef TESTING
+				start = std::clock();
+			#endif
+			stim::vec<int> cost = 
+				stim::cuda::get_cost(texbufferID, GL_TEXTURE_2D, numSamples);
+			cudaDeviceSynchronize();
+			#ifdef TESTING
+				duration_cuda = duration_cuda +
+					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
+				num_cuda = num_cuda + 1.0;
+			#endif
+			current_cost = cost[1];
+			return cost[0];
 		}
  
-		///Entry-point into the cuda code for calculating the cost
-		///	of a given samples array (in texture form) 
 		int
-		getCost()
+		getCost(GLuint tID, int n)
 		{
-			createResource();
-			stim::vec<int> cost = 	get_cost(resource, numSamples);
-			destroyResource();
-//			if (cost[1] >= 80)
-//				exit(0);
+			#ifdef TESTING
+				start = std::clock();
+			#endif
+			stim::vec<int> cost = 
+				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n);
+			cudaDeviceSynchronize();
+			#ifdef TESTING
+				duration_cuda = duration_cuda +
+					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
+				num_cuda = num_cuda + 1.0;
+			#endif
 			current_cost = cost[1];
 			return cost[0];
 		}
  
 	public:
+
+		///ininializes the cuda device and environment.
+		void
+		initCuda()
+		{	
+			stim::cudaSetDevice();
+			//GLint max;
+			//glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
+			//std::cout << max << std::endl;
+		}
+
+		//horizonal rectangle forming the spider.
 		stim::rect<float> hor;
+		//vectical rectangle forming the spider.
 		stim::rect<float> ver;	
  
+		//Testing and Timing variables.
+		#ifdef TESTING
+			std::clock_t start;
+			double duration_sampling = 0.0;
+			double duration_cuda = 0.0;
+			double num_sampling = 0.0;
+			double num_cuda = 0.0;
+		#endif
+
 //--------------------------------------------------------------------------//
 //-----------------------------CONSTRUCTORS---------------------------------//
 //--------------------------------------------------------------------------//
  
  
-		///@param samples, the number of samples this spider is going to use.
-		///best results if samples is can create a perfect root.
+		///@param int samples, the number of samples this spider is going to use.
+		///Best results if samples is can create a perfect root.
 		///Default Constructor
 		gl_spider
-		(int samples = 1089)
+		(int samples = 1089, int samplespos = 441,int samplesmag = 144)
 		{
 			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);
-			//setPosition(0.0,0.0,0.0);
-			//setDirection(0.0,0.0,1.0);
-			//setMagnitude(1.0);
 			numSamples = samples;
+			numSamplesPos = samplespos;
+			numSamplesMag = samplesmag;
 		}
  
-		///temporary constructor for convenience, will be removed in further updates.	
+		///Position constructor: floats.
+		///@param float pos_x, position x.
+		///@param float pos_y, position y.
+		///@param float pos_z, position z.
+		///@param float dir_x, direction x.
+		///@param float dir_y, direction y.
+		///@param float dir_z, direction z.
+		///@param float mag_x, size of the vector.
+		///@param int 	samples, number of templates this spider is going to use.
 		gl_spider
 		(float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z,
-			float mag_x, int numSamples = 1089)
+			float mag_x, int numsamples = 1089, int numsamplespos = 441, int numsamplesmag =144)
 		{
 			p = vec<float>(pos_x, pos_y, pos_z);
 			d = vec<float>(dir_x, dir_y, dir_z);
 			m = vec<float>(mag_x, mag_x, mag_x);
 			S = vec<float>(1.0,1.0,1.0);
 			R = vec<float>(1.0,1.0,1.0);
-			//setPosition(pos_x, pos_y, pos_z);
-			//setDirection(dir_x, dir_y, dir_z);
-			//setMagnitude(mag_x);
-		
+			numSamples = numsamples;
+			numSamplesPos = numsamplespos;
+			numSamplesMag = numsamplesmag;
 		}
-	
+
+		///Position constructor: vecs of floats.
+		///@param stim::vec<float> pos, position.
+		///@param stim::vec<float> dir, direction.
+		///@param float mag, size of the vector.
+		///@param int 	samples, number of templates this spider is going to use.
+		gl_spider
+		(stim::vec<float> pos, stim::vec<float> dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144)
+		{
+			p = pos;
+			d = dir;
+			m = vec<float>(mag, mag, mag);
+			S = vec<float>(1.0,1.0,1.0);
+			R = vec<float>(1.0,1.0,1.0);
+			numSamples = samples;
+			numSamplesPos = samplesPos;
+			numSamplesMag = samplesMag;
+		}
+
+		///destructor	
 		~gl_spider
 		(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);
 		}
  
-		///@param GLuint id texture that is going to be sampled.
+		///@param GLuint id, texture that is going to be sampled.
 		///Attached the spider to the texture with the given GLuint ID.
 		///Samples in the default d acting as the init method.
 		///Also acts an init.	
@@ -518,16 +804,26 @@ class gl_spider
 		attachSpider(GLuint id)
 		{
 			texID = id;
-			GenerateFBO(16, numSamples*8);
+			//GenerateFBO(16, numSamples*8);
+			GenerateFBO(16, numSamples*8, texbufferID, fboID); 
+			GenerateFBO(16, numSamplesPos*8, ptexbufferID, pfboID); 
+			GenerateFBO(16, numSamplesMag*8, mtexbufferID, mfboID); 
+			GenerateFBO(16, 216, btexbufferID, bfboID); 
 			setDims(0.6, 0.6, 1.0);
 			setSize(512.0, 512.0, 426.0);
 			setMatrix();
 			dList = glGenLists(3);
 			glListBase(dList);
-			Bind();
-			genDirectionVectors(5*M_PI/4);
-			genPositionVectors();
-			genMagnitudeVectors();
+			Bind(texbufferID, fboID, numSamples);
+				genDirectionVectors(5*M_PI/4);
+			Unbind();
+			Bind(ptexbufferID, pfboID, numSamplesPos);
+				genPositionVectors();
+			Unbind();
+			Bind(mtexbufferID, mfboID, numSamplesMag);
+				genMagnitudeVectors();
+			Unbind();
+			Bind(btexbufferID, bfboID, 27);
 			DrawCylinder();
 			Unbind();
 		}
@@ -556,7 +852,7 @@ class gl_spider
 			return m;
 		}
  
-		///@param vector pos, the new p.
+		///@param stim::vec<float> pos, the new p.
 		///Sets the p vector to input vector pos.
 		void
 		setPosition(vec<float> pos)
@@ -564,9 +860,9 @@ class gl_spider
 			p = pos;
 		}
  
-		///@param x x-coordinate.
-		///@param y y-coordinate.
-		///@param z z-coordinate.
+		///@param float x x-coordinate.
+		///@param float y y-coordinate.
+		///@param float z z-coordinate.
 		///Sets the p vector to the input float coordinates x,y,z.
 		void
 		setPosition(float x, float y, float z)
@@ -576,7 +872,7 @@ class gl_spider
 			p[2] = z;
 		}
  
-		///@param vector dir, the new d.
+		///@param stim::vec<float> dir, the new d.
 		///Sets the d vector to input vector dir.
 		void
 		setDirection(vec<float> dir)
@@ -584,9 +880,9 @@ class gl_spider
 			d = dir;
 		}
  
-		///@param x x-coordinate.
-		///@param y y-coordinate.
-		///@param z z-coordinate.
+		///@param stim::vec<float> x x-coordinate.
+		///@param stim::vec<float> y y-coordinate.
+		///@param stim::vec<float> z z-coordinate.
 		///Sets the d vector to the input float coordinates x,y,z.
 		void
 		setDirection(float x, float y, float z)
@@ -596,7 +892,7 @@ class gl_spider
 			d[2] = z;
 		}
  
-		///@param vector dir, the new d.
+		///@param stim::vec<float> dir, the new d.
 		///Sets the m vector to the input vector mag.	
 		void
 		setMagnitude(vec<float> mag)
@@ -605,17 +901,19 @@ class gl_spider
 			m[1] = mag[0];
 		}
  
-		///@param mag size of the sampled region.
+		///@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[2] = mag;
 		}
  
-
+		///@param float x, voxel size in the x direction.
+		///@param float y, voxel size in the y direction.
+		///@param float z, voxel size in the z direction.
+		///Sets the voxel sizes in each direction. necessary for non-standard data.
 		void
 		setDims(float x, float y, float z)
 		{
@@ -624,6 +922,18 @@ class gl_spider
 			S[2] = z;
 		}
  
+		///@param stim::vec<float> Dims, voxel size.
+		///Sets the voxel sizes in each direction. necessary for non-standard data.
+		void
+		setDims(stim::vec<float> Dims)
+		{
+			S = Dims;
+		}
+
+		///@param float x, size of the data in the x direction.
+		///@param float y, size of the data in the y direction.
+		///@param float z, size of the data in the z direction.
+		///Sets the data volume sizes in each direction.
 		void
 		setSize(float x, float y, float z)
 		{
@@ -631,10 +941,19 @@ class gl_spider
 			R[1] = y;
 			R[2] = z;
 		}
+
+		///@param stim::vec<float> Dims, size of the volume.
+		///Sets the data volume sizes in each direction.
+		void
+		setSize(stim::vec<float> Siz)
+		{
+			S = Siz;
+		}
  
-		///@param dir, the vector to which we are rotating
-		///given a vector to align to, finds the required
-		///axis and angle for glRotatef
+		///@param stim::vec<float> dir, the vector to which we are rotating.
+		///given a vector to align to, finds the required axis and angle for glRotatef.
+		///rotates from 0.0, 0.0, 1.0 to dir.
+		///return is in degrees for use with glRotatef.
 		stim::vec<float>
 		getRotation(stim::vec<float> dir)
 		{
@@ -655,56 +974,156 @@ class gl_spider
 			}
 			return out;
 		}
+
+		///@param stim::vec<float> pos, the position of the seed to be added.
+		///Adds a seed to the seed list.
+		///Assumes that the coordinates passes are in tissue space.
+		void
+		setSeed(stim::vec<float> pos)
+		{
+			seeds.push(pos);
+		}
+
+		///@param stim::vec<float> dir, the direction of the seed to be added.
+		///Adds a seed to the seed directions list.
+		void
+		setSeedVec(stim::vec<float> dir)
+		{
+			seedsvecs.push(dir);
+		}
+
+		///@param float mag, the size of the seed to be added.
+		///Adds a seed to the seed list.
+		///Assumes that the coordinates passes are in tissue space.
+		void
+		setSeedMag(float mag)
+		{
+			seedsmags.push(mag);
+		}
+
+
+		///@param float x, x-position of the seed to be added.
+		///@param float y, y-position of the seed to be added.
+		///@param float z, z-position of the seed to be added.
+		///Adds a seed to the seed list.
+		///Assumes that the coordinates passes are in tissue space.
+		void
+		setSeed(float x, float y, float z)
+		{
+			seeds.push(stim::vec<float>(x, y, z));
+		}
+
+		///@param float x, x-direction of the seed to be added.
+		///@param float y, y-direction of the seed to be added.
+		///@param float z, z-direction of the seed to be added.
+		///Adds a seed to the seed directions list.
+		void
+		setSeedVec(float x, float y, float z)
+		{
+			seedsvecs.push(stim::vec<float>(x, y, z));
+		}
  
-		///Function to get back the framebuffer Object attached to the spider.
-		///For external access.
-		GLuint
-		getFB()
+		///Method to get the top of the seed positions stack.
+		stim::vec<float> 
+		getLastSeed()
+		{
+			stim::vec<float> tp = seeds.top();
+			return tp;
+		}
+
+		///Method to get the top of the seed direction stack.
+		stim::vec<float> 
+		getLastSeedVec()
+		{
+			stim::vec<float> tp = seedsvecs.top();
+			return tp;
+		}
+
+		///Method to get the top of the seed magnitude stack.
+		float
+		getLastSeedMag()
 		{
-			return fboID;
+			float tp = seedsmags.top();
+			return tp;
 		}
  
-		///Method for controling the buffer and texture binding in order to properly
-		///do the render to texture.
+		///deletes all data associated with the last seed.
 		void
-		Bind()
+		popSeed()
 		{
-			float len = 8.0;
-			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
-			glFramebufferTexture2D(
-				GL_FRAMEBUFFER,
-				GL_COLOR_ATTACHMENT0,
-				GL_TEXTURE_2D,
-				texbufferID,
-				0);
-			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
-			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
-			glDrawBuffers(1, DrawBuffers);
-			glBindTexture(GL_TEXTURE_2D, texbufferID);
-			glClearColor(1,1,1,1);
-			glClear(GL_COLOR_BUFFER_BIT);
-			glMatrixMode(GL_PROJECTION);
-			glLoadIdentity();
-			glMatrixMode(GL_MODELVIEW);
-			glLoadIdentity();
-			glViewport(0,0,2.0*len, numSamples*len);
-			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
-			glEnable(GL_TEXTURE_3D);
-			glBindTexture(GL_TEXTURE_3D, texID);
+			seeds.pop();
+			seedsvecs.pop();
+			seedsmags.pop();
+		}
+		
+		///returns the seeds position stack.
+		std::stack<stim::vec<float> >
+		getSeeds()
+		{
+			return seeds;
+		}
  
-			CHECK_OPENGL_ERROR
+		///returns true if all seed stacks are empty, else false.
+		bool
+		Empty()
+		{
+			//return (seeds.empty() && seedsvecs.empty() && seedsmags.empty());
+			return (seeds.empty() && seedsvecs.empty());
+		}
+
+		///@param std::string file:file with variables to populate the seed stacks.
+		///Adds a seed to the seed list, including the position, direction and magnitude.
+		///Assumes that the coordinates passes are in tissue space.
+		void
+		setSeeds(std::string file)
+		{
+			std::ifstream myfile(file.c_str());
+			string line;
+			if(myfile.is_open())
+           		{
+                   		while (getline(myfile, line))
+                   		{
+                           		float x, y, z, u, v, w, m;
+                           		myfile >> x >> y >> z >> u >> v >> w >> m;
+					setSeed(x, y , z);
+					setSeedVec(u, v, w);
+					setSeedMag(m);
+                   		}
+                   	myfile.close();
+	          	} else {
+        	                  std::cerr<<"failed" << std::endl;
+			}
 		}
  
-		///Method for Unbinding all of the texture resources
+		///Saves the network to a file.
 		void
-		Unbind()
+		saveNetwork(std::string name)
 		{
-			//Finalize GL_buffer
-			glBindTexture(GL_TEXTURE_3D, 0);                      
-			glDisable(GL_TEXTURE_3D);
-			glBindFramebuffer(GL_FRAMEBUFFER,0);
-			glBindTexture(GL_TEXTURE_2D, 0);
+			sk.save(name);
+		}
+
+		///returns a COPY of the entire stim::glObj object.
+		stim::glObj<float>
+		getNetwork()
+		{
+			return sk;
+		}
+
+		///returns a COPY of the entire stim::glnetwork object.
+		stim::glnetwork<T>
+		getGLNetwork()
+		{
+			return nt;
+		}
+		
+		///Function to get back the framebuffer Object attached to the spider.
+		///For external access.
+		GLuint
+		getFB()
+		{
+			return bfboID;
 		}
+
 //--------------------------------------------------------------------------//
 //-----------------------------TEMPORARY METHODS----------------------------//
 //--------------------------------------------------------------------------//
@@ -726,12 +1145,26 @@ class gl_spider
 		int
 		Step()
 		{
-		//	Bind();
-			findOptimalDirection();
-			findOptimalPosition();
-			findOptimalScale();
-			branchDetection();
-		//	Unbind();
+			Bind(texbufferID, fboID, numSamples);
+			CHECK_OPENGL_ERROR
+			#ifdef TESTING
+				start = std::clock();
+			#endif
+				findOptimalDirection();
+			Unbind();
+			Bind(ptexbufferID, pfboID, numSamplesPos);
+				findOptimalPosition();
+			Unbind();
+			Bind(mtexbufferID, mfboID, numSamplesMag);
+				findOptimalScale();
+			Unbind();
+			CHECK_OPENGL_ERROR
+
+			#ifdef TESTING
+				duration_sampling = duration_sampling +
+					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
+				num_sampling = num_sampling + 1.0; 
+			#endif
 			return current_cost;
 		}
  
@@ -742,16 +1175,6 @@ class gl_spider
 			std::cout << cT << std::endl;
 		}
  
-		/* Method for initializing the cuda devices, necessary only
-			there are multiple cuda devices */
-		void
-		initCuda()
-		{	
-			stim::cudaSetDevice();
-			//GLint max;
-			//glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
-			//std::cout << max << std::endl;
-		}
  
 //--------------------------------------------------------------------------//
 //-----------------------------EXPERIMENTAL METHODS-------------------------//
@@ -760,36 +1183,410 @@ class gl_spider
 		void
 		DrawCylinder()
 		{	 
-			 Bind();
 			 glNewList(dList+3, GL_COMPILE);
 			 float z0 = -0.5; float z1 = 0.5; float r0 = 0.5;
 			 float x,y;
-			 float xold = 0.5; float yold = 0.5;
-			 float step = 360.0/numSamples;
-		 	 int j = 0;
+			 float xold = 0.5; float yold = 0.0;
+			 float step = 360.0/numSamples*32;
+			 //float step = 360.0/8.0;
 			 glEnable(GL_TEXTURE_3D);
 			 glBindTexture(GL_TEXTURE_3D, texID);
 			 glBegin(GL_QUAD_STRIP);
+			 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);
 					 glTexCoord3f(x,y,z0); 
-					 glVertex2f(0.0, j*0.1+0.1);
+					 glVertex2f(0.0, j*6.4+6.4);
+//					 glVertex2f(0.0, j*27.0+27.0);
 					 glTexCoord3f(x,y,z1); 
-					 glVertex2f(16.0, j*0.1+0.1);
+					 glVertex2f(16.0, j*6.4+6.4);
+//					 glVertex2f(16.0, j*27.0+27.0);
 					 glTexCoord3f(xold,yold,z1); 
-					 glVertex2f(16.0, j*0.1); 
+					 glVertex2f(16.0, j*6.4); 
+//					 glVertex2f(16.0, j*27.0); 
 					 glTexCoord3f(xold,yold,z0); 
-					 glVertex2f(0.0, j*0.1);
+					 glVertex2f(0.0, j*6.4);
+//					 glVertex2f(0.0, j*27.0);
 					 xold=x;
 					 yold=y;
 					 j++;
 				}
 		      	 glEnd();  
 			 glEndList();
-			 Unbind();
 		}
+		
+
+		///@param min_cost the cost value used for tracing
+		///traces out each seedpoint in the seeds queue to completion in both directions.
+		void
+		trace(int min_cost)
+		{	
+			Bind();
+			rev = stim::vec<float>(0.0,0.0,1.0);
+			bool sEmpty = true;
+			float lastmag = 16.0;;
+			while(!seeds.empty())
+			{
+				//clear the currently traced line and start a new one.
+				cL.clear();
+				cM.clear();
+				sk.Begin(stim::OBJ_LINE);
+				stim::vec<float> curSeed = seeds.top();
+//				std::cout << "The current seeds is " << curSeed << std::endl;
+				stim::vec<float> curSeedVec = seedsvecs.top();
+				float curSeedMag = seedsmags.top();
+				seeds.pop();
+				seedsvecs.pop();
+				seedsmags.pop();
+//				std::cout << "The current seed Vector is " << curSeedVec << std::endl;
+				setPosition(curSeed);
+				setDirection(curSeedVec);
+				cL.push_back(curSeed);
+				cM.push_back(curSeedMag);
+				sk.createFromSelf(GL_SELECT);
+				traceLine(min_cost);
+
+				sk.rev();
+		//		std::cout << "reversed" << std::endl;
+				std::reverse(cL.begin(), cL.end());
+				std::reverse(cM.begin(), cM.end());
+				setPosition(curSeed);
+				setDirection(-rev);
+				setMagnitude(16.0);
+				sk.createFromSelf(GL_SELECT);
+				traceLine(min_cost);
+				sk.End();
+			}
+			Unbind();
+		}
+
+		///@param min_cost the cost value used for tracing
+		///traces the seedpoint passed to completion in one directions.
+		void
+		traceLine(int min_cost)
+		{
+			stim::vec<float> pos;
+			stim::vec<float> mag;
+			int h;
+			bool started = false;
+			bool running = true;
+			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
+			while(running)
+			{
+				int cost = Step();
+				if (cost > min_cost){
+					running = false;
+					break;
+				} else {
+					//Have we found an edge?
+					pos = getPosition();
+					if(pos[0] > size[0] || pos[1] > size[1]
+					 || pos[2] > size[2] || pos[0] < 0
+					 || pos[1] < 0 || pos[2] < 0)
+					{
+//					       std::cout << "Found Edge" << std::endl;
+						running = false;
+						break;
+					}
+					//If this is the first step in the trace,
+					// save the direction
+					//(to be used later to trace the fiber in the opposite direction)
+					if(started == false){
+						rev = -getDirection();
+						started = true;
+					}
+//					std::cout << i << p << std::endl;
+					m = getMagnitude();
+					//Has the template size gotten unreasonable?
+					if(m[0] > 75 || m[0] < 1){
+//						std::cout << "Magnitude Limit" << std::endl;
+						running = false;
+						break;
+					}
+					else
+					{
+						h = selectObject(pos, getDirection(), m[0]);
+						//Have we hit something previously traced?
+						if(h != -1){
+							std::cout << "I hit a line" << h << std::endl;
+							running = false;
+							break;
+						}
+						else {          
+							cL.push_back(stim::vec<float>(p[0], p[1],p[2]));//
+							sk.TexCoord(m[0]);
+							sk.Vertex(p[0], p[1], p[2]);
+							Bind(btexbufferID, bfboID, 27);
+							CHECK_OPENGL_ERROR
+							branchDetection();
+							CHECK_OPENGL_ERROR
+							Unbind();
+							CHECK_OPENGL_ERROR
+						}
+				 	}
+                 		}
+         		}
+		}	
+
+
+		int
+		selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag) 
+		{
+		//Define the varibles and turn on Selection Mode
+
+			float s = 3.0;
+			GLuint selectBuf[2048];
+			GLint hits;
+			glSelectBuffer(2048, selectBuf);
+			glDisable(GL_CULL_FACE);
+			(void) glRenderMode(GL_SELECT);
+		//Init Names stack
+
+			glInitNames();
+			glPushName(1);
+
+			CHECK_OPENGL_ERROR
+		//What would that vessel see in front of it.
+				camSel.setPosition(loc);
+				camSel.setFocalDistance(mag/s);
+				camSel.LookAt((loc[0]+dir[0]*mag/s),
+					 (loc[1]+dir[1]*mag/s),
+					 (loc[2]+dir[2]*mag/s));
+				ps  = camSel.getPosition();
+				ups = camSel.getUp();
+				ds  = camSel.getLookAt();
+				glMatrixMode(GL_PROJECTION);
+				glPushMatrix();
+				glLoadIdentity();
+				glOrtho(-mag/s/2.0, mag/s/2.0, -mag/s/2.0, mag/s/2.0, 0.0, mag/s/2.0);
+				glMatrixMode(GL_MODELVIEW);
+				glPushMatrix();
+				glLoadIdentity();
+				
+				CHECK_OPENGL_ERROR
+				gluLookAt(ps[0], ps[1], ps[2],
+					 ds[0], ds[1], ds[2],
+					 ups[0], ups[1], ups[2]);
+
+//				sk.Render();
+				nt.Render();								
+
+				CHECK_OPENGL_ERROR
+
+
+//				glLoadName((int) sk.numL());
+				glLoadName(nt.sizeE());
+
+//				sk.RenderLine(cL);
+				nt.RenderLine(cL);	
+
+//				glPopName();
+				glFlush();
+
+				glMatrixMode(GL_PROJECTION);
+				glPopMatrix();
+				glMatrixMode(GL_MODELVIEW);
+				CHECK_OPENGL_ERROR
+				glPopMatrix();
+
+		//	glEnable(GL_CULL_FACE);
+			hits = glRenderMode(GL_RENDER);
+			int found_hits = processHits(hits, selectBuf);
+			return found_hits;
+		}
+
+		//Given a size of the array (hits) and the memory holding it (buffer)
+		//returns whether a hit tool place or not.
+		int
+		processHits(GLint hits, GLuint buffer[])
+		{
+			GLuint names, *ptr;
+			//printf("hits = %u\n", hits);
+			ptr = (GLuint *) buffer;
+		//	for (int i = 0; i < hits; i++) { /*  for each hit  */
+			names = *ptr;
+		//		printf (" number of names for hit = %u\n", names);
+			ptr++;
+			ptr++; //Skip the minimum depth value.  
+			ptr++; //Skip the maximum depth value.
+		//		printf ("   the name is ");                                                        
+		//		for (int j = 0; j < names; j++) {     /*  for each name */
+		//			printf ("%u ", *ptr); ptr++;
+		//		}
+		//		printf ("\n");
+	//		}
+
+			
+			if(hits == 0)
+			{
+				return -1;
+			}
+			else
+			{
+//				printf ("%u ", *ptr);
+				return *ptr;
+			}
+		}
+
+		void
+		clearCurrent()
+		{
+			cL.clear();
+			cM.clear();
+		}	
+	
+		void
+		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec<float> spos, 
+				stim::vec<float> smag, stim::vec<float> sdir)
+		{
+                        std::vector<stim::vec<float> > ce = in.first.centerline();                
+                        std::vector<stim::vec<float> > cm = in.first.centerlinemag();
+			//if the fiber is longer than 2 steps (the number it takes to diverge)
+			if(ce.size() > 3)
+			{	
+				//if we did not hit a fiber
+				if(in.second == -1)
+				{
+					spos[0] = spos[0]-sdir[0]*smag[0]/2.;
+					spos[1] = spos[1]-sdir[1]*smag[0]/2.;
+					spos[2] = spos[2]-sdir[2]*smag[0]/2.;
+					int h = selectObject(spos, -sdir, smag[0]);
+					//did we start with a fiber?
+					if(h != -1)
+						nt.addEdge(ce, cm, h, -1);
+					else
+						nt.addEdge(ce, cm, -1, -1);
+				}
+				//if we hit a fiber?
+				else if(in.second != -1)
+				{
+					nt.addEdge(ce,cm,-1, in.second);
+					spos[0] = spos[0]-sdir[0]*smag[0]/2.;
+					spos[1] = spos[1]-sdir[1]*smag[0]/2.;
+					spos[2] = spos[2]-sdir[2]*smag[0]/2.;
+					int h = selectObject(spos, -sdir, smag[0]);
+					//did start with a fiber?
+					if(h != -1){	
+			//			std::cout << "got here double" << smag.str() << std::endl;
+						nt.addEdge(ce,cm, h, in.second);	
+					}
+				}
+			}		
+		}
+
+
+		void
+		printSizes()
+		{
+			std::cout << nt.sizeE() << " edges " << std::endl;
+			std::cout << nt.sizeV() << " nodes " << std::endl;
+			
+		}
+
+		std::pair<stim::fiber<float>, int >
+		traceLine(stim::vec<float> pos, stim::vec<float> mag, int min_cost)
+		{
+			//starting (seed) position and magnitude.
+			stim::vec<float> spos = getPosition();
+			stim::vec<float> smag = getMagnitude();
+			stim::vec<float> sdir = getDirection();	
+
+			Bind();
+			sk.Begin(stim::OBJ_LINE);
+
+
+//			sk.createFromSelf(GL_SELECT);
+			nt.createFromSelf(GL_SELECT);
+
+			cL.push_back(pos);
+			cM.push_back(mag);
+
+//			setPosition(pos);
+//			setMagnitude(mag);
+			int h;
+			bool started = false;
+			bool running = true;
+			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
+			while(running)
+			{
+				int cost = Step();
+				if (cost > min_cost){
+					running = false;
+					sk.End();
+					pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
+					addToNetwork(a, spos, smag, sdir);
+					return a;
+					break;
+				} else {
+					//Have we found the edge of the map?
+					pos = getPosition();
+					if(pos[0] > size[0] || pos[1] > size[1]
+					 || pos[2] > size[2] || pos[0] < 0
+					 || pos[1] < 0 || pos[2] < 0)
+					{
+//					       std::cout << "Found Edge" << std::endl;
+						running = false;
+						sk.End();
+						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
+						addToNetwork(a, spos, smag, sdir);
+						return a;
+						break;
+					}
+					//If this is the first step in the trace,
+					// save the direction
+					//(to be used later to trace the fiber in the opposite direction)
+					if(started == false){
+						rev = -getDirection();
+						started = true;
+					}
+//					std::cout << i << p << std::endl;
+					//Has the template size gotten unreasonable?
+					mag = getMagnitude();
+					if(mag[0] > 75 || mag[0] < 1){
+//						std::cout << "Magnitude Limit" << std::endl;
+						running = false;
+						sk.End();
+						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
+						addToNetwork(a, spos, smag, sdir);
+						return a;
+						break;
+					}
+					else
+					{
+						h = selectObject(p, getDirection(), m[0]);
+						//Have we hit something previously traced?
+						if(h != -1){
+							running = false;
+							sk.End();
+							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
+							addToNetwork(a, spos, smag, sdir);
+							return a;
+							break;
+						}
+						else {  
+							cL.push_back(stim::vec<float>(p[0], p[1],p[2]));
+							cM.push_back(stim::vec<float>(m[0], m[0]));
+//							cM.push_back(m[0]);
+
+							sk.TexCoord(m[0]);
+							sk.Vertex(p[0], p[1], p[2]);
+							Bind(btexbufferID, bfboID, 27);
+							CHECK_OPENGL_ERROR
+							branchDetection();
+							CHECK_OPENGL_ERROR
+							Unbind();
+							CHECK_OPENGL_ERROR
+							
+						}
+				 	}
+                 		}
+         		}
+		}
+
+			
+
 };
 }
 #endif
@@ -31,10 +31,9 @@ template&lt;typename T&gt;
 class gl_texture : public virtual image_stack<T>
 {
 	private:
-		///Method: setTextureType
 		///	Sets the internal texture_type, based on the data
 		///	size. Either 3D, 2D, 1D textures.
-
+		
 		void
 		setTextureType()
 		{
@@ -60,16 +59,14 @@ class gl_texture : public virtual image_stack&lt;T&gt;
  
 	public:
  
-		///Method: Basic Constructor
-		///	Creates an instance of the gl_texture object.
+		///default constructor
 		gl_texture()
 		{
  
 		}
  
-		///Method: Path Constructor
-		///@param string file_path path to the directory with the files.
-		///	Creates an instance of the gl_texture object with a path to the data.
+		///@param string path to the directory with the image files.
+		///Creates an instance of the gl_texture object with a path to the data.
  
 		gl_texture(std::string file_path)
 		{
@@ -77,8 +74,8 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			image_stack<T>::load_images(path.append("/*.jpg"));
 			setTextureType();
 		}
-		///Method:getSize
-		///returns the dimentions of 
+
+		///returns the dimentions of the data in the x, y, z directions. 
 		vec<int>
 		getSize()
 		{
@@ -86,7 +83,6 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			return size;
 		}
  
-		///Method:setTexParam
 		///@param GLint interp     --GL_LINEAR, GL_NEAREST...
 		///@param GLint twrap      --GL_REPEAR, GL_CLAMP_TO_EDGE...
 		///@param GLenum dataType  --GL_UNSIGNED_BYTE, GL_FLOAT16...
@@ -103,7 +99,7 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			type	   = dataType;
 			format	   = dataFormat;
 		}
-		///Method:setDims
+
 		///@param x size of the voxel in x direction
 		///@param y size of the voxel in y direction
 		///@param z size of the voxel in z direction
@@ -116,9 +112,7 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			S[3] = z;
 		}
  
-		///Method:getDims
-		///	get the dimenstions of the voxels.
-
+		///Returns a stim::vec that contains the x, y, z sizes of the voxel.
 		vec<float>
 		getDims()
 		{
@@ -126,11 +120,8 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			return dims;
 		}	
  
-		///Method:setPath
 		///@param file_Path location of the directory with the files
 		///	Sets the path and calls the loader on that path.
-
-
 		void
 		setPath(std::string file_path)
 		{
@@ -139,19 +130,14 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			setTextureType();
 		}
  
-		///Method: getPath
-		///Outputs: string path
-		///	Returns the path associated with an instance of the gl_texture class.
-
+		///	Returns an std::string path associated with an instance of the gl_texture class.
 		std::string
 		getPath()
 		{
 			return path;
 		}
  
-		///Method: getTexture
-		///Outputs: GLuint texID
-		///	Returns the id of the texture create by/associated with the 
+		///	Returns the GLuint id of the texture created by/associated with the 
 		///	instance of the gl_texture class.
  
 		GLuint
@@ -160,7 +146,6 @@ class gl_texture : public virtual image_stack&lt;T&gt;
 			return texID;
 		}
  
-		///Method: createTexture
 		///	Creates a texture and from the loaded data and
 		///	assigns that texture to texID
 		//TO DO :::: 1D textures
-#ifndef RTS_PLANE_H
-#define RTS_PLANE_H
+#ifndef STIM_PLANE_H
+#define STIM_PLANE_H
  
 #include <iostream>
 #include <stim/math/vector.h>
-#include "rts/cuda/callable.h"
+#include <stim/cuda/cudatools/callable.h>
+#include <stim/math/quaternion.h>
  
  
-namespace stim{
-template <typename T, int D> class plane;
+namespace stim
+{
+template<typename T> class plane;
 }
  
-template <typename T, int D>
-CUDA_CALLABLE stim::plane<T, D> operator-(stim::plane<T, D> v);
-
-namespace stim{
-
-template <class T, int D = 3>
-class plane{
-
-	//a plane is defined by a point and a normal
-
-private:
-
-	vec<T, D> P;	//point on the plane
-	vec<T, D> N;	//plane normal
-
-	CUDA_CALLABLE void init(){
-		P = vec<T, D>(0, 0, 0);
-		N = vec<T, D>(0, 0, 1);
-	}
-
-
-public:
-
-	//default constructor
-	CUDA_CALLABLE plane(){
-		init();
-	}
+template<typename T>
+CUDA_CALLABLE stim::plane<T> operator-(stim::plane<T> v);
+
+namespace stim
+{
+
+template <typename T>
+class plane
+{
+	protected:
+		stim::vec<T> P;
+		stim::vec<T> N;
+		stim::vec<T> U;
+
+		///Initializes the plane with standard coordinates.
+		///
+		CUDA_CALLABLE void init()
+		{
+			P = stim::vec<T>(0, 0, 0);
+			N = stim::vec<T>(0, 0, 1);
+			U = stim::vec<T>(1, 0, 0);
+		}
+
+	public:
  
-	CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
-		P = p;
-		N = n.norm();
-	}
-
-	CUDA_CALLABLE plane(T z_pos){
-		init();
-		P[2] = z_pos;
-	}
-
-	//create a plane from three points (a triangle)
-	CUDA_CALLABLE plane(vec<T, D> a, vec<T, D> b, vec<T, D> c){
-		P = c;
-		N = (c - a).cross(b - a);
-		if(N.len() == 0)	//handle the degenerate case when two vectors are the same, N = 0
-			N = 0;
-		else
-			N = N.norm();
-	}
-
-	template< typename U >
-	CUDA_CALLABLE operator plane<U, D>(){
-
-		plane<U, D> result(N, P);
-		return result;
-	}
-
-	CUDA_CALLABLE vec<T, D> norm(){
-		return N;
-	}
-
-	CUDA_CALLABLE vec<T, D> p(){
-		return P;
-	}
-
-	//flip the plane front-to-back
-	CUDA_CALLABLE plane<T, D> flip(){
-		plane<T, D> result = *this;
-		result.N = -result.N;
-		return result;
-	}
-
-	//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
-	CUDA_CALLABLE int face(vec<T, D> v){
-		
-		T dprod = v.dot(N);		//get the dot product between v and N
-
-		//conditional returns the appropriate value
-		if(dprod < 0)
-			return 1;
-		else if(dprod > 0)
-			return -1;
-		else
-			return 0;
-	}
-
-	//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back)
-	CUDA_CALLABLE int side(vec<T, D> p){
-
-		vec<T, D> v = p - P;	//get the vector from P to the query point p
-
-		return face(v);
-	}
-
-	//compute the component of v that is perpendicular to the plane
-	CUDA_CALLABLE vec<T, D> perpendicular(vec<T, D> v){
-		return N * v.dot(N);
-	}
-
-	//compute the projection of v in the plane
-	CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
-		return v - perpendicular(v);
-	}
-
-	CUDA_CALLABLE void decompose(vec<T, D> v, vec<T, D>& para, vec<T, D>& perp){
-		perp = N * v.dot(N);
-		para = v - perp;
-	}
-
-	//get both the parallel and perpendicular components of a vector v w.r.t. the plane
-	CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){
-
-		v_perp = v.dot(N);
-		v_par = v - v_perp;
-	}
-
-	//compute the reflection of v off of the plane
-	CUDA_CALLABLE vec<T, D> reflect(vec<T, D> v){
-
-		//compute the reflection using N_prime as the plane normal
-		vec<T, D> par = parallel(v);
-		vec<T, D> r = (-v) + par * 2;
-
-		/*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
-		std::cout<<str()<<std::endl;
-		std::cout<<"v: "<<v<<std::endl;
-		std::cout<<"r: "<<r<<std::endl;
-		std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
-		std::cout<<"Parallel: "<<par<<std::endl;*/
-		return r;
-
-	}
-
-	CUDA_CALLABLE rts::plane<T, D> operator-()
-	{
-		rts::plane<T, D> p = *this;
-
-		//negate the normal vector
-		p.N = -p.N;
-
-		return p;
-	}
-
-	//output a string
-	std::string str(){
-		std::stringstream ss;
-		ss<<"P: "<<P<<std::endl;
-		ss<<"N: "<<N;
-		return ss.str();
-	}
-
-	///////Friendship
-	//friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
-
-
+		CUDA_CALLABLE plane()
+		{
+			init();
+		}
+
+		CUDA_CALLABLE plane(vec<T> n, vec<T> p = vec<T>(0, 0, 0))
+		{
+			init();
+			P = p;
+			rotate(n.norm());
+		}
+
+		CUDA_CALLABLE plane(T z_pos)
+		{
+			init();
+			P[2] = z_pos;
+		}
+
+		//create a plane from three points (a triangle)
+		CUDA_CALLABLE plane(vec<T> a, vec<T> b, vec<T> c)
+		{
+			init();
+			P = c;
+			stim::vec<T> n = (c - a).cross(b - a);
+			try
+			{
+				if(n.len() != 0)
+				{
+					rotate(n.norm());
+				} else {
+				 	throw 42;
+				}
+			}
+			catch(int i)
+			{
+				std::cerr << "No plane can be creates as all points a,b,c lie on a straight line" << std::endl;
+			}  
+		}
+	
+		template< typename U >
+		CUDA_CALLABLE operator plane<U>()
+		{
+			plane<U> result(N, P);
+			return result;
+
+		}
+
+		CUDA_CALLABLE vec<T> n()
+		{
+			return N;
+		}
+
+		CUDA_CALLABLE vec<T> p()
+		{
+			return P;
+		}
+
+		CUDA_CALLABLE vec<T> u()
+		{
+			return U;
+		}
+
+		///flip the plane front-to-back
+		CUDA_CALLABLE plane<T> flip(){
+			plane<T> result = *this;
+			result.N = -result.N;
+			return result;
+		}
+
+		//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane,     -1 = intersects back)
+		CUDA_CALLABLE int face(vec<T> v){
+			
+			T dprod = v.dot(N);             //get the dot product between v and N
+
+			//conditional returns the appropriate value
+			if(dprod < 0)
+				return 1;
+			else if(dprod > 0)
+				return -1;
+			else
+				return 0;
+		}
+
+		//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = bac    k)
+		CUDA_CALLABLE int side(vec<T> p){
+
+			vec<T> v = p - P;    //get the vector from P to the query point p
+
+			return face(v);
+		}
+
+		//compute the component of v that is perpendicular to the plane
+		CUDA_CALLABLE vec<T> perpendicular(vec<T> v){
+			return N * v.dot(N);
+		}
+
+		//compute the projection of v in the plane
+		CUDA_CALLABLE vec<T> parallel(vec<T> v){
+			return v - perpendicular(v);
+		}
+
+		CUDA_CALLABLE void setU(vec<T> v)
+		{
+			U = (parallel(v.norm())).norm();		
+		}
+
+		CUDA_CALLABLE void decompose(vec<T> v, vec<T>& para, vec<T>& perp){
+			perp = N * v.dot(N);
+			para = v - perp;
+		}
+
+		//get both the parallel and perpendicular components of a vector v w.r.t. the plane
+		CUDA_CALLABLE void project(vec<T> v, vec<T> &v_par, vec<T> &v_perp){
+
+			v_perp = v.dot(N);
+			v_par = v - v_perp;
+		}
+
+		//compute the reflection of v off of the plane
+		CUDA_CALLABLE vec<T> reflect(vec<T> v){
+
+			//compute the reflection using N_prime as the plane normal
+			vec<T> par = parallel(v);
+			vec<T> r = (-v) + par * 2;
+			return r;
+
+		}
+
+		CUDA_CALLABLE stim::plane<T> operator-()
+		{
+			stim::plane<T> p = *this;
+
+			//negate the normal vector
+			p.N = -p.N;
+			return p;
+		}
+
+		//output a string
+		std::string str(){
+			std::stringstream ss;
+			ss<<"P: "<<P<<std::endl;
+			ss<<"N: "<<N<<std::endl;
+			ss<<"U: "<<U;
+			return ss.str();
+		}
+
+
+		CUDA_CALLABLE void rotate(vec<T> n)
+		{
+			quaternion<T> q;
+			q.CreateRotation(N, n);
+			
+			N = q.toMatrix3() * N;
+			U = q.toMatrix3() * U;
+
+		}
+
+		CUDA_CALLABLE void rotate(vec<T> n, vec<T> &X, vec<T> &Y)
+		{
+			quaternion<T> q;
+			q.CreateRotation(N, n);
+			
+			N = q.toMatrix3() * N;
+			U = q.toMatrix3() * U;
+			X = q.toMatrix3() * X;
+			Y = q.toMatrix3() * Y;
+
+		}
  
 };
-
+		
+		
 }
-
-//arithmetic operators
-
-//negative operator flips the plane (front to back)
-//template <typename T, int D>
-
-
-
-
 #endif
+#ifndef RTS_PLANE_H
+#define RTS_PLANE_H
+
+#include <iostream>
+#include <stim/math/vector.h>
+#include "rts/cuda/callable.h"
+
+
+namespace stim{
+template <typename T, int D> class plane;
+}
+
+template <typename T, int D>
+CUDA_CALLABLE stim::plane<T, D> operator-(stim::plane<T, D> v);
+
+namespace stim{
+
+template <class T, int D = 3>
+class plane{
+
+	//a plane is defined by a point and a normal
+
+private:
+
+	vec<T, D> P;	//point on the plane
+	vec<T, D> N;	//plane normal
+
+	CUDA_CALLABLE void init(){
+		P = vec<T, D>(0, 0, 0);
+		N = vec<T, D>(0, 0, 1);
+	}
+
+
+public:
+
+	//default constructor
+	CUDA_CALLABLE plane(){
+		init();
+	}
+	
+	CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
+		P = p;
+		N = n.norm();
+	}
+
+	CUDA_CALLABLE plane(T z_pos){
+		init();
+		P[2] = z_pos;
+	}
+
+	//create a plane from three points (a triangle)
+	CUDA_CALLABLE plane(vec<T, D> a, vec<T, D> b, vec<T, D> c){
+		P = c;
+		N = (c - a).cross(b - a);
+		if(N.len() == 0)	//handle the degenerate case when two vectors are the same, N = 0
+			N = 0;
+		else
+			N = N.norm();
+	}
+
+	template< typename U >
+	CUDA_CALLABLE operator plane<U, D>(){
+
+		plane<U, D> result(N, P);
+		return result;
+	}
+
+	CUDA_CALLABLE vec<T, D> norm(){
+		return N;
+	}
+
+	CUDA_CALLABLE vec<T, D> p(){
+		return P;
+	}
+
+	//flip the plane front-to-back
+	CUDA_CALLABLE plane<T, D> flip(){
+		plane<T, D> result = *this;
+		result.N = -result.N;
+		return result;
+	}
+
+	//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
+	CUDA_CALLABLE int face(vec<T, D> v){
+		
+		T dprod = v.dot(N);		//get the dot product between v and N
+
+		//conditional returns the appropriate value
+		if(dprod < 0)
+			return 1;
+		else if(dprod > 0)
+			return -1;
+		else
+			return 0;
+	}
+
+	//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back)
+	CUDA_CALLABLE int side(vec<T, D> p){
+
+		vec<T, D> v = p - P;	//get the vector from P to the query point p
+
+		return face(v);
+	}
+
+	//compute the component of v that is perpendicular to the plane
+	CUDA_CALLABLE vec<T, D> perpendicular(vec<T, D> v){
+		return N * v.dot(N);
+	}
+
+	//compute the projection of v in the plane
+	CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
+		return v - perpendicular(v);
+	}
+
+	CUDA_CALLABLE void decompose(vec<T, D> v, vec<T, D>& para, vec<T, D>& perp){
+		perp = N * v.dot(N);
+		para = v - perp;
+	}
+
+	//get both the parallel and perpendicular components of a vector v w.r.t. the plane
+	CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){
+
+		v_perp = v.dot(N);
+		v_par = v - v_perp;
+	}
+
+	//compute the reflection of v off of the plane
+	CUDA_CALLABLE vec<T, D> reflect(vec<T, D> v){
+
+		//compute the reflection using N_prime as the plane normal
+		vec<T, D> par = parallel(v);
+		vec<T, D> r = (-v) + par * 2;
+
+		/*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
+		std::cout<<str()<<std::endl;
+		std::cout<<"v: "<<v<<std::endl;
+		std::cout<<"r: "<<r<<std::endl;
+		std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
+		std::cout<<"Parallel: "<<par<<std::endl;*/
+		return r;
+
+	}
+
+	CUDA_CALLABLE rts::plane<T, D> operator-()
+	{
+		rts::plane<T, D> p = *this;
+
+		//negate the normal vector
+		p.N = -p.N;
+
+		return p;
+	}
+
+	//output a string
+	std::string str(){
+		std::stringstream ss;
+		ss<<"P: "<<P<<std::endl;
+		ss<<"N: "<<N;
+		return ss.str();
+	}
+
+	///////Friendship
+	//friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);
+
+
+
+};
+
+}
+
+//arithmetic operators
+
+//negative operator flips the plane (front to back)
+//template <typename T, int D>
+
+
+
+
+#endif
@@ -41,13 +41,14 @@ public:
 		z = u_hat[2]*(T)sin(theta/2);
 	}
  
-	CUDA_CALLABLE void CreateRotation(vec<T> from, vec<T> to){
+	void CreateRotation(vec<T> from, vec<T> to){
  
 		vec<T> r = from.cross(to);			//compute the rotation vector
 		T theta = asin(r.len());				//compute the angle of the rotation about r
 		//deal with a zero vector (both k and kn point in the same direction)
-		if(theta == (T)0)
+		if(theta == (T)0){
 			return;
+		}
  
 		//create a quaternion to capture the rotation
 		CreateRotation(theta, r.norm());
-#ifndef RTS_RECT_H
-#define RTS_RECT_H
+#ifndef STIM_RECT_H
+#define STIM_RECT_H
+
  
 //enable CUDA_CALLABLE macro
 #include <stim/cuda/cudatools/callable.h>
+#include <stim/math/plane.h>
 #include <stim/math/vector.h>
 #include <stim/math/triangle.h>
-#include <stim/math/quaternion.h>
 #include <iostream>
 #include <iomanip>
 #include <algorithm>
+#include <assert.h>
  
 namespace stim{
  
 //template for a rectangle class in ND space
-template <class T>
-struct rect
+template <typename T>
+class rect : plane <T>
 {
 	/*
 		^                   O
 		|                   
 		|                   
-		Y         C         
+		Y         P         
 		|                   
 		|                   
 		O---------X--------->
@@ -28,106 +30,143 @@ struct rect
  
 private:
  
-	stim::vec<T> C;
 	stim::vec<T> X;
 	stim::vec<T> Y;
  
-	CUDA_CALLABLE void scale(T factor){
-		X *= factor;
-		Y *= factor;
-	}
  
  
-	CUDA_CALLABLE void normal(vec<T> n){		//orient the rectangle along the specified normal
-
-		n = n.norm();								//normalize, just in case
-		vec<T> n_current = X.cross(Y).norm();	//compute the current normal
-		quaternion<T> q;							//create a quaternion
-		q.CreateRotation(n_current, n);				//initialize a rotation from n_current to n
-
-		//apply the quaternion to the vectors and position
-		X = q.toMatrix3() * X;
-		Y = q.toMatrix3() * Y;
-	}
-
-	CUDA_CALLABLE void init(){
-		C = vec<T>(0, 0, 0);
-		X = vec<T>(1, 0, 0);
-		Y = vec<T>(0, 1, 0);
-	}
  
 public:
  
-	CUDA_CALLABLE rect(){
+	using stim::plane<T>::n;
+	using stim::plane<T>::P;
+	using stim::plane<T>::N;
+	using stim::plane<T>::U;
+	using stim::plane<T>::rotate;
+
+	///base constructor.
+	CUDA_CALLABLE rect()
+	 : plane<T>()
+	{
 		init();
 	}
  
-	//create a rectangle given a size and position
-	CUDA_CALLABLE rect(T size, T z_pos = (T)0){
+	///create a rectangle given a size and position in Z space.
+	///@param size: size of the rectangle in ND space.
+	///@param z_pos z coordinate of the rectangle.
+	CUDA_CALLABLE rect(T size, T z_pos = (T)0)
+	 : plane<T>(z_pos)
+	{
 		init();			//use the default setup
 		scale(size);	//scale the rectangle
-		C[2] = z_pos;
 	}
  
  
-	//create a rectangle from a center point, normal, and size
-	CUDA_CALLABLE rect(vec<T> c, vec<T> n = vec<T>(0, 0, 1)){
+	///create a rectangle from a center point, normal
+	///@param c: x,y,z location of the center.
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE rect(vec<T> c, vec<T> n = vec<T>(0, 0, 1))
+		: plane<T>()
+	{
 		init();			//start with the default setting
-		C = c;
 		normal(n);		//orient
 	}
  
+	///create a rectangle from a center point, normal, and size
+	///@param c: x,y,z location of the center.
+	///@param s: size of the rectangle.
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE rect(vec<T> c, T s, vec<T> n = vec<T>(0, 0, 1))
+		: plane<T>()
+	{
+		init();			//start with the default setting
+		scale(s);
+		center(c);
+		rotate(n, X, Y);
+	}
+
+	///creates a rectangle from a centerpoint and an X and Y direction vectors.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
 	CUDA_CALLABLE rect(vec<T> center, vec<T> directionX, vec<T> directionY )
+		 : plane<T>((directionX.cross(directionY)).norm(),center)
 	{
-		C = center;
 		X = directionX;
 		Y = directionY;
 	}
  
+	///creates a rectangle from a size, centerpoint, X, and Y direction vectors.
+	///@param size of the rectangle in ND space.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
 	CUDA_CALLABLE rect(T size, vec<T> center, vec<T> directionX, vec<T> directionY )
+		: plane<T>((directionX.cross(directionY)).norm(),center)
 	{	
-		C = center;
 		X = directionX;
 		Y = directionY;
 		scale(size);
 	}
-
-	CUDA_CALLABLE rect(vec<T> size, vec<T> center, vec<T> directionX, vec<T> directionY )
+	
+	///creates a rectangle from a size, centerpoint, X, and Y direction vectors.
+	///@param size of the rectangle in ND space, size[0] = size in X, size[1] = size in Y.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
+	CUDA_CALLABLE rect(vec<T> size, vec<T> center, vec<T> directionX, vec<T> directionY)
+		: plane<T>((directionX.cross(directionY)).norm(), center)
 	{	
-		C = center;
 		X = directionX;
 		Y = directionY;
 		scale(size[0], size[1]);
 	}
-	
-	CUDA_CALLABLE void scale(T factor1, T factor2){
+
+	CUDA_CALLABLE void scale(T factor){
+		X *= factor;
+		Y *= factor;
+	}
+
+	///scales a rectangle in ND space.
+	///@param factor1: size of the scale in the X-direction.
+	///@param factor2: size of the scale in the Y-direction.	
+	CUDA_CALLABLE void scale(T factor1, T factor2)
+	{
 		X *= factor1;
 		Y *= factor2;
 	}
  
+	///@param n; vector with the normal.
+	///Orients the rectangle along the normal n.
+	CUDA_CALLABLE void normal(vec<T> n)
+	{	
+		//orient the rectangle along the specified normal
+		rotate(n, X, Y);
+	}
+
+	///general init method that sets a general rectangle.
+	CUDA_CALLABLE void init()
+	{
+		X = vec<T>(1, 0, 0);
+		Y = vec<T>(0, 1, 0);
+	}
+
 	//boolean comparison
 	bool operator==(const rect<T> & rhs)
 	{
-		if(C == rhs.C && X == rhs.X && Y == rhs.Y)
+		if(P == rhs.P && X == rhs.X && Y == rhs.Y)
 			return true;
 		else
 			return false;
 	}
  
-	/*******************************************
-	Return the normal for the rect
-	*******************************************/
-	CUDA_CALLABLE stim::vec<T> n()
-	{
-        return (X.cross(Y)).norm();
-	}
  
 	//get the world space value given the planar coordinates a, b in [0, 1]
 	CUDA_CALLABLE stim::vec<T> p(T a, T b)
 	{
 		stim::vec<T> result;
 		//given the two parameters a, b = [0 1], returns the position in world space
-		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = this->P - X * (T)0.5 - Y * (T)0.5;
 		result = A + X * a + Y * b;
  
 		return result;
@@ -142,16 +181,16 @@ public:
 	std::string str()
 	{
 		std::stringstream ss;
-		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = P - X * (T)0.5 - Y * (T)0.5;
 		ss<<std::left<<"B="<<std::setfill('-')<<std::setw(20)<<A + Y<<">"<<"C="<<A + Y + X<<std::endl;
 		ss<<std::setfill(' ')<<std::setw(23)<<"|"<<"|"<<std::endl<<std::setw(23)<<"|"<<"|"<<std::endl;
 		ss<<std::left<<"A="<<std::setfill('-')<<std::setw(20)<<A<<">"<<"D="<<A + X;
  
-        return ss.str();
+        	return ss.str();
  
 	}
  
-	//scales the rectangle by a value rhs
+	///multiplication operator scales the rectangle by a value rhs.
 	CUDA_CALLABLE rect<T> operator*(T rhs)
 	{
 		//scales the plane by a scalar value
@@ -164,36 +203,44 @@ public:
  
 	}
  
-	//computes the distance between the specified point and this rectangle
+	///computes the distance between the specified point and this rectangle.
+	///@param p: x, y, z coordinates of the point to calculate distance to.
 	CUDA_CALLABLE T dist(vec<T> p)
 	{
         //compute the distance between a point and this rect
  
-		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = P - X * (T)0.5 - Y * (T)0.5;
  
-        //first break the rect up into two triangles
-        triangle<T> T0(A, A+X, A+Y);
-        triangle<T> T1(A+X+Y, A+X, A+Y);
+		//first break the rect up into two triangles
+		triangle<T> T0(A, A+X, A+Y);
+		triangle<T> T1(A+X+Y, A+X, A+Y);
  
  
-        T d0 = T0.dist(p);
-        T d1 = T1.dist(p);
+		T d0 = T0.dist(p);
+		T d1 = T1.dist(p);
  
-        if(d0 < d1)
-            return d0;
-        else
-            return d1;
+		if(d0 < d1)
+		    return d0;
+		else
+		    return d1;
+	}
+
+	CUDA_CALLABLE T center(vec<T> p)
+	{
+		this->P = p;
 	}
  
+	///Returns the maximum distance of the rectangle from a point p to the sides of the rectangle.
+	///@param p: x, y, z point.
 	CUDA_CALLABLE T dist_max(vec<T> p)
 	{
-		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
-        T da = (A - p).len();
-        T db = (A+X - p).len();
-        T dc = (A+Y - p).len();
-        T dd = (A+X+Y - p).len();
+		vec<T> A = P - X * (T)0.5 - Y * (T)0.5;
+		T da = (A - p).len();
+		T db = (A+X - p).len();
+		T dc = (A+Y - p).len();
+		T dd = (A+X+Y - p).len();
  
-        return std::max( da, std::max(db, std::max(dc, dd) ) );
+		return std::max( da, std::max(db, std::max(dc, dd) ) );
 	}
 };
  
+#ifndef RTS_RECT_H
+#define RTS_RECT_H
+
+//enable CUDA_CALLABLE macro
+#include <stim/cuda/cudatools/callable.h>
+#include <stim/math/vector.h>
+#include <stim/math/triangle.h>
+#include <stim/math/quaternion.h>
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+
+namespace stim{
+
+//template for a rectangle class in ND space
+template <class T>
+struct rect
+{
+	/*
+		^                   O
+		|                   
+		|                   
+		Y         C         
+		|                   
+		|                   
+		O---------X--------->
+	*/
+
+private:
+
+	stim::vec<T> C;
+	stim::vec<T> X;
+	stim::vec<T> Y;
+
+	CUDA_CALLABLE void scale(T factor){
+		X *= factor;
+		Y *= factor;
+	}
+	
+
+
+public:
+
+	///base constructor.
+	CUDA_CALLABLE rect(){
+		init();
+	}
+
+	///create a rectangle given a size and position in Z space.
+	///@param size: size of the rectangle in ND space.
+	///@param z_pos z coordinate of the rectangle.
+	CUDA_CALLABLE rect(T size, T z_pos = (T)0){
+		init();			//use the default setup
+		scale(size);	//scale the rectangle
+		C[2] = z_pos;
+	}
+
+	
+	///create a rectangle from a center point, normal
+	///@param c: x,y,z location of the center.
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE rect(vec<T> c, vec<T> n = vec<T>(0, 0, 1)){
+		init();			//start with the default setting
+		C = c;
+		normal(n);		//orient
+	}
+
+	///create a rectangle from a center point, normal, and size
+	///@param c: x,y,z location of the center.
+	///@param s: size of the rectangle.
+	///@param n: x,y,z direction of the normal.
+	CUDA_CALLABLE rect(vec<T> c, T s, vec<T> n = vec<T>(0, 0, 1)){
+		init();			//start with the default setting
+		C = c;
+		scale(s);
+		normal(n);		//orient
+	}
+
+	///creates a rectangle from a centerpoint and an X and Y direction vectors.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
+	CUDA_CALLABLE rect(vec<T> center, vec<T> directionX, vec<T> directionY )
+	{
+		C = center;
+		X = directionX;
+		Y = directionY;
+	}
+
+	///creates a rectangle from a size, centerpoint, X, and Y direction vectors.
+	///@param size of the rectangle in ND space.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
+	CUDA_CALLABLE rect(T size, vec<T> center, vec<T> directionX, vec<T> directionY )
+	{	
+		C = center;
+		X = directionX;
+		Y = directionY;
+		scale(size);
+	}
+	
+	///creates a rectangle from a size, centerpoint, X, and Y direction vectors.
+	///@param size of the rectangle in ND space, size[0] = size in X, size[1] = size in Y.
+	///@param center: x,y,z location of the center.
+	///@param directionX: u,v,w direction of the X vector.
+	///@param directionY: u,v,w direction of the Y vector.
+	CUDA_CALLABLE rect(vec<T> size, vec<T> center, vec<T> directionX, vec<T> directionY )
+	{	
+		C = center;
+		X = directionX;
+		Y = directionY;
+		scale(size[0], size[1]);
+	}
+
+	///scales a rectangle in ND space.
+	///@param factor1: size of the scale in the X-direction.
+	///@param factor2: size of the scale in the Y-direction.	
+	CUDA_CALLABLE void scale(T factor1, T factor2){
+		X *= factor1;
+		Y *= factor2;
+	}
+
+	///@param n; vector with the normal.
+	///Orients the rectangle along the normal n.
+	CUDA_CALLABLE void normal(vec<T> n){		//orient the rectangle along the specified normal
+
+		n = n.norm();								//normalize, just in case
+		vec<T> n_current = X.cross(Y).norm();	//compute the current normal
+		quaternion<T> q;							//create a quaternion
+		q.CreateRotation(n_current, n);				//initialize a rotation from n_current to n
+
+		//apply the quaternion to the vectors and position
+		X = q.toMatrix3() * X;
+		Y = q.toMatrix3() * Y;
+	}
+
+	///general init method that sets a general rectangle.
+	CUDA_CALLABLE void init(){
+		C = vec<T>(0, 0, 0);
+		X = vec<T>(1, 0, 0);
+		Y = vec<T>(0, 1, 0);
+	}
+
+	//boolean comparison
+	bool operator==(const rect<T> & rhs)
+	{
+		if(C == rhs.C && X == rhs.X && Y == rhs.Y)
+			return true;
+		else
+			return false;
+	}
+
+	/*******************************************
+	Return the normal for the rect
+	*******************************************/
+	CUDA_CALLABLE stim::vec<T> n()
+	{
+        return (X.cross(Y)).norm();
+	}
+
+	//get the world space value given the planar coordinates a, b in [0, 1]
+	CUDA_CALLABLE stim::vec<T> p(T a, T b)
+	{
+		stim::vec<T> result;
+		//given the two parameters a, b = [0 1], returns the position in world space
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+		result = A + X * a + Y * b;
+
+		return result;
+	}
+
+	//parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
+	CUDA_CALLABLE stim::vec<T> operator()(T a, T b)
+	{
+		return p(a, b);
+	}
+
+	std::string str()
+	{
+		std::stringstream ss;
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+		ss<<std::left<<"B="<<std::setfill('-')<<std::setw(20)<<A + Y<<">"<<"C="<<A + Y + X<<std::endl;
+		ss<<std::setfill(' ')<<std::setw(23)<<"|"<<"|"<<std::endl<<std::setw(23)<<"|"<<"|"<<std::endl;
+		ss<<std::left<<"A="<<std::setfill('-')<<std::setw(20)<<A<<">"<<"D="<<A + X;
+
+        return ss.str();
+
+	}
+
+	///multiplication operator scales the rectangle by a value rhs.
+	CUDA_CALLABLE rect<T> operator*(T rhs)
+	{
+		//scales the plane by a scalar value
+
+		//create the new rectangle
+		rect<T> result = *this;
+		result.scale(rhs);
+
+		return result;
+
+	}
+
+	///computes the distance between the specified point and this rectangle.
+	///@param p: x, y, z coordinates of the point to calculate distance to.
+	CUDA_CALLABLE T dist(vec<T> p)
+	{
+        //compute the distance between a point and this rect
+
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+
+        //first break the rect up into two triangles
+        triangle<T> T0(A, A+X, A+Y);
+        triangle<T> T1(A+X+Y, A+X, A+Y);
+
+
+        T d0 = T0.dist(p);
+        T d1 = T1.dist(p);
+
+        if(d0 < d1)
+            return d0;
+        else
+            return d1;
+	}
+
+	CUDA_CALLABLE T center(vec<T> p)
+	{
+		C = p;
+	}
+
+	///Returns the maximum distance of the rectangle from a point p to the sides of the rectangle.
+	///@param p: x, y, z point.
+	CUDA_CALLABLE T dist_max(vec<T> p)
+	{
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
+        T da = (A - p).len();
+        T db = (A+X - p).len();
+        T dc = (A+Y - p).len();
+        T dd = (A+X+Y - p).len();
+
+        return std::max( da, std::max(db, std::max(dc, dd) ) );
+	}
+};
+
+}	//end namespace rts
+
+template <typename T, int N>
+std::ostream& operator<<(std::ostream& os, stim::rect<T> R)
+{
+    os<<R.str();
+    return os;
+}
+
+
+#endif
@@ -71,8 +71,9 @@ struct vec : public std::vector&lt;T&gt;
 	vec( const vec<T>& other){
 		unsigned int N = other.size();
 		resize(N);							//resize the current vector to match the copy
-		for(unsigned int i=0; i<N; i++)		//copy each element
-			at(i) = other[i];
+		for(unsigned int i=0; i<N; i++){	//copy each element
+			at(i) = other[i];
+		}
 	}
  
 	//I'm not sure what these were doing here.
+#ifndef STIM_CYLINDER_H
+#define STIM_CYLINDER_H
+#include <iostream>
+#include <stim/math/circle.h>
+#include <stim/math/vector.h>
+
+
+namespace stim
+{
+template<typename T>
+class cylinder
+{
+	private:
+		stim::circle<T> s;
+		std::vector< stim::vec<T> > pos;
+		std::vector< stim::vec<T> > mags;
+		std::vector< T > L;
+	
+		void
+		init()
+		{
+
+		}
+
+		void
+		init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
+		{
+			pos = inP;
+			mags = inM;
+			L.resize(pos.size()-1);
+			T temp = (T)0;
+			for(int i = 0; i < L.size()-1; i++)
+			{
+				temp += (pos[i] - pos[i+1]).len();
+				L[i] = temp;
+			}
+		}
+		
+		stim::vec<T>
+		d(int idx)
+		{
+			return (pos[idx] - pos[idx+1]).norm();
+			
+		}
+
+		T
+		getl(int j)
+		{
+			for(int i = 0; i < j-1; ++i)
+			{
+				temp += (pos[i] - pos[i+1]).len();
+				L[i] = temp;
+			}
+		}
+
+		int
+		findIdx(T l)
+		{
+			int i = pos.size()/2;
+			while(i > 0 && i < pos.size())
+			{
+				if(L[i] < l)
+				{
+					i = i/2;
+				}
+				else if(L[i] < l && L[i+1] > l)
+				{
+					break;
+				}
+				else
+				{
+					i = i+i/2;
+				}
+			}
+			return i;
+		}
+
+	public:
+		cylinder()
+		{
+
+		}
+
+		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
+		///The higher the number of sides, the more rectangeles compose the surface of the cylinder.
+		///@param inP:  Vector of stim vecs composing the points of the centerline.
+		///@param inM:  Vector of stim vecs composing the radii of the centerline.
+		cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
+		{
+			init(inP, inM);
+		}
+
+
+		///Returns a position vector at the given p-value (p value ranges from 0 to 1).
+		stim::vec<T>
+		p(T pvalue)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0)
+				return;
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
+		}
+
+		stim::vec<T>
+		p(T l, int idx)
+		{
+			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
+		}
+
+		///Returns a radius at the given p-value (p value ranges from 0 to 1).
+		T
+		r(T pvalue)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0)
+				return;
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
+		}
+
+		T
+		r(T l, int idx)
+		{
+			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
+		}
+
+
+		///returns the position of the point with a given pvalue and theta on the surface
+		///in x, y, z coordinates. Theta is in degrees from 0 to 360
+		stim::vec<T>
+		surf(T pvalue, T theta)
+		{
+			if(pvalue < 0.0 || pvalue > 1.0)
+				return;
+			T l = pvalue*L[L.size()-1];
+			int idx = findIdx(l);
+			stim::vec<T> ps = p(l, idx); 
+			T m = r(l, idx);
+			stim::vec<T> dr = d(idx);
+			s = stim::circle<T>(ps, m, dr);
+			return(s.p(theta));
+		}
+
+		std::vector<std::vector<vec<T> > >
+		getPoints(int sides)
+		{
+			if(pos.size() < 2)
+			{
+				return;
+			} else {
+				std::vector<std::vector <vec<T> > > points;
+				points.resize(pos.size());
+				stim::vec<T> d = (pos[0] - pos[1]).norm();
+				s = stim::circle<T>(pos[0], mags[0][0], d);
+				points[0] = s.getPoints(sides);
+				for(int i = 1; i < pos.size(); i++)
+				{
+					d = (pos[i] - pos[i-1]).norm();
+					s.center(pos[i]);
+					s.normal(d);
+					s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]);
+					points[i] = s.getPoints(sides);
+				}
+				return points;
+			}
+		}
+		
+};
+
+}
+#endif
@@ -29,6 +29,8 @@ private:
 	void
 	init()
 	{
+		if(glIsList(dList))
+			glDeleteLists(dList, 1);
 		dList = glGenLists(1);
 		glListBase(dList);
  
@@ -40,16 +42,25 @@ private:
 	}
  
 	void
-	Create()
+	Create(GLenum mode)
 	{
+//		GLuint selectBuf[2048];
+//		GLint hits;
+//		glSelectBuffer(2048, selectBuf);
+		
 		int len = (int) stim::obj<T>::numL();
 		std::vector< stim::vec<float> > line;
 		glNewList(dList, GL_COMPILE);
 	//	glColor3f(0.0, 1.0, 0.0);
-		glLineWidth(2.5);
+		glLineWidth(3.5);
 		for(int i = 0; i < len; i++){
 			line = stim::obj<T>::getL_V(i);
-			glColor3ub(rand()%255, rand()%255, rand()%255);
+			if(mode == GL_SELECT)
+			{
+				glLoadName(i);
+			}
+			glColor3f(0.0, 1.0, 0.05);
+			//glColor3ub(rand()%255, rand()%255, rand()%255);
 			glBegin(GL_LINE_STRIP);
 			for(int j = 0; j < line.size(); j++){
 					glVertex3f(
@@ -71,21 +82,21 @@ public:
 	}
  
 	void
-	createFromSelf()
+	createFromSelf(GLenum mode = GL_RENDER)
 	{
 	//	glPopMatrix();
 		init();
-		Create();
+		Create(mode);
 	//	glPushMatrix();
 	}
  
 	void
-	createFromFile(std::string filename)
+	createFromFile(std::string filename, GLenum mode = GL_RENDER)
 	{
 		stim::obj<T>::load(filename);
 		glPushMatrix();		//Safety Operation to avoid changing the current matrix.
 		init();
-		Create();
+		Create(mode);
 		glPopMatrix();
 		CHECK_OPENGL_ERROR
 	}
+#ifndef STIM_GLNETWORK_H
+#define STIM_GLNETWORK_H
+
+#include <GL/glew.h>
+#include <GL/glut.h>
+#include "network.h"
+#include <stim/visualization/cylinder.h>
+#include <stim/math/vector.h>
+#include <list>
+#include <ANN/ANN.h>
+#include "fiber.h"
+
+
+namespace stim
+{
+template <typename T>
+class glnetwork : public virtual network<T>
+{
+private:
+	using stim::network<T>::E;
+
+	GLuint dList;		///displaylist for the lines.
+	GLuint cList;		///displaylist for the cylinders.
+
+	void init()
+	{
+		///clear lists if there is data in them.
+		///adding points may create errors or uncessary duplicate points.
+		if(glIsList(dList))
+			glDeleteLists(dList, 1);
+		if(glIsList(cList))
+			glDeleteLists(cList, 1);
+		dList = glGenLists(1);		///create the lists
+		cList = glGenLists(1);
+
+		///set up the Line list.
+		glListBase(dList);
+		glMatrixMode(GL_PROJECTION);
+		glLoadIdentity;
+		glMatrixMode(GL_MODELVIEW);
+		glLoadIdentity;
+
+		///set up the cylinder List.
+		glListBase(cList);
+		glMatrixMode(GL_PROJECTION);
+		glLoadIdentity;
+		glMatrixMode(GL_MODELVIEW);
+		glLoadIdentity;
+	}
+
+	void
+	Create(GLenum mode, int sides = 8)
+	{
+		glListBase(dList);
+		glNewList(dList, GL_COMPILE);
+		glLineWidth(3.5);
+		for(int i = 0; i < E.size(); i++)
+		{
+			if(mode == GL_SELECT)
+			{
+//				glLineWidth(3.5);
+				glLoadName(i);
+			}
+			else{
+//				glLineWidth(1.0+1.0*i);
+			}
+			glColor3f(0.0, 1.0-0.05*i, i*0.05);
+			std::vector<stim::vec<T> > line = getEdgeCenterLine(i);
+			glBegin(GL_LINE_STRIP);
+			for(int j = 0; j < line.size(); j++)
+			{
+				glVertex3f(line[j][0],
+					line[j][1],
+					line[j][2]);
+			}
+			glEnd();
+		}
+		glEndList();
+
+		glListBase(cList);
+		glNewList(cList, GL_COMPILE);
+
+		for(int i = 0; i < E.size(); i++)
+		{
+			if(mode == GL_SELECT)
+			{
+				glLoadName(i);
+			}
+			glColor3f(1.0, 1.0, 0.0);
+			std::vector<stim::vec<T> > line = getEdgeCenterLine(i);
+			std::vector<stim::vec<T> > linemag = getEdgeCenterLineMag(i);
+			stim::cylinder<T> cyl(line, linemag);
+			std::vector<std::vector<stim::vec<T > > > p = cyl.getPoints(sides);
+			for(int i = 0; i < p.size()-1; i++)
+			{
+		                for(int j = 0; j < p[0].size()-1; j++)
+                		{
+		                        glColor4f(1.0, 1.0, 0.0, 0.5);
+        		                glEnable(GL_BLEND);
+                		        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+                		        glBegin(GL_QUADS);
+                		                glVertex3f(p[i][j][0], p[i][j][1], p[i][j][2]);
+                        		        glVertex3f(p[i][j+1][0], p[i][j+1][1], p[i][j+1][2]);
+                       		                glVertex3f(p[i+1][j+1][0], p[i+1][j+1][1], p[i+1][j+1][2]     );
+                       			        glVertex3f(p[i+1][j][0], p[i+1][j][1], p[i+1][j][2]);
+                    			glEnd();
+                        		glDisable(GL_BLEND);
+
+                        		glColor4f(1.0, 0.0, 1.0, 1.0);
+					glLineWidth(2.0);
+                       			glBegin(GL_LINES);
+                       	        		glVertex3f(p[i][j][0], p[i][j][1], p[i][j][2]);
+               			                glVertex3f(p[i][j+1][0], p[i][j+1][1], p[i][j+1][2]);
+                        		        glVertex3f(p[i][j][0], p[i][j][1], p[i][j][2]);
+                           			glVertex3f(p[i+1][j][0], p[i+1][j][1], p[i+1][j][2]     );
+                        		glEnd();
+                		}
+
+			}
+
+			
+
+		}
+		glEndList();
+//		CHECK_OPENGL_ERROR
+	}
+	
+public:
+	using stim::network<T>::sizeE;
+	using stim::network<T>::getEdgeCenterLine;
+	using stim::network<T>::getEdgeCenterLineMag;
+	glnetwork()
+	{
+
+	}
+	
+	void
+	createFromSelf(GLenum mode = GL_RENDER, int sides = 8)
+	{
+		init();
+		Create(mode, sides);
+	}
+
+	void
+	Render()
+	{
+		glCallList(dList);
+//		CHECK_OPENGL_ERROR
+	}
+
+	void
+	RenderCylinders()
+	{
+		glCallList(cList);
+//		CHECK_OPENGL_ERROR
+	}
+
+	void
+	RenderLine(std::vector<stim::vec<T> > l)
+	{
+                glColor3f(0.5, 1.0, 0.5);
+                glLineWidth(3.0);
+                glBegin(GL_LINE_STRIP);
+                for(int j = 0; j < l.size(); j++){
+                        glVertex3f(
+                                l[j][0],
+                                l[j][1],
+                                l[j][2]
+                                );
+                }
+                glEnd();
+	}
+	
+};
+}
+
+#endif
@@ -683,14 +683,14 @@ public:
 		l.resize(nP);
  
 		//copy the points from the point list to the stim vector
-		unsigned int pi;
+		unsigned int pie;
 		for(unsigned int p = 0; p < nP; p++){
  
 			//get the index of the geometry point
-			pi = L[i][p][0] - 1;
+			pie = L[i][p][0] - 1;
  
 			//get the coordinates of the current point
-			stim::vec<T> newP = V[pi];
+			stim::vec<T> newP = V[pie];
  
 			//copy the point into the vector
 			l[p] = newP;
@@ -737,14 +737,14 @@ public:
 		l.resize(nP);
  
 		//copy the points from the point list to the stim vector
-		unsigned int pi;
+		unsigned int pie;
 		for(unsigned int p = 0; p < nP; p++){
  
 			//get the index of the geometry point
-			pi = L[i][p][1] - 1;
+			pie = L[i][p][1] - 1;
  
 			//get the coordinates of the current point
-			stim::vec<T> newP = VT[pi];
+			stim::vec<T> newP = VT[pie];
  
 			//copy the point into the vector
 			l[p] = newP;