目標

在使用OpenCV進行影片輸入和相似性測量教程中，我已經介紹了用於檢查兩張影像相似性的PSNR和SSIM方法。正如你所看到的，執行過程相當耗時，尤其是在SSIM的情況下。然而，如果OpenCV針對CPU的實現效能不能滿足你的需求，並且你的系統恰好擁有NVIDIA CUDA GPU裝置，那麼一切都還沒有失去。你可以嘗試將自己的演算法移植或編寫到顯示卡上。

本教程將詳細介紹如何使用OpenCV的GPU模組進行編碼。作為先決條件，你應該已經知道如何處理core、highgui和imgproc模組。所以，我們的主要目標是：

與CPU相比有什麼不同？
為PSNR和SSIM建立GPU程式碼
最佳化程式碼以獲得最大效能

原始碼

你也可以在OpenCV源庫的samples/cpp/tutorial_code/gpu/gpu-basics-similarity/gpu-basics-similarity目錄中找到原始碼和影片檔案，或者從此處下載。完整的原始碼相當長（由於透過命令列引數控制應用程式和效能測量）。因此，為了避免這些部分過於混亂，這裡只提供了函式本身。

PSNR返回一個浮點數，如果兩個輸入影像相似，該值介於30到50之間（越高越好）。

double getPSNR(const Mat& I1, const Mat& I2)
{
    Mat s1;
    absdiff(I1, I2, s1); // |I1 - I2|
s1.convertTo(s1, CV_32F); // 不能對8位值進行平方運算
s1 = s1.mul(s1); // |I1 - I2|^2
 
    Scalar s = sum(s1); // 每通道求和
 
    double sse = s.val[0] + s.val[1] + s.val[2]; // 各通道求和
 
    if( sse <= 1e-10) // 對於小值返回零
        return 0;
    else
    {
        double mse =sse /(double)(I1.channels() * I1.total());
        double psnr = 10.0*log10((255*255)/mse);
        return psnr;
    }
}

double getPSNR_CUDA(const Mat& I1, const Mat& I2)
{
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
gI1.upload(I1);
gI2.upload(I2);
 
gI1.convertTo(t1, CV_32F);
gI2.convertTo(t2, CV_32F);
 
cuda::absdiff(t1.reshape(1), t2.reshape(1), gs);
cuda::multiply(gs, gs, gs);
 
    Scalar s = cuda::sum(gs);
    double sse = s.val[0] + s.val[1] + s.val[2];
 
    if( sse <= 1e-10) // 對於小值返回零
        return 0;
    else
    {
        double mse =sse /(double)(gI1.channels() * I1.total());
        double psnr = 10.0*log10((255*255)/mse);
        return psnr;
    }
}

struct BufferPSNR // 最佳化的CUDA版本
{ // 在CUDA上資料分配非常昂貴。使用緩衝區解決：分配一次，以後重複使用。
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
    cuda::GpuMat buf;
};

double getPSNR_CUDA_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)
{
b.gI1.upload(I1);
b.gI2.upload(I2);
 
b.gI1.convertTo(b.t1, CV_32F);
b.gI2.convertTo(b.t2, CV_32F);
 
cuda::absdiff(b.t1.reshape(1), b.t2.reshape(1), b.gs);
cuda::multiply(b.gs, b.gs, b.gs);
 
    double sse = cuda::sum(b.gs, b.buf)[0];
 
    if( sse <= 1e-10) // 對於小值返回零
        return 0;
    else
    {
        double mse = sse /(double)(I1.channels() * I1.total());
        double psnr = 10.0*log10((255*255)/mse);
        return psnr;
    }
}

SSIM返回影像的MSSIM。這也是一個介於零和一之間的浮點數（越高越好），但是每個通道有一個值。因此，我們返回一個Scalar OpenCV資料結構

Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
    const double C1 = 6.5025, C2 = 58.5225;
/***************************** 初始化 **********************************/
    int d = CV_32F;
 
    Mat I1, I2;
i1.convertTo(I1, d); // 不能在一個位元組大的值上計算
i2.convertTo(I2, d);
 
    Mat I2_2 = I2.mul(I2); // I2^2
    Mat I1_2 = I1.mul(I1); // I1^2
    Mat I1_I2 = I1.mul(I2); // I1 * I2
 
/*************************** 初始化結束 **********************************/
 
    Mat mu1, mu2; // 初步計算
    GaussianBlur(I1, mu1, Size(11, 11), 1.5);
    GaussianBlur(I2, mu2, Size(11, 11), 1.5);
 
    Mat mu1_2 = mu1.mul(mu1);
    Mat mu2_2 = mu2.mul(mu2);
    Mat mu1_mu2 = mu1.mul(mu2);
 
    Mat sigma1_2, sigma2_2, sigma12;
 
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
sigma1_2 -= mu1_2;
 
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
sigma2_2 -= mu2_2;
 
    GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
sigma12 -= mu1_mu2;
 
    Mat t1, t2, t3;
 
t1 = 2 * mu1_mu2 + C1;
t2 = 2 * sigma12 + C2;
t3 = t1.mul(t2); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
t1 = mu1_2 + mu2_2 + C1;
t2 = sigma1_2 + sigma2_2 + C2;
t1 = t1.mul(t2); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
    Mat ssim_map;
    divide(t3, t1, ssim_map); // ssim_map = t3./t1;
 
    Scalar mssim = mean( ssim_map ); // mssim = ssim圖的平均值
    return mssim;
}

Scalar getMSSIM_CUDA( const Mat& i1, const Mat& i2)
{
    const float C1 = 6.5025f, C2 = 58.5225f;
/***************************** 初始化 **********************************/
    cuda::GpuMat gI1, gI2, gs1, tmp1,tmp2;
 
gI1.upload(i1);
gI2.upload(i2);
 
gI1.convertTo(tmp1, CV_MAKE_TYPE(CV_32F, gI1.channels()));
gI2.convertTo(tmp2, CV_MAKE_TYPE(CV_32F, gI2.channels()));
 
vector<cuda::GpuMat> vI1, vI2;
cuda::split(tmp1, vI1);
cuda::split(tmp2, vI2);
    Scalar mssim;
 
    Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(vI2[0].type(), -1, Size(11, 11), 1.5);
 
    for( int i = 0; i < gI1.channels(); ++i )
    {
        cuda::GpuMat I2_2, I1_2, I1_I2;
 
cuda::multiply(vI2[i], vI2[i], I2_2); // I2^2
cuda::multiply(vI1[i], vI1[i], I1_2); // I1^2
cuda::multiply(vI1[i], vI2[i], I1_I2); // I1 * I2
 
/*************************** 初始化結束 **********************************/
        cuda::GpuMat mu1, mu2; // 初步計算
gauss->apply(vI1[i], mu1);
gauss->apply(vI2[i], mu2);
 
        cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
cuda::multiply(mu1, mu1, mu1_2);
cuda::multiply(mu2, mu2, mu2_2);
cuda::multiply(mu1, mu2, mu1_mu2);
 
        cuda::GpuMat sigma1_2, sigma2_2, sigma12;
 
gauss->apply(I1_2, sigma1_2);
cuda::subtract(sigma1_2, mu1_2, sigma1_2); // sigma1_2 -= mu1_2;
 
gauss->apply(I2_2, sigma2_2);
cuda::subtract(sigma2_2, mu2_2, sigma2_2); // sigma2_2 -= mu2_2;
 
gauss->apply(I1_I2, sigma12);
cuda::subtract(sigma12, mu1_mu2, sigma12); // sigma12 -= mu1_mu2;
 
        cuda::GpuMat t1, t2, t3;
 
mu1_mu2.convertTo(t1, -1, 2, C1); // t1 = 2 * mu1_mu2 + C1;
sigma12.convertTo(t2, -1, 2, C2); // t2 = 2 * sigma12 + C2;
cuda::multiply(t1, t2, t3); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
cuda::addWeighted(mu1_2, 1.0, mu2_2, 1.0, C1, t1); // t1 = mu1_2 + mu2_2 + C1;
cuda::addWeighted(sigma1_2, 1.0, sigma2_2, 1.0, C2, t2); // t2 = sigma1_2 + sigma2_2 + C2;
cuda::multiply(t1, t2, t1); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
        cuda::GpuMat ssim_map;
cuda::divide(t3, t1, ssim_map); // ssim_map = t3./t1;
 
        Scalar s = cuda::sum(ssim_map);
mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols);
 
    }
    return mssim;
}

struct BufferMSSIM // 最佳化的CUDA版本
{ // 在CUDA上資料分配非常昂貴。使用緩衝區解決：分配一次，以後重複使用。
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
    cuda::GpuMat I1_2, I2_2, I1_I2;
vector<cuda::GpuMat> vI1, vI2;
 
    cuda::GpuMat mu1, mu2;
    cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
 
    cuda::GpuMat sigma1_2, sigma2_2, sigma12;
    cuda::GpuMat t3;
 
    cuda::GpuMat ssim_map;
 
    cuda::GpuMat buf;
};

Scalar getMSSIM_CUDA_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b)
{
    const float C1 = 6.5025f, C2 = 58.5225f;
/***************************** 初始化 **********************************/
 
b.gI1.upload(i1);
b.gI2.upload(i2);
 
    cuda::Stream stream;
 
b.gI1.convertTo(b.t1, CV_32F, stream);
b.gI2.convertTo(b.t2, CV_32F, stream);
 
cuda::split(b.t1, b.vI1, stream);
cuda::split(b.t2, b.vI2, stream);
    Scalar mssim;
 
    Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(b.vI1[0].type(), -1, Size(11, 11), 1.5);
 
    for( int i = 0; i < b.gI1.channels(); ++i )
    {
cuda::multiply(b.vI2[i], b.vI2[i], b.I2_2, 1, -1, stream); // I2^2
cuda::multiply(b.vI1[i], b.vI1[i], b.I1_2, 1, -1, stream); // I1^2
cuda::multiply(b.vI1[i], b.vI2[i], b.I1_I2, 1, -1, stream); // I1 * I2
 
gauss->apply(b.vI1[i], b.mu1, stream);
gauss->apply(b.vI2[i], b.mu2, stream);
 
cuda::multiply(b.mu1, b.mu1, b.mu1_2, 1, -1, stream);
cuda::multiply(b.mu2, b.mu2, b.mu2_2, 1, -1, stream);
cuda::multiply(b.mu1, b.mu2, b.mu1_mu2, 1, -1, stream);
 
gauss->apply(b.I1_2, b.sigma1_2, stream);
cuda::subtract(b.sigma1_2, b.mu1_2, b.sigma1_2, cuda::GpuMat(), -1, stream);
        //b.sigma1_2 -= b.mu1_2; - 這會導致額外的資料傳輸操作
 
gauss->apply(b.I2_2, b.sigma2_2, stream);
cuda::subtract(b.sigma2_2, b.mu2_2, b.sigma2_2, cuda::GpuMat(), -1, stream);
        //b.sigma2_2 -= b.mu2_2;
 
gauss->apply(b.I1_I2, b.sigma12, stream);
cuda::subtract(b.sigma12, b.mu1_mu2, b.sigma12, cuda::GpuMat(), -1, stream);
        //b.sigma12 -= b.mu1_mu2;
 
        //這裡也會因為呼叫operator*(Scalar, Mat)而導致額外的資料傳輸
cuda::multiply(b.mu1_mu2, 2, b.t1, 1, -1, stream); //b.t1 = 2 * b.mu1_mu2 + C1;
cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
cuda::multiply(b.sigma12, 2, b.t2, 1, -1, stream); //b.t2 = 2 * b.sigma12 + C2;
cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -12, stream);
 
cuda::multiply(b.t1, b.t2, b.t3, 1, -1, stream); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
cuda::add(b.mu1_2, b.mu2_2, b.t1, cuda::GpuMat(), -1, stream);
cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
 
cuda::add(b.sigma1_2, b.sigma2_2, b.t2, cuda::GpuMat(), -1, stream);
cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -1, stream);
 
 
cuda::multiply(b.t1, b.t2, b.t1, 1, -1, stream); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
cuda::divide(b.t3, b.t1, b.ssim_map, 1, -1, stream); // ssim_map = t3./t1;
 
stream.waitForCompletion();
 
        Scalar s = cuda::sum(b.ssim_map, b.buf);
mssim.val[i] = s.val[0] / (b.ssim_map.rows * b.ssim_map.cols);
 
    }
    return mssim;
}

如何實現？——GPU

如上所示，每種操作我們都有三種類型的函式。一個用於CPU，兩個用於GPU。我為GPU製作了兩個版本的原因是為了說明，簡單地將CPU程式碼移植到GPU上實際上可能會使其變慢。如果你想獲得一些效能提升，你需要記住一些規則，我將在後面詳細介紹。

GPU模組的開發旨在儘可能地與CPU對應版本相似。這使得移植過程更容易。在編寫任何程式碼之前，你需要做的第一件事是將GPU模組連結到你的專案，幷包含該模組的標頭檔案。GPU的所有函式和資料結構都位於cv名稱空間的gpu子名稱空間中。你可以透過use namespace關鍵字將其新增到預設名稱空間，或者透過cv::在所有地方顯式標記它以避免混淆。我將採用後者。

#include <opencv2/gpu.hpp> // GPU結構和方法

GPU是“圖形處理單元”的縮寫。它最初是為了渲染圖形場景而構建的。這些場景在某種程度上建立在大量資料之上。然而，這些資料並非以順序方式相互依賴，並且可以進行並行處理。因此，GPU將包含多個較小的處理單元。這些不是最先進的處理器，在一對一測試中，它們會落後於CPU。然而，其優勢在於其數量。近年來，將GPU這些大規模並行能力用於非圖形場景（也包括渲染）的趨勢日益增長。這催生了圖形處理單元上的通用計算（GPGPU）。

GPU有自己的記憶體。當你使用OpenCV從硬碟讀取資料到Mat物件時，這發生在你的系統記憶體中。CPU透過其快取直接處理這些資料，但GPU不能。它必須將計算所需的資訊從系統記憶體傳輸到自己的記憶體中。這是透過上傳過程完成的，並且耗時。最終，結果必須下載回系統記憶體，供CPU檢視和使用。不建議將小型函式移植到GPU，因為上傳/下載時間將大於並行執行所獲得的收益。

Mat物件僅儲存在系統記憶體（或CPU快取）中。要將OpenCV矩陣獲取到GPU，你需要使用其GPU對應項cv::cuda::GpuMat。它的工作方式與Mat類似，但僅限於2D，並且其函式不返回引用（不能混合GPU引用和CPU引用）。要將Mat物件上傳到GPU，你需要建立類的例項後呼叫upload函式。要下載，你可以直接賦值給Mat物件或使用download函式。

Mat I1; // 主記憶體項 - 例如，用imread讀取影像到此
gpu::GpuMat gI; // GPU矩陣 - 目前為空
gI1.upload(I1); // 將資料從系統記憶體上傳到GPU記憶體
 
I1 = gI1; // 下載，gI1.download(I1) 也可以

一旦你的資料上傳到GPU記憶體中，你就可以呼叫OpenCV中支援GPU的函式。大多數函式與CPU上的函式名稱相同，區別在於它們只接受GpuMat輸入。

另一點需要記住的是，並非所有通道數的影像都能在GPU上高效地實現演算法。通常，我發現GPU影像的輸入影像需要是單通道或四通道的，並且元素大小為char或float型別。抱歉，GPU不支援double型別。對於某些函式，傳入其他型別的物件會導致丟擲異常，並在錯誤輸出上顯示錯誤訊息。文件在大多數地方詳細說明了輸入所接受的型別。如果你有三通道影像作為輸入，你可以做兩件事：要麼新增一個新通道（並使用char元素），要麼拆分影像併為每個影像呼叫函式。第一種方法不建議使用，因為它會浪費記憶體。

對於某些函式，元素的（相鄰項的）位置不重要時，快速解決方案是將其重塑為單通道影像。PSNR實現就是這種情況，因為對於absdiff方法，鄰居的值不重要。然而，對於GaussianBlur，這不是一個選項，因此SSIM需要使用split方法。有了這些知識，你可以製作可行的GPU程式碼（就像我的GPU程式碼一樣）並執行它。你會驚訝地發現它可能比你的CPU實現慢。

最佳化

這是因為你將記憶體分配和資料傳輸的成本拋諸腦後。在GPU上，這個成本非常高。另一種最佳化可能性是藉助cv::cuda::Stream引入非同步OpenCV GPU呼叫。

GPU上的記憶體分配成本相當大。因此，如果可能的話，儘可能少地分配新記憶體。如果你建立一個打算多次呼叫的函式，那麼一個好主意是在第一次呼叫時只分配函式的所有區域性引數。為此，你可以建立一個包含所有將使用的區域性變數的資料結構。例如，在PSNR的情況下，這些變數是：
struct BufferPSNR // 最佳化的GPU版本

{ // 在GPU上資料分配非常昂貴。使用緩衝區解決：分配一次，以後重複使用。

gpu::GpuMat gI1, gI2, gs, t1,t2;

gpu::GpuMat buf;

};

然後在主程式中建立它的一個例項
BufferPSNR bufferPSNR;

最後，每次呼叫時都將其傳遞給函式
double getPSNR_GPU_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)

現在你可以這樣訪問這些區域性引數：b.gI1、b.buf等等。只有當新矩陣的大小與前一個不同時，GpuMat才會重新分配記憶體。
避免不必要的功能資料傳輸。一旦你轉向GPU，任何小的資料傳輸都會變得非常重要。因此，如果可能，請進行所有就地計算（換句話說，不要建立新的記憶體物件——原因在上一節中已解釋）。例如，儘管算術運算可能更容易用一行公式表達，但它會更慢。在SSIM的情況下，我需要計算：
b.t1 = 2 * b.mu1_mu2 + C1;

儘管上面的呼叫會成功，但請注意，其中存在一個隱藏的資料傳輸。在進行加法之前，它需要將乘法結果儲存在某個地方。因此，它會在後臺建立一個區域性矩陣，將C1值新增到其中，最後將其賦值給t1。為了避免這種情況，我們使用gpu函式而不是算術運算子：
gpu::multiply(b.mu1_mu2, 2, b.t1); //b.t1 = 2 * b.mu1_mu2 + C1;

gpu::add(b.t1, C1, b.t1);
使用非同步呼叫（cv::cuda::Stream）。預設情況下，每當你呼叫一個GPU函式時，它都會等待呼叫完成並返回結果。然而，可以進行非同步呼叫，這意味著它會呼叫操作執行，為演算法進行昂貴的資料分配，然後立即返回。現在你可以呼叫另一個函式，如果你願意的話。對於MSSIM，這是一個小的最佳化點。在我們的預設實現中，我們將影像拆分為通道，併為每個通道呼叫GPU函式。透過流可以實現小程度的並行化。透過使用流，我們可以在GPU已經在執行給定方法時進行資料分配和上傳操作。例如，我們需要上傳兩張影像。我們將它們一個接一個地排隊，然後呼叫處理它們的函式。函式將等待上傳完成，但在此期間，它會為下一個要執行的函式進行輸出緩衝區分配。
gpu::Stream stream;

stream.enqueueConvert(b.gI1, b.t1, CV_32F); // 上傳

gpu::split(b.t1, b.vI1, stream); // 方法（將流作為最後一個引數傳遞）。

gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2

CV_32F
#define CV_32F
Definition interface.h:78

結果和結論

在配備低端NVIDIA GT220M的Intel P8700筆記型電腦CPU上，效能資料如下：

PSNR CPU時間（10次執行平均）：41.4122毫秒。結果：19.2506
PSNR GPU時間（10次執行平均）：158.977毫秒。結果：19.2506
GPU最佳化初始呼叫：31.3418毫秒。結果：19.2506
PSNR GPU最佳化時間（/ 10次執行）：24.8171毫秒。結果：19.2506
 
MSSIM CPU時間（10次執行平均）：484.343毫秒。結果為B0.890964 G0.903845 R0.936934
MSSIM GPU時間（10次執行平均）：745.105毫秒。結果為B0.89922 G0.909051 R0.968223
MSSIM GPU初始呼叫時間：357.746毫秒。結果為B0.890964 G0.903845 R0.936934
MSSIM GPU最佳化時間（/ 10次執行）：203.091毫秒。結果為B0.890964 G0.903845 R0.936934

在這兩種情況下，我們都實現了與CPU實現相比近100%的效能提升。這可能正是你的應用程式工作所需的改進。你可以在YouTube上觀看此功能的執行時例項。