[경고] 아래 글을 읽지 않고 "라돈 변환"을 보면 바보로 느껴질 수 있습니다.
[그림 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$으로 바꾼 스펙트럼 영역(spectral domain) 혹은 파수 영역(wavenumber domain)의 절편(slice)이다.
식 (7)에 의해 공간 영역
(spatial domain)에서 회전 각도를 $\theta$만큼 돌리면서 구한 사영 함수
$p_\theta(s)$의 푸리에 변환 $P_\theta (\kappa)$는 파수 영역에서도 각도 $\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$
따라서 라돈 변환으로 물체의 영상 $f(x, y)$을 정상적으로 복원하기 위해, $f(x, y)$가 라돈 역변환의 조건을 만족하는지 미리 확인하거나 이를 가정해야 한다.
1. 허프 변환(Hough transform)
[그림 1.1] 허프 변환(Hough transform)에 사용하는 직선의 매개변수 $(r, \theta)$(출처: wikipedia.org)
[그림 4, 7]이 채택하는 라돈 변환의 매개변수 $(s, \theta)$를 재미있게 활용하는 예시가 1959년
허프 42세, 이승만 정부 시절에 나온
허프 변환(Hough transform)이다[5]. 허프
Paul Hough(1917–1989)가 거품 상자
(bubble chamber)의 직선 궤적 추적용으로 개발한 허프 변환은 복잡한 영상에서 직선과 같은 특정 도형을 효과적으로 도출하는 특징 추출
(feature extraction) 기법이다. 원래 허프 변환은 디지털 영상 처리
(digital image processing)에 주로 사용되다가 요즘은 기만체
(decoy)를 판별하는
레이다 신호 처리
(radar signal processing)에도 활용된다[6].
라돈 변환에 쓰이는 [그림 7]의 $(s, \theta)$ 대신, 허프 변환은 매개변수를 [그림 1.1]과 같은 $(r, \theta)$로 선택한다. 보통 일차 함수를 $y$ = $mx+b$로 표기하기 때문에, 기울기 $m$과 $y$절편 $b$로 된 매개변수 $(m, b)$를 선택하면 더 편리하고 직관적일 것 같다. 왜 [그림 1.1]과 같은 매개변수 $(r, \theta)$를 쓸까? 조금만 생각하면 쉽게 답할 수 있다. 초등 수준에서 매개변수 $(m, b)$는 간단해서 좋지만, 직선이 $y$축에 붙는 경우 기울기 $m$이 무한대로 갈 수 있어서 큰 문제가 된다. 이를 해결하기 위해 [그림 1.1]처럼 원점에서 직선까지 거리인 $r$ 및 $x$축과
파란색 법선의 각도 $\theta$로 직선의 방정식을 표현한다. 매개변수 $(r, \theta)$로 된 직선의 방정식을 만들 때는
헤세 표준형(Hesse normal form)을 이용한다. 직선 위의 점은 $\bar r$ = $(x, y)$, 법선 벡터는 $\bar n$ = $(r\cos \theta, r\sin \theta)$이다. 그러면 헤세 표준형에 따라 $(\bar r - \bar n) \cdot \bar n$ = $0$이 나와서 직선의 방정식이 $(r, \theta)$로 표현된다.
(1.1)
만약 $r$ = $0$이면, 이 직선은 원점을 지나서 $\hat n \cdot \bar r$을 만족한다. 여기서 $\hat n$ = $(\cos \theta, \sin \theta)$이다. 그래서 식 (1.1)은 $r$ = $0$에서도 성립함으로써 $x \cos \theta + y \sin \theta$ = $0$이 잘 나온다.
(a) 세 점을 지나는 직선의 추정
(b) 영상에서 두 직선의 결정
[그림 1.2] 허프 변환을 활용한 직선의 특징 추출의 예시(출처: wikipedia.org)
[그림 1.2]는 허프 변환을 이용해서 직선의 특징을 추출하는 간단한 예시를 보여준다. [그림 1.2(a)]와 같이 세 점이 주어진 경우에 가장 적합한 직선의 방정식을 결정한다. 각 세 점 $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$을 식 (1.1)에 대입하고 각도 $\theta$를 15˚ 간격으로 바꾸면서 $r$을 계산한다. 각 세 점에서 $(r, \theta)$가 거의 변하지 않는 $\theta$ = $60^\circ$인 파란색 직선이 우리가 원하는 답이다. 이번에는 [그림 1.2(b)]의 좌측 영상에서 발생한 직선을 인식한다. 검정색으로 된 직선의 점 $(x, y)$를 선택하고 식 (1.1)에서 독립 변수 $\theta$를 바꾸면서 종속 변수 $r$을 찾는다. 이런 방식으로 계산한 $(r, \theta)$를 [그림 1.2(b)]의 우측 영상에 흑백 화소(pixel)로 찍는다. 이때 흑백 영상을 만들기 위해 도수 분포(度數分布, frequency distribution) 개념을 차용한다. 각 점에서 찾은 $(r, \theta)$의 화소 도수가 늘어날수록 더 밝은 색을 할당한다. 결과적으로 [그림 1.2(b)]의 우측 영상에서 도수가 많아서 굉장히 밝은 두 점이 나타난다. 따라서 좌측 영상에 나온 직선은 2개이며, 직선의 방정식은 밝은 두 점이 표현하는 $(r, \theta)$로 결정된다.
[참고문헌]
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(1.1)
