滤波反投影法重建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直接反投影重建法
研究发现，反投影后重建的密度函数f b (x , y ) 与实际的密度函数f (x , y ) 满足如下卷积关系
f b (x , y ) =f (x , y ) **
1r
(1)
上式说明了直接反投影法造成图像模糊的原因。一个以δ函数形式出现的密度函数用直接反投影法重建后成了一个带“长尾巴”的1/r 形式的图像，如图4所示。
图4直接反投影重建法造成的伪像
为了得到真实的重建密度函数，必须设法消除1/r 的影响。一种可能的方法是把式(1)转换到频域。因为时域中两个函数的卷积的傅立叶变换是这两个函数分别的傅立叶变换的乘积，因此式(1)相应的频域表达式为
F b (ρ, β) =F (ρ, β)
1
ρ
(2)
式中F b (ρ, β) 和F (ρ, β) 分别是f b (x , y ) 与f (x , y ) 的傅立叶变换，1/ρ为1/r 的傅立叶变换。从公式(2)可得
f (x , y ) =F -1[ρF b (ρ, β)]
-
(3)
从式(3)可知，为了获得真实的密度函数f (x , y ) ，可以先求出反投影函数f b (x , y ) 的傅立叶变换F b (ρ, β) ，在频域中对F b (ρ, β) 加上权重ρ后求其逆傅立叶变换，就能得到所要的密度函数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 滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限，所以频域中的滤波函数
ρ可表示为
⎧ρ,
H R -L (ρ) =⎨
⎩
ρ≤ρ0
(4)
其滤波函数如图7所示。
图7 R-L 滤波函数
连续的R-L 卷积函数为
2
h R -L (R ) =ρ0[2sinc (2ρ0R ) -sin c 2(ρ0R )]
(5)
离散化的R-L 卷积函数为
⎧1
n =0⎪4T 2,
⎪
h R -L (nT ) =⎨0, n 为偶数
⎪1
⎪-222, n 为奇数⎩πn T
图8（a ）（b ）分别为连续形式、离散形式的R-L 卷积函数。
(6)
图8 R-L 卷积函数，（a ）连续形式，（b ）离散形式
R-L 滤波函数的特点是函数形式简单，重建的图像轮廓清晰。存在的问题是：由于在频域中用矩形函数截断了滤波函数，在相应的空域中造成振荡响应，即所谓的Gibb’s现象。另外，如果投影数据中含有噪声，则重建的图像质量也不够满意。
（2）S-L 滤波函数
与R-L 滤波函数不同的是，S-L 滤波函数在频域中不用矩形函数去截断
ρ滤波函数，
而是用一些比较平滑的窗函数来约束滤波函数。例如，可以令窗函数W (ρ) 为
W (ρ) =sin c (
于是可得S-L 滤波函数为
ρρ) rect () 2ρ02ρ0
H S -L (ρ) =ρsin c (
其滤波函数如图9所示。
ρρ) rect () 2ρ02ρ0
(7)
图9 S-L 滤波函数
对应的卷积函数为
14ρ1-4ρ0R sin(2πρ0R )
h S -L (R ) =(0) 2
2π1-(4ρ0R ) 2
其离散化的卷积函数为
(8)
h S -L (nT ) =-
2
π2T 2(4n 2-1)
(9)
图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 1
图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 interpolation
50
100
150
200
250
300
350
图13 180个投影角、不用滤波、线性插值重建得到的图像
Sinogram 1: 180 projections, Ram-Lak (Ramp) filter, Linear interpolation
[***********]50
50
100
150
200
250
300
350
图14 180个投影角、采用Ram-Lak 滤波、线性插值重建得到的图像
7．图15为180个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像。仔细比较图14与图15，可以看到图15比图14略微清晰，原因是三次插值比线性插值效果要好，但三次插值计算量比线性插值要大很多。
Sinogram 1: 180 projections, Ram-Lak Hamming filter, Cubic interpolation
50
100
150
200
250
300
350
图15 180个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像
8．图16为30个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。与图15相比，可见投影角度数目，即投影数据较少时，图像有线状条纹，比较模糊。
Sinogram 1: 30 projections, Ram-Lak Hamming filter, Cubic interpolation
50
100
150
200
250
300
350
图16 30个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像
9．分别选择不同的投影数据，选用不同的滤波器、插值函数及投影角数目重建图像，并比较重建图像的差异。
10．在MATLAB 命令窗中输入“type iradon”命令，窗口显示逆雷登变换的MATLAB 程序，
读懂该程序，在关键的代码后用中文注释。
11．画出程序中用到的Ram-Lak Cosine，Ram-Lak Hamming，Ram-Lak Hann滤波器的滤波函数及其卷积函数。
五、实验结果
1．采用高阶插值函数其优点是重建图像质量好，其缺点是
2．如果重建后的图像细节很多，即高频分量多，则该选用
六、结果讨论与思考题
1．为什么投影角度数目少，重建图像不清晰？
2．为什么滤波反投影重建图像要比直接反投影重建图像要清晰得多？
滤波反投影法重建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直接反投影重建法
研究发现，反投影后重建的密度函数f b (x , y ) 与实际的密度函数f (x , y ) 满足如下卷积关系
f b (x , y ) =f (x , y ) **
1r
(1)
上式说明了直接反投影法造成图像模糊的原因。一个以δ函数形式出现的密度函数用直接反投影法重建后成了一个带“长尾巴”的1/r 形式的图像，如图4所示。
图4直接反投影重建法造成的伪像
为了得到真实的重建密度函数，必须设法消除1/r 的影响。一种可能的方法是把式(1)转换到频域。因为时域中两个函数的卷积的傅立叶变换是这两个函数分别的傅立叶变换的乘积，因此式(1)相应的频域表达式为
F b (ρ, β) =F (ρ, β)
1
ρ
(2)
式中F b (ρ, β) 和F (ρ, β) 分别是f b (x , y ) 与f (x , y ) 的傅立叶变换，1/ρ为1/r 的傅立叶变换。从公式(2)可得
f (x , y ) =F -1[ρF b (ρ, β)]
-
(3)
从式(3)可知，为了获得真实的密度函数f (x , y ) ，可以先求出反投影函数f b (x , y ) 的傅立叶变换F b (ρ, β) ，在频域中对F b (ρ, β) 加上权重ρ后求其逆傅立叶变换，就能得到所要的密度函数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 滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限，所以频域中的滤波函数
ρ可表示为
⎧ρ,
H R -L (ρ) =⎨
⎩
ρ≤ρ0
(4)
其滤波函数如图7所示。
图7 R-L 滤波函数
连续的R-L 卷积函数为
2
h R -L (R ) =ρ0[2sinc (2ρ0R ) -sin c 2(ρ0R )]
(5)
离散化的R-L 卷积函数为
⎧1
n =0⎪4T 2,
⎪
h R -L (nT ) =⎨0, n 为偶数
⎪1
⎪-222, n 为奇数⎩πn T
图8（a ）（b ）分别为连续形式、离散形式的R-L 卷积函数。
(6)
图8 R-L 卷积函数，（a ）连续形式，（b ）离散形式
R-L 滤波函数的特点是函数形式简单，重建的图像轮廓清晰。存在的问题是：由于在频域中用矩形函数截断了滤波函数，在相应的空域中造成振荡响应，即所谓的Gibb’s现象。另外，如果投影数据中含有噪声，则重建的图像质量也不够满意。
（2）S-L 滤波函数
与R-L 滤波函数不同的是，S-L 滤波函数在频域中不用矩形函数去截断
ρ滤波函数，
而是用一些比较平滑的窗函数来约束滤波函数。例如，可以令窗函数W (ρ) 为
W (ρ) =sin c (
于是可得S-L 滤波函数为
ρρ) rect () 2ρ02ρ0
H S -L (ρ) =ρsin c (
其滤波函数如图9所示。
ρρ) rect () 2ρ02ρ0
(7)
图9 S-L 滤波函数
对应的卷积函数为
14ρ1-4ρ0R sin(2πρ0R )
h S -L (R ) =(0) 2
2π1-(4ρ0R ) 2
其离散化的卷积函数为
(8)
h S -L (nT ) =-
2
π2T 2(4n 2-1)
(9)
图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 1
图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 interpolation
50
100
150
200
250
300
350
图13 180个投影角、不用滤波、线性插值重建得到的图像
Sinogram 1: 180 projections, Ram-Lak (Ramp) filter, Linear interpolation
[***********]50
50
100
150
200
250
300
350
图14 180个投影角、采用Ram-Lak 滤波、线性插值重建得到的图像
7．图15为180个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像。仔细比较图14与图15，可以看到图15比图14略微清晰，原因是三次插值比线性插值效果要好，但三次插值计算量比线性插值要大很多。
Sinogram 1: 180 projections, Ram-Lak Hamming filter, Cubic interpolation
50
100
150
200
250
300
350
图15 180个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像
8．图16为30个投影角、采用Ram-Lak Hamming 滤波、三次插值重建得到的图像。与图15相比，可见投影角度数目，即投影数据较少时，图像有线状条纹，比较模糊。
Sinogram 1: 30 projections, Ram-Lak Hamming filter, Cubic interpolation
50
100
150
200
250
300
350
图16 30个投影角、采用Ram-Lak Hamming滤波、三次插值重建得到的图像
9．分别选择不同的投影数据，选用不同的滤波器、插值函数及投影角数目重建图像，并比较重建图像的差异。
10．在MATLAB 命令窗中输入“type iradon”命令，窗口显示逆雷登变换的MATLAB 程序，
读懂该程序，在关键的代码后用中文注释。
11．画出程序中用到的Ram-Lak Cosine，Ram-Lak Hamming，Ram-Lak Hann滤波器的滤波函数及其卷积函数。
五、实验结果
1．采用高阶插值函数其优点是重建图像质量好，其缺点是
2．如果重建后的图像细节很多，即高频分量多，则该选用
六、结果讨论与思考题
1．为什么投影角度数目少，重建图像不清晰？
2．为什么滤波反投影重建图像要比直接反投影重建图像要清晰得多？