added debugging code to the cuda code that is enabled with -DDebug in cmake list…

…s. Fixed minor errors in the grids and image stack classes

added debugging code to the cuda code that is enabled with -DDebug in cmake list…
…s. Fixed minor errors in the grids and image stack classes
Pavel Govyadinov
1 parent efb0a2b5
Showing 7 changed files with 121 additions and 28 deletions Show diff stats
stim/cuda/branch_detection.cuh
stim/cuda/filter.cuh
stim/cuda/testKernel.cuh
stim/envi/envi_header.h
stim/gl/gl_spider.h
stim/gl/gl_texture.h
stim/grids/image_stack.h
@@ -11,7 +11,7 @@ typedef unsigned int uint;
 std::vector< stim::vec<float> >
-find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
+find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y, int iter = 0)
 {
 	float		sigma		= 2.0;
 	unsigned int	conn		= 7;
@@ -22,7 +22,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	stringstream name;
-	cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, threshold);
+	cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, threshold, iter);
 	cudaDeviceSynchronize();
@@ -26,6 +26,11 @@ namespace stim
 	float* LoG;
 	float* res;
 	float* centers;
+
+//#ifdef DEBUG
+	float* print;
+//#endif
+
 	stim::cuda::cuda_texture tx;
@@ -44,6 +49,13 @@ namespace stim
 		HANDLE_ERROR(
 			cudaMalloc( (void**) &centers, DIM_Y*DIM_X*sizeof(float))
 		);
+
+//#ifdef DEBUG
+		HANDLE_ERROR(
+			cudaMalloc( (void**) &print, DIM_Y*DIM_X*sizeof(float))
+		);
+//#endif
+
 	//	checkCUDAerrors("Memory Allocation, Result");
 	}
@@ -58,6 +70,11 @@ namespace stim
 		HANDLE_ERROR(
 			cudaFree(centers)
 		);
+//#ifdef DEBUG
+		HANDLE_ERROR(
+			cudaFree(print)
+		);
+//#endif
 			free(LoG);
 	}
@@ -89,7 +106,7 @@ namespace stim
 	//Shared memory would be better.
 	__global__
 	void
-	applyFilter(cudaTextureObject_t texIn, unsigned int DIM_X, unsigned int DIM_Y, int kr, int kl, float *res, float* gpuLoG){
+	applyFilter(cudaTextureObject_t texIn, unsigned int DIM_X, unsigned int DIM_Y, int kr, int kl, float *res, float* gpuLoG, float* p){
 	//R = floor(size/2)
 	//THIS IS A NAIVE WAY TO DO IT, and there is a better way)
@@ -101,16 +118,15 @@ namespace stim
 	//	float val = 0;
 		float tu = (x-kr+xi)/(float)DIM_X;
 		float tv = (y-kr+yi)/(float)DIM_Y;
+		int idx = y*DIM_X+x;
 		shared[xi][yi] = gpuLoG[yi*kl+xi]*(255.0-(float)tex2D<unsigned char>(texIn, tu, tv));
 		__syncthreads();
-	
 		//x = max(0,x);
 		//x = min(x, width-1);
 		//y = max(y, 0);
 		//y = min(y, height - 1);
-		int idx = y*DIM_X+x;
 	//	int k_idx;
                 for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
                 {
@@ -135,11 +151,12 @@ namespace stim
                 __syncthreads();
                 if(xi == 0 && yi == 0)
                         res[idx] = shared[0][0];
+
 	}
 	extern "C"
 	float *
-	get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, float threshold)
+	get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, float threshold, int iter = 0)
 	{
 		tx.SetTextureCoordinates(1);
 		tx.SetAddressMode(1, 3);
@@ -153,7 +170,14 @@ namespace stim
 		dim3 numBlocks(DIM_X, DIM_Y);
 		dim3 threadsPerBlock(sizeK, sizeK);
-		applyFilter <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), DIM_X, DIM_Y, floor(sizeK/2), sizeK, res, gpuLoG);
+		applyFilter <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), DIM_X, DIM_Y, floor(sizeK/2), sizeK, res, gpuLoG, print);
+
+		#ifdef DEBUG
+			stringstream name;
+			name.str("");
+			name << "Fiber Cylinder " << iter << ".bmp";
+			stim::gpu2image<float>(res, name.str(), DIM_X, DIM_Y, 0, 255);
+		#endif
 		stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
@@ -53,11 +53,31 @@
 		float valIn             = tex2D<unsigned char>(texIn, x, y);
 		float templa		= templ(x, 32)*255.0;
-		print[idx]              = abs(valIn-templa);             ///temporary
+		//print[idx]              = abs(valIn-templa);             ///temporary
+		print[idx]		= abs(valIn);
 		//print[idx]              = abs(templa);             ///temporary
 	}
+	///Find the difference of the given set of samples and the template
+	///using cuda acceleration.
+	///@param stim::cuda::cuda_texture t	--stim texture that holds all the references
+	///					  to the data.
+	///@param float* result			--a pointer to the memory that stores the result.
+	__global__
+	//void get_diff (float *result)
+	void get_diff2 (cudaTextureObject_t texIn, float *print, int dx)
+	{       
+		int x   = threadIdx.x + blockIdx.x * blockDim.x;
+		int y   = threadIdx.y + blockIdx.y * blockDim.y;
+		int idx = y*dx+x;
+	//	int idx = y*16+x;
+
+		float valIn             = tex2D<unsigned char>(texIn, x, y);
+		print[idx]              = abs(valIn);             ///temporary
+
+	}
+
 	void test(cudaTextureObject_t tObj, int x, int y, std::string nam)
 	{
@@ -86,3 +106,31 @@
 		cleanUP();
 	}
+
+	void test(cudaTextureObject_t tObj, int x, int y, std::string nam, int iter)
+	{
+
+		//Bind the Texture in GL and allow access to cuda.
+
+		//initialize the return arrays.
+
+		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(2, 2);
+		//dim3 threadsPerBlock(8, 108);
+
+
+//		get_diff <<< blocks, threads >>> (tx.getTexture(), print);
+		get_diff2 <<< numBlocks, threadsPerBlock >>> (tObj, print, x);
+
+		cudaDeviceSynchronize();
+		stringstream name;      //for debugging
+		name << nam.c_str();
+		stim::gpu2image<float>(print, name.str(),x,y,0,255);
+	  
+		cleanUP();
+	}
@@ -8,6 +8,7 @@
 #include <vector>
 #include <algorithm>
 #include <stdlib.h>
+#include <cmath>
 //information from an ENVI header file
 //A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
@@ -27,7 +27,6 @@
 #include <stim/cuda/branch_detection.cuh>
 #include "../../../volume-spider/glnetwork.h"
 #include <stim/visualization/cylinder.h>
-#include <stim/cuda/testKernel.cuh>
 #include <iostream>
 #include <fstream>
 #ifdef TIMING
@@ -40,6 +39,9 @@
 	#include <ctime>
 #endif
+#ifdef DEBUG
+	#include <stim/cuda/testKernel.cuh>
+#endif
 namespace stim
 {
@@ -143,8 +145,8 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		#ifdef DEBUG
-			stringstream name;
 			int iter;
+			stringstream name;
 			int iter_pos;
 			int iter_dir;
 			int iter_siz;
@@ -294,6 +296,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			DrawLongCylinder(n, l_template, l_square);	///Draw the cylinder.
 			stim::cylinder<float> cyl(cL, cM);
 			std::vector< stim::vec<float> > result = find_branch(cylinder_texID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);		///find all the centers in cuda
+			
 			stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);			///the borders of the texture.
 			float pval;									//pvalue associated with the points on the cylinder.
 			if(!result.empty())								///if we have any points
@@ -957,7 +960,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			iter_dir = 0;
 			iter_siz = 0;
 #endif
-			stepsize = 3.0;
+			stepsize = 6.0;
 			n_pixels = 16.0;
 			srand(100);	
@@ -1380,20 +1383,31 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		Step()
 		{
 			#ifdef DEBUG
-			std::cerr << "Took a step" << std::endl;
+			std::cerr << "Took a step";
 			#endif
 			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
 			CHECK_OPENGL_ERROR
 				findOptimalDirection();
 			Unbind();
+			#ifdef DEBUG
+			std::cerr << " " << current_cost;
+			#endif
 			Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
 				findOptimalPosition();
 			Unbind();
+			#ifdef DEBUG
+			std::cerr << " " << current_cost;
+			#endif
 			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
 				findOptimalRadius();
 			Unbind();
+			#ifdef DEBUG
+			std::cerr << " " << current_cost;
+			#endif
 			CHECK_OPENGL_ERROR
-
+			#ifdef DEBUG
+			std::cerr << std::endl;
+			#endif
 			return current_cost;
 		}
@@ -1539,9 +1553,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 //				findOptimalDirection();
 //			Unbind();
 //THIS IS EXPERIMENTAL
-			Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
-				findOptimalRadius();
-			Unbind();
+		//	Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
+		//		findOptimalRadius();
+		//	Unbind();
 //THIS IS EXPERIMENTAL
 //				cL.push_back(curSeed);
@@ -30,9 +30,9 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 		GLenum cpu_type;
 		GLenum gpu_type;
 		GLenum format;					//format for the texture (GL_RGBA, GL_LUMINANCE, etc.)
-		using image_stack<T>::R;
+		using image_stack<T,F>::R;
 		//using image_stack<T>::S;
-		using image_stack<T>::ptr;
+		using image_stack<T,F>::ptr;
 		///	Sets the internal texture_type, based on the data dimensions
 		void setTextureType(){
@@ -247,7 +247,7 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 		}
 		///returns the dimentions of the data in the x, y, z directions. 
-		vec<int> getSize(){
+		stim::vec<int> getSize(){
 			stim::vec<int> size(R[1], R[2], R[3]);
 			return size;
 		}
@@ -282,7 +282,7 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 		///@param file_mask specifies the file(s) to be loaded
 		///	Sets the path and calls the loader on that path.
 		void load_images(std::string file_mask){
-			image_stack<T>::load_images(file_mask);				//load images
+			image_stack<T, F>::load_images(file_mask);				//load images
 			guess_parameters();
 		}
@@ -292,13 +292,18 @@ class gl_texture : public virtual image_stack&lt;T, F&gt;
 			return texID;
 		}
-		
+	
 		T* getData(){
 			return ptr;
 		}
-		
-		
+	
+		void setData(T* rts)
+		{
+				
+		}
+
 	};
+
 }
@@ -17,13 +17,13 @@ class image_stack : public virtual stim::grid&lt;T, 4, F&gt;{
 protected:
 	//using stim::grid<T, 4>::S;
-	using stim::grid<T, 4>::R;
-	using stim::grid<T, 4>::ptr;
-	using stim::grid<T, 4>::read;
+	using stim::grid<T, 4, F>::R;
+	using stim::grid<T, 4, F>::ptr;
+	using stim::grid<T, 4, F>::read;
 public:
 	//default constructor
-	image_stack() : grid<T, 4>() {
+	image_stack() : grid<T, 4, F>() {
 	}
@@ -113,7 +113,8 @@ public:
 	/// @param i is the page to be saved
 	void save_image(std::string file_name, unsigned int i){		
 		stim::image<T> I;											//create an image		
-		I.set_interleaved_rgb(&ptr[ i * R[0] * R[1] * R[2] ], R[1], R[2], R[0]);	//retrieve the interlaced data from the image - store it in the grid
+		I.set_interleaved(&ptr[ i * R[0] * R[1] * R[2] ], R[1], R[2], R[0]);	//retrieve the interlaced data from the image - store it in the grid
+//		I.set_interleaved_rgb(&ptr[ i * R[0] * R[1] * R[2] ], R[1], R[2], R[0]);	//retrieve the interlaced data from the image - store it in the grid
 		I.save(file_name);
 	}