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)

댓글 20개 :

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

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

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

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

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

      삭제
  4. 혹시 저랑 연락하면서 이것좀 알려주살수 있나요 이해가 잘 안되서 정말 궁금합니다 보셧으면 연락 주세요 010 7341 8946

    답글삭제
    답글
    1. 궁금한 점은 댓글로 남겨주세요.

      본문이 잘 이해되지 않으면, 푸리에 급수, 푸리에 변환, 적분 방정식 등의 지식이 부족한 게 원인일 수 있어요.
      상단에 있는 기존 글부터 공부해서 정복을 해야 합니다.
      라돈 변환의 더 쉬운 설명은 참고문헌 [4]를 참고할 수도 있어요.

      삭제
  5. 식 19번에 로 제곱과 s 제곱의 위치가 바뀐것 같은데요?

    답글삭제
    답글
    1. 그리고 추가로 가장 마지막 좌표 변환을 통해서 f(0,0)의 값 만을 알던 걸 갑자기 모든 위치의 값을 알 수 있게 할 수 있나요? 아직 세세한 내용을 이해를 못하겠어서 그냥 언뜻 보기엔 수학적 속임수(?) 같은 느낌이 들거든요

      삭제
    2. 1. 식 (19)는 문제 없습니다.

      2. [그림 4]처럼 방위각을 돌리면서 광선을 쏜 값인 $p_\theta(s)$를 알아야 라돈 역변환으로 내부 함수 $f(x, y)$를 구할 수 있어요. 이 과정이 바로 CT나 MRI에 직접 쓰입니다.

      삭제
    3. 지금보니 제가 약간 말투가 공격적이었던 것 같네요;; 그 점 사과드립니다.
      분명 유도과정에선 루트 안 값이 양수가 되기 위해서 적분이 로가 s보다 클 때 수행되었는데 식 19번은 그 위치가 바뀌어 있어서 음수가 아닌가...? 하는 생각이 들었습니다. 혹시 시간이 나신다면 이 점 설명해주시면 정말 감사하겠습니다.

      삭제
    4. 공격적이라고 느낀 적 없습니다 ^^

      식 (19)가 자연스럽지 않으면, 아벨의 적분 방정식인 식 (18)을 보세요. 식 (18)에서는 $x, t$를 바꾸었고요, 식 (19)는 $s, \rho$를 바꾸지 않았어요. 변수 선택은 우리의 자유입니다.

      삭제
    5. 혹시 이 내용은 어디서 공부하신 건가요? 그리고 만약 제가 가능하다면 이 내용을 인용하고 싶은데 어떻게 언급하면 좋을까요?

      삭제
    6. 1. 참고문헌 [1]에 라돈 논문의 번역본이 링크로 설정되어 있어요. 이 IEEE TMI 논문으로 공부했어요.

      2. 혹시 논문 같은 공식 문서를 쓰신다면, 이 블로그를 인용하지 마시고 IEEE TMI 논문을 쓰세요.

      삭제
    7. 친절하신 답변에 정말 감사드립니다. 마지막으로 하나만 더 여쭤볼까 하는데 제가 이 라돈변환 파트를 한 1년동안 분석하면서 이 페이지에서 설명되는 전반적인 내용은 이해가 됐습니다. 그런데 이제 IEEE 논문을 확인해보니까 전혀 새로운 내용이 전개되는 것 같더라구요. 그래서 제가 이해한 내용이 맞을지 확인해주실 수 있으실까요?
      1. 이 페이지에서 라돈 역변환 공식을 유도했다
      2. 논문에서는 5가지의 theorem을 설명하는데 4번째 까지는 이 페이지 내용의 검증이고 5페이지는 저도 잘 모르겠습니다;;
      (제가 이 페이지에서 정말 많은 것을 배울 수 있었습니다. 페이지 설립부터 답변까지 모든 과정과 노고에 감사드립니다.)

      삭제
    8. 3. 이 변환을 통해 방사선을 조사한 함수를 변환해서 물체의 방사선 흡수도(f(x,y))를 나타낼 수 있다.
      4. 하지만 이 실제로는 laminogram, 푸리에 변환을 이용해서 데이터 변조과정을 거쳐 blur이 최소화된 이미징을 실시한다.
      마지막 궁극의 질문 : 그렇다면 실제로 이 변환을 실시하지 않는 이유는 무엇일까?

      삭제
    9. 2번째 질문에 5페이지가 아니라 5번째 theorem을 말하는 것입니다. 감사합니다

      삭제
    10. 1. 논문에 라돈 역변화 공식이 있어요. 참고문헌 [1]의 정리 III이 식 (24)입니다.

      2. 본문에 라돈 논문 전부를 소개하고 있지 않아요. 내용이 쉬운 정리 III까지만 상세히 증명했어요.

      3. 라돈 변환이 CT의 기본 이론입니다. 식 (2)가 이를 나타냅니다.

      4. 식 (24)는 수학적으로 타당하지만(라돈이 유일성까지 증명), 정확하게 적분하기가 어려워요.
      그래서 실제 계산에서는 다양한 편법이 동원됩니다.

      그중 푸리에 변환은 FFT로 인해 빠른 계산이 가능합니다.
      어떤 적분 변환이든지 보간을 잘 해서 FFT를 쓸 수 있다면 만사형통입니다.
      예를 들어, 여기 라돈 변환이든 또 다른 한켈 변환이든 신호 처리할 때는 푸리에 변환이 효율적입니다.
      푸리에 변환을 쓰면 더 많은 자료를 더 빨리 계산할 수 있어요.

      삭제
    11. 여태 의문이었던게 명쾌하게 해답됐네요 감사합니다

      삭제

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