디지털 영상 처리 - 투영에 의한 영상 재구성 2

 안녕하세요. 지난 포스팅의 디지털 영상 처리 - 투영에 의한 영상 재구성 1에서는 CT의 기본적인 원리인 역투영을 통한 3D 재구성에 대해서 설명하고 기본적인 수학 도구인 라돈 변환(Radon Transform)과 용어인 시노그램(Sinogram), 래미노그램(Laminogram)에 대해서 설명하였습니다. 기본적인 알고리즘은 주어졌지만, 아쉽게도 저희가 원하는 퀄리티가 나오지 않았습니다. 오늘은 이러한 재구성 퀄리티를 증가시키는 푸리에-박편 정리(Fourier-Slicing Theorem)를 설명하고 이를 활용한 개선된 재구성 알고리즘을 소개해드리도록 하겠습니다. 

 

1. 푸리에-박편 정리(Fourier-Slicing Theorem)

 일단, 지난 포스팅에서의 $\rho$ 변수에 대한 1D 푸리에 변환은 아래와 같이 정의됩니다. 참고로 $\rho$는 직교좌표계의 직선을 $x\cos{\theta} + y\sin{\theta} = \rho$로 변환하였을 때의 변수이며, 원점에서 직선까지의 가장 짧은 직선거리를 의미합니다. 

 

$$G(\omega, \theta) = \int_{-\infty}^{\infty} g(\rho, \theta) e^{-j2\pi\omega\rho} \; d\rho$$

 

 여기서 $\omega$는 주파수 변수(Frequency Variable), $\theta$는 고정된 각도를 의미합니다. 만약, 푸리에 변환에 대한 자세한 설명을 원하신다면 아래의 링크에 순서대로 정리를 해놓았으니 꼭 참조해보시길 바랍니다. 

 

디지털 영상 처리 - 주파수 공간 필터링 기초 1

 그리고 지난 포스팅에서 저희가 배웠던 라돈 변환을 써보도록 하겠습니다. 

 

$$g(\rho, \theta) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \delta(x\cos{\theta} + y\sin{\theta} - \rho) \; dxdy$$

 

 이 식을 $\rho$에 대한 1D 푸리에 변환 수식에 그대로 대입하면 아래와 같습니다. 

 

$$G(\omega, \theta) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \delta(x\cos{\theta} + y\sin{\theta} - \rho) e^{-j2\pi\omega\rho} \; dxdyd\rho$$

$$ = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \left[\int_{-\infty}^{\infty} \delta(x\cos{\theta} + y\sin{\theta} - \rho) e^{-j2\pi\omega\rho} \; d\rho\right]\; dxdy$$

$$ = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-j2\pi\omega(x\cos{\theta} + y\sin{\theta})}\; dxdy$$

 

 위 식에서 마지막 줄은 임펄스 함수의 푸리에 변환 시 적용되는 선별 특성을 활용한 것입니다. 그 다음으로 $\mu = \omega\cos{\theta}, \nu = \omega\sin{\theta}$라고 놓으면 아래와 같이 변합니다. 

 

$$ G(\omega, \theta)= \left[\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-j2\pi(x\mu + y\nu)}\; dxdy\right]_{\mu=\omega\cos{\theta}, \nu=\omega\sin{\theta}}$$

$$ = \left[F(\mu, \nu)\right]_{\mu=\omega\cos{\theta}, \nu=\omega\sin{\theta}}$$

$$ = F(\omega\cos{\theta}, \omega\sin{\theta})$$

 

 이때, $F(\mu, \nu)$는 $f(x, y)$의 2D 푸리에 변환입니다. 이제 결과를 정리해보면 아래와 같습니다. 

 

$$ G(\omega, \theta)= F(\omega\cos{\theta}, \omega\sin{\theta})$$

 

 위 식을 보통 푸리에-박편 정리(Fourier-Slicing Theorem) 또는 투영-박편 정리(Projection-Slicing Theorem)이라고 합니다. 이 정리는 저희에게 어떤 결과를 알려줄까요? 이 식의 가장 처음은 특정 각도에서의 투영인 $g(\rho, \theta)$에 $\rho$에 대해서 1D 푸리에 변환을 적용한 것이였습니다. 그리고 그 결과는 직선 경로 상의 $f(x, y)$의 2D 푸리에 변환의 박편입니다. 

 

 

 위 그림을 해석하면 더욱 쉽게 이해할 수 있습니다. 특정 각도 $\theta$의 투영 결과인 $g(\rho, \theta)$는 왼쪽 그림의 "Projection"을 의미합니다. 그리고 이것을 1D 푸리에 변환을 하게 되면 주파수 공간에서 투영을 만들었을 때 사용했던 각도 $\theta$인 선을 따라서 $F(\mu, \nu)$를 따라서 값을 추출했다는 의미입니다. 

 

2. 평행-빔 필터링된 역투영을 이용한 재구성(Reconstruction Using Parallel-Beam Filtered Backprojection)

 이전 포스팅에서도 보았지만 직접 역투영을 적용하면 블러링이 너무 심해서 저희가 원하는 정도의 퀄리티가 나오지 않았습니다. 이를 해결하는 방법이 바로 "평행-빔 필터링된 역투영을 이용한 재구성"입니다. 이 과정은 역투영을 계산하기 전에 투영을 필터링함으로써 이루어집니다. 일단 $F(\mu, \nu)$의 2D 역푸리에 변환은 아래와 같습니다. 

 

$$f(x, y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(\mu, \nu) e^{j2\pi(x\mu + y\nu)} \; d\mu d\nu$$

 

 그리고 이전 절과 같이 $\mu = \omega\cos{\theta}, \nu = \omega\sin{\theta}$라고 놓으면 그 미분은 $d\mu d\nu = \omega d\omega d\theta$가 되고 이를 극좌표로 쓰면 아래와 같습니다. 

 

$$f(x, y) = \int_{0}^{2\pi} \int_{0}^{\infty} F(\omega\cos{\theta}, \omega\sin{\theta}) e^{j2\pi(x\cos{\theta} + y\sin{\theta})} \; \omega d\omega d\theta$$

 

 여기서 푸리에-박편 정리 $ G(\omega, \theta)= F(\omega\cos{\theta}, \omega\sin{\theta})$를 위 식에 대입하면 아래와 같이 바뀝니다. 

 

$$f(x, y) = \int_{0}^{2\pi} \int_{0}^{\infty} G(\omega, \theta) e^{j2\pi(x\cos{\theta} + y\sin{\theta})} \; \omega d\omega d\theta$$

 

 그 다음에 설명하지는 않았지만 라돈 변환의 1D 푸리에 변환의 특성인 $G(\omega, \theta + \pi) = -G(\omega, \theta)$를 활용하고 0 ~ 180도의 적분과 180 ~ 360도의 적분으로 나누어서 하면 아래와 같이 표현됩니다. 

 

$$f(x, y) = \int_{0}^{\pi} \int_{0}^{\infty} |\omega| G(\omega, \theta) e^{j2\pi\omega(x\cos{\theta} + y\sin{\theta})} \;  d\omega d\theta$$

 

 앗, 그런데 $\omega$에 대한 적분 항들 중  $x\cos{\theta} + y\sin{\theta} = \rho$를 만족하기 때문에 다시 한번 아래와 같이 쓸 수 있습니다. 

 

$$f(x, y) = \int_{0}^{\pi} \left[\int_{0}^{\infty} |\omega| G(\omega, \theta) e^{j2\pi\omega\rho} \;  d\omega \right]_{\rho = x\cos{\theta} + y\sin{\theta}}d\theta$$

 

 위 식은 $|\omega|$ 항이 곱해진 $G(\omega, \theta)$의 1D 역푸리에 변환이라고 볼 수 있습니다. 이 식을 해석해보면 $|\omega|$라는 함수로 필터링 시킨 것이라고 볼 수 있습니다. 그런데 $|\omega|$ 함수는 $-\infty$ ~ $\infty$ 까지 정의되기 때문에 역푸리에 변환이 정의되지 않습니다. 이러한 현상을 방지하기 위해서 저희는 어떤 필터링 함수를 사용할 것인지 결정해야합니다. 가장 좋은 방법은 박스 필터를 씌우는 것이겠죠. 하지만, 박스 필터는 주파수 공간에서 불연속인 구간을 만들기 때문에 그에 대응되는 영상 공간에서는 무한히 반복하는 물결 파동 현상이 발생할 수 있습니다. 따라서 저희는 어떤 부드러운 필터 함수를 사용해야만 합니다. 오늘 설명하고 있는 투영에 의한 재구성에서는 "해밍 윈도우(Hamming Window)" 또는 "한 윈도우(Hann Window)" 함수를 자주 사용합니다. 이 식은 아래와 같이 정의됩니다. 

 

$$h(\omega) = \begin{cases} c + (c - 1)\cos{\frac{2\pi\omega}{M - 1}} & 0 \le \omega \le (M  - 1)\\ 0 & \text{Otherwise}\end{cases}$$

 

 여기서 $c = 0.54$면 해밍 윈도우, $c = 0.5$면 한 윈도우라고 합니다. 두 윈도우의 가장 큰 차이점은 한 윈도우의 끝 점들이 0인라는 것외에는 없습니다. 실제로 대부분에 응용에서 두 윈도우의 차이를 구별하기는 어렵다고 하네요. 그리고 위 윈도우를 곱한 뒤 역푸리에 변환을 적용해주면 됩니다. 그러면 그냥 $|\omega|$를 곱해서 적용하는 것보다는 훨씬 나은 결과를 얻을 수 있습니다. 

 

 다시 돌아와서 전체 알고리즘을 정리해볼 수 있을 거 같습니다. 

 

1. 각 투영 $g(\rho, \theta)$의 1D 푸리에 변환을 계산한다.
2. 해밍 윈도우에 의해 곱해진 필터 함수 $|\omega|$를 각 푸리에 변환에 곱한다.
3. 각 필터링된 변환 결과에 대해서 1D 역푸리에 변환을 얻는다.
4. STEP3에서 얻은 모든 1D 역변환을 적분한다.(=더한다.)

 위 그림은 해밍 윈도우를 적용했을 때와 $|\omega|$를 적용했을 때의 차이를 보여주고 있습니다. (c) 그림을 보면 (d) 그림에 비해서 물결 파동 현상이 더 심하다는 것을 알 수 있습니다. 

 

 저희는 $f(x, y) = \int_{0}^{\pi} \left[\int_{0}^{\infty} |\omega| G(\omega, \theta) e^{j2\pi\omega\rho} \;  d\omega \right]_{\rho = x\cos{\theta} + y\sin{\theta}}d\theta$을 통해서 개선된 역투영을 할 수 있음을 알게 되었습니다. 하지만, 저희는 주파수 공간과 영상 공간 사이의 관계성을 보여주는 컨볼루션 정리(Convolution Theorem)을 여기에 적용해볼 수 있습니다. 이를 위해서 $s(\rho)$를 $|\omega|$의 1D 역푸리에 변환이라고 하면 저희는 아래와 같이 쓸 수 있습니다. 

 

$$f(x, y) = \int_{0}^{\pi} \left[\int_{0}^{\infty} |\omega| G(\omega, \theta) e^{j2\pi\omega\rho} \;  d\omega \right]_{\rho = x\cos{\theta} + y\sin{\theta}} \; d\theta$$

$$ = \int_{0}^{\pi} \left[ s(\rho) \star g(\rho, \theta) \right]_{\rho = x\cos{\theta} + y\sin{\theta}} \; d\theta$$

$$ = \int_{0}^{\pi} \left[ \int_{-\infty}^{\infty} g(\rho, \theta)s(x\cos{\theta} + y\sin{\theta})\right] \; d\theta$$

 

 마지막 줄은 컨볼루션의 정의를 직접적으로 활용하여 얻은 결과입니다. 위 식에서 의미하는 것은 각도 $\theta$에서의 각 개별 역투영들은 해당 투영 $g(\rho, \theta)$와 램프 필터 $s(\rho)$의 역푸리에 변환을 컨볼루션 시켜서 얻을 수 있습니다. 그리고 마찬가지로 역투영된 영상은 모든 역투영된 개별 영상들을 합하여 얻을 수 있습니다. 일반적으로 컨볼루션 연산 자체가 더 계산에서 효과적이기 때문에 현대 CT 시스템은 해당 접근법을 많이 사용한다고 합니다. 

 

3. 부채-빔 필터링된 역투영을 이용한 재구성(Reconstruction Using Fan-Beam Filtered Backprojection)

 지금까지는 평행-빔에 집중하였으나 최근 CT 시스템은 1개의 소스로부터 여러 개의 투영을 하는 부채-빔 시스템을 활용하고 있습니다. 

 

 위 그림에서 검출기 들은 소스의 반대편에 호를 따라서 배치됩니다. 그리고 $p(\alpha, \beta)$는 부채-빔 투영을 의미한다고 했을 때 $\alpha$는 중심 광선에 대해 측정된 특정 검출기의 각도 위치, $\beta$는 $y$-축에 대한 측정된 광원의 각도 변위를 의미합니다. 부채-빔의 광선은 이전에 언급한 평행-빔 영상화 과정과 마찬가지로 $L(\rho, \theta)$의 형태로 표현될 수 있습니다. 

 

 아마 몇몇 분들은 생각하셨겠지만, 여기서 바로 푸리에-박편 정리를 적용할 수 없습니다. 왜냐하면 푸리에-박편 정리 자체가 평행-빔 시스템을 중점에 두고 만들어진 정리이기 때문이죠. 따라서 이를 어떻게 사용할 수 있을 지 부터 알아보아야 합니다. 일단 위 그림에서 아래의 관계식이 성립한다는 것으로부터 시작됩니다. 

 

$$\theta = \beta + \alpha$$

$$\rho = D \sin{\alpha}$$

 

 여기서 $D$는 광원 중심으로부터 $xy$-평면의 원점까지의 거리를 의미합니다. 여기서 저희가 평면의 원점을 중심으로하는 반지름 $T$인 원형 영역 안에 있는 객체를 포착하려고 한다고 가정하겠습니다. 그러면 원점에서 빔까지의 거리 $\rho$는 $|\rho| > T$에 대해서 $g(\rho, \theta) = 0$이라고 볼 수 있습니다. 이는 $T$보다 더 큰 영역은 아예 고려 대상에서 제외하겠다는 의미입니다. 위와 같은 사실들을 조합하면 저희가 평행-빔을 기반으로한 컨볼루션 복원 수식은 아래와 같이 바뀝니다. 

 

$$ f(x, y) = \frac{1}{2}\int_{0}^{2\pi}  \int_{-T}^{T} g(\rho, \theta)s(x\cos{\theta} + y\sin{\theta} - \rho) \; d\rho d\theta$$

 

 그리고 저희는 현재 $\alpha, \beta$에 따라서 투영결과가 달라지기 때문에 $\alpha, \beta$에 대한 적분을 구해야합니다. 이를 위해서 먼저 좌표계를 극좌표 $(r, \phi)$로 바꾸도록 하겠습니다. 즉, $x = r\cos{\phi}, y = r\sin{\phi}$라고 두겠습니다. 그러면 아래와 같은 식을 얻을 수 있습니다. 

 

$$x\cos{\theta} + y\sin{\theta} = r\cos{\phi}\cos{\theta} + r\sin{\phi}\sin{\theta} = r\cos{\left(\theta - \phi\right)}$$

 

 그러면 아래와 같이 적분식이 변하게 됩니다. 

 

$$ f(x, y) = \frac{1}{2}\int_{0}^{2\pi}  \int_{-T}^{T} g(\rho, \theta)s(r\cos{\left(\theta - \phi\right)} - \rho) \; d\rho d\theta$$

 

 하지만 여전히 적분이 $\rho, \theta$에 대한 것이기 때문에 좌표변환을 위해 $\theta = \beta + \alpha$와 $\rho = D \sin{\alpha}$ 사실을 이용하면 됩니다. 그러면 아래와 같이 식이 변하게 됩니다. 

 

$$ f(x, y) = \frac{1}{2}\int_{-\alpha}^{2\pi - \alpha}  \int_{\arcsin{(-T/D)}}^{\arcsin{(T/D)}} g(D\sin{\alpha}, \alpha + \beta)s(r\cos{\left(\beta + \alpha - \phi\right)} - D\sin{\alpha}) D\cos{\alpha} \; d\alpha d\beta$$

 

 여기서 $d\rho d\theta = D\cos{\alpha} d\alpha d\beta$를 이용하였습니다. 일단, 위 식은 상당히 복잡합니다. 위 식을 간단하게 만들기 위해서 기하학적 요소를 첨가해보도록 하겠습니다. 

 

 

 위 그림을 보면서 설명해보도록 하겠습니다. 먼저, $\beta$의 범위를 고려해보겠습니다. $\beta$는 위 적분 식에 의해서 $\alpha$ ~ $2\pi - \alpha$까지 커버되고 있습니다. 그리고 이 범휘는 360도 모든 범위를 포함하고 있습니다. 그리고 $\beta$에 대한 모든 함수는 주기가 $2\pi$가 된다는 점을 고려하면 나머지 적분들의 한계들은 각각 $0$과 $2\pi$로 대체될 수 있습니다. 

 

 다음으로 $\arcsin{(T/D)}$는 $|\rho| > T$에 대한 최대값 $\alpha_{m}$을 가지고 그 너머에서는 $g = 0$이기 때문에 안쪽 적분의 한계들은 각각 $-\alpha_{m}$, $\alpha_{m}$으로 대체될 수 있습니다. 

 

 그리고 $L(\rho, \theta)$를 고려해보도록 하겠습니다. 이 선을 따라서의 부채-빔 광선의 합은 같은 선을 따라서의 평행-빔 광선의 합과 동일해야합니다. 그리고 이는 $(\alpha, \beta)$와 $(\rho, \theta)$의 대응값에 대해서 모든 광선 합에 대해서 참입니다. 따라서 $p(\alpha, \beta) = g(\rho, \theta)$임을 알 수 있고 $p(\alpha, \beta) = g(D\sin{\alpha}, \alpha + \beta)$가 됩니다.

 

 위의 모든 조건들을 고려하여 다시 수식을 정리하면 아래와 같이 쓸 수 있습니다. 

 

$$ f(x, y) = \frac{1}{2}\int_{0}^{2\pi}  \int_{-\alpha_{m}}^{\alpha_{m}} p(\alpha, \beta)\left[s(r\cos{\left(\beta + \alpha - \phi\right)} - D\sin{\alpha}\right] D\cos{\alpha} \; d\alpha d\beta$$

 

 이제부터는 위 수식이 부채-빔 시스템에서의 역투영 수식의 기본입니다. 이 식을 조금 더 다듬으면 컨볼루션 형태로 바꿀 수 있습니다. 

 

 위 그림에서 $r\cos{\left(\beta + \alpha - \phi\right)}  - D\sin{\alpha} = R\sin{\left(\alpha^{'} - \alpha\right)}$임을 증명이 되었다고 가정하겠습니다. 여기서 $R$은 광원으로부터 부채 광선의 임의의 점, $\alpha^{'}$는 이 광선과 중심 광선 사이의 각도를 의미합니다. 여기서 $R$과 $\alpha^{'}$가 $r, \phi, \beta$에 의해 결정되고 있습니다. 이 식을 다시 적분에 대입하면 아래의 수식을 얻을 수 있습니다. 

 

$$ f(x, y) = \frac{1}{2}\int_{0}^{2\pi}  \int_{-\alpha_{m}}^{\alpha_{m}} p(\alpha, \beta)\left[R\sin{\left(\alpha^{'} - \alpha\right)}\right] D\cos{\alpha} \; d\alpha d\beta$$

 

 그리고 $s(R\sin{\alpha}) = \left(\frac{\alpha}{R\sin{\alpha}}\right)^{2} s(\alpha)$임을 안다고 가정하면 다시 수식을 아래와 같이 고칠 수 있습니다. 

 

$$f(x, y) = \int_{0}^{2\pi} \frac{1}{R^{2}} \left[\int_{\alpha_{m}}^{\alpha_{m}} q(\alpha, \beta) h^{\alpha^{'} - \alpha} \; d\alpha\right] \; d\beta$$

 

 여기서 $h(\alpha) = \frac{1}{2} \left(\frac{\alpha}{\sin{\alpha}}\right)^{2} s(\alpha)$, $q(\alpha, \beta) = p(\alpha, \beta) D\cos{\alpha}$를 의미합니다. 그리고 위 식에서 안쪽 적분은 $q(\alpha, \beta), h(\alpha^{'} - \alpha)$ 사이의 컨볼루션임을 알 수 있습니다. 

 

 영상 복원과 관련된 내용은 여기까지 입니다! 다음 포스팅부터는 컬러 공간에 대한 이야기를 해보도록 하겠습니다.