（应用数学专业论文）基于分块迭代的快速代数重建算法研究.pdf

基于分块迭代的快速代数重建算法研究 摘要 计算机层析成像( 简称c t ) 技术检测精度高、重建的图像具有无影像 重叠、空间分辨率和密度分辨率高、可以直接数字化处理与分析等特点，成 为近十凡年发展起来的一种新的无损检测技术，广泛应用于航空、航天、机 械、船舶、公安、海关等诸多领域。 c t 技术分为硬件和软件两部分。重建算法是软件部分的核心，其包括 变换方法和级数重建法等。变换方法以卷积反投影( f b p ) 方法最为常用， 该算法重建速度快，成像质量较好，但是其要求完全的、等间距的采样数据， 通常重建图像伪影较重。代数重建方法( a r t ) 是级数重建法的典型形式， 其适用于不同方式的采样数据，对不完全数据也可重建图像，但是，计算量 大、重建速度慢，影响了该算法的应用范围。 本文主要针对经典a r t 方法收敛速度慢的缺点，分析了影响它收敛速 度的关键因素，提出了一种基于分块迭代的快速a r t 方法。其基本思想是 对要重建的图像逐级分块，通过迭代使图像收敛速度逐渐加快，逼近于重建 的图像。在此方法中运用了求解投影系数矩阵的快速计算方法，并对投影数 据按相关性最小的原则进行排序。利用分块迭代的快速a r t 方法，对x 射 线工业c t 实采投影数据进行图像重建，并与f b p 方法、经典a r t 方法重 建的图像进行了比较，测试结果表明：该方法重建图像精度高，伪影轻，并 有较高的密度分辨率和空间分辨率，较经典a r t 方法重建速度快。 在图像重建中，投影数据的采集方式、实际物体的模型参数、重建图 像的大小等，都刑重建图像的质量有一定影响。针对上述因素，本文在带噪 声投影数据处理、不完全数据重建、投影数据插值、低分辨率数据重建高分 辨率图像等方面进行了一系列改进试验。数值试验结果表明，相应重建图像 的空问分辨率和密度分辨率得到不同程度的改善。 本文研究结果应用于国产工业c t 机专用成像软件，对研制国产高分辨 率工业c t 机及其成像软件有重要意义。 西安理工大学硕士学位论文第i 页 基于分块选代的快速代数重建算法研究 西安理工大学硕士学位论文 第i i 页 基于分块迭代的快速代数重建算法研究 a b s t r a c t t h ei m a g eo fc o m p u t e dt o m o g r a p h y ( c t ) i so fg o o dq u a l i t ya n dh i 曲 r e s o l u t i o ni tc a r lb ep r o c e s s e da n da n a l y z e dd i g i t a l l y t h e r e f o r e ，c th a sb e c o m ea n e wn o n d e s t r u c t i v ed e t e c t i o nt e c h n i q u e ，w h i c hi s a p p l i e dw i d e l yi na v i a t i o n ， s h i p p i n g , a u t o m o b i l e ，p u b l i cs e c u r i t ya n d o t h e ra r e a s t h ef i l t e rb a c k - p r o j e c t i o n ( f b p ) t e c h n i q u ea n dt h ea l g e b r a i cr e c o n s t r u c t i o n t e c h n i q u e ( a r t ) a r et w oo ft h em a j o ri m a g er e c o n s t r u c t i o nm e t h o d si nc t t h e r e c o n s t r u c t i o ns p e e do ff b pi sf a s t t h ei m a g er e c o n s t r u c t e db yf b pi so fh i g h s p a t i a lr e s o l u t i o n ，b u tw i t i lh e a v ya r t i f a c t s ，w h i l et h ei m a g er e c o n s t r u c t e db ya r t i so f h i 曲- d e n s i t yr e s o l u t i o n a n d i so f s l i g h ta r t i f a c t s i n o r d e r t or e c o n s t r u c t i m a g e ， g e n e r a l l yf b pn e e d sc o m p l e t ep r o j e e t i o nd a t a ，b u ta r t i sa v a i l a b l et oe i t h e r c o m p l e t eo ri n c o m p l e t ep r o j e c f i o nd a t a t h em a j o rd i s a d v a n t a g eo f a r t i si t sl o w c o n v e r g e n ts p e e d i nt h i sp a p e r ，t h ea u t h o rp r o p o s e sam o d i f i e da r t ( m a r t ) a l g o r i t h mb a s e d o nb l o c k i t e r a t i o n ，a f t e ra n a l y z i n gt h ef a c t o r st h a te f f e c tt h ec o n v e r g e n ts p e e di n t h ec o n v e n t i o n a la r t t h ea l g o r i t h mc o n v e r g e sf a s t s oi ti sa v a i l a b l et o r e c o n s t r u c tm o d e r a t ea n dl a r g ei m a g eb yu s i n ge i t h e rp a r a l l e lo rd i v e r g e n t p r o j e c t i o nd a t a ，a sw e l la st h ed a t ac o l l e c t e da ts e r i e sa r b i t r a r ys a m p l e da n g l e sm a d s a m p l e dr a d i i a n di td o e sn o tr e q u i r et h es t o r a g eo fp r o j e c t i o nm a t r i x ，c o m p a r e d w i t ht h ec o n v e n t i o n a la r tt h en u m e r i c a lr e s u l t ss h o wt h a tm a r tr e c o a s t r u c t s i m a g ew i t hf a s ts p e e d ，h i g hr e s o l u t i o na n ds l i g h ta r t i f a c t s k e y w o r d s ：c o m p u t e dt c m o g r a p h y f b pa r t i m a g e r e c o n s t r u c t i o nb l o c k - i t e r a t i o n 西安理工大学硕士学位论文第l i i 页 基于分块选代的快速代数重建算法研究 1 c t 概述 第一章绪论 计算机层析成像( c o m p u t e dt o m o g r a p h y ，简称c t ) 技术是利用具有 某种能量的射线源( x 射线，y 射线) 对物体进行断层扫描。根据在物体外部 所获得的某种物理量( 指物质对射线的衰减系数) 的一维投影数据，运用特 定的重建算法。进而重建物体特定断层的无影像重叠的二维图像。 c t 机有三个组成部分。一是数据采集系统；二是图像重建，主要是图 像重建算法；三是图像后处理和显示系统。 c t 断层图像具有无影像重叠、空间和密度分辨率高、可直接进行数字 化处理等特点，是近十几年发展的新一代非接触性无损检测技术，成为关键 部件检测、机械仿型设计、安全检查的方面强有力的手段，广泛应用于航空、 航天、机械、汽车、船舶、公安等领域。 2 国内外研究情况 自7 0 年代第一台c t $ j l 问世以来，计算机层析成像技术一直是国内外研 究的热点之一，其研究可分为硬件和软件两方面。硬件方面发展出具有低、 中、高或超高能的y 射线、x 射线、电子束、直线加速器的c t 机。软件方面 结合检测系统从数学模型、数据采集、成像方法、软件实现等角度，提高成 像分辨率和检测速度。特别值得注意的是些通过软件方法替代部分硬件功 能的技术，由此达到降低成本，提高效率等目的。 c t 技术的理论由德国数学家r a d o n 建立( 1 9 1 7 ) ，6 0 年代由美国的 a m c o r m a r k 等人进一步发展。在7 0 年代，世界上第一台医用c t 机由英国人 c n h o u n s f i e i d 研制成功。a m c o r m a r k 和c n h o u n s f i e l d 两人由此获1 9 7 9 年 诺贝尔医学奖。在医学方面，西方发达国家已先后研制出具有高分辨率的螺 旋c t 、可超高速成像的电子束c t 等设备。工业方面，7 0 年代末起美国利用 西安理工大学硕士学位论文第1 页 基于分块迭代的快速代数重建算法研究 研制的透射式工业ct 设备对军工产品的关键部件做无损检测，美国科学测 量系统公司的工业ct 机还应用于c a d c a m 方面，进行加工元件的仿型制 造，开展逆向工程学的研究。 8 0 年代中期，国内一些大学、研究所开始从事c t 理论与应用试验的研 究，其中有清华大学、中科院高能所、重庆大学、中国计量科学院、上海交 通大学、北京信息工程学院、东北大学、浙江大学等单位。9 0 年代初我国自 行研制成功第一台透射式y 射线工业ct 样机，清华大学等单位研制成功x 射线微焦点i c t 机，打破了西方国家在对我国的封锁局面。但是，由于国内 工业c t 的关键器件在研制过程中，目前国内用于军事工业和民用工业关键 材料、主要部件的无损检测的工业c t 设备部分依赖于进口。在工业c t 关键 技术f 包括成像软件) 上西方国家对中国一直采取技术封锁和保密，以保持其 在工业c t 上的垄断地位。例如出口到中国的工业c t 机的数据采集和成像系 统被固化加密，无法二次开发和改进，维修也依靠国外厂家。此外，还限制 对华出口某些类型的高分辨率3 2 、j k c t 机。因此发展具有自主知识产权的国 产高分辨率工业c t 机( 包括成像软件) 是非常迫切和必要的。 3 本文主要工作 本文主要针对经典a r t ( 简n ：a r t ) 算法收敛速度慢的缺点，并对影响 a r t 算法收敛速度的主要因素，如投影系数矩阵( 简称系数矩阵) 计算、迭 代格式、收敛准则等进行分析，提出了一种基于分块迭代的快速a r t 算法。 其基本思想是采用对图像逐级分块，通过迭代使图像逐步细化，最终逼近于 重建的图像。该算法中运用了求解系数矩阵非零元的快速计算方法，按相关 性最小的原则对投影数据进行拟正交排序。用实采数据对该算法进行了测 试，结果表明：该算法与经典a r t 方法比较，速度有了较大提高，重建图像 伪影轻，有较高的密度分辨率和空间分辨率。 重建图像的过程中，除考虑重建算法外，投影数据的采集方式、数据量 的大小、实际物体的模型参数、重建图像和投影数据的关系等上述因素对图 像质量也有一定程度的影响。因此，本文结合上述因素进行带有噪声投影数 西安理工大学硕士学位论文第2 页 基于分块迭代的快速代数重建算法研究 据处理、不完全数据重建、投影数据插值、低分辨率数据重建高分辨率图像 等方面进行了一系列探索。数值测试结果表明，它们在不同程度上提高了图 像的质量。 本文提出的算法被集成于工业c t 机成像专用软件，合作完成的有关论 文在2 0 0 0 年国际c t 和三维断面成像会议上报告，并获得了“2 0 0 0 年亚洲 c t 科技进展荣誉奖杯”。在读研期间，本人还参与了由导师主持的国防科 技项目的研究，该项目已通过国家技术鉴定和有关部门组织的验收。 本文其余各章组织如下：第二章介绍c t 成像的基本理论，包括投影理 论、r a d o n 变换和逆变换，以及离散重建算法。第三章主要介绍经典a r t 算法的基本思想、迭代公式、收敛准则和特点。第四章提出了一种分块迭代 迭代的快速a r t 算法，内容包括分块迭代方法、投影系数矩阵的计算、投 影数据的拟正交排序、收敛准则的改进等。第五章探讨提高图像质量的若干 方法，比如噪声投影数据处理、不完全数据重建、投影数据插值、低分辨率 数据重建高分辨率图像。附录介绍作者结合工业c t 机的研制，参与设计开 发工业c t 成像专用软件的有关工作。 西安理工大学硕士学位论文 第3 页 基于分块送代的快速代数重建算法研究 第二章c t 成像的基本理论 1 投影理论 c t 的图像重建，是基于射线与物质的相互作用。在衰减系数为的均 匀介质中，当强度为厶的射线穿过厚度为z 的介质时，其强度衰减，服从比 尔定律( b e e r sl a w ) 1 。 1 = 厶e x p ( 一棚 若非均匀介质，则衰减系数可看作点坐标x 和y 的函数卢( x ，y ) 【2 】， 衰减强度，必须沿吸收路经进行积分，其中是射线穿过介质的路径。 ，2 ，。e x p ( 一l p ( z ，y ) d s ( 21 ) y 7 册疆耳 匆| | 。 0 磁慕 幽1 ( x ，y ) 与物质的种类和射线的能量有关，密度大的物质对射线的吸收 率大。并且，( z ，y ) 的分布与物质分布密度函数成线性关系，将上述方程 进行变形。可以得出一新的函数，其中线积分沿射柬路径计算，可以得到 投影数据，即线积分值正( p ，0 ) 如一l n ( 毒) 2 i 施，y 陟 ( 2 2 ) 丘( p ，0 ) 为投影数据，因为可以唯一确定p ，0 ，改变p ，0 就可得 西安理工大学硕士学位论文第4 页 基于分块选代的快速代数重建算法研究 到各个方向、不同位置上的丘( p ，0 ) 的值【3 。 2 r a d o n 变换及其逆变换 由1 1 知，c t 投影数据是线性衰减系数的积分。图像重建的实质是根据 物体所有可能的投影数据，估算物体内部的密度分布。数学基础就是r a d o n 变换及其逆变换，由r a d o n 在1 9 1 7 年给出了证明 4 】。 ，j 、 一 r 蔚、 矿、 、 、 o i 、r 图2 如图2 所示，直线三的方程为 p = x c o s 0 + y s i n 0 ， 其中p 为直线到原点的距离，0 为赢线的法线与x 轴的交角。r p ，们与 o x y 平面的直线三对应。 在( 22 ) 中，令 f ( x ，力= ( * y ) ， f ( p ，0 ) = 皿( p ，0 ) 那么( 2 2 ) 变为 夕( p ，日) = i ，( x ，y ) ( 2 3 ) r a d o n 变换是一个算子，用鼾来表示，当r a d o n 变换算子作用于一函数 的时候产生另一函数：在( 23 ) 中，称于( p ，口) 是，y ) 的r a d o n 变换， 记为毗，。 西安理工大学硕士学位论文 第5 页 基于分块迭代的快速代数重建算法研究 在图2 中，记向量x = ( x ，_ y ) ，或 通过代换，( 2 3 ) 为 ，( p ，p ) = e ，( p c 。s 0 一j s i n bp s i n 口+ s c o s 0 ) d s ( 2 4 ) 设 # = ( c o s 0 , s i n 0 ) ， f 1 = ( 一s i n 0 , c o s 0 ) ( p ，r ) 为z 在坐标( 参 1 ) 坐标系中的坐标，那么 p # + f # 1 = ( p c o s 0 - t s i n 0 , p s i n 0 + t c o s o ) = ( x ，y ) = x 在上述条件下，( 2 4 ) 改为 于( p ， ) = e ，( p # + r # 1 ) a t ( 2 5 ) 由上式，直线的方程为 p 一# x = 0 根据广义函数( j 函数) 的性质，( 25 ) 也可写为。 夕( p ，f ) = ，( z ) 6 ( p f x ) a x d y ( 2 6 ) r 2 为- f 求f f l f ( x ，j ，) ，r a d o n 在1 9 1 7 年对公式( 2 6 ) 进行求逆运算，这 就是著名的r a d o n 定理，函数f ( x ，y ) 的解为 m 川= 萨- 1p e 南型学咖( 2 ， 西安理工大学硕士学位论文第6 页 口口蜀瞄 j s 一 十目口 s 1o c s p p 1 1 = x y ，、【 基于分块迭代的快速代数重建算法研究 由公式( 2 7 ) ，r a d o n 逆变换由关于投影数据的导数算子、h i l b e r t 算子以 及反投影算子复合而成，其中h i l b e r t 算子为无限区域上的奇异积分算子， 难以直接计算；导数算子对数据噪声有放大作用，反投影算子可引起伪影。 此外实际采集的数据，有时还是不完全的。因此在实际层析成像计算中， r a d o n 逆变换不便用直接数值求解，而采用离散重建算法。 c t 投影数据实际上是沿一组平行射线( 平行束) 或发散射线( 扇形束) 进行采集的，所能得到的只是有限个离散角度，有限位置上的投影数据集( 本 论文都用平行束数据) 。一个重建算法就是通过有限个数对( p ，日) 的测量 值 g 矿】( p 口) ( 根据物理量进行估计) ，来计算，的逼近值。 设定数据采集方法，并且 9 i ，】( p ，0 ) 的测量值由，个数对( p 。，目。) ， ( p ：，0 2 ) ，( p ，研) 确定，数据集为 避厂( n ，目，) ，观厂( p ：，曰：) ， 守i ，( p j ，目，) ) 定义吼，f = 9 v ( p f ，b ) ，1 i ，。在下文中用y ，表示吼。f 的测量值 ( 投影数据) ，用y 表示，维列向量，并且它的第i 个元素是y 。称为投影数 据向量。图像重建就是根据投影数据集】，来估计图像分布，。 目前常用的、发展较为成熟的、从投影重建图像的算法，基本上分为两 大类，一类是变换方法，另一类是级数展，f 法。 变换方法中最常用的是卷积反投影算法( f b p ) ，其特点是重建速度快， 空间分辨率高。但是该算法要求完全的等间隔采样数据，且重建图像伪影较 重；快速傅立叶变换算法( f f t ) 重建速度很快，但是，从极坐标网格到直 角坐标网格内插、利用最近邻域值或双线性内插时，误差较大。另外，小波 变换方法是近几年发展的一种新的算法，该算法有良好的局部性质，主要应 用于局部重建。 级数展开法中有代表性的是代数重建方法( a r t ) 。算法简单，适用于 不同的采样数据，对不完全数据亦可重建图像，也可以结合先验的知识，列 具有参数模型的实物进行成像。还有联合代数重建方法( s i r t ) ，这两种 西安理工大学硕士学位论文第7 页 基于分块迭代的快速代数重建算法研究 方法都是通过迭代来求得线性方程组的解。因此，这两种算法的缺点都是计 算量大，收敛速度慢，对存贮空间要求高，难于重建大的图像。 西安理工大学硕士学位论文第8 页 基于分块迭代的快速代数重建算法研究 第三章代数重建算法 a r t ( a l g e b r ar e c o n s t r u c t i o nt e c h n i q u e s ) 是船数重建法的一种典型重建方 法。事实上，a r t 所表示的最一般的代数迭代方法由k a c z m a r z ( 1 9 3 7 ) 提出， 来求解相容的线性代数方程组。g o r d o n 和h e r m a n 等人进一步提出用a r t 方法解大型稀疏矩阵和不等式【5 】。本章介绍这一经典算法。 1 a r t 的数学模型 a r t 的目的就是由某一断层不同方向的投影值，通过迭代解线性方程 组的方法得到物体被测区域的衰减系数的分布。衰减系数和物体的密度有正 比关系，可以反映图像密度的分布。通过图像显示技术，将密度分布转化为 图像中灰度分布，就可以得到物体该断层的图像。 在c t 中，图像区域就是重建区域，并且在某一点的图像密度值就是 在该点物质组织在某种能量下的相对衰减系数值，假设开始将物体被重建区 域离散化为一幅数字化图像，被重建的物体密度在像素内平均取值，那么， 我们就可以得到每一点的图像值。 如图3 所示，该图像是一t l x ”的数字化图像，每一网格表示一像素， 将象素从1 到，进行编号，并定义 州，= 括蠹鬈猱筹嬲i 于是，图像f ( x ，y ) 的n 的数字化图像f ( x ，y ) ， 为 7 ( x ，y ) = z 。b ( x ，y ) 其中一是厂在筇，个像素内的平均值。 西安理工大学硕士学位论文 第9 页 基于分块迭代的快速代数重建算法研究 上式可间记为 7 = z x 户， ( 3 1 ) j - i 基图像被选定后，就能用基图像的线性组合来表达任何图像，。记 x = ( z 。，x ：，巧) ，称为图像向量，因此，将求解图像近似值，的问题，转 化为，寻求一个图像向量x ，使得7 尽可能接近，。 由上一章r a d o n 变换的定义，可得 m ，= m 。7 = x j o t q ( 3 2 ) j - i 6 ，是用户所定义的函数，通常虢，b j 能用解析方法进行计算。如下图所示， 9 i ，6 ，就是x r a y 射线源检测器对的第f 个位置上的直线( 既由原点到与 t y 轴夹角为饥的直线上距离为p ，的直线。，见图3 所示) 和第，个像素 交线的长度a b 。 t 啦蠢 。 、i 、i ￥ 芽 0 、 x 墉 办 7 1 q j 幽3 以0 表示孵。b ，的计算值，即 0 = 辨，b 若令y ，为投影数据向量的第，个分量，( 3 2 ) 改为 j y 。= 白_ ( 3 3 ) 西安理工大学硕士学位论文 第1 0 页 基于分块选代的快速代数重建算法研究 在公式( 3 3 ) 中，( i = 1 , 2 ，3 ，；j = 1 , 2 3 一，) ，并且有 y 为第i 条射线的投影值； x ，为第j 个网格的图像值： “为第f 条射线经过第，个网格的射线长度。 其中，为射线投影数，j 为像素数，公式( 3 3 ) 可表示为矩阵形式 y = r x ( 3 4 ) 其中，r 为i x l 阶投影矩阵，z 为待求的j x l 阶图像矩阵，五为x j 阶 系数矩阵。 方程( 3 4 ) 为具有，个未知量的，个联立方程组，如对矩阵直接求逆， 也能得到一个解，也就是 并= r 。1 y( 3 5 ) 但是，直接求解上述线性方程组，由于以下原因而不采用，和，数 目过大，导致运算量很大；而且线性方程组不适定，不存在唯一解；如果投 影数据含有噪声或量化误差，方程可能无解。注意到r 是典型的稀疏矩阵， 大多数矗= 0 ，因此，通常通过迭代的方法，求方程解。 2 a r t 的迭代方法： a r t 的基本过程是：给定重建区域的一个图像初值，再将所得的投影值 残差沿其射线方向均匀的反投射回去，不断的对图像进行修正，直到满足所 需要求，然后结束迭代过程。 首先对未知数z ，设定一组初值x ；“，对每一条射线( 如第i 条) ，经过迭 代计算后，求得一组修正值g ，使 x j 为z j + c f ，( 扛1 , 2 3 ，；j = 1 , 2 ，3 ，一，l ，) 并且，满足下面表达式 j o ( x ，+ q ) = y ， 西安理工大学硕士学位论文第1 1 页 基于分块选代的快速代数重建算法研究 那么，g 可以计算出 其中，p 。和a 可有以下公式给出 = 勺x ， a = 儿一p ( 3 6 ) 公式( 3 6 ) 给出了第i 条射线通过的第，各像素的修正值c , j 的迭代公 式，这样就可以对x ，值进行修正。事实上，公式( 3 6 ) 可写为 从而 上式可改写为向量形式 y 。一k ，x 尸 e ( 。 x r ) - x ，+ 嘴 其中，五( 。为上式中的未知向量第次的迭代值 t ，为系数矩阵r 的第i ；个分量( j 维列向量) ； y 。为第i 。条射线的实测投影值； i 。2 为i ，的转置向量； 是i ，与z 。的内积，既第次的计算投影值。 由上述的迭代公式可知，a r t 方法在逐次运算过程中产生一系列的向 量z ( ”，z ( ”，一( ”，当足够大时，这个向量序列收敛于所要求的估计值 西安理工夫学硕士学位论文第1 2 页 士矽 = g 一o i | x j j x 基于分块迭代的快速代数重建算法研究 z + 。 a r t 的具体迭代步骤如如此下( 首先对i = l 进行修正) 首先对各像素值赋初值 x = x ? ，u = l ，2 ，3 ，t ，) 计算第i 个方程的估计值，即 j p 。= x ；o ) j 4 计算残差 i = y i p ? 计算第f 个方程的修正值 y 4 - 名- x ，值进行第一次修正( x ，+ c f ) ，由于每一条射线只是通过 像素的最少一部分，只对射线通过的那部分像素进行修正，而其它未通过的 那部分像素，保持不变： 将第一个方程修正的值x ，代入到第二个方程求p ：，在重复以上 一等步骤，既由第二个方程对各x 值进行修正，如此循环，直到第，个 方程； 经过从i 到，各方程的修正，这是得到的各x ，值叫做完成了第一轮 迭代。一般经过第k 轮迭代以后，各x ，的值和前一轮的值满足给定的阀值 ( 由收敛准则决定) ，就可以认为各z 的值达到了要求，否则，还要进行下 一轮迭代，直到符合收敛要求的结果为止。 对不同的实际图像重建问题又发展了一系列的a r t 方法。为了加快收 敛速度，抑制噪声，提出了阻尼a r t 法。 它的主要迭代公式为： 西安理工大学硕士学位论文第1 3 页 士驴 1 i 咯 基于分块选代的快速代数重建算法研究 _ = 5 ( r + 1 ) = z 肼+ ) x一 ”为阻尼因子，当0 ( 办” 2 时上式收敛e 实际计算中，一般日 彤” l 。其实，阻尼因子的作用是按照迭代的次数合理地分配残差、压 迭代时数据的振荡、防止上述迭代公式发散。 3 a r t 的收敛准则： a 】u 方法迭代运算时，只能采用有限的次数，因此，根据迭代的精度j 求，确定收敛准则。 一般采用的收敛准则有以下几种： 定义残差： 一d 。= 去j 姜c 缈 mv 智、。 其中，缈，= y ，一 。 定义方差： 矿。= 寺( x ，- x ) 2 智 定义熵： 拈一击n 喜审t n 争1 智、x o 上述几种迭代准则，随着t 的增大，d 趋于零，v 趋于极小，s i 于极大。方差 ，r 和熵s r 反映了解的平滑程度。方差极小，熵极大求出的j 一个平滑图像，但是在实际的重建过程中，既要求残差小，又要求方差小， 就要控制迭代次数，也就是要取一个折衷解，如果要用单精度浮点数组作； 精度标准，一般认为d 1 0 4 就可以认为达到要求。因此，按照通常的* 度要求，当d 。 s 停止迭代( s 为给定的精度) 。 西安理工大学硕士学位论文 第1 4 页 基于分块迭代的快速代数重建算法研究 4 a r t 的特点 首先，a r t 方法是一种逐次迭代的算法，它避免直接系数矩阵求逆的 运算，减少了运算量。 其次，在每一个迭代计算中，只用到系数矩阵r 的一行元素，可以节 省大量存贮空间。 再之，可直接对超定方程组求解，而不需转化为正则方程，当方程组 相容时，可以得到精确解。 此外，可以将投影值和残差均匀地反投影，不至于使误差集中在一起 造成解的畸变；a r t 方法还可适用于不同方式的投影数据，对不完全数据 也能进行图象重建。 5 影响a r t 收敛速度的因素 上一节介绍了a r t 的优点，但是，a r t 也有不足之处。由于a r t 是解大 型线性方程组，系数矩阵是稀疏矩阵，而且线性方程组是不适定的，因此， 在迭代过程中，运算量大，对存贮空间要求高，收敛速度慢。本节针对这些 缺点，对影响a r t 收敛速度的关键因素，进行了分析。 第，初值的选取。因为，初值的选取对收敛的速度有很大的影响， 如果它和图像值很接近，那么，根据迭代准则，初值很快收敛于图像值。 第二，系数矩阵的设计和运算。因为，在每次迭代的修正过程中，在 计算投影估计值时，只用到系数矩阵的非零元素，而非零元素占系数矩阵的 很少一部分。因此，设计有效的系数矩阵生成算法非常必要。 第三，a r t 迭代公式的设计。算法的迭代公式是晟核心的内容，一个好 的迭代公式直接影响到算法的收敛速度。因此，人们想出各种各样的迭代方 法来加快收敛速度，最典型的算法有以下几种，正交化重建方法( 系数矩阵 正交化) ，采用非线性最小二乘法，还有多准则迭代收敛法等。这些改进目 的都是用来加快a r t 算法的收敛速度，但对较大的图象，上述难以奏效。 西安理工大学硕士学位论文第1 5 页 基于分块迭代的快速代数重建算法研究 第四，迭代收敛准则的设计。有效的停机准则可以避免不必要的下轮迭 代，加快收敛的速度。 针对上述因素，对a r t 算法进行改进，才能加快收敛速度，可重建高质 量的图象。 西安理工大学硕士学位论文第1 6 页 基于分块选代的快速代数重建算法研究 第四章分块迭代的快速代数重建方法m 1 上一章介绍了a r t 算法的基本思想、迭代公式和特点，并分析了影响a r t 算法收敛速度的主要因素。将在此基础上提出一种基于分块迭代的快速代数 重建算法。通过大量的模拟数据和实采数据的测试，表明该算法收敛速度快， 重建图像伪影轻，并有较高的空间分辨率和密度分辨率。 在经典a r t 算法迭代过程中，开始就对图像所有的像素值进行修正， 事实上，这是不必要的。因为，迭代时的图像初值一般是零或者是图像的某 种平均值，和实际图像值相差很大。开始就对每个像素修正，这种方法并不 理想。因此，本文采用了基于分块的a r t 算法，对重建的图像分块，把每 块图像作为一像素，进行修正，很明显，这种方法减少了重建像素的个数， 缩短重建时间。该方法最显著的特点是分块图像大小可以改变，通过修正后 的图像值来确定分块图像的大小。当计算图像值和原图像值较接近时，分块 图像较小，反之，较大。整个分块迭代是图像逐渐细化、收敛逐渐加快的过 程，这是此算法关键点之一。 1 1 分块基本思想 首先，将图像重建过程分级，每一级图像分成不同的块，相邻两级图像 块之间有一定的关系，图像块大小由级数决定。然后，确定一初始图像，赋 予初值，对投影数据抽样，开始迭代。 首次迭代结束后，图像块数为初始图像块数的四倍，将上次迭代的值 分块赋值( 每一图像块一分为四) ，相邻的像素再进行线性插值，以避免图 像的块效应，再按照图像块大小对投影数据抽样，进行迭代修正。 随着级数增加，下一级图像块数为上一级图像块数的四倍，按照分块 递进关系进行迭代，当级数到最后一级时，分块迭代结束，按照收敛准则， 如果图像值满足一定精度要求，全部迭代结束，否则，进行下一轮分块迭代。 西安理工大学硕士学位论文第1 7 页 基于分块选代的快速代数重建算法研究 1 2 分块迭代格式 设重建图像为n n 个像素( 一般为2 的倍数，这样处理比较方便) ， 设置分块级数为h 。 首次图像块为 1 1 = n 2 一1 u 2 “” 每块图像为2 ( “1 】2c “1 像素大小，图像向量为 x o k ( x f l l x 一，x 孙) ，( m 1 = n 2 ”4 。1 n 2 ”“) ， 然后对图像向量赋予初值，进行运算。 第二级时，图像块数是首次迭代的四倍，为n 2 “2 2 ( “”，图像 向量为首次图像向量分量个数的四倍，为 x 。k ( x 黔x 黔x 茹) ，( m “) - 2 。“1 x n 2 ”。1 ) ， 将首次图像向量迭代的结果赋予第二级，由于每一级向量大小不一致， 所以，有一种特殊的赋值关系，不同级图像块之间的赋值在以下用图例说明。 重复以上的迭代方法，进行第三级、第四级、以及第k 级的迭代，直 到所有的级数迭代完为止。 一般地，图像之间的关系为： 第k 次分块后，图像块数为 1 = n 2 “2 ( 6 “， 图像向量为 “= ( z p ，z 一( k - o ，) ，( m “= n 2 “一1 n 2 “一1 ) ， 第k + 1 次分块后，图像块数为 “= 2 4 删u 2 “”， 图像向量为 五”= ( x x ，x ) ，( m “1 - n 2 ”“1 n 2 “) 图像块数之问的递进关系为 西安理工大学硕士学位论文 第1 8 页 基于分块迭代的快速代数重建算法研究 ，( o + 1 ) = 4 ，( “ 不同级的图像之间赋值时，按照图像块进行，赋值结束后，再进行排序 ( 和系数矩阵对应) 按照上述的公式，有如下赋值关系： 设第k 次第s 行第t 列图像块为x ? ， i = s x u z “。+ f 那么，上述元素所对应的第t + 1 次的元素为x ，”， j = ( 2 xs - 1 ) x 2 “1 + ( 2 t 一1 ) ， ( 2 x s 一1 ) 2 “1 + ( 2 f ) ， ( 2 x s ) x n 2 6 + 1 + ( 2 f ) ， ( 2 j ) 2 “1 + ( 2 t 1 ) 定义x ! ：为图像矩阵在第次分块迭代时，第s 行第f 列图像块值。 由于第女+ 1 次的四个像素值都一样，为了防止方块效应，将其中三个 元素进行邻域插值，图4 、图5 说明了赋值和插值的过程。 扦 第2 t 刊 x 。 x ，e ， 叶t t - lx 2 一i 毪!一l心汝 2 12 t心 潜笏骚【2 x 2，2 1一#彩绣 趁 角 2 s 行 第k 次弟k + l 谈 图4图5 图4 所示第k 次第s 行第f 列元素x ! ：，图5 中，第k + 1 一、“。，- - u 一“2 ( k n + d ，、 x “2 s “- l 、, 2 ，、x 盟l 。、x 芸嚣，的值和第k 次的元素x 砦的值相等，为了使图像平 滑，上述任何一个元素要和相邻的元素进行插值，比如x 器j 。和相邻元素 x 2 ：，、z i ：。、x 剀。进行插值，将插值后的值赋予x 嬲| 。然后进 行第k 次迭代。 西安理工大学硕士学位论文 第1 9 页 基于分块迭代的快速代数重建算法研究 按照上述的递推关系，每一级分块迭代后，当图像的值满足一定的精 度要求后，进行下一级分块迭代。 图6 、图7 、图8 、图9 是用分块的代数重建算法重建的每一级的图像。 图6 图8 图7 图9 从实验的结果看，分块算法收敛速度快，图像质量比较高。 在每一级的分块迭代中。根据图像分块的太小、投影数据的数据量之间 西安理工大学硕士学位论文 第2 0 页 基于分块选代的快速代数重建算法研究 的关系进行投影数据抽样，目的是为了加快收敛速度和提高图像的分辨率。 数据的抽样是基于相关性最小的原则，运用了两种方法，一是对投影数据的 方向进行抽样，二是对每一投影方向的投影数据进行抽样，下面是数据抽样 的理论分析【9 】。 由公式( 3 3 ) ，对于任意投影数据，都可以表示为一个线性方程，考虑 以下两个方程，它们是由不同的投影数据产生的。 。靠寸 札叫， 6 0 ，= ( 1 1 ，1 2 ，3 ，) o j2 u j ? 1 r i ? 2 r i 3 ，”j i ? n ) c o s 6 f 叫= 。j l 卯i l l 田j l ， 其中，卯，是线性方程的系数向量，c o s l z u 是印，出的相关系数。 从几何的角度，如果把每一个线性方程看作一个超平面，那么， c o s t z ，就是超平面单位法向量的内积，当两个向量正交时，很明显，相关 系数为零，而线性相关时，相关系数最大。但是，在用迭代法求解线性方程 组时，迭代次数与相关系数有很大关系，相关系数越小，迭代次数越少，给 定的初值很快逼近原方程的解。 在投影数据抽样时，除满足投影数据之间的相关性比较小之外，还必 须保证每一像素都要被修正，因此，抽样后相邻的投影数据必须有一定的关 系( 相邻射线间距必须小于网格长度) ，否则，可能导致少数像素修正的次 数不够，而使这个区域的图像变模糊，影响图像的分辨率。 投影数据抽样后，数据量明显减少，因此，收敛速度得到了提高。由 于在不同级的重建过程中，数据抽样时，抽样比率依赖于分块图像的大小， 图像块变小时，抽样比率变小。通过实验测试，利用抽样数据进行迭代，在 开始迭代时，只需几次就可以重建物体断层图像的轮廓，重建时间明显缩短， 西安理工大学硕士学位论文第2 l 页 基于分块选代的快速代数重建算法研究 图像质量有所改善，收敛速度明显提高。 2 系数矩阵计算 根据a r t 算法的迭代公式，每次迭代都要计算投影估计值，因此，系 数矩阵的计算直接影响到该算法收敛的速度。 1 1 系数矩阵模型 假设要重建物体的图像区域为n xn 个像素( 一般重建区域选取物体 实际区域的最大半径的外围正方形) ，那么，对任一投影数据，在迭代时系 数矩阵大小为n n ，因此，墓建的系数矩阵模型应和重建图像区域相对应， 为n n 个网格，如图1 0 所示，其中射线( 产生投影数据) 与网格相交时， 系数矩阵元素的值为交线长度，否则，取零。实际上，任何经过该重建区域 的射线，与网格相交的个数不大于2 n 个，每个网格交线长度为系数矩阵 的非零元，非零元素的计算直接影响到系数矩阵的计算。但是，要计算非零 元的个数，必须设计很复杂搜索算法，这样会增加运算的时间，降低求系数 矩阵的效率，因此，只对系数矩阵的非零元的计算设计优化的算法。 图1 0 由公式( 3 4 ) ，可知，x 是2 维列向量，而，代表重建区域图像的像 素个数，并且图像区域是一个方阵，为了计算的方便，对图像区域进行排序， 对应关系为： 西安理工大学硕士学位论文 第2 2 页 基于分块迭代的快速代致重建算法研究 ；置，( i = 1 , 2 ，n ，j = l ，2 ，) ； k = j nq - f ： 而系数矩阵置是m x n 2 维( m 代表投影数据个数) ，r 的每一行向量 的元素应和图像向量的元素一一对应，所以，在计算出每条射线与所有网格 交线的长度后，也要进行转化排序， 其中， 。= 毋( i - 1 , 2 ，n ，= 1 2 n ，= 1 , 2 ，肼) ，k = ，+ i ， o 代表第z 条射线在第k 个位置的值，其中k 是系数矩阵元素的序号，下一 节是o 的计算方法。 1 2 系数矩阵非零元的计算 系数矩阵的计算最关键的部分，是非零元。实际上射线与网格交线长度 的计算。传统的逐行扫描搜索算法，此算法简单，计算量大计算方法是先 按照模型建立直角坐标系。确定射线的方程，把所有的网格定位，然后从第 1 行到第行进行计算，每一行采用单向递进搜索的方法，求出与网格相交 的长度和序号。 本文所采用的是另一种沿射线进行搜索的算法，采用了二叉树遍历算法 的思想，先求出射线与重建区域相交的一条边，计算出第一个网格的长度和 它的序号，再采用取模( 减少加法和乘法运算，尽量利用优化的系统函数) 和内积运算( 避免求线段交点的运算) ，再判断下次相交的网格，然后计算 下一网格长度和序号，按照此算法循环进行，直到与下一条边相交时( 射线 和图像区域只能和两条边相交，或者和某一条边重合) ，整个运算结束。 在存贮网格长度和序号时，我们设计了一种数据结构，这个结构包括两 个元素，网格长度( 系数矩阵的值) 和序号，把序号作为结构的元素是为快 速找到投影值进行修正，由于每条射线与重建区域相交的网格数不大于 2 n 个，所以此结构的存贮空间不大于2 x n 个，由相交的网格数决定，目 的是迭代时不增加额外的运算量。 系数矩阵的值有两种计算方法，一种是当射线与网格相交时，取值为l ， 不相交时，取值为零，是一级精度。另一种方法是计算与网格交线的眨度， 西安理工大学硕士学位论文第2 3 页 基于分块选代帕快速代数重建算法研究 是二级精度。当重建图像是低频部分比较多时，效果一样，但当高频部分占 据多数时，后一种方法比较理想。 在上述算法的程序实现中，对代码进行了优化处理，这样，就加快了运 算速度，减少求系数矩阵所用的时间，在重建5 1 2 x 5 1 2 图像，投影数据为 5 7 6 1 8 8 1 时，直接用此算法计算系数矩阵，第一种方法计算时间为4 秒左 右，第二种方法大约8 秒左右，使系数矩阵的实时运算成为可能，下图是用 不同的系数矩阵，对原木采集的数据重建的图像。 图1 1 图1 2 图1 1 是用第一种方法重建的图像，图1 2 是用第二种方法重建的图像。 从图像上可以看出，图1 2 的原木年轮的纹理和图1 1 相比较，没有图1 1