[訊號處理] 2d Convolution & 簡單的影像處理 (利用 MATLAB )

一般在處理影像的時候我們會把 影像(image) 視為 二維離散訊號， 更一般的說法是將其視為矩陣，亦即給定任一張 (灰階) 影像 我們可以將其分割成 $M \times N$ 矩陣 (像素 pixel)，並且將其記作 $x[m,n]$，其中 $m=0,...,M-1$ 與 $n = 0,...,N-1$ 。

Comments:
如果要處理的影像是彩色的，一個常用的做法是把該影像以 $R,G,B$ 三色分別存成 三個矩陣，最後在做疊合。在此不做贅述。

以下我們利用 MATLAB 來執行 (灰階)影像處理，我們使用的圖檔是 MATLAB內建的圖檔 cameraman.tif ，當然讀者可自行讀入任何自己想要的圖檔。在MATLAB 輸入

MATLAB Code For Loading the Image
x = imread('cameraman.tif' )
imagesc(x)

則會顯示

圖1a: cameraman.tif 原圖 ($256 \times 256$)

Comments:
1. 上述影像訊號 $x[m,n]$ 為 $M \times N = 256 \times 256 $ 矩陣。
2. 一般而言在影像處理領域，我們將圖片最左上角點視為座標原點，對應的影像訊號為 $x[0,0]$ 。
3. 上述影像透過 MATLAB 2017a 讀圖會顯示有一些藍綠等顏色，若要使此圖顯示為灰階(gray scale) 讀者可鍵入 colormap(gray) 來使其變成灰階影像，亦即

MATLAB Code For Gray Scale Colormap
x = imread('cameraman.tif' )
imagesc(x)
colormap(gray)

則我們會得到如下圖

圖1b: 灰階影像

現在我們可對上述影像進行一些常見的基本處理。令 $x[m,n]$ 為輸入影像訊號，且假設 影像 濾波器 為 線性非時變 (Linear Time Invariant, LTI) 故其對應的 impulse response $h[m,n]$ 可以用來完全描述我們的影像濾波器，最後 我們令 影像處理過後的輸出訊號 $y[m,n]$ 可表為 $x[m,n]$ 與 $h[m,n]$ 的 convolution ，差別僅在此時我們的 convolution運算為二維運算，故我們寫成
\[
y[m,n] = h[m,n]**x[m,n]: = \sum\limits_{k = 0}^{M - 1} {\sum\limits_{l = 0}^{N - 1} {h\left[ {k,l} \right]x\left[ {m - k,n - l} \right]} }
\] 其中 $**$ 表示二維 convolution，在 MATLAB 中我們可以使用 conv2 來執行此運算。下圖顯示了上述討論的觀點：

Comments:
1. 上述討論中提及的 脈衝響應 $h[m,n]$ 在影像處理中又被稱為 點擴散函數 (point spread function, PSF)
2. 給定任意 影像 $x[m,n]$ 其影像尺寸為 $M_1 \times N_1$ 且 影像濾波器  $h[m,n]$ 具有 影像尺寸為$M_2 \times N_2$ 則其經過 2d convolution之後的 輸出 $y[m,n]$ 之影像尺寸會變成
$$
(M_1 + M_2 -1) \times (N_1 + N_2 -1)
$$ 故上述 2d convolution 執行之後原圖尺寸會變成
$$
(256+2-1) \times (256 + 2 -1) = 257 \times 257
$$也就是說透過2d convolution之後 影像邊緣 會多跑出一些額外的像素。讀者可參考下方後續的討論。
3. 上述討論中我們假設 影像濾波器為 LTI ，故輸入與輸出關係為  (time-domain) convolution，那麼熟悉訊號與系統的讀者不難做出以下猜想，對輸入 $x$ 與 輸出 $y$ 與 影像濾波器 $h$ 取 Discrete Fourier Transform (DFT) 可得到 對應的 DFT係數 如 \[\begin{gathered}
  X[k,l] = DFT(x[m,n]); \hfill \\
  H[k,l] = DFT(h[m,n]); \hfill \\
  Y[k,l] = DFT(y[m,n]). \hfill \\
\end{gathered} \] 假設 $N,M$ 夠大，則其在 輸入輸出在頻域上對應的關係為 frequency domain  為相乘
$$
Y[k,l] = H[k,l] X[k,l]
$$
一但計算出 $Y[k,l]$ 在對其取 Inverse Digital Fourier Transform (IDFT) 即可得到 $y[m,n]$。在 MATLAB 中，上述討論可以使用 fft2 與 ifft2 來實現，在此不做贅述。


以下我們探討幾種常見的影像濾波器，亦即幾種常見的 $h[.]$：

Image Filtering:
以下為幾種常見的影像濾波器的例子：
1. 影像模糊 (Image blurring)
\[{h_b}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  {1/10}&{1/10}
\end{array}} \right];\]

MATLAB code for Image Blurring
h_b = 1/10 * [1 1; 1 1];
y_b = conv2(x, h_b, 'same')
imagesc( abs(y_b) )

圖2: 模糊效果

讀者可以比較此圖與之前的圖不難發現此圖較原圖為模糊。另外可注意到模糊之後的影像邊緣部分出現額外不屬於原圖的像素。


2. 邊緣偵測：圖像的邊緣可以透過特定的濾波器設計來偵測，比如說我們可以考慮
\[\begin{gathered}
  {h_h}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  { - 1/10}&{ - 1/10}
\end{array}} \right]; \hfill \\
  {h_v}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{ - 1/10} \\
  {1/10}&{ - 1/10}
\end{array}} \right]. \hfill \\
\end{gathered} \]其中 $h_h$ 用以偵測 影像的水平邊緣，$h_v$ 用以偵測影像的垂直邊緣。

MATLAB Code for Horizontal Edge Detection
h_h = 1/10 * [1 1; -1 -1];
y_h = conv2(x, h_h, 'same')
imagesc( abs(y_h) )

我們得到對於影像水平邊緣的偵測如下圖所示

圖3: 水平邊緣偵測

MATLAB Code for Vertical Edge Detection
h_v = 1/10 * [1 -1; 1 -1];
y_v = conv2(x, h_v, 'same')
imagesc( abs(y_v) )

上述code得到對於垂直邊緣的偵測如下圖

圖4: 垂直邊緣偵測

Comments:
1. 讀者可自行改變 $h_v, h_h, h_b$  的矩陣的係數來看看得到的圖有什麼不同。
2. 影像處理有非常多的有趣的問題可以進一步討論，比如說被 模糊後之後的影像是否可以把它還原？如果可以該怎麼做？ 直接用 conv2 或者用 fft2 何者運算較快？影像太大該怎麼壓縮？等等。