CUDA精进之路（一）：图像处理——大图像分块处理（包括求均值、最大值）

引言

在我的第一篇文章中我简单介绍了CUDA以及我的一些个人学习见解，在本文中我将开始正式开始CUDA实践之旅，众做周知CUDA目前应用的领域十分广泛，它能把一些普通的CPU代码提速几十倍甚至几百倍。在本人所从事的图像处理领域，在一些大图像的处理上（4K以上图像），仅仅依靠CPU进行计算已经完全无法满足工程项目所要求的运行时间，这时候我们就需要利用CUDA对代码进行加速。本文以一个8000*1000图像为例，进行代码实践。

任务要求

将一个8000*1000的单通道图分割为40x40 (可调整)的块, 计算每个块内的像素值的和、最大、最小、均值，将计算结果保存在CPU端。

实现思路

由于我的本本的显卡最大支持的thread数为1024，而任务中40*40的数量显然已经超过了1024，故在计算任务的分配上，我综合考虑后选择了通过两次运算（每次计算数为40*20）来完成，运用共享内存保存每个块中传入图像的像素数据，最后利用归约算法对加法进行优化。

实现环境

VS2013 + CUDA7.5 + Opencv2.4.13

实现代码

1）8000*1000实验图像获取

我用了一段简单的Matlab代码来获取实验图像，代码如下：

img = uint8(zeros(8000,1000));
for i=1:8000for j=1:1000img(i,j) = uint8(255*rand(1));end
end
imwrite(img,'test.jpg')

该程序的目的在于产生一个随机产生0-255值的灰度图像，大小为8000*1000，作为我们的实验图像。

（2）CPU代码实现

我先使用Opencv的方式在CPU端实现了任务，得到了计算结果（方便与GPU实现相比较，同时也可验证CUDA代码的准确性），具体的代码如下：

#include <iostream>
#include <opencv2\opencv.hpp>
#include <stdlib.h>
using namespace std;
using namespace cv;int operateMat(Mat img, int startx, int starty, int size, int thresh);int main()
{Mat img = imread("test.jpg",0);int size = 40;int thresh = 0;cout << "Plaese Input Size:";cin >> size;   //分块区域大小cout << "Plaese Input Thresh value:";cin >> thresh; //阈值double time0 = static_cast<double>(getTickCount()); //计时器开始int row = img.rows;int col = img.cols;int length = row / size;int width = col / size;int sumResult[10000];     //区域内求和结果int averageResult[10000]; //区域内均值结果int maxResult[10000]; //区域内最大值结果int minResult[10000]; //区域内最小值结果int threshNumber[10000];  //区域内阈值结果int count = 0;int x = 0;for (int i = 0; i < length; i++){for (int j = 0; j < width; j++){sumResult[count], maxResult[count], minResult[count], threshNumber[count]  = operateMat(img, i*size, j*size, size, 35);averageResult[count] = sumResult[count] / (size * size);count += 1;}}time0 = ((double)getTickCount() - time0) / getTickFrequency(); //计时器结束cout << time0 << "s" << endl; //输出运行时间system("Pause");return 0;
}int operateMat(Mat img, int startx, int starty, int size, int thresh)
{int sum = 0;int max = 0;int min = 255;int average = 0;int threshCount = 0;for (int i = startx; i < startx + size; i++){uchar *data = img.ptr<uchar>(i);for (int j = starty; j < starty + size; j++){sum += data[j];if (max < data[j]){max = data[j];}if (min > data[j]){min = data[j];}if (data[j] > thresh){threshCount += 1;}}}average = sum / (size*size);return sum, max, min, threshCount;
}

（3）CUDA代码
在CUDA实现上，本程序利用了每个块独有的共享内存，将需要计算的40*40的小块中所有的像素数据通过两次导入的方式（每次导入40*20）全部赋值给大小为1600的共享内存，接着通过规约算法优化运算即可得到所需结果。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cuda.h>
#include <device_functions.h>
#include <opencv2\opencv.hpp>
#include <iostream>
using namespace std;
using namespace cv;__global__ void matSum(uchar *dataIn, int *dataOutSum, int *dataOutMax, int *dataOutMin, int imgHeight, int imgWidth)
{//__shared__ int _data[1600];const int number = 2048;extern __shared__ int _sum[];  //小图像块中求和共享数组__shared__ int _max[number];  //小图像块中求最大值共享数组__shared__ int _min[number];  //小图像块中求最小值共享数组int thread = threadIdx.x + threadIdx.y * blockDim.x; //一个block中所有thread的索引值int threadIndex = threadIdx.x + threadIdx.y * imgWidth; //每个小块中存放数据的thread索引值//每个小块中存放数据的block索引值int blockIndex1 = blockIdx.x * blockDim.x + 2 * blockIdx.y * blockDim.y * imgWidth; //40*20的上半block索引值int blockIndex2 = blockIdx.x * blockDim.x + (2 * blockIdx.y + 1) * blockDim.y * imgWidth; //40*20的下半block索引值int index1 = threadIndex + blockIndex1; //每个block中上半部分索引值int index2 = threadIndex + blockIndex2; //每个block中下半部分索引值//将待计算的40*40小图像块中的所有像素值分两次传送到共享数组中_sum[thread] = dataIn[index1]; //将上半部分的40*20中所有数据赋值到共享数组中_sum[thread + blockDim.x * blockDim.y] = dataIn[index2]; //将下半部分的40*20中所有数据赋值到共享数组中_max[thread] = dataIn[index1];_max[thread + blockDim.x * blockDim.y] = dataIn[index2];_min[thread] = dataIn[index1];_min[thread + blockDim.x * blockDim.y] = dataIn[index2];//memcpy(_sum, _data, 1600 * sizeof(int));//memcpy(_max, _data, 1600 * sizeof(int));//memcpy(_min, _data, 1600 * sizeof(int));  在GPU（Device）中用memcpy函数进行拷贝会导致显卡混乱，故不选择此种方式//利用归约算法求出40*40小图像块中1600个像素值中的和、最大值以及最小值for (unsigned int s = number / 2; s > 0; s >>= 1){if (thread < s){ _sum[thread] += _sum[thread + s]; if (_max[thread] < _max[thread + s]) { _max[thread] = _max[thread + s]; }if (_min[thread] > _min[thread + s]) { _min[thread] = _min[thread + s]; }}__syncthreads(); //所有线程同步}if (threadIndex == 0) { //将每个小块中的结果储存到输出中dataOutSum[blockIdx.x + blockIdx.y * gridDim.x] = _sum[0]; dataOutMax[blockIdx.x + blockIdx.y * gridDim.x] = _max[0];dataOutMin[blockIdx.x + blockIdx.y * gridDim.x] = _min[0];}}int main()
{Mat image = imread("test.jpg", 0); //读取待检测图片int sum[5000]; //求和结果数组int max[5000]; //最大值结果数组int min[5000]; //最小值结果数组imshow("src", image);size_t memSize = image.cols*image.rows*sizeof(uchar);int size = 5000 * sizeof(int);uchar *d_src = NULL;int *d_sum = NULL;int *d_max = NULL;int *d_min = NULL;cudaMalloc((void**)&d_src, memSize);cudaMalloc((void**)&d_sum, size);cudaMalloc((void**)&d_max, size);cudaMalloc((void**)&d_min, size);cudaMemcpy(d_src, image.data, memSize, cudaMemcpyHostToDevice);int imgWidth = image.cols;int imgHeight = image.rows;dim3 threadsPerBlock(40, 20); //每个block大小为40*20dim3 blockPerGrid(25, 200); //将8000*1000的图片分为25*200个小图像块double time0 = static_cast<double>(getTickCount()); //计时器开始matSum << <blockPerGrid, threadsPerBlock, 4096 * sizeof(int) >> >(d_src, d_sum, d_max, d_min, imgHeight, imgWidth);time0 = ((double)getTickCount() - time0) / getTickFrequency(); //计时器结束cout << "The Run Time is :" << time0 << "s" << endl; //输出运行时间cudaMemcpy(sum, d_sum, size, cudaMemcpyDeviceToHost);cudaMemcpy(max, d_max, size, cudaMemcpyDeviceToHost);cudaMemcpy(min, d_min, size, cudaMemcpyDeviceToHost);cout << "The sum is :" << sum[0] << endl;cout << "The max is :" << max[0] << endl;cout << "The min is :" << min[0] << endl;waitKey(0);cudaFree(d_src);cudaFree(d_sum);cudaFree(d_max);cudaFree(d_min);return 0;
}

（4）运算结果比对

CPU与GPU代码运算结果一致，在性能上CPU端代码耗时约0.1S,GPU端代码为0.2ms，加速500倍~