基于CUDA的图像亮度直方图统计

算法：
1、先计算原始图像每个像素的亮度：u = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b)。

2、用一个256大小的数组统计每个亮度的点的数量。
C++实现方法：

memset(out, 0, sizeof(out)); unsigned long offset; unsigned long p; unsigned char r, g, b, u; offset = 0; for(y = 0; y < h; y ++) for(x = 0; x < w; x ++) { p = ((unsigned long *)in)[offset ++]; b = p & 0xff; g = (p >> 8) & 0xff; r = (p >> 16) & 0xff; // 计算亮度，结果0~255 u = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b); // 统计 out[u] ++; }

对1280×1024大小的图像，在QX6600（2.4GHz）上单线程运行时间为17.5ms。



CUDA实现的方法比较复杂。因为如果每个线程处理一个像素然后累加到统计数组中的话，多个线程同时累加一个地址会造成数据错误，需要用原子操作进行排队，大量的线程排队操作会造成计算单元空闲，无法发挥大部分的计算性能。因此从GPU内部架构考虑，最好的方法是分两步。第一步把整个图像分成很多个部分，每个部分用一个线程块（含很多个线程）独立统计，统计结果放在线程块共有的shared memory中（同样需要用原子操作，但是shared
memory在GPU内部，速度很快，而且每个SM有独立的shared memory，因此这样操作并行度很高，能使计算单元满载运行）；第二步则把各个部分独立统计的结果再累加起来得到总的统计结果。代码比较复杂：

#define THREAD_N 128 #define LOOP_N 64 __global__ void histKernel(unsigned char *in, unsigned long *out) { __shared__ unsigned long smem[256]; // shared memory const unsigned long tid = threadIdx.x; const unsigned long bid = blockIdx.x; unsigned long offset = __umul24((__umul24(bid, THREAD_N) + tid), LOOP_N); // (bid * THREAD_N + tid) * LOOP_N int i; unsigned char r, g, b, u; unsigned long p; smem[tid] = smem[tid + 128] = 0; __syncthreads(); // 每个线程块有THREAD_N(128)个线程，每个线程处理LOOP_N(64)个点，统计结果存储在每个线程块的smem[256]中 for(i = 0; i < LOOP_N; i ++) { p = *(unsigned long *)&in[offset << 2]; b = p & 0xff; g = (p >> 8) & 0xff; r = (p >> 16) & 0xff; offset ++; // 计算亮度(0~255) u = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b); // 用原子操作统计，防止线程同时进行“读-修改-写”操作时造成冲突 atomicAdd((int *)&smem[u], 1); } __syncthreads(); // 线程块统计计算完成后，汇总各线程块的统计结果 // 把结果从smem[256]累加到global memory中 // 128字交替访存，以满足各线程的合并访问要求以及防止shared memory的bank conflict，提高效率 atomicAdd((int *)&out[tid], smem[tid]); atomicAdd((int *)&out[tid + 128], smem[tid + 128]); } extern "C" float histCall(unsigned char *in, unsigned long *out, int w, int h) { unsigned char *device_src = 0; unsigned long *device_dest = 0; cudaMalloc((void **)&device_src, w * h * sizeof(unsigned long)); cudaMalloc((void **)&device_dest, 256 * sizeof(unsigned long)); cudaMemset(device_dest, 0, 256 * sizeof(unsigned long)); unsigned int timer = 0; cutCreateTimer(&timer); cudaMemcpy(device_src, in, w * h * sizeof(unsigned long), cudaMemcpyHostToDevice); cutStartTimer(timer); // 各线程块并行统计，每个线程块处理THREAD_N*LOOP_N个点，然后汇总累加至最终结果 histKernel<<<w * h / THREAD_N / LOOP_N, THREAD_N>>>(device_src, device_dest); cudaThreadSynchronize(); cutStopTimer(timer); cudaMemcpy(out, device_dest, 256 * sizeof(unsigned long), cudaMemcpyDeviceToHost); float ms = cutGetTimerValue(timer); cutDeleteTimer(timer); cudaFree(device_dest); cudaFree(device_src); return ms; }

对1280×1024大小的图像，在Geforce GTX 285上计算时间为0.54ms，比QX6600上单线程提速约30倍。



使用CUDA编程时要合理利用GPU内部的shared memory以及熟悉其优化规则，以达到最高效率。