title: Convolution using CUDA with Sobel Filter
body: <p>So I have implemented (or at least tried to) a Sobel filter in CUDA and my code is as follows. When I execute this file, i get half of a correctly Sobel filtered image and the other half is black. I am not able to upload the pictures because they are in .pgm format. So what the code does is read in a grayscale image which is in .pgm format and convolves Sobel filter mask with it using shared memory concept. I used a 1024 by 1024 .pgm image as input and it return a Sobel filtered image with edges that has half of it cut horizontally so it is black in the bottom half. Can someone please help me out here. Also, I am a little curious about the code, I do not really understand what the second batch loading does, so can you please explain that too.</p>
  
  <p><strong>sobel.cu</strong></p>
  
  <pre><code>/* sobel.cu */
  
  #include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;float.h&gt;
  #include &lt;time.h&gt;
  #include "mypgm.h"
  
  #define Mask_width  3
  #define Mask_radius Mask_width/2
  #define TILE_WIDTH 16
  #define w (TILE_WIDTH + Mask_width - 1)
  #define clamp(x) (min(max((x), 0.0), 1.0))
  
  
  __global__ void convolution(float *I, const float* __restrict__ M, float *P, int width, int height) {
      __shared__ float N_ds[w][w];
      int k;
  
      // First batch loading
      int dest = threadIdx.y * TILE_WIDTH + threadIdx.x,
          destY = dest / w, destX = dest % w,
          srcY = blockIdx.y * TILE_WIDTH + destY - Mask_radius,
          srcX = blockIdx.x * TILE_WIDTH + destX - Mask_radius,
          src = srcY * width + srcX;
      if (srcY &gt;= 0 &amp;&amp; srcY &lt; height &amp;&amp; srcX &gt;= 0 &amp;&amp; srcX &lt; width)
          N_ds[destY][destX] = I[src];
      else
          N_ds[destY][destX] = 0;
      for (int iter = 1; iter &lt;= (w*w) / (TILE_WIDTH*TILE_WIDTH); iter++)
      {
          // Second batch loading
          dest = threadIdx.y * TILE_WIDTH + threadIdx.x + TILE_WIDTH * TILE_WIDTH;
          destY = dest / w, destX = dest % w;
          srcY = blockIdx.y * TILE_WIDTH + destY - Mask_radius;
          srcX = blockIdx.x * TILE_WIDTH + destX - Mask_radius;
          src = srcY * width + srcX;
          if (destY &lt; w) {
              if (srcY &gt;= 0 &amp;&amp; srcY &lt; height &amp;&amp; srcX &gt;= 0 &amp;&amp; srcX &lt; width)
                  N_ds[destY][destX] = I[src];
              else
                  N_ds[destY][destX] = 0;
          }
      }
      __syncthreads();
  
      float accum = 0;
      int y, x;
      for (y = 0; y &lt; Mask_width; y++)
          for (x = 0; x &lt; Mask_width; x++)
              accum += N_ds[threadIdx.y + y][threadIdx.x + x] * M[y * Mask_width + x];
  
      y = blockIdx.y * TILE_WIDTH + threadIdx.y;
      x = blockIdx.x * TILE_WIDTH + threadIdx.x;
      if (y &lt; height &amp;&amp; x &lt; width)
          P[y * width + x] = accum;
  
      __syncthreads();
  
  }
  
  void sobel_filtering()
  /* Spatial filtering of image data */
  /* Sobel filter (horizontal differentiation */
  /* Input: image1[y][x] ---- Outout: image2[y][x] */
  
  {
      /* Definition of Sobel filter in horizontal direction */
      float weight[3][3] = { { -1,  0,  1 },
                { -2,  0,  2 },
                { -1,  0,  1 } };
      float pixel_value;
  
      int x, y, i, j;  /* Loop variable */
      float * deviceInputImageData;
      float * deviceOutputImageData;
      float * deviceMaskData;
  
      cudaMalloc((void **)&amp;deviceInputImageData, x_size1 * y_size1 * sizeof(float));
      cudaMalloc((void **)&amp;deviceOutputImageData, x_size1 * y_size1 * sizeof(float));
      cudaMalloc((void **)&amp;deviceMaskData, 3 * 3 * sizeof(float));
  
      cudaMemcpy(deviceInputImageData, image1, x_size1 * y_size1 * sizeof(float), cudaMemcpyHostToDevice);
      cudaMemcpy(deviceMaskData, weight, 3 * 3 * sizeof(float), cudaMemcpyHostToDevice);
  
      /* Maximum values calculation after filtering*/
      printf("Now, filtering of input image is performed\n\n");
      x_size2 = x_size1;
      y_size2 = y_size1;
      for (y = 0; y &lt; y_size2; y++) {
          for (x = 0; x &lt; x_size2; x++) {
              image2[y][x] = 0;
          }
      }
  
      dim3 dimGrid(ceil((float)x_size1 / TILE_WIDTH), ceil((float)y_size1 / TILE_WIDTH));
      dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);
      convolution&lt;&lt;&lt;dimGrid, dimBlock&gt;&gt;&gt;(deviceInputImageData, deviceMaskData, deviceOutputImageData, x_size1, y_size1);
  
  
      cudaMemcpy(image2,
          deviceOutputImageData,
          x_size2 * y_size2 * sizeof(float),
          cudaMemcpyDeviceToHost);
  
      cudaFree(deviceInputImageData);
      cudaFree(deviceOutputImageData);
      cudaFree(deviceMaskData);
  
  }
  
  
  int main()
  {
      load_image_data();   /* Input of image1 */
  
      clock_t begin = clock();
      sobel_filtering();   /* Sobel filter is applied to image1 */
      clock_t end = clock();
      double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
      printf("\n\nTiming result of multiplication of matrix-vector: %f\n", time_spent);
      save_image_data();   /* Output of image2 */
      return 0;
  }
  </code></pre>
  
  <p><strong>mypgm.h</strong></p>
  
  <pre><code>/* pgm file IO headerfile ------ mypgm.h */
  
  /* Constant declaration */
  
  //#define MAX_IMAGESIZE  1024
  
  #define MAX_IMAGEWIDTH  3840
  #define MAX_IMAGEHEIGHT  2160
  #define MAX_BRIGHTNESS  255 /* Maximum gray level */
  #define GRAYLEVEL       256 /* No. of gray levels */
  #define MAX_FILENAME    256 /* Filename length limit */
  #define MAX_BUFFERSIZE  256
  
  /* Global constant declaration */
  /* Image storage arrays */
  float image1[MAX_IMAGEWIDTH][MAX_IMAGEHEIGHT],
  image2[MAX_IMAGEWIDTH][MAX_IMAGEHEIGHT];
  int x_size1, y_size1, /* width &amp; height of image1*/
  x_size2, y_size2; /* width &amp; height of image2 */
  
  /* Prototype declaration of functions */
  void load_image_data( ); /* image input */
  void save_image_data( ); /* image output*/
  void load_image_file(char *); /* image input */
  void save_image_file(char *); /* image output*/
  
  
  /* Main body of functions */
  
  void load_image_data()
  /* Input of header &amp; body information of pgm file */
  /* for image1[ ][ ]，x_size1，y_size1 */
  {
      char file_name[MAX_FILENAME];
      char buffer[MAX_BUFFERSIZE];
      FILE *fp; /* File pointer */
      int max_gray; /* Maximum gray level */
      int x, y; /* Loop variable */
  
      /* Input file open */
      printf("\n-----------------------------------------------------\n");
      printf("Monochromatic image file input routine \n");
      printf("-----------------------------------------------------\n\n");
      printf("     Only pgm binary file is acceptable\n\n");
      printf("Name of input image file? (*.pgm) : ");
      scanf("%s", file_name);
      fp = fopen(file_name, "rb");
      if (NULL == fp) {
          printf("     The file doesn't exist!\n\n");
          exit(1);
      }
      /* Check of file-type ---P5 */
      fgets(buffer, MAX_BUFFERSIZE, fp);
      if (buffer[0] != 'P' || buffer[1] != '5') {
          printf("     Mistaken file format, not P5!\n\n");
          exit(1);
      }
      /* input of x_size1, y_size1 */
      x_size1 = 0;
      y_size1 = 0;
      while (x_size1 == 0 || y_size1 == 0) {
          fgets(buffer, MAX_BUFFERSIZE, fp);
          if (buffer[0] != '#') {
              sscanf(buffer, "%d %d", &amp;x_size1, &amp;y_size1);
          }
      }
      /* input of max_gray */
      max_gray = 0;
      while (max_gray == 0) {
          fgets(buffer, MAX_BUFFERSIZE, fp);
          if (buffer[0] != '#') {
              sscanf(buffer, "%d", &amp;max_gray);
          }
      }
      /* Display of parameters */
      printf("\n     Image width = %d, Image height = %d\n", x_size1, y_size1);
      printf("     Maximum gray level = %d\n\n", max_gray);
      if (x_size1 &gt; MAX_IMAGEWIDTH || y_size1 &gt; MAX_IMAGEHEIGHT) {
          printf("     Image size exceeds %d x %d\n\n",
              MAX_IMAGEWIDTH, MAX_IMAGEHEIGHT);
          printf("     Please use smaller images!\n\n");
          exit(1);
      }
      if (max_gray != MAX_BRIGHTNESS) {
          printf("     Invalid value of maximum gray level!\n\n");
          exit(1);
      }
      /* Input of image data*/
      for (y = 0; y &lt; y_size1; y++) {
          for (x = 0; x &lt; x_size1; x++) {
              image1[y][x] = (unsigned char)fgetc(fp);
          }
      }
      printf("-----Image data input OK-----\n\n");
      printf("-----------------------------------------------------\n\n");
      fclose(fp);
  }
  
  void save_image_data()
  /* Output of image2[ ][ ], x_size2, y_size2 in pgm format*/
  
  {
      char file_name[MAX_FILENAME];
      FILE *fp; /* File pointer */
      int x, y; /* Loop variable */
  
      /* Output file open */
      printf("-----------------------------------------------------\n");
      printf("Monochromatic image file output routine\n");
      printf("-----------------------------------------------------\n\n");
      printf("Name of output image file? (*.pgm) : ");
      scanf("%s", file_name);
      fp = fopen(file_name, "wb");
      /* output of pgm file header information */
      fputs("P5\n", fp);
      fputs("# Created by Image Processing\n", fp);
      fprintf(fp, "%d %d\n", x_size2, y_size2);
      fprintf(fp, "%d\n", MAX_BRIGHTNESS);
      /* Output of image data */
      for (y = 0; y &lt; y_size2; y++) {
          for (x = 0; x &lt; x_size2; x++) {
              fputc(image2[y][x], fp);
          }
      }
      printf("\n-----Image data output OK-----\n\n");
      printf("-----------------------------------------------------\n\n");
      fclose(fp);
  }
  
  void load_image_file(char *filename)
  /* Input of header &amp; body information of pgm file */
  /* for image1[ ][ ]，x_size1，y_size1 */
  {
      char buffer[MAX_BUFFERSIZE];
      FILE *fp; /* File pointer */
      int max_gray; /* Maximum gray level */
      int x, y; /* Loop variable */
  
      /* Input file open */
      fp = fopen(filename, "rb");
      if (NULL == fp) {
          printf("     The file doesn't exist!\n\n");
          exit(1);
      }
      /* Check of file-type ---P5 */
      fgets(buffer, MAX_BUFFERSIZE, fp);
      if (buffer[0] != 'P' || buffer[1] != '5') {
          printf("     Mistaken file format, not P5!\n\n");
          exit(1);
      }
      /* input of x_size1, y_size1 */
      x_size1 = 0;
      y_size1 = 0;
      while (x_size1 == 0 || y_size1 == 0) {
          fgets(buffer, MAX_BUFFERSIZE, fp);
          if (buffer[0] != '#') {
              sscanf(buffer, "%d %d", &amp;x_size1, &amp;y_size1);
          }
      }
      /* input of max_gray */
      max_gray = 0;
      while (max_gray == 0) {
          fgets(buffer, MAX_BUFFERSIZE, fp);
          if (buffer[0] != '#') {
              sscanf(buffer, "%d", &amp;max_gray);
          }
      }
      if (x_size1 &gt; MAX_IMAGEWIDTH || y_size1 &gt; MAX_IMAGEHEIGHT) {
          printf("     Image size exceeds %d x %d\n\n",
              MAX_IMAGEWIDTH, MAX_IMAGEHEIGHT);
          printf("     Please use smaller images!\n\n");
          exit(1);
      }
      if (max_gray != MAX_BRIGHTNESS) {
          printf("     Invalid value of maximum gray level!\n\n");
          exit(1);
      }
      /* Input of image data*/
      for (y = 0; y &lt; y_size1; y++) {
          for (x = 0; x &lt; x_size1; x++) {
              image1[y][x] = (float)fgetc(fp);
          }
      }
      fclose(fp);
  }
  
  void save_image_file(char *filename)
  /* Output of image2[ ][ ], x_size2, y_size2 */
  /* into pgm file with header &amp; body information */
  {
      FILE *fp; /* File pointer */
      int x, y; /* Loop variable */
  
      fp = fopen(filename, "wb");
      /* output of pgm file header information */
      fputs("P5\n", fp);
      fputs("# Created by Image Processing\n", fp);
      fprintf(fp, "%d %d\n", x_size2, y_size2);
      fprintf(fp, "%d\n", MAX_BRIGHTNESS);
      /* Output of image data */
      for (y = 0; y &lt; y_size2; y++) {
          for (x = 0; x &lt; x_size2; x++) {
              fputc(image2[y][x], fp);
          }
      }
      fclose(fp);
  }
  </code></pre>