라돈 변환(Radon Transform)

[경고] 아래 글을 읽지 않고 "라돈 변환"을 보면 바보로 느껴질 수 있습니다.

1. 푸리에 급수의 시작

2. 형보다 나은 아우: 푸리에 변환

3. 한켈 변환

4. 아벨의 적분 방정식

[그림 1] 라돈 변환과 컴퓨터 단층 촬영(출처: wikipedia.org)

푸리에 급수(Fourier series)와 푸리에 변환(Fourier transform)은 현대 수학을 만든 시초라 할 수 있다. 푸리에 변환에 버금가는 대단한 기여를 한 적분 변환중 하나는 라돈 변환(Radon transform)이다. 라돈 변환은 수학자 라돈Johann Radon(1887–1956)이 1917년라돈 30세, 일제 식민지 시절에 제안했다[1]. 사영(projection)을 이용해 원래 물체의 영상을 복원하는 라돈 변환의 대표적 응용은 [그림 1]에 소개한 컴퓨터 단층 촬영(computed tomography, CT)이다. CT 기술은 EMI의 공학자인 하운스필드Godfrey Hounsfield(1919–2004)에 의해 1972년하운스필드 53세, 박정희 정부 시절에 최초로 상용화되었다.

[그림 2] 미국을 방문한 비틀즈(출처: wikipedia.org)

음반 회사로 유명한 EMI는 비틀즈(The Beatles) 때문에 엄청난 매출과 수입을 올렸다. 비틀즈가 번 돈은 EMI가 만들던 CT 기술에도 투입되어서 인류를 위한 의료 기술 개발에 훌륭하게 쓰였다. 비틀즈의 노래를 사랑한 팬들은 알게 모르게 CT 기술에 대한 간접 지원을 한 것이었다. 돈을 매개로 한 예술과 과학의 놀라운 만남이다.

[그림 3] CT 주사기(scanner)의 내부 구조(출처: wikipedia.org)

[그림 3]은 전형적인 CT 주사기(scanner)의 내부 구조를 보여준다. X선관(X-ray tube) T에서 나온 X선 X는 공동[빈 구멍] 위치에 놓일 물체에 일부 흡수되면서 X선 검출기(X-ray detector) D까지 거의 직선으로 진행한다. 더 세부적으로 보면, X선은 직진성이 매우 높아서 송신기에서 수신기까지 가는 경로는 거의 직선으로 가정할 수 있고, 물체에 의해 흡수되기 때문에 물체의 종류에 따라 투과율이 달라진다. 예를 들어 전기장의 거리별 감쇠(attenuation)는 $|E(z)|$ = $|E_0| e^{-\alpha z}$로 표현되므로, [그림 3, 4]처럼 X선관에서 발생한 X선이 물체를 거쳐 검출기에 수신되는 복사 선속(radiant flux) 혹은 복사 전력(radiant power) $I_\theta(s)$는 다음처럼 공식화한다.

                  (1)

여기서 $I_0$는 X선관의 복사 선속, $\alpha(s, t; \theta)$는 물체를 통과할 때 생기는 위치별 감쇠 상수, $p_\theta(s)$는 사영 함수(projection function)이다. 식 (1)을 정리해서 사영 함수를 물체의 감쇠 상수로 표기할 수도 있다.

                  (2)

여기서 $I_\theta(s)/I_0$는 X선 검출기의 수신 전력 비율이다. 따라서 라돈 변환의 핵심인 사영 함수는 복사 선속에 대한 로그 함수 비율과 등가이다. 즉 라돈이 제안한 라돈 변환은 [그림 3]에 보여준 물체에 대한 X선(X-ray) 투과와 밀접히 연결된다.

[그림 4] 라돈 변환의 좌표계(출처: wikipedia.org)

(1) 물체(object)

(2) 광선의 시작선(starting line of rays)
(3) 광선의 종료선(ending line of rays)

(4) 투과 광선(transmission ray)

(5) 자료 원(datum circle)

(6) 원점(origin)

(7) 사영 함수(projection function) $p_\theta(s)$

[그림 4]는 라돈 변환을 정의할 때 사용하는 좌표계를 표현한다. [그림 3, 4]를 비교하면 수학적인 라돈 변환을 물리적 실재로 구현한 결과가 CT 기술임을 잘 알 수 있다. [그림 4]에서 주어진[혹은 고정된] $s, \theta$에 대해 점 $(x, y)$를 새로운 좌표축 $t$로 표현하면 다음과 같다.

                  (3)

여기서 $s$ = $x \cos \theta + y \sin \theta$, $t$ = $-x \sin \theta + y \cos \theta$, $x^2 + y^2$ = $s^2 + t^2$이다. 식 (3)을 간편하게 2차원 회전 행렬(rotation matrix)처럼 쓸 수도 있다.

                  (4)

여기서 점 $(s, t)$를 위한 좌표축 $s, t$는 기준 좌표축 $x, y$를 $\theta$만큼 회전시켜 얻는다. 시작선에서 나와서 종료선까지 적분하는 사영 함수(projection function) $P_\theta(s)$는 다음과 같은 선 적분(line integral)으로 정의한다.

                  (5)

여기서 $f(x, y)$는 고정점 $z$ = $z_0$에 대한 물체의 2차원 절편(slice)이며, $(x, y)$는 식 (3)의 정의를 사용한다. 사영 함수 $p_\theta(s)$는 공동도(空洞圖) 혹은 사이노그램(sinogram)이라고도 한다. 공동도는 구멍(sinus, opening)에 대한 그림(gram)이란 뜻이다. 변수 $t$가 변하는 범위를 전체 실수로 확장하면 식 (5)는 라돈 변환이 된다.

                  (6)

각도를 $\theta + \pi$로 바꾸면, $p_{\theta + \pi}(s)$ = $p_\theta(-s)$ 혹은 $p_\theta(s)$ = $p_{\theta + \pi}(-s)$가 항상 성립한다.

[그림 5] 푸리에 절편 정리(출처: wikipedia.org)

사영 함수는 푸리에 변환 관점으로 관찰해야 숨겨진 내면을 볼 수 있다. 푸리에 절편 정리(Fourier slice theorem) 혹은 사영-절편 정리(projection-slice theorem)에 따라 사영 함수의 푸리에 변환을 구한다. 2차원에 대한 푸리에 절편 정리는 다음과 같다. 

                  (7)

식 (7)에 제시한 푸리에 절편 정리에 의해, 사영 함수 $p(x)$의 푸리에 변환 $\mathfrak{F}[p(x)]$는 차원 하나가 축소된 푸리에 변환 $P(\xi)$ = $F(\xi, 0)$을 만든다. 이 개념을 $N$차원 공간으로 일반화해서 식 (7)을 $N$차원 푸리에 절편 정리로 확장한다.

                  (8)

여기서 $p(\cdot)$는 좌표 성분 $x_{N-m+1}, \cdots, x_N$에 대한 공간 영역의 사영(projection), $\mathfrak{F}_{N-m}[\cdot]$는 $x_1, \cdots, x_{N-m}$에 대한 푸리에 변환, $F(\cdot)$는 $\xi_{N-m+1}, \cdots, \xi_N$ 성분을 $0$으로 바꾼 파수 영역의 절편(slice)이다.

식 (7)에 의해 공간 영역(spatial domain)에서 회전 각도를 $\theta$만큼 돌리면서 구한 사영 함수 $p_\theta(s)$의 푸리에 변환 $P_\theta (\kappa)$는 파수 영역(spectral domain)에서도 각도 $\theta$만큼 회전한다. 따라서 다음 관계식을 보간(interpolation)해서 $f(x, y)$의 2차원 푸리에 변환 $F(\xi, \eta)$를 근사한다.

                  (9)

2차원 푸리에 절편 정리인 식 (9)를 이용해서 파수 영역의 2차원 푸리에 변환 $F(\xi, \eta)$를 쉽게 구할 수 있다고 오해할 수 있다. 하지만 파수 영역의 절편 함수(slice function) $P_\theta (\kappa)$는 원점 $(\xi, \eta)$ = $(0, 0)$를 지나는 직선이라서, 원점 부근의 표본[낮은 파수 성분]이 과대 평가되고 원점을 벗어난 표본[높은 파수 성분]은 과소 평가되는 심각한 문제가 있다.

[그림 6] 복소 영역에서 절편 함수의 분포(출처: wikipedia.org)

파수 영역의 절편 함수[파란색 선]가 가진 한계를 [그림 6]이 분명히 보여준다. 원점에는 직선이 잘 모이지만, 반지름이 커지면 직선간의 간격이 넓어져서 문제가 된다. 그래서 [그림 6]과 같은 파수 영역의 절편 함수를 적절히 보간(interpolation)해서 근사적으로 절편 표본간의 간격을 비슷하게 만든다. 그 다음에 보간한 $F(\xi, \eta)$의 푸리에 역변환을 취해서 물체의 2차원 영상 $f(x, y)$를 다음처럼 근사적으로 복원한다.

                  (10)

여기서 $\operatorname{Int}[{\bf A}]$는 집합(set) $\bf A$의 원소로 만드는 보간을 의미한다.

식 (10)은 근사가 많이 들어간 영상 복원법이라서, 수학적으로 더 엄밀하고 세련되게 이론을 전개할 필요가 있다. 이를 위해 원통 좌표계에서 사용하는 푸리에 변환인 한켈 변환(Hankel transform)을 식 (10)에 적용해서 정확한 라돈 역변환(inverse Radon transform)을 정의한다. 모든 각도 $\theta$에 대해 식 (9)를 계산한 후 복원한 물체의 영상 $f(x, y)$는 다음과 같다[3].

                  (11)

사영 함수의 푸리에 변환 $P_\theta (\kappa)$는 $\theta$에 대해 주기적이므로 푸리에 급수로 전개한다. 이 결과를 식 (11)에 넣고 한켈 변환처럼 정리해서 라돈 역변환을 얻는다.

                  (12)

여기서 $x$ = $\rho \cos \phi$, $y$ = $\rho \sin \phi$, $J_n (\cdot)$는 제$n$차 제1종 베셀 함수(Bessel function of the first kind)이다. 파수 영역에서 정의된 $P_n(\kappa)$는 $\theta$에 대한 평균값으로 정확히 구한다.

                  (13)

한켈 변환으로 라돈 역변환을 유도할 수 있지만, 권위 있는 원류는 라돈이 제안한 방법[1]이다. 방사선 원소인 라돈(Radon, Rn) 혹은 라듐(radium, Ra)과 이름이 비슷해서 위험하다는 오해를 받을 때도 있지만, 라돈은 인류에 기여한 아름다운 수학자이다. 라돈이 제안한 라돈 역변환은 사영 함수 $p_\theta(s)$로부터 시작한다. 디랙 델타 함수(Dirac delta function)를 도입해서 사영 함수를 2차원 공간에서 다시 표현한다.

                  (14)

[그림 4]에서 보여준 좌표축의 회전 각도 $\theta$에 대한 평균 사영 함수 $\bar p(s)$는 다음과 같다.

                  (15)

데카르트 좌표계(Cartesian coordinate system)를 원통 좌표계(circular cylindrical coordinate system)로 바꾸어서 디랙 델타 함수의 적분을 구한다.

                  (16)

여기서 $x$ = $\rho \cos \phi$, $y$ = $\rho \sin \phi$, $\rho > |s|$, $\theta_1$과 $\theta_2$는 $\delta(\cdot)$의 입력 변수에 대한 근이며 $\rho  \cos (\theta_{1,2} - \phi)$ = $s$가 성립한다. 식 (16)을 식 (15)에 넣고 $f(x, y)$의 평균값 $\bar f(\rho)$가 피적분 함수에 나타나도록 정리한다.

                  (17)

식 (17)은 아벨의 적분 방정식(Abel's integral equation)이므로, 다음과 같은 해석적인 해법으로 쉽게 해결된다.

          (18)

여기서 $\alpha$ = $1/2$이다. 따라서 $\bar f(\rho)$는 다음처럼 공식화된다.

                  (19)

여기서 식 (18)의 오른쪽 식에 나온 극한은 $0$이라 가정한다. 결국 반지름 $\rho$를 $0$으로 보내면, 함수값 $f(0, 0)$를 사영 함수의 적분으로 정확히 표현할 수 있다.

                  (20)

[그림 7] 일반적인 라돈 역변환을 위한 좌표 변환

임의의 위치에서 $f(x, y)$를 완벽하게 복원하기 위해, [그림 7]처럼 원점 $(0, 0)$을 $(x_0, y_0)$으로 보내는 좌표 변환에 따라 평균값 $\bar f(\rho)$와 $\bar p(s)$를 다음과 같이 일반화한다.

                  (21)

여기서 $\bar f(0, 0; r)$ = $\bar f(r)$, $\bar p(0, 0; q)$ = $\bar p(q)$, $q$는 $(x_0, y_0)$를 기준으로 양수나 음수가 된다. 식 (15), (17)에 식 (21)을 적용해서 식 (19)를 변형한다.

                  (22)

최종적으로 반지름 $r$이 $0$으로 가는 극한에 의해 식 (22)는 물체의 2차원 영상 $f(x, y)$가 된다.

                  (23)

식 (23)의 적분에는 함수 미분소가 있어서 다소 복잡하므로, 부분 적분을 적용해서 다음처럼 단순화한다.

                  (24)

[라돈 역변환의 조건]

라돈 역변환이 성립하기 위한 조건은 다음과 같다.

• 물체의 영상 $f(x, y)$는 연속이다.
• 식 (17)의 성립을 위해 다음 이중 적분은 수렴한다: 
  
  
• 식 (22)를 만족하기 위해 임의의 점 $(x, y)$에서 다음 극한이 성립한다: $\lim_{r \to \infty} \bar f(x, y; r)$ = $0$ 

[참고문헌]

[1] J. Radon, "Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten (On the determination of functions by their integral values along certain manifolds)," Berichte über die Verhandlungen der Königlich-Sächsischen Akademie der Wissenschaften zu Leipzig (Reports on the Proceedings of the Royal Saxonian Academy of Sciences in Leipzig), vol. 69, pp. 262–277, Apr. 1917. (In German, 방문일 2020-11-28)

[2] 박미경, 이의택, "Computed Tomography에 의한 절편 영상의 추정", 전자통신동향분석, 제12권, 제6호, pp. 50–57, 1997년 12월.

[3] E. Won, "Derivation of the inverse Radon transformation," Eunil Won's Home Page, Nov. 2007. (방문일 2020-12-05)

[4] 한화택, 미적분의 쓸모: 미래를 예측하는 새로운 언어, 더퀘스트, 2021.