GPU上的相似度檢測(PNSR 和 SSIM)

學習目標

在?OpenCV的視頻輸入和相似度測量?教程中我們已經學習了檢測兩幅圖像相似度的兩種方法：PSNR和SSIM。正如我們所看到的，執行這些算法需要相當長的計算時間，其中SSIM（結構相似度）的算法代價相當高昂。然而，OpenCV現在支持Nvidia的CUDA加速，如果你有一塊支持CUDA的的Nvidia顯卡。您可以將算法改為使用GPU計算從而大幅提高效率。

本教程將提供一個很好的示例來演示如何讓OpenCV來使用GPU這些操作。當然在這之前，你應該了解一下如何操作core，highgui 和imgproc 模塊。因此我們的目標是:

• GPU與CPU有什么不同?
• 為PSNR何SSIM編寫GPU 代碼
• 優化代碼以獲得最佳性能

源代碼

你也可以在OPENCV源庫文件夾中?samples/cpp/tutorial_code/gpu/gpu-basics-similarity/gpu-basics-similarity?找到源代碼以及一些視頻文件或者?從這里下載?。 完整的源碼很長（因為其中包含了通過命令行參數進行應用程序和性能測量等無關代碼）。因此在這里只出現一部分關鍵代碼。

如果兩個輸入圖像的相似，PSNR 返回一個浮點數在30和50之間，(數值越高，相符程度越高).

123456789 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

double getPSNR(const Mat& I1, const Mat& I2) {Mat s1; absdiff(I1, I2, s1); // |I1 - I2|s1.convertTo(s1, CV_32F); // cannot make a square on 8 bitss1 = s1.mul(s1); // |I1 - I2|^2Scalar s = sum(s1); // sum elements per channeldouble sse = s.val[0] + s.val[1] + s.val[2]; // sum channelsif( sse <= 1e-10) // for small values return zeroreturn 0;else{double mse =sse /(double)(I1.channels() * I1.total());double psnr = 10.0*log10((255*255)/mse);return psnr;} }double getPSNR_GPU_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);gpu::absdiff(b.t1.reshape(1), b.t2.reshape(1), b.gs);gpu::multiply(b.gs, b.gs, b.gs);double sse = gpu::sum(b.gs, b.buf)[0];if( sse <= 1e-10) // for small values return zeroreturn 0;else{double mse = sse /(double)(I1.channels() * I1.total());double psnr = 10.0*log10((255*255)/mse);return psnr;} }struct BufferPSNR // Optimized GPU versions { // Data allocations are very expensive on GPU. Use a buffer to solve: allocate once reuse later.gpu::GpuMat gI1, gI2, gs, t1,t2;gpu::GpuMat buf; };double getPSNR_GPU(const Mat& I1, const Mat& I2) {gpu::GpuMat gI1, gI2, gs, t1,t2; gI1.upload(I1);gI2.upload(I2);gI1.convertTo(t1, CV_32F);gI2.convertTo(t2, CV_32F);gpu::absdiff(t1.reshape(1), t2.reshape(1), gs); gpu::multiply(gs, gs, gs);Scalar s = gpu::sum(gs);double sse = s.val[0] + s.val[1] + s.val[2];if( sse <= 1e-10) // for small values return zeroreturn 0;else{double mse =sse /(double)(gI1.channels() * I1.total());double psnr = 10.0*log10((255*255)/mse);return psnr;} }

SSIM 返回圖像的結構相似度指標。這是一個在0-1之間的浮點數(越接近1，相符程度越高), 這個算法對每個通道有一個值，所以最終會產生一個OpenCV的?Scalar?數據結構:

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208

}Scalar getMSSIM( const Mat& i1, const Mat& i2) { const double C1 = 6.5025, C2 = 58.5225;/***************************** INITS **********************************/int d = CV_32F;Mat I1, I2; i1.convertTo(I1, d); // cannot calculate on one byte large valuesi2.convertTo(I2, d); Mat I2_2 = I2.mul(I2); // I2^2Mat I1_2 = I1.mul(I1); // I1^2Mat I1_I2 = I1.mul(I2); // I1 * I2/*************************** END INITS **********************************/Mat mu1, mu2; // PRELIMINARY COMPUTINGGaussianBlur(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;/ FORMULA 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 = average of ssim mapreturn mssim; }Scalar getMSSIM_GPU( const Mat& i1, const Mat& i2) { const float C1 = 6.5025f, C2 = 58.5225f;/***************************** INITS **********************************/gpu::GpuMat gI1, gI2, gs1, t1,t2; gI1.upload(i1);gI2.upload(i2);gI1.convertTo(t1, CV_MAKE_TYPE(CV_32F, gI1.channels()));gI2.convertTo(t2, CV_MAKE_TYPE(CV_32F, gI2.channels()));vector<gpu::GpuMat> vI1, vI2; gpu::split(t1, vI1);gpu::split(t2, vI2);Scalar mssim;for( int i = 0; i < gI1.channels(); ++i ){gpu::GpuMat I2_2, I1_2, I1_I2; gpu::multiply(vI2[i], vI2[i], I2_2); // I2^2gpu::multiply(vI1[i], vI1[i], I1_2); // I1^2gpu::multiply(vI1[i], vI2[i], I1_I2); // I1 * I2/*************************** END INITS **********************************/gpu::GpuMat mu1, mu2; // PRELIMINARY COMPUTINGgpu::GaussianBlur(vI1[i], mu1, Size(11, 11), 1.5);gpu::GaussianBlur(vI2[i], mu2, Size(11, 11), 1.5);gpu::GpuMat mu1_2, mu2_2, mu1_mu2; gpu::multiply(mu1, mu1, mu1_2); gpu::multiply(mu2, mu2, mu2_2); gpu::multiply(mu1, mu2, mu1_mu2); gpu::GpuMat sigma1_2, sigma2_2, sigma12; gpu::GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);sigma1_2 -= mu1_2;gpu::GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);sigma2_2 -= mu2_2;gpu::GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);sigma12 -= mu1_mu2;/ FORMULA gpu::GpuMat t1, t2, t3; t1 = 2 * mu1_mu2 + C1; t2 = 2 * sigma12 + C2; gpu::multiply(t1, t2, t3); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))t1 = mu1_2 + mu2_2 + C1; t2 = sigma1_2 + sigma2_2 + C2; gpu::multiply(t1, t2, t1); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))gpu::GpuMat ssim_map;gpu::divide(t3, t1, ssim_map); // ssim_map = t3./t1;Scalar s = gpu::sum(ssim_map); mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols);}return mssim; } struct BufferMSSIM // Optimized GPU versions { // Data allocations are very expensive on GPU. Use a buffer to solve: allocate once reuse later.gpu::GpuMat gI1, gI2, gs, t1,t2;gpu::GpuMat I1_2, I2_2, I1_I2;vector<gpu::GpuMat> vI1, vI2;gpu::GpuMat mu1, mu2; gpu::GpuMat mu1_2, mu2_2, mu1_mu2; gpu::GpuMat sigma1_2, sigma2_2, sigma12; gpu::GpuMat t3; gpu::GpuMat ssim_map;gpu::GpuMat buf; }; Scalar getMSSIM_GPU_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b) { int cn = i1.channels();const float C1 = 6.5025f, C2 = 58.5225f;/***************************** INITS **********************************/b.gI1.upload(i1);b.gI2.upload(i2);gpu::Stream stream;stream.enqueueConvert(b.gI1, b.t1, CV_32F);stream.enqueueConvert(b.gI2, b.t2, CV_32F); gpu::split(b.t1, b.vI1, stream);gpu::split(b.t2, b.vI2, stream);Scalar mssim;for( int i = 0; i < b.gI1.channels(); ++i ){ gpu::multiply(b.vI2[i], b.vI2[i], b.I2_2, stream); // I2^2gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2gpu::multiply(b.vI1[i], b.vI2[i], b.I1_I2, stream); // I1 * I2gpu::GaussianBlur(b.vI1[i], b.mu1, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream);gpu::GaussianBlur(b.vI2[i], b.mu2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream);gpu::multiply(b.mu1, b.mu1, b.mu1_2, stream); gpu::multiply(b.mu2, b.mu2, b.mu2_2, stream); gpu::multiply(b.mu1, b.mu2, b.mu1_mu2, stream); gpu::GaussianBlur(b.I1_2, b.sigma1_2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream);gpu::subtract(b.sigma1_2, b.mu1_2, b.sigma1_2, stream);//b.sigma1_2 -= b.mu1_2; - This would result in an extra data transfer operationgpu::GaussianBlur(b.I2_2, b.sigma2_2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream);gpu::subtract(b.sigma2_2, b.mu2_2, b.sigma2_2, stream);//b.sigma2_2 -= b.mu2_2;gpu::GaussianBlur(b.I1_I2, b.sigma12, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream);gpu::subtract(b.sigma12, b.mu1_mu2, b.sigma12, stream);//b.sigma12 -= b.mu1_mu2;//here too it would be an extra data transfer due to call of operator*(Scalar, Mat)gpu::multiply(b.mu1_mu2, 2, b.t1, stream); //b.t1 = 2 * b.mu1_mu2 + C1; gpu::add(b.t1, C1, b.t1, stream);gpu::multiply(b.sigma12, 2, b.t2, stream); //b.t2 = 2 * b.sigma12 + C2; gpu::add(b.t2, C2, b.t2, stream); gpu::multiply(b.t1, b.t2, b.t3, stream); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))gpu::add(b.mu1_2, b.mu2_2, b.t1, stream);gpu::add(b.t1, C1, b.t1, stream);gpu::add(b.sigma1_2, b.sigma2_2, b.t2, stream);gpu::add(b.t2, C2, b.t2, stream);gpu::multiply(b.t1, b.t2, b.t1, stream); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2)) gpu::divide(b.t3, b.t1, b.ssim_map, stream); // ssim_map = t3./t1;stream.waitForCompletion();Scalar s = gpu::sum(b.ssim_map, b.buf); mssim.val[i] = s.val[0] / (b.ssim_map.rows * b.ssim_map.cols);}return mssim; }

如何在GPU上實現相同算法?

你可以總共有三種類型的函數來針對每個操作。其中包括一個CPU和2個GPU函數。之所以做兩個GPU函數是為了說明一個常識：簡單的移植你的CPU程序到GPU反而會讓速度變慢。如果你想要一些性能上的改善，你將需牢記住一些規則，我在后面詳細說明。

開發的GPU模塊盡可能多的與CPU對應，這樣才能方便移植。在你編寫任何代碼之間要做的第一件事是連GPU模塊到你的項目中，包括模塊的頭文件。所有GPU的函數和數據結構都在以?cv?命名空間的子空間?gpu?內。你可以將其添加為默認?use namespace?關鍵字, 也可以顯示的通過cv::來避免混亂。我在后面會做。

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

GPU代表**圖形處理單元**。最開始是為渲染各種圖形場景而建立，這些場景是基于大量的矢量數據建立的。由于矢量圖形的特殊性，數據不需要以串行的方式一步一步執行的，而是并行的方式一次性渲染大量的數據。從GPU的結構上來說，不像CPU基于數個寄存器和高速指令集，GPU一般有數百個較小的處理單元。這些處理單元每一個都比CPU的核心慢很多很多。然而，它的力量在于它的數目眾多，能夠同時進行大量運算。在過去的幾年中已經有越來越多的人常識用GPU來進行圖形渲染以外的工作。這就產生了GPGPU（general-purpose computation on graphics processing units)的概念。

圖像處理器有它自己的內存，一般稱呼為顯存。當你從硬盤驅動器讀數據并產生一個?Mat?對象的時候，數據是放在普通內存當中的（由CPU掌管）CPU可以直接操作內存, 然而GPU不能這樣，對于電腦來說GPU只是個外設，它只能操作它自己的顯存，當計算時，需要先讓CPU將用于計算的信息轉移到GPU掌管的顯存上。 這是通過一個上傳過程完成，需要比訪問內存多很多的時間。而且最終的計算結果將要回傳到你的系統內存處理器才能和其他代碼交互，由于傳輸的代價高昂，所以注定移植太小的函數到GPU上并不會提高效率。

Mat對象僅僅存儲在內存或者CPU緩存中。為了得到一個GPU能直接訪問的opencv 矩陣你必須使用GPU對象?GpuMat?。它的工作方式類似于2維 Mat，唯一的限制是你不能直接引用GPU函數。（因為它們本質上是完全不同的代碼，不能混合引用）。要傳輸*Mat*對象到*GPU*上并創建GpuMat時需要調用上傳函數，在回傳時，可以使用簡單的重載賦值操作符或者調用下載函數。

Mat I1; // 內存對象，可以用imread來創建 gpu::GpuMat gI; // GPU 矩陣 - 現在為空 gI1.upload(I1); //將內存數據上傳到顯存中I1 = gI1; //回傳, gI1.download(I1) 也可以

記住：一旦你的數據被傳到GPU顯存中，你就只能調用GPU函數來操作，大部分gpu函數名字與原來CPU名字保持一致，不同的是他們只接收?GpuMat?輸入。可以在:?在線文檔?中找到一些說明，或者查閱OpenCV參考手冊。

另外要記住的是：并非所有的圖像類型都可以用于GPU算法中。很多時候，GPU函數只能接收單通道或者4通道的uchar或者float類型（CV_8UC1,CV_8UC4,CV_32FC1,CV_32FC4），GPU函數不支持雙精度。如果你試圖向這些函數傳入非指定類型的數據時，這些函數會拋出異常或者輸出錯誤信息。這些函數的文檔里大都說明了接收的數據類型。如果函數只接收4通道或者單通道的圖像而你的數據類型剛好是3通道的話，你能做的就是兩件事：1.增加一個新的通道（使用char類型）或 2.將圖像的三個通道切分為三個獨立的圖像，為每個圖像調用該函數。不推薦使用第一個方法，因為有些浪費內存。

對于不需要在意元素的坐標位置的某些函數，快速解決辦法就是將它直接當作一個單通道圖像處理。這種方法適用于PSNR中需要調用的absdiff方法。然而，對于GaussianBlur就不能這么用，對于SSIM，就需要使用分離通道的方法。知道這些知識就已經可以讓你的GPU開始運行代碼。但是你也可能發現GPU“加速”后的代碼仍然可能會低于你的CPU執行速度。

優化

因為大量的時間都耗費在傳輸數據到顯存這樣的步驟了，大量的計算時間被消耗在內存分配和傳輸上，GPU的高計算能力根本發揮不出來，我們可以利用?gpu::Stream?來進行異步傳輸，從而節約一些時間.

＃. GPU上的顯存分配代價是相當大的。因此如果可能的話，應該盡可能少的分配新的內存。如果你需要創建一個調用多次的函數，就應該在一開始分配全部的內存，而且僅僅分配一次。在第一次調用的時候可以創建一個結構體包含所有將使用局部變量。 對于PSNR例子:

struct BufferPSNR // 優化版本{ // 基于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中，任何小的數據傳輸都是一次大的開銷，所以應該盡量避免不必要的數據傳輸。因此應該盡可能的在現有對象上計算（換句話說，不創造新的顯存對象，上面已經解釋過）。例如，雖然使用算術表達式會讓一個公式更淺顯易懂但是在GPU代碼里，這個步驟是很緩慢的。例如在SSIM例子中需要計算：

b.t1 = 2 * b.mu1_mu2 + C1;

在這個表達式的調用過程中必然會有一個隱性內存分配過程，調用乘法的結果必然要用一個臨時對象儲存，然后才能和*C1*相加。如果直接用表達式寫GPU代碼，就必然會創建一個臨時變量來儲存乘積的矩陣，然后加上?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);

使用異步調用 (the?gpu::Stream)。默認情況下，無論何時你調用一個GPU函數，系統都將等待調用完成并返回后繼結果，這就意味著在GPU計算的時候CPU沒干活。如果使異步調用，就會意味著它將調用后立即返回，并在后臺異步執行，也就是說雖然CPU需要的數據還沒從GPU返回，但是并不妨礙CPU進行其它的計算或者調用另一個函數。這是一個小的MSSIM可優化點。在我們的默認實現中，我們將圖像分裂成多個通道，然后為每個通道調用GPU函數。這本來應該是一個平行操作，我們可以建立一個流，通過使用流我們可以異步進行數據分配、傳輸等各種操作，并且給GPU執行操作方法以便GPU在合適的時候自動調用。例如我們需要傳輸兩幅圖像。我們依次調用現成的函數處理它們。該函數還是要等到上傳數據完畢才會開始工作，但是在執行第二個函數的時候，原來的輸出緩存將直接作為第二個操作的輸出參數而不需要傳輸過程。

gpu::Stream stream;stream.enqueueConvert(b.gI1, b.t1, CV_32F); // 上傳gpu::split(b.t1, b.vI1, stream); // 分離方法(在最后的參數里傳遞stream). gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2

結果和結論

以下是一個Intel P8700的筆記本CPU，與一個低端的Nvidia GT220M的性能比較，結果如下:

Time of PSNR CPU (averaged for 10 runs): 41.4122 milliseconds. With result of: 19.2506 Time of PSNR GPU (averaged for 10 runs): 158.977 milliseconds. With result of: 19.2506 Initial call GPU optimized: 31.3418 milliseconds. With result of: 19.2506 Time of PSNR GPU OPTIMIZED ( / 10 runs): 24.8171 milliseconds. With result of: 19.2506Time of MSSIM CPU (averaged for 10 runs): 484.343 milliseconds. With result of B0.890964 G0.903845 R0.936934 Time of MSSIM GPU (averaged for 10 runs): 745.105 milliseconds. With result of B0.89922 G0.909051 R0.968223 Time of MSSIM GPU Initial Call 357.746 milliseconds. With result of B0.890964 G0.903845 R0.936934 Time of MSSIM GPU OPTIMIZED ( / 10 runs): 203.091 milliseconds. With result of B0.890964 G0.903845 R0.936934

在這個案例中，使用GPU加速比單純使用CPU提高了近100%的性能。這對構建高速實時系統甚有幫助。你可以在并不存在于天朝人民概念中的`YouTube 上看到實例視頻 ｀<https://www.youtube.com/watch?v=3_ESXmFlnvY>`_.

from:?http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/gpu/table_of_content_gpu/table_of_content_gpu.html#table-of-content-gpu

總結

以上是生活随笔為你收集整理的OpenCV之gpu 模块. 使用GPU加速的计算机视觉：GPU上的相似度检测(PNSR 和 SSIM)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。