Lab5: Finished hopefully

Lab5/filter.cu

@@ -55,8 +55,8 @@ __global__ void filter_1(unsigned char *image, unsigned char *out, const unsigne
 
   while(local_index < LOCAL_MEM_SIZE)
   {
-    int global_x = max(0, blockOffsetX - kernelsizex + local_index % FILTER_W);
-	  int global_y = max(0, blockOffsetY - kernelsizey + local_index / FILTER_W);
+    int global_x = min(max(0, blockOffsetX - kernelsizex + local_index % FILTER_W), imagesizex - 1);
+	  int global_y = min(max(0, blockOffsetY - kernelsizey + local_index / FILTER_W), imagesizey - 1);
     unsigned global_index = global_y * imagesizex + global_x;
 
     local_memory[3*local_index+0] = image[3*global_index+0];
@@ -65,6 +65,7 @@ __global__ void filter_1(unsigned char *image, unsigned char *out, const unsigne
 
     local_index += threadsPerBlock;
   }
+  __syncthreads();
 
   unsigned x = blockIdx.x * blockDim.x + threadIdx.x;
   unsigned y = blockIdx.y * blockDim.y + threadIdx.y;
@@ -97,7 +98,7 @@ __global__ void filter_1(unsigned char *image, unsigned char *out, const unsigne
   }
 }
 
-__global__ void filter_separable(unsigned char *image, unsigned char *out, const unsigned int imagesizex, const unsigned int imagesizey, const int kernelsize, const int stride)
+__global__ void filter_separable(unsigned char *image, unsigned char *out, const unsigned int imagesizex, const unsigned int imagesizey, const int kernelsize, const bool horizontal)
 {
   __shared__ unsigned char local_memory [STATIC_SHARED_MEM_SIZE_SEP * 3];
   // map from blockIdx to pixel position
@@ -106,14 +107,28 @@ __global__ void filter_separable(unsigned char *image, unsigned char *out, const
   int blockOffsetX = blockIdx.x * blockDim.x;
   int blockOffsetY = blockIdx.y * blockDim.y;
 
-  int local_index = threadIdx.x;
+  /* Threadidx.x or threadIdx.y will always be zero depending on direction */
+  int local_index = threadIdx.x + threadIdx.y;
+
 
   int FILTER_SIZE = (kernelsize * 2 + 1) + (2 * kernelsize);
 
   while(local_index < FILTER_SIZE)
   {
-    int global_x = max(0, blockOffsetX - kernelsize + local_index);
-    int global_y = max(0, blockOffsetY);
+    int global_x = 0;
+    int global_y = 0;
+
+    if(horizontal)
+    {
+      global_x = min(max(0, blockOffsetX - kernelsize + local_index), imagesizex - 1);
+      global_y = blockOffsetY;
+    }
+    else
+    {
+      global_x = blockOffsetX;
+      global_y = min(max(0, blockOffsetY - kernelsize + local_index), imagesizey - 1);//max(0, blockOffsetY);
+    }
+
     unsigned global_index = global_y * imagesizex + global_x;
 
     local_memory[3*local_index+0] = image[3*global_index+0];
@@ -122,10 +137,12 @@ __global__ void filter_separable(unsigned char *image, unsigned char *out, const
 
     local_index += threadsPerBlock;
   }
-
+  __syncthreads();
   unsigned x = blockIdx.x * blockDim.x + threadIdx.x;
   unsigned y = blockIdx.y * blockDim.y + threadIdx.y;
 
+  /* Since one of threadIdx.x or threadIdx.y is 0 we can just add them */
+  local_index = threadIdx.x + threadIdx.y;
   if (x < imagesizex && y < imagesizey) // If inside image
   {
     int d;
@@ -135,21 +152,90 @@ __global__ void filter_separable(unsigned char *image, unsigned char *out, const
    	sumx=0;sumy=0;sumz=0;
    	for(d=-kernelsize;d<=kernelsize;d++)
     {
-			unsigned dx = kernelsize + threadIdx.x + d;
+			unsigned dx = kernelsize + local_index + d;
 
       sumx += local_memory[dx*3+0];
       sumy += local_memory[dx*3+1];
       sumz += local_memory[dx*3+2];
  		}
-    out[(y*imagesizex+x)*3+0] = local_memory[kernelsize + threadIdx.x + 0];
-    out[(y*imagesizex+x)*3+1] = local_memory[kernelsize + threadIdx.x + 1];
-    out[(y*imagesizex+x)*3+2] = local_memory[kernelsize + threadIdx.x + 2];
-    // out[(y*imagesizex+x)*3+0] = sumx/divby;
-  	// out[(y*imagesizex+x)*3+1] = sumy/divby;
-  	// out[(y*imagesizex+x)*3+2] = sumz/divby;
+
+
+    out[(y*imagesizex+x)*3+0] = sumx/divby;
+  	out[(y*imagesizex+x)*3+1] = sumy/divby;
+  	out[(y*imagesizex+x)*3+2] = sumz/divby;
   }
 }
 
+__global__ void filter_gaussian(unsigned char *image, unsigned char *out, const unsigned int imagesizex, const unsigned int imagesizey, const int kernelsize, const bool horizontal)
+{
+  __shared__ unsigned char local_memory [STATIC_SHARED_MEM_SIZE_SEP * 3];
+  // map from blockIdx to pixel position
+  unsigned threadsPerBlock  = blockDim.x;
+
+  int blockOffsetX = blockIdx.x * blockDim.x;
+  int blockOffsetY = blockIdx.y * blockDim.y;
+
+  /* Threadidx.x or threadIdx.y will always be zero depending on direction */
+  int local_index = threadIdx.x + threadIdx.y;
+
+
+  int FILTER_SIZE = (kernelsize * 2 + 1) + (2 * kernelsize);
+
+  while(local_index < FILTER_SIZE)
+  {
+    int global_x = 0;
+    int global_y = 0;
+
+    if(horizontal)
+    {
+      global_x = min(max(0, blockOffsetX - kernelsize + local_index), imagesizex - 1);
+      global_y = blockOffsetY;
+    }
+    else
+    {
+      global_x = blockOffsetX;
+      global_y = min(max(0, blockOffsetY - kernelsize + local_index), imagesizey - 1);//max(0, blockOffsetY);
+    }
+
+    unsigned global_index = global_y * imagesizex + global_x;
+
+    local_memory[3*local_index+0] = image[3*global_index+0];
+    local_memory[3*local_index+1] = image[3*global_index+1];
+    local_memory[3*local_index+2] = image[3*global_index+2];
+
+    local_index += threadsPerBlock;
+  }
+  __syncthreads();
+
+  unsigned x = blockIdx.x * blockDim.x + threadIdx.x;
+  unsigned y = blockIdx.y * blockDim.y + threadIdx.y;
+
+  /* Since one of threadIdx.x or threadIdx.y is 0 we can just add them */
+  local_index = threadIdx.x + threadIdx.y;
+  if (x < imagesizex && y < imagesizey) // If inside image
+  {
+    int d;
+    unsigned int sumx, sumy, sumz;
+    float divby[5] = {1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0};
+   	sumx=0;sumy=0;sumz=0;
+
+   	for(d=-kernelsize;d<=kernelsize;d++)
+    {
+			unsigned dx = kernelsize + local_index + d;
+
+      sumx += local_memory[dx*3+0] * divby[kernelsize+d];
+      sumy += local_memory[dx*3+1] * divby[kernelsize+d];
+      sumz += local_memory[dx*3+2] * divby[kernelsize+d];
+    //  printf("Divby factor %f kernelsize+d %d\n", divby[kernelsize+d], kernelsize+d);
+ 		}
+
+    out[(y*imagesizex+x)*3+0] = sumx;// / old_divby;
+  	out[(y*imagesizex+x)*3+1] = sumy;// / old_divby;
+  	out[(y*imagesizex+x)*3+2] = sumz;/// old_divby;
+  }
+
+}
+
 __global__ void filter_old(unsigned char *image, unsigned char *out, const unsigned int imagesizex, const unsigned int imagesizey, const int kernelsizex, const int kernelsizey)
 {
   int x = blockIdx.x * blockDim.x + threadIdx.x;
@@ -180,9 +266,98 @@ __global__ void filter_old(unsigned char *image, unsigned char *out, const unsig
   out[(y*imagesizex+x)*3+2] = sumz/divby;
   }
 }
+
+__global__ void filter_median(unsigned char *image, unsigned char *out, const unsigned int imagesizex, const unsigned int imagesizey, const int kernelsizex, const int kernelsizey)
+{
+
+  __shared__ unsigned char local_memory [STATIC_SHARED_MEM_SIZE_PX * 3];
+  // map from blockIdx to pixel position
+  unsigned threadNumInBlock = threadIdx.x + blockDim.x * threadIdx.y;
+  unsigned threadsPerBlock  = blockDim.x * blockDim.y;
+
+  int blockOffsetX = blockIdx.x * blockDim.x;
+  int blockOffsetY = blockIdx.y * blockDim.y;
+
+  int local_index = threadNumInBlock;
+
+  int FILTER_W = (kernelsizex * 2 + 1) + (2 * kernelsizex);
+  int FILTER_H = (kernelsizey * 2 + 1) + (2 * kernelsizey);
+
+  unsigned LOCAL_MEM_SIZE = FILTER_W * FILTER_H;
+
+  while(local_index < LOCAL_MEM_SIZE)
+  {
+    int global_x = min(max(0, blockOffsetX - kernelsizex + local_index % FILTER_W), imagesizex - 1);
+	  int global_y = min(max(0, blockOffsetY - kernelsizey + local_index / FILTER_W), imagesizey - 1);
+    unsigned global_index = global_y * imagesizex + global_x;
+
+    local_memory[3*local_index+0] = image[3*global_index+0];
+    local_memory[3*local_index+1] = image[3*global_index+1];
+    local_memory[3*local_index+2] = image[3*global_index+2];
+
+    local_index += threadsPerBlock;
+  }
+  __syncthreads();
+  unsigned x = blockIdx.x * blockDim.x + threadIdx.x;
+  unsigned y = blockIdx.y * blockDim.y + threadIdx.y;
+
+  if (x < imagesizex && y < imagesizey) // If inside image
+  {
+    int dy, dx;
+    unsigned int sumx[256] = {0}, sumy[256] = {0}, sumz[256] = {0};
+
+    unsigned local_x = threadIdx.x;
+    unsigned local_y = threadIdx.y;
+
+   	for(dy=-kernelsizey;dy<=kernelsizey;dy++)
+   		for(dx=-kernelsizex;dx<=kernelsizex;dx++)
+   		{
+        unsigned yy = kernelsizey + local_y + dy;
+  			unsigned xx = kernelsizex + local_x + dx;
+        unsigned mem_access = yy * FILTER_W + xx;
+
+        ++sumx[local_memory[mem_access*3+0]];
+        ++sumy[local_memory[mem_access*3+1]];
+        ++sumz[local_memory[mem_access*3+2]];
+   		}
+
+    unsigned num_x = 0;
+    unsigned num_y = 0;
+    unsigned num_z = 0;
+    int index_x = 0;
+    int index_y = 0;
+    int index_z = 0;
+
+    unsigned half_filter_size = ((kernelsizex*2+1) * (kernelsizey*2+1))/2;
+    while(num_x < half_filter_size)
+    {
+      // printf("Thread: %d Num_x %d index_x %d\n",y* blockDim.x + x, num_x, index_x);
+      num_x += sumx[index_x++];
+    }
+    while(num_y < half_filter_size)
+    {
+      num_y += sumy[index_y++];
+    }
+    while(num_z < half_filter_size)
+    {
+      num_z += sumz[index_z++];
+    }
+
+    // while(num_x < half_filter_size && num_y < half_filter_size && num_z < half_filter_size)
+    // {
+    //   num_x += sumx[index_x++];
+    //   num_y += sumy[index_y++];
+    //   num_z += sumz[index_z++];
+    // }
+
+    out[(y*imagesizex+x)*3+0] = index_x;
+  	out[(y*imagesizex+x)*3+1] = index_y;
+  	out[(y*imagesizex+x)*3+2] = index_z;
+  }
+}
 // Global variables for image data
 
-unsigned char *image, *pixels, *dev_bitmap, *dev_input;
+unsigned char *image, *pixels, *dev_bitmap, *dev_input, *dev_bitmap_tmp;
 unsigned int imagesizey, imagesizex; // Image size
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -199,41 +374,75 @@ void computeImages(int kernelsizex, int kernelsizey)
 	pixels = (unsigned char *) malloc(imagesizex*imagesizey*3);
 	cudaMalloc( (void**)&dev_input, imagesizex*imagesizey*3);
 	cudaMemcpy( dev_input, image, imagesizey*imagesizex*3, cudaMemcpyHostToDevice );
+
+  cudaMalloc( (void**)&dev_bitmap_tmp, imagesizex*imagesizey*3);
 	cudaMalloc( (void**)&dev_bitmap, imagesizex*imagesizey*3);
 
-	dim3 grid(imagesizex,imagesizey);
-  ResetMilli();
-	filter_old<<<grid,1>>>(dev_input, dev_bitmap, imagesizex, imagesizey, kernelsizex, kernelsizey); // Awful load balance
-  cudaThreadSynchronize();
-  printf("[Un-Optimized] %f\n", GetSeconds());
 
-  dim3 gridsize(imagesizex/(kernelsizex*2+1), imagesizey/(kernelsizey*2+1));
+  /* Given filter */
+	// dim3 grid(imagesizex,imagesizey);
+  // ResetMilli();
+	// filter_old<<<grid,1>>>(dev_input, dev_bitmap, imagesizex, imagesizey, kernelsizex, kernelsizey); // Awful load balance
+  // cudaThreadSynchronize();
+  // printf("[Un-Optimized] %f\n", GetSeconds());
+
+  /* Our first block filter */
+  dim3 gridsize(imagesizex/(kernelsizex*2+1)+1, imagesizey/(kernelsizey*2+1)+1);
   dim3 blocksize(kernelsizex*2+1, kernelsizey*2+1);
   ResetMilli();
   filter_1<<<gridsize, blocksize>>>(dev_input, dev_bitmap, imagesizex, imagesizey, kernelsizex, kernelsizey); // Nice load balance
   cudaThreadSynchronize();
   printf("[Optimized]    %f\n", GetSeconds());
 
-  dim3 gridsize_sep(imagesizex/(kernelsizex*2+1), imagesizey);
-  dim3 blocksize_sep(kernelsizex*2+1);
-  // unsigned char *tmp = dev_input;
+  /* Separable filter */
+  dim3 gridsize_sep_horizontal(imagesizex/(kernelsizex*2+1) + 1, imagesizey);
+  dim3 gridsize_sep_vertical(imagesizex, imagesizey/(kernelsizex*2+1) + 1);
+
+  dim3 threads_horizontal(kernelsizex*2+1);
+  dim3 threads_vertical(1, kernelsizex*2+1);
+  // // unsigned char *tmp = dev_input;
   ResetMilli();
-  /* X-direction */
-  filter_separable<<<gridsize_sep, blocksize_sep>>>(dev_input, dev_bitmap, imagesizex, imagesizey, kernelsizex, 1);
+  // /* X-direction */
+  filter_separable<<<gridsize_sep_horizontal, threads_horizontal>>>(dev_input, dev_bitmap_tmp, imagesizex, imagesizey, kernelsizex, true);
   cudaThreadSynchronize();
   /* Y-direction */
-  // dim3 gridsize_sep(imagesizex, imagesizey/(kernelsizex*2+1));
-  // filter_separable<<<gridsize_sep, blocksize_sep>>>(tmp, dev_bitmap, imagesizex, imagesizey, kernelsizex, 1);
+  filter_separable<<<gridsize_sep_vertical,   threads_vertical>>>(dev_bitmap_tmp, dev_bitmap, imagesizex, imagesizey, kernelsizex, false);
   cudaThreadSynchronize();
   printf("[Separable]    %f\n", GetSeconds());
 
+
+  /* Gaussian filter i.e. 2 x 5*1 filters run horizontal and vertical with gaussian weights */
+  if(kernelsizex == 2)
+  {
+    ResetMilli();
+    /* X-direction */
+    filter_gaussian<<<gridsize_sep_horizontal, threads_horizontal>>>(dev_input, dev_bitmap_tmp, imagesizex, imagesizey, kernelsizex, true);
+    cudaThreadSynchronize();
+    /* Y-direction */
+    filter_gaussian<<<gridsize_sep_vertical,  threads_vertical>>>(dev_bitmap_tmp, dev_bitmap, imagesizex, imagesizey, kernelsizex, false);
+    cudaThreadSynchronize();
+
+    printf("[Gaussian]    %f\n", GetSeconds());
+  }
+
+
+
+  /* Median filter */
+  ResetMilli();
+  filter_median<<<gridsize, blocksize>>>(dev_input, dev_bitmap, imagesizex, imagesizey, kernelsizex, kernelsizey); // Nice load balance
+  cudaThreadSynchronize();
+  printf("[Median]    %f\n", GetSeconds());
+
 //	Check for errors!
   cudaError_t err = cudaGetLastError();
   if (err != cudaSuccess)
       printf("Error: %s\n", cudaGetErrorString(err));
 	cudaMemcpy( pixels, dev_bitmap, imagesizey*imagesizex*3, cudaMemcpyDeviceToHost );
+
+  cudaFree( dev_bitmap_tmp );
 	cudaFree( dev_bitmap );
 	cudaFree( dev_input );
+
 }
 
 // Display images
@@ -279,7 +488,7 @@ int main( int argc, char** argv)
 	glutDisplayFunc(Draw);
 
 
-  for(int i = 2; i < 10; i++)
+  for(int i = 2; i < 5; i++)
   {
      printf("Frame is %i x %i\n", i*2+1, i*2+1);
 	   computeImages(i, i);