CT图像重建
昆明理工大学信息工程与自动化学院学生实验报告
(2009—2010学年第一学期)
一、实验目的与意义
医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及X光仪器,CT仪器,MRI 仪器和核医学仪器的工作原理及成像方法。其中CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。
CT技术是医学成像系统中的一种重要手段。它通过特定的算法,利用计算机的高速运算功能,可以在短时间内快速呈现人体断层图像。让学生练习CT图像的重建有助于学生理解CT算法的内容,熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。
二、实验算法原理
1、MATLAB处理数字图像的基本函数;
2、X-CT三维图像重建的基本算法。
CT图象重建有四种基本的算法:矩阵法,迭代法,傅立叶算法,反投影算法.我们采用的方法为卷积反投影. 卷积反投影有:平行光束投影的卷积反投影算法, 等角扇形光来投影的重建算法.
1).平行光束投影的卷积反投影算法
从投影重建三维物体的图像,就是重建一个个横断面。这样三堆图像的重建就归结为二维图象的重建。二维图像的重建问题可以从数学上描述如下。
假定),(y x g 表示一个二维的未知函数,通过),(y x g 的直线称为光钱(见图2.1)。沿光线),(y x g 的积分称作光线积分。沿相同方向的一组光线积分,就构成一个投影。图2.1中垂直于直线'
CC (与X 轴夹角为θ)的光线所形成。
图2.1 ),(y x g 在θ方向的投影)(t P θ
的投影)(t P θ,称之为),(y x g 在θ方向的投影。光线积分和投影在数学上可以定义如下:
在图2.1中直线AB 的方程为:
1sin cos t Y X =+θθ (2.1) 其中1t 是AB 到原点的距离,),(y x g 沿AB 的积分为:
dxdy t y x y x g ds y x g t P AB
)sin cos (),(),()(11-+==??+∞
∞
-θθδθ (2.2)
对于给定的θ,),(y x g 在θ方向的投影)(t P θ是t 的函数。如果),(y x g 在各个方向的投影已知,),(y x g 就可以唯一确定。下面就讨论卷积反投影重建算法。
假定投影方向θ,如图2.2,将坐标),(y x 旋转θ角(逆时针方向)形成坐标系),(s t 。),(y x g 在),(s t 坐标系中为),(s t g 。
图2.2 傅立叶切片定理示意图
坐标系),(s t 与),(y x 之间的关系为:
???
? ?????? ?
?-=???? ??y x s t θθθθcos sin sin cos (2.3)
显然
()ds s t g t P ?
+∞
∞
-=
),(θ (2.4)
令)(w S θ为)(t P θ的傅立叶变换则 dt wt j t P w S )2ex
p()()(πθθ-=?
+∞
∞- dsdt wt j s t g )2ex p(),(π-=
?
+∞
∞
- (2.5)
将上式变换到),(y x 坐标系中,注意到变换的可比行列式
1cos sin sin cos =-=????????=?θθθθy
s y t
x s t t (2.6) 从而得到:
dxdy y x j y x g w S ?+-=??
+∞∞-+∞
∞-)]sin cos (2ex p[),()(θθπωθ
dxdy vy ux j y x g )](2ex p[),(+-=??
+∞∞-+∞
∞
-π (2.7)
其中
??
?==θ
ωθ
ωsin cos v u (2.8)
若令),(y x g 的傅立叶变换为),(v u G ,由(2.8)可知
),(),()(θωωθG v u G S == (2.9) 若),(y x g 的傅立叶变换为),(v u G 的极坐标表示。这说明),(y x g 在θ方向的投影)(t P θ
傅立叶变换)(w S θ等于),(v u G 在与u 轴成θ角的直线上的值。这就是著名的傅立叶投影切片定理。可见在整个),(v u 平面),(v u G 可以利用各个方向的投影来得到,从而),(y x g 也可以通过求),(v u G 的傅立叶反交换的办法求得: dudv vy ux j v u G y x g )](2ex p[),(),(+=??
+∞∞-+∞
∞
-π (2.10)
变换到极坐标中 ???==θ
ωθ
ωsin cos v u , ω=?
得到
θωθθπωθωπd d y x j G y x g )]sin cos (2ex p[),(),(20
+=??
∞
(2.11)
经推导得 ?
???
????=+∞∞-π
θθωπωωω0
)2exp()(),(d d t j S y x g (2.12) 其中
θθsin cos y x t += 若令
ωπωωωθθd t j S t Q )2ex p()()(?+∞
∞
-= (2.13)
则
?
+=
π
θθθθ0
)sin cos (),(d y x Q y x g (2.14)
(2.13)式右端是两频谱函数)(w S θ和)(ωH 的乘积的傅立叶反变换。)(w S θ是投影)(t P θ 傅立叶变换。若)(ωH 的傅立叶反变换为)(t h ,则根据卷积定理有: ?
+∞
∞
--=τττθθd t h P t Q )()()( (2.15)
或
)()()(t h t P t Q *=θθ 其中
ωπωωd t j t h )2ex p()(?
+∞
∞
-=
(2.16)
当图像的频谱是有限带宽时,则上式变为 ωπωωωω
d t j t h )2ex p()(0
?+-=
(2.17)
由于图象及其频谱都是离散采样的, 假定图象采样间隔为τ, 则根据采样定理τω2/10=。为了进行数学处理,只需知道h (t)在有限带宽上的离散采样点的值.这样我们有
??
?
??-=2222/10)4/(1)(τπττn n h
(2.18)
其中n 为正负整数。 (2.18)的离散形式为 ∑∞
-∞
=-=m m n h m P n Q τ
ττ
τθθ)()()( (2.19)
假定)(τθm P 在1,1,0-??????=N m 之外的值为0,则上式变为 []∑-=-=1
)()()(N m m n h m P n Q ττττθθ (2.20)
或
∑---=-=1
)
1()(][)(N N m m h m n P n Q τττ
τθθ (2.21)
其中1,2,1,0-??????=N n 从而可见为确定)(t P θ的N 个采样点上的)(τθn Q 的值,需要使用)(τn h 的2N — 1个点上的值,从n=一(N — 1)到(N — 1)。
为求得)(τθn Q ,利用傅立叶变换计算卷积是比较快的方法,为清除循环卷积的周期交叠效应,实际上)(τn h 取2N 个点,)(τθm P 补0,使之有(2N —1)个元素,则)(t P θ在N 个采样点上就避免了交叠,如果使用以2为基的FFT(快速傅立叶变换)算法, )(τθm P 和)(τn h 都必需朴0至(2N 一1)个元素,(2N 一1)为大于等于2N —l 的最小的2的整数幂。计算)(τθn Q 的过程可以写为
]0)((]0)([([)(ττττθθn h FFT n P FFT FFT n Q ??= (2.22)
其中FFT 和IFFT 分别表示快速傅立叶变换和反变换, 光滑窗是在滤波过程中加入的光滑因子,例如引用汉明窗 ,有时可以改进重建效果。对于各个θ方向的投影, 得到)(τθn Q 之后就可以由(2.22)来计算
),(y x g 。重建步骤可以归纳为:
第一步:卷积,也称滤波,由(2.22)对每个θ方向计算)(τθn Q 。 第二步。反投影,由(2.14)的近似形式
∑=+=
M
i i i
i
y x Q M
y x g 1
)sin cos (),(θθ
π
θ (2.23)
来计算),(y x g 的近似值),(?y x g
。M 为投影个数i θ为投影方向角,他们均匀的分布在0~π的范围内。 当计算)sin cos (i i y x Q i θθθ+时,i i y x t θθsin cos +=,不一定在)(τθn Q 的整离散点上,这就需要插值求得,预先将)(τθn Q 插值加密,即最靠近的点,可以提高计算速度。 2).等角扇形光来投影的重建算法
几乎所有的快遗CT 设备都是用的扇形光束。这里叙述的是等角度光束投影,如图2.3,测量投影数据的探测器等间距地分布在1D 2D 弧上,弧的半径为2D , D 为光源到图像中心的距离。在下文中,),(φr f 图象在极坐标中的表示。)(γβR 表示在方向角为β的投影中位量角为γ的光线产生的投影数据。通过中心的光线其γ为0。L 表示从光源到像素),(φr 的距离。
图2.3 等扇形束投影重建算法中的变量
)sin(2),,(22φβφγβ-++=Dr r D L (2.24)
γ表示在方向角为β的投影中通过像素),(φr 的光线的位置角
)
sin()
cos(tan tan
),,(11
'
φβγφβγφγβγ-+-==--D E S E P (2.25) 图像),(φr f 和扇形投影)(γβR 有下述关系
βγγγγγγγγφγγβπ
d h rR D L r f )())
sin((21)(cos 1),('
'
'20
2---?=??
- (2.26) 其中r 和φ是投影中光线的最大位置角,从而可以得到这种重建算法的执行步骤: 第一步:投影的修改,假定投影的抽样间隔为α,抽样数据)(αβn R 通过下式进行修改, αααββn D n R n R cos )()(1'
= (2.27)
第二步:卷积(滤波)将第一步修改了的投影与响应函数)(γg 进行卷积 )()sin 1(21)(2
'
γγ
γh g =
(2.28) αααγγγ
ββd m n g n R Q )()()(-=?- (2.29)
其离散形式为: