1、CT、MRI图像重建算法上海理工大学 聂升东泰山医学院 邱建峰常用的重建算法 二维傅利叶变换法（中心切片理论）普通反投影 反投影法 平行线束 滤波反投影 去伪影重建 （雷登变换）扇形线束 卷积反投影 投影重建 迭代法（数值迭代）拟和逼近法（算术拟和）两个重要工具迪拉克函数-抽样卷积-滤波 傅利叶变换是将任意一周期信号或非周期信号变换成用其自身频率特性表述的一种形式，使信号的变化与频率变化之间建立内在的联系，从分析频率特性的角度来揭示信号本身的变化观律。如图所示的矩形波信号，经过傅里叶变换成频率变化的形式。二维傅里叶变换法将各个投影进行，再把各角度上的变换结果汇集起来，在变换成求得的傅里叶变换的
2、频域曲面，再改为空间直角坐标。按公式进行后即可得到重建图像。二维傅里叶变换法是最理想的图像重建方法之一。但该方法需要进行正、反两次傅里叶变换，计算量比较大，在实际应用中不易实现。二维傅里叶变换法理论依据：中心切片理论优点：算法精确缺点：复杂耗时，运算量庞大定义密度函数f（x，y）在某一方向上的投影函数g （R）的一维傅里叶变换函数 G （p），是密度函数f（x，y）的二维傅里叶变换F（p，）在p，平面上沿同一方向过原点直线上的值。中心切片理论中心切片理论的证明特殊情况 =900一般情况中心切片的应用指出了投影重建图像的可能性二维傅里叶变换法如果我们在不同角度下取得足够多的投影函数数据，并进行傅
3、里叶变换，那么变换后的数据将充满整个（u，v）平面。得到频域函数F（u，v）的全部值，将其进行傅里叶反变换，得到原始的密度函数f（x，y）即所要的图像。填充（u，v）平面需要按照角度填充，即筛选角度，完成这一步骤需要做极坐标转换。空间域雷登空间极坐标频率域二维傅里叶变换法算法：描述/过程图验证：sinc（2R）缺点：插值/运算量大反投影法(back projection)又称总和法，是利用投影数值近似地复制出值的二维分布。基本原理是将所测得的投影值按其原路径平均的分配到每一点上，各个方向上投影值反投影后，在影像处进行叠加，从而推体出原图像。反投影重建算法的一般步骤：原像取投影反投影重建重建后图
4、像反投影重建算法设被测人体断面上器官或组织的吸收系数分布为，X线束扫描时在某一角度方向的投影表示为：则在角度的反投影可表示为：式中的是由沿反方向进行反投影所产生的吸收系数，-函数起筛选角度作用。将上式全部角度上的反投影值相加，即对应从0变化到所有反投影值加在一起，可得到图像重建的吸收系数分布为：dxdyRyxyxfRPsincos,dRRyxRPyxbsincos,00sincos,dRRyxRpddyxbyxfb投影重建算法的基本内容：“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和（的平均值）”。式中xk表示象素的值，pk,i为经过象素的第i条射线投影。可以这样理解：其
5、中f(x,y,z)是身体组织密度。,11pnkk ipixpn(,)lpf x y z dl反投影重建算法的物理概念算法举例113711152pxxxx222610146pxxxx3356787pxxxx4491011121pxxxx551611165pxxxx664710133pxxxx算法举例根据反投影算法x1=p5=5x6=p2+p3+p5=18平均化处理，除以投影线数目xi=xi/6000005200100000056237181271108136250.8310.3300.51.16321.160.061.661.330.160.510.330.83反投影重建后原像素值再除以投影线数
6、，平均化断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和的平均值123456星状伪迹反投影重建后，原来为0的点不再为0，形成伪迹00000520010000000.8310.330.51.16321.160.061.661.330.160.510.330.83原像素值再除以投影线数，平均化星状伪迹我们考虑孤立点源反投影重建，中心点A经n条投影线投影后，投影值均为1：1/n1/n1/n1/n11/n1/n1/n1/n121(.)1Anfpppn000010000星状伪迹产生星状伪迹的原因在于：内插内插方法 1.紧邻内插:在 时，若则以 中的内容 直接代替 。否则，用 2.线形内
7、插：如果 ，则 cos()mprm()(1)()*cos()/mmmpndpndpndrmnddmmcos()cos()(1)rmndrmnd (,)mnd cos()mprm()mpnd(1)mpnd cos()mprmcos()(1)ndrmnd 在反投影法中，反投影重建图像的吸收系数与实际图像之间有关系：ryxfyxfb1,校正模糊失真的步骤可先将fb(x，y)作二维傅里叶变换，然后将变换结果用加权，则得真正图像的二维傅里叶变换式，有：在此基础上再行二维傅里叶反变换，按公式获原图像吸收系数分布面函数f（x，y）。显然仍需进行二维傅里叶正、反变换，因此并未减少计算量。,bFF直接直接反投反
8、投2DIFT2DFTX频域频域空域空域普通反投影法仍需二维傅里叶变换先反投再修整滤波反投影算法去除星状伪迹的两种方法：取投影取投影一维滤波器反投影重建二维滤波器反投影重建输入图像（原图像）输入图像（原图像）输出图像(同原图像)输出图像(同原图像)星状尾迹的去除（b）（a）第一种方法中的二维滤波器较难实现。第二种方法交换了滤波与反投影的顺序，此时投影数据是一维的，易于实现。此方法的理论基础仍是中心切片定理先修整再反投滤波反投影算法滤波反投影法也是解析法中的一种。这种方法为了消除模糊因子1/r的影响，并将，既可校正失真，又可简化计算量，提高图像重建速度。采用卷积计算的滤波反投影法当前应用最为广泛，
9、故也称卷积反投影法(convolution back projection，CBP)。反投反投1DFTX 频域频域空域空域滤波反投影法1D IFT卷积反投影法滤波反投影仍需1D傅利叶变换能否更简单？解决方法：卷积、滤波卷积反投影法在空间域进行滤波先滤波再反投，但不需要傅利叶变换 11111FRgRgFF反投反投X C（R）空域空域卷积反投影法成像方法是在后投影之前，对所有的投影数据进行卷积，使结果的图像无所谓的“星月样”晕伪影。成像的过程可分成三步：是获取全部的投影数据并作预处理。是将所得数据与卷积因子相乘，其间须通过大量的数学运算，同时采用的算法还须考虑图像的分辨率和噪声等。，其后投影数据值
10、中正负值相互抵消，根据系统显示的不同矩阵大小，各投影滤过的原始数据被投影成像并通过监视器显示。要考虑的问题不可积，其逆变换不存在不可积，其逆变换不存在需要近似函数代替其逆变换需要近似函数代替其逆变换卷积函数卷积函数待建图像为a(x,y),它的2D傅氏变换为:12(,)(,)AA 12(,)(,)()rAAF p x ()(,)PP 待建图像其中1212(,)(,)(,)a ra x yFA 2cos()0|(,)jrdPed d 2cos()|(,)jrPed cos()cos()()(,)|(,)|cos(),rxrrrrxrrh xp xg xg r 上式的物理意义是：|()rF h x(
11、,)rg xcos()rxr(,)rp x可得0(,)cos(),a rg rd(,)r(,)r称“滤波反投影方程”，它体现了滤波（卷积）反投影算法的三个步骤：1 把在固定视角 下测得的投影 经过滤波，得到滤波后投影 ；2 对每一个 ，把 反投射于满足 的射线上的所有各点 ；3 将步骤2中的反投影值对所有 进行累加（积分），得重建后的图像。ii(,)rip x(,)rig x(,)rig x(,)rcos()rixr 00(,)cos(),a rg rd 滤波函数的选择为了醒目起见，可将式(4.39)中的m 作为常数，于是有：及同理，把 看成常数，由 及 得：)()*()*()(rrp ndh
12、 ndxp x()*()rp ndh x(,)(,)*()rrrg xp xh x(,)(,)rrg xp x)()*()(rrrp xh xp x常用滤波函数举例1.RL滤波函数(i)系统函数式中，B=1/2d，且：()R LH()()2R LHWrectB1,0,2BrectBB(ii)相应的冲激响应 ()R Lrhx2()rBjaxR LrBhxed2222sin(2)sin()rrBcx BBcx B(iii)采样序列采样间隔为d，对应的最高不失真空间频率为1/(2d)的离散形式如下：偶数 奇数 如将 的离散表示进行线性内插，则得另一连续的空域函数 ，是 的一次近似。()R Lhnd(
13、)R Lrhx22221,04()0,1,R Lndhndnnnd()R Lhnd()R Lhnd()R Lrhx2.SL滤波函数选取sinc函数作为窗函数。于是有：相应的冲激响应 为：2()sin(/2)sin/22S LBHcBrectBB()S Lrhx22()sin2rBjxS LrBBhxedB221 4sin(4)1 42()21(4)rrrBxBxBBx如令 y=4BXr，则得：221sin()1 42()()21S LyyBhyy采样序列对 进行均匀采样，得：将这一离散形式进行线性内插，则得连续的空域函数 ：()S Lhnd()S Lhnd2222(),0,1,2,(41)S
14、Lhndndn ()S Lrhx2()()()*()S LrS Lrrrnhxhxxndx对经线性内插后的函数 进行傅立叶变换得：用SL滤波函数重建的图像其振荡响应减小，对含噪声数据的重建质量也较RL滤波函数情况为好，但在低频段不及RL滤波函数的重建质量高。()S Lrhx21sin(/2)()sin2(/2)S Lrddhxdd F重排算法 把一个视图中采得的扇形数据重新组合成平行的射线投影数据，然后用上一节的算法重建。直接重建算法 不必数据重排，只需适当加权即可运用与上一节类似的算法重建等角射线型等距射线型12D D0S：X射线源：检测器阵所在弧线：中心射线00S D 扇形位置由该中心射线
15、与y轴交角 确定，同一扇形中的射线 由规定，射线的绝对位置由（，）唯一确定。0S E重建算法的推导扇束算法平行束算法结论扇束情况下的重建算法较为复杂，但实质没有改变。可采用平行束情况下的算法实现，只需加以适当地修正即可。反投影重建实际实验CTsimMatlab-radonCTsim仿真CT扫描CTsim软件是由Kevin M.Rosenberg等学者编写开发的CT重建仿真程序，其源代码开放。用户可选择herman head、sheep-logan、unit pulse三种体模，进行滤波反投影重建和福利叶变换重建。1、选择要进行扫描的体模。打开CTsim软件，按file-create phant
16、om选择体模sheep-logan。2、对体模进行平行线束（平移-旋转式）扫描投影。按process-projections，在projection parameter中选择geometry-parallel,trace level中选择projections，其它扫描参数使用默认参数。软件进行体模扫描演示，扫描完成后，将各角度投影值按顺序排列在二维平面上。3、选择reconstruct-filtered backprojection，在filtered backprojection parameter中选择filter method为FFT。软件进行滤波反投影重建，选用的滤波方法为快速福里叶
17、变换，获得重建的sheep-logan体模图像。4、比较重建图像与原体模的差异。5、在trace level中选择plot，重复1-3步，观察扫描过程中，每个投影值的变化。6、在projection parameter中选择geometry-parallel（扇形束扫描），重复1-3步，观察扫描过程和重建影像。其扫描过程、投影图像、重建图像。比较投影图像、重建图像与平行线束扫描获得的图像的区别。7、选用unit pulse体模，重复上述过程，观察重建图像中的伪影。8、更换其它体模，进行上述方法，进行重建，观察重建方法与重建图像。MATLAB仿真迭代重建迭代重建同样使用反投影，但不进行滤波。方法
18、：反投影重建图像，对新图像进行投影，将新投影与真实投影比较，进行修正。通常用两投影差反投来改善图像。并进行迭代求最佳。迭代法代数重建算法（ART）最有效、对噪声敏感迭代最小平方技术（ILST）常用ART迭代一、两次，然后转入ISLT同步迭代重建算法 迭代法准确运算量大 通常迭代十次左右，而滤波反投影法只花迭代一次的时间核素成像常用算法MRI的重建算法在MRI设备中，傅立叶成像（fourier imaging）占主流修改的劳特伯投影重建修改：脉冲NMR 使用两个梯度线圈 投影数据进行内插MRI的重建算法90度RF脉冲后，外加梯度场，梯度法线方向投影就是FID信号。用s（t，MRI的重建算法MRI
19、的重建算法对傅立叶变换来讲：CT中，R与K是一对共轭分量；MRI中，w与t是一对共轭分量。因此，用连续波NMR得到的频域投影P（MRI的重建算法脉冲磁共振中：时域投影s（t，MRI的重建算法MRI的重建算法zwx不同方向梯度存在时观察到的谱P（W，）P（Wk，）经一维福利叶变换后给出截面tztxxwzFTFT-21974年NMR谱学家Ernst小组用二维傅立叶成像法充作了劳伯特试验。采集数据方法改变。MRI的重建算法T1T2一次采集64个数据重复64次，但使用64个T1时间形成正交二维时域，但真实时间是一维的MRI的重建算法t2t1Kx=gxt1Kz=gzt2投影重建低频分量重建精度高，高频低二维傅立叶变换法高低频精度相同，图像质量好福利叶成像中，滤波函数即峰形函数自动包含在多维傅立叶变换中，因此算法体现经济性傅立叶成像在效果上是全三维数据采集，信息量增益明显。通过选择性面激发或选择性线激发简化到二维或一维。MRI的重建算法
