[Bio] #3 CT(Computed Tomography), matlab으로 CT 구현하기(라돈 역변환)

 

우리가 흔히 씨티(CT) 촬영이라고 말하는

컴퓨터단층촬영(Computed Tomography scan)에 대해 알아보자.

 

 

---

CT scanner

CT란?

CT scanner의 구조를 보면 알 수 있듯이,

여러 방향에서 엑스선 투과기법으로

인체 단면을 촬영하여 복원하는 의료영상촬영법이다.

 

출처: 삼성서울병원 홈페이지

 

위의 사진처럼 인체의 횡단면 하나에 대해,

전 방향에서 X선을 송, 수신하여 여러장의 영상을 촬영해

횡단면 하나의 영상으로 합한다.

 

X선을 투과하는 모습(왼)과 sinogram(오)

 

왼쪽 사진의 흰색 그래프를 보면,

투과한 X선들이 한 평행한 센서에 감지되는 것을 알 수 있다.

이때 센서에 맺힌 X선 흡광도를 여러 방향에 대해 합한 뒤

그래프로 시각화한 것이 오른쪽 사진의 사이노그램(sinogram)이다.

 

 

 

여기서 짚고 넘어갈 개념이 '라돈 변환(Radon transform)'이다.

 

2차원 Radon Transform

 

Rf(α,s)의 기하적인 의미는,

함수 f를 원점에서 s(거리)만큼 떨어져 있고 α(각)와 수직인 모든 점에서 적분한 것이다.

즉 특정 각도와 거리에 대한 투과도 값을 갖는 projection 데이터가 된다.

 

sinogram은 이러한 라돈 변환(Radon transform)을 물리적으로 구현한 결과이며,

우리는 sinogram의 역 라돈 변환을 통해 원본 데이터를 재구성할 수 있다.

 

 

---

matlab 내장 함수로 팬텀 이미지의 sinogram을 뽑고 다시 역변환 해보자.

 

obj = phantom(256);
imagesc(obj); colormap gray; axis image; axis off;

Shepp Logan phantom image

위와 같이 256*256 size의 phantom 이미지를 확인할 수 있다.

phantom(n)은 n*n phantom 배열을 반환하는 함수이며,

imagesc는 배열을 이미지로 시각화해주는 함수다.

 

thetas = linspace(0,180,256);
[sinogram, x] = radon(obj, thetas);
imagesc(sinogram);

 

각을 0부터 180도까지 일정하게 256개로 나눠

built-in radon function으로 해당 각들에 대해 projection해주면 sinogram을 얻을 수 있다.

 

이때 [sinogram, x] = radon(obj, thetas); 에서

sinogram은 라돈변환의 결과로 얻어지는 projection 데이터를 갖는 2차원 배열이다.

x는 라돈 변환 결과의 거리 값 데이터인데, 여기선 필요없다.

 

이를 통해 얻은 sinogram은 다음과 같다.

 

phantom 이미지의 sinogram

 

 

CT의 원리대로라면

우리가 내장 phantom 함수를 이용해 얻은 phantom이미지가 진짜 인체가 될 것이고,

X선을 투과해 얻은 가공 전의 데이터가 sinogram이다.

이 sinogram을 역변환 하기 전까지 우리는 인체 내부 이미지를 알 수 없다.

 

이제 역변환을 통해 우리가 CT 촬영으로 보게 될 이미지를 얻어보자.

 

---

matlab에 내장된 라돈 역변환 함수는 다음과 같이 사용할 수 있다.

I = iradon(R,theta)
I = iradon(R,theta,interp,filter,frequencyScaling,outputSize)

filter를 사용해 엄밀히 변환하면 더 좋겠지만,

일단 간단히 첫번째 방식으로 역변환해보자.

 

I = iradon(sinogram,thetas);
imagesc(I); colormap gray;

내장 iradon function을 사용해 역변환

 

오!

sinogram의 역변환으로 원본 데이터와 거의 일치하는 이미지를 얻었다.

 

 

 

 

 

이번에는 직접 라돈 역변환 코드를 짜보자.

여기서 Back-projection 방법을 사용할 것이다.

 

이 방법은 쉽게 말해 sinogram의 각 열 데이터(한 방향에 대한 projection 데이터)를 추출해

해당 각만큼 회전시킨 뒤 다시 합하는 것이다.

 

[row, col] = size(sinogram);
backproj = zeros(size(sinogram));
for i=1:1:length(thetas)
    backproj = backproj + imrotate(repmat(sinogram(:, i), [1,col]), -thetas(i), 'bilinear', 'crop');
end
backproj=backproj';
figure;
imagesc(backproj); colormap gray; axis image

 

반복문 안을 보자.

sinogram(:, i)를 col번 반복하여 2차원 배열을 만들고, 각에 대해 회전시켜 합하고 있다.

이렇게 하면 sinogram의 열 데이터가 행 방향으로 복사되고, 그 이미지를 해당 각에 맞게 회전시켜 합하게 된다.

 

(그냥 backproj를 출력하면 90도 돌아가 있길래 전치시켜주었다.)

 

back-projection 방식으로 역변환한 이미지는 다음과 같다.

 

 

오!

 

희뿌연 느낌이 섞이고 왜곡, blur 효과가 가미된 듯하긴 하지만..

원본 이미지를 알아볼 수 있다.

 

 

 

---

 

 

앞으로 바이오 쪽을 더 공부하거나 연구할 일이 있을지 모르겠지만..

전에 공부했던 내용들이 한번 들어서는 잊어버릴 것 같은 것들이라,

종종 생각날 때 정리하러 와야겠다!

 

 

 

 

2023.07.25 - [Study/Bio Technology] - [Bio] 바이오전기전자 #2 Matlab으로 DFT 구현하기

[Bio] 바이오전기전자 #2 Matlab으로 DFT 구현하기