OpenCL 學習step by step (7) 灰度圖Histogram計算(1)

histogram翻譯成中文就是直方圖，在計算機圖像處理和視覺技術中，通常用histogram來進行圖像匹配，從而完成track，比如meanshift跟蹤算法中，經常要用到圖像的直方圖。

灰度圖的histogram計算，首先要選擇bin(中文可以稱作槽)的數量，對於灰度圖，像素的範圍通常是[0-255],所以bin的數目就是256，然後我們循環整幅圖像，統計出每種像素值出現的次數，放到對應的bin中。比如bin[0]中放的就是整幅圖像中灰度值為0的像素個數，bin[1]中放的就是整幅圖像中灰度值為1的像素個數……

下面的直方圖就是灰度圖lenna對應的直方圖。

灰度圖直方圖的cpu計算特別簡單，定義一個數組hostBin[256],初始化所有數組元素為0，然後循環整幅圖像，得到直方圖，代碼如下：

//cpu求直方圖

void cpu_histgo()

    {

    int i, j;

    for(i = 0; i < height; ++i)

        {

        for(j = 0; j < width; ++j)

            {

            //printf("data: %d\n", data[i * width + j] );

            hostBin[data[i * width + j]]++;

            //printf("hostbin %d=%d\n", data[i * width + j], hostBin[data[i * width + j]]);

            }

        }

    }

如何使用opencl,來計算灰度圖，就沒有那麼簡單了。我們知道gpu的優勢是並行計算，如何把圖像分塊，來並行計算直方圖，是我們討論的重點。下面是一副512*512圖像的thread，workgroup劃分：

我們設定圖像的寬度是bins的整數倍，即256的倍數，高度是workgroup size(本程序中，設置為128)的倍數，如果圖像高寬不是bins和workgroup size的倍數，則我們通過下面的公式把圖像的寬度和高度變成它們的倍數：

//width是binSize的整數倍，height是groupsize的整數倍

width = (width / binSize ? width / binSize: 1) * binSize;

height = (height / groupSize ? height / groupSize: 1) * groupSize;

則512*512的圖像可以分為8個work group，每個workgroup包括128個thread，每個thread計算256個像素的直方圖，並把結果放到該thread對應的local memroy空間，在kenrel代碼結束前，合併一個workgroup中所有thread的直方圖，生成一個workgroup塊的直方圖，最後在host端，合併8個workgroup塊的直方圖，產生最終的直方圖。

OpenCL的memory對象主要有3個，dataBuffer用來傳入圖像數據，而minDeviceBinBuf大小是workgroup number *256, 即每個workgroup對應一個bin，另外一個kernel函數的第二個參數，它的大小為workgroup size*256, 用於workgroup中的每個線程存放自己256個像素的直方圖結果。

//創建2個OpenCL內存對象

dataBuf = clCreateBuffer(

    context,

    CL_MEM_READ_ONLY,

    sizeof(cl_uchar) * width  * height,

    NULL,

    0);

//該對象存放每個block塊的直方圖結果

midDeviceBinBuf = clCreateBuffer(

    context,

    CL_MEM_WRITE_ONLY,

    sizeof(cl_uint) * binSize * subHistgCnt,

    NULL,

    0);

   …

    status = clSetKernelArg(kernel, 1, groupSize * binSize * sizeof(cl_uchar), NULL); //local memroy size, lds for amd

下面看看kernel代碼是如何計算workgroup塊的直方圖。

__kernel

void histogram256(__global const uchar* data,

                  __local uchar* sharedArray,

                  __global uint* binResult)

{

    size_t localId = get_local_id(0);

    size_t globalId = get_global_id(0);

    size_t groupId = get_group_id(0);

    size_t groupSize = get_local_size(0);

下面這部分代碼初始化每個thread對應的local memory，也就是對應的256個bin中計數清零。sharedArray大小是workgroup size * 256 = 128 * 256

//初始化共享內存

    for(int i = 0; i < BIN_SIZE; ++i)

        sharedArray[localId * BIN_SIZE + i] = 0;

通過barrier設置workgroup中所有thread的同步點，保證所有thread都完成初始化操作。

    barrier(CLK_LOCAL_MEM_FENCE);

下面的代碼，計算thread中，256個像素的直方圖，比如對於workgroup 0中的thread 0它計算的256個像素為綠條的部分像素，注意：每個thread的包含的像素並不是連續的。

    

//計算thread直方圖

    for(int i = 0; i < BIN_SIZE; ++i)

    {

        uint value = (uint)data[groupId * groupSize * BIN_SIZE + i * groupSize + localId];

        sharedArray[localId * BIN_SIZE + value]++;

    }

通過fence，保證每個thread都完成各自的直方圖計算。

    barrier(CLK_LOCAL_MEM_FENCE); 

下面是合併各個thread的直方圖形成整個workgroup像素塊的直方圖，每個thread合併2個bin，比如thread 0，合併bin0和bin128。

   //合併workgroup中所有線程的直方圖，產生workgroup直方圖

    for(int i = 0; i < BIN_SIZE / groupSize; ++i)

    {

        uint binCount = 0;

        for(int j = 0; j < groupSize; ++j)

            binCount += sharedArray[j * BIN_SIZE + i * groupSize + localId];

        binResult[groupId * BIN_SIZE + i * groupSize + localId] = binCount;

    }

}

最終在host端，我們還要把每個workgroup塊的直方圖合併成得到整個圖像的直方圖，主要代碼如下：

// 合併子塊直方圖值

for(i = 0; i < subHistgCnt; ++i)

    {

    for( j = 0; j < binSize; ++j)

        {

        deviceBin[j] += midDeviceBin[i * binSize + j];

        }

    }

完整的代碼請參考：

工程文件gclTutorial7

代碼下載：http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

原文作者：邁克老狼