计算机断层扫描算法 (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也是机器学习和图像处理中很重要的库.
网友评论