#ifndef STIM_GL_SPIDER_H
#define STIM_GL_SPIDER_H

//#include <GL/glew.h>
#include <GL/glut.h>
#include <cuda.h>
#include <cuda_gl_interop.h>
#include <cudaGL.h>
#include <math.h>
#include <stim/gl/gl_texture.h>
#include <stim/visualization/camera.h>
#include <stim/gl/error.h>
#include <stim/math/vector.h>
#include <stim/math/vec3.h>
#include <stim/math/rect.h>
#include <stim/math/matrix.h>
#include <stim/math/constants.h>
#include <stim/cuda/spider_cost.cuh>
#include <stim/cuda/cudatools/glbind.h>
#include <stim/cuda/arraymath.cuh>
#include <stim/cuda/cuda_texture.cuh>
#include <stim/cuda/cudatools.h>
#include <stim/cuda/ivote.cuh>
#include <stim/visualization/glObj.h>
#include <vector>
#include <stim/cuda/branch_detection.cuh>
#include "../../../volume-spider/fiber.h"
#include "../../../volume-spider/glnetwork.h"
#include <stim/visualization/cylinder.h>
#include <stim/cuda/testKernel.cuh>

//#include <stim/cuda/testKernel.cuh>

#include <iostream>
#include <fstream>
#ifdef TIMING
	#include <ctime>
	#include <cstdio>	
#endif

#ifdef TESTING
	#include <iostream>
	#include <cstdio>
	#include <ctime>
#endif


namespace stim
{

template<typename T>
class gl_spider : public virtual gl_texture<T>
{
	//doen't use gl_texture really, just needs the GLuint id.
	//doesn't even need the texture iD really.
	private:

		#ifdef TIMING
			double branch_time;// = 0;
			double direction_time;// = 0;
			double position_time;// = 0;
			double size_time;// = 0;
			double cost_time;// = 0;
			double network_time;// = 0;
			double hit_time;// = 0;
		#endif

		//
		stim::vec3<float> p;  	//vector designating the position of the spider.
		stim::vec3<float> d;	//vector designating the orientation of the spider
						//always a unit vector.
		stim::vec<float> m;	//magnitude of the spider vector.
						//mag[0] = length.
						//mag[1] = width.
		std::vector<stim::vec3<float> > dV;	//A list of all the direction vectors.
		std::vector<stim::vec3<float> > pV;	//A list of all the position vectors.
		std::vector<stim::vec3<float> > mV;	//A list of all the size vectors.

		stim::matrix<float, 4> cT;		//current Transformation matrix
							//From tissue space to texture space.
		GLuint texID;
		stim::vec<float> S;			//Size of a voxel in the volume.
		stim::vec<float> R;			//Dimensions of the volume.


		//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;// = 5.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::vec3<float> > seeds;		//seed positions.
		std::stack< stim::vec3<float> > seedsvecs;	//seed directions.
		std::stack< float > 		seedsmags;	//seed magnitudes.

		std::vector< stim::vec3<float> > cL;	//Positions of line currently being traced.
		std::vector< stim::vec3<float> > cD;	//Direction of line currently being traced.
		std::vector< stim::vec<float> > cM;	//Magnitude of line currently being traced.

		stim::glnetwork<float> nt;		//object for storing the network.

		stim::vec<float> rev;			//reverse vector;
		stim::camera camSel;
		stim::vec3<float> ps;
		stim::vec3<float> ups;
		stim::vec3<float> ds;

		//static const float t_length = 16.0;
		float t_length;

		//cuda texture variables that keep track of the binding.
		stim::cuda::cuda_texture t_dir;
		stim::cuda::cuda_texture t_pos;
		stim::cuda::cuda_texture t_mag;
		

//--------------------------------------------------------------------------//
//-------------------------------PRIVATE METHODS----------------------------//
//--------------------------------------------------------------------------//

		/// Method for finding the best scale for the spider.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. 
		void
		findOptimalDirection()
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif
			setMatrix();			//create the transformation matrix.
			glCallList(dList);		//move the templates to p, d, m.
			glFinish();
//			glFlush();
			#ifdef TIMING
				direction_time += gpuStopTimer();
			#endif
			#ifdef TESTING
//				test(texbufferID, GL_TEXTURE_2D,2*t_length,numSamples*t_length, "Final_Cost_Direction.bmp");
			#endif

			int best = getCost(t_dir.getTexture(), t_dir.getAuxArray() ,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();	//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 (direction) for the spider.
		/// Not sure if necessary since the next p (position) for the spider
		/// will be at d * m.
		void
		findOptimalPosition()
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif
			setMatrix();			//create the transformation matrix.
			glCallList(dList+1);		//move the templates to p, d, m.
			glFinish();
//			glFlush();
			#ifdef TIMING
				position_time += gpuStopTimer();
			#endif

			#ifdef TESTING
//				test(ptexbufferID, GL_TEXTURE_2D,2*t_length, numSamplesPos*t_length, "Final_Cost_Position.bmp");
			#endif
			int best = getCost(t_pos.getTexture(), t_pos.getAuxArray(), numSamplesPos);		//find min cost.
//			std::cerr << best << std::endl;
			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.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. */
		void
		findOptimalScale()
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif
			setMatrix();			//create the transformation.
			glCallList(dList+2);		//move the templates to p, d, m.
			glFinish();
//			glFlush();
			#ifdef TIMING
				size_time += gpuStopTimer();
			#endif
			#ifdef TESTING
//				test(mtexbufferID, GL_TEXTURE_2D, 2*t_length, numSamplesMag*t_length, "Final_Cost_Position.bmp");
			#endif
			int best = getCost(t_mag.getTexture(), t_mag.getAuxArray(), numSamplesMag);		//get best cost.
			setMagnitude(m[0]*mV[best][0]);	//adjust the magnitude.
		}


		///subject to change.
		///finds branches.
		///depreciated
		void
		branchDetection()
		{
			setMatrix();
			glCallList(dList+3);
			std::vector< stim::vec<float> > result = find_branch(
					btexbufferID, GL_TEXTURE_2D, 16, 216);
			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
			if(!result.empty())
			{
				for(int i = 1; i < result.size(); i++)
				{
					stim::vec<float> cylp(
						0.5 * cos(2*M_PI*(result[i][1])),
						0.5 * sin(2*M_PI*(result[i][1])),
						result[i][0]-0.5,
						1.0);
					cylp = cT*cylp;
					
					stim::vec3<float> vec(
						cylp[0]*S[0]*R[0],
						cylp[1]*S[1]*R[1],
						cylp[2]*S[2]*R[2]);
					stim::vec3<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],
						   -p[1] + cylp[1]*S[1]*R[1],
						   -p[2] + cylp[2]*S[2]*R[2]);
						seeddir = seeddir.norm();
					float seedm = m[0];
// Uncomment for global run 
					if(
					 !(vec[0] > size[0] || vec[1] > size[1]
					 || vec[2] > size[2] || vec[0] < 0
					 || vec[1] < 0 || vec[2] < 0))
					{
						setSeed(vec);
						setSeedVec(seeddir);
						setSeedMag(seedm);
					}
				}
			}    
			
		}


		///finds all the branches in the a given fiber.
		///using LoG method.
		void
		branchDetection2(int n = 8, int l_template = 8, int l_square = 8)
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif

			if(cL.size() < 4){}
			else{
			setMatrix(1);
			DrawLongCylinder(n, l_template, l_square);
			stim::cylinder<float> cyl(cL, cM);
			std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
			float pval;
			if(!result.empty())
			{
				for(int i = 0; i < result.size(); i++)
				{
					int id = result[i][2];
					if(fmod(result[i][2], id) != 0 && id != 0)
					{
						
						pval = ((cyl.getl(id+1)-cyl.getl(id))*
							(fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);
					}
					else if(id == 0)
					{
						pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);
					}
					else
					{
						pval = (cyl.getl(id)/cyl.getl(cL.size()-1));
					}
					stim::vec3<float> v = cyl.surf(pval, result[i][0]);
					stim::vec3<float> di = cyl.p(pval);
					float rad = cyl.r(pval);
					if(
					 !(v[0] > size[0] || v[1] > size[1]
					 || v[2] > size[2] || v[0] < 0
					 || v[1] < 0 || v[2] < 0))
					{
						setSeed(v);
						setSeedVec((v-di).norm());
						setSeedMag(rad);
					}
				}
			}
			}
			#ifdef TIMING
				branch_time += gpuStopTimer();
			#endif
		}


		float uniformRandom()
		{
			return (  (float)(rand()))/(  (float)(RAND_MAX));
		}

		float normalRandom()
		{
			float u1 = uniformRandom();
			float u2 = uniformRandom();
			return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));
		}

		stim::vec3<float> uniformRandVector()
		{
			stim::vec3<float> r(uniformRandom(), uniformRandom(), 1.0);
			return r;
		}

		stim::vec3<float> normalRandVector()
		{
			stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);
			return r;
		}

		
//--------------------------------------------------------------------------//
//---------------------TEMPLATE CREATION METHODS----------------------------//
//--------------------------------------------------------------------------//



		///@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 = M_PI/2)
		{

			//Set up the vectors necessary for Rectangle creation.
			stim::vec3<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			stim::vec3<float> pos(0.0,0.0,0.0);
			stim::vec3<float> mag(1.0, 1.0, 1.0);
			stim::vec3<float> dir(0.0, 0.0, 1.0);

			float PHI[2], Z[2], range;
			PHI[0] = solidAngle/2;
			PHI[1] = asin(0);
			
			Z[0] = cos(PHI[0]);
			Z[1] = cos(PHI[1]);
			
			range = Z[0] - Z[1];

			float z, theta, phi;
			glNewList(dList, GL_COMPILE); 		
			for(int i = 0; i < numSamples; i++)
			{
				z = uniformRandom()*range + Z[1];
				theta = uniformRandom()*2*M_PI;
				phi = acos(z);
				stim::vec3<float> sph(1, theta, phi);
				stim::vec3<float> cart = sph.sph2cart();
				dV.push_back(cart);
				if(cos(Y.dot(cart)) < 0.087)
				{
						Y[0] = 0.0; Y[1] = 1.0;
				}else{
					Y[0] = 1.0; Y[1] = 0.0;
				}
                                hor = stim::rect<float>(mag,
					 pos, cart,
                       			((Y.cross(cart)).cross(cart)).norm());
       				ver = stim::rect<float>(mag,
					 pos, cart,
                       			hor.n());
				UpdateBuffer(0.0, 0.0+i*t_length);
			}
			glEndList();
		}

		///@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.
			stim::vec3<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			stim::vec3<float> pos(0.0,0.0,0.0);
			stim::vec3<float> mag(1.0, 1.0, 1.0);
			stim::vec3<float> dir(0.0, 0.0, 1.0);

			//Set up the variable necessary for vector creation.
			glNewList(dList+1, GL_COMPILE);
			for(int i = 0; i < numSamplesPos; i++)
			{
				stim::vec3<float> temp = uniformRandVector();
				temp = temp*delta*2.0 - delta/2.0;
				temp[2] = 0.0;
				pV.push_back(temp);
                		hor = stim::rect<float>(mag,
					 temp, dir,
                      			((Y.cross(d)).cross(d))
					.norm());
                		ver = stim::rect<float>(mag,
					 temp, dir,
                        		hor.n());
				UpdateBuffer(0.0, 0.0+i*t_length);
			}
			glEndList();
		}

		///@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)
		{
			
			//Set up the vectors necessary for Rectangle creation.
			stim::vec3<float> Y(1.0, 0.0, 0.0);		//orthogonal vec.
			stim::vec3<float> pos(0.0, 0.0, 0.0);
			stim::vec3<float> mag(1.0, 1.0, 1.0);
			stim::vec3<float> dir(0.0, 0.0, 1.0);

			//Set up the variable necessary for vector creation.
			int dim = (sqrt(numSamplesMag)-1)/2;
			float min 	= 1.0-delta;
			float max 	= 1.0+delta;
			float step	= (max-min)/(numSamplesMag-1);
			float factor;
			stim::vec3<float> temp(0.0,0.0,0.0);

			glNewList(dList+2, GL_COMPILE);
			for(int i = 0; i < numSamplesMag; i++){
				//Create linear index
				factor = (min+step*i)*mag[0];
				temp = factor;
				mV.push_back(temp);	
				hor = stim::rect<float>(temp,
					 pos, dir, 
       	       				((Y.cross(d)).cross(d))
					.norm());
               			ver = stim::rect<float>(temp,
					 pos, dir,
                       			hor.n());
				UpdateBuffer(0.0, 0.0+i*t_length);
			CHECK_OPENGL_ERROR
			}
			glEndList();
		}

		///@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)
		{	
			stim::vec3<float>p1; 
        	        stim::vec3<float>p2; 
	                stim::vec3<float>p3; 
                	stim::vec3<float>p4;	
			p1 = hor.p(1,1);
			p2 = hor.p(1,0);
			p3 = hor.p(0,0);
			p4 = hor.p(0,1);
			glBegin(GL_QUADS);
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x,v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+t_length, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+t_length, v_y+t_length);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x, v_y+t_length);
			 glEnd();

			 p1 = ver.p(1,1);
			 p2 = ver.p(1,0);
			 p3 = ver.p(0,0);
			 p4 = ver.p(0,1);
		 	 glBegin(GL_QUADS);
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x+t_length, v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+2.0*t_length, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+2.0*t_length, v_y+t_length);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x+t_length, v_y+t_length);
			glEnd(); 
		}
		


//--------------------------------------------------------------------------//
//--------------------------------GL METHODS--------------------------------//
//--------------------------------------------------------------------------//

		///@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)
		{
			glDeleteFramebuffers(1, &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, NULL);   
//			delete[] texels;
			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
			glBindTexture(GL_TEXTURE_2D, 0);
		}

		///@param uint width sets the width of the buffer.
		///@param uint height sets the height of the buffer.
		///Function for setting up the 2D buffer that stores the samples.
		void
		GenerateFBO(unsigned int width, unsigned int height)
		{
			glGenFramebuffers(1, &fboID);
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
//			int numChannels = 1;
//			unsigned char* texels = new unsigned char[width * height * numChannels];
			glGenTextures(1, &texbufferID);
			glBindTexture(GL_TEXTURE_2D, texbufferID);

			//Textures repeat and use linear interpolation, luminance format.
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,
				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, NULL);   
//			delete[] texels;
			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
			glBindTexture(GL_TEXTURE_2D, 0);
			CHECK_OPENGL_ERROR
		}


		///IF type == 0
		///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).

		///IF type == 1
		///Method for using the gl manipulation to set up a matrix 
		///To transform from tissue space into texture space.
		///All transformation happen in glMatrixMode(GL_TEXTURE).
		///All transformation happen in glMatrixMode(GL_TEXTURE).
		void setMatrix(int type = 0)
		{
			if(type == 0)
			{
				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();

				//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();
				
				CHECK_OPENGL_ERROR
				//revert back to default gl mode.
				glMatrixMode(GL_MODELVIEW);
			}
			else if(type == 1)
			{
				glMatrixMode(GL_TEXTURE);
				glLoadIdentity();
				glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
				glMatrixMode(GL_MODELVIEW);
			}
		}

		///Method for controling the buffer and texture binding.
		///Clears the buffer upon binding.
		void
		Bind()
		{
			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*t_length, numSamples*t_length);
			gluOrtho2D(0.0,2.0*t_length,0.0,numSamples*t_length);
			glEnable(GL_TEXTURE_3D);
			glBindTexture(GL_TEXTURE_3D, texID);

			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);
			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
			findOptimalDirection();
			findOptimalPosition();
			findOptimalScale();
			Unbind();
			Bind(btexbufferID, bfboID, 27);
			branchDetection();
			Unbind();
			return current_cost;
		}

		


//--------------------------------------------------------------------------//
//--------------------------------CUDA METHODS------------------------------//
//--------------------------------------------------------------------------//
		

		///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()
		{
			stim::vec<int> cost = 
//				stim::cuda::get_cost(texbufferID, GL_TEXTURE_2D, numSamples);
			cudaDeviceSynchronize();
			current_cost = cost[1];
			return cost[0];
		}

//		int
//		getCost(GLuint tID, int n)
//		{
//			#ifdef TIMING
//				gpuStartTimer();
//			#endif
//			stim::vec<int> cost = 
//				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n, 2*t_length, t_length);
//			#ifdef TIMING
//				cost_time += gpuStopTimer();
//			#endif
//			current_cost = cost[1];
//			return cost[0];
//		}

		int
		getCost(cudaTextureObject_t tObj, float* result, int n)
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif
			stim::vec<int> cost = 
				stim::cuda::get_cost(tObj, result, n, 2*t_length, t_length);
			#ifdef TIMING
				cost_time += gpuStopTimer();
			#endif
			current_cost = cost[1];
			return cost[0];
		}

	public:

		///ininializes the cuda device and environment.
		void
		initCuda()
		{	
			stim::cudaSetDevice();
		}

		//horizonal rectangle forming the spider.
		stim::rect<float> hor;
		//vectical rectangle forming the spider.
		stim::rect<float> ver;	

		//Timing variables.
		#ifdef TESTING
			std::clock_t start;
		#endif

//--------------------------------------------------------------------------//
//-----------------------------CONSTRUCTORS---------------------------------//
//--------------------------------------------------------------------------//


		///@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 samplespos = 441,int samplesmag = 144)
		{
//			std::cout << "I ran this constructor" << std::endl;
			p = stim::vec3<float>(0.0, 0.0, 0.0);
			d = stim::vec3<float>(0.0, 0.0, 1.0);
			m = stim::vec<float>(1.0, 1.0);
			S = stim::vec3<float>(1.0, 1.0, 1.0);
			R = stim::vec3<float>(1.0, 1.0, 1.0);
//			std::cout << samples << std::endl;
			numSamples = samples;
//			std::cout << numSamples << std::endl;
			numSamplesPos = samplespos;
			numSamplesMag = samplesmag;
		}

		///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, int numsamplespos = 441, int numsamplesmag =144)
		{
			p = stim::vec3<float>(pos_x, pos_y, pos_z);
			d = stim::vec3<float>(dir_x, dir_y, dir_z);
			m = stim::vec<float>(mag_x, mag_x, mag_x);
			S = stim::vec3<float>(1.0,1.0,1.0);
			R = stim::vec3<float>(1.0,1.0,1.0);
			numSamples = numsamples;
			numSamplesPos = numsamplespos;
			numSamplesMag = numsamplesmag;
		}

		///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::vec3<float> pos, stim::vec3<float> dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144)
		{
			p = pos;
			d = dir;
			m = vec<float>(mag, mag, mag);
			S = vec3<float>(1.0,1.0,1.0);
			R = vec3<float>(1.0,1.0,1.0);
			numSamples = samples;
			numSamplesPos = samplesPos;
			numSamplesMag = samplesMag;
		}

		///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.
		///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.	
		void
		attachSpider(GLuint id)
		{
#ifdef TIMING
			branch_time = 0;
			direction_time = 0;
			position_time = 0;
			size_time = 0;
			cost_time = 0;
			network_time = 0;
			hit_time = 0;
#endif
			stepsize = 2.5;
			t_length = 16.0;

			srand(100);	
			texID = id;
			//GenerateFBO(16, numSamples*8);
			GenerateFBO(t_length*2, numSamples*t_length, texbufferID, fboID); 
			std::cout << numSamples << std::endl;
			CHECK_OPENGL_ERROR
			GenerateFBO(t_length*2, numSamplesPos*t_length, ptexbufferID, pfboID); 
			std::cout << numSamplesPos << std::endl;
			CHECK_OPENGL_ERROR
			GenerateFBO(t_length*2, numSamplesMag*t_length, mtexbufferID, mfboID); 
			std::cout << numSamplesMag << std::endl;
			CHECK_OPENGL_ERROR
			GenerateFBO(16, 216, btexbufferID, bfboID); 
			CHECK_OPENGL_ERROR
			t_dir.MapCudaTexture(texbufferID, GL_TEXTURE_2D);
			t_dir.Alloc(numSamples);
			t_pos.MapCudaTexture(ptexbufferID, GL_TEXTURE_2D);
			t_pos.Alloc(numSamplesPos);
			t_mag.MapCudaTexture(mtexbufferID, GL_TEXTURE_2D);
			t_mag.Alloc(numSamplesMag);
//			setDims(0.6, 0.6, 1.0);
//			setSize(1024.0, 1024.0, 1024.0);
			setMatrix();
			dList = glGenLists(3);
			glListBase(dList);
			Bind(texbufferID, fboID, numSamples, t_length);
				genDirectionVectors(5*M_PI/4);
			Unbind();
			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
				genPositionVectors();
			Unbind();
			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
				genMagnitudeVectors();
			Unbind();
			Bind(btexbufferID, bfboID, 27);
			DrawCylinder();
			Unbind();
		}
		
//--------------------------------------------------------------------------//
//-----------------------------ACCESS METHODS-------------------------------//
//--------------------------------------------------------------------------//
		///Returns the p vector.
		vec3<float>
		getPosition()
		{
			return p;
		}
	
		///Returns the d vector.
		vec3<float>
		getDirection()
		{
			return d;
		}

		///Returns the m vector.
		stim::vec<float>
		getMagnitude()
		{
			return m;
		}
	
		///@param stim::vec<float> pos, the new p.
		///Sets the p vector to input vector pos.
		void
		setPosition(stim::vec3<float> pos)
		{
			p = pos;
		}

		///@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)
		{
			p[0] = x;
			p[1] = y;
			p[2] = z;
		}
		
		///@param stim::vec<float> dir, the new d.
		///Sets the d vector to input vector dir.
		void
		setDirection(stim::vec3<float> dir)
		{
			d = dir;
		}
		
		///@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)
		{
			d[0] = x;
			d[1] = y;
			d[2] = z;
		}
			
		///@param stim::vec<float> dir, the new d.
		///Sets the m vector to the input vector mag.	
		void
		setMagnitude(stim::vec<float> mag)
		{
			m[0] = mag[0];
			m[1] = mag[0];
		}
		
		///@param float mag, size of the sampled region.
		///Sets the m vector to the input mag for both templates.
		void
		setMagnitude(float mag)
		{
			m[0] = mag;
			m[1] = mag;
		}
		
		///@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)
		{
			S[0] = x;
			S[1] = y;
			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::vec3<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)
		{
			R[0] = x;
			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::vec3<float> Siz)
		{
			S = Siz;
		}
		
		///@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::vec3<float> dir)
		{
			stim::vec<float> out(0.0,0.0,0.0,0.0);
			stim::vec<float> from(0.0,0.0,1.0);
			out[0] = acos(dir.dot(from))*180/M_PI;
			if(out[0] < 1.0){
				out[0] = 0.0;
				out[1] = 0.0;
				out[2] = 0.0;
				out[3] = 1.0;
			} else {
				stim::vec<float> temp(0.0, 0.0, 0.0);;
				stim::vec<float> dir1(dir[0], dir[1], dir[2]);
				temp = (from.cross(dir1)).norm();
				out[1] = temp[0];
				out[2] = temp[1];
				out[3] = temp[2];
			}
			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::vec3<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::vec3<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)
		{
//			std::cout << "sMag: " << mag << std::endl;
			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)
		{
//			std::cout << "sPos: " << x << " " << y << " " << z << std::endl;
			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)
		{
//			std::cout << "sDir: " << x << " " << y << " " << z << std::endl;
			seedsvecs.push(stim::vec<float>(x, y, z));
		}
		
		///Method to get the top of the seed positions stack.
		stim::vec3<float> 
		getLastSeed()
		{
			stim::vec3<float> tp = seeds.top();
			return tp;
		}

		///Method to get the top of the seed direction stack.
		stim::vec3<float> 
		getLastSeedVec()
		{
			stim::vec3<float> tp = seedsvecs.top();
			return tp;
		}

		///Method to get the top of the seed magnitude stack.
		float
		getLastSeedMag()
		{
			float tp = seedsmags.top();
			return tp;
		}

		///deletes all data associated with the last seed.
		void
		popSeed()
		{
			seeds.pop();
			seedsvecs.pop();
			seedsmags.pop();
		}
		
		///returns the seeds position stack.
		std::stack<stim::vec3<float> >
		getSeeds()
		{
			return seeds;
		}

		///sets the number of direction templates.
		void
		setNumberDirectionTemplates(int n)
		{
			numSamples = n;
		}

		///sets the number of position templates.
		void
		setNumberPositionTemplates(int n)
		{
			numSamplesPos = n;
		}

		///sets the number of position templates.
		void
		setNumberSizeTemplates(int n)
		{
			numSamplesMag = n;
		}

		#ifdef TIMING
		///Returns the timings at the moment the method is called.
		///In the following order: Branch, Direction, Position, Size, Cost, Network, Hit_detetion.
		std::vector<double> 
		getTimings()
		{
			std::vector <double> ret;
			ret.resize(7);
			ret[0] = branch_time;
			ret[1] = direction_time;
			ret[2] = position_time;
			ret[3] = size_time;
			ret[4] = cost_time;
			ret[5] = network_time;
			ret[6] = hit_time;

			return ret;
		}
		#endif

		///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;
			}
		}
		
		///Saves the network to a file.
		void
		saveNetwork(std::string name)
		{
			stim::glObj<float> sk;
			for(int i = 0; i < nt.sizeE(); i++)
			{
				std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
                 		std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
				sk.Begin(stim::OBJ_LINE);
				for(int j = 0; j < ce.size(); j++)
				{
					sk.TexCoord(cm[j][0]);
					sk.Vertex(ce[j][0], ce[j][1], ce[j][2]);
				}
				sk.End();
			}	
			sk.save(name);
		}

		///Depreciated, but might be reused later()
		///returns a COPY of the entire stim::glObj object.
		stim::glObj<float>
		getNetwork()
		{

		}

		///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----------------------------//
//--------------------------------------------------------------------------//

		///temporary Method necessary for visualization and testing.
		void
		Update()
		{
			vec3<float> Y(1.0,0.0,0.0);
			if(cos(Y.dot(d))< 0.087){
				Y[0] = 0.0; Y[1] = 1.0;}
			hor = stim::rect<float>(m, p, d.norm(),
				((Y.cross(d)).cross(d)).norm());
			ver = stim::rect<float>(m, p, d.norm(),
				 hor.n());
		}


		int
		Step()
		{
			Bind(texbufferID, fboID, numSamples, t_length);
			CHECK_OPENGL_ERROR
				findOptimalDirection();
			Unbind();
			Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
				findOptimalPosition();
			Unbind();
			Bind(mtexbufferID, mfboID, numSamplesMag, t_length);
				findOptimalScale();
			Unbind();
			CHECK_OPENGL_ERROR

			return current_cost;
		}


		void
		printTransform()
		{
			std::cout << cT << std::endl;
		}


//--------------------------------------------------------------------------//
//-----------------------------EXPERIMENTAL METHODS-------------------------//
//--------------------------------------------------------------------------//

		void
		MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
		{
			float PHI[2], Z[2], range;
			PHI[0] = asin(solidAngle/2);
			PHI[1] = asin(0);
			
			Z[0] = cos(PHI[0]);
			Z[1] = cos(PHI[1]);
			
			range = Z[0] - Z[1];

			std::vector<stim::vec3<float> > vecsUni;
			for(int i = 0; i < numSamplesPos; i++)
			{
				stim::vec3<float> a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0);
				a[0] = a[0]-0.4;
				a[1] = a[1]-0.4;
				vecsUni.push_back(a);
			}

			stringstream name;
			for(int i = 0; i < numSamplesPos; i++)
				name << vecsUni[i].str() << std::endl;
			
			std::ofstream outFile;
			outFile.open("New_Pos_Vectors.txt");
			outFile << name.str().c_str();
		}
		void
		DrawCylinder()
		{	 
			 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.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*6.4+6.4);
//					 glVertex2f(0.0, j*27.0+27.0);
					 glTexCoord3f(x,y,z1); 
					 glVertex2f(16.0, j*6.4+6.4);
//					 glVertex2f(16.0, j*27.0+27.0);
					 glTexCoord3f(xold,yold,z1); 
					 glVertex2f(16.0, j*6.4); 
//					 glVertex2f(16.0, j*27.0); 
					 glTexCoord3f(xold,yold,z0); 
					 glVertex2f(0.0, j*6.4);
//					 glVertex2f(0.0, j*27.0);
					 xold=x;
					 yold=y;
					 j++;
				}
		      	 glEnd();  
			 glEndList();
		}

///need to return the cylinder.
		void
		DrawLongCylinder(int n = 8, int l_template = 8,int l_square = 8)
		{
			int cylLen = cL.size()-1;
			GenerateFBO(n*l_square, cylLen*l_template, btexbufferID, bfboID);
			Bind(btexbufferID, bfboID, cylLen, l_template*l_square/2.0);
			stim::cylinder<float> cyl(cL, cM);
			std::vector<std::vector<stim::vec3<float> > > p = cyl.getPoints(n);
			for(int i = 0; i < p.size()-1; i++)	///number of circles
			{
				for(int j = 0; j < p[0].size()-1; j++)		///points in the circle
				{
					glBegin(GL_QUADS);
						glTexCoord3f(p[i][j][0], p[i][j][1], p[i][j][2]);
						glVertex2f(j*l_square,  i*(float)l_template); 

						glTexCoord3f(p[i][j+1][0], p[i][j+1][1], p[i][j+1][2]);
						glVertex2f(j*l_square+l_square, i*(float)l_template); 

						glTexCoord3f(p[i+1][j+1][0], p[i+1][j+1][1], p[i+1][j+1][2]);
						glVertex2f(j*l_square+l_square, i*(float)l_template+(float)l_template); 

						glTexCoord3f(p[i+1][j][0], p[i+1][j][1], p[i+1][j][2]);
						glVertex2f(j*l_square,i*(float)l_template+(float)l_template); 
					glEnd();
				}
			}
			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)
		{	
			stim::vec3<float> curSeed; 
			stim::vec3<float> curSeedVec;			
			float curSeedMag;
			while(!Empty())
			{
				//clear the currently traced line and start a new one.
				cL.clear();
				cM.clear();
				cD.clear();
				curSeed = seeds.top();
				curSeedVec = seedsvecs.top();
				curSeedMag = seedsmags.top();
				seeds.pop();
				seedsvecs.pop();
				seedsmags.pop();
//				std::cout << "The current seed Vector is " << curSeedVec << std::endl;
				setPosition(curSeed);
				setDirection(curSeedVec);
				setMagnitude(curSeedMag);
//				cL.push_back(curSeed);
//				cM.push_back(curSeedMag);
//				cD.push_back(curSeedMag);
				pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost);
			}
		}

		int
		selectObject(stim::vec3<float> loc, stim::vec3<float> dir, float mag) 
		{
		//Define the varibles and turn on Selection Mode
			
			#ifdef TIMING
				gpuStartTimer();
			#endif

//			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/stepsize);
				camSel.LookAt((loc[0]+dir[0]*mag/stepsize),
					 (loc[1]+dir[1]*mag/stepsize),
					 (loc[2]+dir[2]*mag/stepsize));
				ps  = camSel.getPosition();
				ups = camSel.getUp();
				ds  = camSel.getLookAt();
				glMatrixMode(GL_PROJECTION);
				glPushMatrix();
				glLoadIdentity();
				glOrtho(-mag/stepsize/2.0, mag/stepsize/2.0, -mag/stepsize/2.0, mag/stepsize/2.0, 0.0, mag/stepsize/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);
			#ifdef TIMING
				hit_time += gpuStopTimer();
			#endif

			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 *ptr;
			ptr = (GLuint *) buffer;
			ptr++;
			ptr++; //Skip the minimum depth value.  
			ptr++; //Skip the maximum depth value.

			
			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::vec3<float> spos, 
				stim::vec<float> smag, stim::vec3<float> sdir)
		{
			#ifdef TIMING
				 double s = std::clock();
			#endif
			
                        std::vector<stim::vec3<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 && h < nt.sizeE())
						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 && h < nt.sizeE()){	
			//			std::cout << "got here double" << smag.str() << std::endl;
						nt.addEdge(ce,cm, h, in.second);	
					} else { nt.addEdge(ce,cm, -1, -1);}
				}
			}		

			#ifdef TIMING
				double nt = (std::clock() - s) / (double) CLOCKS_PER_SEC;
				network_time += nt * 1000.0;
			#endif
		}

//		void
//		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos, 
//				stim::vec<float> smag, stim::vec3<float> sdir)
//		{
//
//		}


		void
		printSizes()
		{
			std::cout << nt.sizeE() << " edges " << std::endl;
			std::cout << nt.sizeV() << " nodes " << std::endl;
			
		}

		std::pair<stim::fiber<float>, int >
		traceLine(stim::vec3<float> pos, stim::vec<float> mag, int min_cost)
		{
			//starting (seed) position and magnitude.
			stim::vec3<float> spos = getPosition();
			stim::vec<float> smag = getMagnitude();
			stim::vec3<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::vec3<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();
					branchDetection2();
					pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
					addToNetwork(a, spos, smag, sdir);
//					std::cout << "the cost of " << cost << " > " << min_cost << std::endl;
					return a;
					break;
				} else {
					//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)
					{
						running = false;
						branchDetection2();
						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
						addToNetwork(a, spos, smag, sdir);
//						std::cout << "I hit and edge" << std::endl;
						return a;
						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;
					}
					//Has the template size gotten unreasonable?
					mag = getMagnitude();
					if(mag[0] > 75 || mag[0] < 1){
						running = false;
						branchDetection2();
						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
						addToNetwork(a, spos, smag, sdir);
//						std::cout << "The templates are too big" << std::endl;
						return a;
						break;
					}
					else
					{
						h = selectObject(p, getDirection(), m[0]);
						//Have we hit something previously traced?
						if(h != -1){
//							std::cout << "I hit the fiber " << h << std::endl;
							running = false;
							branchDetection2();
							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
							addToNetwork(a, spos, smag, sdir);
							return a;
							break;
						}
						else {  
							cL.push_back(stim::vec3<float>(p[0], p[1],p[2]));
							cM.push_back(stim::vec<float>(m[0], m[0]));
//							Bind(btexbufferID, bfboID, 27);
							Unbind();
							CHECK_OPENGL_ERROR
							
						}
				 	}
                 		}
         		}
		}
};
}
#endif