[訊號處理] 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 何者運算較快？影像太大該怎麼壓縮？等等。

[MATLAB] 如何使用 Latex 數學符號來標示圖形的 x,y 軸

一般在 MATLAB 使用圖形常會加入 x,y 軸來幫助讀者了解圖形內容，有時候我們想要在 x-label 顯示 數學式子來使其更為簡潔：比如說我想要在 x-軸加入 $\hat{K}$ 則只要在 MATLAB 中使用以下指令即可：

xlabel('$$\hat{K}$$','Interpreter','Latex')

注意到上述指令中， "＄＄...＄＄" 符號用來告知 MATLAB 我們要使用 Latex 語法，下圖為執行上述指令 x-label 後，在圖形上顯示的樣子：

更為進階的 latex語法也完全支援：

xlabel('$$\hat{K}, \tilde{x}, \int_\mathcal{X} f_X(x)dx$$','Interpreter','Latex')

上述指令執行後結果如下圖：

[數值方法] Forward Euler's Method 與 其改良型 求解 ODE

給定 初始值問題(Initial Value Problem, IVP)
\[
y'(t) = f(t,y(t)),\;\;\;\; y(0) = y_0
\] 一般而言，上述 一般的 IVP 問題並沒有辦法寫下解析解，但我們可以退而求其次詢問是否有合適的數值方法求解上述 初始值問題，最為基本的想法是利用所謂的 Forward Euler method:

Forward Euler Method:
\[
y_{n+1} = y_n + h \cdot f(t_n, y_n)
\]其中 $t_n = n h$。令終止時間 $T$ 在 $n$ 步之後到達，則 $T = n h$。

Comment:
1. 上式中 $h$ 稱為 迭代步長 (step size)。
2. $t_{n+1} = t_n + h$
3. 儘管 Forward Euler Method 非常簡便，但其數值誤差相當大，以下我們看個實際例子：


Example: 
考慮 IVP $ y'=t-y $ 且 $ y(0) = 0 $ ，
(a) 試求解上述 IVP
(b) 現在令 $h := 0.1$ ，試求用 Forward Euler Method 求解 $y_2$
(Note: $y_2 = y(t_2) = y(2h)$)
(c) 比較 (a) 與 (b)

Solution (a):
首先改寫原式：
\[y'\left( t \right) =  - y\left( t \right) + t
\]由上式可知此為線性常係數 ODE 其解可立即求得
\[\begin{array}{l}
y\left( t \right) = {e^{ - 1}}y\left( 0 \right) + \int_0^t {{e^{ - 1\left( {t - s} \right)}}sds} \\
 \Rightarrow y\left( t \right) = 0 + \underbrace {{e^{ - t}}\int_0^t {{e^s}sds} }_{ = {e^{ - t}}\left( {{e^t}(t - 1) + 1} \right)}\\
 \Rightarrow y\left( t \right) = (t - 1) + {e^{ - t}}
\end{array}\]
Solution (b):
由 $y_0 := y(0)=0$ 出發，則我們可以計算
\[\begin{array}{*{20}{l}}
{{y_1} = {y_0} + h \cdot f({t_0},{y_0})}\\
{ \Rightarrow {y_1} = 0 + \left( {0.1} \right) \cdot \left( {0 \cdot h - {y_0}} \right)}\\
{ \Rightarrow {y_1} = \left( {0.1} \right) \cdot \left( {0 - 0} \right) = 0}
\end{array}\]同理，我們接著計算 $y_2$
\[\begin{array}{l}
{y_2} = {y_1} + h \cdot f({t_1},{y_1})\\
 \Rightarrow {y_2} = {y_1} + h \cdot \left( {1 \cdot h - {y_1}} \right)\\
 \Rightarrow {y_2} = 0 + \left( {0.1} \right) \cdot \left( {0.1 - 0} \right)\\
 \Rightarrow {y_2} = 0.01
\end{array}
\]
Solution (c)
由 (a) 可知
\[\begin{array}{l}
y\left( {2h} \right) = (2h - 1) + {e^{ - 2h}}\\
 \Rightarrow y\left( {0.2} \right) =(0.2 - 1) + {e^{ - 0.2}} \approx {\rm{0.01873}}
\end{array}
\]但由 (b) 可知我們得到的 $y_2 = y(2h) =0.01$ 亦即具有誤差
\[
|y(2h) - y_{2}| = |0.01873 - 0.01| = 0.0873
\]

以下為簡單的 MATLAB code 實現 Forward Euler Method:

t0 = 0; %Initial time

y0 = 0; %Initial value

h = 0.1; %Step size

tfinal = 1; %Final time

F = @(t,y) t - y;

y = y0;

yout = y;

for t = 0: h: tfinal-h

    y = y + h * F(t,y); %Euler Method

    yout =[yout; y];

end

t_interval = 0: h: tfinal

plot(t_interval,  yout, 'o-')

Improved Euler Method
為了改良上述誤差過大的問題，我們可採用以下 Predictor-Corrector 架構來改善 Forward Euler Method，亦即我們首先透過 Forward Euler Method 建構

Predictor Equation
\[
\tilde{y}_{n+1} := y_n + h \cdot f(t_n,y_n)
\]接著再修正上述的預測結果:
Corrector  Equation:
\[
y_{n+1} = y_n + \frac{h}{2}(f(t_n,y_n) + f(t_{n+1}, \tilde{y}_{n+1}))
\]
我們以下再次使用前述例子來看看改良型 Euler 是否有得到更加精準的結果：

Example (Revisit): 
考慮 IVP $y'=t-y$ 且 $y(0)=0$ ，
(a) 令 $h := 0.1$ ，試求用 Improved Euler Method 求解 $y_2$

Solution:
由 $y_0 = 0$ 開始，為了計算 $y_1$ 我們先透過 predictor 計算 $\tilde{y}_{1} $:
\[
\tilde{y}_{n+1} := y_n + h \cdot f(t_n,y_n) = 0
\]接著用Corrector  Equation:
\[\begin{array}{l}
{y_1} = {y_0} + \frac{h}{2}\left( {f({t_0},{y_0}) + f({t_1},{{\tilde y}_1})} \right)\\
 \Rightarrow {y_1} = 0 + \frac{h}{2}\left( {\left( {{t_0} - {y_0}} \right) + \left( {{t_1} - {{\tilde y}_1}} \right)} \right)\\
 \Rightarrow {y_1} = \frac{h}{2}\left( {\left( {0 \cdot h - 0} \right) + \left( {1 \cdot h - 0} \right)} \right)\\
 \Rightarrow {y_1} = \frac{h}{2}h = \frac{1}{2}{\left( {0.1} \right)^2} = 0.005
\end{array}
\]有了 $y_1$ 我們再透過 Predictor 計算 $\tilde{y}_{2} $:
\[\begin{array}{l}
{{\tilde y}_2} = {y_1} + h \cdot f({t_1},{y_1})\\
 \Rightarrow {{\tilde y}_2} = 0.005 + h \cdot \left( {{t_1} - {y_1}} \right)\\
 \Rightarrow {{\tilde y}_2} = 0.005 + h \cdot \left( {1 \cdot h - 0.005} \right)\\
 \Rightarrow {{\tilde y}_2} \approx {\rm{0.01450}}
\end{array}\]再用 Corrector 計算 $y_2$
\[\begin{array}{l}
{y_2} = {y_1} + \frac{h}{2}(f({t_1},{y_1}) + f({t_2},{{\tilde y}_2}))\\
 \Rightarrow {y_2} = 0.005 + \frac{h}{2}(\left( {{t_1} - {y_1}} \right) + \left( {{t_2} - {{\tilde y}_2}} \right))\\
 \Rightarrow {y_2} = 0.005 + \frac{{0.1}}{2}(\left( {1 \cdot h - 0.005} \right) + \left( {2 \cdot h - {\rm{0.01450}}} \right))\\
 \Rightarrow {{y}_2} \approx {\rm{0.0190}}
\end{array}
\]現在回憶前述例子中，我們知道 $0.01873$ 此表明此例中，improved Euler 確實大幅降低計算誤差。

Comments:
1. 實際數值應用中，大多數數值方法採用所謂 Runge-Kutta Method 來獲取更加精確或者穩定的數值解，但未免離題，在此不做贅述。

[MATLAB] 如何對 ezplot 擷取圖形資料

MATLAB 中所提供的 symbolic toolbox 是非常強大的符號運算工具，其中的 ezplot 幫助我們能對未定符號進行快速畫圖，現在假定我們非常滿意圖形的結果，並且想試圖擷取 ezplot中的 x軸與 y軸的 圖形資料該怎麼辦？

以下我們用一個例子說明如何擷取圖形資料:

syms x %使用 symbolic toolbox 定義變數 $x$ 

f = ezplot( x^2, [-1,1] ); %使用 ezplot 繪製 $x^2$ 且限定 區間 $[-1,1]$

f_xaxis = get(f, 'xdata' ); %從 ezplot 提取 圖形 $x$ 軸資料

f_yaxis = get(f,  'ydata');  %從 ezplot 提取 圖形 $y$ 軸資料

那麼上述中的 f_xaxis 即為我們原函數在ezplot 的x軸資料， f_yaxis 為我們原函數在ezplot 的y軸資料。一旦擷取完畢，要做剩餘的運算會變得相當容易，比如說要找ezplot 圖中的最大值，我們可以直接針對剛剛截取到的 y軸資料使用

max(f_yaxis)

NOTE:
由於我們是從 ezplot 的 x軸與 y軸資料來擷取最大值，故此法不保證取得的結果為 "真正" 的極值。僅是一種簡便的方法。如果讀者目的在於要找極值，則建議使用 fmincon 或者 fminunc 等函數會更為精準。

[MATLAB] 將 symbolic expression 轉成 latex 程式碼

一般而言在 MATLAB 使用中不免會碰到使用 symbolic toolbox 情況，但有時表示式非常繁雜，如果又想要把該表示方程式轉寫成 latex 貼到論文中該怎麼辦？

MATLAB 提供一個非常方便的功能

latex(.)

可以幫助我們直接轉換 MATLAB symbolic expression 成為 LATEX code

[隨機分析] Euler-Maruyama 法求 隨機微分方程 數值解 (利用 MATLAB)

在此我們介紹 Euler-Maruyama Method 求解 隨機微分方程   (Stochastic Differential Equation, SDE)，現在 考慮  SDE 的積分型式 可寫為：
\[
X(t) = X_0 + \int_0^t f(X(s))ds + \int_0^t g(X(s))dW(s), \;\; 0 \le t \le T
\]其中 $f, g$ 為 純量函數 且 初始值 $X_0$ 為隨機變數；另外 第二項積分為 Ito integral。上式若有解則其解 $X(t)$ 對任意 $t$ 皆為隨機變數。

一般而言，上式可改寫為較為簡潔的 隨機微分方程的微分形式：
\[
dX(t) = f(X(t))dt + g(X(t)) dW(t),\;\; X(0)=X_0, \; 0 \le t \le T \ \ \ \ \ (*)
\]
Comments: 
1. 讀者需注意我們不可寫 $dW(t)/dt$ 因為 Browian motion 為 處處連續但處處不可微 (with probability 1)
2. 若 $g = 0$ 且 $X_0$ 為常數，則 SDE 退化成一般的 ODE；亦即
\[
dX(t) = f(X(t))dt, \;\; X(0)=X_0, \; 0 \le t \le T
\] (此時即可用 Euler method 求數值解。)
3. 關於 SDE 何時有解，讀者可參考BLOG 相關系列文章 ：
[隨機分析] Uniqueness and Existence theorem for S.D.E. (1)- Uniqueness
[隨機分析] Uniqueness and Existence theorem for S.D.E. (2) - Picard Iteration for SDE
[隨機分析] Uniqueness and Existence theorem for S.D.E. (3) - An intermediate result (Upper bound) of Picard Iteration
[隨機分析] Uniqueness and Existence theorem for S.D.E. (4) - The Existence of Solution



那麼現在回歸主題，在此我們的目的是：
對 SDE 在關心的區間 $t \in [0,T]$ 之間進行數值求解 (透過 MATLAB )。

第一步首先將此區間 $[0,T]$ 離散化，亦即 對某些 $N$ 而言，定義 $\Delta t := T/N$ 且 對 $j=1,...,N$ ，定義 $\tau_j := j \Delta t$。 並且稱 數值近似解 $X_j :=X(\tau_j) $。

利用 Euler-Maruyama (EM) 法，我們可將前述 SDE 在離散化區間內 改寫如下\[{X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,N
\]
Comment: 注意到上式若寫成積分形式即為
\[X({\tau _j}) = X({\tau _{j - 1}}) + \int_{{\tau _{j - 1}}}^{{\tau _j}} f (X(s))ds + \int_{{\tau _{j - 1}}}^{{\tau _j}} g (X(s))dW(s)
\]

Example
利用 Euler-Maruyama 法求解下列線性 SDE
\[
dX(t) = \mu X(t) dt + \sigma X(t) dW(t),\;\; X(0) = X_0
\]其中 $\mu, \sigma $ 為 常數 (亦即在此例我們選 $f(X(t)) = \mu X(t)$ 且 $g(X(t)) = \sigma X(t)$)。此例為 Geometric Brownian Motion (GBM) 模型，且此 SDE 有解析解如下：
\[X\left( t \right) = {X_0}{e^{\left( {\mu  - \frac{{{\sigma ^2}}}{2}} \right)t + \sigma W\left( t \right)}}\]以下我們將利用 MATLAB 實現 上述 解析解 與透過 EM 法 求數值解並比較其差異。


現在我們利用 MATLAB 實現 Euler-Maruyama 法：

首先計算 discretized Brownian path over $[0,1]$
定義 stepsize $\Delta t = R \delta t$ (一般泛取 $R=1$)且 由於 EM 法\[{X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,L\]上式還需要計算 $W(\tau_j) - W(\tau_{j-1})$，故我們計算：
\[\begin{array}{l}
W({\tau _j}) - W({\tau _{j - 1}}) = W\left( {j\Delta t} \right) - W\left( {\left( {j - 1} \right)\Delta t} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} = W\left( {jR\delta t} \right) - W\left( {\left( {j - 1} \right)R\delta t} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} = \sum\limits_{k = jR - R + 1}^{jR} {d{W_k}}
\end{array}\]上式計算出現在下列程式碼 (第14行)


我們將 MATLAB 程式碼附上如下： (考慮 $\mu=2$, $\sigma =1$ 且 $X_0=1$。 $R=1$)
(line 10 : line 14: 產生 standard Brownian Motion)
(line 16: 解析解 for GBM SDE)
(line 19: line 25: EM-method)

執行結果如下：

上圖中 紅線為 EM 法結果，藍線為 解析解結果。讀者可調整 $R$ 值來檢驗兩者差距。


ref: Desmond J. Higham, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM REVIEW Vol 43, No. 3, pp. 525-546, 2001.

[MATLAB] cdfplot 圖形變色的小技巧

這次要介紹一個 繪製機率常用 的指令：

cdfplot

此指令用來繪製 empirical cumulative probability distribution function, cdf:

用法如下：

cdfplot(X)


其中 X 為產生的隨機變數

下圖即為使用cdfplot所產生的圖形的一個例子：

另外要介紹一個小技巧：
如果想要呈現兩張以上的 cdf 圖形，則使用不同顏色來區別會是一個好方法，但是 cdfplot指令並不支援使用者改變顏色，故我們可以使用 set 指令來幫助我們 (亦可使用 plot tool 來改變顏色)：

假設有兩組隨機變數 X, Y, 那麼使用不同顏色的 cdfplot 方法如下:

    cdfplot(Ｘ) %繪製藍色線條 X 的cdf

    hold on
 
    Y_color = cdfplot(Y) %繪製紅色線條 Y 的 cdf 並記做 Y_color

    set(Y_color ,  'Color' ,  'red'); %設定顏色

    legend('X', 'Y') %對應資料名稱

    grid on  %開啟網格

那麼執行後圖形如下：(下圖的 X 與 Y 為透過 randn.m 指令產生 1,000 組 隨機變數 繪製而成)