Iterative Attenuation Correction and Image Reconstruction Using Time-Of-Flight Positron Emission Tomography

양전자방출단층촬영기의 비행시간정보를 이용한 반복적 감쇠보정 및 영상재구성

Abstract

In this paper, an iterative method is proposed to perform attenuation correction and image reconstruction simultaneously for positron emission tomography, by using the time-of-flight information. Numerical simulation results are presented to demonstrate an improved performance of the proposed method in attenuation correction and image reconstruction.

1. 서 론

양전자방출단층촬영기(Positron Emission Tomography, PET)는 인체에 주사 또는 흡입된 방사성의 약품의 체내분포를 영상의 형태로 표현하는 의료기기로, 악성종양, 뇌신경계 및 심장 질환의 진단에 중요한 정보를 제공한다[1]. PET는 방사성동위원소에서 양전자가 방출직후에 전자와 결합하여 쌍소멸하며 거의 180도 반대방향으로 두 개의 감마선을 방출하는 현상을 이용한다. 구체적으로, 두 개의 검출기 b1,b2가 각각 t1,t2 시간(이를 비행시간(Time-Of-Flight, TOF)정보라 한다)에 동일한 양전자방출에서 기인한 감마선을 검출하였다면, 그 두 개의 검출기를 연결하는 어느 중간부분에 양전자가 하나 존재했다는 것을 의미하게 된다. 이러한 관측정보 (b1의 위치, t1, b2의 위치, t2)가 나열된 데이터를 list mode 데이터라고 하는데, 여기서, 주어진 검출기 쌍 b=(b1,b2)에 대응하는 list mode 데이터의 수 nb(이를 Raw Data라 하자)는 두 개의 검출기 b1,b2를 잇는 반응선상에 분포한 양전자의 밀도를 의미하게 된다. PET는 검출기를 여러 곳에 배치하여 얻은 이러한 데이터를 이용하여 양전자의 분포를 영상의 형태로 재구성한다[2].

비행시간정보를 이용하여 PET의 영상재구성을 향상시키는 것은 PET의 개발초기부터 논의되었으나, 당시의 PET의 낮은 시간분해능하에서는 어려운 일이었다. 현재 각급 병원에 있는 PET 기기는 약 3ns의 시간분해능을 가지고 있어 9cm이내로 양전자의 위치를 추적할 수 있다. 수년 내에 병원에 활발히 보급될 것으로 예상되는 PET 기기는 100∼300 ps의 시간분해능을 가져 하나의 list mode 데이터를 가지고 3∼9mm이내로 양전자 위치를 추적할 수 있을 것으로 예상된다. 이러한 PET의 시간분해능의 향상은 비행시간정보를 이용한 PET의 영상재구성 향상의 가능성을 높여주고 있다[1].

PET는 관측데이터를 영상재구성에 사용하기에 앞서 우연계수보정, 산란보정, 감쇠보정 등의 데이터의 정확성을 높이는 전처리 과정이 필요하다[2]. 모든 보정이 중요하지만, 본 논문에서는 감쇠보정만 고려하고자 한다. 감쇠는 감마선이 검출기에 의해 측정되지 못하고 인체 조직 내에 흡수되는 현상인데, 보정하는 방법의 하나는 X-ray CT(Computerized Tomography) 또는 MRI(Magnetic Resonance Imaging) 기기를 이용하여 감쇠를 일으키는 인체장기에 대한 해부학적 영상을 동시에 측정하는 것이다. 그러나, 부가 장비의 사용은 PET기기의 가격을 높이는 단점이 있다[1].

논문 [3]의 연구결과는 PET의 TOF 정보는 영상 재구성의 향상과 함께 감쇠보정에 쓰일 수 있는 가능성을 보여주고 있다. 그러나, 제안한 TOF PET 데이터를 이용한 감쇠보정 방법은 실제 데이터의 다양한 물리적 특성을 고려하지 않고 수학적 함수의 해석적 특성을 고려하여, PET의 시간분해능의 향상에 의한 장점을 충분히 활용할 수 없는 단점이 있다.

이 논문에서는 PET의 TOF 정보를 이용한 기존의 감쇠보정 방법[3]의 문제점을 살펴보고 그 문제점을 극복하면서 영상재구성을 동시에 수행하는 반복법을 제안하고자 한다.

 

2. 이 론

PET 영상재구성은 list mode 데이터를 분석하는 방법에 따라 다음 2가지 방법으로 나눌 수 있다.

A) Raw Data nb를 그대로 사용하는 경우

이 경우에는 관측된 list mode 데이터를 Raw Data nb로 재배열하고, 재구성하고자 하는 방사성의 약품의 체내분포영상 fυ(여기서, υ는 화소 인덱스로 fυ는 화소 υ에서의 방사성의약품의 체내분포밀도를 의미)와의 nb와의 선형관계

를 정의하는 시스템 행렬 Pb,υ를 미리 계산한 후에, 반복적 영상재구성을 적용하여 fυ를 추정한다.

이 방법은 관측된 데이터와 재구성하려는 영상과의 관계를 정확히 구현할 수 있다는 장점이 있으나, 관측된 데이터와 재구성하고하는 영상과의 관계를 결정하는 PET의 설계적 특징과 물리적 현상을 많이 고려하면 할수록 영상재구성에 더 많은 계산이 필요로 한다는 단점이 있다[2].

B) Raw Data nb를 재배열한 데이터 lc를 사용하는 경우

이 경우에는 관측된 list mode 데이터(b1,t1,b2,t2)를 검출기 쌍 (b1,b2)를 잇는 반응선을 기준으로 재배열하여 얻은 데이터 lc를 사용한다. 이 때, 반응선은 각도, 원점으로부터의 거리, 높이, 횡단면에 대한 기울어짐을 기준으로 분류하여 재배열한다. 재배열된 데이터 lc와 재구성하고자 하는 fυ와의 선형관계

를 이용하여 방사성의약품의 체내분포영상 fυ를 직접적(Filtered Backprojection, FBP) 또는 반복적으로 추정한다. 여기서, 시스템 행렬 Kc,υ는 미리 계산하여 저장한 것을 사용하기보다는 필요할 때 마다 계산하여 사용하는 방식을 이용한다[2].

이 논문에서 제안하고자 하는 반복적 감쇠보정 및 영상재구성 알고리즘은 (2)의 방법에 따라 재배열된 데이터에 적용되어 2차원 또는 3차원 영상재구성에 사용되는 것이다. 그러나, 설명의 편이를 위해 지금부터의 설명은 2차원 영상재구성에 한정하고, 이산화된 인덱스 c,υ 대신 연속변수 (θ,u)(여기서, θ)는 각도변수를 u는 원점으로부터의 거리를 나타내는 변수이다), (x,y)를 사용하고자 한다. 이 논문에서 사용하고자 하는 수학적 기호들과 그 의미를 Table 1에 정리하였다.

Table 1.List of mathematical notations and their meanings

Table 1의 일부 수학적 기호들은 수식을 통해 설명하는 것이 더 분명하다. 예를 들어,

로 설명할 수 있다[4]. 한편, TOF 정확도를 가리키는 점확산함수 w가

의 형태의 가우스 함수라고 가정한다. 이 조건을 (4)에 적용하면

을 얻을 수 있음은 자명하다.

실제 관측데이터 m(θ,u,t)는 보통 관측 잡음, 감쇠, 산란, 재배열과정에서의 오차 등 다양한 원인에 의해 왜곡된다[2]. 그러나, 만약, 왜곡이 감쇠에 의해서만 이루어지는 경우에는 실제 관측데이터 m(θ,u,t)와 이상적인 데이터 y(θ,u,t)의 관계는

로 기술할 수 있다[4]. 동일한 이유로 M(θ,u)와 Y(θ,u)의 관계는

로 기술할 수 있다.

이 논문에서 제안하고자 하는 반복적 감쇠보정 방법은 반복법은 물론, FBP와 같은 직접법에도 사용될 수 있으나, 설명의 편이를 위해 FBP를 이용하여 제안하는 방법을 설명하고자 한다. FBP 영상재구성은 다음과 같은 수학적 정리에 근거를 두고 있다[4].

여기서, 푸리에 변환의 실제 계산은 FFT(Fast Fourier Transform)를 이용하여 계산한다[4].

PET의 TOF를 이용한 감쇠보정과 관련된 가장 대표적인 연구 [3]의 핵심내용은

의 공식을 이용하여 감쇠률 A(θ,u)의 편미분을 계산할 수 있다는 것이다. 여기서,

이다.

 

3. 제안한 방법

논문 [3]에서는 (12), (13), (14)을 이용하여 감쇠보정할 수 있음을 수학적으로 증명하였다. 그러나, 이 방법은 잡음에 오염된 데이터에 적용하는 경우에는 수치미분의 불안정성 때문에 만족 할 만한 결과를 얻기 힘들다. 이러한 문제점을 극복하기 위해 본 논문에서는 미분과정을 푸리에 변환을 이용하여 계산하는데, 특정주파수이하의 정보만 절단하여 사용하는 방법을 제안한다. 구체적으로,

와 같이 수치미분은 푸리에 해석으로 대체하고 잡음에 의한 영향을 주파수영역 |ξ|>Ω에 대한 정보를 제거함으로써 극복하고자 한다. 동일한 방법으로, 편미분 과 2차 편미분 은

의 방법으로 계산한다. 변수 θ에 대한 주기함수에 대한 편미분 은

와 같이 절단 푸리에 급수전개를 이용하여 계산한다. 여기서, 2π 주기함수 m(θ,u,t)가

으로 표현됨을 이용하였다.

식 (15), (16), (17), (18)의 결과를 (12), (13), (14)에 대입하여 을 계산한 후에 다음의 최소자승법을 이용하여 A를 계산한다.

여기서, (20)의 은 (12), (13), (14)에 의해 관측데이터로부터 계산된 결과이고. 미분연산자를 표현하는 행렬 P는 중앙차분법에 의하여 이산화 하였다[5]. 식 (20)의 정규화항 λ∥A∥2은 최소자승법이 잡음에 강건하게 만들기 위해 제안되었다. 식 (21)의 제한조건은 감쇠계수 µ(x,y)에 대한 일관성 조건[4]을, 그리고, (22)는 감쇠계수 µ(x,y)가 x2 + y2 < r02에서만 존재한다는 가정을 이용하였다.

이 논문의 시뮬레이션에서는 A는 (20)에서 유도되는 정규방정식(Normal Equation)에 공역기울기법(Conjugate Gradient Method, CGM)을 5회 적용하여 계산하였다[5,6,7].

지금까지 설명한 내용은 m(θ,u,t)에서 A(θ,u)를 추정하는 과정으로, 이것이 정확히 수행되었다면, (9)를 이용하여 감쇠보정된 데이터 Y(θ,u)를 A(θ,u)로부터 계산하고, 최종적으로 (10), (11)을 이용하여 영상 f(x,y)를 재구성할 수 있다. 그러나, 이 일련의 과정들은 주어진 관측데이터 m(θ,u,t)의 부정확성, 수치미분, 수치적분 과정의 불안정성 등에 의한 어려움을 내포하고 있다. 이러한 문제점들을 해결하기 위해, 이 논문에서는 감쇠보정과 영상재구성을 교차적으로 진행하는 다음의 반복법을 제안한다.

Iterative Attenuation Correction and Image Reconstruction

Starting with all-one image f0:

for n = 0, 1, 2,...,N

Compute

Compute An+1/2(θ,u) by using CGM on (20) with An(θ,u) as initial guess

Compute Yn+1(θ,u) = exp(A(θ,u))M(θ,u)

Compute fn+1(x,y) by using FBP on Yn+1(θ,u)

end

 

4. 실험 결과

Fig. 1은 시뮬레이션에 사용된 영상 f(x,y)로 그 크기는 201 × 201 이다. 영상 f(x,y)의 화소의 크기는 2 × 2 mm2 라고 가정하였다.

Fig. 1.Test images f(x,y) of size 201 × 201, (a) Hot Rod, (b) Cold Rod.

Fig. 2(a)는 시뮬레이션에 사용된 감쇠계수 µ(x,y)로 그 크기는 201 × 201 이고, Fig. 2(b)는 감쇠계수 µ(x,y)를 식 (3)에 적용하여 계산한 감쇠률 A(θ,u)를 보여주고 있다. Fig. 2(b)의 감쇠률 A(θ,u)의 크기는 200 × 201로 하였는데, 이는 각도 θ의 수를 200개(즉, 각도를 0.9° 씩 변화시켜가며 관측)로, 그리고, 원점으로부터의 거리 u의 수를 201개로 함을 의미한다.

Fig. 2.(a) Test attenuation coefficients µ(x,y) of size 201 × 201, (b) attenuation rate A(θ,u) of size 200 × 201 computed by (3) using µ(x,y).

TOF 정확도를 가리키는 점확산함수 w의 σ(식(6) 참조))가 5mm라고 가정하였다. 또한, 시간변수 t의 수가 101라고 가정하였고, 관측환경의 불확실성을 모사하기 위해 관측데이터 m(θ,u,t)의 1%를 표준편차로 갖는 정규잡음을 추가하였다. Fig. 3(a)는 Fig. 1(a)의 ‘Hot Rod’에 대한 관측데이터 M(θ,u)를, Fig. 3(b)는 Fig. 1(b)의 ‘Cold Rod’에 대한 관측데이터 M(θ,u)이다.

Fig. 3.Observed data M(θ,u) (of size 200 × 201) (a) of ‘Hot Rod’ in Fig. 1(a) and (b) of ‘Cold Rod’ in Fig. 1(b).

Fig. 4는 Fig. 3(a)의 관측데이터 M(θ,u)에 FBP를 적용하여 재구성한 영상이다. FBP 시뮬레이션에는 절단주파수(Cutoff Frequency)를 0.5로 한 Hann 필터를 적용하였다. Fig. 4(a)는 감쇠보정을 하지 않은 관측데이터 M(θ,u)에 FBP를 적용한 결과이고, Fig. 4(b)는 제안한 반복적 감쇠보정 및 영상재구성 알고리즘을 10회 적용한 결과이다. 원영상 Fig. 1(a)와의 비교를 통해 제안한 방법의 효과가 분명함을 알 수 있다.

Fig. 4.FBP Reconstruction using observed data M(θ,u) in Fig. 3(a). (a) reconstructed image without attenuation correction, (b) reconstructed image with the attenuation correction by the proposed method.

Fig. 5는 Fig. 3(b)의 관측데이터 M(θ,u)에 FBP를 적용하여 재구성한 영상이다. Fig. 5(a)는 감쇠보정을 하지 않은 관측데이터 M(θ,u)에 FBP를 적용한 결과이고, Fig. 5(b)는 제안한 반복적 감쇠보정 및 영상재구성 알고리즘을 10회 적용한 결과이다. 원영상 Fig. 1(b)와의 비교를 통해 제안한 방법의 효과가 분명함을 알 수 있다.

Fig. 5.FBP Reconstruction using observed data M(θ,u) in Fig. 3(b). (a) reconstructed image without attenuation correction, (b) reconstructed image with the attenuation correction by the proposed method.

 

5. 결 론

제안한 방법은 PET의 TOF 정보를 이용하여 감쇠보정과 영상재구성을 동시에 반복적으로 진행하는 것으로, 그 핵심은 정확히 추정된 감쇠정보가 향상된 영상재구성을 제공하듯이 정확한 영상재구성이 정확한 감쇠보정을 유도할 수 있다는 사실을 이용한 것이다. 제안한 방법의 가능성은 4장의 실험 결과를 통해 확인할 수 있었다.

TOF 정보를 이용한 영상재구성의 응용방향의 하나는 TOF 정보를 이용한 영상재구성의 향상이 관측에 필요한 방사능 피폭량의 감소를 유도하는 것이다[8,9]. 본 논문에서 제안한 방법이 별도의 CT 또는 MRI 측정 없이 TOF 정보를 이용한 감쇠보정과 저선량 방사능 피폭으로도 가능한 영상재구성 알고리즘을 개발에 기여하기를 기대한다.