滤波反投影法重建CT图像实验指导书一、实验目的1.了解傅立叶.docx

滤波反投影法重建CT图像实验指导书
一、实验目的
1.了解傅立叶变换法、直接反投影法重建CT图像的原理;
2.掌握滤波反投影法重建CT图像的原理和基木方法。

二、实验器材
装有MATLAB程序的PC机，滤波反投影法图像重建演示软件，投影数据。

三、实验原理
CT图像重建问题实际上就是如何从投影数据中解算出成像平面上各像索点的衰减系数。

图像重建的算法有多种，如反投影法、傅立叶变换法、迭代法、滤波反投彩法等。

在介绍算法前，有必要先介绍从投影重建图像的重要依据，即中心切片定理。

1.中心切片定理
密度函数/(x,)0在某一•方向上的投影函数g&(/?)的一•维傅立叶变换函数g&(p)是原密度两数f(x,y)的二维傅立叶变换函数F(p,&)在(Q,&)平面上沿同一•方向且过原点的直线上的值。

图1屮心切片定理
2.傅立叶变换法
如果在不同角度下取得足够多的投影函数数据，并作傅立叶变换，根据小心切片定理，变换后
的数据将充满整个(仏叭平而。

一旦频域函数F(w,v)或F(/?,0)的全部值都
得到后，将真做傅立叶反变换，就能得到原始的密度函数.f(x,y)，即所要重建的图像。

上述图像重建算法称为傅立叶变换法，图2给出了傅立叶变换重建方法的流程图。

图屮指出，对于每次测得的投影数据先作一维傅立叶变换。

根据中心切片定理，可将此变换结果看成二维频率域屮同样角度下过原点的直线上的值。

在不同投影角卜•所得的一维变换函数可在频域中构成完整的二维傅立叶变换函数，将此二维变换函数做一次逆变换，就得到了所要求的空间域屮的密度函数。

为了在二维逆变换屮采用快速傅立叶变换算法，通常在逆变换前要将极坐标形式的频域函数变换成直角坐标形式的数据。

图2傅立叶变换重建图像的过程
采用傅立叶变换法重建图像吋，投影函数的一维傅立叶变换在频域中为极处标形式, 把极坐标形式的数据通过插补运算转换为直角坐标形式的数据时，计算工作量较大。

此外, 在极坐标形式的频域数据屮，离原点较远的频率较高的部分数据比较稀疏，当这些位登上的数据转换到肯角坐标下时，需经插补，这将引入一定程度的谋差。

即，在重建的图像屮, 高频分虽可能冇较大火真。

3.直接反投影法
总接反投影法把每次测得的投影数据“原路”反投彫到投影线的各个像素上。

即指定投影线上所冇各点的值等于所测得的投影值。

如图3所示的例了中，被探查物体只是在原点位置上的一个致密点。

在一个角度上测最到其投彩值时，就把这个值赋给投彩线上的所侑点。

于是，从不同角度进行反投影后的重建图像是由以原点为中心的一系列辐射线。

显然，原点位置上的分布密度最高；愈往四周，密度愈低。

这当然对以说是粗略地把图像恢复出来了。

但问题是除了密度最人的屮心点，四周出现了逐渐变浅的云晕状阴影。

图3直接反投影重建法
研究发现，反投影后重建的密度函数f b(x,y)与实际的密度函数/(兀』)满足如下卷积
关系
东(兀刃= /Uy)**~ (1)
r
上式说明了直接反投影法造成图像模糊的原因。

一个以》函数形式出现的密度函数用肓接反投影法重建后成了一个带“长尾巴”的1/厂形式的图像，如图4所示。

图4宜接反投影重建法造成的伪像
为了得到真实的重建密度函数，必须设法消除1/广的影响。

一种町能的方法是把式(1)转换到频域。

因为时域屮两个两数的卷积的傅立叶变换是这两个函数分别的傅立叶变换的乘积, 因此式(1)相应的频域表达式为
休(p,0) = F 9,0)丄(2)
P
式中F b(p, J3)和F(p,0)分别是f h(x,y)与f(x,y)的傅立叶变换,1/p为1/厂的傅立叶变换。

从公式(2)可得
/(x,y) = L[0F；(p,0)] (3)
从式(3)可知，为了获得真实的密度函数/(x,y),可以先求出反投影函数人(兀,y)的傅立叶变换％9，0),在频域中对耳(00)加上权重Q后求其逆傅立叶变换，就能得到所要的密度函数/(兀，刃，如图5所示。

用这样的方法重建图像当然是可行的，但它还是没有避免计算二维傅立叶变换的问
题。

两次二维傅立叶变换所花费的时间述是相当可观的。

直接
反投影
2D
F(Q,0)
图5直接反投影法的图像校正
4.滤波反投影法
直接反投影法重建的图像是模糊的。

虽然这种模糊的图像可以经过修正后再现原始的密度函数，但修正的过程是很费时的。

这就是先反投影，后修疋重建方法存在的问题。

滤波反投彩法采用先修正、后反投影的做法，同样可得到原始的密度函数。

其基木方法是：在某一投影角下取得了投影函数（一维函数）后，对此一维投影函数作滤波处理，得到一个经过修正的投影函数；然后再将此修正后的投影函数作反投影运算，得到所需的密度函数。

这一方法中要解决的主要问题是如何修止投影函数才能使Z在反投影后能重建原密度函数。

图6 给出了滤波反投影法重建图像的流程图。

滤波反投影法重建图像有以卞儿个步骤：
（1）对某一角度卜•的投影函数作一维傅立叶变换;
（2）对（1）的变换结果乘上一维权重因了|°|；
（3）对（2）的加权结果作一维逆傅立叶变换；
（4）用（3）中得出的修正过的投影函数做直接反投影；
（5）改变投影角度，重复（1）〜（4）的过程，直到完成全部180度的反投影。

5•滤波函数
滤波函数的选取是滤波反投影法的关键问题。

在实际系统中选择滤波函数时还要考虑 许多其他
因素，包括系统的带宽、信噪比与分辨率等。

下面介绍两屮不同类型的滤波函数。

(1) R-L 滤波函数
R-L 滤波函数的基木出发点是认为实际的二维图像函数总有一个频率上限，所以频域屮 的滤波函数制可表示为
\p\ § A ) 0
其滤波函数如图7所示。

连续的R-L 卷积函数为
h R _L (R)=犹[2 sin c(2p Q R) - sin c 2 (p 0/?)]
(5)
离散化的R-L 卷积函数为
1 4T 7
為2严，汹奇数
图
8 (a) (b)分别为连续形式、离散形式的R 丄卷积函数。

图8 R 丄卷积函数，Q)连续形式，(b)离散形式
(〃)=
0, 斤为偶数
R ・L 滤波函数的特点是函数形式简单，重建的图像轮廓清晰。

存在的问题是：由于在频 域中用短形函数截断了滤波函数，在相应的空域中造成振荡响应，即所谓的Gibb ，s 现象。

另外，如果投彫数据中含有噪声，则重建的图像质量也不够满意。

(2) S 丄滤波函数
• R -L 滤波函数不同的是，S-L 滤波函数在频域屮不用矩形函数公截断|"|滤波函数， 而是用一些比较平滑的窗函数來约束滤波函数。

例如，可以令窗函数W(Q )为
于是可得S ・L 滤波两数为
其滤波函数如图9所示。

对应的卷积函数为
其离散化的卷积函数为
图10 (a) (b)分别为连续形式、离散形式的S ・L 卷积函数。

H —(P )Tp|sin
(7)
4炖)2「4几肌皿2切/?)
7T
1一(4炖貯
(M)
2 ^•2r 2(4n 2 -1)
W(p) = sin
rect
用S-L滤波函数重建的图像屮振荡相丿应较小，对含噪声的数据重建出来的图像质暈也较R-L滤波函数重建的图像质量要好。

但是,S-L滤波函数重建的图像在高频响应方面不如R-L 滤波函数好，这是因为S-L滤波函数在高频段偏离了理想的滤波函数|p| o
滤波反投影法无需等待所冇数据采集完毕后再作处理。

当扫描系统在作机械运动时，每当在一个角度上获得了投影函数后，计算机马上就可对具进行傅立叶变换等处理。

这样，一旦扫描结束示只要再作数毫秒的处理就可得到重建的图像。

四.实验步龙
1.运行MATLAB程序，打开stan_imgreclab.m文件，运行该M文件，显示滤波反投影法图像重建演示软件主界而，如图11所示。

图10 S-L卷积函数，（a）连续形式，（b）离散形式
力S
(a) (b) 图11滤波反投影法图像重建演示软件主界面
2.鼠标点击主界血屮“Open”按钮，弹出文件选择对话框，选择sinogram 1.mat数据后，显示正弦图，如图12所示。

投影数据是以正弦图的方式存放在sinogram I.mat数据文件中。

横坐标为投影角，纵坐标为某一投影角下的投影值。

所以匸弦图可以看成由不同角度卜•采集的投影函数值一行行蒂放起来的数据集。

Sirwyar 1
50
100
150
200 >
250
300-
350
400
450
500
20406080100120160180
图12 sinogram 1.mat数据文件屮存放的投影数据
3.打开Filter 下面的下拉框,可选择Ram-Lak （Ramp）, Shepp-Logan, Ram-Lak Cosine, Ram-Lak Hamming, Ram-Lak Hann, Special 滤波器，也町以不川滤波器，选择None。

4.打开Interpolation F面的下拉框,可选择不同的插值方式,有Linear （线性）,Nearest （最近邻）,Spline （样条）,Cubic （三次）插值方式。

5.Number of Projections Fifil的文本框可输入投影数据的数目，其最大值不能超过180。

6•图13为18（）个投影角、不用滤波、线性插值重建得到的图像，可以看到该图像非常模糊。

图14为同样条件下，采用Ram-Lak滤波后重建得到的图像，可见为一清晰的头部断面图像。

因此，滤波反投影重建算法要优于直接反投影重建算法。

50
100
150
200
250
300
350
图13 180个投影
角、不用滤波、线性插值重建得到的图像
图14 18()个投影角、采用Ram-Lak 滤波、线性插值重建得到的图像
7. 图15为180个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。

仔细 比较图14与图15,可以看到图15比图14略微清晰，原因是三次插值比线性插值效果要好, 但三次插值计算罐比线性插值要人很多。

50 100 150 200 250 300 350
Sinogram 1: 180 projectio ns, Ram-Lak (Ramp) filter, Lin ear i nterpolati on
50
100
150
200
250
300
350
50
100
150
200
250
300
350
图15 180个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像
8. 图16为3()个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。

与图 15相比，可见投影角度数目，即投影数据较少时，图像有线状条纹，比较模糊。

图16 30个投彩角、采丿IJ Ram-Lak Hamming 滤波、三次插值重建得到的图像
9. 分别选择不同的投影数据，选用不同的滤波器、插值函数及投影角数冃重建图像，并比 较重建图像的差异。

1().在MATLAB 命令窗中输入“type iradon 命令，窗口显示逆雷登变换的MATLAB 程序,
100
150
200
250
300
350
50 100 150 200 250 300 350
Sinogram 1: 30 projections, Ram-Lak Hamming filter, Cubic interpolation
50
100
150
200
250
300
350
50 100 150 200 250 300 350
读懂该程序，在关键的代码后用中文注释。

11.画出程序屮丿IJ到的Ram-Lak Cosine, Ram-Lak Hamming, Ram-Lak Hann 滤波器的滤波函数及其港积函数。

五、实验结果
1.采用高阶插值函数其优点是重建图像质量好，其缺点是_____________________________ 0
2.如果重建后的图像细节很多，即高频分聚多，则该选用_________________ 滤波函数较好。

六、结果讨论与思考题
1.为什么投彫角度数目少，重建图像不清晰？
2.为什么滤波反投影重建图像耍比直接反投影重建图像要清晰得多?。