日志

[原创]基于CUDA的图像亮度直方图统计

已有 2151 次阅读2009-8-7 23:20 |个人分类:高性能计算|系统分类:兴趣爱好| 直方图, 统计

算法：

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 ++;
    }

对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, 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;
}