李侑青的高等影像處理教室

影像處理與圖形識別常使用空間遮罩運算 (Mask Operation)，

或稱為旋積 (Convolution)。在基礎影像處理課程裡，

大家可能學過平均濾波器 (Average Filter)、邊偵測用的 Sobel Operator，

更複雜一點的遮罩，

如高斯平滑化 (Gaussian Smoothing)、LOG (Laplasian of Gaussian)。

這些遮罩影像處理計算，都可以用數學上的旋積計算來表達。

這篇文章要分享的是目前快速計算旋積的幾種常用方法。

假設影像大小為 N × N，遮罩大小為 M × M。

對影像中的每個點進行遮罩運算，則每個點需要 M² 次計算，

整張影像就需要 N²M² 次計算。

這種方法唯一的好處是程式容易撰寫，只要四個迴圈，容易實踐。

相信不少人在初學影像處理課程時，都寫過類似的程式。

然而這種計算方式太過費時，對於目前動輒千萬畫素的影像來說，

無法做到即時影像處理。

所幸，在基礎影像處理課程裡面，

有提供傅立葉變換 (Fourier Transform) 這項好用工具。

傅立葉變換能將影像從空間域 (Spatial Domain) 轉換至頻率域 (Frequency Domain)，

我們稍微複習一下一維的正轉換與反轉換公式：

正轉換是反轉換是  

並且，一個訊號經由正轉換與反轉換之後，還是原來的訊號，

也就是： f(t) = F^-{ F{ f(t) } }，左邊公式的意思，

就是把訊號 f(t) 丟給 F{} 轉換，

然後再把 F{} 轉換完的結果丟給 F^-{} 做反轉換，

所得的結果和原始訊號 f(t) 一樣。

二維的傅立葉變換相當於在兩個正交方向上各做一次一維的傅立葉變換，

故不冗餘述說。並且在影像處理實務上，

都是直接使用快速傅立葉變換 (Fast Fourier Transform, FFT) 計算。

更多關於傅立葉變換的介紹，

可以至 Wiki 查詢 Fourier Transform，或閱讀影像處理相關書籍。

這裡我將直接運用傅立葉變換應用於影像處理的特性：

空間域的旋積等於頻率域的相乘。

假設有一個影像 f(x, y) 以及一個遮罩 g(x,y)，

影像的遮罩計算以數學表達就是旋積計算： f(x,y)*g(x,y)，

利用傅立葉變換可以得到：  

公式就是利用空間域的旋積等於頻率域的相乘這一特性來計算旋積。

先把影像和遮罩個別轉換到頻率域，然後從低頻到高頻彼此一一交乘，

再將交乘的結果做反轉換，就可以得到整張影像的旋積計算結果。

影像與遮罩經過傅立葉變換在頻率域計算旋積，有什麼好處嗎？

下面跟大家說明。

對於一個 N² 大小的影像，快速傅立葉變換需要 N² log N 次計算。

而影像與遮罩在頻率域相乘，共有 N2 個點，因此需要 N2 次計算。

整體來說，影像的正轉換、遮罩的正轉換，頻率域相乘、反轉換，

各別需要 N² log N、N² log N、N²、N² log N 次計算。

因此整個計算的複雜度是： O( N² log N )。

相較於在空間域直接做遮罩計算的複雜度：O( N²M² )，當 M² >> log N 時，

在頻率域計算旋積的速度會遠快於直接在空間域進行遮罩計算。

目前影像處理或圖形識別相關的國際期刊論文在討論快速旋積計算時，

基本上都是使用上述的方法。唯一的差別在於快速傅立葉 (FFT) 的改進方法。

而一維快速傅立葉變換的速度一直卡在 N log N 這個時間量級，

頂多再乘上一個分數，稍微再快一點點，

目前尚未有人證明這是否已經是 FFT 的理論速度極限。

對於通用性的快速遮罩計算（或說快速旋積計算），一般都是使用上述方法。

這也表示一般的遮罩計算速度受到快速傅立葉變換的限制。

在某些影像處理的應用裡，遮罩的數值可能不是許多個，而是同一個數值。

比如說加總遮罩，遮罩內的每個數值都是 1，

或者是平均濾波器 (Average Filter) 的遮罩，相當於加總遮罩乘上某個常數。

這類數值單一的遮罩計算，可以使用一些技巧，讓計算速度更快。

以加總遮罩計算來說，我們需要求得影像某個位置周圍的加總值。

如下圖所示，藍色框是影像範圍，我們想要求紅色遮罩區域的加總值。

那麼計算該區加總值的過程，可以用下圖的過程來替代。

一、假設已有第一張圖裡，紅框區域的加總值

二、減去第二張圖裡，上方綠色區域的加總值

三、再減去第三張圖裡，左方綠色區域的加總值

四、再加上最後一張圖裡，左上角綠色交疊區域的加總值

以上過程假設我們可以很快地得到從左上角到某點的矩形區域加總值，

那麼任意區域的加總值，可以用兩個減法與一個加法就算出來。

換句話說，無論是加總遮罩計算還是平均濾波器計算，

一個 N² 大小的影像，我們都可以在 O( N² ) 時間計算完畢。

相較於頻率域旋積的時間複雜度 O( N² log N )，又快了 log N 倍。

那麼現在的疑問在於，這個假設成立嗎？

有什麼方法可以很快得到包含左上角的矩形區域加總值呢？

有的。在演算法上，有個稱為 Scan 的演算法可以使用，

也有人將這演算法稱為 Prefix Sum。

意思是運用遞迴的方式，逐項計算每一項與之前所有項的總和。

以一維來說，假設有一個陣列 f [x]，

我們想計算第 x 項到第 1 項的總和：s [x] = f [1] + ... + f [x]

對於每個位置 x，我們不必每次都從第 1 項算到第 x 項，

而是可以遞迴地計算 s [x] = f [x] + s [x-1]

相同的原理，對於二維陣列來說，

從第 [1, 1] 位置到第 [u, v] 這矩形區間的總和，

也可以遞迴地計算：s [u, v] = f [u, v] + s [u-1, v] + s [u, v-1] - s [u-1, v-1]

得到所有位置的加總表。對於 N² 大小的影像，需要 3N² 次計算。

那麼我們要計算右下角是 [u, v]、左上角是 [u-X+1, v-Y+1] 的矩形區域總和，

就可以利用前面得到的加總表快速算得：

因此每個位置的加總遮罩計算，只需要 3 次加減法就能得到。

整體來說，二維平均濾波或加總遮罩計算，只需要 O( N² ) 次計算即可。

並且相對於頻率域的旋積計算，更省去了浮點數的乘法與指數運算，

速度又更快許多。

大家現在可能會覺得，似乎 Scan 演算法就只能夠用在平均遮罩運算，

下面我就要舉一個影像處理的期刊論文，讓大家看看它的應用。

論文名稱：Fast normalized cross correlation for defect detection

來自於期刊：Pattern Recognition Letters, no.24, pp.2625-2631, 2003.

這篇論文在講的是，假設有一張標準影像 t (x, y)，

我們現在輸入一張相同大小的待測影像 f (x, y)，

它要計算每個影像點 (x, y) 的相似程度。

計算方式是把 (x, y) 周圍 M 範圍內的所有點拿來計算相關係數，M = 2m+1，

公式為：  

公式裡面，待測影像 f 在 (x, y) 的平均值為  

類似地，標準影像 t 在 (x, y) 的平均值為  

如果我們直接使用公式去計算兩張影像各點的相關係數，

則時間複雜度為 O (N²M² )。

從這邊我們可能還看不出來要如何加速計算相關係數。

因此，需要把公式做一些整理。

首先是分母的部分，它需要計算 f 在 (x, y) 局部的變異數，

我們可以把它展開來，得到  

再把  f 在 (x, y) 的平均值帶進去展開，得到  

用相同方式整理標準影像的局部變異數，可得到  

到此為止，分母的部分整理完畢。

現在我們再來整理分子的部分。首先，將公式裡的乘法展開

我們看一下最右邊一項，它是待測影像的平均值減去標準影像的某個加總。

由於在局部內，各點減去平均值的總和為零，（請讀者自行證明了，很簡單）

所以最右一項的值為零，因此公式就剩下兩個項目：  

再將標準影像的局部平均值帶入整理：    

綜合分母與分子項目，整個相關係數的計算變成

當中，每個局部加總項的計算，我們都可以用 Scan 演算法來加速。

所以，計算一個相關係數，只要 3 次減法、3 次乘法、1 次除法，

對於整張影像來說，計算時間複雜度為 O( N² )。

大家看一下，從 O( N²M² ) 加速到 O( N² )，

只要 M = 10 就相當於加速 100 倍！

所以大家可別瞧不起加總運算的加速，

這篇論文讓我們看到相關係數計算的加速，再推廣一點思考，

是不是關於統計方面的運算，都能用 Scan 演算法加速呢？

除了平均濾波器這類單值的遮罩之外，

高斯平滑化 (Gaussian Smoothing)、

高斯的拉普拉斯運算 (Laplasian of Gausian, LoG)，

這些高斯相關的遮罩在影像處理時也很常被使用到。

那麼，有沒有快速的高斯平滑演算法呢？

很可惜，目前沒有；現有的加速方式，就是使用頻率域的加速計算。

不過，既然是做平滑化、模糊化的計算，有時候精確性就不那麼重要。

若能允許一些誤差存在，那麼，

用近似演算法 (Approximation) 是加速遮罩運算的好選擇。

思考觀念與 Scan 演算法類似。

在 Scan 的遞迴裡面，當前的值是由新的點加上過去的值。

而近似方式也可以用遞迴來完成，當前的值是由新的點加上過去值的線性組合。

為了方便解說，我以一維信號來舉例。我要把信號 f (x) 和遮罩 g (x) 做旋積計算，

則計算結果 y (x) = f (x) * g (x)

每個點 x 都要做一次遮罩運算，太過耗時，因此用近似方式遞迴完成

y (x) = f (x) + a₁ y (x-1) + a₂ y (x-2) + ... + a_n y (x-n)

則長度為 N 的信號與 M 的遮罩做旋積，原本需要 O( NM ) 次計算，

現在每個點的計算次數固定，整個旋積用近似方式計算只要 O( N ) 的時間。

舉例來說，讓我們再看一篇論文。

論文名稱是：Recursive implementation of the Gaussian Filter

來自期刊：Signal Processing, no.44, pp.139-151, 1995.

該篇論文提出高斯平滑濾波器的近似遞迴公式：

y (t) = B x (t) + ( b₁ y (t-1) + b₂ y (_t-2) + b₃ y (t-3)) / b₀

其中，

B = 1 – (b₁ + b₂ + b₃) / b₀

b₀ = 1.57825 + 2.44413 q + 1.4281 q² + 0.422205 q³

b₁ = 2.44413 q + 2.85619 q² + 1.26661 q³

b₂ = -1.4281 q² – 1.26661 q³

b₃ = 0.422205 q³

然後 s 就是高斯函數公式裡面的標準差，控制著高斯濾波器的長相。

該論文用大量的實驗，尋找適合近似高斯遮罩的線性組合參數，

實驗顯示近似效果良好，雖然有一點點誤差存在，

但相較於旋積計算，用近似的方式能獲得速度上的優勢。

因此，如果某個影像處理的應用可以容許一些誤差，

大家不妨可以思考去設計線性組合來近似遮罩運算。

總上結論，在這篇文章裡首先我介紹了：

傳統的空間遮罩演算法，時間複雜度為 O ( N²M² )

然後介紹傳統的頻率域加速演算法，時間複雜度為 O ( N² log N )

接著介紹遮罩值唯一時可以用 Scan 算法加速，時間複雜度為 O ( N² )

最後介紹使用線性組合來近似遮罩運算，時間複雜度為 O ( N² )

這些是影像處理在快速旋積計算的基本概念。

李侑青,資訊工程博士

影像處理,Image Processing,旋積,Convolution,遮罩運算,Mask Operation

快速傅立葉變換,Fast Fourier Transform,高斯模糊,高斯平滑,Gaussian Smoothing