高斯模糊演算法的實現和優化
前兩年我發過一文：Win32下的C++高斯模糊演算法例項，裡面給出了一個高斯模糊的實現，並寫了粗略的簡介。
不過當時內容講得非常簡單，而且附帶的例子演算法是有缺陷的：
- 一是對圖片的邊角採用“跳過”的方式處理，導致模糊後的圖片有黑邊；
- 二是演算法本身採用的是二維矩陣，效率上不如一維高斯模糊好。
鑑於此，這裡重新整理並試圖完整的講述一下高斯模糊。
一、高斯模糊是什麼
模糊演算法，不論是使用哪種演算法，目的都是為了讓圖片看起來不如原來那麼清晰。
清晰的圖片，畫素間的過渡會較為乾脆利落，簡而言之，就是畫素之間的差距比較大。
而模糊的本質，其實就是使用某種演算法把影象畫素和畫素之間的差距縮小，讓中間點和周圍點變得差不多；即，讓中間點取一個範圍內的平均值。
模糊到了極致，比如用於計算模糊的取值區域為整張圖片，就會得到一張全圖所有畫素顏色都差不多的圖片：
左邊是原圖，右邊是徹底模糊之後的對比圖
計算模糊的取值區域就是濾鏡區域，那麼關鍵就是，我們採用什麼曲線函式來生成平均值了。
假設濾鏡區域為正方形，且半徑為r，我們可以用如下曲線函式來計算圖片上某一個點的值：
這是一個非常簡單的直線函式，求得的點值即為算術平均值，圖片上某個點的值僅和模糊半徑r有關，與座標的位置無關。
由此我們可以得到用於模糊影象的濾鏡演算法：
void Average(filter_t& kernel, long radius) { long diamet = Diamet(radius); // (r * 2) + 1 kernel.set(radius, diamet * diamet); // kernel size is d * d double average = 1.0 / (double)kernel.size(); for(long n = 0, i = 0; i < diamet; ++i) for(long j = 0; j < diamet; ++j, ++n) kernel[n] = average; }
然後使用如下演算法遍歷整張圖片，並使用濾鏡處理畫素：
void Blur2D(bitmap_t& bitmap, filter_t& kernel) { for(long inx = 0, y = 0; y < bitmap.h(); ++y) { for(long x = 0; x < bitmap.w(); ++x, ++inx) { double r = 0.0, g = 0.0, b = 0.0; for (long n = 0, j = -kernel.radius(); j <= kernel.radius(); ++j) { long j_k = Edge(j, y, bitmap.h()); for (long i = -kernel.radius(); i <= kernel.radius(); ++i, ++n) { long i_k = Edge(i, x, bitmap.w()); long inx_k = inx + j_k * bitmap.w() + i_k; r += bitmap[inx_k].r * kernel[n]; g += bitmap[inx_k].g * kernel[n]; b += bitmap[inx_k].b * kernel[n]; } } bitmap[inx].r = Clamp<bitmap_t::channel_t>(r); bitmap[inx].g = Clamp<bitmap_t::channel_t>(g); bitmap[inx].b = Clamp<bitmap_t::channel_t>(b); } } }
上面的演算法，第7-19行，把濾鏡範圍內的畫素收集起來，並按照濾鏡給出的值將所有畫素合併計算為中心點畫素的值。
其中，Edge為邊緣處理演算法：
template <typename T>
T Edge(T i, T x, T w)
{
T i_k = x + i;
if (i_k < 0) i_k = -x;
else if (i_k >= w) i_k = w - 1 - x;
else i_k = i;
return i_k;
}它將超出邊界的濾鏡範圍檔回邊界之內；
Clamp為畫素通道的限制演算法：
template <typename T1, typename T2>
T1 Clamp(T2 n) { return (T1)(n > (T1)~0 ? (T1)~0 : n); }每個畫素通道的最大值不能超過255，若賦予的值超出了通道最大值，則將通道值限制在最大值上。
使用上面的演算法，設定濾鏡半徑為5，處理出來的圖片效果如下：
左邊是原圖，右邊是對比圖
模糊成功了，但效果有些不盡人意，圖片的一些邊緣細節模糊得並不柔和，感覺像沒戴眼鏡的近視眼看起來的效果。
從上面的演算法裡我們可以看出，濾鏡中填充的每一個值，可以看做是一個權重。在處理影象的時候，通過它計算出濾鏡範圍內每個畫素的權重值，最後將它們相加，得到的就是濾鏡中心點的畫素值。
使用上面的模糊演算法，在處理圖片的時候，對所有的畫素點是一視同仁的。但實際上，離中心點越遠的點，重要性應該越低才更合理。
也就是說，越靠近濾鏡邊緣的畫素，權重更小才會更符合實際。
高斯分佈，即正態分佈曲線，形狀大概如下圖：
它就是一個可以計算出符合上面要求的權重分佈的函式，對應的二維形式如下：
使用高斯分佈曲線作為濾鏡演算法的模糊演算法，稱之為高斯模糊。
二、高斯模糊演算法實現
在前文的闡述裡，我們已經有了高斯模糊演算法的曲線函數了。通過這個二維的曲線函式，我們可以得到與之對應的程式演算法如下：
void Gauss(filter_t& kernel, long radius)
{
long diamet = Diamet(radius); // (r * 2) + 1
kernel.set(radius, diamet * diamet); // kernel size is d * d
double sigma = (double)radius / 3.0;
double sigma2 = 2.0 * sigma * sigma;
double sigmap = sigma2 * PI;
for(long n = 0, i = -kernel.radius(); i <= kernel.radius(); ++i)
{
long i2 = i * i;
for(long j = -kernel.radius(); j <= kernel.radius(); ++j, ++n)
kernel[n] = exp(-(double)(i2 + j * j) / sigma2) / sigmap;
}
}演算法的第6-8行，先行計算出了函式裡一些與座標無關的值，然後在第14行，使用高斯二維函式計算出濾鏡裡每個點的權重。
其中第6行σ的計算，參考正態分佈曲線圖，可以知道 3σ 距離以外的點，權重已經微不足道了。反推即可知道當模糊半徑為r時，取σ為 r/3 是一個比較合適的取值。
設定濾鏡半徑為5，通過與前文一樣的遍歷演算法，我們可以得到下面的模糊效果：
左邊是原圖，右邊是對比圖
這個效果不夠明顯，我們可以設更大一點的濾鏡，比如10，處理出來的效果如下：
左邊是原圖，右邊是對比圖
好了，模糊成功了。效率怎麼樣呢？
在我的電腦上（i3M330，2.13G），使用 半徑10 的濾鏡處理 400 × 649 的圖片，Debug版本需要大約 8s 左右，Release版本則為916ms。
這個速度就算在Release下也偏慢了。
研究一下我們的演算法：
當濾鏡半徑為r時，演算法的時間複雜度是O(x × y × (2r)²)。
能否降低一點呢？
三、高斯模糊演算法優化
引用維基百科裡的原文：
“除了圓形對稱之外，高斯模糊也可以在二維影象上對兩個獨立的一維空間分別進行計算，這叫作線性可分。這也就是說，使用二維矩陣變換得到的效果也可以通過在水平方向進行一維高斯矩陣變換加上豎直方向的一維高斯矩陣變換得到。”
這個特徵讓我們可以使用一維的高斯函式，在橫縱兩個方向上分別處理一次影象，得到和二維高斯函式一樣的效果。
使用一維高斯函式，演算法的時間複雜度就變為O(2 × x × y × 2r)，兩相比較，演算法效率高出r倍。
一維形式的高斯函式如下：
從而得到新的程式演算法如下：
void Gauss(filter_t& kernel, long radius)
{
kernel.set(radius, Diamet(radius));
static const double SQRT2PI = sqrt(2.0 * PI);
double sigma = (double)radius / 3.0;
double sigma2 = 2.0 * sigma * sigma;
double sigmap = sigma * SQRT2PI;
for(long n = 0, i = -kernel.radius(); i <= kernel.radius(); ++i, ++n)
kernel[n] = exp(-(double)(i * i) / sigma2) / sigmap;
}同樣，演算法的第7-9行，先行計算出座標無關的值，然後在第12行中使用高斯函式計算出權重。
為了配合一維濾鏡，遍歷演算法需改為如下形式：
void Blur1D(bitmap_t& bitmap, filter_t& kernel)
{
Normalization(kernel);
buffer_t buff(bitmap);
for(long inx = 0, y = 0; y < bitmap.h(); ++y)
{
for(long x = 0; x < bitmap.w(); ++x, ++inx)
{
for(long n = 0, i = -kernel.radius(); i <= kernel.radius(); ++i, ++n)
{
long i_k = Edge(i, x, bitmap.w());
long inx_k = inx + i_k;
buff[inx].r += bitmap[inx_k].r * kernel[n];
buff[inx].g += bitmap[inx_k].g * kernel[n];
buff[inx].b += bitmap[inx_k].b * kernel[n];
}
}
}
for(long inx = 0, x = 0; x < bitmap.w(); ++x)
{
for(long y = 0; y < bitmap.h(); ++y)
{
inx = y * bitmap.w() + x;
double r = 0.0, g = 0.0, b = 0.0;
for(long n = 0, i = -kernel.radius(); i <= kernel.radius(); ++i, ++n)
{
long i_k = Edge(i, y, bitmap.h());
long inx_k = inx + i_k * bitmap.w();
r += buff[inx_k].r * kernel[n];
g += buff[inx_k].g * kernel[n];
b += buff[inx_k].b * kernel[n];
}
bitmap[inx].r = Clamp<bitmap_t::channel_t>(r);
bitmap[inx].g = Clamp<bitmap_t::channel_t>(g);
bitmap[inx].b = Clamp<bitmap_t::channel_t>(b);
}
}
}演算法中：
第3行，用於對濾鏡做“歸一化”處理；
第5行，對待處理的圖片生成一塊快取區，用於儲存處理的中間結果；
第15-17行，將第一次橫方向的處理結果儲存在快取區中；
第32-34行，將第二次縱方向的處理結果儲存在臨時變數中；
最後，第36-38行，將最終的處理結果賦值給圖片畫素的每個通道。
關於第3行的歸一化處理，演算法如下：
void Normalization(filter_t& kernel)
{
double sum = 0.0;
for(int n = 0; n < kernel.size(); ++n)
sum += kernel[n];
if (Equal(sum, 1.0)) return;
for(int n = 0; n < kernel.size(); ++n)
kernel[n] = kernel[n] / sum;
}目的是讓濾鏡的權重總值等於1。否則的話，使用總值大於1的濾鏡會讓影象偏亮，小於1的濾鏡會讓影象偏暗。
第5行及後面的影象處理演算法採用了快取區。
這是因為若不使用快取區快取中間的影象結果，那麼第一遍處理完的影象就只能儲存在bitmap裡了。而bitmap的畫素，每通道都是0-255的整數，如果把中間結果儲存在bitmap的記憶體空間裡，第二遍遍歷的處理精度就會大大降低，最後處理出來的效果質量也將大打折扣。
同樣使用半徑為10的濾鏡，處理後的影象效果如下：
從左往右，第一張為原圖，第二張為一維高斯模糊的效果，第三張為二維高斯模糊的效果
雖然理論上一維和二維的處理應該是沒有差別的，但是由於計算精度的原因，實際處理後的效果還是有差別的。
但是差別很細微，就上圖而言，一維處理之後的圖片，模糊效果更加平滑一些。
使用一維演算法，處理同樣圖片的時間，Debug下只需要846ms，Release下僅僅78ms，差不多都是二維演算法的10倍左右。效率的差距和前面計算的時間複雜度差距一致。
四、其它模糊演算法
從上面的文字裡，我們已經知道了高斯模糊演算法其實就是利用高斯函式或者說曲線，生成模糊濾鏡，然後對影象畫素做處理的演算法。
那麼我們可以考慮使用其它曲線來作為我們的模糊曲線麼？
當然是可以的。
文章的最開始就使用了最簡單的算術平均值曲線作為模糊曲線，只是實際的效果較為一般。
我們可以再嘗試其他的曲線函式，比如直線函式：
求絕對值是為了讓函式曲線保持“中間高兩邊低”
使用其他函式得到的濾鏡演算法和處理的圖片效果本文就不再一一贅述了。感興趣的朋友可以下載文後所附的程式碼。
附件裡的模糊演算法實現在filter.h中，VC和gcc編譯器下均編譯通過。測試平臺基於win32下的Qt 5.0.1，編譯器為MinGW 4.7。