오늘 갑자기 CT 영상이 어떻게 복원되는지가 궁금해졌습니다. 그래서 CT 영상 복원에 대한 가설을 간단하게 세우고 직접 실험을 해보았습니다.

Hypothesis

CT에 대해 아래와 같은 사실을 알고 있습니다.

• CT는 도넛 모양의 구조물이 회전하며 그 구멍 안의 물체를 촬영한다.
• CT는 한 번에 한 단면만을 촬영한다.
• CT는 X-Ray를 이용한다.
• ray가 물체를 통과할 때 그 세기가 지수적으로 감소한다.

거기에 다음과 같은 추가적인 추론을 할 수 있습니다.

• CT는 빔을 한 쪽에서 방사하고 다른 쪽에서 검출한다.
• 만약 빔이 한 가닥이라면 이는 속이 찬 원통과 속이 빈 원통을 구분할 수 없다. 그러므로 빔은 여러 가닥이다.

Modeling

이를 바탕으로 먼저 다음과 같이 균질한 물체에서 X-ray의 감쇄를 모델링했습니다.

$$I = I_0 e^{-\mu x}$$

이때 $I_0$는 물체를 통과하기 전의 빛의 세기, $I$는 물체를 통과한 후의 빛의 세기, $\mu$는 물체의 투과계수, $x$는 물체의 두께입니다.

계산의 편의를 위하여 $I_0$를 1로 두면 아래와 같이 변형할 수 있습니다.

$$I = e^{-\mu x}$$

이로부터 X-ray가 물체를 통과하면서

• 미소 면적 (그러므로 균질하다고 가정할 수 있는) $\Delta S_0$, $\Delta S_1$, $\Delta S_2$, ... 를 지날 때,
• 각 미소 면적의 감쇄율을 $\mu_0$, $\mu_1$, $\mu_2$, ... 라고 하고
• ray가 각 미소 면적을 지나는 길이를 $\Delta x_0$, $\Delta x_1$, $\Delta x_2$, ... 라고 하면,

물체를 통과한 후의 빛의 세기 $I$는 다음과 같이 표현할 수 있습니다.

$$I = e^{-\mu_0 \Delta x_0} e^{-\mu_1 \Delta x_1} e^{-\mu_2 \Delta x_2} ... \\ = e^{-\sum_{i=0}^n \mu_i \Delta x_i}$$

이는 계산하기 어려우므로 양변에 로그를 취하면 다음과 같이 변형할 수 있습니다.

$$\ln I =- \sum\_{i=0}^n \mu_i \Delta x_i$$

여기서 $n$과 $\mu_i$는 ray마다 다른 값을 가지므로 여러 ray들을 한 번에 행렬로 표현할 수 없습니다. 이를 해결하기 위해 다음과 같은 방법을 사용합니다.

• 측정할 어떤 단면을 미소 면적으로 쪼갠 후 (쪼개는 방법은 적당한 테셀레이션이면 충분합니다. 사각형 그리드 형태로 쪼개는 것이 계산이 편할 것입니다.)
• 각 미소 면적에 인덱스를 할당합니다. (인덱스를 할당하는 순서는 상관없습니다.)
• 선택 행렬 $W$를 정의합니다. $W_{i,j}$는 $i$번째 ray가 $j$번째 미소 면적을 지나는 거리입니다. 예컨대 각 변의 길이가 $d$인 정사각형 그리드에 대하여 이 값은 $0 \leq W_{i,j} \leq \sqrt{2}d$을 가질 수 있습니다.

그러면 위 수식은 다음과 같은 형태로 표현됩니다.

$$\ln I*i = - \sum*{j=0}^n \mu*j W*{i,j}$$

이는 다음과 같이 행렬 연산으로 표현할 수 있습니다.

$$\ln I = -W\mu$$

이때 $\mu$는 각 미소 단면의 투과계수를 원소로 가지는 벡터이고, $W$는 선택 행렬입니다. $I$는 측정된 각 ray의 빛의 세기를 원소로 가지는 벡터입니다.

• $I$는 실제로 측정된 알려진 값입니다.
• $W$는 계산에 의해 얻어지는 알려진 값입니다.
• $\mu$는 알고자 하는 값입니다.

단면을 미소 면적으로 나눈 개수를 n, ray의 개수를 m이라고 할 때, 각 행렬은 다음과 같은 크기를 가집니다.

$$\mu = (n, 1) \\ W = (m, n) \\ I = (m, 1)$$

이때 일반적으로 $n \neq m$이므로 $W$는 square matrix가 아닙니다. 그러므로 이 행렬 연산은 exact solution을 가지지 않으며, pseudo-inverse를 이용하여 해를 구해야 합니다. 그러면 아래와 같은 형태로 표현할 수 있습니다.

$$\mu = -W^+ \ln I$$

계산의 편의를 위해 $D = -\ln I$라고 하면

$$W^+ D= \mu\\(W^+ = (W^T W)^{-1} W^T)$$

Implementation

• 미소 면적은 64x64, 각 변의 길이가 1인 정사각형 그리드로 간주합니다.
• 계산상의 편의를 위해 $i$번째 ray와 $j$번째 미소 면적 사이의 거리를 $d_{i,j}$라고 할 때, $W_{i,j} = max(0, 1 - \frac{d_{i,j}}{d})$로 근사했습니다. 그리드 길이를 1로 두었으므로 이는 $W_{i,j} = max(0, 1 - d_{i,j})$입니다.
• CT를 묘사하기 위해 m개의 평행한 ray들이 180도 회전하며 촬영한다고 가정했습니다. 이때 한 스텝은 $\pi/256$ rad 입니다.
• $W^+$를 직접 계산하지 않고 곧바로 행렬 방정식 $W^+ D= \mu$를 푸는 것이 조금 더 빠르므로 그렇게 구현했습니다.
• $\mu$ 행렬은 원래는 $[0,1]$의 값을 가져야 하나 오차로 인해 범위를 넘는 값이 나옵니다. 이를 방지하기 위해 $\mu$ 행렬의 원소가 $[0,1]$의 값을 가지도록 clip합니다.

코드는 Python, Numpy, OpenCV를 사용하여 구현했으며 아래 리포지토리에 있습니다.

https://github.com/unknownpgr/ct-reconstruction

• config.py는 해상도 관련 설정을 담고 있습니다.
• calculate.py는 앞서 설명한 행렬 연산을 수행합니다.
• visualize.py는 결과를 시각화합니다.

Result

• 원본 이미지
• 재구성된 이미지

원본 이미지와 대단히 흡사한 이미지가 얻어진 것을 확인할 수 있습니다.

Discussion

위 과정에서 가장 오랜 시간이 걸리는 작업은 pseudo-inverse를 구하는 것입니다. Pseduo-inverse 연산의 computational time complexity는 $O(m^2n)$인데 ($m$은 ray의 개수, $n$은 미소 면적의 개수) $W$ 행렬의 크기는 $(23000, 4096)$ 정도로 대단히 크기 때문입니다. 그러나 ray를 방사하는 방법만 정해지면 $W$ 행렬은 변하지 않습니다. 따라서 ray를 방사하는 방법을 미리 정해두고 $W$ 행렬 및 그 pseudo-inverse를 미리 계산해두면 이후에는 시간을 크게 절약할 수 있을 것입니다.