引言
在我的第一篇文章中我简单介绍了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:8000
for j=1:1000
img(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*20
dim3 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倍~