深入理解CUDA图像并行处理

在CUDA编程中，一个CUDA Kernel是由众多线程（threds）组成，而这些线程又可以被组织成一个或多个block块。在同一线程块中，线程ID是从0开始连续编号的，可以通过内置变量threadIdx来获取：

// 获取本线程的索引，blockIdx 指的是线程块的索引，blockDim 指的是线程块的大小，threadIdx 指的是本线程块中的线程索引
int tid = blockIdx.x * blockDim.x + threadIdx.x;

以对图像的归一化处理为例，需要对图片中的每一个像素点的三个通道值分别除以255，相比于使用CPU进行串行计算，我们可以使用CUDA核函数创建更多的线程和线程块来充分利用GPU的并行处理能力：

// 计算需要的线程总量（高度 x 宽度）:640*640=409600
int jobs = dst_height * dst_width;
// 一个线程块包含256个线程
int threads = 256;
// 计算线程块的数量（向上取整）
int blocks = ceil(jobs / (float)threads);

// 调用kernel函数
preprocess_kernel<<<blocks, threads>>>(
  img_buffer_device, dst, dst_width, dst_height, jobs); // 函数的参数

这里我们定义每个线程块的线程数量为256，线程块的数量为ceil(jobs / (float)threads)，总的线程总量要大于等图片的像素数量。当启动Kernel函数时，GPU上的每个线程都会执行相同的程序代码，从而实现更高效的并行计算，函数具体实现如下：

// 一个线程处理一个像素点
__global__ void preprocess_kernel(
    uint8_t *src, float *dst, int dst_width,
    int dst_height, int edge)
{
  int tid = blockDim.x * blockIdx.x + threadIdx.x;
  if (tid >= edge)
    return;

  int dx = tid % dst_width; // 计算当前线程对应的目标图像的x坐标
  int dy = tid / dst_width; // 计算当前线程对应的目标图像的y坐标

  // normalization（对原图中(x,y)坐标的像素点3个通道进行归一化）
  float c0 = src[dy * dst_width * 3 + dx * 3 + 0] / 255.0f;
  float c1 = src[dy * dst_width * 3 + dx * 3 + 1] / 255.0f;
  float c2 = src[dy * dst_width * 3 + dx * 3 + 2] / 255.0f;

  // bgr to rgb
  float t = c2;
  c2 = c0;
  c0 = t;

  // rgbrgbrgb to rrrgggbbb
  // NHWC to NCHW
  int area = dst_width * dst_height;
  float *pdst_c0 = dst + dy * dst_width + dx;
  float *pdst_c1 = pdst_c0 + area;
  float *pdst_c2 = pdst_c1 + area;
  *pdst_c0 = c0;
  *pdst_c1 = c1;
  *pdst_c2 = c2;
}

其中tid是本线程的索引，dst_width和dst_height是图像的宽和高，edge是图片的像素数量，每一个线程处理一个像素点。由于线程索引是从0开始计数的，我们要确保tid不能超过图片的像素数量edge：

int tid = blockDim.x * blockIdx.x + threadIdx.x;
  if (tid >= edge)
    return;

由于图像数据以行优先（row-major）顺序连续存储在内存中，每个像素由3个字节表示（BGR）。为了获取每个线程所处理的像素点在内存中的起始位置，我们可以先计算当前线程所对应图像的x和y坐标即dx和dy：

int dx = tid % dst_width; // 计算当前线程对应的目标图像的x坐标
int dy = tid / dst_width; // 计算当前线程对应的目标图像的y坐标

然后获取当前线程所处理的像素点在内存中的起始位置：dy * dst_width * 3 + dx * 3，*3是因为每个像素点有3个通道值，在内存中的排列方式为:BGRBGRBGR...，最后再/255对原图中(x,y)坐标的像素点3个通道值进行归一化：

// normalization
float c0 = src[dy * dst_width * 3 + dx * 3 + 0] / 255.0f;
float c1 = src[dy * dst_width * 3 + dx * 3 + 1] / 255.0f;
float c2 = src[dy * dst_width * 3 + dx * 3 + 2] / 255.0f;

• dy * dst_width * 3：定位到第dy行的起始位置
• dx * 3：在当前行中定位到第dx个像素的起始位置
• + 0, + 1, + 2：分别访问B、G、R三个通道的值除以255

交换变量做BGR到RGB的通道转换：

// bgr to rgb
float t = c2;
c2 = c0;
c0 = t;

目标图像（RGB）像素点在内存中的排列方式为RRR...GGG...BBB，当前像素点R通道的值在目标图像中内存地址为(dst + dy * dst_width + dx)，G通道的值在目标图像中内存地址为(dst + dy * dst_width + dx) + area，加上1个通道的偏移量area，以此类推，完成对图像的通道转换：

// NHWC to NCHW
// rgbrgbrgb to rrrgggbbb
int area = dst_width * dst_height;
float *pdst_c0 = dst + dy * dst_width + dx;
float *pdst_c1 = pdst_c0 + area;
float *pdst_c2 = pdst_c1 + area;
*pdst_c0 = c0;
*pdst_c1 = c1;
*pdst_c2 = c2;