滤波反投影法重建CT 图像实验指导书一、 实验目的1． 了解傅立叶变换法、直接反投影法重建CT 图像的原理；2． 掌握滤波反投影法重建CT 图像的原理和基本方法。
二、 实验器材装有MATLAB 程序的PC 机，滤波反投影法图像重建演示软件，投影数据。
三、 实验原理CT 图像重建问题实际上就是如何从投影数据中解算出成像平面上各像素点的衰减系数。
图像重建的算法有多种，如反投影法、傅立叶变换法、迭代法、滤波反投影法等。
在介绍算法前，有必要先介绍从投影重建图像的重要依据，即中心切片定理。
1． 中心切片定理密度函数(,)f x y 在某一方向上的投影函数()g R θ的一维傅立叶变换函数()g θρ是原密度函数(,)f x y 的二维傅立叶变换函数(,)F ρθ在(,)ρθ平面上沿同一方向且过原点的直线上的值。
图1 中心切片定理2．傅立叶变换法如果在不同角度下取得足够多的投影函数数据，并作傅立叶变换，根据中心切片定理，变换后的数据将充满整个(,)u v 平面。
一旦频域函数(,)F u v 或(,)F ρβ的全部值都得到后，将其做傅立叶反变换，就能得到原始的密度函数(,)f x y ，即所要重建的图像。
上述图像重建算法称为傅立叶变换法，图2给出了傅立叶变换重建方法的流程图。
图中指出，对于每次测得的投影数据先作一维傅立叶变换。
根据中心切片定理，可将此变换结果看成二维频率域中同样角度下过原点的直线上的值。
在不同投影角下所得的一维变换函数可在频域中构成完整的二维傅立叶变换函数，将此二维变换函数做一次逆变换，就得到了所要求的空间域中的密度函数。
为了在二维逆变换中采用快速傅立叶变换算法，通常在逆变换前要将极坐标形式的频域函数变换成直角坐标形式的数据。
图2 傅立叶变换重建图像的过程采用傅立叶变换法重建图像时，投影函数的一维傅立叶变换在频域中为极坐标形式，把极坐标形式的数据通过插补运算转换为直角坐标形式的数据时，计算工作量较大。
此外，在极坐标形式的频域数据中，离原点较远的频率较高的部分数据比较稀疏，当这些位置上的数据转换到直角坐标下时，需经插补，这将引入一定程度的误差。
即，在重建的图像中，高频分量可能有较大失真。
3．直接反投影法直接反投影法把每次测得的投影数据“原路”反投影到投影线的各个像素上。
即指定投影线上所有各点的值等于所测得的投影值。
如图3所示的例子中，被探查物体只是在原点位置上的一个致密点。
在一个角度上测量到其投影值时，就把这个值赋给投影线上的所有点。
于是，从不同角度进行反投影后的重建图像是由以原点为中心的一系列辐射线。
显然，原点位置上的分布密度最高；愈往四周，密度愈低。
这当然可以说是粗略地把图像恢复出来了。
但问题是除了密度最大的中心点，四周出现了逐渐变浅的云晕状阴影。
图3直接反投影重建法研究发现，反投影后重建的密度函数(,)b f x y 与实际的密度函数(,)f x y 满足如下卷积关系 1(,)(,)(1)b f x y f x y r =** 上式说明了直接反投影法造成图像模糊的原因。
一个以δ函数形式出现的密度函数用直接反投影法重建后成了一个带“长尾巴”的1/r 形式的图像，如图4所示。
图4直接反投影重建法造成的伪像为了得到真实的重建密度函数，必须设法消除1/r 的影响。
一种可能的方法是把式(1)转换到频域。
因为时域中两个函数的卷积的傅立叶变换是这两个函数分别的傅立叶变换的乘积，因此式(1)相应的频域表达式为1(,)(,)(2)b F F ρβρβρ=式中(,)b F ρβ和(,)F ρβ分别是(,)b f x y 与(,)f x y 的傅立叶变换，1/ρ为1/r 的傅立叶变换。
从公式(2)可得1(,)[(,)](3)b f x y F ρρβ-=从式(3)可知，为了获得真实的密度函数(,)f x y ，可以先求出反投影函数(,)b f x y 的傅立叶变换(,)b F ρβ，在频域中对(,)b F ρβ加上权重ρ后求其逆傅立叶变换，就能得到所要的密度函数(,)f x y ，如图5所示。
用这样的方法重建图像当然是可行的，但它还是没有避免计算二维傅立叶变换的问题。
两次二维傅立叶变换所花费的时间还是相当可观的。
图5 直接反投影法的图像校正4．滤波反投影法直接反投影法重建的图像是模糊的。
虽然这种模糊的图像可以经过修正后再现原始的密度函数，但修正的过程是很费时的。
这就是先反投影，后修正重建方法存在的问题。
滤波反投影法采用先修正、后反投影的做法，同样可得到原始的密度函数。
其基本方法是：在某一投影角下取得了投影函数（一维函数）后，对此一维投影函数作滤波处理，得到一个经过修正的投影函数；然后再将此修正后的投影函数作反投影运算，得到所需的密度函数。
这一方法中要解决的主要问题是如何修正投影函数才能使之在反投影后能重建原密度函数。
图6给出了滤波反投影法重建图像的流程图。
图6 滤波反投影法滤波反投影法重建图像有以下几个步骤：（1）对某一角度下的投影函数作一维傅立叶变换；（2）对（1）的变换结果乘上一维权重因子 ；（3）对（2）的加权结果作一维逆傅立叶变换；（4）用（3）中得出的修正过的投影函数做直接反投影；（5）改变投影角度，重复（1）~（4）的过程，直到完成全部180度的反投影。
5．滤波函数滤波函数的选取是滤波反投影法的关键问题。
在实际系统中选择滤波函数时还要考虑许多其他因素，包括系统的带宽、信噪比与分辨率等。
下面介绍两中不同类型的滤波函数。
（1）R-L 滤波函数R-L 滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限，所以频域中的滤波函数ρ可表示为 0,()(4)0R L H ρρρρ-⎧≤=⎨⎩其滤波函数如图7所示。
图7 R-L 滤波函数连续的R-L 卷积函数为22000()[2sin (2)sin ()](5)R L h R c R c R ρρρ-=-离散化的R-L 卷积函数为22221,04()0,(6)1,R L n T h nT n n n Tπ-⎧=⎪⎪=⎨⎪⎪-⎩为偶数为奇数图8（a ）（b ）分别为连续形式、离散形式的R-L 卷积函数。
图8 R-L 卷积函数，（a ）连续形式，（b ）离散形式R-L 滤波函数的特点是函数形式简单，重建的图像轮廓清晰。
存在的问题是：由于在频域中用矩形函数截断了滤波函数，在相应的空域中造成振荡响应，即所谓的Gibb’s 现象。
另外，如果投影数据中含有噪声，则重建的图像质量也不够满意。
（2）S-L 滤波函数与R-L 滤波函数不同的是，S-L 滤波函数在频域中不用矩形函数去截断ρ滤波函数，而是用一些比较平滑的窗函数来约束滤波函数。
例如，可以令窗函数()W ρ为00()sin ()()22W c rect ρρρρρ= 于是可得S-L 滤波函数为 00()sin ()()(7)22S L H c rect ρρρρρρ-=其滤波函数如图9所示。
图9 S-L 滤波函数对应的卷积函数为200020414sin(2)1()()(8)21(4)S L R R h R R ρρπρπρ--=-其离散化的卷积函数为 2222()(9)(41)S L h nT T n π-=--图10（a ）（b ）分别为连续形式、离散形式的S-L 卷积函数。
图10 S-L卷积函数，（a）连续形式，（b）离散形式用S-L滤波函数重建的图像中振荡相应较小，对含噪声的数据重建出来的图像质量也较R-L滤波函数重建的图像质量要好。
但是，S-L滤波函数重建的图像在高频响应方面不如R-L滤波函数好，这是因为S-L滤波函数在高频段偏离了理想的滤波函数 。
滤波反投影法无需等待所有数据采集完毕后再作处理。
当扫描系统在作机械运动时，每当在一个角度上获得了投影函数后，计算机马上就可对其进行傅立叶变换等处理。
这样，一旦扫描结束后只要再作数毫秒的处理就可得到重建的图像。
四、实验步骤1．运行MA TLAB程序，打开start_imgreclab.m文件，运行该M文件，显示滤波反投影法图像重建演示软件主界面，如图11所示。
图11 滤波反投影法图像重建演示软件主界面2．鼠标点击主界面中“Open ”按钮，弹出文件选择对话框，选择sinogram1.mat 数据后，显示正弦图，如图12所示。
投影数据是以正弦图的方式存放在sinogram1.mat 数据文件中。
横坐标为投影角，纵坐标为某一投影角下的投影值。
所以正弦图可以看成由不同角度下采集的投影函数值一行行叠放起来的数据集。
Sinogram 12040608010012014016018050100150200250300350400450500图12 sinogram1.mat 数据文件中存放的投影数据3．打开Filter 下面的下拉框，可选择Ram-Lak (Ramp)，Shepp-Logan ，Ram-Lak Cosine ，Ram-Lak Hamming ，Ram-Lak Hann ，Special 滤波器，也可以不用滤波器，选择None 。
4．打开Interpolation 下面的下拉框，可选择不同的插值方式，有Linear （线性）, Nearest （最近邻）,Spline （样条）,Cubic （三次）插值方式。
5．Number of Projections 下面的文本框可输入投影数据的数目，其最大值不能超过180。
6．图13为180个投影角、不用滤波、线性插值重建得到的图像，可以看到该图像非常模糊。
图14为同样条件下，采用Ram-Lak 滤波后重建得到的图像，可见为一清晰的头部断面图像。
因此，滤波反投影重建算法要优于直接反投影重建算法。
Sinogram 1: 180 projections, None filter, Linear interpolation5010015020025030035050100150200250300350图13 180个投影角、不用滤波、线性插值重建得到的图像Sinogram 1: 180 projections, Ram-Lak (Ramp) filter, Linear interpolation5010015020025030035050100150200250300350图14 180个投影角、采用Ram-Lak 滤波、线性插值重建得到的图像7．图15为180个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。
仔细比较图14与图15，可以看到图15比图14略微清晰，原因是三次插值比线性插值效果要好，但三次插值计算量比线性插值要大很多。
Sinogram 1: 180 projections, Ram-Lak Hamming filter, Cubic interpolation5010015020025030035050100150200250300350图15 180个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像8．图16为30个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。