fixed a stask bug, added support for branch detection and minimized the number o…

…f call to branch detection as well as well as branching points added in the border regions (where the branch outpoint is outside of the search space)

fixed a stask bug, added support for branch detection and minimized the number o…
…f call to branch detection as well as well as branching points added in the border regions (where the branch outpoint is outside of the search space)
Pavel Govyadinov
1 parent 84eff8b1
Showing 10 changed files with 132 additions and 74 deletions Show diff stats
stim/cuda/branch_detection.cuh
stim/cuda/ivote/update_dir.cuh
stim/cuda/ivote/vote.cuh
stim/cuda/sharedmem.cuh
stim/cuda/spider_cost.cuh
stim/cuda/templates/conv2.cuh
stim/cuda/templates/conv2sep.cuh
stim/cuda/templates/gaussian_blur.cuh
stim/cuda/testKernel.cuh
stim/gl/gl_spider.h
@@ -117,7 +117,6 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	atan_2d(cpuTable, rmax);
 	cudaMemcpy(gpuTable, cpuTable, bytes_table, cudaMemcpyHostToDevice);
-	test(texbufferID, texType);
 	t.MapCudaTexture(texbufferID, texType);
 	cudaDeviceSynchronize();
@@ -125,55 +124,41 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 		gpuI, sigma, x, y, t.getTexture(), t.getArray()
 		);
 	cudaDeviceSynchronize();
-	stim::gpu2image<float>(gpuI, "Blur.jpg", 16, 8*27, 0, 255);	
 	stim::cuda::gpu_gradient_2d<float>(
 		gpuGrad, gpuI, x, y
 		);
 	cudaDeviceSynchronize();
-//	stim::gpu2image<float>(gpuGrad, "Grad.jpg", 16, 8*27, 0, 1);	
-	stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y, M_PI);
+	stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y);
 	cudaDeviceSynchronize();
-//	stim::gpu2image<float>(gpuGrad, "Cart_2_polar.jpg", 16, 8*27, 0, 1);	
 	cudaDeviceSynchronize();
 	for (int i = 0; i < iter; i++)
 	{
 		stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
-		name << "Vote" << i << ".bmp";
-		stim::gpu2image<float>(gpuVote, name.str(), 16, 8*27, 0, 255);
-		name.str("");
-		name.clear();	
 	cudaDeviceSynchronize();
-//		std::cout << "got here7_a_" << i << std::endl;
 		stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);
 	cudaDeviceSynchronize();
-//		std::cout << "got here7_b_" << i << std::endl;
 		phi = phi - dphi;
 	}
 	cudaDeviceSynchronize();
 	stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, final_t, conn, x, y);
-//	std::cout << "got here_sdkfj" << std::endl;
 	cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
-	stim::gpu2image<float>(gpuCenters, "Centers.jpg", 16, 8*27, 0, 255);	
 	for(int i = 0; i < pixels; i++)
 	{
 		int ix = (i % x);
 		int iy = (i / x);
-//		std::cout << i << " : " << ix <<" : " <<  iy << std::endl;
-		if((cpuCenters[i] == 1) && (ix > 4) && (ix < x-4))
+		if((cpuCenters[i] == 1) && (ix > 5) && (ix < x-5))
 		{
-			std::cout << ix << " : " << iy << std::endl; 
 			float x_v = (float) ix;
 			float y_v = (float) iy;
 			output.push_back(stim::vec<float>((x_v/(float)x),
 							  (y_v/(float)y), 0.0));	
-			std::cout << x_v/16.0 << " : " << y_v/216.0 << std::endl; 
 		}
 	}
@@ -164,6 +164,8 @@ namespace stim{
 			//free allocated memory
 			cudaFree(gpuDir);
+			cudaDestroyTextureObject(texObj);
+
 		}
 		template<typename T>
@@ -211,4 +213,4 @@ namespace stim{
 	}
 }
-#endif
 \ No newline at end of file
+#endif
@@ -124,6 +124,9 @@ namespace stim{
 			cuda_vote <<< blocks, threads,share_bytes >>>(gpuVote, texObj, gpuTable, phi, rmax, x , y);
+			cudaDestroyTextureObject(texObj);
+			cudaFreeArray(cuArray);
+
 		}
@@ -169,4 +172,4 @@ namespace stim{
 	}
 }
-#endif
 \ No newline at end of file
+#endif
@@ -59,7 +59,7 @@ namespace stim{
 					//perform the copy
 					if(sx < X && sy < Y)
-						dest[sy * X + sx] = tex2D<D>(src, tx, ty);
+						dest[sy * X + sx] = abs(255 - tex2D<D>(src, tx, ty));
 				}
 			}
 		}
@@ -155,8 +155,6 @@ namespace stim{
 		t.UnmapCudaTexture();
 		cleanUP();
 		ret[0] = idx; ret[1] = (int) output[idx];
-		std::cout << output[idx] << std::endl;
-		
 		free(output);
 		return ret;
 	}
@@ -107,6 +107,8 @@ namespace stim{
 			//call the kernel to do the multiplication
 			//cuda_conv2 <<< blocks, threads >>>(img, mask, copy, w, h, M);
 			cuda_conv2 <<< blocks, threads >>>(img, mask, copy, texObj, w, h, M);
+			cudaDestroyTextureObject(texObj);
+			cudaFreeArray(cuArray);
 		}
@@ -143,4 +145,4 @@ namespace stim{
 }
-#endif
 \ No newline at end of file
+#endif
@@ -213,6 +213,8 @@ namespace stim{
 			//free allocated memory
 			cudaFree(cuArray);
+			cudaDestroyTextureObject(texObj);
+
 		}
 		/// Applies a Gaussian blur to a 2D image stored on the CPU
@@ -43,6 +43,7 @@ namespace stim{
 			stim::cuda::tex_conv2sep(out, x, y, texObj, cuArray, gpuKernel0, kwidth, gpuKernel0, kwidth);
 			HANDLE_ERROR(cudaFree(gpuKernel0));
+			free(kernel0);
 		}
@@ -58,6 +59,7 @@ namespace stim{
 			//copy the kernel to the GPU
 			T* gpuKernel0;
+			HANDLE_ERROR(cudaMalloc(&gpuKernel0, kwidth*sizeof(T)));
 			HANDLE_ERROR(cudaMemcpy(gpuKernel0, kernel0, kwidth * sizeof(T), cudaMemcpyHostToDevice));
 			//perform the gaussian blur as a separable convolution
@@ -39,7 +39,7 @@
 		int y   = threadIdx.y + blockIdx.y * blockDim.y;
 		int idx = y*16+x;
-		float valIn             = tex2D<unsigned char>(texIn, x, y)/255.0;
+		float valIn             = tex2D<unsigned char>(texIn, x, y);
 		print[idx]              = abs(valIn);             ///temporary
@@ -80,7 +80,7 @@
 		cudaDeviceSynchronize();
 		stringstream name;      //for debugging
 		name << "FromTex.bmp";
-		stim::gpu2image<float>(print, name.str(),16,216,0,1);
+		stim::gpu2image<float>(print, name.str(),16,216,0,255);
 		tx.UnmapCudaTexture();
 		cleanUP();
@@ -70,13 +70,14 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		GLuint btexbufferID;
 		int numSamples;				//The number of templates in the buffer.
-		float stepsize = 6.0;			//Step size.
+		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::vec<float> rev;			//reverse vector;
@@ -84,6 +85,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		stim::vec<float> ps;
 		stim::vec<float> ups;
 		stim::vec<float> ds;
+		std::vector<stim::vec<float> > last3;
@@ -147,42 +149,43 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			setMatrix();
 			glCallList(dList+3);
-		//	test(btexbufferID, GL_TEXTURE_2D, 8*27);
 			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++)
 				{
-					std::cout << result[i] << std::endl;
-//					std::cout <<"["<< result[i][0]*16.0 << ", "
-//					<< result[i][1]*16.0 << ", " << 
-//					   result[i][2]*216.0 <<"]"  <<  std::endl;
-					
 					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);
-					stim::vec<float> cylv(
-						0.5 * cos(2*M_PI*(result[i][1]))*S[0]*R[0],
-						0.5 * sin(2*M_PI*(result[i][1]))*S[1]*R[1],
-						(result[i][0]-0.5)*S[2]*R[2],
-						0.0);
-					std::cout << cylp << std::endl;
 					cylp = cT*cylp;
-					cylv = cT*cylv;
-					std::cout << cylp[0]*S[0]*R[0] << ", "
-					<< cylp[1]*S[1]*R[1] << ", " << 
-					   cylp[2]*S[2]*R[2] << std::endl;
-					setSeed(cylp[0]*S[0]*R[0],
+					stim::vec<float> vec(
+						cylp[0]*S[0]*R[0],
 						cylp[1]*S[1]*R[1],
 						cylp[2]*S[2]*R[2]);
-					cylv.norm();
-					setSeedVec(cylv[0],
-						   cylv[1],
-						   cylv[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;
+					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
+					 && 
+					 !(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);
+					}
 				}
 			}    
@@ -207,7 +210,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///Method for populating the vector arrays with sampled vectors.
 		///uses the default d vector <0,0,1>
 		void
-		genDirectionVectors(float solidAngle = 3*M_PI/2)
+		genDirectionVectors(float solidAngle = 5/M_PI*4)
 		{
 			//ofstream file;
 			//file.open("dvectors.txt");
@@ -608,6 +611,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			#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;
@@ -837,11 +841,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
 		void
-		setSeedVec(stim::vec<float> pos)
+		setSeedVec(stim::vec<float> dir)
+		{
+			seedsvecs.push(dir);
+		}
+
+		void
+		setSeedMag(float mag)
 		{
-			seedsvecs.push(pos);
+			seedsmags.push(mag);
 		}
+
 		///@param x, y, z: variables for the x, y and z coordinate of the seed
 		///Adds a seed to the seed list.
 		///Assumes that the coordinates passes are in tissue space.
@@ -864,6 +875,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			return tp;
 		}
+		stim::vec<float> 
+		getLastSeedVec()
+		{
+			stim::vec<float> tp = seedsvecs.top();
+			return tp;
+		}
+
+		float
+		getLastSeedMag()
+		{
+			float tp = seedsmags.top();
+			return tp;
+		}
+
+		void
+		popSeed()
+		{
+			seeds.pop();
+			seedsvecs.pop();
+		//	seedsmags.pop();
+		}
+
 		std::stack<stim::vec<float> >
 		getSeeds()
 		{
@@ -894,9 +927,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
                                   		((float) y), 
                                   		((float) z)));
 					seedsvecs.push(stim::vec<float>(
-						((float) x),
-						((float) y),
-						((float) z)));
+						((float) u),
+						((float) v),
+						((float) w)));
                    		}
                    	myfile.close();
 	          	} else {
@@ -956,12 +989,34 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			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;
+		}
+
+		int
+		StepP()
+		{
+			Bind();
+			CHECK_OPENGL_ERROR
+			#ifdef TESTING
+				start = std::clock();
+			#endif
+			findOptimalDirection();
+			findOptimalPosition();
+			findOptimalScale();
+			Unbind();
+			CHECK_OPENGL_ERROR
 			Bind(btexbufferID, bfboId, 27);
 			CHECK_OPENGL_ERROR
 			branchDetection();
 			CHECK_OPENGL_ERROR
 			Unbind();
 			CHECK_OPENGL_ERROR
+
 			#ifdef TESTING
 				duration_sampling = duration_sampling +
 					 (std::clock() - start) / (double) CLOCKS_PER_SEC;
@@ -1032,33 +1087,34 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{	
 			Bind();
 			rev = stim::vec<float>(0.0,0.0,1.0);
+			bool sEmpty = true;
+			float lastmag = 16.0;;
 			while(!seeds.empty())
 			{
 				//clear the currently traced line and start a new one.
 				cL.clear();
 				sk.Begin(stim::OBJ_LINE);
 				stim::vec<float> curSeed = seeds.top();
+//				std::cout << "The current seeds is " << curSeed << std::endl;
 				stim::vec<float> curSeedVec = seedsvecs.top();
+				seeds.pop();
+				seedsvecs.pop();
+//				std::cout << "The current seed Vector is " << curSeedVec << std::endl;
 				setPosition(curSeed);
 				setDirection(curSeedVec);
-				setMagnitude(16.0);
 				cL.push_back(curSeed);
 				sk.createFromSelf(GL_SELECT);
 				traceLine(min_cost);
 				sk.rev();
+		//		std::cout << "reversed" << std::endl;
 				std::reverse(cL.begin(), cL.end());
 				setPosition(curSeed);
-				setDirection(rev);
+				setDirection(-rev);
 				setMagnitude(16.0);
 				sk.createFromSelf(GL_SELECT);
 				traceLine(min_cost);
-				
-				//temporary glObj rendering code.
-
 				sk.End();
-				seeds.pop();
-				seedsvecs.pop();
 			}
 			Unbind();
 		}
@@ -1072,21 +1128,23 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			stim::vec<float> mag;
 			bool h;
 			bool started = false;
+			bool running = true;
 			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
-			while(1)
+			while(running)
 			{
 				int cost = Step();
 				if (cost > min_cost){
-//				       std::cout << "Min Cost" << std::endl;
+					running = false;
 					break;
 				} else {
 					//Have we found an edge?
 					pos = getPosition();
 					if(pos[0] > size[0] || pos[1] > size[1]
 					 || pos[2] > size[2] || pos[0] < 0
-					 || p[1] < 0 || p[2] < 0)
+					 || pos[1] < 0 || pos[2] < 0)
 					{
 //					       std::cout << "Found Edge" << std::endl;
+						running = false;
 						break;
 					}
 					//If this is the first step in the trace,
@@ -1101,6 +1159,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					//Has the template size gotten unreasonable?
 					if(m[0] > 75 || m[0] < 1){
 //						std::cout << "Magnitude Limit" << std::endl;
+						running = false;
 						break;
 					}
 					else
@@ -1108,12 +1167,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						h = selectObject(pos, getDirection(), m[0]);
 						//Have we hit something previously traced?
 						if(h){
-//						std::cout << "Hit a Previous Line" << std::endl;
+						running = false;
 						break;
 						}
 						else {          
 							sk.TexCoord(m[0]);
 							sk.Vertex(p[0], p[1], p[2]);
+							Bind(btexbufferID, bfboId, 27);
+							CHECK_OPENGL_ERROR
+							branchDetection();
+							CHECK_OPENGL_ERROR
+							Unbind();
+							CHECK_OPENGL_ERROR
 						}
 				 	}
                  		}
@@ -1132,7 +1197,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			glSelectBuffer(2048, selectBuf);
 			glDisable(GL_CULL_FACE);
 			(void) glRenderMode(GL_SELECT);
-
 		//Init Names stack
 			glInitNames();
@@ -1151,7 +1215,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				glMatrixMode(GL_PROJECTION);
 				glPushMatrix();
 				glLoadIdentity();
-				glOrtho(-mag/s, mag/s, -mag/s, mag/s, 0.0, mag/s/4.0);
+				glOrtho(-mag/s, mag/s, -mag/s, mag/s, 0.0, mag/s/2.0);
 				glMatrixMode(GL_MODELVIEW);
 				glPushMatrix();
 				glLoadIdentity();
@@ -1173,7 +1237,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				CHECK_OPENGL_ERROR
 				glPopMatrix();
-			glEnable(GL_CULL_FACE);
+		//	glEnable(GL_CULL_FACE);
 			hits = glRenderMode(GL_RENDER);
 			bool found_hits = processHits(hits, selectBuf);
 			return found_hits;
@@ -1185,19 +1249,19 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		processHits(GLint hits, GLuint buffer[])
 		{
 			GLuint names, *ptr;
-			printf("hits = %u\n", hits);
+			//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);
+			//	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");
+			//	printf ("   the name is ");                                                        
+			//	for (int j = 0; j < names; j++) {     /*  for each name */
+			//		printf ("%u ", *ptr); ptr++;
+			//	}
+			//	printf ("\n");
 			}
 			if(hits == 0)
 				return 0;