使用 CUDA 进行图像处理

• 当下生活在高清摄像头的时代，这种摄像头能捕获高达1920*1920像素的高解析度画幅。想要实施的处理这么多的数据，往往需要几个TFlops地浮点处理性能，这些要求CPU也无法满足
• 通过在代码中使用CUDA，可以利用GPU提供的强大地计算能力
• CUDA支持多维地Grid和块，因此可以根据图像地尺寸、数据量大小，合理的分配块和线程进行图像处理
• 简单图像处理过程地特定编程模式：

for(int i=0;i<image_height;i++)
{
    for(int j=0;j<image_width;j++)
    {
        //Pixel Processing code for pixel located at(i,j)
    }
}

• 将像素处理映射到CUDA地一批线程上：

int i = blockidx.y * blockDim.y + threadIdx.y
int j = blockidx.x * blockDim.x + threadIdx.x

1. 在GPU上通过CUDA进行直方图统计

• 首先介绍CPU版本的直方图统计，实现如下：

int h_a[1000] = Random values between 0 and 15

//假设图像取值范围在【0-15】，定义数组并初始化
int histogram[16];
for(int i=0;i<16;i++)
{
    histogram[i] = 0;
}
//统计每个值的个数
for(int i=0;i<1000;i++)
{
    histogram[h_a[i]]+=1;
}

• 下面写一个同样功能的GPU代码，我们将使用3种不同的方法写这个代码，前两种方法的内核代码如下：

__global__ void histogram_without_atomic(int* d_b, int* d_a)
{
	int tid = threadIdx.x + blockDim.x * blockIdx.x;
	int item = d_a[tid];
	if (tid < SIZE)
	{
		d_b[item]++;
	}

}

__global__ void histogram_atomic(int* d_b, int* d_a)
{
	int tid = threadIdx.x + blockDim.x * blockIdx.x;
	int item = d_a[tid];
	if (tid < SIZE)
	{
		atomicAdd(&(d_b[item]), 1);
	}
}

• 第一个函数是最简单方式实现的直方图统计，每个线程读取 1 个元素值。使用线程ID作为输入数组的索引获取该元素的数值，然后此值再将对应的d_b结果数组中的索引位置处进行 +1 操作。最后d_b数组应该包含输入数据中0-15之间每个值的频次，这种方式将得出错误的结果，因为对相同的存储器位置将有大量的线程试图同时进行不安全的修改，其运行结果如下：
  
• 第二个函数用原子操作实现统计，避免多线程并行时的资源占用导致的计算异常问题，其计算结果如下：
  
• main函数如下：

int main()
{
	//定义设备变量并分配内存
	int h_a[SIZE];
	for (int i = 0; i < SIZE; i++) {

		h_a[i] = i % NUM_BIN;
	}
	int h_b[NUM_BIN];
	for (int i = 0; i < NUM_BIN; i++) {
		h_b[i] = 0;
	}

	// 声明GPU指针变量
	int* d_a;
	int* d_b;

	// 分配GPU变量内存
	cudaMalloc((void**)&d_a, SIZE * sizeof(int));
	cudaMalloc((void**)&d_b, NUM_BIN * sizeof(int));

	// transfer the arrays to the GPU
	cudaMemcpy(d_a, h_a, SIZE * sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, NUM_BIN * sizeof(int), cudaMemcpyHostToDevice);

	// 进行直方图统计
	//histogram_without_atomic << <((SIZE + NUM_BIN - 1) / NUM_BIN), NUM_BIN >> > (d_b, d_a);
	histogram_atomic << <((SIZE+NUM_BIN-1) / NUM_BIN), NUM_BIN >> >(d_b, d_a);

	// copy back the sum from GPU
	cudaMemcpy(h_b, d_b, NUM_BIN * sizeof(int), cudaMemcpyDeviceToHost);
	printf("Histogram using 16 bin without shared Memory is: \n");
	for (int i = 0; i < NUM_BIN; i++) {
		printf("bin %d: count %d\n", i, h_b[i]);
	}

	// free GPU memory allocation
	cudaFree(d_a);
	cudaFree(d_b);
	return 0;
}

• 当我们试图测量使用了原子操作的该代码的性能的时候，你会发现相比CPU的性能，对于很大规模的数组，GPU的实现更慢。这就引入了一个问题：我们真的应当使用CUDA进行直方图统计吗？如果必须能否将这个计算更快些？
• 这两个问题的答案都是：YES 。如果我们在一个块中用共享内存进行直方图统计，最后再将每个块的部分统计结果叠加到全局内存上的最终结果上去。这样就能加速该操作。这是因为整数加法满足交换律。我需要补充的是：只有当原始数据就在GPU的显存上的时候，才应当考虑使用GPU计算，否则完全不应当 cudaMemcpy 过来再计算，因为仅 cudaMemcpy 的时间就将等于或者大于 CPU 计算的时间，用共享内存进行直方图统计的内核函数代码实现如下：

#include <stdio.h>
#include <cuda_runtime.h>

#define SIZE 1000
#define NUM_BIN 256

__global__ void histogram_shared_memory(int* d_b, int* d_a)
{
	int tid = threadIdx.x + blockDim.x * blockIdx.x;
	int offset = blockDim.x * gridDim.x;
	__shared__ int cache[256];
	cache[threadIdx.x] = 0;
	__syncthreads();

	while (tid < SIZE)
	{
		atomicAdd(&(cache[d_a[tid]]), 1);
		tid += offset;
	}
	__syncthreads();
	atomicAdd(&(d_b[threadIdx.x]), cache[threadIdx.x]);
}

• 我们要为当前的每个块都统计一次局部结果，所以需要先将共享内存清空，然后用类似之前的方式在共享内存中进行直方图统计。这种情况下，每个块只会统计部分结果存储在各自的共享内存中，并非像以前那样直接统计为全局内存上的总体结果。
• 本例中，块中256个线程进行共享内存上的256个元素的访问，而原本的代码则在全局内存上的16个元素位置上进行访问。因为共享内存本身要比全局内存具有更高效的并行访问性能，同时将16个统一的竞争访问的位置放宽到了每个共享内存上的256个竞争位置，这两个因素共同缩小了原子操作累计统计的时间。
• 最终还需要进行一次原子操作，将每个块的共享内存上的部分统计结果累加到全局内存上的最终统计结果。因为整数加法满足交换律，我们不需要担心每个块执行的顺序。
• main函数如上一个类似：

int main()
{
	// generate the input array on the host
	int h_a[SIZE];
	for (int i = 0; i < SIZE; i++) {
		//h_a[i] = bit_reverse(i, log2(SIZE));
		h_a[i] = i % NUM_BIN;
	}
	int h_b[NUM_BIN];
	for (int i = 0; i < NUM_BIN; i++) {
		h_b[i] = 0;
	}

	// declare GPU memory pointers
	int* d_a;
	int* d_b;

	// allocate GPU memory
	cudaMalloc((void**)&d_a, SIZE * sizeof(int));
	cudaMalloc((void**)&d_b, NUM_BIN * sizeof(int));

	// transfer the arrays to the GPU
	cudaMemcpy(d_a, h_a, SIZE * sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, NUM_BIN * sizeof(int), cudaMemcpyHostToDevice);

	// launch the kernel
	histogram_shared_memory << <SIZE / 256, 256 >> > (d_b, d_a);

	// copy back the result from GPU
	cudaMemcpy(h_b, d_b, NUM_BIN * sizeof(int), cudaMemcpyDeviceToHost);
	printf("Histogram using 16 bin is: ");
	for (int i = 0; i < NUM_BIN; i++) {
		printf("bin %d: count %d\n", i, h_b[i]);
	}

	// free GPU memory allocation
	cudaFree(d_a);
	cudaFree(d_b);

	return 0;
}

• 执行结果：