我的CUDA學習之旅1——大圖像分塊處理程序

01-26

前言

在我的第一篇文章中本人簡單介紹了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)); endendimwrite(img,test.jpg)

只是一個隨機產生0-255值的程序，可產生一張隨機的灰度圖：

（2）CPU代碼實現

我先使用Opencv的方式在CPU端實現了任務，得到了計算結果（方便與GPU實現相比較，同時也可驗證CUDA代碼的準確性），具體的代碼如下：

#include <iostream>#include <opencv2opencv.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代碼

#include "cuda_runtime.h"#include "device_launch_parameters.h"#include <cuda.h>#include <device_functions.h>#include <opencv2opencv.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倍~