safety commit, I broke something and couldn't go back because I didn't commit

Pavel Govyadinov
1 parent d06d3dbd
Showing 3 changed files with 129 additions and 45 deletions Show diff stats
stim/cuda/branch_detection.cuh
stim/cuda/testKernel.cuh
stim/gl/gl_spider.h
@@ -111,6 +111,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
+	test(texbufferID, texType, x, y);
 	std::vector<stim::vec<float> >  output;
 	initCuda(bytes_table, bytes_ds); 
@@ -123,6 +124,8 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	stim::cuda::tex_gaussian_blur2<float>(
 		gpuI, sigma, x, y, t.getTexture(), t.getArray()
 		);
+	stim::gpu2image<float>(gpuI, "Blur.jpg", x,y , 0, 255);
+//	stim::gpu2image<float>(t.getArray(), "ORIGINAL.jpg", x,y , 0, 255);
 	cudaDeviceSynchronize();
@@ -146,12 +149,14 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	cudaDeviceSynchronize();
 	stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, final_t, conn, x, y);
+	stim::gpu2image<float>(gpuCenters, "Centers.jpg", x, y, 0, 1);
 	cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
+	stim::cpu2image<float>(cpuCenters, "CentersXPU.jpg", x, y, 0, 1);
 	for(int i = 0; i < pixels; i++)
 	{
 		int ix = (i % x);
 		int iy = (i / x);
-		if((cpuCenters[i] == 1) && (ix > 4) && (ix < x-4))
+		if((cpuCenters[i] != 0))
 		{
 			float x_v = (float) ix;
@@ -14,9 +14,9 @@
 	///Initialization function, allocates the memory and passes the necessary
 	///handles from OpenGL and Cuda.
 	///@param DIM_Y			--integer controlling how much memory to allocate.
-	void initArray()
+	void initArray(int x, int y)
 	{
-			cudaMalloc( (void**) &print, 216*16*sizeof(float));     ///temporary
+			cudaMalloc( (void**) &print, x*y*sizeof(float));     ///temporary
 	}
 	///Deinit function that frees the memery used and releases the texture resource
@@ -48,12 +48,14 @@
 	{       
 		int x   = threadIdx.x + blockIdx.x * blockDim.x;
 		int y   = threadIdx.y + blockIdx.y * blockDim.y;
-		int idx = y*16+x;
+		int idx = y*64+x;
+//		int idx = y*32+x;
+//		int idx = y*16+x;
 		float valIn             = tex2D<unsigned char>(texIn, x, y);
 		float templa		= templ(x);
-		//print[idx]              = abs(valIn);             ///temporary
-		print[idx]              = abs(templa);             ///temporary
+		print[idx]              = valIn;             ///temporary
+		//print[idx]              = abs(templa);             ///temporary
 	}
@@ -64,7 +66,7 @@
 	///@param GLenum texType	--either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D
 	///				  may work with other gl texture types, but untested.
 	///@param DIM_Y, the number of samples in the template.
-	void test(GLint texbufferID, GLenum texType)
+	void test(GLint texbufferID, GLenum texType, int x, int y)
 	{
 		//Bind the Texture in GL and allow access to cuda.
@@ -72,16 +74,12 @@
 		//initialize the return arrays.
-		initArray();
-		
-		int x = 16;
-		int y = 27*8;
-		y = 8* 1089; 
+		initArray(x,y);
+		dim3 numBlocks(1, y);
+		dim3 threadsPerBlock(x, 1);
 		int max_threads = stim::maxThreadsPerBlock();
 		//dim3 threads(max_threads, 1);
 		//dim3 blocks(x / threads.x + 1, y);	
-		dim3 numBlocks(1, 1089);
-		dim3 threadsPerBlock(16, 8);
 		//dim3 numBlocks(2, 2);
 		//dim3 threadsPerBlock(8, 108);
@@ -92,7 +90,7 @@
 		cudaDeviceSynchronize();
 		stringstream name;      //for debugging
 		name << "FromTex.bmp";
-		stim::gpu2image<float>(print, name.str(),16,1089*8,0,1.0);
+		stim::gpu2image<float>(print, name.str(),x,y,0,255);
 		tx.UnmapCudaTexture();
 		cleanUP();
@@ -23,6 +23,7 @@
 #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>
@@ -218,6 +219,37 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
+
+		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);
+			std::cerr << "the number of points is " << result.size() << std::endl;
+			if(!result.empty())
+			{
+				for(int i = 0; i < result.size(); i++)
+				{
+					std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1] << std::endl;
+					stim::vec<float> v = cyl.surf(result[i][0], result[i][1]*360);
+					std::cout <<  v[0] << " ," << v[1] << " ," << v[2] << std::endl;
+					stim::vec<float> di = cyl.p(result[i][0]);
+					std::cout <<  di[0] <<  " ," << di[1] << " ," << di[2] << std::endl;
+					float rad = cyl.r(result[i][0]);
+					std::cout << rad << std::endl;
+					setSeed(v);
+					setSeedVec(stim::vec<float>(-di[0]+v[0], -di[1]+v[1], -di[2]+v[2]).norm());
+					setSeedMag(cyl.r(result[i][0]));
+				}
+			}
+			std::cout << "I ran the new branch detection" << std::endl;
+			}
+		}
+
 //--------------------------------------------------------------------------//
 //---------------------TEMPLATE CREATION METHODS----------------------------//
@@ -459,6 +491,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		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;
@@ -506,45 +539,60 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
+		///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()
+		void setMatrix(int 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();
+			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);
+				//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()
+		Bind(float len = 8.0)
 		{
-			float len = 8.0;
 			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
 			glFramebufferTexture2D(
 				GL_FRAMEBUFFER,
@@ -576,9 +624,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///@param GLuint &framebufferID, framebuffer used for storage.
 		///@param int nSamples, number of rectanges to create. 
 		void
-		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples)
+		Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
 		{
-			float len = 8.0;
 			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
 			glFramebufferTexture2D(
 				GL_FRAMEBUFFER,
@@ -1216,6 +1263,36 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		      	 glEnd();  
 			 glEndList();
 		}
+
+		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
@@ -1515,6 +1592,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				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;
@@ -1529,6 +1607,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 //					       std::cout << "Found Edge" << std::endl;
 						running = false;
 						sk.End();
+		branchDetection2();
 						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
 						addToNetwork(a, spos, smag, sdir);
 						return a;
@@ -1548,6 +1627,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 //						std::cout << "Magnitude Limit" << std::endl;
 						running = false;
 						sk.End();
+		branchDetection2();
 						pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
 						addToNetwork(a, spos, smag, sdir);
 						return a;
@@ -1560,6 +1640,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						if(h != -1){
 							running = false;
 							sk.End();
+		branchDetection2();
 							pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
 							addToNetwork(a, spos, smag, sdir);
 							return a;
@@ -1574,7 +1655,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 							sk.Vertex(p[0], p[1], p[2]);
 							Bind(btexbufferID, bfboID, 27);
 							CHECK_OPENGL_ERROR
-							branchDetection();
+//							branchDetection();
 							CHECK_OPENGL_ERROR
 							Unbind();
 							CHECK_OPENGL_ERROR