codebase / stimlib

#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 <stack>
#include <stim/cuda/branch_detection.cuh>
#include "../../../volume-spider/glnetwork.h"
#include <stim/visualization/cylinder.h>
#include <stim/cuda/testKernel.cuh>
#include <iostream>
#include <fstream>
#ifdef TIMING
	#include <ctime>
	#include <cstdio>	
#endif

#ifdef TESTING
	#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;	//normalized direction of travel
		float m;			//size of the spider in tissue space.

		std::vector<stim::vec3<float> > dV;	//A list of all the direction vectors.
		std::vector<stim::vec3<float> > pV;	//A list of all test positions (relative to p)
		std::vector<float> mV;	//A list of all the size vectors.

		stim::matrix<float, 4> cT;			//current Transformation matrix (tissue)->(texture)
		GLuint texID;						//OpenGL ID for the texture to be traced
		stim::vec3<float> S;					//Size of a voxel in the volume.
		stim::vec3<float> R;					//Dimensions of the volume.


		//GL and Cuda variables
		GLuint dList;				//ID of the starting display lists (series of 4)
									//dList + 0 = direction template rectangles
									//dList + 1 = position template rectangles
									//dList + 2 = size template rectangles
									//dList + 3 = branch detection cylinder around the fiber

		GLuint fboID;				//framebuffer ID for direction templates
		GLuint texbufferID;			//texture ID for direction templates
		GLuint direction_buffID;				//framebuffer ID, position templates
		GLuint direction_texID;		//texture ID, position templates

		GLuint position_buffID;				//framebuffer ID, position templates
		GLuint position_texID;		//texture ID, position templates

		GLuint radius_buffID;				//framebuffer ID, radius templates
		GLuint radius_texID;		//texture ID, radius templates

		GLuint cylinder_buffID;				//framebuffer ID, cylinder (surrounding fiber)
		GLuint cylinder_texID;		//texture ID, cylinder

		int numSamples;				//The number of templates in the buffer.
		int numSamplesPos;
		int numSamplesMag;

		float length;				//this will be a function of the radius
		float stepsize;				//this will be a function of the length

		int current_cost;			//variable to store the cost of the current step		
		
		//Tracing variables.
		std::stack< stim::vec3<float> > seeds;		//seed positions
		std::stack< stim::vec3<float> > seedsvecs;	//seed directions
		std::stack< float > 		seedsmags;		//seed magnitudes

		std::vector< stim::vec3<float> > cL;	//centerline up to the current point
		std::vector< stim::vec3<float> > cD;	//directions up to the current point (debugging)
		std::vector< float > cM;		//radius up to the current point

		stim::glnetwork<float> nt;				//network object holding the currently traced centerlines
		stim::glObj<float> sk;					//OBJ file storing the network (identical to above)

		//consider replacing with two seed points facing opposite directions
		stim::vec<float> rev;					//reverse vector

		//selection mode - detecting fiber intersections
		stim::camera camSel;					//camera for selection mode (detecting collisions)
		stim::vec3<float> ps;					//position for the selection camera
		stim::vec3<float> ups;					//up direction for the selection camera
		stim::vec3<float> ds;					//direction for the selection camera

		float n_pixels;							//length of the template (in pixels)

		//cuda texture variables that keep track of the binding.
		stim::cuda::cuda_texture t_dir;				//cuda_texture object used as an interface between OpenGL and cuda for direction vectors.
		stim::cuda::cuda_texture t_pos;				//cuda_texture object used as an interface between OpenGL and cuda for position vectors.
		stim::cuda::cuda_texture t_mag;				//cuda_texture object used as an interface between OpenGL and cuda for size vectors.
		

		#ifdef DEBUG
			stringstream name;
			int iter;
			int iter_pos;
			int iter_dir;
			int iter_siz;
		#endif

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

		/// 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();	//Timer for profiling
			#endif
			setMatrix();			//create the transformation matrix.
			glCallList(dList);		//move the templates to p, d, m.
			glFinish();			//flush the pipeline
			#ifdef TIMING
				direction_time += gpuStopTimer();	//profiling
			#endif

			int best = getCost(t_dir.getTexture(), t_dir.getAuxArray() ,numSamples);		//find min cost.
			#ifdef DEBUG
				name.str("");
				name << "Final_Cost_Direction_fiber_"<< iter << "_" << iter_dir << ".bmp";	
				test(t_dir.getTexture(), n_pixels*2.0, numSamples*n_pixels, name.str());
				iter_dir++;
			#endif
			stim::vec<float> next(		///calculate the next vector.
			dV[best][0]*S[0]*R[0],
			dV[best][1]*S[1]*R[1],
			dV[best][2]*S[2]*R[2],
			0);
			next = (cT*next).norm();	///transform the next vector into Tissue space.
			setPosition(	p[0]+next[0]*m/stepsize,
					p[1]+next[1]*m/stepsize,
					p[2]+next[2]*m/stepsize);	
			setDirection(next[0], next[1], next[2]);		//move forward and change direction.
		}

		/// Method for finding the best d (direction) for the spider.
		/// Not sure if necessary since the next p (position) for the spider
		/// will be at d * m.
		void
		findOptimalPosition()
		{
			#ifdef TIMING
				gpuStartTimer();	//timer for profiling
			#endif
			setMatrix();			//create the transformation matrix.
			glCallList(dList+1);		//move the templates to p, d, m.
			glFinish();			//flush the pipeline
//			glFlush();
			#ifdef TIMING
				position_time += gpuStopTimer();		///timer for profiling
			#endif

			int best = getCost(t_pos.getTexture(), t_pos.getAuxArray(), numSamplesPos);		//find min cost.
			#ifdef DEBUG
				name.str("");
				name << "Final_Cost_Position_" << iter << "_" << iter_pos << ".bmp";	
				test(t_pos.getTexture(), n_pixels*2.0, numSamplesPos*n_pixels, name.str());
				iter_pos++;
			#endif
			stim::vec<float> next(		//find next position.
 				pV[best][0],
				pV[best][1],
				pV[best][2],
				1);
			next = cT*next;			//transform the next position vector into tissue space.
			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();			//flush the drawing pipeline.
			#ifdef TIMING
				size_time += gpuStopTimer();
			#endif
			int best = getCost(t_mag.getTexture(), t_mag.getAuxArray(), numSamplesMag);		//get best cost.
			#ifdef DEBUG
				name.str("");
				name << "Final_Cost_Size_" << iter << "_" <<  iter_siz << ".bmp";	
				test(t_mag.getTexture(), n_pixels*2.0, numSamplesMag*n_pixels, name.str());
				iter_siz++;
			#endif
			setMagnitude(m*mV[best]);	//adjust the magnitude.
		}




		///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();		///timer for performance analysis
			#endif

			if(cL.size() < 4){}		///if the size of the fiber is less then 4 we do nothing.
			else{
			setMatrix(1);			///finds the current transformation matrix
			DrawLongCylinder(n, l_template, l_square);	///Draw the cylinder.
			stim::cylinder<float> cyl(cL, cM);
			std::vector< stim::vec<float> > result = find_branch(cylinder_texID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);		///find all the centers in cuda
			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);			///the borders of the texture.
			float pval;									//pvalue associated with the points on the cylinder.
			if(!result.empty())								///if we have any points
			{
				for(int i = 0; i < result.size(); i++)					///for each point
				{
					int id = result[i][2];
					if(fmod(result[i][2], id) != 0 && id != 0)			///if the remainer is odd
					{
						
						pval = ((cyl.getl(id+1)-cyl.getl(id))*
							(fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);		///calculate pvalue
					}
					else if(id == 0)									///if the point is on the edge
					{
						pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);			
					}
					else
					{
						pval = (cyl.getl(id)/cyl.getl(cL.size()-1));					///if the point is somewhere on the surface of the cylinder other than the edge
					}
					stim::vec3<float> v = cyl.surf(pval, result[i][0]);					///find the coordinates of the point at pval on the surface in tissue space.
					stim::vec3<float> di = cyl.p(pval);							///find the coord of v in tissue space projected on the centerline.
					float rad = cyl.r(pval)/2;								///find the radius at the pvalue's location
					if(
					 !(v[0] > size[0] || v[1] > size[1]
					 || v[2] > size[2] || v[0] < 0
					 || v[1] < 0 || v[2] < 0))								///if the v point is INSIDE the volume
					{
						setSeed(v);									///add a seedpoint's position.
						setSeedVec((v-di).norm());							///add a seedpoints direction
						setSeedMag(rad);								///add the starting radius.
					}
				}
			}
			}
			#ifdef TIMING
				branch_time += gpuStopTimer();									///timer for performance.
			#endif
		}


		float uniformRandom()
		{
			return (  (float)(rand()))/(  (float)(RAND_MAX));							///generates a random number between 0 and 1 using the uniform distribution.
		}

		float normalRandom()
		{
			float u1 = uniformRandom();
			float u2 = uniformRandom();
			return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));							///generate a random number using the normal distribution between 0 and 1.
		}

		stim::vec3<float> uniformRandVector()
		{
			stim::vec3<float> r(uniformRandom(), uniformRandom(), 1.0);						///generate a random vector using the uniform distribution between 0 and 1.
			return r;
		}

		stim::vec3<float> normalRandVector()
		{
			stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);						///generate a random vector using the normal distribution between 0 and 1.
			return r;
		}

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



		///@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 = stim::PI/2)
		{

			//Set up the vectors necessary for Rectangle creation.
			stim::vec3<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			stim::vec3<float> pos(0.0,0.0,0.0);		//center point of a rectangle
			float mag = 1.0;				//size of the generated rectangle.
			stim::vec3<float> dir(0.0, 0.0, 1.0);		//normal of the rectangle

			float PHI[2], Z[2], range;			
			PHI[0] = solidAngle/2;				///Project the solid angle into spherical coordinates
			PHI[1] = asin(0);
			
			Z[0] = cos(PHI[0]);				///Project the z into spherical coordinates
			Z[1] = cos(PHI[1]);
			
			range = Z[0] - Z[1];				///The range the possible values can be.

			float z, theta, phi;
			glNewList(dList, GL_COMPILE); 			///create a display list of all the direction templates.
			for(int i = 0; i < numSamples; i++)		///for each sample
			{
				z = uniformRandom()*range + Z[1];	///generate a z coordinate
				theta = uniformRandom()*stim::TAU;	///generate a theta coordinate
				phi = acos(z);				///generate a phi from the z.
				stim::vec3<float> sph(1, theta, phi);	///combine into a vector in spherical coordinates.
				stim::vec3<float> cart = sph.sph2cart();///convert to cartesian.	
				dV.push_back(cart);			///save the generated vector for further use.
				#ifdef DEBUG
//					std::cout << cart << std::endl;
				#endif
				if(cos(Y.dot(cart)) < 0.087)		///make sure that the Y is not parallel to the new vector.
				{
						Y[0] = 0.0; Y[1] = 1.0;
				}else{
					Y[0] = 1.0; Y[1] = 0.0;
				}

                                hor = stim::rect<float>(mag,		///generate a rectangle with the new vectro as a normal.
					 pos, cart,
                       			((Y.cross(cart)).cross(cart)).norm());

				#ifdef DEBUG
			//		std::cout << hor.n() << std::endl;
				#endif
       				ver = stim::rect<float>(mag,		///generate another rectangle that's perpendicular the first but parallel to the cart vector.
					 pos, cart,
                       			hor.n());
				UpdateBuffer(0.0, 0.0+i*n_pixels);	///Put the necessary points into the diplaylist.
			}
			glEndList();					///finilize the display list.
		}

		///@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);		//center point of a rectangle
			float mag = 1.0;				///size of each rectangle
			stim::vec3<float> dir(0.0, 0.0, 1.0);		///normal of the rectangle plane.

			//Set up the variable necessary for vector creation.
			glNewList(dList+1, GL_COMPILE);			///generate a new display list.
			pV.push_back(pos);
			hor = stim::rect<float>(mag,                    ///generate a rec     tangle with the new vector as a normal.  
                	        pos, dir,
	                       	((Y.cross(d)).cross(d))
	                        .norm());
                        ver = stim::rect<float>(mag,                    ///generate anoth     er rectangle that's perpendicular the first but parallel to the cart vector.
        	                pos, dir,
	                        hor.n()); 
                        UpdateBuffer(0.0, 0.0+0*n_pixels);
			for(int i = 1; i < numSamplesPos; i++)		///for the number of position samples
			{
				stim::vec3<float> temp = uniformRandVector();	///generate a random point on a plane.
				temp = temp*delta*2.0 - delta;		///scale the point and center it around 0.0
				temp[2] = 0.0;
				pV.push_back(temp);				///save the point for further use.
                		hor = stim::rect<float>(mag,			///generate a rectangle with the new vector as a normal.
					 temp, dir,
                      			((Y.cross(d)).cross(d))
					.norm());
                		ver = stim::rect<float>(mag,			///generate another rectangle that's perpendicular the first but parallel to the cart vector.
					 temp, dir,
                        		hor.n());
				UpdateBuffer(0.0, 0.0+i*n_pixels);		///sample the necessary points and put them into a display list.
			}
			glEndList();					///finilize the display list.
			#ifdef DEBUG
				for(int i = 0; i < numSamplesPos; i++)
					std::cout << pV[i] << std::endl;
			#endif
		}

		///@param float delta, How much the rectangles are allowed to expand.
		///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);		//center of the future rect.
			float mag = 1.0;				///size of the rectangle
			stim::vec3<float> dir(0.0, 0.0, 1.0);		///normal of the rectangle plane.

			//Set up the variable necessary for vector creation.
			float min 	= 1.0-delta;			///smallest size
			float max 	= 1.0+delta;			///largers size.
			float step	= (max-min)/(numSamplesMag-1);	///the size variation from one rect to the next.
			float factor;
			glNewList(dList+2, GL_COMPILE);
			for(int i = 0; i < numSamplesMag; i++){		///for the number of position samples
				//Create linear index
				factor = (min+step*i)*mag;		///scaling factor
				mV.push_back(factor);			///save the size factor for further use.
				hor = stim::rect<float>(factor,		///generate a rectangle with the new vector as a normal.
					 pos, dir, 
       	       				((Y.cross(d)).cross(d))
					.norm());
               			ver = stim::rect<float>(factor,		///generate another rectangle that's perpendicular the first but parallel to the cart vector.
					 pos, dir,
                       			hor.n());
				UpdateBuffer(0.0, 0.0+i*n_pixels);	///sample the necessary points and put them into a display list.
			CHECK_OPENGL_ERROR
			}
			glEndList();					///finilize the displaylist.
		}

		///@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; 				///first point.
        	        stim::vec3<float>p2; 				///second point.
	                stim::vec3<float>p3; 				///third point.
                	stim::vec3<float>p4;				///fourth point.
			p1 = hor.p(1,1);				///generate the top right point from the horizontal template.
			p2 = hor.p(1,0);				///generate the bottom right point from the horizonatal template.
			p3 = hor.p(0,0);				///generate the bottom left point from the horizontal template.
			p4 = hor.p(0,1);				///generate the top left point from the horizonatal template.
			glBegin(GL_QUADS);				///generate the Quad from the 4 points.
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x,v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+n_pixels, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+n_pixels, v_y+n_pixels);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x, v_y+n_pixels);
			 glEnd();					///finish the quad.

			 p1 = ver.p(1,1);				///generate the top right point from the vertical template.
			 p2 = ver.p(1,0);                               ///generate the bottom right point from the vertical template.
			 p3 = ver.p(0,0);                               ///generate the bottom left point from the vertical template.
			 p4 = ver.p(0,1);                               ///generate the top left point from the vertical template.
		 	 glBegin(GL_QUADS);				///generate the Quad from the 4 points.
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x+n_pixels, v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+2.0*n_pixels, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+2.0*n_pixels, v_y+n_pixels);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x+n_pixels, v_y+n_pixels);
			glEnd(); 					///finish the quad.
		}
		


//--------------------------------------------------------------------------//
//--------------------------------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);				///clear the framebuffer.
			glGenFramebuffers(1, &framebufferID);					///generate a clean buffer.
			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);			///bind the new buffer.
//			int numChannels = 1;
//			unsigned char* texels = new unsigned char[width * height * numChannels];
			glGenTextures(1, &textureID);						///generate a texture that will attach to the buffer.
			glBindTexture(GL_TEXTURE_2D, textureID);

			//Textures repeat and use linear interpolation, luminance format.
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);		///Set up the texture to repeat at edges.
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);		
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);	///Set up the texture to use Linear interpolation
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,				///Create the texture with no data.
				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, NULL);   
//			delete[] texels;
			glBindFramebuffer(GL_FRAMEBUFFER, 0); 					///Bind the frontbuffer
			glBindTexture(GL_TEXTURE_2D, 0);					///Unbind the texture.
		}


		///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,
					 m,
					 m);
				//get and store the current transformation matrix for later use.
				glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
				cT.set(curTrans);
				
				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, direction_buffID);//set up GL buffer		
			glFramebufferTexture2D(
				GL_FRAMEBUFFER,
				GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_2D,
				direction_texID,
				0);
			glBindFramebuffer(GL_FRAMEBUFFER, direction_buffID);
			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
			glDrawBuffers(1, DrawBuffers);
			glBindTexture(GL_TEXTURE_2D, direction_texID);
			glClearColor(1,1,1,1);
			glClear(GL_COLOR_BUFFER_BIT);
			glMatrixMode(GL_PROJECTION);
			glLoadIdentity();
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glViewport(0,0,2.0*n_pixels, numSamples*n_pixels);
			gluOrtho2D(0.0,2.0*n_pixels,0.0,numSamples*n_pixels);
			glEnable(GL_TEXTURE_3D);
			glBindTexture(GL_TEXTURE_3D, texID);

			CHECK_OPENGL_ERROR
		}

		///Method for controling the buffer and texture binding.
		///Clears the buffer upon binding.
		///@param GLuint &textureID, texture to be bound.
		///@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);	///Bind the framebuffer.
			glFramebufferTexture2D(					///associate it with the texture
				GL_FRAMEBUFFER,
				GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_2D,
				textureID,
				0);
			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);	///Bind the framebuffer.
			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};		///generate the drawbuffer.
			glDrawBuffers(1, DrawBuffers);				///set the drawbuffer.
			glBindTexture(GL_TEXTURE_2D, textureID);		///Bind the texture passed.
			glMatrixMode(GL_PROJECTION);				///clear out the draw matrices
			glLoadIdentity();
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glViewport(0,0,2.0*len, nSamples*len);			///set up viewport
			gluOrtho2D(0.0,2.0*len,0.0,nSamples*len);		///Set up ortho
			glEnable(GL_TEXTURE_3D);				
			glBindTexture(GL_TEXTURE_3D, texID);			///bind the main texture (return to original state).

			CHECK_OPENGL_ERROR
		}
		
		///Unbinds all texture resources.
		void
		Unbind()
		{
			//Finalize GL_buffer
			glBindTexture(GL_TEXTURE_3D, 0);                      ///Bind the front buffer.
			CHECK_OPENGL_ERROR
			glBindTexture(GL_TEXTURE_2D, 0);                      ///Bind the default GL texture.
			CHECK_OPENGL_ERROR
			glBindFramebuffer(GL_FRAMEBUFFER, 0);		      ///Bind the defautl framebuffer.
			CHECK_OPENGL_ERROR
			glDisable(GL_TEXTURE_3D);			      ///Turn off texturing.
			CHECK_OPENGL_ERROR
		}

		


//--------------------------------------------------------------------------//
//--------------------------------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(cudaTextureObject_t tObj, float* result, int n)
		{
			#ifdef TIMING
				gpuStartTimer();
			#endif
			stim::vec<int> cost = 
				stim::cuda::get_cost(tObj, result, n, 2*n_pixels, n_pixels);	///call the cuda function with the appropriate texture buffer.
			#ifdef TIMING
				cost_time += gpuStopTimer();
			#endif
			current_cost = cost[1];							///current cost.
			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)
		{
			p = stim::vec3<float>(0.0, 0.0, 0.0);
			d = stim::vec3<float>(0.0, 0.0, 1.0);
			m = 1.0;
			S = stim::vec3<float>(1.0, 1.0, 1.0);
			R = stim::vec3<float>(1.0, 1.0, 1.0);
			numSamples = samples;
			numSamplesPos = samplespos;
			numSamplesMag = samplesmag;
		#ifdef DEBUG
			iter = 0;
			iter_pos = 0;
			iter_dir = 0;
			iter_siz = 0;
		#endif
		}

		///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 = mag_x;
			S = stim::vec3<float>(1.0,1.0,1.0);
			R = stim::vec3<float>(1.0,1.0,1.0);
			numSamples = numsamples;
			numSamplesPos = numsamplespos;
			numSamplesMag = numsamplesmag;
		#ifdef DEBUG
			iter = 0;
			iter_pos = 0;
			iter_dir = 0;
			iter_siz = 0;
		#endif
		}

		///Position constructor: vecs of floats.
		///@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 = mag;
			S = vec3<float>(1.0,1.0,1.0);
			R = vec3<float>(1.0,1.0,1.0);
			numSamples = samples;
			numSamplesPos = samplesPos;
			numSamplesMag = samplesMag;
		#ifdef DEBUG
			iter = 0;
			iter_pos = 0;
			iter_dir = 0;
			iter_siz = 0;
		#endif
		}

		///destructor	
		~gl_spider
		(void)
		{
			Unbind();
			glDeleteTextures(1, &direction_texID);
			glDeleteBuffers(1, &direction_buffID);
			glDeleteTextures(1, &position_texID);
			glDeleteBuffers(1, &position_buffID);
			glDeleteTextures(1, &radius_texID);
			glDeleteBuffers(1, &radius_buffID);     
			glDeleteTextures(1, &cylinder_texID);
			glDeleteBuffers(1, &cylinder_buffID);
		}

		///@param GLuint id, texture that is going to be sampled.
		///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
#ifdef DEBUG
			iter = 0;
			iter_pos = 0;
			iter_dir = 0;
			iter_siz = 0;
#endif
			stepsize = 3.0;
			n_pixels = 16.0;

			srand(100);	
			texID = id;
			GenerateFBO(n_pixels*2, numSamples*n_pixels, direction_texID, direction_buffID); 
			CHECK_OPENGL_ERROR
			GenerateFBO(n_pixels*2, numSamplesPos*n_pixels, position_texID, position_buffID); 
			CHECK_OPENGL_ERROR
			GenerateFBO(n_pixels*2, numSamplesMag*n_pixels, radius_texID, radius_buffID); 
			CHECK_OPENGL_ERROR
			GenerateFBO(16, 216, cylinder_texID, cylinder_buffID); 
			CHECK_OPENGL_ERROR
			t_dir.MapCudaTexture(direction_texID, GL_TEXTURE_2D);
			t_dir.Alloc(numSamples);
			t_pos.MapCudaTexture(position_texID, GL_TEXTURE_2D);
			t_pos.Alloc(numSamplesPos);
			t_mag.MapCudaTexture(radius_texID, GL_TEXTURE_2D);
			t_mag.Alloc(numSamplesMag);
			setMatrix();
			dList = glGenLists(3);
			glListBase(dList);
			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
				genDirectionVectors(5*stim::PI/4);
			Unbind();
			Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
				genPositionVectors(0.2);
			Unbind();
			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
				genMagnitudeVectors();
			Unbind();
			Bind(cylinder_texID, cylinder_buffID, 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.
		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 float mag, size of the sampled region.
		///Sets the m vector to the input mag for both templates.
		void
		setMagnitude(float mag)
		{
			m = 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::vec3<float> from(0.0,0.0,1.0);
			out[0] = acos(dir.dot(from))*180/stim::PI;
			#ifdef DEBUG
				std::cout << "out is " << out << std::endl;
				std::cout << "when rotating from " << from << " to " << dir << std::endl;
			#endif
			if(out[0] < 0.01){
				out[0] = 0.0;
				out[1] = 0.0;
				out[2] = 0.0;
				out[3] = 1.0;
			}
			else if(out[0] < -180.0+0.01)
			{
				out[0] = 180.0;
				out[1] = 1.0;
				out[2] = 0.0;
				out[3] = 0.0;
			} else {
				stim::vec3<float> temp(0.0, 0.0, 0.0);;
				temp = (from.cross(dir)).norm();
				out[1] = temp[0];
				out[2] = temp[1];
				out[3] = temp[2];
			}
			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)
		{
			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));
		}
		
		///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()
		{
			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 cylinder_buffID;
		}

//--------------------------------------------------------------------------//
//-----------------------------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()
		{
			#ifdef DEBUG
			std::cerr << "Took a step" << std::endl;
			#endif
			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
			CHECK_OPENGL_ERROR
				findOptimalDirection();
			Unbind();
			Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
				findOptimalPosition();
			Unbind();
			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
				findOptimalScale();
			Unbind();
			CHECK_OPENGL_ERROR

			return current_cost;
		}


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


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

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

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

			stringstream name;
			for(int i = 0; i < numSamplesPos; i++)
				name << vecsUni[i].str() << std::endl;
			
			std::ofstream outFile;
			outFile.open("New_Pos_Vectors.txt");
			outFile << name.str().c_str();
		}
		void
		DrawCylinder()
		{	 
			 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*stim::TAU/360.0);
					 y=r0*sin(i*stim::TAU/360.0);
					 glTexCoord3f(x,y,z0); 
					 glVertex2f(0.0, j*6.4+6.4);
//					 glVertex2f(0.0, j*27.0+27.0);
					 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.
///SOMETHING MIGHT BE GOING ON HERE IN GENERATE BUFFER.
		void
		DrawLongCylinder(int n = 8, int l_template = 8,int l_square = 8)
		{
			int cylLen = cL.size()-1;
			GenerateFBO(n*l_square, cylLen*l_template, cylinder_texID, cylinder_buffID);
			Bind(cylinder_texID, cylinder_buffID, cylLen, l_template*l_square/2.0);
			stim::cylinder<float> cyl(cL, cM);
			std::vector<std::vector<stim::vec3<float> > > p = cyl.getPoints(n);
			for(int i = 0; i < p.size()-1; i++)	///number of circles
			{
				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();
				setPosition(curSeed);
				setDirection(curSeedVec);
				setMagnitude(curSeedMag);

				#ifdef DEBUG
					std::cout << "The new seed " << curSeed << curSeedVec << curSeedMag << std::endl;
				#endif
				
//			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
//			CHECK_OPENGL_ERROR
//				findOptimalDirection();
//			Unbind();
//THIS IS EXPERIMENTAL
			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
				findOptimalScale();
			Unbind();
//THIS IS EXPERIMENTAL

//				cL.push_back(curSeed);
//				cM.push_back(curSeedMag);
//				cD.push_back(curSeedMag);
				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(std::vector<stim::vec3<float> > L, std::vector<float > M)
		{
			if(L.size() > 3)
			{
				sk.Begin(stim::OBJ_LINE);
				for(int i = 0; i < L.size(); i++)
				{
					sk.TexCoord(M[i]);
					sk.Vertex(L[i][0], L[i][1], L[i][2]);
				}
				sk.End();

				nt.addEdge(L,M);
				#ifdef DEBUG
					iter++;
				#endif
			}
		}


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

		void
		traceLine(stim::vec3<float> pos, float mag, int min_cost)
		{
			//starting (seed) position and magnitude.
			stim::vec3<float> spos = getPosition();
			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;
					branchDetection2();
					addToNetwork(cL, cM);
					#ifdef DEBUG
					std::cerr << "the cost of " << cost << " > " << min_cost << std::endl;
					#endif
					break;
				} else {
					//Have we found the edge of the map?
					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();
						addToNetwork(cL, cM);
						#ifdef DEBUG
						std::cerr << "I hit and edge" << std::endl;
						#endif
						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 > 75 || mag < 1){
						running = false;
						branchDetection2();
						addToNetwork(cL, cM);
						#ifdef DEBUG
						std::cerr << "The templates are too big" << std::endl;
						#endif
						break;
					}
					else
					{
						h = selectObject(p, getDirection(), m);
						//Have we hit something previously traced?
						if(h != -1){
							#ifdef DEBUG
							std::cerr << "I hit the fiber " << h << std::endl;
							#endif
							running = false;
							branchDetection2();
							addToNetwork(cL, cM);
							break;
						}
						else {  
							cL.push_back(stim::vec3<float>(p[0], p[1],p[2]));
							cM.push_back(m);
//							Unbind();
							CHECK_OPENGL_ERROR
						}
				 	}
                 		}
         		}
				#ifdef DEBUG
				std::cout << "I broke out!" << std::endl;
				#endif
		}
};
}
#endif