昆明理工大学信息工程与自动化学院学生实验报告（ 2009—2010学年 第 一 学期 ）课程名称：医学成像系统与放射治疗装置 开课实验室： 3208 2008 年 12 月24 日一、实验目的与意义医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及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 )2exp()()(πθθ-=⎰+∞∞- dsdt wt j s t g )2ex p(),(π-=⎰+∞∞- （2.5）将上式变换到),(y x 坐标系中，注意到变换的可比行列式1cos sin sin cos =-=∂∂∂∂∂∂∂∂=⎰θθθθys y tx 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)的近似形式∑=+=Mi i iiy x Q My 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),('''202---⋅=⎰⎰- （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）其离散形式为： ααααββ)()()(1'm n g m R n Q Mm -=∑- （2.30）这里也和上节一样可以加进一光滑窗，改进重建效应。

第三步：反投影，就是执行（2.30）的积分 βγφγβφβπd Q L r f )(),,(1),(202⎰= （2.31）近似的有：)(),,(12),(12γφγβπβQ L my x f Mi ∑==（2.32）其),(y x f 是图象),(y x f 的近似图像，φcos r x =，φsin r y =， M 是投影数，假定投影均匀地分布在0~π2内。