라돈변환과 FFT 그리고 CT

FFT의 실생활 적용

푸리에 변환(Fourier transform, FT)은 시간이나 공간에 대한 함수를 시간 또는 공간 주파수 성분으로 분해하는 변환을 말한다.

푸리에 급수/변환의 목적은 어떤 함수 f(t)에 대해 주어진 시간 축에서 다수의 사인파 또는 코사인파로 분해하는 것이다

FFT는 우리가 인식하지 못할 뿐, 음악,영상,의료등 많은 영역에서 사용되고 있다.

이중 의료 영상 중 컴퓨터 단층촬영 (CT)에 도입해 볼 것이다.

컴퓨터 단층촬영 (CT)

컴퓨터단층촬영, CT(Computed Tomography)는 X선 발생장치가 있는 원통형의 기계를 사용해 내부를 촬영하는 기술, 또는 그러한 기술을 사용한 의학 검사이다.

위 사진과 같이 X선을 몸에 투영시켜 적분을 통해 감쇠로 인한 값을 저장하여 사진으로 복원한다

즉, CT촬영 결과의 원본은 우리가 흔히 아는 사진이 아닌 Sinogram으로 이루어져 있다

Sinogram이 어떻게 얻는지; 우리가 흔히 아는 CT사진이 되는지 알아보자

Radon Transform

수학에서 라돈 변환은 평면에서 정의된 함수를 2차원 공간에서 정의된 함수로 변환하는 적분변환이다
특정 선에서의 값은 해당 선위의 함수의 선적분과 같다
정리하자면, 라돈변환의 핵심 연산은 선적분이며 CT촬영에서 몸에 X선을 여러 각도로 투과시키면서 선적분을 하고, 그 X선 투과결과가 Sinogram이다.

$R f(\theta, s) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \, \delta(x \cos\theta + y \sin\theta - s) \, dx \, dy$

Inverse Radon Transform

라돈변환으로 Sinogram을 얻는다하더라도 직관적으로 이를 분석하기 어렵다.
따라서 Sinogram을 다시 역변환을 해서 우리의 몸 단층의 이미지를 얻는다.
단순히 Sinogram을 다방면으로 적분을 해서 원본의 이미지를 얻을것이라 생각할 수도 있지만, 그렇지 않다
적분만으로 역변환을 할 경우, 이미지가 퍼져서 흐릿해보일 수 있기 때문에 다음과 같은 과정을 거쳐야 한다.

Frequency domain filtering(주파수 도메인 필터링)

BackProjection(백프로젝션)

Frequency Domain Filtering (주파수 도메인 필터링)

위에 말했듯이 라돈변환은 저주파가 강조되서 단순히 적분만 하면 흐릿해 지기 때문에 먼저 필터링 과정을 거친다.
크게 보면 신호나 이미지르르 더 정확히 해석하고 원하는 정보만 추출하기 위해서이다.

필터링 과정

각도 𝜃마다 $R_f(\theta, s)$ 을 추출한다

s에 대해 푸리에 변환을 사용해서 주파수 변환한다

주파수로 바꾸면 원하는 부분만 쉽게 조작할 수 있다

주파수에 $H(\omega)$ 곱한다 (대표적으로 Ram-lak filter)

다시 푸리에 역변환을 한다

Fourier Transform 사용 이유

시간 영역에서 필터링은 복잡하지만 주파수 영역에서는 단순 곱셈으로 나온다

\[x(t) * h(t) \xrightarrow{\mathcal{F}} X(\omega) \cdot H(\omega)\]

Ram-Lak Filter

위에 주파수에 곱하는 필터 함수에서 대표적으로 사용하는 것이다.

$H(\omega) = |\omega| \quad \text{(high-pass filter)}$

저주파 성분 억제 (라돈변환은 저주파가 강조된다)

고주파 강조 (경계를 살려서 흐림을 방지한다)

Backprojection (백프로젝션)

백프로젝션(Backprojection)이란, 라돈 변환(선적분)으로 얻은 Sinogram 데이터를 이미지 평면으로 되쏘아 원래 이미지를 복원하는 과정이다.
라돈변환은 여러 각도에서 x선방향으로 선적분 한 것이기 때문에 선적분 값을 그 직선을 따라 뿌리고 모든 각도로 반복해서 누적시키면 원래의 이미지가 복원된다.

$f(x, y) \approx \int_0^\pi \text{FilteredProjection}(\theta, s)\, d\theta \quad \text{where } s = x \cos \theta + y \sin \theta$

(x,y) : 복원 할 이미지의 좌표

𝜃 : 각도

s : 점위 해당되는 투영 직선 위치

모든 각도별로 sinogram에서 해당 좌표값을 다 더한다

위에서 말했듯이 그냥 backproject을 할 경우 흐릿해지기 때문에 필터링을 거친다

이제 이를 코드로 구현해 보자

Implement in Code

source code
Python 소스코드 파일첨부 | CT Medical Image
run
32.53s

Radon Transfrom

과정 요약

라돈변환은 이미지를 여러각도로 회전 시켜야 하는데 회전 과정에서 이미지가 잘릴 수 도 있다.

때문에 잘리는것을 방지하기위해 잘리지 않는 최소크기(대각선) 만큼 패딩을 해준다

원래는 CT가 회전하면서 투영하는 것을 이미지를 회전하면서 투영하는 것으로 구현한다.

회전된 이미지에 대해 선적분을 한다.

결과를 sinogram으로 쌓고 반환한다.

여기서 선적분은 실제로 선적분 연산을 하지는 않는다. 이미지가 픽셀로 조밀하기 때문에 수치근사한다

디지털 이미지는 연속이 아닌 discrete하기 때문에 적분대신 픽셀 줄 합산을 한다

실제로도 CT 스캐너도 무한한 해상도로 선적분하는 것이 아닌 이산데이터를 받아서 수치근사한다

# 이미지는 흑백이미지를 담는다 (라돈변환은 색보다 밝기정보에 민감)
# 또한 흑백은 2D (64x64)지만 RGB는 3D (64x64x3) 배열이여서 메모리 절약도 가능
def radon_transform(image, angles):

    # h : 높이 w : 너비
    h, w = image.shape

    # 대각선 길이 계산
    diag_len = int(np.ceil(np.sqrt(h**2 + w**2)))

    # zero-pedding (이미지 회전시 손상 방지)
    pad_width = ((diag_len - h) // 2, (diag_len - w) // 2)
    padded = np.pad(image, (pad_width, pad_width), mode='constant', constant_values=0)

    sinogram = []

    for theta in angles:
        # 이미지 회전 (reshape = False : 이미지 크기유지, order = 1: 1차 선형보간)
        # 선형보간 : 두점사이의 직선의 방정식 (order = 1 : 1차)
        rotated = rotate(padded, angle=theta, reshape=False, order=1)

        # 각도에서 수직 방향으로 선적분 (픽셀 합)
        projection = np.sum(rotated, axis=0)
        sinogram.append(projection)

    # 결과값 방향이 90도 돌아가있어서 전치행렬로 바꿔줌
    return np.array(sinogram).T

Fast Fourier Transform frequency

FFT 결과에 해당하는 주파수 배열 생성 함수
해상도를 높이는 Ram-Lak Filter는 주파수 필터이기 때문에 주파수로 바꿔줘야한다.

# n : 샘플 개수, d : 샘플 간격
def fftfreq(n, d=1.0):

    # 주파수 간격 = 1 / 전체 구간 길이(n * d)
    Hz = 1.0 / (n * d)

    results = np.arange(n)

    # FFT, DFT의 결과는 양수 음수가 공존해있다
    # 음부부분도 표현하기 위해 뒤쪽 음수부분을 양수로 바꿔준다
    negative_n = (n + 1) // 2
    results[negative_n:] -= n

    # 실제 주파수 = 주파수인덱스 * 단주파수
    return results * val

Ram-Lak Filter

역 라돈 변환에서 첫번째 과정인 필터링 과정에서 필수적인 함수이다

과정 요약

앞서 만든 fft_freq()로 주파수 축을 생성한다

주파수 축으로 Ram-Lak Filter를 생성한다

FFT 결과 주파수 축을 곱해 주파수 필터링을 한다

배열을 다시 원래 크기로 만들어 준 후 리턴한다

# 필터링 시킬 sinogram과 주파수 변환에 필요한 FFT, fft_freq를 넣는다
def Ram_Lak(sinogram, FFT, fft_freq):

    # 한 방향에서의 픽셀
    n = sinogram.shape[0]

    # 주파수 축 생성
    # 푸리에 변환 결과의 인덱스의 Hz
    Hz = fft_freq(n)

    # 주파수의 크기만 필요하기 때문에 절대치 씌운다
    # 고역통과 필터 (high-pass filter)
    ramp = np.abs(Hz)

    # 필터링 결과 배열
    filtered = np.zeros_like(sinogram, dtype=np.float64)

    # 각도별로 반복
    for i in range(sinogram.shape[1]):

        # i 번째 각도 선적분 결과
        signal = sinogram[:, i]
        signal_len = len(signal)

        # FFT를 위한 zero-pedding 수행
        padded_len = 2 ** int(np.ceil(np.log2(signal_len)))
        signal_padded = np.pad(signal, (0, padded_len - signal_len))

        # FFT 수행 (주파수 도메인으로 변환)
        spectrum = FFT(signal_padded)

        # 결과값을 담을 필터도 zero-pedding으로 길이 맞춰준다
        ramp_padded = np.pad(ramp, (0, len(spectrum) - len(ramp)))

        # Ram-lak filter의 목적
        # 주파수 필터링은 단순한 곱셈으로 이루어진다
        # 시간 도메인에서 필터링은 복잡하지만 주파수는 필터와 곱하기만 하면 된다
        spectrum_filtered = spectrum * ramp_padded

        # 역 푸리에 변환
        inverse_fft = np.conj(FFT(np.conj(spectrum_filtered))) / len(spectrum_filtered)

        # 실수 부만 추출, 원래 길이로 되돌린다
        filtered[:, i] = np.real(inverse_fft[:signal_len])

    return filtered

Inverse Radon Transform

def inverse_radon(sinogram, theta, FFT, fft_freq, size = None):
    
    h, w = sinogram.shape

    # 사이즈를 지정하지 않았을경우 sinogram크기에 맞춘다
    if size is None:
        size = h

    # sinogram을 주파수 도메인 필터링 한다
    filtered = Ram_Lak(sinogram, FFT, fft_freq)

    # 역변환 결과 배열 생성
    inverse_transform = np.zeros((size, size), dtype = np.float32)

    # 역변환 결과 담을 좌표 생성
    # 이미지를 중심에 오게하기 위해 좌표를 중심으로 맞춘다
    mid = size // 2
    x = np.arange(size) - mid
    y = np.arange(size) - mid
    # 이미지 중심에 오게하는 (x, y)좌표 쌍 생성
    X, Y = np.meshgrid(x, y)

    # 마찬가지로 각도별로 Backprojection 수행한다
    for i, angle in enumerate(theta):

        # sinogram 한줄 추출한다
        sinogram = filtered[:, i]

        # 역변환 이므로 각도는 반대로 돌린다
        line = np.tile(sinogram, (size, 1))

        # 라돈변환과 동일하게 누적시킨다
        rotated = rotate(line, angle = -angle, reshape=False, order=1)
        inverse_transform += rotated

    # 정규화(Nomalization)
    # 각도별로 누적된 값이 커지기 때문에 너무 밝아지거나 왜곡될 수 있다
    # 보정 계수 :pi / 2N(각도 수)
    inverse_transform *= np.pi / (2 * len(theta))

    return inverse_transform

Data Loading

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import rotate
from skimage.transform import radon

# dataset 파일 확인을 위해 주피터 노트북에서 실행
import numpy as np

data = np.load('dataset.npz', allow_pickle=True)
image = data['image']

print("type:", type(image))
print("shape:", getattr(image, 'shape', 'No shape'))
print("dtype:", image.dtype if isinstance(image, np.ndarray) else 'Not ndarray')

475개의 이미지를 갖고있다

data = np.load('dataset.npz', allow_pickle=True)

## 첫번째 이미지 사용
image = data['image'][0]

# 컬러를 흑백으로 변환
if image.ndim == 3 and image.shape[2] == 3:
    image = np.mean(image, axis=2)

# 0도부터 180도까지 균등하게 각도 생
theta = np.linspace(0., 180., max(image.shape), endpoint=False)

Use Function

sinogram = radon_transform(image, theta)

reconstructed = inverse_radon(sinogram, theta, FFT, fft_freq)

Visualization

# 라돈 역변환 적용
reconstructed = inverse_radon(sinogram, theta, FFT, fft_freq)

plt.figure(figsize=(18, 6))

plt.subplot(1, 3, 1)
plt.title("원본 이미지")
plt.imshow(image, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.title("Sinogram")
plt.imshow(sinogram, cmap='gray', aspect='auto')
plt.xlabel("Angle")
plt.ylabel("Detector pixel")

plt.subplot(1, 3, 3)
plt.title("Inverse Radon Transfrom")
plt.imshow(reconstructed, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()

output

dataset 0번째 이미지

dataset 1번째 이미지