2008 年在一個 PS 討論群裡，有網友不解 Photoshop 的高斯模糊中的半徑是什麼含義，因此當時我寫了這篇文章：

　　在那篇文章中，主要講解了高斯模糊中的半徑的含義，是二維正態分佈的方差的平方根，並且給出了演算法的理論描述。現在我又打算把該演算法用 c++ 實現出來，於是有了下面的這個 DEMO。

　　起初我是按照演算法理論直接實現，即使用了二維高斯模板，結果發現處理時間很長，對一個圖片竟然能達到大約數分鐘之久。這樣肯定是不對的，所以我百度了一下，發現這個問題應該採用分別進行兩次一維高斯模糊就可以了^[1]，這樣演算法的時間複雜度的一個係數，就從 O ( σ ^2 ) 降低到了 O ( σ )。這樣演算法於是速度提高到了毫秒級。下表給出分別用二維模糊的原始方法，和兩次一維模糊累加的方法的演算法成本比較：

　　其中：σ ：方差平方根（Photoshop 中的高斯模糊半徑）；n = w * h （圖片的畫素數量）。具體時間和圖片大小和高斯半徑的大小有關，一個粗略的大概情況為，演算法（1）的耗時為分鐘級，演算法（2）的耗時為毫秒到秒級。可見演算法（2）比演算法（1）速度更快，但相比演算法（1）來說演算法（2）具有較高的空間需求。

　　注：當然上面的空間複雜度並不是絕對的，例如，可以通過對影象進行序列的切片處理，既可減小演算法（2）的空間需求。

　　兩種演算法在高斯半徑為常數條件下，都是關於圖片大小的線性演算法，區別在於常數係數的大小不同，前者是高斯半徑（模板尺寸）的平方級，後者是高斯半徑（模板尺寸）的線性級別。這個改進，非常類似於我此前有一篇部落格中給出的，對一個油畫效果濾鏡的演算法改進，也是通過把常數係數，從模板尺寸的平方級別降低到線性級別，使演算法速度獲得提高的。

　　在理論上，高斯模板是無邊界和無限擴充套件的一個二維曲面，在實現時，就必須對這個曲面截斷為有限大的二維模板。那麼在哪裡截斷呢？根據下圖所示的一維正態分佈貢獻：

　　圖1. 正態分佈的貢獻比

　　此圖來自參考資料 [1]，根據資料文中敘述，此圖實際來源於（Maybe blocked by the GFW）

　　從圖 1 中可以看到，在 3σ 以外的貢獻比例非常小，為 0.1 %，因此我們截斷模板時，對模板邊界定義為 3σ ；

　　int r = ( int ) ( sigma * 3 + 0.5 );  // 完整模板的邏輯尺寸：( 2 * r + 1 ) * ( 2 * r + 1 );

　　二維高斯模板的計算公式是：

　　下圖給出了二維模板的視覺化結果。採用的視覺化方法是，根據上面的公式和模板邊界，生成二維高斯模板，然後取一個縮放因子 f = 255 / 模板中心點的資料，以此縮放因子把模板資料等比縮放，然後繪製成灰度圖片，這樣中心點的亮度就被提高到最亮。視覺化效果中，每個單元格對應著一個模板資料，單元格大小為 8 * 8 或者 16 * 16 畫素。

　　圖 2. 二維模板的視覺化結果

　　圖 2 中，左側是人們常見的 3 x 3 模板（σ ≈  0.849），圍繞中心點的 3 x 3 的浮點資料為：

　　採用演算法（2），要完成高斯模糊，對圖片分別進行兩個方向的一維高斯模糊即可。例如，先對圖片進行水平方向的模糊，得到中間結果，然後再對這個中間結果進行垂直方向的模糊，即得到最終結果。下圖是一個演示圖，給出了原圖在兩個方向上分別單獨進行一維高斯模糊的結果，以及最終的結果：

　　圖3. 演算法（2）中一維模糊的中間結果

　　僅在這個圖片的例子中，我把我寫的演算法的處理結果，在 Photoshop 中開啟和 Photoshop 自帶的高斯模糊的處理結果做差值對比，發現兩者是相同的。

　　我實現的 DEMO 程式（Windows 平臺）的介面如下所示：

　　圖4. DEMO 程式的主視窗 UI

　　通過點選選單 - 視覺化 - 二維高斯模板可以在右側的檢視中生成一個灰度圖片，即二維高斯模板的視覺化結果。

　　在程式介面的客戶區下方有一個控制面板，可以選擇高斯模糊的演算法引數，高斯半徑的意義和 Photoshop 中的高斯模糊半徑的意義相同，都是演算法中的 σ。

　　演算法引數中：

　　（1）支援多執行緒處理。根據我的觀察，執行緒數設定為和 CPU 核心數相同是比較合適的。執行緒數比 CPU 核心數更多，也是沒有什麼意義的，因為演算法執行時，CPU 已經滿負荷運轉了。開啟更多執行緒，也不能再提高速度了。

　　假設 CPU 核數為 p，開啟的多執行緒數量 >= p，則演算法速度大約為單執行緒處理的 p 倍。（當 CPU 滿負荷時，執行緒數量取得更大，也沒有提高速度的意義了）

　　注：此處的 CPU 核心數應該為 CPU 的物理核心數，而非模擬出來的多核數目。

　　（2）浮點型別：支援 float 和 double。它是高斯模板的資料的型別，也是進行畫素加權累加時的資料型別，根據我的觀察，float 和 double 的速度相差不大。基本相同。

　　（3）高斯半徑：即 σ。演算法的常數係數為 O(σ)。很顯然，σ 的值取得越大，演算法耗時將會越長。在 DEMO 中，其允許範圍和 Photoshop 的要求一致，是 [0.1, 250]。

　　在實現演算法時，我也嘗試了對 255 個灰度值 * 模板資料的結果進行快取和查表處理，但是發現不能有效提高速度，所以最終我放棄了這種方法。這可能是因為，演算法的計算只是一個浮點乘法，對資料的讀取動作，並不能做到比浮點乘法更快。所以這裡採用快取也就顯得沒有必要了。

　　在本 DEMO 中，濾鏡處理是放在 UI 執行緒中進行的，這使得在濾鏡處理時間較長時（例如高斯半徑取值很大，圖片也很大），介面會有些卡，可以把濾鏡處理動作放在一個新建的後臺執行緒中執行。這是比較容易實現的。

　　對演算法的使用方法：

　　在 C++ 程式中，使用我寫的這個演算法是非常簡單的，例如：

#include "GaussBlurFilter.hpp"

CGaussBlurFilter<double> _filter;
_filter.SetSigma(3.5); // 設定高斯半徑
_filter.SetMultiThreads(true, 4); // 開啟多執行緒，使用者建議的執行緒數為 4；

// lpSrcBits / lpDestBits: 畫素資料的起始地址，必須以 4 bytes 對齊，
// 注意：不論高度為正或者負，lpBits 都必須為所有畫素中地址值最低的那個畫素的地址。
// bmWidth, bmHeight: 影象寬度和高度（畫素），高度允許為負值；
// bpp: 位深度，支援 8（灰度）, 24（真彩色）， 32
_filter.Filter(lpSrcBits, lpDestBits, bmWidth, bmHeight, bpp);

　　需要注意的是，在多執行緒處理中，我使用了 Windows API （例如 CreateThread）等，這使得 GaussBlurFilter.hpp 目前只能用在 Windows 平臺，如果要在其他平臺使用，應當修改和多執行緒有關的 API 函式呼叫。

　　高度值可以為正也可以為負，但畫素資料的地址 lpBits 都必須是所有畫素中，地址值最小的那個畫素的地址。即，假設圖片左上角點的座標為原點，如果圖片高度為正數（bottom - up），則 lpBits 是左下角畫素 （col = 0，row = height - 1）的地址。如果圖片高度為負數（top-down），則 lpBits 是左上角畫素（col = 0，row = 0） 的地址。影象資料的掃描行寬度必須以 4 Bytes 對齊，即通過下面的公式計算掃描行寬度：

　　int stride = ( bmWidth * bpp + 31 ) / 32 * 4; //掃描行寬度，對齊到 4 Bytes

（上式為程式語言表達，非數學表達，即利用了整數除法對小數部分的截斷性。）

　　bpp：影象的畫素位深度。只支援 8 （灰度索引影象），24，32 這幾個值。對於 32 bpp 的影象來說，最後一個畫素通道是表徵畫素的不透明度，也就是 alpha，對於 alpha 如何參與到演算法中，我想了下，有多種處理方法，但都好像沒有什麼容易理解的物理意義，所以在程式碼裡我忽略了 alpha 通道。

　　【相關下載】：

　　（1）Demo 可執行檔案（包含 GaussBlurFilter.hpp）：GaussBlurDemo_Bin.zip

　　（2）Demo 完整原始碼（包含 GaussBlurFilter.hpp 和 可執行檔案）：GaussBlurDemo_Src.zip

　　【參考資料】