前言
cuda对图像处理过程进行加速是很常见的操作，而且图像处理算法中sobel算子又是最为简单的算法，本篇博客就是使用cuda opencv c++实现对图像进行sobel边缘提取加速操作，cpu和gpu处理速度进行一个对比试验
具体流程
首先贴上代码
#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;
//Sobel算子边缘检测核函数
__global__ void sobelInCuda(unsigned char* dataIn, unsigned char* dataOut, int imgHeight, int imgWidth)
{
int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
int index = yIndex * imgWidth + xIndex;
//printf("xIndex : %d, yIndex: %d\n", xIndex, yIndex);
//printf("index : %d\n", index);
int Gx = 0;
int Gy = 0;
if (xIndex > 0 && xIndex < imgWidth - 1 && yIndex > 0 && yIndex < imgHeight - 1)
{
Gx = dataIn[(yIndex - 1) * imgWidth + xIndex + 1] + 2 * dataIn[yIndex * imgWidth + xIndex + 1] + dataIn[(yIndex + 1) * imgWidth + xIndex + 1]
- (dataIn[(yIndex - 1) * imgWidth + xIndex - 1] + 2 * dataIn[yIndex * imgWidth + xIndex - 1] + dataIn[(yIndex + 1) * imgWidth + xIndex - 1]);
Gy = dataIn[(yIndex - 1) * imgWidth + xIndex - 1] + 2 * dataIn[(yIndex - 1) * imgWidth + xIndex] + dataIn[(yIndex - 1) * imgWidth + xIndex + 1]
- (dataIn[(yIndex + 1) * imgWidth + xIndex - 1] + 2 * dataIn[(yIndex + 1) * imgWidth + xIndex] + dataIn[(yIndex + 1) * imgWidth + xIndex + 1]);
dataOut[index] = (abs(Gx) + abs(Gy)) / 2;
}
}
//Sobel算子边缘检测CPU函数
void sobel(Mat srcImg, Mat dstImg, int imgHeight, int imgWidth)
{
int Gx = 0;
int Gy = 0;
for (int i = 1; i < imgHeight - 1; i++)
{
uchar* dataUp = srcImg.ptr<uchar>(i - 1);
uchar* data = srcImg.ptr<uchar>(i);
uchar* dataDown = srcImg.ptr<uchar>(i + 1);
uchar* out = dstImg.ptr<uchar>(i);
for (int j = 1; j < imgWidth - 1; j++)
{
Gx = (dataUp[j + 1] + 2 * data[j + 1] + dataDown[j + 1]) - (dataUp[j - 1] + 2 * data[j - 1] + dataDown[j - 1]);
Gy = (dataUp[j - 1] + 2 * dataUp[j] + dataUp[j + 1]) - (dataDown[j - 1] + 2 * dataDown[j] + dataDown[j + 1]);
out[j] = (abs(Gx) + abs(Gy)) / 2;
}
}
}
int main()
{
Mat grayImg = imread("lena.jpg", 0);
int imgHeight = grayImg.rows;
int imgWidth = grayImg.cols;
Mat gaussImg;
//高斯滤波
GaussianBlur(grayImg, gaussImg, Size(3, 3), 0, 0, BORDER_DEFAULT);
//Sobel算子CPU实现
Mat dst(imgHeight, imgWidth, CV_8UC1, Scalar(0));
double time1 = static_cast<double>(cv::getTickCount());
sobel(gaussImg, dst, imgHeight, imgWidth);
double time2 = static_cast<double>(cv::getTickCount());
std::cout << "cpu Time use: " << 1000 * (time2 - time1) / cv::getTickFrequency() << "ms" << std::endl;//输出运行时间
cv::imwrite("dst_cpu.jpg", dst);
//CUDA实现后的传回的图像
Mat dstImg(imgHeight, imgWidth, CV_8UC1, Scalar(0));
//创建GPU内存
unsigned char* d_in;
unsigned char* d_out;
cudaMalloc((void**)&d_in, imgHeight * imgWidth * sizeof(unsigned char));
cudaMalloc((void**)&d_out, imgHeight * imgWidth * sizeof(unsigned char));
//将高斯滤波后的图像从CPU传入GPU
cudaMemcpy(d_in, gaussImg.data, imgHeight * imgWidth * sizeof(unsigned char), cudaMemcpyHostToDevice);
dim3 threadsPerBlock(32, 32);
dim3 blocksPerGrid((imgWidth + threadsPerBlock.x - 1) / threadsPerBlock.x, (imgHeight + threadsPerBlock.y - 1) / threadsPerBlock.y);
//记录起始时间
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
//调用核函数
sobelInCuda << <blocksPerGrid, threadsPerBlock >> > (d_in, d_out, imgHeight, imgWidth);
//将图像传回GPU
cudaMemcpy(dstImg.data, d_out, imgHeight * imgWidth * sizeof(unsigned char), cudaMemcpyDeviceToHost);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
float elapsedTime;
cudaEventElapsedTime(&elapsedTime, start, stop);
printf("gpu time cost: %3.5f ms", elapsedTime);
cv::imwrite("dst_gpu.jpg", dstImg);
//释放GPU内存
cudaFree(d_in);
cudaFree(d_out);
cudaEventDestroy(start);
cudaEventDestroy(stop);
return 0;
}
当然，这里的代码是参考别人。下面着重讲讲对这部分代码的理解吧
首先，Sobel算子：是离散微分算子（discrete differentiation operator），用来计算图像灰度的近似梯度，梯度越大越有可能是边缘。sobel边缘提取就是使用sobel算子对图像进行卷积操作，这里sobel算子选用的是一个3*3的矩阵
sobel算子cpu的实现，这样没有太多可说的，基本上就是对一个二维矩阵进行矩阵乘积操作，然后再将矩阵乘积结果回填到目标像素位置
而sobel算子的gpu实现，则相对来说比较复杂一些。首先创建了block和grid的尺寸，即代码中的threadsPerBlock和blocksPerGrid，使用cudaMalloc分配图像尺寸大小的内存空间，用来作为核函数的输入变量和输出变量，最后调用核函数完成并行操作，即sobel边缘提取
再看看核函数内部的实现
int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
这两行代码主要是获取线程块内的坐标位置，这个位置和图像矩阵位置是一一对应的
int index = yIndex * imgWidth + xIndex;
而这行代码，则是实现图像矩阵坐标和data_in， data_out这两个指针里面存储的一维数组的映射关系，(xIndex, yIndex)该坐标像素值即为每次线程运行时，做sobel运算时的中心坐标，而实际求Gx和Gy时，还需要获取该中心坐标周围的8个坐标值来进行运算，所以是这样进行书写的代码
Gx = dataIn[(yIndex - 1) * imgWidth + xIndex + 1] + 2 * dataIn[yIndex * imgWidth + xIndex + 1] + dataIn[(yIndex + 1) * imgWidth + xIndex + 1]
- (dataIn[(yIndex - 1) * imgWidth + xIndex - 1] + 2 * dataIn[yIndex * imgWidth + xIndex - 1] + dataIn[(yIndex + 1) * imgWidth + xIndex - 1]);
Gy = dataIn[(yIndex - 1) * imgWidth + xIndex - 1] + 2 * dataIn[(yIndex - 1) * imgWidth + xIndex] + dataIn[(yIndex - 1) * imgWidth + xIndex + 1]
- (dataIn[(yIndex + 1) * imgWidth + xIndex - 1] + 2 * dataIn[(yIndex + 1) * imgWidth + xIndex] + dataIn[(yIndex + 1) * imgWidth + xIndex + 1]);
运算的结果存回到dataOut
整个测试中输入图片是女神lena
cpu和gpu输出的结果为
可以看出，gpu和cpu的结果图是一模一样的
cpu和gpu的时间测试
上一部分叙述了cpu和gpu版本的sobel算子边缘提取，现在来测试一下两个版本算法的运行时间
cpu的版本使用的时间测试工具为opencv中自带的计时工具，我一直在使用，比较靠谱。测试方式也比较简单，就是在运行sobel函数前后加上 static_cast<double>(cv::getTickCount())函数，然后前后两次count数相减再除于count的频次就可以了，这里我将count乘以1000，这样测试得到的时间基本单位就是ms
gpu版本使用的时间测试工具是cuda自带的工具cudaEventElapsedTime，在计算时延之前还需要进行cudaEventRecord操作，避免没有gpu线程没有完全执行结束
最终测试的结果如下：
可以看出，对于同一个输入图片，gpu的运行时间是cpu大概1/5的样子
这里我对gpu进行计时的时候将图像从gpu传到cpu的操作也进行了计时，如果将这部分除去，gpu的时间会更短
可以看出，这次sobel算子边缘提取，真正在gpu上进行运算的时间仅为0.02ms，绝大部分时间是将数据从gpu拷贝到cpu的时间，这个时间会占据整个流程的85%左右
总结
整个探索过程还是比较有意义，既熟悉了sobel算子的原理，又熟悉了gpu并行运算的整体流程，还做cpu与gpu运算版本时间对比。
