通常，我們談的高斯模糊，都知道其是可以行列分離的演算法，現在也有著各種優化演算法實作，而且其速度基本是和引數大小無關的，但是，在我們實際的應用中，我們可能會發現，有至少50%以上的場景中，我們并不需要大半徑的高斯，反而是微小半徑的模糊更有用武之地（比如Canny的預處理、簡單去噪等），因此，小半徑的高斯是否能進一步加速就值的研究，正因為如此，一些商業軟體都提供了類似的功能，比如在halon中，直接的高斯模糊可以用smooth_image實作，但是你在其幫助檔案中搜索gauss關鍵字后，你會發現有以下兩個函式：
gauss_filter — Smooth using discrete gauss functions.
gauss_image — Smooth an image using discrete Gaussian functions.
兩個函式的功能描述基本是一個意思，但是在gauss_image函式的注釋下有這么一條：
gauss_image is obsolete and is only provided for reasons of backward compatibility. New applications should use the operator gauss_filter instead.
即這個函式已經過時，提供他只是為了向前兼容，新的應用建議使用gauss_filter 函式，那我們再來看下halcon中其具體的描述：
Signature
gauss_filter(Image : ImageGauss : Size : )
Description
The operator gauss_filter smoothes images using the discrete Gaussian, a discrete approximation of the Gaussian function,
The smoothing effect increases with increasing filter size. The following filter sizes (Size) are supported (the sigma value of the gauss function is indicated in brackets):
3 (0.600)
5 (1.075)
7 (1.550)
9 (2.025)
11 (2.550)
For border treatment the gray values of the images are reflected at the image borders. Notice that, contrary to the operator gauss_image, the relationship between the filter mask size and its respective value for the sigma parameter is linear.
可見gauss_filter的Size只能取3、5、7、9、11這五個值，括號里給出了對應的sigma值，
這種小半徑的模糊的優化其實在我博客里有講過好幾個這方面的，但是前面講述的基本都是直徑不超過5，半徑不大于2的，比如這里的3和5就可以直接用那種方法處理，
我們先來看看這個權重怎么計算：
簡單的例子，比如Size=3時，其實就是一個3*3的卷積，這個3*3的卷積核可以用下述方式計算出來：
const int Diameter = 3, Radius = 1; float Sum = 0, Delta = 0.600, Weight[Diameter][Diameter]; for (int Y = 0; Y < Diameter; Y++) { for (int X = 0; X < Diameter; X++) { Weight[X][Y] = exp(-((X - Radius) * (X - Radius) + (Y - Radius) * (Y - Radius)) / (2 * Delta * Delta)); Sum += Weight[X][Y]; } } // 0.027682 0.111015 0.027682 // 0.111015 0.445213 0.111015 // 0.027682 0.111015 0.027682 for (int Y = 0; Y < Diameter; Y++) { for (int X = 0; X < Diameter; X++) { Weight[X][Y] /= Sum; } }
注意卷積和要中心對稱，
要計算這個3*3的卷積，可以直接用浮點資料計算，很明顯，這里我們可以直接硬計算，而無需其他什么優化技巧，但是為了不必要的浮點計算，我們很明顯可以把這個卷積核定點話，弄成整數，比如整體乘以16384倍，得到下面的卷積核：
// 454 1819 454 // 1819 7292 1819 // 454 1819 454
這樣，就可以直接借助整數的乘法和一個移位來得到最終的結果值，
另外，還有一個特點，就是借助于SIMD執行還可以是實作一次性進行4個整數的計算，如果在厲害一點，還可以使用_mm_madd_epi16這個特別的SIMD指令，一次性實作8位整數的計算，效率大大的提高，
有興趣的可以找找我以前的博文，
當半徑大于3時，在使用直接卷積就帶來了一定的性能問題，比如直徑為7時，每個點的計算量有49次了，這個時候即使借助于SSE也會發現，其耗時和優化后的任意核的高斯相比已經不具有任何優勢了，當半徑進一步加大時，反而超過了任何核心的優化效果，這個時候我們就需要使用另外一個特性了，即高斯卷積的行列分離特性，我們以9*9的高斯卷積核為例：
使用類似上面的代碼，可以得到這個時候的歸一化的高斯卷積核如下：
// 0.000825 0.001936 0.003562 0.005135 0.005801 0.005135 0.003562 0.001936 0.000825 // 0.001936 0.004545 0.008363 0.012056 0.013619 0.012056 0.008363 0.004545 0.001936 // 0.003562 0.008363 0.015386 0.022181 0.025057 0.022181 0.015386 0.008363 0.003562 // 0.005135 0.012056 0.022181 0.031977 0.036124 0.031977 0.022181 0.012056 0.005135 // 0.005801 0.013619 0.025057 0.036124 0.040809 0.036124 0.025057 0.013619 0.005801 // 0.005135 0.012056 0.022181 0.031977 0.036124 0.031977 0.022181 0.012056 0.005135 // 0.003562 0.008363 0.015386 0.022181 0.025057 0.022181 0.015386 0.008363 0.003562 // 0.001936 0.004545 0.008363 0.012056 0.013619 0.012056 0.008363 0.004545 0.001936 // 0.000825 0.001936 0.003562 0.005135 0.005801 0.005135 0.003562 0.001936 0.000825
我們看到這個卷積核其轉置后的結果和原型一模一樣，這樣的卷積核就具有行列可分離性，我們要得到其可分離的卷積列向量或者行列量，可以通過歸一化其第一行或者第一列的得到，比如將上述核心第一行歸一化后得到：
0.028714 0.067419 0.124039 0.178822 0.202011 0.178822 0.124039 0.067419 0.028714
用matlab驗證下：
>> a=[0.028714 0.067419 0.124039 0.178822 0.202011 0.178822 0.124039 0.067419 0.028714] a = 0.0287 0.0674 0.1240 0.1788 0.2020 0.1788 0.1240 0.0674 0.0287 >> a'*a ans = 0.0008 0.0019 0.0036 0.0051 0.0058 0.0051 0.0036 0.0019 0.0008 0.0019 0.0045 0.0084 0.0121 0.0136 0.0121 0.0084 0.0045 0.0019 0.0036 0.0084 0.0154 0.0222 0.0251 0.0222 0.0154 0.0084 0.0036 0.0051 0.0121 0.0222 0.0320 0.0361 0.0320 0.0222 0.0121 0.0051 0.0058 0.0136 0.0251 0.0361 0.0408 0.0361 0.0251 0.0136 0.0058 0.0051 0.0121 0.0222 0.0320 0.0361 0.0320 0.0222 0.0121 0.0051 0.0036 0.0084 0.0154 0.0222 0.0251 0.0222 0.0154 0.0084 0.0036 0.0019 0.0045 0.0084 0.0121 0.0136 0.0121 0.0084 0.0045 0.0019 0.0008 0.0019 0.0036 0.0051 0.0058 0.0051 0.0036 0.0019 0.0008
和前面計算的核是一致的，
這個時候的策略就需要改變了，不能直接計算，我們分配一個臨時的中國記憶體，考慮到精度問題，建議中間記憶體至少使用short型別，我們對原始資料先進行行方向的一維卷積，并取適當的移位資料，將這個中間結果保留在臨時的記憶體中，然后在對臨時記憶體記性列方向的卷積，保存到目標中，考慮到卷積時邊緣部分會超出邊界，所以還可以使用一個臨時擴展的記憶體，提前把邊界位置的內容設計好并填充進去，計算時，就可以連續訪問了，
其實這里也有兩種選擇，即先只計算那些領域不會超出邊界的中心像素（使用SIMD優化），然后再用普通的C代碼組防邊界溢位的普通演算法，但是測驗發現，這些普通的C代碼的耗時占整體的比例有點夸張了，還不如前面的做個臨時擴展記憶體來的快速和方便，
同樣的道理，水平和垂直方向的一維卷積也應該用定點化來實作，同樣的可借助于_mm_madd_epi16指令，
我們測驗了halcon的gauss_filter 的速度，測驗代碼如下所示：
HalconCpp::HObject hoImage0; HalconCpp::ReadImage(&hoImage0, "D:\\1.bmp"); int max_iter = 100; int multi = 3; timepoint tb; long long tp; static int w[] = { 3, 5, 7, 9, 11, 13, 15, 21, 31, 41, 51 }; static int h[] = { 3, 5, 7, 9, 11, 13, 15, 21, 31, 41, 51 }; HalconCpp::SetSystem("parallelize_operators", "false"); HalconCpp::HObject hoImageT; for (int i = 0; i < 5; i++) { tb = time_now(); for (int k = 0; k < max_iter; ++k) HalconCpp::GaussFilter(hoImage0, &hoImageT, w[i]); tp = time_past(tb); std::cout << "GaussFilter " << w[i] << " use time:" << tp / max_iter/1000 << "ms" << std::endl; }
測驗的結果如下：
為了測驗的公平，我們關閉了halcon的多執行緒優化方面的功能，即使用了如下的陳述句：
HalconCpp::SetSystem("parallelize_operators", "false");
我也對我優化后的演算法進行了速度測驗，主要耗時如下表所示：
和halcon相比，基本在同一個數量級別上，
不過halcon的smooth_image似乎非常的慢，即使我不用其"gauss"引數，同樣的圖片，其耗時如下所示：
呼叫代碼為：
HalconCpp::SmoothImage(hoImage0, &hoImageT, "deriche1", w[i]);
但是說明一點，smoothimage的deriche1引數耗時和sigma值無關，
最近關于高斯模糊方面的我寫了不少文章，都綜合在我的SSE優化DEMO里，最近我也把這個DEMO做了更好的分類管理，如下圖所示：
有興趣的朋友可以從：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar?t=1660121429 下載，
如果想時刻關注本人的最新文章，也可關注公眾號：
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/504723.html
標籤：其他
下一篇：leetcode 114. Flatten Binary Tree to Linked List 二叉樹展開為鏈表(簡單)