影像復原

影像復原的目的是在預先定義好的意義上改善一幅影像，不同於影像增強主要是一個主觀的程式，影像復原大致為一個客觀的程式。修復是利用退化現象的某種先驗知識，試圖把已經退化的影像加以重建或修復。

## 影像退化 / 復原程式的模型

如下圖，退化程式可以被模式化成一個退化函數（Degradation function），連同加成性雜訊（Noise）η(x,y) 共同作用在一輸入影像 f (x,y) 上，產生一退化影像 g (x,y)：

$$
g(x,y)=H[f(x,y)]+\eta (x,y)
$$

因此就可以利用對退化函數 H 以及雜訊 η(x,y) 的了解，來獲得一個原始影像的估測 f ̂(x,y)。如果 H 為一個線性空間不變數（linear spatially invariant）的程式，則可以證明退化影像在空間域（spatial domain）為

$$
g(x,y)=h(x,y)*f(x,y)+\eta (x,y)
$$

其中 h (x,y) 是退化函數的空間表示，\* 代表迴旋積（convolution），因此我們可以將上式中的模型寫成等效的頻率域（frequency domain）表示式：

$$
G(u,v)=H(u,v)F(u,v)+N(u,v)
$$

其中大寫字母的各項是對於迴旋積方程式中各項的傅立葉轉換。而影像的復原可以大致分為兩個情形進行分析，第一為假定 H 是一個恆等運算子，而我們只處理由雜訊所造成的退化。而第二個情形是檢視在 H 和 η 都存在的情況進行影像的復原。

## 雜訊模型

雜訊模型的特性與效應對影像的復原是很重要的，而雜訊模型大致可以分為兩個基本類型：空間域中的雜訊（由雜訊機率密度所描述）和頻率領域中的雜訊（由雜訊的各種傅立葉性質所描述）

### 常見的雜訊模型

* 高斯（Gaussian）

$$
p(z)={\frac {1}{{\sqrt {2\pi }}\sigma }}e^{-(z-\mu )^{2}/2\sigma ^{2}}
$$

* 雷利（Rayleigh）

 $p(z)={\frac {2}{b}}(z-a)e^{-(z-a)^{2}/b}$ ,for z≥a

* Erlang,Gamma(a,b)

 $p(z)={\frac {a^{b}z^{(}b-1)}{(b-a)!}}e^{(}-az)$ ,for z≥0

* 指數型態（Exponential）

 $p(z)=e^{-az}$ ,for z≥0

* 鹽與胡椒（salt-and-pepper）

$$
p(z)=P_{a}\sigma (z-a)+P_{b}\sigma (z-b)
$$

### 估計雜訊參數

通常藉由分析影像的傅立葉頻譜來估測週期性雜訊的參數。週期性雜訊有產生頻率尖波的傾向，即使用肉眼也通常能夠檢視。在雜訊間波很顯著的情況，或者當對有關干擾的頻率有些了解時，自動化的操作分析是可能的。

對於在空間域的雜訊，從感測器規格可大略知道 PDF 的參數，但是從樣本影像來估測這些參數是必要的。因此變成了從雜訊的平均值和變異數來估測解 PDF 所需的參數 a 和 b 的一個問題。設 z\_i 為在一影像中表示強度為準的一個離散隨機變數，並假設 p (zi)，i=0,1,2,3,…,L-1 為相對應的正規化的直方圖，其中 L 是可能的強度值數目。一個值方圖的成分 p (zi) 是強度值 zi 發生機率的一個估測，且值方圖可以視為強度 PDF 的一個近似。

## 空間域濾波

### 平均濾波器

* 算術平均濾波器 Arithmetic Mean Filter
* 幾何平均濾波器 Geometric Mean Filter
* 調和平均濾波器 Harmonic Mean Filter
* 反調和平均濾波器 Contraharmonic Mean Filter

### 排序統計值（Order-Statistics）濾波器

* 中值濾波器 Median Filter
* 最大值最小值濾波器 Max and Min Filter
* 中間點濾波器 Midpoint Filter
* Alpha 修整平均濾波器 Alpha-trimmed Mean Filter

### 適應性空間（Adaptive）濾波器

#### 適應性中值濾波器（adaptive median filter）

Zmin=Sxy 中的最小強度值

Zmax=Sxy 中的最大強度值

Zmed=Sxy 中的中間強度值

Zxy= 在座標 (x,y) 處的強度值

_Level A:_

A1 = zmed – zmin

A2 = zmed – zmax

If A1 > 0 AND A2 < 0, Go to level B

Else increase the window size

If window size ≤ Smax repeat level A

Else output zxy

_Level B:_

B1 = zxy – zmin

B2 = zxy – zmax

If B1 > 0 AND B2 < 0, output zxy

Else output zmed

## 以頻率域濾波降低週期性雜訊

週期性雜訊經常以在傅立葉頻譜中可見的脈衝狀串集呈現。對這些成分濾波的主要方法是經由帶陷濾波。n 階發特沃斯帶陷濾波器的轉移函數為： $H(u,v)={\cfrac {1}{1+\left[{\cfrac {D_{0}^{2}}{D_{1}(u,v)D_{1}(u,v)}}\right]^{n}}}$

其中

$$
D_{1}(u,v)=[(u-M/2-u_{0})^{2}+(v-N/2-v_{0})^{2}]^{\frac {1}{2}}
$$

且

$$
D_{2}(u,v)=[(u-M/2+u_{0})^{2}+(v-N/2+v_{0})^{2}]^{\frac {1}{2}}
$$

其中 (u0,v0) 且依對稱性 (-u0,-v0) 是「凹陷」的位置，而 D0 是他們半徑的量測。

### 帶斥濾波器（Bandreject filters）

* 理想帶斥濾波器

$$
H(u,v)={\begin{cases}1&{\mbox{if }}D(u,v)<D_{0}-{\frac {W}{2}}\\0&{\mbox{if }}D_{0}-{\frac {W}{2}}\leqslant D(u,v)\leqslant D_{0}+{\frac {W}{2}}\\1&{\mbox{if }}D(u,v)>D_{0}+{\frac {W}{2}}\\\end{cases}}
$$

* 發特沃斯帶陷濾波器 Butterworth Bandreject Filter

$$
H(u,v)={\cfrac {1}{1+\left[{\cfrac {D(u,v)W}{D^{2}(u,v)D_{0}^{2}(u,v)}}\right]^{2n}}}
$$

* 高斯帶斥濾波器 Gaussian Bandreject Filter

$$
H(u,v)=1-e^{-{\frac {1}{2}}\left[{\cfrac {D^{2}(u,v)D_{0}^{2}(u,v)}{D(u,v)W}}\right]^{2}}
$$

### 帶通濾波器（Bandpass filters）

$$
H_{BP}(u,v)=1-H_{BR}(u,v)
$$

### 凹口型濾波器（Notch Filter）

$$
H_{NP}(u,v)=1-H_{NR}(u,v)
$$

* 理想凹口型濾波器 Ideal Notch Reject Filter

$$
H(u,v)={\begin{cases}0&{\mbox{if }}D_{1}(u,v)\leqslant D_{0}{\mbox{ or }}D_{1}(u,v)\leqslant D_{0}\\1&{\mbox{otherwise }}\end{cases}}
$$

* 發特沃斯凹口型濾波器 Butterworth Notch Reject Filter

$$
H(u,v)={\frac {1}{1+\left[{\cfrac {D_{0}^{2}}{D_{1}(u,v)D_{2}(u,v)}}\right]^{n}}}
$$

* 高斯凹口型濾波器 Gaussian Notch Reject Filter

$$
H(u,v)=1-e^{-{\frac {1}{2}}\left[{\cfrac {D_{1}(u,v)D_{2}(u,v)}{D_{0}^{2}(u,v)}}\right]^{2}}
$$

## 直接反濾波

我們復原一張退化影像所能採取的最簡單方法是形成形式如下的一個估測： ${\hat {F}}(u,v)={\frac {G(u,v)}{H(u,v)}}$ 然後藉由 ${\hat {F}}(u,v)$ 的反傅立葉轉換獲得域個影像的相對應估測，這個方法被稱為反濾波（inverse filtering），由影像復原模型，我們可以將我們的估測表示成： ${\hat {F}}(u,v)=F(u,v)+{\frac {N(u,v)}{H(u,v)}}$ 由此式可知，即使我們確切的知道 H (u,v)，我們仍無法復原 F (u,v)，因為雜訊分量是一個它的傅立葉轉換 N (u,v) 未知的隨機函數。此外通常實際上有一個問題是函數 H (u,v) 有許多零點。即使 N (u,v) 這一項可忽視，將他除以 H (u,v) 幾乎為零的值會主宰復原的估測。

試圖反濾波的典型方法是形成比值 ${\hat {F}}(u,v)={\frac {G(u,v)}{H(u,v)}}$ ，然後限制獲得這個反濾波的頻率範圍到「接近」原點的頻率。此想法是 H (u,v) 中零點比較不可能再接近原點處發生，因為通常轉換的大小在該區域中有其最高值。有許多基調的變形，其中在 H 是零或靠近零的 (u,v) 處特別處理。這種方法有時稱為虛擬反（pseudoinverse）濾波。

## Wiener 濾波

Wiener 濾波尋求使以下統計誤差函數最小化的估測 ${\hat {f}}$  :

$$
e^{2}=E{(f-{\hat {f}}^{2})}
$$

其中 E 是期望值運算子而 f 是未退化的影像。此表示式在頻率域中的解為： ${\hat {F}}(u,v)=\left[{\frac {1}{H(u,v)}}{\frac {\left|H(u,v)\right|^{2}}{\left|H(u,v)\right|^{2}+S_{\eta }(u,v)/S_{f}(u,v)}}\right]G(u,v)$

其中

 $H(u,v)$ = 退化函數

$$
\left|H(u,v)\right|^{2}=H*(u,v)H(u,v)
$$

 $H*(u,v)=H(u,v)$ 的共軛複數

 $S_{\eta }(u,v)$ = 雜訊方功率頻譜

 $S_{f}(u,v)$ = 未退化影像的功率頻譜

比值 $S_{\eta }$ (u,v)/ $S_{f}$ (u,v) 稱為雜訊對訊號功率比，可以看出對所有 u 和 v 的相關值，如果雜訊功率頻譜為零，則此比值成為零，而 Wiener 濾波器簡化成在反濾波器。

## 參考資料

* Gonzalez, Rafael C., Richard Eugene Woods, and Steven L. Eddins. Digital image processing using MATLAB.. Pearson Education India. 2004.
* Wiener, Norbert (1949), _Extrapolation, Interpolation, and Smoothing of Stationary Time Series_. New York: Wiley. ISBN 0-262-73005-7
* Brown, Robert Grover and Patrick Y.C. Hwang (1996) _Introduction to Random Signals and Applied Kalman Filtering_. 3 ed. New York: John Wiley & Sons. ISBN 0-471-12839-2