디지털 영상 처리 - 제한된 최소 제곱 필터링

 안녕하세요. 지난 포스팅의 디지털 영상 처리 - 최소 평균 오차(Wiener) 필터링 구현하기에서는 MATLAB에 이미 구현된 deconvwnr 함수를 이용하였습니다. 그러나 deconvwnr의 치명적인 단점이 있습니다. 바로 어떤 커널로 블러링이 되었는 지 알아야 디컨볼루션을 진행할 수 있었습니다. 실제로 다시 한번 Wiener 필터링 수식 상기해보도록 하겠습니다.

 

$$\hat{F}(\mu, \nu) = \left[\frac{1}{H(\mu, \nu)} \cdot \frac{\left|H^{2}(\mu, \nu)\right|}{\left|H(\mu, \nu)\right|^{2} + \frac{S_{\eta}(\mu, \nu)}{S_{f}(\mu, \nu)}}\right]G(\mu, \nu)$$

 

 여기서 저희는 열화되지 않은 영상과 노이즈의 전력 스펙트럼 $S_{f}(\mu, \nu), S_{\eta}(\mu, \nu)$을 알아야한다는 문제가 있습니다. 물론 모르기 때문에 적당하게 $\frac{S_{\eta}(\mu, \nu)}{S_{f}(\mu, \nu)} = K$의 상수라고 두고 진행할 수도 있습니다. 이와 같은 근사법은 좋은 결과가 나오기도 하지만 항상 해결책은 아닙니다. 따라서 노이즈에 대한 최소한의 정보만을 가지고 영상을 복원할 수 있는 방법을 생각해보도록 하겠습니다. 

 

 보통 부가된 노이즈에 대한 평균과 분산은 열화된 영상으로부터 충분히 계산될 수 있습니다. 따라서 이제부터 설명드릴 방법은 노이즈의 평균과 분산만을 가지고 영상을 복원하는 방법에 대해서 설명하도록 하겠습니다. 일단 열화 과정을 모델링 및 컨볼루션 정리를 통해 저희는 아래의 식을 정의하도록 하겠습니다. 

 

$$\mathbf{g} = \mathbf{H}\mathbf{f} + \mathbf{\eta}$$

 

 여기서 $\mathbf{g}$를 만드는 방법에 대해서 설명해보도록 하겠습니다. 일단 어떤 열화된 영상 $g(x, y)$가 $M \times N$ 크기를 가진다고 하겠습니다. 그러면 $g(x, y)$의 첫번째 행의 $N$ 개의 값을 하나의 벡터 $\mathbf{g}_{1}$라고 하겠습니다. 이와 같이 정의하면 $i$번째 행을 벡터로 나타내면 $\mathbb{g}_{i}$가 됩니다. 이렇게 정의된 각 행 벡터들을 일렬로 나열한 것이 $\mathbf{g}$입니다. 따라서 $\mathbf{g}$은 $MN \times 1$의 크기를 가진다는 것을 알 수 있죠.

 

 그리고 $\mathbf{f}$와 $\mathbf{\eta}$는 $\mathbf{g}$와 동일한 크기를 가져야하기 때문에 둘 다 $MN \times 1$의 크기를 가집니다. 그렇다면 $\mathbf{H}$는 어떨까요? 행렬곱을 아신다면 아주 간단하게 $MN \times MN$이 되야함을 알 수 있습니다. 

 

 이 구조를 보시면 영상 복원이라는 것은 단순한 행렬곱으로 이루어진다는 것을 관찰할 수 있습니다. 하지만, 그 연산 비용은 만만치 않죠. 예를 들어서 $M = N = 512$의 영상이 주어졌다고 가정하겠습니다. 그러면 $\mathbf{g}$, $\mathbf{f}$, $\mathbf{\eta}$은 $262,144 \times 1$의 크기를 가지고 $\mathbf{H}$은 $262,144 \times 262,144$ 크기를 가집니다. 꽤나 많은 연산이 소요될 것으로 예상할 수 있습니다. 문제는 이뿐만이 아닙니다. 기본적으로 열화행렬 $H$의 특징 중 하나는 노이즈에 매우 민감하게 반응한다는 것입니다. 따라서 저희는 이 두 가지 문제를 완화시키면서 최적의 영상 복원을 하는 것으로 목표가 바뀌게 됩니다. 

 

 한편 $H$의 노이즈 민감성을 줄이기 위한 한 가지 대응책은 Laplacian을 적용하는 것입니다. 따라서 저희는 아래와 같이 정의된 기준 함수 $C$를 보도록 하겠습니다. 

 

$$C = \sum_{x = 0}^{M - 1} \sum_{y = 0}^{N - 1} \left[\nabla^{2} f(x, y)\right]^{2}$$

 

 그리고 이 식은 $\|\mathbf{g} - \mathbf{H}\mathbf{\hat{f}}\|^{2} = \|\mathbf{\eta}\|^{2}$이라는 제한조건이 걸립니다. 아마도 이를 최적화 문제로 다시 쓰면 아래와 같이 쓸 수 있습니다. 

 

$$\begin{aligned} \min_{f(x, y)} \quad& C = \sum_{x = 0}^{M - 1} \sum_{y = 0}^{N - 1} \left[\nabla^{2} f(x, y)\right]^{2} \\ \textrm{s.t.}  \quad& \|\mathbf{g} -\mathbf{H}\mathbf{\hat{f}}\|^{2} = \|\mathbf{\eta}\|^{2} \end{aligned}$$

 

위 수식의 최적해는 이미 아래의 수식으로 구할 수 있다고 알려져 있습니다. 

 

$$\hat{F}(\mu, \nu) = \left[\frac{H^{*}(\mu, \nu)}{\left|H(\mu, \nu)\right|^{2} + \gamma \left|P(\mu, \nu)\right|^{2}}\right]G(\mu, \nu)$$

 

 여기서 $\gamma$는 위의 제한식을 맞추기 위해 도입된 파라미터이고 $P(\mu, \nu)$는 아래의 푸리에 변환을 의미합니다. 

 

$$p(x, y) = \begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end{bmatrix}$$

 

 위의 $p(x, y)$를 보면 Laplacian 커널과 동일하다는 것을 알 수 있습니다. 남은 문제는 $\gamma$의 선택입니다. 어떻게 선택해야 좋을까요? 가장 간단한 방법은 대화식으로 바꾸어가면 테스트해보는 방법입니다.

 

 

 위의 그림은 제한된 최소 제곱 필터링을 이용해서 얻은 결과입니다. 여기서 $\gamma$는 시각적으로 가장 좋은 결과가 되도록 수작업으로 선택하였다고 합니다. 이전 포스팅의 Wiener 필터링과 비교해보면 구조가 많이 살아있는 것을 볼 수 있습니다. 

 

하지만 여전히 $\gamma$를 수동으로 정하는 것에 불만이 있으실 겁니다. 이를 위해서 반복적으로 $\gamma$를 얻는 방법도 존재합니다. 

 

 가장 먼저 나머지 벡터 $\mathbf{r} = \mathbf{g} - \mathbf{H}\mathbf{\hat{f}}$라고 정의하도록 하겠습니다. 그리고 제한된 최소 제곱 필터링의 최적해로부터 $\hat{F}(\mu, \nu)$은 $\gamma$의 함수이기 때문에 $\mathbf{r}$ 역시 $\gamma$의 함수입니다. 

 

 한편, 1973년 Hunt의 논문에 의해서 위와 같이 정의된 $\mathbf{r}$의 norm이 $\gamma$에 대해 단조 증가 함수임이 증명되었습니다. 즉, $\phi(\gamma) = \mathbf{r}^{T}\mathbf{r} = \|r\|^{2}$이 단조 증가 함수입니다. 일전에 저희가 미리 도입한 제한식을 이용하면 $\mathbf{r}$과 $\mathbf{\eta}$ 사이에는 아래와 같은 관계식이 존재합니다. 

 

$$\|\mathbf{r}\|^{2} = \|\eta\|^{2} \pm a$$

 

 여기서 $a$는 정확도 계수(accuracy factor)로서 $a = 0$이라면 100% 정확하게 위의 제한식을 만족한다는 것을 의미합니다. 여기서 $\phi(\gamma)$가 단조 증가함수라는 조건에 의해서 저희가 원하는 $\gamma$를 찾는 것은 굉장히 쉽습니다. 

 

STEP1. $\gamma$의 초깃값을 저장한다.

STEP2. $\|\mathbf{r}\|^{2}$를 계산한다.

STEP3. $\|\mathbf{r}\|^{2} = \|\eta\|^{2} \pm a$를 만족하면 중지한다. $\|\mathbf{r}\|^{2} < \|\eta\|^{2} - a$이라면 $\gamma$를 증가시킨다. $\|\mathbf{r}\|^{2} > \|\eta\|^{2} + a$이라면 $\gamma$를 감소시킨다. 

 

 위의 알고리즘을 어느정도 반복하면 최적의 추정을 얻을 수 있습니다. 물론 이 방법이 아닌 Newton-Rapson 알고리즘과 같은 훨씬 더 좋은 알고리즘을 고려해볼 수도 있습니다.