2020년 11월 28일 토요일

라돈 변환(Radon Transform)

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


[그림 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)

댓글 6개 :

  1. 답글
    1. 세상을 바꾼 변환이라 많이 어려워요. 다만, 여기에 소개한 내용은 라돈 논문의 도입부입니다. 너무 양이 많아서 수학적으로 더 세밀한 부분은 생략했어요.

      삭제
  2. 역라돈변환이 적분인이유도 저기에 설명되어있는건가용

    답글삭제
    답글
    1. 네. 참고문헌 [1]에 잘 증명되어 있어요.

      삭제
  3. radon transform 설명에 mr 영상이라니..

    답글삭제
    답글
    1. 라돈 변환은 제가 감동받은 적분 변환 중의 하나입니다. 난공불락 같은 완결된 푸리에 급수를 원통 좌표계로 확장한 신기한 방식이 라돈 변환입니다.

      삭제

욕설이나 스팸글은 삭제될 수 있습니다. [전파거북이]는 선플운동의 아름다운 인터넷을 지지합니다.