我的书签
添加书签
移除书签4.2 直方图
直方图是用来计数或可视化离散数据的频度(比如，出现次数)，直方图多用于图像处理。例如本节的例子，我们将创建一个直方图，其统计像素值不超过256-bit的图像。图4.1中，左边的图为输入，右边的图为产生好的直方图。例子中我们将使用局部内存，以及在OpenCL内核中使用局部和全局原子操作。
图4.1 从一副256-bit图产生的直方图。直方图中显示了相关像素的频度。
理论上，直方图算法本身很简单。例子中，每个值都由自己对应的位置，直方图的实现如下：
1. int histogtam[HIST_BINS]
3. main(){
4. for (each input value){
5. histogtam[value]++;
6. }
7. }
与第3章的向量相加不同，直方图不那么容易并行，在多线程的状态下期中的自加操作将会导致竞争。不过，虽然效率有些低，但也可以使用原子加操作来完成直方图中的数据的自加操作。下面的伪代码，就完全能用多线程完成直方图的计算。
1. int histogtam[HIST_BINS]
3. createHistogram(){
4. for (each of my value){
5. atomic_add(histogtam[value], 1);
6. }
7. }
9. main(){
10. for (number of threads){
11. spawn_thread(createHistogram);
12. }
13. }
这个实现并不高效，因为其对每个像素都进行了原子操作。更高效的一种方式是使用局部直方图，其只统计该区域内的像素。当一个线程计算完成自己的局部直方图后，就会自动的加到全局直方图中。
1. int histogram[HIST_BINS]
3. createHistogram(){
5. int localHistorgram[HIST_BINS];
7. for (each of my value){
8. localHistorgram[value]++;
9. }
11. for (each bin){
12. atomic_add(histogram[bin], localHistorgram[bin]);
13. }
14. }
16. main(){
17. for (number of threads){
18. spawn_thread(createHistogram);
19. }
20. }
很多情况下，OpenCL的多线程和CPU多线程在原理上是一样的——不过二者的线程粒度有些不同。如将第一种多线程的方式换成OpenCL中的工作项来做，因为每个工作项都要对共享的全局内存使用原子操作，其执行效率也很低。第8章中，我们将了解到，全局变量的访问延迟要远高于寄存器和局部内存。如多线程的实现，如果有太多的线程访问同一个位置，那么这个程序的效率将会大大降低。
不过，我们也不想在每个工作项中备份直方图。GPU工作项将私有数据存储到寄存器，当寄存器被占满，多于的私有变量会存储在全局内存中，这对于性能来说有弊无益。
最好的办法就是在每个工作组中创建一份局部积分图。局部内存中的数据，每个工作组中的所有工作项都可以共享访问。局部内存一般会分布在GPU的片上内存中，其访问速度要比访问全局内存快的多。如同第二种CPU多线程算法，当工作组完成局部积分图时，其会传递给全局内存，并使用原子加操作将对应位置上的数据原子加到全局内存中。不过，这种实现方式也有问题：对局部内存的访问上存在条件竞争。这里需要你对目标设备的架构有所了解。对于很多GPU来说，原子操作访问局部内存的效率很高。在AMD Radeon GPU上，原子单元位于片上暂存式存储器中。因此，局部内存上的原子操作的效率要比全局原子操作的效率高很多。下面的例子中，我们将使用到局部原子操作来生成局部直方图。
OpenCL内核的参考代码如下：
1. #define HIST_BINS 256
3. __kernel
4. void histogram(__global int *data,
5. int numData,
6. __global int *histogram){
8. __local int localHistorgram[HIST_BINS];
10. int lid = get_local_id(0);
11. int gid = get_glaobal_id(0);
13. /* Initialize local histogram to zero */
14. for (int i = lid; i < HIST_BINS; i += get_local_size(0)){
15. localHistorgram[i] = 0;
16. }
18. /* Wait nutil all work-items within
19. * the work-group have completed their stores */
20. barrier(CLK_LOCAL_MEM_FENCE);
22. /* Compute local histogram */
23. for (int i = gid; i < numData; i += get_glaobal_size(0)){
24. atomic_add(&localHistorgram[data[i]], 1);
25. }
27. /* Wait nutil all work-items within
28. * the work-group have completed their stores */
29. barrier(CLK_LOCAL_MEM_FENCE);
31. /* Write the local histogram out to
32. * the global histogram */
33. for (int i = lid; i < HIST_BINS; i += get_glaobal_size(0)){
34. atomic_add(&histogram[i], localHistorgram[i]);
35. }
36. }
代码清单4.1 计算直方图的OpenCL内核代码
代码清单4.1的实现中包含如下5步：
- 初始化局部直方图内的值为0 (第14行)
- 同步工作项，确保相应的数据全部更新完毕 (第23行)
- 计算局部直方图 (第26行)
- 再次同步工作项，确保相应的数据全部更新完毕 (第35行)
- 将局部直方图写入到全局内存中 (第39行)
1,3,5展示了如何在OpenCL中对内存共享区域(全局或局部内存)进行读写。当我们需要工作项需要访问不同的内存位置时，我们可以以工作项的唯一标识ID为基准，然后加上所有工作项的数量作为跨度(例如，工作组内以工作组中工作项的数量，计算对应工作项所要访问的局部内存位置。或以NDRange中的尺寸，访问全局内存)。第1步中，我们以工作组的尺寸为跨距，用来将局部直方图初始化为0。这就允许我们更加灵活的对工作组尺寸进行设置及配置，并且能保证现有模块功能的正确性。第3步中使用同样的方式读取全局内存中的数据，以及第5步中也以相同的方式将局部内存中的数据写出。
第2和第4步使用栅栏对两步间的操作进行同步，其指定的内存栅栏，将同步工作组中的所有工作项。栅栏和内存栅栏将在第7章详细讨论。现在就能确保工作组中的所有工作项都要到达该栅栏处，只要有线程没有达到，已达到的线程就不能执行下面的操作。局部内存栅栏就是用来保证所有工作项都到达栅栏处，以代表局部直方图更新完毕。
为了让全局直方图得到正确的结果，我们也需要对全局积分图进行初始化。可以在数组创建之后，直接使用主机端API clEnqueueFillBuffer()对数据进行初始化。clEnqueueFillBuffer()的参数列表如下:
1. cl_int
2. clEnqueueFillBuffer(
3. cl_command_queue command_queue,
4. cl_mem buffer,
5. const void *pattern,
6. size_t offset,
7. size_t size,
8. cl_uint num_events_in_wait_list,
9. const cl_event *event_wait_list,
10. cl_event *event)
该API类似于C中的memset()函数。buffer参数就是要初始化的数组对象，具体的值由pattern指定。与memset()不同，pattern可以指定为任意的OpenCL支持类型，比如：标量、整型向量或浮点类型。pattern_size用来指定pattern所占空间。size参数用来指定数组内初始化的字节数，其值必须是pattern_size的整数倍。offset参数用来指定数组起始初始化的位置或偏移。
除了初始化直方图数组部分，主机端代码与第3章向量相加的代码很类似。代码清单4.2提供了直方图统计主机端的完整代码。内核代码在代码清单4.1中，将其存成名为histogram.cl以供4.2中的代码使用。代码中有些工具函数在4.6节中会详细介绍。这些函数用于读写BMP文件，这里提供其源码的在线地址(http://booksite.elsevier.com/9780128014141)
1. /* System includes */
2. #include <stdio.h>
3. #include <stdlib.h>
4. #include <sting.h>
6. /* OpenCL includes */
7. #include <CL/cl.h>
9. /* Utility functions */
10. #include "utils.h"
11. #include "bmp_utils.h"
13. static const int HIST_BINS = 256;
15. int main(int argc, char *argv[]){
17. /* Host data */
18. int *hInputImage = NULL;
19. int *hOutputHistogram = NULL;
21. /* Allocate space for the input image and read the
22. * data from disk */
23. int imageRows;
24. int imageCols;
25. hInputImage = readBmp("../../Images/cat.bmp", &imageRows, &imageCols);
26. const int imageElements = imageRows * imageCols;
27. const size_t imageSize = imageElements * sizeof(int);
29. /* Allocate space for the histogram on the host */
30. const int histogramSize = HIST_BINS * sizeof(int);
31. hOutputHistogram = (int *)malloc(histogramSize);
32. if (!hOutputHistogram){ exit(-1); }
34. /* Use this check the output of each API call */
35. cl_int status;
37. /* Get the first platform */
38. cl_platform_id platform;
39. status = clGetPlatformIDs(1, &platform, NULL);
40. check(status);
42. /* Get the first device */
43. cl_device_id device;
44. status = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
45. check(status);
47. /* Create a command-queue and associate it with the device */
48. cl_command_queue cmdQueue;
49. context = clCreateContext(NULL, 1, &device, NULL, NULL, &status);
50. check(status);
52. /* Create a buffer object for the output histogram */
53. cl_mem bufOutputHistogram;
54. bufOutputHistogram = clCreateBuffer(context, CL_MEM_WRITE_ONLY, histogramSize, NULL, &status);
55. check(status);
57. /* Write the input image to the device */
58. status = clEnqueueWriteBuffer(cmdQueue, bufInputImage, CL_TRUE, 0, imageSize, hInputImage, 0, NULL, NULL);
59. check(status);
61. /* Initialize the output histogram with zero */
62. int zero = 0;
63. status = clEnqueueFillBuffer(cmdQueue, bufOutputHistogram, &zero, sizeof(int), 0, histogramSize, 0, NULL, NULL);
64. check(status);
66. /* Create a program with source code */
67. char *programSource = readFile("histogram.cl");
68. size_t prograSourceLen = strlen(programSource);
69. cl_program program = clCreateProgramWithSouce(context, 1, (const char **)&programSource, &prograSourceLen, &status);
70. check(status);
72. /* Build (compile) the program for the device */
73. status = clBuildProgram(program, 1, &device, NULL, NULL, NULL);
74. if (status != CL_SUCCESS){
75. printCompilerError(program, device);
76. exit(-1);
77. }
79. /* Create the kernel */
80. cl_kernel kernel;
81. kernel = clCreateKernel(program, "histogram", &status);
82. check(status);
84. /* Set the kernel arguments */
85. status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufInputImage);
86. status |= clSetKernelArg(kernel, 1, sizeof(int), &imageElements);
87. status |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufOutputHistogram);
89. /* Define the index space and work-group size */
90. size_t globalWorkSize[1];
91. globalWorkSize[0] = 1024;
93. size_t localWorkSize[1];
94. localWorkSize[0] = 64;
96. /* Enqueue the kernel for execution */
97. status = clEnqueueNDRangeKernel(cmdQueue, kernel, 1, NULL, globalWorkSize, localWorkSize, 0, NULL, NULL);
98. check(status);
100. /* Read the output histogram buffer to the host */
101. status = clEnqueuReadBuffer(cmdQueue, bufOutputHistogram, CL_TRUE, 0, histogramSize, hOutputHistogram, 0, NULL, NULL);
102. check(status);
104. /* Free OpenCL resources */
105. clReleaseKernel(kernel);
106. clReleaseProgram(program);
107. clReleaseCommandQueue(cmdQueue);
108. clReleaseMemObject(bufInputImage);
109. clReleaseMemObject(bufOutputHistogram);
110. clReleaseContext(context);
112. /* Free host resource */
113. free(hInputImage);
114. free(hOutputHistogram);
115. free(programSource);
117. return 0;
118. }
代码清单4.2 直方图统计的主机端代码。注意，check(cl_int status)是用来检查之前执行命令的状态是否为CL_SUCCESS。
