發佈於 作者: adahbingee

Image Deconvolution: Frequency domain and Wiener filtering

Model

上一篇提到基本的幾個做法，本篇主要介紹其中 frequency domain的解法，並且解析一些實作上遇到的問題，繼續從 convolution model 開始

\[\begin{equation}
\label{eq:model}
I=L\otimes{P}
\end{equation}\]

Naive method

回來看式\(\eqref{eq:model}\)，其中\(I\)和\(P\)已知求\(L\)。左右各做FFT，\(\otimes\)在 frequency domain 以\(\cdot\)表 inner product

\[
\begin{equation}
F(I)=F(L)\cdot{F(P)}
\end{equation}
\]

求\(L\)很簡單，除掉就好，注意這邊是複數除法參考 Appendix A

\[
\begin{equation}
F(L)=\frac{F(I)}{F(P)}
\label{eq:naive}
\end{equation}
\]

最後再 inverse FFT 返回影像

\[
\begin{equation}
L=F^{-1}{(F(L))}
\end{equation}
\]

看起來似乎很好實做，算出 PSF 與 \(I\) 的 spectrum ，然後點對點做複數除法，最後再行 IFFT 求解即可。(前提是要小心 spectrum 對齊和用 zero padding 處理 warp-around error 的問題)

Divide by zero problem

但如果\(F(P)\)中有0項或很小的數字，\(F(I)\)項除以\(F(P)\)的值將會被放大到近無窮大。如果\(F(I)\)是完美的就沒什麼問題，但\(F(I)\)含有雜訊的成分時，雜訊就會被放大很多，屆時返回的影像只會是充滿雜訊的圖片。在 deblur 的例子中，\(P\)基本上是一個low pass filter，也就是高頻的 magnitude 接近0，除上一個接近0的數字將會不合理的放大雜訊，這不是我們要的結果! 怎麼辦?

左圖為正確 deblur 影像，右圖為 naive method 返回的 deblur 影像

一個直觀的想法，把那些太大的 noise 頻率壓下來，乘上一個 scalar term ，於是有了下面的式子

\[
\begin{equation}
\hat{F}(L)=\frac{F(I)}{F(P)}\left(\frac{|F(P)|^2}{|F(P)|^2+\lambda}\right),\lambda{\ge}0
\label{eq:wiener}
\end{equation}
\]

當\(\lambda\)接近0時，後面的乘項會越接近1；反之則越小。而\(|F(P)|^2\)和\(\lambda\)的影響力關係，假設\(\lambda\)固定，當\(|F(P)|^2\)越小，則\(\lambda\)影響力越大，恰好可以壓抑那些高頻低 magnitude 的部分。

舉例\(\lambda=0.5\)當\(|F(P)|^2=1\)和\(|F(P)|^2=2\)時分別 scalar 為0.67和0.8，是不是有點感覺了?這麼美麗的式子怎來的─ Wiener deconvolution。

如 \(\eqref{eq:e1}\) 所示，可以看做 Naive method 變形，使用修改過的 PSF \(F(G)\)，部分 magnitude 被墊高使其不會除0。當 \(\lambda=0\) 時退化為原本的 Naive method。

Wiener deconvolution

Wiener deconvolution假設一個additive noise \(N\)被應用在式\(\eqref{eq:model}\)

\[
\begin{equation}
I=L\otimes{P}+N
\end{equation}
\]

在 frequency domain 表示

\[
\begin{equation}
F(I)=F(L)\cdot{F(P)}+F(N) \label{eq:addmodel}
\end{equation}
\]

\(\hat{F}(L)\)表示我們估測的latent image，不可避免的已經混入 noise ，其中\(F(G)=F(P)^{-1}\)

\[
\begin{equation}
\hat{F}(L)=\frac{F(I)}{F(P)}=F(I){\cdot}F(G) \label{eq:e1}
\end{equation}
\]

我們的目標就是找一個 \(F(G)\) 最小化\(\left|F(L)-\hat{F}(L)\right|^2\)的期望值\(E\)，這也就是 Wiener filtering 的核心思想

\[
\begin{equation}
\arg\min_{F(G)}{E\left|F(L)-\hat{F}(L)\right|^2}
\end{equation}
\]

下面就是具體做法推導: 把 quadratic form 乘開後，稍微整理一下，一次微分找最小值，最後得到 \(F(G)\)。對於作者這外行人來說需要巨大勇氣!要複習一下複數乘法，微分連鎖律

\[
\epsilon (L)=E\left|F(L)-\hat{F}(L)\right|^2
\]

將\(\hat{F}(L)\)以\(\eqref{eq:e1}\)替換

\[
\begin{align*}
\epsilon(L)&=E\left|F(L)-\hat{F}(L)\right|^2 \\
&=E\left|F(L)-F(I){\cdot}F(G)\right|^2
\end{align*}
\]

將\(F(I)\)以式\(\eqref{eq:addmodel}\)替換

\[
\begin{align*}
\epsilon(L)&=E\left|F(L)-\hat{F}(L)\right|^2 \\
&=E\left|F(L)-F(I){\cdot}F(G)\right|^2 \\
&=E\left|F(L)-\left[F(L)\cdot{F(P)}+F(N)\right]\cdot{F(G)}\right|^2 \\
&=E\left|F(L)-F(L)\cdot{F(P)}\cdot{F(G)}+F(N)\cdot{F(G)}\right|^2 \\
&=E\left|F(L)\left[1-F(P)\cdot{F(G)}\right]+F(N)\cdot{F(G)}\right|^2
\end{align*}
\]

這裡我們需要一點勇氣，把平方展開，注意 complex 。並且假設 noise 與 latent image 是 independent 的，所以可以忽略所有\(F(N)\cdot{F(L)}\)項。其中 \(*\) 代表 complex conjugate

\[
\begin{align*}
\epsilon(L)&=E\left|F(L)-\hat{F}(L)\right|^2 \\
&=E\left|F(L)-F(I){\cdot}F(G)\right|^2 \\
&=E\left|F(L)-\left[F(L)\cdot{F(P)}+F(N)\right]\cdot{F(G)}\right|^2 \\
&=E\left|F(L)-F(L)\cdot{F(P)}\cdot{F(G)}+F(N)\cdot{F(G)}\right|^2 \\
&=E\left|F(L)\cdot\left[1-F(P)\cdot{F(G)}\right]+F(N)\cdot{F(G)}\right|^2 \\
&=E\left\{\left[1-F(P)\cdot{F(G)}\right]^{*}\cdot\left[1-F(P)\cdot{F(G)}\right]\cdot{F(L)^{*}}\cdot{F(L)}+F(N)\cdot{F(N)^*}\cdot{F(G)}\cdot{F(G)^*}\right\} \\
&=E\left\{\left[1-F(P)\cdot{F(G)}\right]^{*}\cdot\left[1-F(P)\cdot{F(G)}\right]\cdot{\left|F(L)\right|^2}+\left|F(N)\right|^2\cdot{F(G)}\cdot{F(G)^*}\right\}
\end{align*}
\]

一次微分找這個 quadratic form 的最小值為0，得到\(F(G)\)

\[
\begin{align*}
\frac{d\epsilon(L)}{dF(G)}=0 \\
-F(P){\cdot}\left[1-F(P)\cdot{F(G)}\right]^{*}\cdot{\left|F(L)\right|^2}+\left|F(N)\right|^2\cdot{F(G)^*} = 0 \\
-F(P)\cdot{\left|F(L)\right|^2}+\left|F(P)\right|^2{\cdot}F(G)^{*}\cdot{\left|F(L)\right|^2}+\left|F(N)\right|^2\cdot{F(G)^*} = 0 \\
-F(P)\cdot{\left|F(L)\right|^2}+F(G)^{*}\left(\left|F(P)\right|^2{\cdot}\left|F(L)\right|^2+\left|F(N)\right|^2\right) = 0 \\
\end{align*}
\]

\[
\begin{align*}
F(G)^{*}&=\frac{F(P)\cdot{\left|F(L)\right|^2}}{\left|F(P)\right|^2{\cdot}\left|F(L)\right|^2+\left|F(N)\right|^2} \\
F(G)&=\frac{F(P)^{*}\cdot{\left|F(L)\right|^2}}{\left|F(P)\right|^2{\cdot}\left|F(L)\right|^2+\left|F(N)\right|^2} \\
F(G)&=\frac{F(P)^{*}}{\left|F(P)\right|^2+\frac{\left|F(N)\right|^2}{\left|F(L)\right|^2}} \\
F(G)&=\frac{F(P)}{F(P)}\frac{F(P)^{*}}{\left|F(P)\right|^2+\frac{\left|F(N)\right|^2}{\left|F(L)\right|^2}} \\
F(G)&=\frac{1}{F(P)}\frac{\left|F(P)\right|^2}{\left|F(P)\right|^2+\frac{\left|F(N)\right|^2}{\left|F(L)\right|^2}}
\end{align*}
\]

到此就推導完成 \(F(G)\) ，將其帶入\(\eqref{eq:e1}\)，得到\(\hat{F}(L)\)

\[
\begin{equation}
\hat{F}(L)=\frac{F(I)}{F(P)}\left(\frac{|F(P)|^2}{|F(P)|^2+\lambda}\right),\lambda=\frac{\left|F(N)\right|^2}{\left|F(L)\right|^2}
\end{equation}
\]

How to choice \(\lambda\) – the signal to noise ratio

基本上 \(\lambda\) 是與頻率成正比，並且 \(\lambda=0\) 時 \(\eqref{eq:wiener}\) 會退化成 naive from \(\eqref{eq:naive}\)

PSF Anchor Problem

這邊由於 PSF 的 anchor 在中間，所以如果沒有修正整張圖會有半個 PSF 的circular shift，在 spatial domain 要先做前處理，由中間切割，上下左右對調。matlab function: psf2otf

psf2otf

// return frequency domain transformed point spread function with centering anchor by circular shift
static Mat psf2otf(const Mat &psf, const Size &size, Point2i anchor = Point2i(-1, -1)) {
    Mat paddedPsf(size.height, size.width, CV_64FC1, Scalar(0));

    // circular shift psf anchor to position (0, 0)
    if (anchor.x < 0 || anchor.y < 0 || anchor.y >= psf.rows || anchor.x >= psf.cols) {
        anchor.x = psf.cols / 2;
        anchor.y = psf.rows / 2;
    }

    // psf sub-region
    // A B
    // C O
    Mat A = psf(Range(0, anchor.y), Range(0, anchor.x));
    Mat B = psf(Range(0, anchor.y), Range(anchor.x, psf.cols));
    Mat C = psf(Range(anchor.y, psf.rows), Range(0, anchor.x));
    Mat O = psf(Range(anchor.y, psf.rows), Range(anchor.x, psf.cols));

    // target sub-region
    // O C
    // B A
    Mat tA = paddedPsf(Range(paddedPsf.rows-A.rows, paddedPsf.rows), Range(paddedPsf.cols-A.cols, paddedPsf.cols));
    Mat tB = paddedPsf(Range(paddedPsf.rows-B.rows, paddedPsf.rows), Range(0, B.cols));
    Mat tC = paddedPsf(Range(0, C.rows), Range(paddedPsf.cols-C.cols, paddedPsf.cols));
    Mat tO = paddedPsf(Range(0, O.rows), Range(0, O.cols));

    A.copyTo(tA);
    B.copyTo(tB);
    C.copyTo(tC);
    O.copyTo(tO);

    Mat otf;
    cv::dft(paddedPsf, otf, DFT_COMPLEX_OUTPUT);

    return otf;
}

Border Boundary Artifacts

參考[1]，模糊的圖像在拍攝時是有混合了來自於不在畫面的影像 (超過FOV取景範圍)，由於圖片在 frequency domain 是無限重複延伸的，所以在 deconvolution 時會取到對邊的資訊來計算，如果畫面邊界一邊含有高頻資訊，對邊卻沒有，就會在對邊產生不合理的 ringing。

這篇論文旨在單邊建立一個低頻平滑漸層連接兩個邊界 (ABC區塊，寬度僅需 PSF 大小)，可以使結果減少畫面邊界的 ringing；舉例來說 B 區塊左邊是較接近 O 區塊右側， B 區塊右邊較接近 O 區塊左側，依此類推。畫面中物體強 edge 附近的 ringing 是其他原因 (ill-posed problem)。

RBA2

RBA

實作上我是一個傻瓜懶人招，直接 double size 原圖作 mirror repeat，也得到類似的效果XD，注意這個做法很浪費資源。

Reference

  1. Non-blind Image Deconvolution with Adaptive Regularization, Jong-Ho Lee et al. PCM 2010
  2. https://en.wikipedia.org/wiki/Wiener_deconvolution

Appendix A.

在複數空間的乘法、除法

\[(x+yi)=(a+bi)(c+di)\]

\[(a+bi)=\frac{x+yi}{c+di}\]

\[a=\frac{cx+dy}{c^2+d^2}\]

\[b=\frac{-dx+cy}{c^2+d^2}\]

在〈Image Deconvolution: Frequency domain and Wiener filtering〉中有 1 則留言

  1. 自動引用通知: Image Deconvolution: Optimization-based algorithms | Computer Vision Notes

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *