ARM NEON的快速高斯模糊圖像濾波器

Question

我正在試圖製作高斯模糊圖像濾波器的移動快速版本。

我讀過其他的問題，比如：Fast Gaussian blur on unsigned char image- ARM Neon Intrinsics- iOS Dev

對於我的目的，我只有一個固定大小（7×7）固定西格馬（2）高斯濾波器需要。因此，在針對ARM NEON進行優化之前，我正在C++中實現一維高斯內核，並直接在移動環境（Android with NDK）中將性能與OpenCV gaussianBlur（）方法進行比較。這樣就會導致更簡單的代碼進行優化。

但是結果是我的實現比OpenCV4Android版本慢了10倍。我讀過OpenCV4 Tegra已經優化了GaussianBlur的實現，但我不認爲標準的OpenCV4Android有這樣的優化，那麼爲什麼我的代碼如此之慢？

這裏是我的執行（注：近邊界應用過濾器時reflect101用於像素反射）：

Mat myGaussianBlur(Mat src){ 
    Mat dst(src.rows, src.cols, CV_8UC1); 
    Mat temp(src.rows, src.cols, CV_8UC1); 
    float sum, x1, y1; 

    // coefficients of 1D gaussian kernel with sigma = 2 
    double coeffs[] = {0.06475879783, 0.1209853623, 0.1760326634, 0.1994711402, 0.1760326634, 0.1209853623, 0.06475879783}; 
    //Normalize coeffs 
    float coeffs_sum = 0.9230247873f; 
    for (int i = 0; i < 7; i++){ 
     coeffs[i] /= coeffs_sum; 
    } 

    // filter vertically 
    for(int y = 0; y < src.rows; y++){ 
     for(int x = 0; x < src.cols; x++){ 
      sum = 0.0; 
      for(int i = -3; i <= 3; i++){ 
       y1 = reflect101(src.rows, y - i); 
       sum += coeffs[i + 3]*src.at<uchar>(y1, x); 
      } 
      temp.at<uchar>(y,x) = sum; 
     } 
    } 

    // filter horizontally 
    for(int y = 0; y < src.rows; y++){ 
     for(int x = 0; x < src.cols; x++){ 
      sum = 0.0; 
      for(int i = -3; i <= 3; i++){ 
       x1 = reflect101(src.rows, x - i); 
       sum += coeffs[i + 3]*temp.at<uchar>(y, x1); 
      } 
      dst.at<uchar>(y,x) = sum; 
     } 
    } 

    return dst; 
} 

Answer 1

這是實現@保羅R的所有建議和@ SH1之後的代碼，總結如下：

1）只使用整數運算（精確適量）

2）中的值相加在應用乘法之前距離掩模中心相同距離的像素，以減少乘法的次數。

3）僅適用水平濾波器，以利用存儲的由矩陣

4）單獨的週期從周圍與圖像內的邊緣的行不使反射功能的不必要呼叫。我完全刪除了反射的功能，包括它們沿着邊緣的循環。另外，作爲個人觀察，爲了在不調用（緩慢）函數「round」或「cvRound」的情況下改善舍入，我已經將臨時和最終像素結果都加到0.5f（= 32768 in integer精度）來減少與OpenCV相比的誤差/差異。

現在性能比OpenCV慢15到6倍左右。

但是，得到的矩陣與OpenCV的高斯模糊所得到的矩陣並不完全相同。這不是由於算術長度（足夠）以及刪除錯誤仍然存​​在。請注意，這是兩個版本產生的矩陣之間的最小差異，在像素強度的0和2（絕對值）之間。 OpenCV使用的係數與getGaussianKernel具有相同大小和sigma的係數相同。

Mat myGaussianBlur(Mat src){ 

Mat dst(src.rows, src.cols, CV_8UC1); 
Mat temp(src.rows, src.cols, CV_8UC1); 
int sum; 
int x1; 

double coeffs[] = {0.070159, 0.131075, 0.190713, 0.216106, 0.190713, 0.131075, 0.070159}; 
int coeffs_i[7] = { 0 }; 
for (int i = 0; i < 7; i++){ 
     coeffs_i[i] = (int)(coeffs[i] * 65536); //65536 
} 

// filter horizontally - inside the image 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = src.ptr<uchar>(y); 
    for(int x = 3; x < (src.cols - 3); x++){ 
     sum = ptr[x] * coeffs_i[3]; 
     for(int i = -3; i < 0; i++){ 
      int tmp = ptr[x+i] + ptr[x-i]; 
      sum += coeffs_i[i + 3]*tmp; 
     } 
     temp.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 
// filter horizontally - edges - needs reflect 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = src.ptr<uchar>(y); 
    for(int x = 0; x <= 2; x++){ 
     sum = 0; 
     for(int i = -3; i <= 3; i++){ 
      x1 = x + i; 
      if(x1 < 0){ 
       x1 = -x1; 
      } 
      sum += coeffs_i[i + 3]*ptr[x1]; 
     } 
     temp.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = src.ptr<uchar>(y); 
    for(int x = (src.cols - 3); x < src.cols; x++){ 
     sum = 0; 
     for(int i = -3; i <= 3; i++){ 
      x1 = x + i; 
      if(x1 >= src.cols){ 
       x1 = 2*src.cols - x1 - 2; 
      } 
      sum += coeffs_i[i + 3]*ptr[x1]; 
     } 
     temp.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 

// transpose to apply again horizontal filter - better cache data locality 
transpose(temp, temp); 

// filter horizontally - inside the image 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = temp.ptr<uchar>(y); 
    for(int x = 3; x < (src.cols - 3); x++){ 
     sum = ptr[x] * coeffs_i[3]; 
     for(int i = -3; i < 0; i++){ 
      int tmp = ptr[x+i] + ptr[x-i]; 
      sum += coeffs_i[i + 3]*tmp; 
     } 
     dst.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 
// filter horizontally - edges - needs reflect 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = temp.ptr<uchar>(y); 
    for(int x = 0; x <= 2; x++){ 
     sum = 0; 
     for(int i = -3; i <= 3; i++){ 
      x1 = x + i; 
      if(x1 < 0){ 
       x1 = -x1; 
      } 
      sum += coeffs_i[i + 3]*ptr[x1]; 
     } 
     dst.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 
for(int y = 0; y < src.rows; y++){ 
    uchar *ptr = temp.ptr<uchar>(y); 
    for(int x = (src.cols - 3); x < src.cols; x++){ 
     sum = 0; 
     for(int i = -3; i <= 3; i++){ 
      x1 = x + i; 
      if(x1 >= src.cols){ 
       x1 = 2*src.cols - x1 - 2; 
      } 
      sum += coeffs_i[i + 3]*ptr[x1]; 
     } 
     dst.at<uchar>(y,x) = (sum + 32768)/65536; 
    } 
} 

transpose(dst, dst); 

return dst; 
} 

Answer 2

如果這是專門爲8幅圖像，那麼你真的不想浮點係數，尤其不是雙精度。你也不想爲x1，y1使用浮點數。您應該使用整數作爲座標，並且您可以使用係數的固定點（即整數）來保留整數域中的所有濾波算術，例如，

Mat myGaussianBlur(Mat src){ 
    Mat dst(src.rows, src.cols, CV_8UC1); 
    Mat temp(src.rows, src.cols, CV_16UC1); // <<< 
    int sum, x1, y1; // <<< 

    // coefficients of 1D gaussian kernel with sigma = 2 
    double coeffs[] = {0.06475879783, 0.1209853623, 0.1760326634, 0.1994711402, 0.1760326634, 0.1209853623, 0.06475879783}; 
    int coeffs_i[7] = { 0 }; // <<< 
    //Normalize coeffs 
    float coeffs_sum = 0.9230247873f; 
    for (int i = 0; i < 7; i++){ 
     coeffs_i[i] = (int)(coeffs[i]/coeffs_sum * 256); // <<< 
    } 

    // filter vertically 
    for(int y = 0; y < src.rows; y++){ 
     for(int x = 0; x < src.cols; x++){ 
      sum = 0; // <<< 
      for(int i = -3; i <= 3; i++){ 
       y1 = reflect101(src.rows, y - i); 
       sum += coeffs_i[i + 3]*src.at<uchar>(y1, x); // <<< 
      } 
      temp.at<uchar>(y,x) = sum; 
     } 
    } 

    // filter horizontally 
    for(int y = 0; y < src.rows; y++){ 
     for(int x = 0; x < src.cols; x++){ 
      sum = 0; // <<< 
      for(int i = -3; i <= 3; i++){ 
       x1 = reflect101(src.rows, x - i); 
       sum += coeffs_i[i + 3]*temp.at<uchar>(y, x1); // <<< 
      } 
      dst.at<uchar>(y,x) = sum/(256 * 256); // <<< 
     } 
    } 

    return dst; 
} 

Answer 3

這個問題的一大部分是，這個算法過於精確，正如@PaulR指出的那樣。通常最好保持你的係數表不比你的數據更精確。在這種情況下，由於您似乎正在處理uchar數據，因此您可以使用大致8位的係數表。

保持這些權重很小將在您的NEON實施中特別重要，因爲您的算術範圍越窄，您可以一次處理的路線就越多。

除此之外，突出的第一個主要放緩是在主循環內具有圖像邊緣反射代碼。這會使大部分工作效率降低，因爲在這種情況下通常不需要做任何特別的事情。

如果您在邊緣附近使用特殊版本的循環，並且在安全的情況下使用不稱爲reflect101()函數的簡化內部循環，它可能會更好。

其次（與原型代碼更相關）是可以在應用加權函數之前將窗口的翅膀一起添加，因爲表格在兩側包含相同的係數。

sum = src.at<uchar>(y1, x) * coeffs[3]; 
for(int i = -3; i < 0; i++) { 
    int tmp = src.at<uchar>(y + i, x) + src.at<uchar>(y - i, x); 
    sum += coeffs[i + 3] * tmp; 
} 

這可以爲每個像素節省6次乘法，而且這是圍繞控制溢出條件進行一些其他優化的一個步驟。

然後還有一些與內存系統有關的其他問題。

兩步法原則上是好的，因爲它可以幫助您避免執行大量重新計算。不幸的是，它可以將有用的數據從L1緩存中提取出來，這會讓所有事情變得相當慢。這也意味着，當你將結果寫入內存時，你會量化中間和，這會降低精度。

當您將此代碼轉換爲NEON時，您需要關注的一件事是試圖將您的工作集保留在寄存器文件中，但在完全使用之前不會丟棄計算結果。

當人們做使用兩遍時，通常中間數據被轉置 - 也就是說，一列輸入成爲一行輸出。

這是因爲CPU真的不喜歡跨輸入圖像的多行獲取少量數據。如果您收集一堆水平像素並對其進行過濾，則效率會更高（因爲緩存的工作方式）。如果臨時緩衝區被轉置，那麼第二遍也將收集在一起的一組水平點（這將在原始方向上垂直），並且它再次轉換它的輸出以便它以正確的方式出現。

如果您優化以保持工作集本地化，那麼您可能不需要這種換位技巧，但值得了解，以便您可以爲自己設置健康的基準性能。不幸的是，像這樣的本地化確實會迫使你回到非最優的內存提取，但是更廣泛的數據類型可以減輕懲罰。