CUDA的骨骼化加速

0. 简介

CUDA作为并行加速的利器，目前已经被越来越多的利用在图像处理领域中，很多时候我们需要做耗时的图像工作时，使用GPU加速是一个很好地解决方法。充分利用CUDA架构特有的常量存储器和共享存储器对普通并行算法进行改进.讨论了如何根据程序和显卡设备的固有属性来分配线程以达到最高的GPU占用率，从而得到最优的加速效果。

1. CUDA加速

在基于GPU的并行系统中, CPU和GPU各司其职, CPU负责复杂的流控制等需要串行处理的部分, 而密集型数据的并行计算部分则交由GPU完成.有别于原始的串行算法, 本文中CPU只负责简单的体素点划分, 而复杂的势场力计算则交由GPU完成.计算时将各个内部点和边界点平均分配给每个线程, 多线程并行执行.具体的程序流程见下图：

2. Opencv 骨骼化

在使用 Opencv时候我们可以直接调用cv::Mat完成骨骼化的工作。但是我们可以看到这部分内部存在for循环中，非常耗时，所以这部分可以尝试使用CUDA来完成加速。

/**
 * 执行一次细化迭代， 通常不会直接从代码中调用此函数。
 *
 * @param  im   范围为0-1的二值图像
 * @param  iter  0=even, 1=odd
 */
void thinningIteration(cv::Mat& im, int iter, int& diff)
{
    cv::Mat marker = cv::Mat::zeros(im.size(), CV_8UC1);

    for (int i = 1; i < im.rows-1; i++)
    {
        for (int j = 1; j < im.cols-1; j++)
        {
            if(im.at<uchar>(i,j) == 1) {
                uchar p2 = im.at<uchar>(i-1, j);
                uchar p3 = im.at<uchar>(i-1, j+1);
                uchar p4 = im.at<uchar>(i, j+1);
                uchar p5 = im.at<uchar>(i+1, j+1);
                uchar p6 = im.at<uchar>(i+1, j);
                uchar p7 = im.at<uchar>(i+1, j-1);
                uchar p8 = im.at<uchar>(i, j-1);
                uchar p9 = im.at<uchar>(i-1, j-1);

                int A  = (p2 == 0 && p3 == 1) + (p3 == 0 && p4 == 1) + 
                         (p4 == 0 && p5 == 1) + (p5 == 0 && p6 == 1) + 
                         (p6 == 0 && p7 == 1) + (p7 == 0 && p8 == 1) +
                         (p8 == 0 && p9 == 1) + (p9 == 0 && p2 == 1);
                int B  = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
                int m1 = iter == 0 ? (p2 * p4 * p6) : (p2 * p4 * p8);
                int m2 = iter == 0 ? (p4 * p6 * p8) : (p2 * p6 * p8);

                if (A == 1 && (B >= 2 && B <= 6) && m1 == 0 && m2 == 0){
                    marker.at<uchar>(i,j) = 1;
                    diff = 1;
                }
            }
        }
    }

    im &= ~marker;
}

/**
 * 用于细化给定二值图像的函数
 *
 * @param  im  范围为0-255的二进制图像
 */
void thinning(Mat& inimg, Mat& outimg)
{
    int diff = 0;
    inimg.copyTo(outimg);
    do {
        // cout << "diff = " << diff << endl;
        diff = 0;
        thinningIteration(outimg, 0, diff);
        thinningIteration(outimg, 1, diff);
    } 
    while (diff != 0); 
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

3. CUDA加速骨骼化

// Thinning.h
// 实现 Zhang-sune 细化算法

#ifndef __THINNING_H__
#define __THINNING_H__

#include "Image.h"
#include "ErrorCode.h"

class Thinning {

protected:

    // 成员变量：highPixel（高像素）
    // 图像内高像素的像素值，可由用户定义。
    unsigned char highPixel;
    
    // 成员变量：lowPixel（低像素）
    // 图像内低像素的像素值，可由用户定义。
    unsigned char lowPixel;

    // 成员变量：imgCon（图像与坐标集之间的转化器）
    // 当参数为坐标集时，实现坐标集与图像的相互转化。
    // ImgConvert imgCon;
    
public:
    // 宏：DEF_BLOCK_X 和 DEF_BLOCK_Y
    // 定义了默认的线程块尺寸。
    unsigned int DEF_BLOCK_X;
    unsigned int DEF_BLOCK_Y;
    // 构造函数：Thinning
    // 无参数版本的构造函数，所有的成员变量皆初始化为默认值。
    __host__ __device__
    Thinning()
    {
        this->highPixel = 255;  // 高像素值默认为 255。
        this->lowPixel = 0;     // 低像素值默认为 0。
        DEF_BLOCK_X = 32;
        DEF_BLOCK_Y = 8; 
    }

    __host__ int thinZS(Image *inimg, Image *outimg);
};

#endif

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

详细的代码可见链接，这里来介绍一下下面的核心部分，对应了骨骼化操作。下面的_thinZS1Ker方法对应了上面Opencv中的thinningIteration算法，只是将传统的Image图像封装成了ImageCuda 。

// ...............
// 直接并行化
// 线程数，处理多少个点有多少线程数
__host__ int Thinning::thinZS(Image *inimg, Image *outimg)
{
    // 局部变量，错误码。
     int errcode;  
     cudaError_t cudaerrcode; 

     // 检查输入图像，输出图像是否为空。
     if (inimg == NULL || outimg == NULL)
         return NULL_POINTER;

     // 声明所有中间变量并初始化为空。
     Image *tempimg = NULL;
     int *devchangecount = NULL;

     // 记录细化点数的变量，位于 host 端。
     int changeCount;

     // 记录细化点数的变量，位于 device 端。并为其申请空间。
     cudaerrcode = cudaMalloc((void **)&devchangecount, sizeof (int));
     if (cudaerrcode != cudaSuccess) {
         return CUDA_ERROR;
     }

     // 生成暂存图像。
     errcode = ImageBasicOp::newImage(&tempimg);
     if (errcode != NO_ERROR)
         return errcode;
     errcode = ImageBasicOp::makeAtCurrentDevice(tempimg, inimg->width, 
                                                 inimg->height);
     if (errcode != NO_ERROR) {
         return errcode;
     }

     // 将输入图像 inimg 完全拷贝到输出图像 outimg ，并将 outimg 拷贝到 
     // device 端。
     errcode = ImageBasicOp::copyToCurrentDevice(inimg, outimg);
     if (errcode != NO_ERROR) {
         // FAIL_THIN_IMAGE_FREE;
         return errcode;
     }
     
     // 提取输出图像
     ImageCuda outsubimgCud;
     errcode = ImageBasicOp::roiSubImage(outimg, &outsubimgCud);
     if (errcode != NO_ERROR) {
         // FAIL_THIN_IMAGE_FREE;
         return errcode;
     }

     // 提取暂存图像
     ImageCuda tempsubimgCud;
     errcode = ImageBasicOp::roiSubImage(tempimg, &tempsubimgCud);
     if (errcode != NO_ERROR) {
         // FAIL_THIN_IMAGE_FREE;
         return errcode;
     }

     // 计算调用 Kernel 函数的线程块的尺寸和线程块的数量。
     dim3 gridsize, blocksize;
     blocksize.x = DEF_BLOCK_X;
     blocksize.y = DEF_BLOCK_Y;
     gridsize.x = (outsubimgCud.imgMeta.width + blocksize.x - 1) / blocksize.x;
     gridsize.y = (outsubimgCud.imgMeta.height + blocksize.y - 1) / blocksize.y;

     // 赋值为 1，以便开始第一次迭代。
     changeCount = 1;
     // int iter_num = 0;
     // 开始迭代，当不可再被细化，即记录细化点数的变量 changeCount 的值为 0 时，
     // 停止迭代。 
     while (changeCount > 0) {
        // iter_num ++;
         // 将 host 端的变量赋值为 0 ，并将值拷贝到 device 端的 devchangecount。
         changeCount = 0;
         cudaerrcode = cudaMemcpy(devchangecount, &changeCount, sizeof (int),
                                  cudaMemcpyHostToDevice);
         if (cudaerrcode != cudaSuccess) {
             return CUDA_ERROR;
         }

         // copy ouimg to tempimg 
         cudaerrcode = cudaMemcpyPeer(tempimg->imgData, tempsubimgCud.deviceId, 
                                      outimg->imgData, outsubimgCud.deviceId, 
                                      outsubimgCud.pitchBytes * outimg->height);
        
         if (cudaerrcode != cudaSuccess) {
             return CUDA_ERROR;
         }

         // 调用核函数，开始第一步细化操作。
         _thinZS1Ker<<<gridsize, blocksize>>>(tempsubimgCud, outsubimgCud, devchangecount);
         if (cudaGetLastError() != cudaSuccess) {
             // 核函数出错，结束迭代函数，释放申请的变量空间。
             // FAIL_THIN_IMAGE_FREE;
             return CUDA_ERROR;
         }

         // copy ouimg to tempimg 
         cudaerrcode = cudaMemcpyPeer(tempimg->imgData, tempsubimgCud.deviceId, 
                                      outimg->imgData, outsubimgCud.deviceId, 
                                      outsubimgCud.pitchBytes * outimg->height);
        
         if (cudaerrcode != cudaSuccess) {
             return CUDA_ERROR;
         }

         // 调用核函数，开始第二步细化操作。
         _thinZS2Ker<<<gridsize, blocksize>>>(tempsubimgCud, outsubimgCud, devchangecount);
         if (cudaGetLastError() != cudaSuccess) {
             // 核函数出错，结束迭代函数，释放申请的变量空间 。
             // FAIL_THIN_IMAGE_FREE;
             return CUDA_ERROR;
         }     
        
         // 将位于 device 端的 devchangecount 拷贝到 host 端上的 changeCount 
         // 变量，进行迭代判断。
         cudaerrcode = cudaMemcpy(&changeCount, devchangecount, sizeof (int),
                                  cudaMemcpyDeviceToHost);
         if (cudaerrcode != cudaSuccess) {
             // FAIL_THIN_IMAGE_FREE;
             return CUDA_ERROR;
         }

    }
    // cout << "thinZS iter_num = " << iter_num << endl;
    // 细化结束后释放申请的变量空间。
    cudaFree(devchangecount);
    ImageBasicOp::deleteImage(tempimg);
    return NO_ERROR;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

thinZS函数是一个并行化函数，通过cudaMalloc申请图像空间，然后使用makeAtCurrentDevice将输入图像存入中间变量tempimg，并通过roiSubImage转化为ImageCuda 类型tempsubimgCud。然后进入while循环开始迭代，当不可再被细化，即记录细化点数的变量 changeCount 的值为 0 时，停止迭代。 并在最后细化结束后释放申请的变量空间。
主程序如下：

…详情请参照古月居