引言

在我的第一篇文章中我簡單介紹了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倍~