在前面的文章中，我们有讲积分图的基本原理、算法层面优化加速，以及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);
-
}
好啦，本文就写到这里，下一篇文章更精彩哦，敬请期待！
微信公众号如下，欢迎扫码关注，欢迎私信技术交流：
转载：https://blog.csdn.net/shandianfengfan/article/details/113449531