ART算法图像重构ppt课件

1、 图像重建是图像处理中的一个重要分支，广泛地应用于物体内部结构图像的检测和观察中，它是一种无损检测技术。 关于图像处理的一些基本内容，如对图像的几何处理，图像的增强，还有复原等，均是从图像到图像，即输入的原始数据是图像，处理后输出的仍是图像。而图像重建是从数据到图像。 图像重建的三种常用检测模型：发射模型、反射模型 计算机层析成像（Computed Tomography, CT）是通过对物体进行不同角度的射线投影测量来获取物体横截面信息的成像技术。 CT的核心技术是由投影数据来重建图像的理论，其实质是由扫描所得到的的投影数据来求出成像平面上每个点的衰减系数值。当强度为 的x-ray通过吸收率为

2、(x,y)的均匀吸收物体，由于均匀吸收，则I必是指数下降，则有0I0exp(y)dybaII0ln(y)dybaII0(x)(x)ln(x,y)dySIpI这里s表示射线经过的体内距离长度0 (x) (x, )(x,y)ln()sIpdyI , x y (x, )p经坐标系旋转变换后可得：而对于任意角度扫描，需要用旋转坐标来描述问题，建立置于扫描系统之上的旋转坐标系 ，即让射线束与旋转坐标系 的 轴平行：所以角每旋转1度就可以取一组投影数据，可得到180组不同的投影。CT就是在收集各角度的投影数据后，利用重建算法处理得到物体的图像。 , x y y是离散值，是测出值！ Radon变换是计算图像

3、在某一指定角度射线方向上的投影的变换方法。二维函数f(x,y)的投影是其在确定方向上的线积分，如下图所示，二维函数f(x,y)在水平方向的线积分就是f(x,y)在y轴上的投影，二维函数f(x,y)在垂直方向的线性积分就是f(x,y)在x轴上的投影。由此，可以沿任意角度 计算函数的投影，计算图像f(x,y)在任意角度的Radon变换。 密度函数在某一方向上的投影函数的一维傅立叶变换函数密度函数在某一方向上的投影函数的一维傅立叶变换函数是原密度函数的二维傅立叶变换函数在平面上沿同一方向是原密度函数的二维傅立叶变换函数在平面上沿同一方向且过原点的直线上值。且过原点的直线上值。1、在不同的角度下取得足

4、够多的投影数据(Radon变换)2、将这些投影数据做一维的Fourier变换，那么变换后的这些数据将充满整个(u,v)平面。（许多过原点成不同夹角的直线）3、也就是说，F(u,v)的全部值都为已知，那么我们将其做一次二维的Fourier逆变换就可以得到原始的衰减系数函数f(x,y)( , )( , )exp 2 ()f x yF u vjuxvy dudv2( cossin )0200( , )( , )( , )( cossin)dR(R) ( cossin)dRjxyjRf x yFed ddFedxyRdgxyR ( , )( , )exp 2 ()f x yF u vjuxvy dud

5、vcosusinv作坐标变换，令：可得出：2(R)( , )jRgFed 表示对投影函数的Fourier变换进行滤波变换，其中 是滤波函数。所以要实现对投影数据实现图像重建，可以采取两步：首先将投影数据和响应脉冲滤波器进行卷积，然后由式对不同旋转角求和，就能实现图像重建。这就是卷积法进行图像重建的基本思路和方法。0( , )(R)(R)( cossin)dRf x ydghxyR所以：式中h(R)为滤波函数纠的空域形式 基本原理是将所测得的投影值按其原路径平均的分配到每一点上，各个方向上投影值反投影后，在影像处进行叠加，从而推体出原图像。 而滤波却是要投影函数的一维Fourier加上权重因子。

6、113711152pxxxx222610146pxxxx3356787pxxxx4491011121pxxxx551611165pxxxx664710133pxxxx根据反投影算法x1=p5 = 5x6=p2+p3+p5=18平均化处理，除以投影线数目xi=xi/6反投影重建后原像素值再除以投影线数，平均化断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和的平均值123456 反投影重建后，原来为0的点不再为0，形成伪迹原像素值再除以投影线数，平均化我们考虑孤立点源反投影重建，中心点A经n条投影线投影后，投影值均为1：p1=p2=.=pn=1因此重建后而其他点均为1/n这类伪迹

7、称为星状伪迹121(.)1Anfpppn 产生星状伪迹的原因在于：反投影重建的本质是把取自有限物体空间的射线投影均匀地回抹(反投影)到射线所及的无限空间的各点之上，包括原先像素值为零的点（其实就是投影数据少产生的！）(a)孤立点源(b)反投影重建图像及星状伪迹 滤波反投影法采用先修正、后反投影的做法，其基本方法是：在某一投影角下取得了投影函数（一维函数）后，对此一维投影函数作滤波处理，得到一个经过修正的投影函数；然后再将此修正后的投影函数作反投影运算，得到所需的密度函数。 滤波反投影法重建图像有以下几个步骤： （1）对某一角度下的投影函数作一维傅立叶变换； （2）对（1）的变换结果乘上一维权重

8、因子； （3）对（2）的加权结果作一维逆傅立叶变换； （4）用（3）中得出的修正过的投影函数做直接反投影； （5）改变投影角度，重复（1）（4）的过程，直到完成全部180度的反投影。滤波函数的选取是滤波反投影法的关键问题 （1）R-L滤波函数 由于在频域中用矩形函数截断了滤波函数，在相应的空域中造成振荡响应 ，重建的图像质量也不够满意 0( )(/ 2)Hrect对应的频域形式为：理想的滤波函数 它是在高频的权重很大，低频的权重很小，所以高频噪声就会很大，所以我们才要对其进行修正（2）S-L滤波函数 与R-L滤波函数不同的是，S-L滤波函数它的关键是把频域的陡峭截止改成缓慢截止。 用S-L滤波

9、函数重建的图像中振荡相应较小，对含噪声的数据重建出来的图像质量也较R-L滤波函数重建的图像质量要好。但是，S-L滤波函数重建的图像在高频响应方面不如R-L滤波函数好，这是因为S-L滤波函数在高频段偏离了理想的滤波函数 对应的频域形式为：00( )sin ()()22S LHcrect %P=imread(lena.jpg); P = phantom(256); %P= rgb2gray(O); R = radon(P,0:179); I0 = iradon(R,0:179,linear,Ram-Lak);I1 = iradon(R,0:179,linear,Shepp-Logan);I2=ir

10、adon(R,0:179,linear,cosine);I3 = iradon(R,0:179,linear,none);subplot(2,3,1), imshow(P), title(Original)subplot(2,3,2), imshow(I0,), title(FBP R-L)subplot(2,3,3), imshow(I1,), title(FBP S-L)subplot(2,3,4), imshow(I2,), title(FBP cosine)subplot(2,3,5), imshow(I3,), title(Unfiltered BP)图像的细节对应的是高频部分，轮廓

11、对应的是图像的低频部分，所以因为没有滤波，细节部分恢复的不好，呈现很“模糊”的情况 一个典型实例： 在matlab图像处理工具箱中，有一个phantom函数，可以用来创建头部的剖视图，首先创建一个头部的256256剖视图，然后分别计算3组不同的Radon变换，第一组采用30个投影，第二组采用90个投影，第三组采用180个投影，用以比较采用不同组数的投影参数重建的图像与原始图像的差别。 由测试结果可以看出：第一组采用30个投影，效果较差；第二组采用90个投影，效果较好；第三组采用180个投影，效果很好，与原始的图像非常接近。 这说明可以通过增加投影的数目，来提高重建图像的质量。 滤波反投影算法要

12、求投影数据必须完全，分布必须均匀。具体地说就是，平移采集投影数据时，应覆盖全部物体区域；相邻射线间均为d。 然而实际应用中，有时无法测到大量的投影数据。例如做CT时，为了避免心脏器官受辐射过久，为了减少剂量，投影数据采集不足； 迭代重建算法就能解决上述滤波反投影算法无能为力的场合啦！ART（代数重建法）是一个迭代的过程，它是一开始就在离散域中进行的,首先把图像离散化，即将欲重建的未知图像离散成一个J=n*n重建图像网格。根据成像的物理过程和相应的数学模型建立待重建图像和投影数据之间的代数方程组，那么图像重建问题就可以转化为解线性方程组问题。x1x2x3x4x5X6x7x8x9x1x2x3x4x

13、5X6x7x8x91号 2号 3号 4号=nn划分为 J个有限像素21ijS ABCr称为加权因子重建模型一ijiijjii=ijjijjjjr xPPr x于是，像素 对射线i的射线投影的贡献为P因为 号射线还经过其他像素，所以 号射线总的射线投影为ijijr1,jr =0,ii第二种表示射线和p 及定义 的方法是：认为像素值集中在像素的中心。即号射线（宽度 ）经过 号像素中心其他ijijr1,jr =0,i第三种定义 的方法是：令射线的宽度为零，间距依旧是 。那么号射线（宽度0）经过 号像素内任一点其他x1x2x3x4x5X6x7x8x911111221993221122229913589

14、33113223991458944114224997551152259912661162269923457.pr xr xr xxpr xr xr xxxxxxpr xr xr xxxxxxpr xr xr xxpr xr xr xxxpr xr xr xxxxxp711722799567888118228999.r xr xr xxxxxpr xr xr xxx1x2x3x4x5X6x7x8x91号5号4号8号以模型三为例，不失一般性如图为一3*3像素的图像的x1,x2,x9为相应的像素值，则各射线和为：9i112812811121314151617181921222324252627282

15、9818283848586878889=1,2,80010000001100110011001=ijjjTTPr xiPRxPPPPxxxxrrrrrrrrrrrrrrrrrrRrrrrrrrrrR上式可写成更简洁的形式：或者用矩阵表示为：式中本例中10011000000100110000000011110000000011110000000001pPRxxJ是 维矢量，称为图像矢量，是I维矢量，称为测量矢量，R是I J维矩阵，称为投影矩阵。我们的任务就是根据测量所得的P和已知的投影矩阵R（规定了像素布置与射线的几何结构，R就可以求出来）即可求出图像矢量x。-1-1RRR p求矩阵 的逆矩阵，

16、也即是x=2449JI=2566.25 10=J n nI此方法虽很直观，但是不实用，因为像素数目 及射线数目 都极大，直接求逆矩阵的运算量很大。例如，假设图像划分为像素，取I=J，则R的维数就为I J （6.25 10 ） 3.9 10而且，在很多情况下，投影数据的个数 小于未知像素J，于是有无穷多个解。1jj1,2,IjijiiTxr pjJxR p第二种办法就是用反投影算法重建（累加法），即把经过 号像素的所有射线和累加，得到 号像素值。写成矩阵的形式为此方法会产生严重的伪迹。不是我们所期望的。但此方法的思想有助于我们理解迭代重建算法。 基本原理是将所测得的投影值按其原路径平均的分配到每

17、一点上，各个方向上投影值反投影后，在影像处进行叠加，从而推体出原图像。 113711152pxxxx222610146pxxxx3356787pxxxx4491011121pxxxx551611165pxxxx664710133pxxxx根据反投影算法x1=p5 = 5x6=p2+p3+p5=18平均化处理，除以投影线数目xi=xi/6反投影重建后原像素值再除以投影线数，平均化断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和的平均值123456 反投影重建后，原来为0的点不再为0，形成伪迹原像素值再除以投影线数，平均化我们考虑孤立点源反投影重建，中心点A经n条投影线投影后，投

18、影值均为1：p1=p2=.=pn=1因此重建后而其他点均为1/n这类伪迹称为星状伪迹121(.)1Anfpppn 产生星状伪迹的原因在于：反投影重建的本质是把取自有限物体空间的射线投影均匀地回抹(反投影)到射线所及的无限空间的各点之上，包括原先像素值为零的点（其实就是投影数据少产生的！）(a)孤立点源(b)反投影重建图像及星状伪迹(0)(1)(1)(2)(k)(k 1)( )(0)( )( )(1)(xxxxxxx1,2,;kTiikTkikikkkARTRxpxpiIxxr xpxx基本思想：先给定一个初始图像求一次近似图像，再据求二次近似图像，如此继续，直至满足预定条件而后止。在根据求时需

19、加一校正值。满足线性不等式而对于图像重建这一具体的情况，我们可以将上式写成r而上式的解，也即ART的迭代公式如下：任意；( )( )2(mod )+1Tkkikikikikkpr xrrikI其他其中我们来分析一下其几何意义和物理意义(从一个例子入手)121212112212=2=0.40.2r =r =p =2.4p =3.00.10.50.40.12.40.20.53.0TTJxxr xpr xpxx已知像素数，射线投影数I 2，求图像矢量x。由可得121112212( )(0)x:0.40.12.4:0.20.53.0819kLxxLxxx 以x 为横坐标，为纵坐标，对应于直线 与直线L

20、为简单起见，令。取(0)(0)0112111(0)(1)(0)11121(1)180.40.14.12.490.40.40.10.170.1+80.442.4-4.1=+=90.180.17TTiTTxxprr rpr xxxrrxH 那么rr由图可见，刚好落在直线上x( )(1)22222(1)(2)(1)222222(2)240.20.54.880.20.20.50.290.5+40.280/ 293.0-4.8=+=80.5142/ 2903.29.0TkTikTTxxrr rpr xxxrrxHp 那么rr由图可见，刚好落在直线上( )(2)11(2)(3)11(2)(3)(4)(5)(2)80/ 290.40.11.592.4142/ 2980/ 29=142/ 29=TkTikTxxpxpxxxxxx那么rr因为r，故有仿照上述方法，可以求得也就是说x收敛于。所以就是所求的解啦！(0)(1)(0)11121+Tpr xxxrr迭代公式：(0)11121(0)11111TTpr xrrpr xrrr校正项改写成(0)1H代表x 到的距离1111rrHr是 的