title: CUDA - median filter implementation does not produce desired results
body: <p>I have been trying to implement the algorithm for the median filter presented in the Wiki article: <a href="http://en.wikipedia.org/wiki/Median_filter#2D_median_filter_pseudo_code" rel="nofollow">http://en.wikipedia.org/wiki/Median_filter#2D_median_filter_pseudo_code</a></p>
  
  <p>To the best of my knowledge, I know that what I have implemented is right. However, when I view the results, I cannot seem to get an output which is similar to the output produced by the <code>median blur</code> function in OpenCV. At the present moment, I am not concerned about speeding up my code by using shared memory or texture memory. I would just like to get things working first. The size of my input image is <code>1024 x 256</code> pixels. </p>
  
  <p><strong>What am I doing wrong? Is there a thread leak in my code?</strong> I know that I should use shared memory to prevent the global read because currently, I am reading data from the global memory a lot. </p>
  
  <p><a href="http://snag.gy/OkXzP.jpg" rel="nofollow">http://snag.gy/OkXzP.jpg</a> -- the first image is the input, second image is my algorithm result, and the third is the openCV <code>medianblur</code> function result. Ideally, I would like my algorithm to output the same result as the <code>medianblur</code> function. </p>
  
  <p>This is all the code that I have written:</p>
  
  <p><strong>kernel implementation</strong></p>
  
  <pre><code>#include "cuda.h"
  #include "cuda_runtime_api.h"
  #include "device_launch_parameters.h"
  #include "device_functions.h"
  #include "highgui.h"
  //#include "opencv2/core/imgproc.hpp"
  //#include "opencv2/core/gpu.hpp"
  #include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;string.h&gt;
  #include &lt;math.h&gt;
  
  // includes, project
  
  #include "cufft.h"
  #include "cublas_v2.h"
  #include "CUDA_wrapper.h"   // contains only func_prototype for function take_input()
  
  
  // define the threads and grids for CUDA
  #define BLOCK_ROWS 32
  #define BLOCK_COLS 16
  
  // define kernel dimensions
  #define KERNEL_DIMENSION 3
  #define MEDIAN_DIMENSION 3
  #define MEDIAN_LENGTH 9
  
  // this is the error checking part for CUDA
  #define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
  inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
  {
     if (code != cudaSuccess) 
     {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
        if (abort) exit(code);
     }
  }
  
  // create two vars for the rows and cols of the image
      int d_imgRows;
      int d_imgCols; 
  
  __global__ void FilterKernel (unsigned short *d_input_img, unsigned short *d_output_img, int d_iRows, int d_iCols)
  
  {  
      unsigned short window[BLOCK_ROWS*BLOCK_COLS][KERNEL_DIMENSION*KERNEL_DIMENSION];
  
      unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;
      unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;   
  
      unsigned int tid = threadIdx.y*blockDim.y+threadIdx.x;
  
      if(x&gt;d_iCols || y&gt;d_iRows)
          return;
  
      window[tid][0]= (y==0||x==0) ? 0.0f : d_input_img[(y-1)*d_iCols+(x-1)];
      window[tid][1]= (y==0) ? 0.0f : d_input_img[(y-1)*d_iCols+x]; 
      window[tid][2]= (y==0||x==d_iCols-1) ? 0.0f : d_input_img[(y-1)*d_iCols+(x+1)];
      window[tid][3]= (x==0) ? 0.0f : d_input_img[y*d_iCols+(x-1)];
      window[tid][4]= d_input_img[y*d_iCols+x];
      window[tid][5]= (x==d_iCols-1) ? 0.0f : d_input_img[y*d_iCols+(x+1)];
      window[tid][6]= (y==d_iRows-1||x==0) ? 0.0f : d_input_img[(y+1)*d_iCols+(x-1)];
      window[tid][7]= (y==d_iRows-1) ? 0.0f : d_input_img[(y+1)*d_iCols+x];
      window[tid][8]= (y==d_iRows-1||x==d_iCols-1) ? 0.0f : d_input_img[(y+1)*d_iCols+(x+1)];
  
     __syncthreads();
  
      // Order elements 
      for (unsigned int j=0; j&lt;9; ++j)
      {
          // Find position of minimum element
          int min=j;
          for (unsigned int l=j+1; l&lt;9; ++l)
              if (window[tid][l] &lt; window[tid][min])
                  min=l;
  
          // Put found minimum element in its place
          const unsigned char temp=window[tid][j];
          window[tid][j]=window[tid][min];
          window[tid][min]=temp;
  
          __syncthreads();
      }
  
      d_output_img[y*d_iCols + x] = (window[tid][4]);
  
  }
  
  void take_input(const cv::Mat&amp; input, const cv::Mat&amp; output)
  {
  
      unsigned short *device_input; 
      unsigned short *device_output;
  
      size_t d_ipimgSize = input.step * input.rows;
      size_t d_opimgSize = output.step * output.rows;
  
      gpuErrchk( cudaMalloc( (void**) &amp;device_input, d_ipimgSize) ); 
      gpuErrchk( cudaMalloc( (void**) &amp;device_output, d_opimgSize) ); 
  
      gpuErrchk( cudaMemcpy(device_input, input.data, d_ipimgSize, cudaMemcpyHostToDevice) );
  
      dim3 Threads(BLOCK_ROWS, BLOCK_COLS);  // 512 threads per block
      dim3 Blocks((input.cols + Threads.x - 1)/Threads.x, (input.rows + Threads.y - 1)/Threads.y);    
  
      //int check = (input.cols + Threads.x - 1)/Threads.x;
      //printf( "blockx %d", check);
  
      FilterKernel &lt;&lt;&lt; Blocks, Threads &gt;&gt;&gt; (device_input, device_output, input.rows, input.cols);
  
      gpuErrchk(cudaDeviceSynchronize());
  
      gpuErrchk( cudaMemcpy(output.data, device_output, d_opimgSize, cudaMemcpyDeviceToHost) );
  
      //printf( "num_rows_cuda %d", num_rows);
      //printf("\n");
  
      gpuErrchk(cudaFree(device_input));
      gpuErrchk(cudaFree(device_output));
  
  }
  </code></pre>
  
  <p><strong>main function</strong> </p>
  
  <pre><code>#pragma once
  #include&lt;iostream&gt;
  #include&lt;opencv2/core/core.hpp&gt;
  #include&lt;opencv2/highgui/highgui.hpp&gt;
  #include&lt;opencv2/imgproc/imgproc.hpp&gt;
  #include&lt;opencv2/gpu/gpu.hpp&gt;
  
  #include &lt;CUDA_wrapper.h&gt;
  
  using std::cout;
  using std::endl;
  
  int main()
  {   
  
      //Read the image from harddisk, into a cv::Mat
      //IplImage *img=cvLoadImage("image.jpg");
      //cv::Mat input(img);
      cv::Mat input = cv::imread("C:/Users/OCT/Documents/Visual Studio 2008/Projects/MedianFilter/MedianFilter/pic1.bmp",CV_LOAD_IMAGE_GRAYSCALE);
  
      //IplImage* input = cvLoadImage("G:/Research/CUDA/Trials/OCTFilter/Debug/pic1.bmp");
      if(input.empty())
      {
          cout&lt;&lt;"Image Not Found"&lt;&lt;endl;
          getchar();
          return -1;
      }
  
      cv::Mat output(input.rows,input.cols,CV_8UC1);
  
      // store the different details of the input image like img_data, rows, cols in variables 
      int Rows = input.rows;
      int Cols = input.cols;
      unsigned char* Data = input.data;
  
      cout&lt;&lt;"image rows "&lt;&lt;Rows&lt;&lt;endl;
      cout&lt;&lt;"image cols "&lt;&lt;Cols&lt;&lt;endl;
      cout&lt;&lt;"\n"&lt;&lt;endl;
      cout&lt;&lt;"data "&lt;&lt;(int)Data&lt;&lt;endl;
      cv::waitKey(0);
  
      // call the device function to take the image as input
      take_input(input, output);
  
      cv::Mat dest; 
  
      medianBlur ( input, dest, 3 );
  
      //Show the input and output
      cv::imshow("Input",input);
      cv::imshow("Output",output);
      cv::imshow("Median blur",dest);
  
      //Wait for key press
      cv::waitKey();
  }
  </code></pre>