计算机断层扫描算法 (computed tomography)

很早我就对CT机有种莫名其妙的兴趣,看到CT拍到的胶片觉得很神奇,一直想明白其中的原理.
这个周末莫名其妙的开始在YouTube上看介绍computed tomography的算法,便开始看算法,写代码,经过多半天的努力,突然豁然开朗.
对计算机断层扫描算法也有了了解,通过手撸代码实现了computed tomography的功能,并初步了解了Radon transform算法和Reconstruction with the Filtered Back Projection (FBP)算法,当然这个也绕不开傅里叶变换.

0.CT原理

我找了一些动画来理解CT的工作原理. 这里只是简单介绍,主要是为了了解成像算法.跟深入的大家可以到维基百科上查查看.

0ctwork.gif

最终计算机通过断层扫描技术形成如下的图像:

1ct.gif

1.原理

简单说,CT通过X光旋转传感器形成断层信号,步进后继续采集形成新的断层信号,(注意这里是信号,还以一维的信号),然后每层信号通过Radon transform算法,形成Sinogram(正弦)图.
Sinogram通过FBP(反向重构算法)形成断层图像.
貌似很简单,比大象放入冰箱还少一步.

便于理解我放一张图:

2ct.png

2.Radon算法

• Radon transform Radon算法实现切片的Sinogram(正弦)图.
  在计算机断层扫描中，断层重建问题是从一组投影[1] _获得断层摄影切片图像。 通过绘制一组穿过2D拍照区域的平行光线，将对象的对比度沿每个光线的积分分配给投影中的单个像素，即可形成投影。 2D对象的单个投影是一维的。 为了能够对对象进行计算机断层成像重建，必须获取多个投影，每个投影对应于射线相对于对象之间的不同角度。 在几个角度上投影的集合称为正弦图，它是原始图像的线性变换。
  为了好理解,我画了一张中间有个圆的图片,模拟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()

我又换了几组不同的圆,大家可以感受到扫描一圈形成的正弦图和实物断层的关系.

3ct.png

4ct.png

5ct.png

最后我们换一张真正的CT图片.结果如下:

6ct.png

最后我们看一张动画了解这个扫描过程.

tomo1.gif

FBP算法

滤波后的反投影的数学基础是傅立叶切片定理 。 它使用傅立叶空间中投影和插值的傅立叶变换获得图像的2D傅立叶变换，然后将其反转以形成重建图像。 滤波后的反投影是执行Radon逆变换的最快方法之一。 FBP的唯一可调参数是滤波器，该滤波器应用于傅立叶变换的投影。 它可用于抑制重建中的高频噪声。

• 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

最后使用动画直观了解一下Radon逆变换获取切片图像的过程.

tomo2.gif

3.说明

这个笔记只是简单把算法走了一遍,有关Radon和Radon逆变换及傅里叶变换,还是很难理解的,我这里使用了skimage中封装好的函数.有关skimage也是机器学习和图像处理中很重要的库.