배경지식

---

이번 장에서는 단층 영상을 frequency filtering 을 통해서 유의미한 결과물로 재조합하는 방법을 다룹니다.

X-ray 를 이용한 CT, 감마선을 이용한 PET는 대표적인 단층 촬영 영상 방식입니다.

CT vs PET

위의 두 영상은 같은 환자의 같은 부위를 촬영한 것입니다. CT 는 일반적으로 우리가 병원에서 자주 촬영하는 X-ray 와 원리가 유사하기 때문에 그 결과물도 비슷합니다. 병변(lesion) 의 유무나 위치를 알 수는 있지만 그것이 암 조직인지는 확신할 수가 없습니다.

PET 는 해부학적으로 디테일한 결과물을 보장하지는 않습니다. 위 사진에서 병변 부위를 제외한 모든 부분이 흐리게 표현됩니다. 하지만 병변 부위는 매우 밝고 선명한데, 이 곳에서 분자 활동이 활발하게 일어남을 의미하고, 이것으로 암 조직이라는 것을 확신할 수 있는 것입니다. 암조직은 끊임 없이 분열하니까!

CT 는 해부학적 구조를 잘 보여주기 때문에 anatomical imaging 로 분류되고, PET 는 molecular activity 라는 기능적인 측면을 잘 보여주기 때문에 functional imaging 로 분류됩니다.

질병의 조기 진단을 위해서 이러한 의료 영상을 얻는 것은 상당히 중요한데, CT 와 PET 모두 같은 원리의 이미지 재조합 알고리즘을 사용합니다.

CT

위 figure 는 1세대 CT의 원리를 설명하는 그림입니다. 앨런 코맥과 고드프리 하운스필드가 이 아이디어로 노벨 생리의학상을 공동 수상했습니다.

X-ray source 로부터 방출된 X-ray 는 인체를 통과하게 되는데, 밀도가 높은 부위를 지날 수록 감쇠 작용(attenuation)이 일어나 에너지 손실이 발생합니다. 그러므로 대상이 두껍거나 밀도가 높을수록 반대 쪽의 센서에는 낮은 X-ray value 가 검출될 것입니다. X-ray source 부터 센서까지의 선을 Line of Response(LOR) 이라고 부릅니다.

최근의 개선된 CT 방식입니다. 요즘은 센서를 많이 둬서 한 번에 많은 LOR 을 얻습니다. 특정 각도에서의 LOR 들의 완전한 집합을 projection 이라고 부릅니다.

여러 각도에서 얻은 projection 을 누적시켜 2D image 를 만든 것을 sinogram 이라고 부릅니다.

Sinogram 얻고자 하는 유의미한 결과물(단층 영상)

Sinogram 은 그 자체로서는 의미가 없습니다. 이것을 유의미한 단층 영상으로 바꾸는 과정이 이미지 재구성(image reconstruction)입니다. CT 기술의 발명 이후 50년 동안 많은 이미지 재구성 방식이 고안되었습니다. 이 방식들은 크게 수학적으로 계산하는 analytic method 와 루프를 계속해서 돌려서 정확도를 올리는 iterative method 로 나뉘는데, analytic method 가 가장 기본적인 방식이며 LTI system 을 이용할 수 있습니다. 최근에는 많은 계산량 때문에 iterative method 가 주목받고 있지만 여전히 analytic method 도 많이 사용됩니다. 오늘 다룰 방식도 analytic method 입니다.

CT 시스템의 기본 개념

---

1G CT

Source 에서 방출되는 X-ray 의 세기를 $I_{0}$ , detector 에서 받는 X-ray 의 세기를 $I_{d}$ 라고 합시다. Detector 가 받는 X-ray 는 수학적으로 어떻게 모델링할까요? X-ray의 에너지는 여러 개의 에너지로 구성된 스펙트럼입니다. 그것을 완벽하게 계산하기에는 수학적으로 상당히 복잡합니다. 그래서 검출되는 에너지는 mean value 를 사용합니다. 그것이 effective energy, $\bar{E}$ 입니다.

위 식에서 s 는 X-ray 축입니다. 축의 각도가 바뀌면 계산 결과도 달라지겠죠. Source 부터 detector 까지의 path 의 밀도들을 모두 적분해서 $e^{-1}$ 의 지수로 사용합니다. 거기에 X-ray 의 원래 에너지인 $I_{0}$ 을 곱하면 $I_{d}$ 가 나옵니다.

가장 중요한 것은 Projection 으로 표현되는 $g_{d}$ 를 구하는 것입니다. 앞서 적분 식에서도 볼 수 있었던 $\mu$ 는 선형 감쇠 계수(linear attenuation coefficient) 입니다. 이 $\mu$ 는 $s$ 에 대한 함수입니다. 결국 $I_{d}$ 를 구한 다음 $g_{d}$ 를 구해야하고, 그것이 X-ray 혹은 CT 에서 얻는 projection 영상입니다. $g_{d}$ 는 $\mu$ 에 대한 선적분입니다.

Line impulse

Line impulse 는 특정 선에 위치한 신호를 의미하는 개념으로, 임펄스 함수($\delta$)가 특정 선을 따라 집중된 형태입니다. 위 그림에서 line impulse $\delta_{\ell}(x, y)$ 는 특정 선 $L(\ell, \theta)$ 에만 값이 있고, 그 외의 모든 위치에서는 0이 됩니다. 여기서 수식

는 선 $L$ 을 따라가며 정의된 임펄스 함수입니다. 이 수식은 $x \cos \theta + y \sin \theta = \ell$ 이라는 조건을 만족하는 점에 대해 값이 있고, 그 외에는 0입니다. 즉, $x \cos \theta + y \sin \theta = \ell$ 라는 식은 평면 상에서 특정 각도 $\theta$ 에 위치한 선을 나타내며, 이 선에 위치하는 모든 점에서 임펄스 함수가 활성화 됩니다. 이 임펄스 함수는 선 위에서만 값을 가지며, 그 외 모든 영역에서는 0이 됩니다.

선 L 은 아래의 수식으로 표현됩니다.

$L(\ell, \theta)=\{(x, y)|x \cos \theta + y \sin \theta = \ell\}$

Line impulse 는 아래의 수식으로 표현됩니다.

$\delta_{\ell}(x,y)=\delta(x\cos\theta+y\sin\theta-\ell)$

X-ray 재조합

우리는 target 의 밀도 분포 함수인 $f(x, y)$ 를 구하는 것이 목적입니다. 하지만 CT 촬영으로 알 수 있는 것은 target 을 지난 X-ray 의 projection 을 log scale 로 conversion 한 $g(\ell, \theta)$ 입니다.

각각의 LOR 은 $L(\ell, \theta)=\{(x, y)|x \cos \theta + y \sin \theta = \ell\}$ 조건을 만족합니다. 그러므로 특정 $\ell$ 과 특정 $\theta$ 에 대한 line impulse 를 밀도 분포 함수인 $f(x, y)$ 에 취해주고 그것을 모든 $x$, $y$ 에 대해 중적분을 해주면 $g(\ell, \theta)$ 입니다. 고정된 $\theta$ 에서 계산한, 즉 특정 $\theta$ 에서 LOR 의 $g$ 가 바로 projection 입니다. 수학적으로는 Radon transform 이라고 부릅니다. 모든 $\ell, \theta$ 에서의 projection 을 구해서 누적시킨 것이 sinogram 입니다.

Sinogram

밀도가 존재하는 한 점을 촬영해서 얻은 sinogram 영상입니다. $\theta$ 의 방향을 유의하면서 (a) 를 봅시다. $\theta = 0$ 일 때는 회전 중심으로부터 밀도가 존재하는 부분이 왼쪽에 위치합니다. 그러므로 (c) 의 1번 부분이 중심으로부터 왼쪽에 위치합니다. $\theta = -\pi/4$ 일 때는 회전 중심으로부터 밀도가 존재하는 부분까지의 거리가 더 좁아졌기 때문에 (c) 의 2번 부분도 1번 부분에 비해 좀 더 중심으로 이동합니다. $\theta=-\pi/2$ 일 때는 회전 중심으로부터 밀도가 존재하는 부분의 위치가 오른쪽으로 바꼈습니다. 그러므로 (c) 의 3번 부분도 중심축으로부터 왼쪽이 아닌 오른쪽에 위치합니다.

원리를 이해했다면 위 그림도 쉽게 이해할 수 있을 것입니다.

위 그림에서는 각도에 따라 밀도를 가지는 부분들이 projection에 따라 가지런히 정렬되어 있을 때도, 산재되어 있을 때도 있습니다. 정렬되어 있을 때는 결과 영상에서 검은색과 흰색이 번갈아가며 나오는 형태를 보입니다. $\theta = \pi/2$ 부분을 주목해봅시다.

Back-projection

Back-projection 은 CT 스캔 후 이미지를 재구성하는 중요한 과정입니다. Sinogram 데이터만으로는 물체의 형태를 알 수 없기 때문에, 데이터를 반대 방향으로 투영하여 원래의 밀도 분포를 복원할 필요가 있습니다. Back-projection 을 통해 각 투영 데이터를 반대 방향으로 투사하여, 모든 각도에서의 데이터를 겹쳐서 원래의 밀도 분포를 복원합니다.

모든 각도에서 back-projection 을 수행하면서 projection 과 밀도가 존재하는 부분이 겹치는 곳이 누적될수록 더 밝게 표현됩니다. 하지만 가장 오른쪽 그림을 보면 알 수 있듯이 밀도가 존재하는 부분의 주변부도 약간은 겹칩니다. 이 때문에 주변부가 흐릿하게 보이는 문제가 발생합니다.

우리가 구하고자하는 밀도 분포 함수인 $f(x, y)$ 와 back-projection 으로 재구성한 밀도 분포 함수인 $f_{b}(x, y)$ 사이의 관계를 수학적으로 표현할 수 있으면 이 문제를 해결할 수 있지 않을까요?

Projection-Slice Theorem(Central Section Theorem)

1-D Fourier transform

Projection $g(\ell, \theta)$ 의 1D Fourier transform 을에 $G(\rho, \theta)$ 라고 합시다. 모든 $\theta$ 에서의 projection 을 Radon transform 식으로 표현할 수 있기 때문에 FT 식에 $g(\ell, \theta)$ 에 Radon transform 식을 대입합니다. 이렇게 하면 다음과 같은 1D FT 식을 얻을 수 있습니다.

2-D Fourier transform

2D FT 를 수행하고 $u, v$ 를 극좌표로 변환하면 다음과 같은 2D FT 식을 얻을 수 있습니다.

이는 1-D FT 의결과와 정확히 같습니다.

Back-projection 기법 1: Fourier method

Projection-slice theorem 의 원리를 활용해 여러 각도에에서 투영 데이터를 수집한 뒤, 각 투영에 대해 1D FT을 수행하면, 전체 2D FT 을 효과적으로 추정할 수 있습니다. 이는 각도별 1D FT 를 구한 후, 그 결과값들을 조합하여 전체 2D Fourier 공간을 복원할 수 있다는 의미입니다. 이 복원된 Fourier 공간을 역변환하면 원래 이미지인 $f(x, y)$ 가 재구성됩니다. 하지만 이 방법은 실제 CT 에서는 널리 사용되고 있지는 않습니다:

• 투영 데이터는 극좌표계로 수집되는데, 이를 2D Foureir 평면의 직교 좌표계에 맞추려면 interpolation 이 필요합니다. 이 interpolation 에서 데이터의 손실이 유발됩니다.
• 이 방법에는 조정 가능한 매개변수가 없기 때문에, 노이즈를 처리하기가 어렵습니다. 노이즈가 포함된 데이터를 그대로 사용해야 하기 때문에 재구성된 이미지 품질이 떨어질 수 있습니다.
• 2D IFT 를 수행하는데 시간이 많이 소요됩니다. 이는 빠른 처리가 필요한 CT 시스템에는 적합하지 않습니다.

Back-projection 기법 2: Filtered Back-projection

Filtered back-projection 은 Fourier method 와는 달리 interpolation 문제를 겪지 않습니다. 또한 2D IFT 과정이 없으므로 시간 효율도 더 좋습니다. 필터를 통해 노이즈를 제어할 수 있기 때문에 더 선명하고 안정적인 이미지를 얻을 수 있다는 장점도 있습니다. 위 계산을 거치면 결국 $f(x, y)$ 는 다음과 같은 수식으로 구할 수 있게 됩니다.

1. 1D 주파수 필터 적용: 각 투영 데이터 $G(\rho, \theta)$ 에 대해 1D 주파수 필터 $H(\rho)=|\rho|$ 를 적용합니다.
2. 1D IFT: 필터링된 투영 $|\rho|G(\rho, \theta)$ 에 대해 1D IFT 을 수행합니다. 이 과정은 주파수 영역에서의 필터링된 데이터를 다시 공간 영역으로 변환하여, 각 투영을 필터링된 상태로 준비하는 단계입니다.
3. Back-projection: 필터링된 투영 데이터를 공간 좌표계 $(x, y)$ 에 back-projection 합니다. 즉, 각 투영을 원래 각도 $\theta$ 에 맞춰 $(x, y)$ 평면에 되돌려서 적층하는 과정입니다.
4. Superposition(중첩 또는 합산): 마지막으로, 모든 각도에서 역투영된 데이터를 합산하여 최종 이미지를 만듭니다. 결과적으로 원래의 밀도 분포 $f(x, y)$ 가 복원됩니다.

FBP step1: Filtering

필터링 후 IFT 를 거치면 +, - 값을 가지는 profile 을 얻을 수 있습니다. 이 때 사용되는 필터는 ramp 라는 HPF 입니다.

FBP step2: Back-projection

FBP step3: Summation

과연 Ramp filter 가 최선일까?

대부분의 signal 은 저주파 성분이 많고 고주파 성분이 적습니다. 실제로는 위 그래프보다도 훨씬 극단적으로 저주파 성분이 많습니다. 반면에 noise 는 모든 주파수 영역대에서 골고루 분포한니다.

이러한 signal 에 왼쪽 그래프와 같은 HPF 를 곱하면 오른쪽 그래프 처럼 signal 분포가 변할 것입니다. 이 과정에서 noise 도 HPF 의 영향을 받기 때문에 고주파 영역의 noise 세기가 커져 signal 보다 noise 가 더 큰 영역도 존재할 수 있게 됩니다. 영상은 선명해지지만 noise 가 너무 많은 부분이 존재할 수 있게 됩니다.

이 문제를 해결하기 위해서 Butterworth filter 나 Gaussian filter 같은 LPF 를 HPF 와 같이 적용하면 BPF 같은 효과가 발생합니다. 저주파 영역대와 고주파 영역대가 모두 억제되면서 고주파 영역대의 noise 세기가 줄어들게됩니다. 다만 LPF 가 개입했으므로 해상도는 약간 떨어지게 됩니다.

FBP Demo

위 figure 는 Noisless target 과 그 필터링 결과입니다. 필터를 적용하지 않은 영상은 굉장히 흐릿해서 종양의 유무나 위치를 파악하기 어렵습니다. Ramp filter 만 적용한 영상은 굉장히 선명하지만 discretize 되어있기 때문에 테두리에 블럭 모양의 artifact 가 생겼습니다. Hanning filter 라는 BPF 를 적용한 영상은 Ramp filter 를 적용한 영상에 비해 선명도는 떨어지지만 필터를 적용하지 않은 영상보다는 훨씬 선명하며, artifact 문제도 꽤나 해결되었습니다.

BPF 의 진가는 노이즈가 많은 환경에서 발휘합니다. BPF 를 적용한 영상이 다른 것들에 비해 훨씬 종양을 파악하기 편합니다.

FBP의 장단점

장점

• 직관적이고 단순하다.
• 하나의 매개변수만으로 조정이 가능하다.
• 오랜 시간 동안 많이 사용되어 경험이 풍부하다.
• 1D 필터를 사용하기 때문에 계산이 빠르고 메모리 효율이 좋다.

단점

• 하나의 매개변수를 사용하기 때문에 구현이 너무 단순하다.
• Nonlinear noise 에 취약하며 streak artifact 를 유발한다.
• 다른 매개변수를 포함시킬 수 없다.

PET, SPECT 에서는 이러한 단점 때문에 FBP 가 더이상 사용되지 않으며, 여러가지 기술의 비교 대상으로만 사용됩니다. 하지만 CT 에서는 여전히 FBP 가 사용됩니다.

출처: 고려대학교 Medical Information Processing Lab.