毕业论文-基于matlab实现的ct重建算法仿真比较研究

毕业论文-基于matlab实现的ct重建算法仿真比较研究内容摘要：

个小方块，每个小块的厚度为 d ，块内视为均匀介质。 这样，每个小块有一个均匀的线性衰减系数，如图所示， 设各个方块的线性衰减系数分别为 nuuuu ,..., 321。 X射线通过每一个方块的衰减为 dueII 101  ，射线通过第二个方块的衰减为dueII 212  ，通过第 n 个方块的衰减为 dunn neIII 。 将上述各式代入整理，有： 图 X射线通过非均匀介质时的体素划分 duuuu neII ).. .(0 321  （ ） 即 dpuuuu n /...321  （ ） 式中， )ln( 0IIp 是通过射线强度检测得到的检测值，是一个已知量，称为投影。 这样，上式中，等号右边可视为一个常数，等待求 解的是等式左边的每一个方块的衰减系数。 由于未知量太多，这个方程没有唯一解。 考虑到围绕受检体的一次扫描可以得到一个有关衰减系数的方程，如果将整个感兴趣断层都划分为大小相等的小方块，每个小方块由于一个固定的衰减系数，这样，上述单束射线环绕扫描得到的每个穿透强度值都可以形成一个关于衰减系数的方程。 只要扫描线不重复且足够多，这些衰减系数的确定值肯定可以通过解线性方程组的方法解第二章 CT 图像重建的基本原理 8 出唯一实数解。 把每个方块的衰减系数值用灰度表示，就可以重建出以衰减系数为特征的断层图像。 这种图像重建方法称为方程组法。 值得说明的是，当上图中 0d 时，投影 p 就是 dlluIIp  )()ln( 0 （ ） 式中，衰减系数 )(lu 是随路径 l 连续变化的函数。 上式说明，入射 X 射线强度与出射 X 射线强度之比经对数运算后，表示沿射线路径上衰减系数的线积分，而投影 p 与射线穿越介质的路径长度成正比。 然而，一般二维断层图像至少应划分出 160 160 个像素，如要用方程组法重建这个断层图像，就要 25600 个独立的方程联立求解，这是一个相当费时的任务。 尽管 1967 年世界首台 CT 试验机采用的就是这个图像重建方法，但当断层被划为更小更多的体素时，即使采用今天最先进的计算机技术，解一组这样的方程组也不是件容易的事。 所以 ，这种方法在后来的实际应用中已不再使用，而是采用下一章将要介绍的反投影解析法。 第三章 CT 图像解析重建算法 9 第三章 CT 图像解析重建算法 通过上一章已经知道 CT 图像重建的基本原理，但由于使用最简单直接的方程组法不易求得断层的各个像素的衰减系数。 这一章介绍现在 CT 最常用的图像重建算法 —— 反投影解析法。 主要介绍了由投影重建图像的系统模型，引进Radon 变换和傅立叶变换，它们对应于 X 射线 CT 的二维重建问题；且给出了Radon 变换在图像重建中的具体形式，即截面函数沿着直线的线积分等于其 Radon 变换，而图像重建过程即是将截面函数沿许多 不同角度下直线的线积分所产生的投影数据进行逆 Radon 变换从而得到截面函数。 另外由逆 Radon 变换的推导给出了两种常用的图像解析重建法：滤波反投影法和卷积反投影法。 反投影重建法（ BP） 反投影是一个求和的过程，它把所有角度上的数据都加在一起。 它把投影域中的数据沿着原来的投影路径“涂抹”回去，但不改变数据的值 ]3[。 直接反投影重建是指，断层平面中某像素的线性衰减系数可由平面内所有经过该像素射线的投影和得出，而整幅图像重建 可以看做为所有方向上扫描投影的累加。 通俗地讲，所谓直接反投影，就是在对某体面一个方向的扫描完成后，以得到的投影为灰度值沿着扫描路径经过的体素回抹到体素对应的像素上。 改变方向后的多次扫描形成多次回抹，同一像素上多次回抹的灰度累加即完成图像重建。 图 断层的像素和射线模型 第三章 CT 图像解析重建算法 10 图 反投影重建示例 如图 所示，其中（ a）为原图像像素，（ b）为反投影重建后图像，（ c）为（ b）中像素值除以投影线数后得到的平均值。 有上述例子可以看出原图中像素值为零的像素，经反投影重建后不再为零。 即有伪迹。 直接反投 影图像产生星状伪迹的原因在于，该方法是把取自有限物体空间的投影均匀地回抹（反投影）到了射线所及的无限空间的各个像素上，包括原来像素值为 0的点。 图 反投影重建的等效系统 设 ),( yxu 为处于 0,0  yx 处的一个点源，即二维断面 ),( yx 上的一个冲击函数 ),( yx ，这时的输出 ),( yxf 就是系统的冲击响应 ),( yxh。 已知 221),( yxyxh   （ ） 可见，反投影重建过程后得到的图像不是原来的图像 ),( yx ，系统存 在星状伪迹。 根据信号与系统理论，对于图 有： ),(**),(),( yxhyxuyxf  （ ） 式中，两个星号表示二维卷积 ]4[。 第三章 CT 图像解析重建算法 11 如果想要去除星状伪迹，可在系统的输出端再加上一个滤波器，设其时域函数为 ),( yxq ，传递函数为 ),( 21 Q ，为了使该加了滤波器的系统输出的图像等于系统输入的原图 ),( yxu ，也就是要求 ),(),(**),( yxyxhyxq  （ ） 即 ),(1**),(22 yxyxyxq   （ ） 对该式取二维傅立叶变换，得到： 11),(2221  yxQ  （ ） 即 222121 ),(  Q （ ） 这是一个二维滤波器。 即使加上某种近似，这种滤波器也不容易实现。 如果21 和 的值取无穷，则是一个不可能实现的理想滤波器。 为此提出了滤波反投影图像重建法。 滤波反投影图像重建法（ FBP） 为了消除直接反投影法产生的星状伪迹，提出了滤波反投影重建的方法，这种方法是在直接反 投影重建方法的基础上增加一个滤波器。 增加滤波器的思路有两种，如图 所示。 第三章 CT 图像解析重建算法 12 图 两种不同的去除星状伪迹反投影图像重建思路 所谓滤波反投影重建法是把获得的投影函数先做卷积处理，即人为设计一种一维滤波器，用它对所得的投影函数进行改造，然后把这些改造过的投影函数进行反投影和累加等处理，就可以达到在一定程度上消除星状伪迹的目的。 在介绍滤波反投影重建法之前先补充滤波反投影重建法要用到的数学知识和投影定理 —— 中心切片定理。 Radon 变换与傅立叶变换 为了有效和快速地对图像进行处理，常需要将原定 义在图像空间的图像以某种形式转换到另外一些空间，并利用在这些空间的特有性质方便地进行一定的加工，最后再转换回图像空间以得到所需的效果 ]7[。 通常 正变换：图像空间到其他空间 反变换：其他空间到图像空间 Radon 变换 Radon 变换 1917 年由奥地利数学家 Radon 提出，是 CT 图像重建思想的数学表达。 Radon 从理论上证明了某物理量的二维分布函数可由该函数在其定义域内的所有方向上的线积分完全确定。 Radon 变换的意义在 于，只要知道了一个未知二维分布函数的所有方向上的线积分，那么就能够求得该二维分布函数。 所谓CT 断层的图像重建，就是求取能够反应断层内部结构和组成的某物理量的二维分布。 断层扫描获得的投影实际上可视为被测物体断层内部结构和组成的不同方向上的线积分，所以 Radon 变换的正变换和逆变换正好对应了 CT 技术的射线投影和图像二维分布函数的重建过程。 第三章 CT 图像解析重建算法 13 对 ),( yxf 的 Radon 变换 ),( pRf 定义为沿由 p 和  定义的直线 l 的线积分。 若设 ),( yxf 表示物体的一个二维断层分布，通过 Radon 变换的到 f 沿不同方向直线的线积分相当于得到物体不同方向的投影。 图 定义 Radon 变换的坐标系统 上述线积分可写为：  dlyxfpR f ),(),(  （ ） 如果借助 Delta 函数，上述线积分还可写为：     d x d yyxpyxfpR f )s inc os(),(),(  （ ） 上式通常称为 ),( yxf 的 Radon 变换。 由上述可知， ),( yxf 关于某直线的 Radon 变换就是 ),( yxf 沿该直线上的一维投影。 实际上当将  也视作可变参数时， ),( pRf 就是一个二元函数，此时 Radon变换是空间域 ),( yx 到变换域 ),( p 的映射。 域 ),( p 中的每个点对应于空间域),( yx 中的一条直线，这里 ),( p 实际上就是《图像分析》中直线的 Hough 变换 ]10[。 傅立叶变换的定义   dedeftf ii  ])([21)( （ ） 其中 第三章 CT 图像解析重建算法 14 dtetfF ti  )()( （实自变量的复值函数） （ ） 称为 )(tf 的 Fourier 变换，记为 )]([ tfF。   deF ti)(21  称为 )(F 的 Fourier 逆变换，记为 )]([1 fF。 中心切片定理 中心切片定理的含义是平行投影的一维傅立叶变换等同于原始物体的二维傅立叶变换的一个切片。 将已知投影数据通过一个简单的二维傅立叶反变换可以得到物体截面的一个评估。 图 中心切片定理示意图 通过上图可知：某断层（或它对应的图像） ),( yxf 在视角为  时得到的平行投影（函数）的一维傅立叶变换，等于 ),( yxf 二维傅立叶变换 ),( 21 F 过原点的一个垂直切片，且切片与 轴 1 相交成  角 ]9[。 数学表达如下： 首先定义代表截面的函数的二维傅里叶变换： d x d yeyxfvuF vyuxj  )(2),(),(  （ ） 类似的定义某个角度  下的一条投影 )(tP 的傅里叶变换： 第三章 CT 图像解析重建算法 15 dtetPS tj   2)()(  （ ） 举个最简单的例子：沿着  =0 的直线方向，物体投影的傅里叶变换等于频域里面  = 0 的情形 d x d yeyxfuF uxj  2),()0,( （ ） 因为相位因子不再依赖于 y ，在此可以将积分分成两部分， dxedyyxfuF uxj  2]),([)0,( （ ） 从平行投影的定义看，可以将上式括号中的部分看成是沿着常 量 x 的积分 dyyxfxP  ),()(0 （ ） 将其代入 ()得到 dxexPuF uxj  20 )()0,( （ ） 等式右边表示的是投影 0P 的一维傅里叶变换；因此可以得到垂直方向投影和投影函数的二维变换的关系式： )()0,( 0 uSuF   （ ） 这是傅里叶切片定理最简单的形式。 很明显物体与坐标系统间的位置方向是独立的。 例如，在图 中，如果 ),(st ,坐标系统为 ),( yx ,坐标系统旋转  ，式()定义的投影的傅里叶变换等于物体沿着一旋转了  角的直线的二维傅里叶变换。 从而引出了成像系统中著名的中心（傅里叶）切片定理 ]2[。 即某图像 ),( yxf 在视角  时的平行投影的傅里叶变换给出了图像二维傅里叶变换 ),( vuF 在与 u 轴夹角为  的一个切片，此切片通过原点。 换句话说，的 )(tP傅里叶变换在数值上对应于 ),( vuF 沿着图中所示直线 BB 的值。 通过将 ),( yx 坐标系统旋转得到 ),(st 坐标系统后，傅里叶切片定理可以得到更广泛的应用。 在 ),(st 系统中沿着常量为 t 的直线进行的投影可表示成：  dsstftP ),()( （ ） 第三章 CT 图像解析重建算法 16 其傅里叶变换   dtetPS tj   2)()( （ ） 将以上两式合并，可以发现    dtedsstfS tj   2]),([)( （ ） 借助旋转坐标定义可以得到关于 ),( yx 坐标系统的变换形式    d x dyeyxfS yxj ))s i n ()c o s ((2),()(   （ ） 等式右边表示一个空间域的傅 里叶变换， )s in (),c os (   vu 则 ))s in ()c o s ((),()(  ，FFS  （ ） 这个。