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/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/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;
  		GLuint ptexbufferID;
  
  		GLuint mfboID;
  		GLuint mtexbufferID;
  

  		GLuint bfboId;
  		GLuint btexbufferID;
  

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

  		float stepsize = 4.0;			//Step size.

  		int current_cost;

  		
  		
  		//Tracing variables.
  		std::stack< stim::vec<float> > seeds;	//Variables for tracing
  		std::stack< stim::vec<float> > seedsvecs;

  		std::stack< float > seedsmags;

  		std::vector< stim::vec<float> > cL;	//Line currently being traced.

  		stim::glObj<float> sk;

  		stim::glnetwork<float> nt;
  

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

  		std::vector<stim::vec<float> > last3;

  
  		/// 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();

  			glCallList(dList);

  			int best = getCost();

  			stim::vec<float> next(
   			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();

  			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]);

  		}
  

  		/// Method for finding the best d for the spider.
  		/// Not sure if necessary since the next p for the spider
  		/// will be at d * m.

  		void
  		findOptimalPosition()
  		{

  			setMatrix();

  			glCallList(dList+1);

  			int best = getCost();

  			stim::vec<float> next(
   			pV[best][0],
  			pV[best][1],
  			pV[best][2],
  			1);

  			next = cT*next;	

  			setPosition(
  					next[0]*S[0]*R[0],

  					next[1]*S[1]*R[1],

  					next[2]*S[2]*R[2]
  				   );

  		}
  	

  		/// 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();

  			glCallList(dList+2);

  			int best = getCost();

  			setMagnitude(m[0]*mV[best][0]);

  		}
  

  		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]/2.0;

  // Uncomment for global run 
  /*					stim::vec<float> lSeed = getLastSeed();

  					if(sqrt(pow((lSeed[0] - vec[0]),2)
  					 + pow((lSeed[1] - vec[1]),2) + 
  					 pow((lSeed[2] - vec[2]),2)) > m[0]/4.0

  					 &&  */

  					if(

  					 !(vec[0] > size[0] || vec[1] > size[1]
  					 || vec[2] > size[2] || vec[0] < 0
  					 || vec[1] < 0 || vec[2] < 0))
  					{
  						setSeed(vec);
  						setSeedVec(seeddir);

  						setSeedMag(seedm);

  					}

  				}
  			}    

  			
  		}
  

  		void

  		Optimize()

  		{

  			/*find the optimum d and scale */ 

  		}

  		
  		
  //--------------------------------------------------------------------------//
  //---------------------TEMPLATE CREATION METHODS----------------------------//
  //--------------------------------------------------------------------------//
  
  		///@param solidAngle, the size of the arc to sample.
  		///Method for populating the vector arrays with sampled vectors.

  		///uses the default d vector <0,0,1>

  		void

  		genDirectionVectors(float solidAngle = 5/M_PI*4)

  		{

  			//ofstream file;
  			//file.open("dvectors.txt");

  			//Set up the vectors necessary for Rectangle creation.

  			vec<float> Y(1.0,0.0,0.0);
  			vec<float> 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;
  			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0));
  			float dp  	= p0/(2.0*((float)dim + 1.0));

  			
  			glNewList(dList, GL_COMPILE);
  			//Loop over the space
  			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 solidAngle, the size of the arc to sample.
  		///Method for populating the buffer with the sampled texture.
  		///uses the default vector <0,0,0>

  		void

  		genPositionVectors(float delta = 0.4)

  		{

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

  			vec<float> temp(0,0,0);

  			int dim = (sqrt(numSamples)-1)/2;
  			stim::rect<float> samplingPlane =

  				 stim::rect<float>(p, d);

  			samplingPlane.scale(mag[0]*delta, mag[0]*delta);

  			float step = 1.0/(dim);

  			//Loop over the samples, keeping the original p sample
  			//in the center of the resulting texture.
  			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 solidAngle, the size of the arc to sample.
  		///Method for populating the buffer with the sampled texture.

  		///uses the default m <1,1,0>

  		void

  		genMagnitudeVectors(float delta = 0.70)

  //		genMagnitudeVectors(float delta = 0.50)

  		{
  			

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

  			int dim = (sqrt(numSamples)-1)/2;
  			float min 	= 1.0-delta;
  			float max 	= 1.0+delta;
  			float step	= (max-min)/(numSamples-1);
  			float factor;

  			vec<float> temp(0.0,0.0,0.0);

  			glNewList(dList+2, GL_COMPILE);

  			for(int i = 0; i < numSamples; 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 v_x x-coordinate in buffer-space,
  		///@param v_y y-coordinate in buffer-space.

  		///Samples the texturespace and places a sample in the provided coordinates
  		///of bufferspace.

  		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 width sets the width of the buffer.
  		///@param 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, GLuint &textureID, GLuint &framebufferID)
  		{
  			glGenFramebuffers(1, &framebufferID);
  			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
  			int numChannels = 1;
  			unsigned char* texels = new unsigned char[width * height * numChannels];
  			glGenTextures(1, &textureID);
  			glBindTexture(GL_TEXTURE_2D, textureID);
  			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 width sets the width of the buffer.
  		///@param 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);
  			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

  		}

  		///Method for using the gl manipulation to alighn 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).
  		void setMatrix()

  		{

  			float curTrans[16];

  			stim::vec<float> rot = getRotation(d);

  			glMatrixMode(GL_TEXTURE);
  			glLoadIdentity();

  			glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);

  			glTranslatef(p[0],
  				     p[1],
  				     p[2]);

  			glRotatef(rot[0], rot[1], rot[2], rot[3]);

  			glScalef(m[0],

  				 m[0],

  				 m[0]);

  			glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
  			cT.set(curTrans);

  		//	printTransform();

  			CHECK_OPENGL_ERROR
  			glMatrixMode(GL_MODELVIEW);

  		}

  
  		///Method for controling the buffer and texture binding in order to properly
  		///do the render to texture.
  		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 in order to properly
  		///do the render to texture.
  		///@param GLuint tbID
  		void
  		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples)
  		{
  			float len = 8.0;
  			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
  			CHECK_OPENGL_ERROR
  
  			glFramebufferTexture2D(
  				GL_FRAMEBUFFER,
  				GL_COLOR_ATTACHMENT0,
  				GL_TEXTURE_2D,
  				textureID,
  				0);
  			CHECK_OPENGL_ERROR
  
  			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
  			CHECK_OPENGL_ERROR
  
  			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
  			glDrawBuffers(1, DrawBuffers);
  			CHECK_OPENGL_ERROR
  
  			glBindTexture(GL_TEXTURE_2D, textureID);
  			CHECK_OPENGL_ERROR
  
  			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
  		}

  		///Method for Unbinding all of the 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
  		}

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

  		///Entry-point into the cuda code for calculating the cost
  		///	of a given samples array (in texture form) 

  		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];

  		}