計(jì)算機(jī)斷層掃描算法 (computed tomography)

很早我就對(duì)CT機(jī)有種莫名其妙的興趣,看到CT拍到的膠片覺(jué)得很神奇,一直想明白其中的原理.
這個(gè)周末莫名其妙的開始在YouTube上看介紹computed tomography的算法,便開始看算法,寫代碼,經(jīng)過(guò)多半天的努力,突然豁然開朗.
對(duì)計(jì)算機(jī)斷層掃描算法也有了了解,通過(guò)手?jǐn)]代碼實(shí)現(xiàn)了computed tomography的功能,并初步了解了Radon transform算法和Reconstruction with the Filtered Back Projection (FBP)算法,當(dāng)然這個(gè)也繞不開傅里葉變換.

0.CT原理

我找了一些動(dòng)畫來(lái)理解CT的工作原理. 這里只是簡(jiǎn)單介紹,主要是為了了解成像算法.跟深入的大家可以到維基百科上查查看.

0ctwork.gif

最終計(jì)算機(jī)通過(guò)斷層掃描技術(shù)形成如下的圖像:

1ct.gif

1.原理

簡(jiǎn)單說(shuō),CT通過(guò)X光旋轉(zhuǎn)傳感器形成斷層信號(hào),步進(jìn)后繼續(xù)采集形成新的斷層信號(hào),(注意這里是信號(hào),還以一維的信號(hào)),然后每層信號(hào)通過(guò)Radon transform算法,形成Sinogram(正弦)圖.
Sinogram通過(guò)FBP(反向重構(gòu)算法)形成斷層圖像.
貌似很簡(jiǎn)單,比大象放入冰箱還少一步.

便于理解我放一張圖:

2ct.png

2.Radon算法

• Radon transform Radon算法實(shí)現(xiàn)切片的Sinogram(正弦)圖.
  在計(jì)算機(jī)斷層掃描中，斷層重建問(wèn)題是從一組投影[1] _獲得斷層攝影切片圖像。 通過(guò)繪制一組穿過(guò)2D拍照區(qū)域的平行光線，將對(duì)象的對(duì)比度沿每個(gè)光線的積分分配給投影中的單個(gè)像素，即可形成投影。 2D對(duì)象的單個(gè)投影是一維的。 為了能夠?qū)?duì)象進(jìn)行計(jì)算機(jī)斷層成像重建，必須獲取多個(gè)投影，每個(gè)投影對(duì)應(yīng)于射線相對(duì)于對(duì)象之間的不同角度。 在幾個(gè)角度上投影的集合稱為正弦圖，它是原始圖像的線性變換。
  為了好理解,我畫了一張中間有個(gè)圓的圖片,模擬CT斷層掃描,生成正弦圖.

import numpy as np
import matplotlib.pyplot as plt

from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage import io,data

#image = shepp_logan_phantom()

#image = rescale(image, scale=0.4, mode='reflect', multichannel=False)
image = io.imread('ct10.jpg')
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))

ax1.set_title("Original")
ax1.imshow(image, cmap=plt.cm.Greys_r)

theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta, circle=True)
ax2.set_title("Radon transform\n(Sinogram)")
ax2.set_xlabel("Projection angle (deg)")
ax2.set_ylabel("Projection position (pixels)")
ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
           extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')

fig.tight_layout()
plt.show()

我又換了幾組不同的圓,大家可以感受到掃描一圈形成的正弦圖和實(shí)物斷層的關(guān)系.

3ct.png

4ct.png

5ct.png

最后我們換一張真正的CT圖片.結(jié)果如下:

6ct.png

最后我們看一張動(dòng)畫了解這個(gè)掃描過(guò)程.

tomo1.gif

FBP算法

濾波后的反投影的數(shù)學(xué)基礎(chǔ)是傅立葉切片定理 。 它使用傅立葉空間中投影和插值的傅立葉變換獲得圖像的2D傅立葉變換，然后將其反轉(zhuǎn)以形成重建圖像。 濾波后的反投影是執(zhí)行Radon逆變換的最快方法之一。 FBP的唯一可調(diào)參數(shù)是濾波器，該濾波器應(yīng)用于傅立葉變換的投影。 它可用于抑制重建中的高頻噪聲。

• CT掃描后獲得的正弦圖

ct02.jpg

• 使用FBP算法
  從正弦圖使用Radon逆變換,獲得切片圖像.

import numpy as np
import matplotlib.pyplot as plt

from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage import io,data

image = io.imread('ct02.jpg')
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta, circle=True)
io.imsave('ct03.jpg',sinogram)

• 獲得切片圖像

ct03.jpg

最后使用動(dòng)畫直觀了解一下Radon逆變換獲取切片圖像的過(guò)程.

tomo2.gif

3.說(shuō)明

這個(gè)筆記只是簡(jiǎn)單把算法走了一遍,有關(guān)Radon和Radon逆變換及傅里葉變換,還是很難理解的,我這里使用了skimage中封裝好的函數(shù).有關(guān)skimage也是機(jī)器學(xué)習(xí)和圖像處理中很重要的庫(kù).