#ifndef STIM_GL_SPIDER_H #define STIM_GL_SPIDER_H #include #include #include #include #include #include #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 #include #include #include #include #include #include #include "../../../volume-spider/fiber.h" #include "../../../volume-spider/glnetwork.h" #include //#include //#include #include #include #ifdef TESTING #include #include #include #endif namespace stim { template class gl_spider : public virtual gl_texture { //doen't use gl_texture really, just needs the GLuint id. //doesn't even need the texture iD really. private: // stim::vec p; //vector designating the position of the spider. stim::vec d; //vector designating the orientation of the spider //always a unit vector. stim::vec m; //magnitude of the spider vector. //mag[0] = length. //mag[1] = width. std::vector > dV; //A list of all the direction vectors. std::vector > pV; //A list of all the position vectors. std::vector > mV; //A list of all the size vectors. stim::matrix cT; //current Transformation matrix //From tissue space to texture space. GLuint texID; stim::vec S; //Size of a voxel in the volume. stim::vec 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 > seeds; //seed positions. std::stack< stim::vec > seedsvecs; //seed directions. std::stack< float > seedsmags; //seed magnitudes. std::vector< stim::vec > cL; //Positions of line currently being traced. std::vector< stim::vec > cD; //Direction of line currently being traced. std::vector< stim::vec > cM; //Magnitude of line currently being traced. stim::glnetwork nt; //object for storing the network. stim::vec rev; //reverse vector; stim::camera camSel; stim::vec ps; stim::vec ups; stim::vec 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 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 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 > result = find_branch( btexbufferID, GL_TEXTURE_2D, 16, 216); stim::vec 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 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 vec( cylp[0]*S[0]*R[0], cylp[1]*S[1]*R[1], cylp[2]*S[2]*R[2]); stim::vec 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 cyl(cL, cM); std::vector< stim::vec > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template); stim::vec 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 v = cyl.surf(pval, result[i][0]); stim::vec 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 Y(1.0,0.0,0.0); //orthogonal vec. vec pos(0.0,0.0,0.0); vec mag(1.0, 1.0, 1.0); vec dir(0.0, 0.0, 1.0); //Set up the variable necessary for vector creation. vec d_s = d.cart2sph().norm(); vec 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(mag, pos, temp, ((Y.cross(temp)).cross(temp)).norm()); ver = stim::rect(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 Y(1.0,0.0,0.0); //orthogonal vec. vec pos(0.0,0.0,0.0); vec mag(1.0, 1.0, 1.0); vec dir(0.0, 0.0, 1.0); //Set up the variable necessary for vector creation. vec temp(0,0,0); int dim = (sqrt(numSamplesPos)-1)/2; //number of position vectors. stim::rect samplingPlane = //plane from which we pull position samples stim::rect(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(mag, temp, dir, ((Y.cross(d)).cross(d)) .norm()); ver = stim::rect(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 Y(1.0,0.0,0.0); //orthogonal vec. vec pos(0.0,0.0,0.0); vec mag(1.0, 1.0, 1.0); vec 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 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(temp, pos, dir, ((Y.cross(d)).cross(d)) .norm()); ver = stim::rect(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::vecp1; stim::vecp2; stim::vecp3; stim::vecp4; 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 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 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 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 hor; //vectical rectangle forming the spider. stim::rect 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(0.0, 0.0, 0.0); d = vec(0.0, 0.0, 1.0); m = vec(1.0, 1.0); S = vec(1.0, 1.0, 1.0); R = vec(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(pos_x, pos_y, pos_z); d = vec(dir_x, dir_y, dir_z); m = vec(mag_x, mag_x, mag_x); S = vec(1.0,1.0,1.0); R = vec(1.0,1.0,1.0); numSamples = numsamples; numSamplesPos = numsamplespos; numSamplesMag = numsamplesmag; } ///Position constructor: vecs of floats. ///@param stim::vec pos, position. ///@param stim::vec 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 pos, stim::vec dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144) { p = pos; d = dir; m = vec(mag, mag, mag); S = vec(1.0,1.0,1.0); R = vec(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 getPosition() { return p; } ///Returns the d vector. vec getDirection() { return d; } ///Returns the m vector. vec getMagnitude() { return m; } ///@param stim::vec pos, the new p. ///Sets the p vector to input vector pos. void setPosition(vec 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 dir, the new d. ///Sets the d vector to input vector dir. void setDirection(vec dir) { d = dir; } ///@param stim::vec x x-coordinate. ///@param stim::vec y y-coordinate. ///@param stim::vec 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 dir, the new d. ///Sets the m vector to the input vector mag. void setMagnitude(vec 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 Dims, voxel size. ///Sets the voxel sizes in each direction. necessary for non-standard data. void setDims(stim::vec 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 Dims, size of the volume. ///Sets the data volume sizes in each direction. void setSize(stim::vec Siz) { S = Siz; } ///@param stim::vec 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 getRotation(stim::vec dir) { stim::vec out(0.0,0.0,0.0,0.0); stim::vec 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 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 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 pos) { seeds.push(pos); } ///@param stim::vec dir, the direction of the seed to be added. ///Adds a seed to the seed directions list. void setSeedVec(stim::vec 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(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(x, y, z)); } ///Method to get the top of the seed positions stack. stim::vec getLastSeed() { stim::vec tp = seeds.top(); return tp; } ///Method to get the top of the seed direction stack. stim::vec getLastSeedVec() { stim::vec 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 > 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 sk; for(int i = 0; i < nt.sizeE(); i++) { std::vector > cm = nt.getEdgeCenterLineMag(i); std::vector > 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 getNetwork() { // return sk; } ///returns a COPY of the entire stim::glnetwork object. stim::glnetwork 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 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(m, p, d.norm(), ((Y.cross(d)).cross(d)).norm()); ver = stim::rect(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 cyl(cL, cM); std::vector > > 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(0.0,0.0,1.0); bool sEmpty = true; float lastmag = 16.0;; stim::vec curSeed; stim::vec 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, int> a = traceLine(p, m, min_cost); } } int selectObject(stim::vec loc, stim::vec 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, int> in, stim::vec spos, stim::vec smag, stim::vec sdir) { std::vector > ce = in.first.centerline(); std::vector > 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, int > traceLine(stim::vec pos, stim::vec mag, int min_cost) { //starting (seed) position and magnitude. stim::vec spos = getPosition(); stim::vec smag = getMagnitude(); stim::vec 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 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, int> a(stim::fiber (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, int> a(stim::fiber (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, int> a(stim::fiber (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, int> a(stim::fiber (cL, cM), h); addToNetwork(a, spos, smag, sdir); return a; break; } else { cL.push_back(stim::vec(p[0], p[1],p[2])); cM.push_back(stim::vec(m[0], m[0])); // Bind(btexbufferID, bfboID, 27); Unbind(); CHECK_OPENGL_ERROR } } } } } }; } #endif