CT 3 图像重建
1. 基本术语
1.1 投影(projection)
投影又称为“射线和”或者线积分
p
(
s
,
θ
)
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
δ
(
x
c
o
s
θ
+
y
s
i
n
θ
−
s
)
d
x
d
y
p(s,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(xcos\theta \space+\space ysin\theta -s)dxdy
p(s,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ + ysinθ−s)dxdy
f
(
x
,
y
)
f(x,y)
f(x,y)为平面x-y上定义的密度函数
δ
\delta
δ为狄拉克函数
\quad
1.2 弦图(sinogram)
用于描述投影的方法之一 弦图就是在重建之前采集一层CT的投影分布图，而弦图对应的空间称为弦空间 弦空间的纵轴代表投影角度，横轴表示每个投影的射线（探测器单元）
弦图可用于检测CT系统是否正常 当某一探测器单元失效，弦图内则含有一条竖直的直线
\quad
1.3 Radon变换
揭示了函数和投影的关系，假定函数为
f
(
x
,
y
)
f(x,y)
f(x,y)，并对其在不同的角度取投影，则
p
(
t
,
θ
)
p(t,\theta)
p(t,θ)和
f
(
x
,
y
)
f(x,y)
f(x,y)的关系式为：
p
(
t
,
θ
)
=
∫
−
∞
∞
∫
−
∞
∞
f
(
x
,
y
)
δ
(
x
c
o
s
θ
+
y
s
i
n
θ
−
t
)
d
x
d
y
p(t,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(xcos\theta \space+\space ysin\theta -t)dxdy
p(t,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ + ysinθ−t)dxdy 逆变换：
f
(
x
,
y
)
=
1
2
π
2
∫
0
π
∫
−
∞
∞
1
x
2
+
y
2
c
o
s
(
t
g
−
1
y
x
−
θ
)
−
t
∂
p
(
t
,
θ
)
∂
t
d
t
d
θ
f(x,y)=\frac{1}{2\pi^2}\int_0^\pi\int_{-\infty}^{\infty}\frac{1}{\sqrt{x^2+y^2}cos(tg^{-1}\frac{y}{x}-\theta)-t}\frac{\partial p(t,\theta)}{\partial t}dtd\theta
f(x,y)=2π21∫0π∫−∞∞x2+y2
cos(tg−1xy−θ)−t1∂t∂p(t,θ)dtdθ
也就是说，一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定
\quad
2 平行束CT图像经典重建算法
由弦图求出CT图像
2.1 直接矩阵求解法
在实际中直接矩阵求解法很难实现
\quad
2.2 迭代法
迭代法的基本步骤：
给出初始矩阵用初始矩阵形成投影把待重建物体投影与模拟投影进行比较若误差满足要求，迭代停止
迭代重建算法的例子可以点击这里参考，出自公众号名为“CT技术共享”，作者为张琼阁老师。 若连接打不开，可在公众号内搜索“CT图像的重建算法（三）迭代重建”
比较精确，但速度慢，多用于核医学设备的图像重建和低剂量CT的图像重建
\quad
2.3 傅里叶变换重建算法
2.3.1 中心切片定律
又称为傅里叶中心切片定理，投影切片定理 一个物体的1D投影的傅立叶变换精确地等于物体2D傅立叶变换在同一角度的直线。 换言之：图像沿某一方向的投影，经过1D 傅立叶变换之后，对应2D傅立叶变换平面的一条线
\quad
2.3.2 傅里叶重建方法
收集CT扫描个角度投影每一投影都计算1DFT规整2D坐标FT平面通过2D反FT算回原图像
\quad
2.3.3 傅里叶重建法的局限性
2D频域上的点不是成矩阵排列的，在做傅立叶逆变换之前需将样本插值转换为笛卡尔坐标（直角坐标）表示频域上
F
(
μ
,
ν
)
F(\mu,\nu)
F(μ,ν) 中某一元素值的改变将导致整幅图像强度的改变，同时还产生明显的阴影伪影难以实现目标重建，逆傅立叶变换的尺寸反比于ROI的尺寸，对于很小的ROI，矩阵太大难以处理对断层的投影作正交变换是一维的，但在求物体图像的逆变换却是二维的，因此，必须将数据都存储起来，等到全部数据完整之后才能进行二维逆变换，
这就要求 硬件内存大，等待的时间长，难于实现实时的图像重建要求
\quad
\quad
2.4 反投影法(Back projection)
由投影重建图像的算法很多，而反投影法是其中最简单、最粗略，也是最基本的算法
原理：“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和（的平均值）”
最后的运算：从每个图像单元中减去背景值（背景强度等于某投射角情况下各投影值之和），再将各吸收系数除以最大公约数，得到最后结果。
其他资料的最后运算：在求出累加值后，再给累加得到的各个像素除以反投影的次数，也就是除以经过像素的射线数。
反投影重建的本质是把取自有限物体空间的射线投影均匀地回抹(反投影)到射线所及的无限空间的各点之上，包括原先像素值为零的点，
→
\rightarrow
→形状伪影
反投影法有时也称为叠加法或总和法。它实现简单而不需要很复杂的数学运算，因此计算速度比较快
但图像边缘会产生模糊
\quad
\quad
2.5 滤波反投影法(FBP)
2.5.1 反投影法和滤波反投影法的区别：
滤波运算或卷积运算的引入（所以FBP也被称为卷积方法） 在反投影前先滤波，即先对1D投影进行滤波，再进行反投影
2.5.2 FBP公式
\quad
P
(
ω
,
θ
)
P(\omega,\theta)
P(ω,θ)表示对应于
θ
\theta
θ角度的单位投影的傅里叶变换
P
(
ω
,
θ
)
∣
ω
∣
P(\omega,\theta)|\omega|
P(ω,θ)∣ω∣ 在空间域表示单位投影被一频域响应为
∣
ω
∣
|\omega|
∣ω∣ 的函数做滤波运算
\quad
2.5.3
高分辨kernels具有高空间分辨力，但噪声大
Shepp-Logan滤波器 平滑图像，损失了部分高频信息
Hamming滤波器 降低了高频噪声
骨滤过器和软组织滤过器 根据诊断需求可选用不同的滤波函数
平滑用于观察软组织 锐利用于观察高分辨力影像
2.5.4 补零问题
原始滤波运算包含一个非周期卷积运算,变到频域后就是周期卷积，直接计算将产竹干涉伪影，即所谓的warp-around效应或interperiod interference。 因此必须在傅立叶变换和滤波操作之前给每一个投影补0，才能避免伪影产生。 0的数目不小于原始投影的样本数减去1(N-1)
FBP实现步骤（平行束）
\quad
\space
3 算法比较
\quad
\space
算法步骤特点迭代法给出初始矩阵用初始矩阵形成投影把待重建物体投影与模拟投影进行比较 若误差满足要求，迭代停止计算精确;耗时长; 改进的迭代算法;是图像重建的一个热点傅里叶变换重建算法对每次测得的投影数据先作1D傅里叶变换在不同投影角度下所得的一维变换函数可在频域中构成完整的二维傅里叶变换函数将此二维变换函数进行逆变换，就得到了所要求的空间域中的密度函数要插值，计算量大；高频部分可能会有明显的失真BP将所测得的投影值按其原路径平均地分配到每一点上;各个方向上投影值反投影后，在影像处进行叠加，从而推断出原图像中心处吸收系数值最大，离中心越远值越低，产生图像的边缘失锐;反投影法会造成影像边缘的不清晰FBP在某一投影角下取得投影函数（一维函数）后，对其作滤波处理再将此修正后的投影函数作反投影运算，得出所需的密度函数只需作一维的傅里叶变换。由于避免了费时的二维傅里叶变换，滤波反投影法明显地缩短了图像重建的时间