[DIP] Deconvolution: Richardson-Lucy deconvolution algorithm

Deconvolution : 
Richardson-Lucy deconvolution algorithm

Image Restoration의 대표적인 예이기도 함.

 

Blurring을 결정하는 PSF (blur kenel이라고도 불림)가 알려진 경우 사용되는 non-blind deblurring algorithm의 대표적 기법.

• The algorithm is based on a PSF (Point Spread Function), where PSF is described as the impulse response of the optical system.
• The blurred image is sharpened through a number of iterations, which needs to be hand-tuned.
• 단점이라면, iteration을 많이 수행할수록 noise성분이 강화된다는 점임.

Degradation and Restoration for Image

 

RL deconv를 사용할 때,
측정된 이미지(=PSF에 의해 blur된 image)와
적용된 PSF는 알려져 있음.

 

Richardson-Lucy deconvolution을 수식으로 표현하면 다음과 같음.

$$\hat{W}_{i,r+1} = \hat{W}_{i,r}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}\hat{W}_{j,r}}
\quad r = 0,1,2,\dots$$

• $W$ : 원본 이미지, deconvolution의 결과물로 얻고자 하는 이미지.
• $\hat{W}$ : doconvolution의 결과 이미지. estimated image
• $H$ : 촬영된(측정된) 이미지. PSF로 blurring이 이루어진 이미지.
• $i,j, k$ : 영상 pixel좌표
• $r$ : iteration number
• $P_{j,k}$ : PSF.
• $W_j$가 $H_k$로 관측될 확률은 $P(H_k|W_j)=P_{j,k}$임.

가운데의 Noisy data는 PSF로 blurring이 이루어지고 Poission Noise가 추가된 것임. 맨 오른쪽에서 RL deconvolution의 결과를 보면, 선명해지긴 했지만, 영상의 경계부분에서 noise가 심해지고, 전체적인 noise들이 증폭된 것을 확인할 수 있음.

---

유도

측정된 영상의 $k$ pixel 의 intensity가 $H_k$ 일 경우,
원래 영상의 $i$ pixel의 intensity가 $W_i$일 확률은 다음과 같음.

$$P(W_i|H_k) = 
  \frac{P(H_k|W_i)P(W_i)}{\sum_{j} P(H_k|W_j)P(W_j)} \tag{1}$$

 

이는 원본 이미지 $W$에서 광자(photon)가 나오고, 이를  촬영된 이미지 $H$에서 측정하는 것으로 생각할 수 있음.

• Photon으로 생각하면,
  $P(W_i|H_k)$는 $H_k$에 해당하는 검출기에서 측정된 photon 중 얼마 만큼이 $W_i$로부터 온 것인지를 의미.
• $P(H_k|W_i)P(W_i)$는
  $W_i$에서 발산되는 photon이 $H_k$에 해당하는 검출기에서 검출될 확률을 의미함.
• $\sum_{j} P(H_k|W_j)P(W_j)$은
  $W$(=원래영상)에서 나오는 photon이 $H_k$에 해당하는 검출기에서 검출될 확률임.

결국, 위 식은 $H$의 $k$에 해당하는 pixel에서 검출된 모든 photon들 중에서 : $H_k$
$W$의 $i$에서 방출된 광자가 얼마인지($W_i$)를 나타내는 확률이 됨: $P(W_i|H_k)$

 

이는 다음과 같이 Bayesian Theorem을 따름.

$P(H_k|W_i)$는 likelihood로 PSF에 해당하며,
이 PSF에 의해 $W_i$의 intensity가 여러 $H_k$에 나누어지게 된다.

 

고로 다음이 성립함.

$$P(W_i) = 
  \sum_{k} P(W_i\cap H_k)=\sum_{k} \color{red}{P(W_i|H_k)}P(H_k) \tag{2}$$

위의 식은 조건부확률의 정의로부터 유도된 것이다.

 

$$P(W_i|H_k)=P(W_i\cap H_k)/P(H_k)$$

모든 $H_k$에 대해 $W_i$로부터 얻어질 확률을 다 더하면 전체 $W_i$의 확률이 구해짐 (← 식(2)).

 

식 1을 식 2에 대입하면 다음과 같음.

$$P(W_i) = 
  \sum_{k} \color{red}{\frac{P(H_k|W_i)\color{blue}{P(W_i)}}{\sum_{j} P(H_k|W_j)P(W_j)}}P(H_k)$$

 

$P(W_i)$가 결국 구해야 하는 값(=원본이미지)으로,
이는 iterative method 등을 통해 구하게 되며 다음과 같은 update 식으로 정리할 수 있음: $r$ 은 iteration number.

$$ P_{r+1}(W_i) = 
  \color{blue}{P_r(W_i)}\sum_{k} \frac{P(H_k|W_i)P(H_k)}{\sum_{j} P(H_k|W_j)P_r(W_j)} \quad
  r = 0,\ 1,\  2,\  \dots$$

초기값으로는 일반적으로 균일한 값을 가지는 영상을 사용함.

 

$W$의 모든 광자가 $H$에 전달된다고 가정 하에 각각의 확률을 정리하면 다음과 같음.

$$\begin{aligned}P(H_k|W_i) &= P_{i,k}\\P(W_i) &= \frac{W_i}{\sum_l W_l} \\
P(H_k) &= \frac{H_k}{\sum_l H_l} = \frac{H_k}{\sum_l W_l} 
\end{aligned}$$

 

위 식을 update식에 대입하면 RL deconvolution의 update식이 유도됨.

$$\begin{aligned}
    \frac{W_{i,r+1}}{\sum_l W_l} &= 
    \frac{W_{i,r}}{\sum_l W_l}
    \sum_{k} \frac{P_{i,k} \frac{H_k}{\sum_l W_l}}
    {\sum_j P_{j,k}\frac{W_{j,r}}{\sum_l W_l}} \\
    \iff \hat{W}_{i,r+1} &= \hat{W}_{i,r}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}\hat{W}_{j,r}}
    \quad r = 0,\ 1,\  2,\  \dots 
\end{aligned}$$

---

Convolution Form

위의 식은 영상시스템이 Linear Shift Invariant라고 가정한 경우 (PSF가 모든 pixel에 대해 동일)이며,
이를 convolution ($\otimes P$)과 correlation ($\otimes P^*$)으로 표현하면 다음과 같음.

 

$$\hat{W}_{r+1}=\hat{W}_{r}\left( \frac{H}{\hat{W}_{r}\otimes P}\otimes P^*\right)$$

• $\otimes$ : convolution (영상이므로 2-D convolution)
• $P^*$ : flipped PSF (horizontally and vertically)이며 결국 $\otimes P^*$는 PSF와 cross correlation을 의미함.

 

---

References

• Richardson, William Hadley. 1972, JOSA, 62, 55–59
• Lucy, L. B. 1974, Astronomical Journal, 79, 745–754

---

관련 ipynb파일

https://gist.github.com/dsaint31x/5c9018946cb775cf94ef6c9aed21942e

DIP_deconv_richardson_lucy.ipynb

---