CUDA实现图像二值化

发表于 | 分类于 | 阅读次数

着重介绍了阈值二值化及ostu二值化的实现，算法本身很简单，主要是将算法并行化，有点费事，同时也注意到了对于时序性较强的算法，cuda移植后性能并不会得到大幅度提升，提升往往在10倍左右，值得深思。

算法实现：DeltaCV 使用例子：Samples

阈值二值化

阈值二值化，这里是指阈值是人为设定的，简单的二值化，鲁棒性低，但计算量小。

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
/*
 * Compare 'threshold()' funciton in OpenCV
 * When:
 *      thresholdMin = thresholdMax and valMin = 0  ==> THRESH_BINARY
 *      thresholdMin = thresholdMax and valMax = 0  ==> THRESH_BINARY_INV
 *      thresholdMax = valMax and thresholdMin = 0  ==> THRESH_TRUNC
 *      thresholdMax = 255 and valMin = 0  ==> THRESH_TOZERO
 *      thresholdMin = 0 and valMax = 0  ==> THRESH_TOZERO_INV
 */

__global__ void thresholdBinarization(unsigned char* dataIn,
                                      unsigned char* dataOut,
                                      short int imgRows,
                                      short int imgCols,
                                      unsigned char thresholdMin,
                                      unsigned char thresholdMax,
                                      unsigned char valMin,
                                      unsigned char valMax)
{
    int xIndex = threadIdx.x + __umul24(blockIdx.x, blockDim.x);
    int yIndex = threadIdx.y + __umul24(blockIdx.y, blockDim.y);

    int tid = __umul24(yIndex,imgCols)+xIndex;

    unsigned char val=dataIn[tid];
    unsigned char res = val;

    if(xIndex < imgCols && yIndex < imgRows)
    {

        if(val>thresholdMax)
        {
            res = valMax;
        }

        if(val<=thresholdMin)
        {
            res = valMin;
        }

        dataOut[tid] = res;
    }
}

这个算法非常简单，设置与图像像素数一样的thread，thread读取每个像素值，直接进行判断，这里使用分支语句的影响很小，貌似可以被编译器优化，不必太担心效率。 而且通过适当的参数设置（见代码中的注释部分），可以得到与opencv中threshold函数一样的效果。

OSTU二值化

OSTU又称为自适应阈值二值化，它的主要思想是迭代，从灰度级0到255，每次按照当前的灰度值将图像分为出背景和前景，如果当前两个部分（背景和前景）的类间方差最大，则说明此时的阈值是最优的。

求直方图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
__global__ void getHist(unsigned char* dataIn, unsigned int* hist)
{
    int xdx = threadIdx.x + __umul24(blockIdx.x, blockDim.x);
    int ydx = threadIdx.y + __umul24(blockIdx.y, blockDim.y);

    int tid = xdx + ydx*gridDim.x*blockDim.x;

    if(tid < 256)
    {
        hist[tid]=0;
    }
    __syncthreads();
    atomicAdd(&hist[dataIn[tid]],1);
}

这里需要注意的是，由于我们面向视频处理，像atomicAdd这类的操作，之前需要将整个数组清空，因为上一帧的处理结果还在内存中，并没有删除。并且还需要同步一下线程，确保整个hist都已被清空。

计算类间方差

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
/*
 * w_A:前景像素级占的比例，Pi的和
 * u_A:前景像素级的均值，i×Pi的和除以w_A
 */
__global__ void ostu(unsigned int* hist,
                     float* sum_Pi,
                     float* sum_i_Pi,
                     float* u_0,
                     float* varance,
                     short int imgRows,
                     short int imgCols)
{

    //清空相关数组，清空上一次的计算结果
    if(blockIdx.x == 0)
    {
        sum_Pi[threadIdx.x] = 0.0f;
        sum_i_Pi[threadIdx.x] = 0.0f;
        if(threadIdx.x==0)
        {
            u_0[0] = 0.0f;
        }
    }
    __syncthreads();

    //计算整幅图的平均灰度、前景的概率、前景的平均灰度值
    unsigned int current_val = hist[threadIdx.x];
    if(blockIdx.x==0)
    {
        atomicAdd(&u_0[0],current_val*threadIdx.x);//sum(i*Pi)+sum(j*Pj)
    }
    else
    {
        if(threadIdx.x < blockIdx.x)
        {
            atomicAdd(&sum_Pi[blockIdx.x-1],current_val);//sum()
            atomicAdd(&sum_i_Pi[blockIdx.x-1],current_val*threadIdx.x);//sum(i*Pi)
        }
    }
    __syncthreads();
    //now we get sum_Pi[256] and sum_i_Pi[256] and w_0
    //下面开始计算类间方差
    int imgSize = imgRows*imgCols;
    if(blockIdx.x>0)
    {
        float f_sum_pi = sum_Pi[blockIdx.x-1]/imgSize;
        float f_sum_pj = 1-f_sum_pi;
        if(f_sum_pj==0)
        {
            varance[blockIdx.x-1]=0;
        }else
        {
            float temp = (u_0[0]/imgSize-sum_i_Pi[blockIdx.x-1]/(f_sum_pi*imgSize));
            varance[blockIdx.x-1] = temp*temp*f_sum_pi/f_sum_pj;
        }
    }
}

最后面计算类间方差的公式是化简后得到公式，参考。上述程序会得到一个大小为256的类间方差数组，里面的值对应每个灰度级作为二值化阈值时对应的类间方差。

寻找最优阈值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
__global__ void tree_max(float* varance,int* thres)
{
        __shared__ float var[256];
        __shared__ int varId[256]; // 512*4/1024 = 2KB
        varId[threadIdx.x] = threadIdx.x;
        var[threadIdx.x] = varance[threadIdx.x];
        for (int i = 1; i < 256; i*=2)
        {
            if(threadIdx.x%(2*i)==0)
            {
                if(var[threadIdx.x]<var[threadIdx.x+i])
                {
                    var[threadIdx.x] = var[threadIdx.x+i];
                    varId[threadIdx.x] = varId[threadIdx.x+i];
                }
            }
            __syncthreads();
        }
        if(threadIdx.x == 0)
        {
            thres[0] = varId[0];
        }
}

这一步是借用并行求和的思想，来并行求最大值，由于我们最终是要得到类间方差最大时的阈值，也就是类间方差最大那个值的索引，所以我们又定义了一个varId变量来保存每一次树型对比时获得的最大类间方差对应的索引，则对比完后，varId[0]保存的就是最大类间方对应的索引。最后我们再使用之前的阈值二值化进行二值化即可。

整体流程及thread、block设置

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
void ostu_gpu(unsigned char* dataIn,
              unsigned char* dataOut,
              unsigned int* hist,
              float* sum_Pi,
              float* sum_i_Pi,
              float* u_0,
              float* varance,
              int* thres,
              short int imgRows,
              short int imgCols)
{

    dim3 tPerBlock_hist(32,32);
    dim3 bPerGrid_hist((imgCols+32-1)/32,(imgRows+32-1)/32);
    //求直方图
    getHist_gpu(dataIn,hist,tPerBlock_hist,bPerGrid_hist);

    dim3 tPerBlock_ostu(256,1);
    dim3 bPerGrid_ostu(257,1);
    //计算类间方差
    ostu<<<bPerGrid_ostu,tPerBlock_ostu>>>(hist,sum_Pi,sum_i_Pi,u_0,varance,imgRows,imgCols);
    //寻找最优阈值
    tree_max<<<1,256>>>(varance,thres);
    //阈值值化
    thresholdBinarization_inner<<<bPerGrid_hist,tPerBlock_hist>>>(dataIn,dataOut,imgRows,imgCols,thres);
}