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
代碼下載：