使用CUDA和Sobel滤波器进行卷积

1 回答

• 2

  您只看到部分图像的原因是您的主机图像缓冲区大小与设备图像缓冲区大小不匹配 .

  在主机上，您具有如下定义的图像缓冲区：

  #define MAX_IMAGEWIDTH  3840
  #define MAX_IMAGEHEIGHT  2160
  ...
  float image1[MAX_IMAGEWIDTH][MAX_IMAGEHEIGHT],
  image2[MAX_IMAGEWIDTH][MAX_IMAGEHEIGHT];
  

  然后，您继续加载尺寸为1024x1024的PGM图像 . 然后，您创建大小为1024x1024的设备存储：

  cudaMalloc((void **)&deviceInputImageData, x_size1 * y_size1 * sizeof(float));
  cudaMalloc((void **)&deviceOutputImageData, x_size1 * y_size1 * sizeof(float));
  

  其中 x_size 和 y_size1 这里是由您的PGM加载例程定义的，对于1024x1024图像，每个例程将为1024 .

  然后，当您从主机到设备进行复制时（在设备 - >主机副本上也会发生类似的问题）：

  cudaMemcpy(deviceInputImageData, image1, x_size1 * y_size1 * sizeof(float), cudaMemcpyHostToDevice);
  

  您将从主机缓冲区复制连续的字节到设备缓冲区 . 这意味着将最大为 MAX_IMAGEWIDTH 的整个宽度的每个主机缓冲区行复制到设备 . 但这不是你想要的 . 您只希望将复制到 x_size1 的每个主机缓冲行复制到设备 .

  有几种方法可以解决这个问题 . 我非常确定最简单的方法就是将 MAX_IMAGEWIDTH 和 MAX_IMAGEHEIGHT 设置为您打算使用的图像的实际值 . 当我这样做时，我得到了一个看似合理的过滤结果 .

  由于这限制了您处理单个图像大小，因此更好的方法是在读取PGM标头数据后动态定义主机图像缓冲区的大小 .

  或者您也可以使用涉及 cudaMemcpy2D 的方法，但这似乎是不必要的复杂 .

  关于第二个问题，“第二批加载”的原因是因为第一批加载仅加载共享内存直到threadblock的维度，因此它加载一个16x16“补丁”的共享内存，每个线程一个元素 . 但是我们需要加载完整的共享内存阵列，因此我们必须进行额外的“批量”加载以填充与过滤器宽度和高度相关的晕区 .

  这是一个似乎对我正常工作的修改过的文件，演示了动态分配主机映像缓冲区的方法：

  #include <stdio.h>
  #include <stdlib.h>
  #include <float.h>
  #include <time.h>
  
  #define MAX_IMAGEWIDTH  2048
  #define MAX_IMAGEHEIGHT 2048
  #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, *image2;
  int x_size1, y_size1, /* width & height of image1*/
  x_size2, y_size2; /* width & 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 & 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", &x_size1, &y_size1);
          }
      }
      /* input of max_gray */
      max_gray = 0;
      while (max_gray == 0) {
          fgets(buffer, MAX_BUFFERSIZE, fp);
          if (buffer[0] != '#') {
              sscanf(buffer, "%d", &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 > MAX_IMAGEWIDTH || y_size1 > 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);
      }
      image1 = (float *)malloc(x_size1*y_size1*sizeof(float));
      /* Input of image data*/
      for (y = 0; y < y_size1; y++) {
          for (x = 0; x < x_size1; x++) {
              image1[y*x_size1+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 < y_size2; y++) {
          for (x = 0; x < x_size2; x++) {
              fputc(image2[y*x_size2+x], fp);
          }
      }
      printf("\n-----Image data output OK-----\n\n");
      printf("-----------------------------------------------------\n\n");
      fclose(fp);
  }
  
  #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];
  
      // 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 >= 0 && srcY < height && srcX >= 0 && srcX < width)
          N_ds[destY][destX] = I[src];
      else
          N_ds[destY][destX] = 0;
      for (int iter = 1; iter <= (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 < w) {
              if (srcY >= 0 && srcY < height && srcX >= 0 && srcX < width)
                  N_ds[destY][destX] = I[src];
              else
                  N_ds[destY][destX] = 0;
          }
      }
      __syncthreads();
  
      float accum = 0;
      int y, x;
      for (y = 0; y < Mask_width; y++)
          for (x = 0; x < 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 < height && x < width)
          P[y * width + x] = accum;
  
  }
  
  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 } };
  
      int x, y;  /* Loop variable */
      float * deviceInputImageData;
      float * deviceOutputImageData;
      float * deviceMaskData;
  
      cudaMalloc((void **)&deviceInputImageData, x_size1 * y_size1 * sizeof(float));
      cudaMalloc((void **)&deviceOutputImageData, x_size1 * y_size1 * sizeof(float));
      cudaMalloc((void **)&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;
      image2 = (float *)malloc(x_size2*y_size2*sizeof(float));
      for (y = 0; y < y_size2; y++) {
          for (x = 0; x < x_size2; x++) {
              image2[y*x_size2+x] = 0;
          }
      }
  
      dim3 dimGrid(ceil((float)x_size1 / TILE_WIDTH), ceil((float)y_size1 / TILE_WIDTH));
      dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);
      convolution<<<dimGrid, dimBlock>>>(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;
  }
  

  注意上面的代码（您的PGM加载例程）有（IMO）缺陷，因为它需要在PGM文件的同一行指定x和y大小，但据我所知，这不是P5的要求PGM文件 . 如果您传递一个有效的P5 PGM文件，该文件在文件的不同行上指定了x和y大小，它将挂起 . 我没有尝试解决这个问题，因为它似乎不是你要问的问题 .

  回复于 2024-04-27T06:31:53+08:00