Merge remote-tracking branch 'origin/Branch_Detection'

Pavel Govyadinov
2 parents c230039b c0f3e9f6
Showing 4 changed files with 214 additions and 138 deletions Show diff stats
stim/cuda/cost.h
stim/gl/gl_spider.h
stim/visualization/glObj.h
stim/visualization/obj.h
@@ -3,17 +3,16 @@
 #include <cuda_runtime.h>
 #include <cublas_v2.h>
 #include <stdio.h>
-#include <stim/visualization/colormap.h>
+#include "../visualization/colormap.h"
 #include <sstream>
-#include <stim/math/vector.h>
-#include <stim/cuda/devices.h>
-#include <stim/cuda/threads.h>
+#include "../math/vector.h"
+#include "../cuda/devices.h"
+#include "../cuda/threads.h"
 ///Cost function that works with the gl-spider class to find index of the item with min-cost.
 typedef unsigned char uchar;
 texture<uchar, cudaTextureType2D, cudaReadModeElementType> texIn;
 float *result;
-float* v_dif;
 cudaArray* srcArray;
 bool testing = false;
@@ -26,30 +25,12 @@ inline void checkCUDAerrors(const char *msg)
 		}
 }
-///Finds the sum of all the pixes in a gives template element.
-///Returns the abosolute value.
-///@param *diff, a pointer to the memory block that holds the pixel-differences.
-float get_sum(float *diff)
-{
-
-	cublasStatus_t ret;
-	cublasHandle_t handle;
-	ret = cublasCreate(&handle);
-	
-	ret = cublasSetVector(20*10, sizeof(*diff), diff, 1, v_dif, 1);
-	float out;
-	ret = cublasSasum(handle, 20*10, v_dif, 1, &out);
-//	cublasDestroy(ret);
-	cublasDestroy(handle);
-	return out;
-}
-
 ///A virtual representation of a uniform template.
 ///Returns the value of the template pixel.
 ///@param x, location of a pixel.
 __device__ float Template(int x)
 {
-	if(x < 20/6 || x > 20*5/6 || (x > 20*2/6 && x < 20*4/6)){
+	if(x < 16/6 || x > 16*5/6 || (x > 16*2/6 && x < 16*4/6)){
 		return 1.0;
 	}else{
 		return 0.0;
@@ -63,15 +44,47 @@ __device__ float Template(int x)
 __global__
 void get_diff (float *result)
 {	
-	//cuPrintf("Hello");
+	//float* shared = SharedMemory();
+	__shared__ float shared[16][8];
 	int x 	= threadIdx.x + blockIdx.x * blockDim.x;
 	int y 	= threadIdx.y + blockIdx.y * blockDim.y;
-	int idx = y*20+x;	
+	int x_t = threadIdx.x;
+	int y_t = threadIdx.y;
+	//int idx = y*16+x;
+	int g_idx = blockIdx.y;
 	float valIn		= tex2D(texIn, x, y)/255.0;
 	float valTemp		= Template(x);
-	result[idx]		= abs(valIn-valTemp);
-	//result[idx]		= abs(valIn);
+	shared[x_t][y_t]	= abs(valIn-valTemp);
+
+	__syncthreads();
+
+	for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
+	{
+		__syncthreads();
+		if (x_t < step)
+		{
+			shared[x_t][y_t] += shared[x_t + step][y_t];
+		}
+	__syncthreads();
+	}
+	__syncthreads();
+
+	for(unsigned int step = blockDim.y/2; step >= 1; step >>= 1)
+	{
+		__syncthreads();
+		if(y_t < step)
+		{
+			shared[x_t][y_t] += shared[x_t][y_t + step];
+		}
+	__syncthreads();
+	}
+	__syncthreads();
+	if(x_t == 0 && y_t == 0)
+		result[g_idx] = shared[0][0];
+
+
+//	//result[idx]		= abs(valIn);
 }
@@ -82,12 +95,6 @@ void get_diff (float *result)
 ///@param DIM_Y, integer controlling how much memory to allocate.
 void initArray(cudaGraphicsResource_t src, int DIM_Y)
 {
-	//cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar> ();
-	//cudaMallocArray(&result, &channelDesc, DIM_X, DIM_Y, 0);
-	//HANDLE_ERROR(
-	//	cudaGraphicsGLRegisterImage(&src,
-	//		fboID,
-	//		GL_TEXTURE_2D,
 	HANDLE_ERROR(
 		cudaGraphicsMapResources(1, &src)	
 	);
@@ -97,10 +104,8 @@ void initArray(cudaGraphicsResource_t src, int DIM_Y)
 	HANDLE_ERROR(
 		cudaBindTextureToArray(texIn, srcArray)
 		);
-	cudaMalloc( (void**) &result, 20*DIM_Y*sizeof(float));
+	cudaMalloc( (void**) &result, DIM_Y*sizeof(float));
         checkCUDAerrors("Memory Allocation Issue 1");	
-	cudaMalloc((void **) &v_dif, 20*10*sizeof(float));
-        checkCUDAerrors("Memory Allocation Issue 2");	
 	//HANDLE_ERROR(
 	//	cudaBindTextureToArray(texIn, ptr, &channelDesc)
 	//	    );		
@@ -117,9 +122,6 @@ void cleanUP(cudaGraphicsResource_t src)
 		cudaGraphicsUnmapResources(1,&src)
 	);
 	HANDLE_ERROR(
-		cudaFree(v_dif)
-	);
-	HANDLE_ERROR(
 		cudaUnbindTexture(texIn)
 	);
 }
@@ -151,25 +153,33 @@ stim::vec&lt;int&gt; get_cost(cudaGraphicsResource_t src, int DIM_Y)
 //	name << "sample_" << inter << "_" << idx << ".bmp";
 //	stim::gpu2image<float>(v_dif, name.str(), 20,10,0,1);
-	float output[DIM_Y];
+	//float output[DIM_Y];
+	float *output;
+	output = (float* ) malloc(DIM_Y*sizeof(float));
 	stim::vec<int> ret(0, 0);
 	float mini = 10000000000000000.0;
-	int idx;
-	initArray(src, DIM_Y*10);
-	dim3 grid(20/2, DIM_Y*10/2);
-	dim3 block(2, 2);
-
-	get_diff <<< grid, block >>> (result);
-	for (int i = 0; i < DIM_Y; i++){
-		output[i] = get_sum(result+(20*10*i));
-		if(output[i] <= mini){
+	int idx = 0;
+	initArray(src, DIM_Y*8);
+	dim3 numBlocks(1, DIM_Y);
+	dim3 threadsPerBlock(16, 8);
+
+
+	get_diff <<< numBlocks, threadsPerBlock >>> (result);
+	cudaMemcpy(output, result, DIM_Y*sizeof(float), cudaMemcpyDeviceToHost);
+
+	for( int i = 0; i<DIM_Y; i++){
+//		std::cout << output[i] << std::endl;
+		if(output[i] < mini){
 			mini = output[i];
 			idx = i;
 		}
-	}	
-
-	output[idx] = get_sum(result+(20*10*idx));
+	}
+  
+//	std::cout << "hello" << std::endl;
+	//output[idx] = get_sum(result+(16*8*idx));
 	cleanUP(src);
 	ret[0] = idx; ret[1] = (int) output[idx];
+	std::cout << output[idx] << std::endl;
+	free(output);
 	return ret;
 }
@@ -121,7 +121,7 @@ class gl_spider
 			setMatrix();
 			glCallList(dList+3);
-//			int best = getCost();
+		//	int best = getCost();
 		}
@@ -185,7 +185,7 @@ class gl_spider
        				 	ver = stim::rect<float>(mag,
 						 pos, temp,
                        				hor.n());
-					UpdateBuffer(0.0, 0.0+idx*10.0);
+					UpdateBuffer(0.0, 0.0+idx*8.0);
 					CHECK_OPENGL_ERROR
 				}
 			}
@@ -233,7 +233,7 @@ class gl_spider
                 			ver = stim::rect<float>(mag,
 						 temp, dir,
                         			hor.n());
-					UpdateBuffer(0.0, 0.0+idx*10.0);
+					UpdateBuffer(0.0, 0.0+idx*8.0);
 				CHECK_OPENGL_ERROR
 				}
 			}
@@ -244,7 +244,8 @@ class gl_spider
 		///Method for populating the buffer with the sampled texture.
 		///uses the default m <1,1,0>
 		void
-		genMagnitudeVectors(float delta = 0.5)
+		genMagnitudeVectors(float delta = 0.70)
+//		genMagnitudeVectors(float delta = 0.50)
 		{
 			//Set up the vectors necessary for Rectangle creation.
@@ -274,7 +275,7 @@ class gl_spider
                			ver = stim::rect<float>(temp,
 					 pos, dir,
                        			hor.n());
-				UpdateBuffer(0.0, 0.0+i*10.0);
+				UpdateBuffer(0.0, 0.0+i*8.0);
 			CHECK_OPENGL_ERROR
 			}
 			glEndList();
@@ -286,7 +287,7 @@ class gl_spider
 		void
 		UpdateBuffer(float v_x, float v_y)
 		{	
-			float len = 10.0;
+			float len = 8.0;
 			stim::vec<float>p1; 
         	        stim::vec<float>p2; 
 	                stim::vec<float>p3; 
@@ -338,13 +339,13 @@ class gl_spider
 					p2[1],
 					p2[2]
 					);
-				glVertex2f(v_x+2*len, v_y);
+				glVertex2f(v_x+2.0*len, v_y);
 				glTexCoord3f(
 					p3[0],
 					p3[1],
 					p3[2]
 					);
-				glVertex2f(v_x+2*len, v_y+len);
+				glVertex2f(v_x+2.0*len, v_y+len);
 				glTexCoord3f(
 					p4[0],
 					p4[1],
@@ -383,47 +384,6 @@ class gl_spider
 			glBindTexture(GL_TEXTURE_2D, 0);
 		}
-		///Method for controling the buffer and texture binding in order to properly
-		///do the render to texture.
-		void
-		Bind()
-		{
-			float len = 10.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 Unbinding all of the texture resources
-		void
-		Unbind()
-		{
-			//Finalize GL_buffer
-			glBindTexture(GL_TEXTURE_3D, 0);                      
-			glDisable(GL_TEXTURE_3D);
-			glBindFramebuffer(GL_FRAMEBUFFER,0);
-			glBindTexture(GL_TEXTURE_2D, 0);
-		}
 		///Method for using the gl manipulation to alighn templates from
 		///Template space (-0.5 0.5) to Texture space (0.0, 1.0),
@@ -558,7 +518,7 @@ class gl_spider
 		attachSpider(GLuint id)
 		{
 			texID = id;
-			GenerateFBO(20, numSamples*10);
+			GenerateFBO(16, numSamples*8);
 			setDims(0.6, 0.6, 1.0);
 			setSize(512.0, 512.0, 426.0);
 			setMatrix();
@@ -704,6 +664,47 @@ class gl_spider
 			return fboID;
 		}
+		///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 Unbinding all of the texture resources
+		void
+		Unbind()
+		{
+			//Finalize GL_buffer
+			glBindTexture(GL_TEXTURE_3D, 0);                      
+			glDisable(GL_TEXTURE_3D);
+			glBindFramebuffer(GL_FRAMEBUFFER,0);
+			glBindTexture(GL_TEXTURE_2D, 0);
+		}
 //--------------------------------------------------------------------------//
 //-----------------------------TEMPORARY METHODS----------------------------//
 //--------------------------------------------------------------------------//
@@ -729,7 +730,7 @@ class gl_spider
 			findOptimalDirection();
 			findOptimalPosition();
 			findOptimalScale();
-	//		branchDetection();
+			branchDetection();
 			Unbind();
 			return current_cost;
 		}
@@ -776,9 +777,9 @@ class gl_spider
 					 glTexCoord3f(x,y,z0); 
 					 glVertex2f(0.0, j*0.1+0.1);
 					 glTexCoord3f(x,y,z1); 
-					 glVertex2f(20.0, j*0.1+0.1);
+					 glVertex2f(16.0, j*0.1+0.1);
 					 glTexCoord3f(xold,yold,z1); 
-					 glVertex2f(20.0, j*0.1); 
+					 glVertex2f(16.0, j*0.1); 
 					 glTexCoord3f(xold,yold,z0); 
 					 glVertex2f(0.0, j*0.1);
 					 xold=x;
 #ifndef STIM_GLOBJ_H
 #define STIM_GLOBJ_H
-#include <stim/visualization/obj.h>
 #include <GL/glew.h>
 #include <GL/glut.h>
 #include <stim/math/vector.h>
+#include <stim/visualization/obj.h>
 namespace stim
@@ -13,23 +13,25 @@ namespace stim
  *  Assist with the visualization the segmented vessels.
  *  Uses
 */
-
-class objGl : public virtual stim::obj<T>
+template <typename T>
+class glObj : public virtual stim::obj<T>
 {
 private:
 	using stim::obj<T>::L;
 	using stim::obj<T>::P;
 	using stim::obj<T>::F;
-	using vec<T>::size;
-	using vec<T>::at;
+//	using stim::obj<T>::numL();
+//	using std::vector<T>::size;
+//	using std::vector<T>::at;
 	GLuint dList;
 	void
-	init
+	init()
 	{
-		dList = glGenList(1);
+		dList = glGenLists(1);
 		glListBase(dList);
+
 		glMatrixMode(GL_PROJECTION);
 		glLoadIdentity;
 		glMatrixMode(GL_MODELVIEW);
@@ -40,15 +42,17 @@ private:
 	void
 	Create()
 	{
-		int len = (int) numL();
+		int len = (int) stim::obj<T>::numL();
 		std::vector< stim::vec<float> > line;
 		glNewList(dList, GL_COMPILE);
-		glColor3f(0.5, 1.0, 0.5);
+	//	glColor3f(0.0, 1.0, 0.0);
+		glLineWidth(2.5);
 		for(int i = 0; i < len; i++){
-			line = getL_V(i);
-			glBegin(GL_LINES);
+			line = stim::obj<T>::getL_V(i);
+			glColor3ub(rand()%255, rand()%255, rand()%255);
+			glBegin(GL_LINE_STRIP);
 			for(int j = 0; j < line.size(); j++){
-					glVectex3f(
+					glVertex3f(
 						line[j][0],
 						line[j][1],
 						line[j][2]
@@ -57,16 +61,33 @@ private:
 			glEnd();
 		}
 		glEndList();
+		CHECK_OPENGL_ERROR
 	}
 public:
-	objGl(std::string filename)
+	glObj()
 	{
-		stim::obj::load(filename);
-		glPopMatrix();		//Safety Operation to avoid changing the current matrix.
+
+	}
+
+	void
+	createFromSelf()
+	{
+	//	glPopMatrix();
 		init();
 		Create();
-		glPushMatrix();
+	//	glPushMatrix();
+	}
+
+	void
+	createFromFile(std::string filename)
+	{
+		stim::obj<T>::load(filename);
+		glPushMatrix();		//Safety Operation to avoid changing the current matrix.
+		init();
+		Create();
+		glPopMatrix();
+		CHECK_OPENGL_ERROR
 	}
@@ -74,7 +95,44 @@ public:
 	Render()
 	{
 		glCallList(dList);
+		CHECK_OPENGL_ERROR
 	}
+	void
+	RenderLine(int i)
+	{
+		std::vector < stim::vec<T> > line;
+		int len = (int) stim::obj<T>::numL();
+		line = stim::obj<T>::getL_V(i);
+		glColor3f(0.5, 1.0, 0.5);
+		glLineWidth(3.0);
+		glBegin(GL_LINE_STRIP);
+		for(int j = 0; j < line.size(); j++){
+			glVertex3f(
+				line[j][0],
+				line[j][1],
+				line[j][2]
+				);
+		}
+		glEnd();
+	}
+
+	void
+	RenderLine(std::vector< stim::vec<T> > l)
+	{
+		glColor3f(0.5, 1.0, 0.5);
+		glLineWidth(3.0);
+		glBegin(GL_LINE_STRIP);
+		for(int j = 0; j < l.size(); j++){
+			glVertex3f(
+				l[j][0],
+				l[j][1],
+				l[j][2]
+				);
+		}
+		glEnd();
+	}
+
+};
 }
-}
+#endif
@@ -352,26 +352,27 @@ public:
 	/// This function terminates drawing of a primitive object, such as a line, face, or point set
 	void End(){
-
 		//copy the current object to the appropriate list
-		switch(current_type){
+		if(current_geo.size() != 0)
+		{
+			switch(current_type){
-		case OBJ_NONE:
-			std::cout<<"STIM::OBJ error, objEnd() called before objBegin()."<<std::endl;
-			break;
+			case OBJ_NONE:
+				std::cout<<"STIM::OBJ error, objEnd() called before objBegin()."<<std::endl;
+				break;
-		case OBJ_POINTS:
-			P.push_back(current_geo);
-			break;
+			case OBJ_POINTS:
+				P.push_back(current_geo);
+				break;
-		case OBJ_LINE:
-			L.push_back(current_geo);
-			break;
+			case OBJ_LINE:
+				L.push_back(current_geo);
+				break;
-		case OBJ_FACE:
-			F.push_back(current_geo);
+			case OBJ_FACE:
+				F.push_back(current_geo);
+			}
 		}
-
 		//clear everything
 		current_type = OBJ_NONE;
 		current_geo.clear();
@@ -380,6 +381,12 @@ public:
 		geo_flag_vt = false;
 		geo_flag_vn = false;
 	}
+//temporary convenience method
+	void rev(){
+		if(current_geo.size() != 0)
+		//	current_geo.reverse(current_geo.begin(), current_geo.end());
+			std::reverse(current_geo.begin(), current_geo.end());
+	}
 	//output the OBJ structure as a string
 	std::string str(){
@@ -642,7 +649,6 @@ public:
 		//get the number of points in the specified line
 		unsigned int nP = L[i].size();
-
 		//create a vector of points
 		std::vector< stim::vec<T> > l;
@@ -688,6 +694,7 @@ public:
 		return l;
 	}
+
 	/// Returns a vector containing the list of texture coordinates associated with each point of a line.
 	/// @param i is the index of the desired line