积分图的一种CUDA并行运算
在前面的文章中，我们有讲积分图的基本原理、算法层面优化加速，以及SSE指令优化加速：
https://blog.csdn.net/shandianfengfan/article/details/109571376
https://blog.csdn.net/shandianfengfan/article/details/109571418
最近在项目中需要使用CUDA来优化非局部均值滤波算法，涉及到积分图的并行计算。因此本文我们来讲一下积分图的并行计算思路，并使用C++和CUDA来实现。
1. 前缀和的概念
我们举个简单的例子来说明前缀和，假设有一组数据：
其前缀和的计算过程如下：
以上计算过程概括成下面的公式：
所以，计算数组中每一个元素的前缀和，就是计算该元素与其前面的所有元素之和。
2. 分行、列并行计算
假设有一张m行n列的图像，我们知道，先分别计算其所有行的前缀和，得到行前缀和图像，再对行前缀和图像分别计算所有列的前缀和，即得到原图像的积分图。如下图所示：
以上计算过程为两个步骤：1. 计算行前缀和；2. 计算列前缀和。所以并行计算思路为：
1. 开启m个线程，每个线程分别对应原图像的一行，负责计算该行的前缀和，所有线程并行计算。该步骤计算结束之后，得到了m*n的行前缀和图像。
2. 开启n个线程，每个线程分别对应行前缀和图像的一列，负责计算该列的前缀和，所有线程并行计算。该步骤计算结束之后，则得到m*n的积分图像。
3. 分块并行计算
这是本文要实现的计算方法，确切地说，该方法也属于分行、列并行计算的方法，只不过它在分行列计算地基础上，进一步分块计算，看起来步骤变多了，实际上计算速度加快不少。
同样假设有一张m行n列的图像，现在要使用分块的思路并行计算其积分图，主要有6个步骤：
步骤一：将每一行数据分成每4个数据的数据块，开启m*(n/4)个线程，每个线程计算一个小块中4个数据的前缀和，所有线程并行计算。
步骤二：开启m个线程，每个线程对应一行数据，负责计算该行中所有小块中最后一个数据的前缀和。
步骤三：开启m*(n/4)个线程，每个线程对应一个小块。对于每一行，从第2个小块开始，将每小块中的数（除最后一个数外）加上它的前一小块的最后一个数，即可得到所有行的前缀和。
步骤四：将每一列数据分成每4个数据的数据块，开启(m/4)*n个线程，每个线程计算一个小块中4个数据的前缀和，所有线程并行计算。
步骤五：开启n个线程，每个线程对应一列数据，负责计算该列中所有小块中最后一个数据的前缀和。
步骤六：开启(m/4)*n个线程，每个线程对应一个小块。对于每一列，从第2个小块开始，将每小块中的数（除最后一个数外）加上它的前一小块的最后一个数，即可得到所有列的前缀和，也即原图的积分图像。
注意，这种分块的方法要求图像的行与列都是4的整数倍，如果不是则需要边缘填充使其为4的整数倍。
下面分步骤上代码。
步骤一核函数代码：
__global__ void cuda_integal0_rc(float *src, float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //col/4
int idy = threadIdx.y + blockDim.y * blockIdx.y; //row
if (idx < col / 4 && idy < row)
{
int offset = idy*col + idx * 4; //每个小块的首地址
integral[offset] = src[offset];
integral[offset + 1] = src[offset + 1] + integral[offset];
integral[offset + 2] = src[offset + 2] + integral[offset + 1];
integral[offset + 3] = src[offset + 3] + integral[offset + 2];
}
}
步骤二核函数代码：
__global__ void cuda_integal1_rc(float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //row
if (idx < row)
{
int offset = idx*col;
for (int i = 7; i < col; i += 4)
{
integral[offset + i] += integral[offset + i - 4];
}
}
}
步骤三核函数代码：
__global__ void cuda_integal2_rc(float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //col/4
int idy = threadIdx.y + blockDim.y * blockIdx.y; //row
if (idx > 0 && idx < col / 4 && idy < row)
{
int offset = idy*col + idx * 4; //每个小块的首地址
float pre_sum = integral[idy*col + (idx*4 - 1)]; //前一个小块的最后一个数
integral[offset] += pre_sum;
integral[offset + 1] += pre_sum;
integral[offset + 2] += pre_sum;
}
}
步骤四核函数代码：
__global__ void cuda_integal3_rc(float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //col
int idy = threadIdx.y + blockDim.y * blockIdx.y; //row/4
if (idx < col && idy < row / 4)
{
int offset = idy*4*col + idx; //每个小块的首地址
int offset1 = offset + col;
int offset2 = offset1 + col;
int offset3 = offset2 + col;
integral[offset1] += integral[offset];
integral[offset2] += integral[offset1];
integral[offset3] += integral[offset2];
}
}
步骤五核函数代码：
__global__ void cuda_integal4_rc(float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //col
if (idx < col)
{
int col_4 = 4 * col;
for (int i = 7; i < row; i += 4)
{
int offset = i*col + idx;
integral[offset] += integral[offset - col_4];
}
}
}
步骤六核函数代码：
__global__ void cuda_integal5_rc(float *integral, int row, int col)
{
int idx = threadIdx.x + blockDim.x * blockIdx.x; //col
int idy = threadIdx.y + blockDim.y * blockIdx.y; //row/4
if (idx < col && idy > 0 && idy < row / 4)
{
int offset = idy*4*col + idx; //每个小块的首地址
float pre_sum = integral[(idy*4 - 1)*col + idx]; //前一个小块的最后一个数
integral[offset] += pre_sum;
integral[offset + col] += pre_sum;
integral[offset + col + col] += pre_sum;
}
}
主体函数代码：
void integal_cuda_rc(Mat src, Mat &dst)
{
Mat tmp;
src.convertTo(tmp, CV_32F);
int row_r = tmp.rows % 4;
int col_r = tmp.cols % 4;
if (row_r || col_r) //如果行与列不是4的整数倍，则边缘填充
{
copyMakeBorder(tmp, tmp, 0, row_r, 0, col_r, BORDER_REFLECT); //扩充边缘
}
const int row = tmp.rows;
const int col = tmp.cols;
float *integral;
cudaMalloc((void**)&integral, row*col*sizeof(float));
dim3 k0block(16, 16); //row*(col/4)
dim3 k0grid((col/4 + k0block.x - 1) / k0block.x, (row + k0block.y - 1) / k0block.y);
dim3 k1block(16, 1); //row
dim3 k1grid((row + k1block.x - 1) / k1block.x, 1);
dim3 k2block = k0block; //row*(col/4)
dim3 k2grid = k0grid;
dim3 k3block(16, 16); //(row/4)*col
dim3 k3grid((col + k3block.x - 1) / k3block.x, (row / 4 + k3block.y - 1) / k3block.y);
dim3 k4block(16, 1); //col
dim3 k4grid((col + k4block.x - 1) / k4block.x, 1);
dim3 k5block = k3block;
dim3 k5grid = k3grid;
float *src_cuda;
cudaMalloc((void**)&src_cuda, row*col * sizeof(float));
cudaMemcpy(src_cuda, (float *)tmp.data, row*col * sizeof(float), cudaMemcpyHostToDevice);
cuda_integal0_rc << <k0grid, k0block >> >(src_cuda, integral, row, col); //步骤一
cuda_integal1_rc << <k1grid, k1block >> >(integral, row, col); //步骤二
cuda_integal2_rc << <k2grid, k2block >> >(integral, row, col); //步骤三
cuda_integal3_rc << <k3grid, k3block >> >(integral, row, col); //步骤四
cuda_integal4_rc << <k4grid, k4block >> >(integral, row, col); //步骤五
cuda_integal5_rc << <k5grid, k5block >> >(integral, row, col); //步骤六
cudaMemcpy((float *)tmp.data, integral, row*col * sizeof(float), cudaMemcpyDeviceToHost);
cudaFree(integral);
tmp(Rect(0, 0, src.cols, src.rows)).copyTo(dst);
}
好啦，本文就写到这里，下一篇文章更精彩哦，敬请期待！
微信公众号如下，欢迎扫码关注，欢迎私信技术交流：
