stim/gl/gl_spider.h

#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/rect.h"
#include "stim/math/matrix.h"
#include "stim/cuda/spider_cost.cuh"
#include <stim/cuda/cudatools/glbind.h>
#include <stim/cuda/arraymath.cuh>
#include <stim/cuda/cudatools.h>
#include <stim/cuda/ivote.cuh>
#include <stim/visualization/glObj.h>
#include <vector>
#include <stim/cuda/branch_detection.cuh>
#include "../../../volume-spider/fiber.h"
#include "../../../volume-spider/glnetwork.h"
#include <stim/visualization/cylinder.h>
//#include <stim/cuda/testKernel.cuh>
//#include <stim/cuda/testKernel.cuh>
#include <iostream>
#include <fstream>
#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:
		//
		stim::vec<float> p;  	//vector designating the position of the spider.
		stim::vec<float> d;	//vector designating the orientation of the spider
						//always a unit vector.
		stim::vec<float> m;	//magnitude of the spider vector.
						//mag[0] = length.
						//mag[1] = width.
		std::vector<stim::vec<float> > dV;	//A list of all the direction vectors.
		std::vector<stim::vec<float> > pV;	//A list of all the position vectors.
		std::vector<stim::vec<float> > mV;	//A list of all the size vectors.
		stim::matrix<float, 4> cT;		//current Transformation matrix
							//From tissue space to texture space.
		GLuint texID;
		stim::vec<float> S;			//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::vec<float> > 	seeds;		//seed positions.
		std::stack< stim::vec<float> > 	seedsvecs;	//seed directions.
		std::stack< float > 		seedsmags;	//seed magnitudes.
		std::vector< stim::vec<float> > cL;	//Positions of line currently being traced.
		std::vector< stim::vec<float> > cD;	//Direction of line currently being traced.
		std::vector< stim::vec<float> > cM;	//Magnitude of line currently being traced.
		stim::glnetwork<float> nt;		//object for storing the network.
		stim::vec<float> rev;			//reverse vector;
		stim::camera camSel;
		stim::vec<float> ps;
		stim::vec<float> ups;
		stim::vec<float> ds;
		
//--------------------------------------------------------------------------//
//-------------------------------PRIVATE METHODS----------------------------//
//--------------------------------------------------------------------------//
		/// Method for finding the best scale for the spider.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. 
		void
		findOptimalDirection()
		{
			setMatrix();			//create the transformation matrix.
			glCallList(dList);		//move the templates to p, d, m.
			int best = getCost(texbufferID,numSamples);		//find min cost.
			stim::vec<float> next(		//find next vector.
 			dV[best][0]*S[0]*R[0],
			dV[best][1]*S[1]*R[1],
			dV[best][2]*S[2]*R[2],
			0);
			next = (cT*next).norm();	//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()
		{
			setMatrix();			//create the transformation matrix.
			glCallList(dList+1);		//move the templates to p, d, m.
			int best = getCost(ptexbufferID, numSamplesPos);		//find min cost.
			std::cerr << best << std::endl;
			stim::vec<float> next(		//find next position.
 				pV[best][0],
				pV[best][1],
				pV[best][2],
				1);
			next = cT*next;			//find next position.
			setPosition(
					next[0]*S[0]*R[0],
					next[1]*S[1]*R[1],
					next[2]*S[2]*R[2]
				   );			//adjust position.
		}
	
		/// Method for finding the best scale for the spider.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. */
		void
		findOptimalScale()
		{
			setMatrix();			//create the transformation.
			glCallList(dList+2);		//move the templates to p, d, m.
			int best = getCost(mtexbufferID, numSamplesMag);		//get best cost.
			setMagnitude(m[0]*mV[best][0]);	//adjust the magnitude.
		}
		///subject to change.
		///finds branches.
		///depreciated
		void
		branchDetection()
		{
			setMatrix();
			glCallList(dList+3);
			std::vector< stim::vec<float> > result = find_branch(
					btexbufferID, GL_TEXTURE_2D, 16, 216);
			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
			if(!result.empty())
			{
				for(int i = 1; i < result.size(); i++)
				{
					stim::vec<float> cylp(
						0.5 * cos(2*M_PI*(result[i][1])),
						0.5 * sin(2*M_PI*(result[i][1])),
						result[i][0]-0.5,
						1.0);
					cylp = cT*cylp;
					
					stim::vec<float> vec(
						cylp[0]*S[0]*R[0],
						cylp[1]*S[1]*R[1],
						cylp[2]*S[2]*R[2]);
					stim::vec<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],
						   -p[1] + cylp[1]*S[1]*R[1],
						   -p[2] + cylp[2]*S[2]*R[2]);
						seeddir = seeddir.norm();
					float seedm = m[0];
// 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)
		{
			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::vec<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::vec<float> v = cyl.surf(pval, result[i][0]);
					stim::vec<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);
					}
				}
			}
			}
		}
		
//--------------------------------------------------------------------------//
//---------------------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 = 5/M_PI*4)
		{
			//Set up the vectors necessary for Rectangle creation.
			vec<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			vec<float> pos(0.0,0.0,0.0);
			vec<float> mag(1.0, 1.0, 1.0);
			vec<float> dir(0.0, 0.0, 1.0);
			//Set up the variable necessary for vector creation.
			vec<float> d_s = d.cart2sph().norm();
			vec<float> temp(0,0,0);
			int dim = (sqrt(numSamples)-1)/2;
			float p0  	= -M_PI;				 //phi angle in spherical coordinates.
			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0)); //step size in Theta.
			float dp  	= p0/(2.0*((float)dim + 1.0));		 //step size in Phi.
			
			glNewList(dList, GL_COMPILE); 		
			//Loop over the above defined space creating distinct vectors. 
			int idx = 0;
			for(int i = -dim; i <= dim; i++){
				for(int j = -dim; j <= dim; j++){
					
					//Create linear index
					idx = (j+dim)+(i+dim)*((dim*2)+1);	
					temp[0] = d_s[0]; 			//rotate vector
					temp[1] = d_s[1]+dp*(float) i;
					temp[2] = d_s[2]+dt*(float) j;
					
					temp = (temp.sph2cart()).norm();	//back to cart
					dV.push_back(temp);
				 	if(cos(Y.dot(temp))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
					else{Y[0] = 1.0; Y[1] = 0.0;}
                                       	hor = stim::rect<float>(mag,
						 pos, temp,
                       				((Y.cross(temp)).cross(temp)).norm());
       				 	ver = stim::rect<float>(mag,
						 pos, temp,
                       				hor.n());
					UpdateBuffer(0.0, 0.0+idx*8.0);
					CHECK_OPENGL_ERROR
				}
			}
			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.
			vec<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			vec<float> pos(0.0,0.0,0.0);
			vec<float> mag(1.0, 1.0, 1.0);
			vec<float> dir(0.0, 0.0, 1.0);
			//Set up the variable necessary for vector creation.
			vec<float> temp(0,0,0);
			int dim = (sqrt(numSamplesPos)-1)/2;	//number of position vectors.
			stim::rect<float> samplingPlane =	//plane from which we pull position samples
				 stim::rect<float>(p, d);
			samplingPlane.scale(mag[0]*delta, mag[0]*delta);
			float step = 1.0/(dim);			//step size.
			//Loop over the samples, keeping the original p samples in the center of the resulting texture to create a large number of position vectors.
			int idx;
			glNewList(dList+1, GL_COMPILE);
			for(int i = -dim; i <= dim; i++){
				for(int j = -dim; j <= dim; j++){
					//Create linear index
					idx = (j+dim)+(i+dim)*((dim*2)+1);	
					temp = samplingPlane.p(
							0.5+step*i,
								 0.5+step*j
										);
					pV.push_back(temp);
                			hor = stim::rect<float>(mag,
						 temp, dir,
                      				((Y.cross(d)).cross(d))
						.norm());
                			ver = stim::rect<float>(mag,
						 temp, dir,
                        			hor.n());
					UpdateBuffer(0.0, 0.0+idx*8.0);
				CHECK_OPENGL_ERROR
				}
			}
			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.
			vec<float> Y(1.0,0.0,0.0);		//orthogonal vec.
			vec<float> pos(0.0,0.0,0.0);
			vec<float> mag(1.0, 1.0, 1.0);
			vec<float> dir(0.0, 0.0, 1.0);
			//Set up the variable necessary for vector creation.
			int dim = (sqrt(numSamplesMag)-1)/2;
			float min 	= 1.0-delta;
			float max 	= 1.0+delta;
			float step	= (max-min)/(numSamplesMag-1);
			float factor;
			vec<float> temp(0.0,0.0,0.0);
			glNewList(dList+2, GL_COMPILE);
			for(int i = 0; i < 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*8.0);
			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)
		{	
			float len = 8.0;
			stim::vec<float>p1; 
        	        stim::vec<float>p2; 
	                stim::vec<float>p3; 
                	stim::vec<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+len, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+len, v_y+len);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x, v_y+len);
			 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+len, v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+2.0*len, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+2.0*len, v_y+len);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x+len, v_y+len);
			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, texels);   
			delete[] texels;
			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
			glBindTexture(GL_TEXTURE_2D, 0);
			CHECK_OPENGL_ERROR
		}
		///@param uint width sets the width of the buffer.
		///@param uint height sets the height of the buffer.
		///Function for setting up the 2D buffer that stores the samples.
		void
		GenerateFBO(unsigned int width, unsigned int height)
		{
			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, texels);   
			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(float len = 8.0)
		{
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
			glFramebufferTexture2D(
				GL_FRAMEBUFFER,
				GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_2D,
				texbufferID,
				0);
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
			glDrawBuffers(1, DrawBuffers);
			glBindTexture(GL_TEXTURE_2D, texbufferID);
			glClearColor(1,1,1,1);
			glClear(GL_COLOR_BUFFER_BIT);
			glMatrixMode(GL_PROJECTION);
			glLoadIdentity();
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glViewport(0,0,2.0*len, numSamples*len);
			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
			glEnable(GL_TEXTURE_3D);
			glBindTexture(GL_TEXTURE_3D, texID);
			CHECK_OPENGL_ERROR
		}
		
		///Method for controling the buffer and texture binding.
		///Clears the buffer upon binding.
		///@param GLuint &textureID, texture to be bound.
		///@param GLuint &framebufferID, framebuffer used for storage.
		///@param int nSamples, number of rectanges to create. 
		void
		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
		{
			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
			glFramebufferTexture2D(
				GL_FRAMEBUFFER,
				GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_2D,
				textureID,
				0);
			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
			glDrawBuffers(1, DrawBuffers);
			glBindTexture(GL_TEXTURE_2D, textureID);
//			glClearColor(1,1,1,1);
//			glClear(GL_COLOR_BUFFER_BIT);
			glMatrixMode(GL_PROJECTION);
			glLoadIdentity();
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glViewport(0,0,2.0*len, nSamples*len);
			gluOrtho2D(0.0,2.0*len,0.0,nSamples*len);
			glEnable(GL_TEXTURE_3D);
			glBindTexture(GL_TEXTURE_3D, texID);
			CHECK_OPENGL_ERROR
		}
		
		///Unbinds all texture resources.
		void
		Unbind()
		{
			//Finalize GL_buffer
			glBindTexture(GL_TEXTURE_3D, 0);                      
			CHECK_OPENGL_ERROR
			glBindTexture(GL_TEXTURE_2D, 0);                      
			CHECK_OPENGL_ERROR
			glBindFramebuffer(GL_FRAMEBUFFER, 0);
			CHECK_OPENGL_ERROR
			glDisable(GL_TEXTURE_3D);
			CHECK_OPENGL_ERROR
		}
		///Makes the spider take a step.
		///starting with the current p, d, m, find the next optimal p, d, m.
		///Performs the branch detection on each step.
		int
		StepP()
		{
			Bind();
			CHECK_OPENGL_ERROR
			#ifdef TESTING
				start = std::clock();
			#endif
			findOptimalDirection();
			findOptimalPosition();
			findOptimalScale();
			Unbind();
			Bind(btexbufferID, bfboID, 27);
			branchDetection();
			Unbind();
			#ifdef TESTING
				duration_sampling = duration_sampling +
					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
				num_sampling = num_sampling + 1.0; 
			#endif
			return current_cost;
		}
		
//--------------------------------------------------------------------------//
//--------------------------------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()
		{
			#ifdef TESTING
				start = std::clock();
			#endif
			stim::vec<int> cost = 
				stim::cuda::get_cost(texbufferID, GL_TEXTURE_2D, numSamples);
			cudaDeviceSynchronize();
			#ifdef TESTING
				duration_cuda = duration_cuda +
					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
				num_cuda = num_cuda + 1.0;
			#endif
			current_cost = cost[1];
			return cost[0];
		}
		int
		getCost(GLuint tID, int n)
		{
			#ifdef TESTING
				start = std::clock();
			#endif
			stim::vec<int> cost = 
				stim::cuda::get_cost(tID, GL_TEXTURE_2D, n);
			cudaDeviceSynchronize();
			#ifdef TESTING
				duration_cuda = duration_cuda +
					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
				num_cuda = num_cuda + 1.0;
			#endif
			current_cost = cost[1];
			return cost[0];
		}
	public:
		///ininializes the cuda device and environment.
		void
		initCuda()
		{	
			stim::cudaSetDevice();
			//GLint max;
			//glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
			//std::cout << max << std::endl;
		}
		//horizonal rectangle forming the spider.
		stim::rect<float> hor;
		//vectical rectangle forming the spider.
		stim::rect<float> ver;	
		//Testing and Timing variables.
		#ifdef TESTING
			std::clock_t start;
			double duration_sampling = 0.0;
			double duration_cuda = 0.0;
			double num_sampling = 0.0;
			double num_cuda = 0.0;
		#endif
//--------------------------------------------------------------------------//
//-----------------------------CONSTRUCTORS---------------------------------//
//--------------------------------------------------------------------------//
		///@param 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 = vec<float>(0.0, 0.0, 0.0);
			d = vec<float>(0.0, 0.0, 1.0);
			m = vec<float>(1.0, 1.0);
			S = vec<float>(1.0, 1.0, 1.0);
			R = vec<float>(1.0, 1.0, 1.0);
			numSamples = samples;
			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 = vec<float>(pos_x, pos_y, pos_z);
			d = vec<float>(dir_x, dir_y, dir_z);
			m = vec<float>(mag_x, mag_x, mag_x);
			S = vec<float>(1.0,1.0,1.0);
			R = vec<float>(1.0,1.0,1.0);
			numSamples = numsamples;
			numSamplesPos = numsamplespos;
			numSamplesMag = numsamplesmag;
		}
		///Position constructor: vecs of floats.
		///@param stim::vec<float> pos, position.
		///@param stim::vec<float> dir, direction.
		///@param float mag, size of the vector.
		///@param int 	samples, number of templates this spider is going to use.
		gl_spider
		(stim::vec<float> pos, stim::vec<float> dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144)
		{
			p = pos;
			d = dir;
			m = vec<float>(mag, mag, mag);
			S = vec<float>(1.0,1.0,1.0);
			R = vec<float>(1.0,1.0,1.0);
			numSamples = samples;
			numSamplesPos = samplesPos;
			numSamplesMag = samplesMag;
		}
		///destructor	
		~gl_spider
		(void)
		{
			Unbind();
			glDeleteTextures(1, &texbufferID);
			glDeleteBuffers(1, &fboID);
			glDeleteTextures(1, &ptexbufferID);
			glDeleteBuffers(1, &pfboID);
			glDeleteTextures(1, &mtexbufferID);
			glDeleteBuffers(1, &mfboID);     
			glDeleteTextures(1, &btexbufferID);
			glDeleteBuffers(1, &bfboID);
		}
		///@param GLuint id, texture that is going to be sampled.
		///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)
		{
			texID = id;
			//GenerateFBO(16, numSamples*8);
			GenerateFBO(16, numSamples*8, texbufferID, fboID); 
			GenerateFBO(16, numSamplesPos*8, ptexbufferID, pfboID); 
			GenerateFBO(16, numSamplesMag*8, mtexbufferID, mfboID); 
			GenerateFBO(16, 216, btexbufferID, bfboID); 
			setDims(0.6, 0.6, 1.0);
			setSize(512.0, 512.0, 426.0);
			setMatrix();
			dList = glGenLists(3);
			glListBase(dList);
			Bind(texbufferID, fboID, numSamples);
				genDirectionVectors(5*M_PI/4);
			Unbind();
			Bind(ptexbufferID, pfboID, numSamplesPos);
				genPositionVectors();
			Unbind();
			Bind(mtexbufferID, mfboID, numSamplesMag);
				genMagnitudeVectors();
			Unbind();
			Bind(btexbufferID, bfboID, 27);
			DrawCylinder();
			Unbind();
		}
		
//--------------------------------------------------------------------------//
//-----------------------------ACCESS METHODS-------------------------------//
//--------------------------------------------------------------------------//
		///Returns the p vector.
		vec<float>
		getPosition()
		{
			return p;
		}
	
		///Returns the d vector.
		vec<float>
		getDirection()
		{
			return d;
		}
		///Returns the m vector.
		vec<float>
		getMagnitude()
		{
			return m;
		}
	
		///@param stim::vec<float> pos, the new p.
		///Sets the p vector to input vector pos.
		void
		setPosition(vec<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(vec<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(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::vec<float> Dims)
		{
			S = Dims;
		}
		///@param float x, size of the data in the x direction.
		///@param float y, size of the data in the y direction.
		///@param float z, size of the data in the z direction.
		///Sets the data volume sizes in each direction.
		void
		setSize(float x, float y, float z)
		{
			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::vec<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::vec<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);;
				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::vec<float> pos)
		{
			seeds.push(pos);
		}
		///@param stim::vec<float> dir, the direction of the seed to be added.
		///Adds a seed to the seed directions list.
		void
		setSeedVec(stim::vec<float> dir)
		{
			seedsvecs.push(dir);
		}
		///@param float mag, the size of the seed to be added.
		///Adds a seed to the seed list.
		///Assumes that the coordinates passes are in tissue space.
		void
		setSeedMag(float mag)
		{
			seedsmags.push(mag);
		}
		///@param float x, x-position of the seed to be added.
		///@param float y, y-position of the seed to be added.
		///@param float z, z-position of the seed to be added.
		///Adds a seed to the seed list.
		///Assumes that the coordinates passes are in tissue space.
		void
		setSeed(float x, float y, float z)
		{
			seeds.push(stim::vec<float>(x, y, z));
		}
		///@param float x, x-direction of the seed to be added.
		///@param float y, y-direction of the seed to be added.
		///@param float z, z-direction of the seed to be added.
		///Adds a seed to the seed directions list.
		void
		setSeedVec(float x, float y, float z)
		{
			seedsvecs.push(stim::vec<float>(x, y, z));
		}
		
		///Method to get the top of the seed positions stack.
		stim::vec<float> 
		getLastSeed()
		{
			stim::vec<float> tp = seeds.top();
			return tp;
		}
		///Method to get the top of the seed direction stack.
		stim::vec<float> 
		getLastSeedVec()
		{
			stim::vec<float> tp = seedsvecs.top();
			return tp;
		}
		///Method to get the top of the seed magnitude stack.
		float
		getLastSeedMag()
		{
			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::vec<float> >
		getSeeds()
		{
			return seeds;
		}
		///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::vec< 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 bfboID;
		}
//--------------------------------------------------------------------------//
//-----------------------------TEMPORARY METHODS----------------------------//
//--------------------------------------------------------------------------//
		///temporary Method necessary for visualization and testing.
		void
		Update()
		{
			vec<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);
			CHECK_OPENGL_ERROR
			#ifdef TESTING
				start = std::clock();
			#endif
				findOptimalDirection();
			Unbind();
			Bind(ptexbufferID, pfboID, numSamplesPos);
				findOptimalPosition();
			Unbind();
			Bind(mtexbufferID, mfboID, numSamplesMag);
				findOptimalScale();
			Unbind();
			CHECK_OPENGL_ERROR
			#ifdef TESTING
				duration_sampling = duration_sampling +
					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
				num_sampling = num_sampling + 1.0; 
			#endif
			return current_cost;
		}
		void
		printTransform()
		{
			std::cout << cT << std::endl;
		}
//--------------------------------------------------------------------------//
//-----------------------------EXPERIMENTAL METHODS-------------------------//
//--------------------------------------------------------------------------//
		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::vec<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)
		{	
//			rev = stim::vec<float>(0.0,0.0,1.0);
			bool sEmpty = true;
			float lastmag = 16.0;;
			stim::vec<float> curSeed; 
			stim::vec<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::vec<float> loc, stim::vec<float> dir, float mag) 
		{
		//Define the varibles and turn on Selection Mode
			float s = 3.0;
			GLuint selectBuf[2048];
			GLint hits;
			glSelectBuffer(2048, selectBuf);
			glDisable(GL_CULL_FACE);
			(void) glRenderMode(GL_SELECT);
		//Init Names stack
			glInitNames();
			glPushName(1);
			CHECK_OPENGL_ERROR
		//What would that vessel see in front of it.
				camSel.setPosition(loc);
				camSel.setFocalDistance(mag/s);
				camSel.LookAt((loc[0]+dir[0]*mag/s),
					 (loc[1]+dir[1]*mag/s),
					 (loc[2]+dir[2]*mag/s));
				ps  = camSel.getPosition();
				ups = camSel.getUp();
				ds  = camSel.getLookAt();
				glMatrixMode(GL_PROJECTION);
				glPushMatrix();
				glLoadIdentity();
				glOrtho(-mag/s/2.0, mag/s/2.0, -mag/s/2.0, mag/s/2.0, 0.0, mag/s/2.0);
				glMatrixMode(GL_MODELVIEW);
				glPushMatrix();
				glLoadIdentity();
				
				CHECK_OPENGL_ERROR
				gluLookAt(ps[0], ps[1], ps[2],
					 ds[0], ds[1], ds[2],
					 ups[0], ups[1], ups[2]);
//				sk.Render();
				nt.Render();								
				CHECK_OPENGL_ERROR
//				glLoadName((int) sk.numL());
				glLoadName(nt.sizeE());
//				sk.RenderLine(cL);
				nt.RenderLine(cL);	
//				glPopName();
				glFlush();
				glMatrixMode(GL_PROJECTION);
				glPopMatrix();
				glMatrixMode(GL_MODELVIEW);
				CHECK_OPENGL_ERROR
				glPopMatrix();
		//	glEnable(GL_CULL_FACE);
			hits = glRenderMode(GL_RENDER);
			int found_hits = processHits(hits, selectBuf);
			return found_hits;
		}
		//Given a size of the array (hits) and the memory holding it (buffer)
		//returns whether a hit tool place or not.
		int
		processHits(GLint hits, GLuint buffer[])
		{
			GLuint names, *ptr;
			//printf("hits = %u\n", hits);
			ptr = (GLuint *) buffer;
		//	for (int i = 0; i < hits; i++) { /*  for each hit  */
			names = *ptr;
		//		printf (" number of names for hit = %u\n", names);
			ptr++;
			ptr++; //Skip the minimum depth value.  
			ptr++; //Skip the maximum depth value.
		//		printf ("   the name is ");                                                        
		//		for (int j = 0; j < names; j++) {     /*  for each name */
		//			printf ("%u ", *ptr); ptr++;
		//		}
		//		printf ("\n");
	//		}
			
			if(hits == 0)
			{
				return -1;
			}
			else
			{
//				printf ("%u ", *ptr);
				return *ptr;
			}
		}
		void
		clearCurrent()
		{
			cL.clear();
			cM.clear();
		}	
	
		void
		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec<float> spos, 
				stim::vec<float> smag, stim::vec<float> sdir)
		{
                        std::vector<stim::vec<float> > ce = in.first.centerline();                
                        std::vector<stim::vec<float> > cm = in.first.centerlinemag();
			//if the fiber is longer than 2 steps (the number it takes to diverge)
			if(ce.size() > 3)
			{	
				//if we did not hit a fiber
				if(in.second == -1)
				{
					spos[0] = spos[0]-sdir[0]*smag[0]/2.;
					spos[1] = spos[1]-sdir[1]*smag[0]/2.;
					spos[2] = spos[2]-sdir[2]*smag[0]/2.;
					int h = selectObject(spos, -sdir, smag[0]);
					//did we start with a fiber?
					if(h != -1)
						nt.addEdge(ce, cm, h, -1);
					else
						nt.addEdge(ce, cm, -1, -1);
				}
				//if we hit a fiber?
				else if(in.second != -1)
				{
					nt.addEdge(ce,cm,-1, in.second);
					spos[0] = spos[0]-sdir[0]*smag[0]/2.;
					spos[1] = spos[1]-sdir[1]*smag[0]/2.;
					spos[2] = spos[2]-sdir[2]*smag[0]/2.;
					int h = selectObject(spos, -sdir, smag[0]);
					//did start with a fiber?
					if(h != -1 && 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);}
				}
			}		
		}
		void
		printSizes()
		{
			std::cout << nt.sizeE() << " edges " << std::endl;
			std::cout << nt.sizeV() << " nodes " << std::endl;
			
		}
		std::pair<stim::fiber<float>, int >
		traceLine(stim::vec<float> pos, stim::vec<float> mag, int min_cost)
		{
			//starting (seed) position and magnitude.
			stim::vec<float> spos = getPosition();
			stim::vec<float> smag = getMagnitude();
			stim::vec<float> sdir = getDirection();	
			Bind();
//			sk.Begin(stim::OBJ_LINE);
//			sk.createFromSelf(GL_SELECT);
			nt.createFromSelf(GL_SELECT);
			cL.push_back(pos);
			cM.push_back(mag);
//			setPosition(pos);
//			setMagnitude(mag);
			int h;
			bool started = false;
			bool running = true;
			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
			while(running)
			{
				int cost = Step();
				if (cost > min_cost){
					running = false;
//					sk.End();
					branchDetection2();
					pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
					addToNetwork(a, spos, smag, sdir);
					return a;
					break;
				} else {
					//Have we found the edge of the map?
					pos = getPosition();
					if(pos[0] > size[0] || pos[1] > size[1]
					 || pos[2] > size[2] || pos[0] < 0
					 || pos[1] < 0 || pos[2] < 0)
					{
						running = false;
						branchDetection2();
						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
						addToNetwork(a, spos, smag, sdir);
						return a;
						break;
					}
					//If this is the first step in the trace,
					// save the direction
					//(to be used later to trace the fiber in the opposite direction)
					if(started == false){
						rev = -getDirection();
						started = true;
					}
					//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);
						return a;
						break;
					}
					else
					{
						h = selectObject(p, getDirection(), m[0]);
						//Have we hit something previously traced?
						if(h != -1){
							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::vec<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