라돈 변환

1. 개요[편집]

라돈 변환 Radon Transform
제안	Johann Radon (1917) — CT보다 54년 빨랐다
정의	$p(\theta,s)=\iint f(x,y)\,\delta(x\cos\theta+y\sin\theta-s)\,dx\,dy$
자료 형태	사이노그램(sinogram) — 점 하나가 사인 곡선 하나
핵심 정리	중앙 슬라이스 정리: 투영의 1D 푸리에 = 2D 푸리에의 방사 절단
표준 역변환	필터보정역투영(FBP) — 램프 필터 $\lvert\nu\rvert$
병적 정도	소볼레프 척도 1/2계 — 약하게 잘못된 문제
응용	CT · PET/SPECT · 전자 토모그래피 · 탄성파 $\tau$-$p$ 변환

물체를 자르지 않고 속을 보려면, 사방에서 그림자를 찍어 두고 그 그림자들로 역산하면 된다.

라돈 변환은 평면 위의 함수 $f(x,y)$ 를, 평면의 모든 직선에 대한 선적분 값의 모임으로 바꾸는 적분변환이다. 직선을 각도 $\theta$ 와 원점으로부터의 부호 있는 거리 $s$ 로 매개변수화하면

$$p(\theta,s) \;=\; \mathcal{R}f(\theta,s) \;=\; \iint f(x,y)\,\delta\bigl(x\cos\theta + y\sin\theta - s\bigr)\,dx\,dy$$

가 된다. 델타 함수가 하는 일은 “그 직선 위만 골라 더하라”는 것뿐이다. 물리적으로는 각도 $\theta$ 방향에서 평행하게 쏜 X선 하나가 물체를 지나며 겪은 총 감쇠량이고, 그래서 $p(\theta,\cdot)$ 를 투영(projection)이라 부른다.

요한 라돈이 이 변환과 그 역공식을 발표한 것이 1917년이다. 응용을 염두에 둔 것이 전혀 아니었고 순수 적분기하학 논문이었는데,¹ 반세기 뒤 코맥이 독립적으로 재발견하고 하운스필드가 EMI에서 실물 스캐너를 만들면서 의학을 바꿔 놓았다. 두 사람은 1979년 노벨 생리의학상을 받았고, 라돈은 1956년에 이미 세상을 떠난 뒤였다.

2. 사이노그램은 왜 사인 곡선인가[편집]

$p(\theta,s)$ 를 $\theta$ 를 세로축, $s$ 를 가로축으로 그린 2차원 영상이 사이노그램이다. 이름의 유래는 계산 한 줄이다. 원 영상이 점 하나, 즉 $f = \delta(x-x_0,\,y-y_0)$ 라면

$$p(\theta,s) = \delta\bigl(s - (x_0\cos\theta + y_0\sin\theta)\bigr) = \delta\bigl(s - r\cos(\theta - \phi)\bigr), \qquad r=\sqrt{x_0^2+y_0^2},\;\; \phi = \operatorname{atan2}(y_0,x_0)$$

즉 점 하나가 사이노그램에서는 사인 곡선 하나를 그린다. 진폭이 원점으로부터의 거리, 위상이 방위각이다. 일반 영상은 이 사인 곡선들의 중첩이므로, 사이노그램을 보면 물체가 회전축에서 얼마나 벗어나 있는지, 금속처럼 감쇠가 극단적인 물체가 어디 있는지를 눈으로 읽을 수 있다. CT 기사들이 재구성 전에 사이노그램을 먼저 보는 이유다.

같은 사인 곡선이 허프 변환에도 나온다. 우연이 아니고, 아래에서 다시 다룬다.

3. 중앙 슬라이스 정리[편집]

라돈 변환 이론의 전부라고 해도 좋은 정리다. 투영 $p(\theta,\cdot)$ 를 $s$ 에 대해 1차원 푸리에 변환한 것을 $P_\theta(\nu)$, 원 영상의 2차원 푸리에 변환을 $F(u,v)$ 라 하면

$$P_\theta(\nu) \;=\; F(\nu\cos\theta,\;\nu\sin\theta)$$

한 방향의 투영을 1차원 푸리에 변환하면, 그것이 2차원 푸리에 평면에서 원점을 지나는 각도 $\theta$ 짜리 직선 위의 값이라는 뜻이다. 증명은 정의를 대입하고 적분 순서만 바꾸면 세 줄이다.

$$P_\theta(\nu) = \int p(\theta,s)e^{-i2\pi\nu s}ds = \iint f(x,y)\,e^{-i2\pi\nu (x\cos\theta+y\sin\theta)}dx\,dy = F(\nu\cos\theta,\nu\sin\theta)$$

이 정리가 주는 것은 “측정하면 무엇을 아는가”에 대한 완전한 답이다. 각도 하나를 측정할 때마다 푸리에 변환 평면의 직선 하나를 채운다. 각도를 촘촘히 돌리면 평면이 채워지고, 2차원 역푸리에로 $f$ 를 얻는다. 존재 정리이자 구성적 알고리즘이며, 동시에 왜 각도가 부족하면 안 되는지도 즉시 설명한다 — 안 잰 각도는 푸리에 평면에 빈 부채꼴로 남는다.

직접 구현(그리딩, direct Fourier method)도 물론 가능하고 고속 푸리에 변환을 쓰면 가장 빠르다. 다만 극좌표 표본을 직교 격자로 보간해야 하는데, 표본이 저주파에 조밀하고 고주파에 성기게 분포하므로 보간 오차가 고주파에서 커지고 그게 영상 전체에 퍼지는 아티팩트가 된다. 그리딩 커널 설계가 그래서 별도 연구 주제이고, 실무의 기본값은 다음 절의 FBP다.

4. 그냥 역투영하면 왜 뭉개지는가[편집]

가장 순진한 발상은 이렇다. 각 투영값을 그것이 나온 직선을 따라 도로 뿌리고(smear back), 모든 각도에 대해 더한다.

$$b(x,y) = \int_0^\pi p\bigl(\theta,\;x\cos\theta+y\sin\theta\bigr)\,d\theta$$

이것이 역투영(backprojection)이고, 수학적으로는 라돈 변환의 수반 연산자($\mathcal{R}^{*}$)이지 역연산자가 아니다. 실제로 계산해 보면

$$b = \mathcal{R}^{*}\mathcal{R}f = f * \frac{1}{r}, \qquad r=\sqrt{x^2+y^2}$$

즉 원 영상을 $1/r$ 이라는 꼬리가 긴 핵으로 흐린 것이 나온다. 2차원에서 $1/r$ 의 푸리에 변환이 $1/\lvert\nu\rvert$ 이므로, 주파수 영역에서 보면 모든 성분이 $\lvert\nu\rvert$ 로 나눠져 있다. 저주파는 과대평가되고 고주파는 죽는다. 그래서 순수 역투영 영상은 물체 형상이 어렴풋이 보이는 안개 덩어리가 된다.

원인은 표본 밀도다. 중앙 슬라이스 정리에 따라 각 투영이 채우는 것은 원점을 지나는 직선이고, 직선들은 원점 근처에서 빽빽하고 멀수록 성기다. 각도를 아무리 촘촘히 돌려도 이 불균형은 남는다. 저주파를 여러 번 세고 고주파를 적게 세는 셈이니, 결과가 저주파로 편향되는 게 당연하다.

5. 필터보정역투영[편집]

고칠 곳이 정확히 한 군데라는 게 아름답다. 2차원 역푸리에 적분을 극좌표로 바꾸면 야코비안이 튀어나온다.

$$f(x,y) = \iint F(u,v)\,e^{i2\pi(ux+vy)}\,du\,dv = \int_0^\pi\!\!\int_{-\infty}^{\infty} F(\nu\cos\theta,\nu\sin\theta)\,e^{i2\pi\nu s}\,\lvert\nu\rvert\,d\nu\,d\theta$$

$du\,dv = \lvert\nu\rvert\,d\nu\,d\theta$ 의 그 $\lvert\nu\rvert$ 가 앞 절에서 잃어버린 바로 그 인자다. 중앙 슬라이스 정리로 $F$ 를 $P_\theta$ 로 바꿔 쓰면 필터보정역투영(filtered backprojection, FBP) 공식이 나온다.

$$f(x,y) = \int_0^\pi \tilde{p}\bigl(\theta,\;x\cos\theta + y\sin\theta\bigr)\,d\theta, \qquad \tilde{P}_\theta(\nu) = \lvert\nu\rvert\,P_\theta(\nu)$$

즉 역투영하기 전에 각 투영을 램프 필터 $\lvert\nu\rvert$ 로 한 번 필터링한다. 이 한 줄이 40년 넘게 임상 CT의 표준 알고리즘이었다. 순서를 바꿔 “역투영 후 2차원 필터링”을 해도 수학적으로는 같지만, 필터링을 1차원으로 하는 쪽이 압도적으로 싸고 투영이 도착하는 대로 처리할 수 있어 스트리밍에도 유리하다. 비용은 각도 $M$, 검출기 칸 $N$, 영상 $N \times N$ 일 때 필터링 $O(MN\log N)$, 역투영 $O(MN^2)$ 으로 역투영이 지배한다.

램프 필터의 공간영역 임펄스응답은 발산하는 물건이라 실제로는 나이퀴스트에서 자른 이산판을 쓴다. 라마찬드란과 락슈미나라야난이 유도한 Ram-Lak 커널은 검출기 간격 $\Delta s$ 에 대해 $h(0) = 1/(4\Delta s^2)$, 짝수 $n$ 에서 0, 홀수 $n$ 에서 $-1/(\pi^2 n^2 \Delta s^2)$ 인 예쁜 닫힌 형식을 갖는다. 문제는 $\lvert\nu\rvert$ 가 고주파를 그대로 통과시켜 잡음을 그대로, 아니 증폭해서 넘긴다는 것이다. 그래서 램프에 창을 씌운다.

창	성격	결과
Ram-Lak	램프를 나이퀴스트에서 직절단	해상도 최대, 잡음 최대, 링잉
Shepp-Logan	램프에 sinc 형태를 곱함	고주파를 완만히 죽임 — 실무 기본값
코사인 · 해밍 · 한	더 강한 저역 감쇠	잡음 억제 우선, 해상도 손실

이 선택이 곧 해상도와 잡음의 저울이며, 정칙화 파라미터를 고르는 문제와 정확히 같은 구조다(티호노프 정규화·역문제 참고). 저선량 CT에서 영상이 뿌옇게 보이는 이유의 상당 부분은 창을 세게 먹였기 때문이다.²

6. 얼마나 재야 하는가, 덜 재면 어떻게 되는가[편집]

표본화. 검출기 간격 $\Delta s$ 가 정하는 나이퀴스트 주파수까지를 각도 방향으로도 똑같이 채우려면, 검출기 칸이 $N$ 개일 때 각도는 대략 $M \approx \frac{\pi}{2}N$ 개가 필요하다는 경험칙이 있다. 512칸이면 800뷰 남짓. 임상 스캐너의 뷰 수가 그 근처인 것이 우연이 아니다.

각도 부족(sparse-view). 각도를 줄이면 푸리에 평면에 방사형 빈틈이 생기고, 그 결과가 고대비 구조에서 뻗어 나오는 줄무늬 아티팩트(streak)다. 금속 임플란트 주변이 방사형으로 지저분해지는 것이 이것이고(금속 아티팩트는 빔 경화까지 겹쳐 더 나쁘다), 대책은 아래의 반복법이나 압축센싱형 정칙화다.

제한각(limited angle). 기하학적 이유로 $[0,\pi)$ 를 다 돌 수 없는 경우가 있다. 전자 토모그래피에서 시료 홀더 때문에 $\pm 60^\circ \sim \pm 70^\circ$ 밖에 기울이지 못하는 것이 대표적이다. 이때 푸리에 평면에는 빠진 쐐기(missing wedge)가 남고, 그 방향의 정보가 아예 없으므로 재구성 결과가 빔 방향으로 늘어나고 그 방향 경계면이 사라진다. 이건 알고리즘으로 완전히 메울 수 있는 문제가 아니다 — 데이터에 없는 것을 만들어 내려면 사전 정보(양수성, 희소성, 형상 모형)를 넣는 수밖에 없다.

이론적으로 라돈 변환은 소볼레프 척도에서 1/2계만큼 매끄럽게 만드는 연산자다. 따라서 역변환은 잘못된 문제이긴 하나 정도가 약하고(mildly ill-posed), 미분방정식의 후향 문제나 디컨볼루션처럼 지수적으로 발산하지는 않는다. 그래서 완전한 데이터만 있으면 FBP 같은 명시적 공식이 실용적으로 작동한다. 재미있는 여담으로, 2차원 라돈 역변환은 비국소적이다 — 한 점의 값을 알려면 그 점을 지나지 않는 직선의 적분값까지 필요하다. 홀수 차원에서는 역공식이 국소 미분만으로 쓰이지만 짝수 차원에서는 힐베르트 변환이 끼어들기 때문이며, 그래서 관심 영역만 스캔해 그 부분만 복원하는 “내부 문제”(interior problem)가 2차원에서는 추가 가정 없이 유일해가 없다.

7. 반복 재구성[편집]

FBP는 “데이터가 충분하고 잡음이 적다”는 전제 위의 명시적 공식이다. 전제가 깨지면 문제를 아예 거대한 선형계로 놓는 편이 낫다.

$$A\mathbf{x} = \mathbf{b}$$

$\mathbf{x}$ 는 화소값 벡터, $\mathbf{b}$ 는 측정된 투영값, $A_{ij}$ 는 $i$ 번 광선이 $j$ 번 화소를 지나는 길이다. $A$ 는 극도로 희소하지만 크기가 $10^6 \times 10^6$ 을 넘고, 대개 과소결정이거나 모순된 계다.

• ART(대수적 재구성법). 본질은 1937년 카츠마시 방법이다. 광선 하나(즉 행 하나)를 골라, 현재 추정값을 그 행이 정의하는 초평면에 사영한다.

  $$\mathbf{x}^{k+1} = \mathbf{x}^{k} + \lambda\,\frac{b_i - \mathbf{a}_i^{\mathsf T}\mathbf{x}^{k}}{\lVert \mathbf{a}_i\rVert^2}\,\mathbf{a}_i$$

  1970년 고든·벤더·허먼이 CT 맥락에서 재발견해 대수적 재구성법이라는 이름을 붙였다. 수렴이 빠르지만 행 순서에 민감하고, 잡음이 있으면 해가 없는 계를 억지로 맞추느라 소금-후추 잡음이 낀다.

• SIRT. 한 번에 모든 행을 반영해 갱신한다. 사실상 $\lVert A\mathbf{x}-\mathbf{b}\rVert^2$ 에 대한 란트베버 반복 내지 경사하강이며, 느리지만 잡음에 훨씬 부드럽다. SART(1984)는 뷰 단위로 갱신하는 절충안이다.

• 통계적 재구성. 방출 단층촬영(PET·SPECT)에서는 계수가 푸아송 분포를 따르므로 최소제곱이 아니라 우도를 최대화해야 하고, 그 표준 알고리즘이 기댓값 최대화 알고리즘 기반 MLEM과 그 가속판 OSEM이다. X선 CT의 MBIR도 같은 계보로, 잡음 모형 + 사전항(대개 에지 보존형)을 명시적으로 넣는다.

• 정칙화. 사전항으로 전변분 잡음제거의 TV 항을 쓰는 것이 sparse-view CT의 국룰이 됐다. 뼈처럼 조각별 상수에 가까운 영상에 특히 잘 듣는다. 압축센싱 이론이 “왜 뷰를 절반으로 줄여도 되는가”에 대한 정당화를 제공한다.

정리하면 FBP는 빠르고 예측 가능하며 선형이라 잡음 전파를 손으로 계산할 수 있고, 반복법은 느리지만 잡음 모형·기하 제약·사전 정보를 넣을 수 있다. 저선량이 화두가 된 2010년대 이후 임상 스캐너가 반복법으로 대거 이동했고,³ 지금은 그 자리를 다시 학습 기반 재구성이 넘보고 있다. 다만 재구성망이 없는 병변을 그려 넣는(hallucination) 위험 때문에, 계측·진단 영역에서는 “데이터 일관성 항을 반드시 유지한다”는 하이브리드 설계가 표준이다.

8. 허프 변환과의 관계[편집]

이진 영상 $f$ 에 대해 $\mathcal{R}f(\theta,s)$ 는 직선 $(\theta,s)$ 위에 놓인 켜진 화소의 개수다. 그런데 허프 변환의 누적기 $A[\theta,s]$ 에 쌓이는 표의 수도 정확히 같은 값이다. 두 변환은 같은 배열을 만든다. 차이는 계산 순서에 있다.

• 허프는 흩뿌리기(scatter)다. 화소 하나가 자기를 지나는 모든 직선에 표를 던진다.
• 라돈은 긁어모으기(gather)다. 직선 하나가 자기 위 화소를 전부 더한다.

역투영이 라돈 변환의 수반이듯, 이 둘도 서로 수반 관계이며 그래서 결과가 같다. 실질적인 차이는 목적이다. 허프는 봉우리 위치만 원하므로 그레이스케일 값도, 역변환도 필요 없고 대신 기울기 방향으로 투표 범위를 좁히거나 임의 형상으로 일반화하는 확장이 자유롭다. 라돈은 배열 전체가 데이터이고 되돌리는 것이 목적이라 표본화·잡음·역문제 이론이 따라붙는다. 같은 사인 곡선 그림을 보면서 한쪽은 봉우리를, 다른 쪽은 영상을 본다.

9. 어디에 쓰이나[편집]

• X선 CT. 원조 응용. 실제 스캐너는 평행빔이 아니라 부채빔·원뿔빔이라 좌표 변환(rebinning)이나 전용 공식이 필요하고, 3차원 원뿔빔의 사실상 표준이 FDK 근사 알고리즘이다. 원뿔빔의 정확한 재구성을 위한 궤적 조건(투이 조건)은 원 궤도 하나로는 만족되지 않으며, 나선형 궤도가 그래서 쓰인다.
• 전자 토모그래피·크라이오 전자현미경. 앞서 말한 빠진 쐐기 문제가 이 분야의 상수다.
• 탄성파 탐사. $\tau$-$p$(slant stack) 변환이 곧 라돈 변환이고, 다중 반사파 제거와 속도 분석에 쓴다. 지진 토모그래피는 광선이 직선이 아니라 굴절되는 만큼 비선형 역문제가 된다.
• 영상처리 일반. 직선·줄무늬 검출, 방향 질감 분석, 스캔 문서의 기울기 보정. 사이노그램에서 봉우리 각도를 읽으면 그게 곧 기울어진 각도다.
• 수치해석 실습. 셰프-로건 팬텀처럼 타원들의 합으로 만든 전산 단층촬영 시험 영상은 라돈 변환의 해석해가 존재해서, 재구성 알고리즘의 검증(V&V)에 이상적인 벤치마크가 된다. 타원 하나의 투영은 닫힌 형식이므로 이산화 오차를 참값과 직접 비교할 수 있다.

10. 관련 문서[편집]

• 푸리에 변환 · 고속 푸리에 변환 · 웨이블릿 변환
• 허프 변환 · 에지 검출 · 특징점 검출
• 역문제 · 티호노프 정규화 · 압축센싱 · 조건수
• 최소자승법 · 켤레기울기법 · 기댓값 최대화 알고리즘
• 전변분 잡음제거 · 비등방성 확산 · 이미지 피라미드
• 볼륨 렌더링 · 라돈 측도
• 전산 단층촬영 · 대수적 재구성법 · 지진 토모그래피

11. Footnotes[편집]

1. Radon, J. (1917). “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten.” 제목을 옮기면 “어떤 다양체를 따른 적분값으로 함수를 결정하는 것에 관하여”다. 응용 이야기는 한 줄도 없다. 순수수학이 반세기 묵혔다가 병원 장비가 되는 이 서사는 응용수학 강의 첫 시간의 단골 소재이며, 참고로 라돈 측도의 그 라돈과 동일 인물이다. ↩

2. 램프 필터를 처음 구현해 보면 “왜 투영값에서 음수가 나오지” 하고 당황하게 된다. 램프는 미분에 가까운 고역 강조 필터라 필터링된 투영은 당연히 음수를 갖는다. 그 음수가 역투영 과정에서 이웃한 양수 기여를 상쇄해 $1/r$ 꼬리를 정확히 지우는 것이 FBP의 작동 원리다. 음수가 싫다고 클리핑하면 영상 전체가 뿌옇게 뜬다 — 실제로 자주 나오는 초보 버그. ↩

3. “저선량 CT”의 저선량은 대체로 광자 수가 적다는 뜻이고, 광자 수가 적으면 푸아송 잡음이 커진다. 여기서 알고리즘이 하는 일은 잡음을 없애는 게 아니라 잡음과 해상도의 교환비를 유리하게 바꾸는 것이다. 물리적으로 안 찍힌 광자는 어떤 알고리즘으로도 되살아나지 않는다. 이 당연한 사실이 마케팅 자료에서는 종종 흐려진다. ↩