（生物医学工程专业论文）ct迭代重建算法的加速方法研究.pdf

东南大学硕士学位论文 由实验结果可以看到，第一，与传统的系统矩阵计算方式相比，基于n u f f t 的傅里叶方法 在前向投影和反投影算子的实现上有明显的速度优势。第二，在采用优化k a i s e r b e s s e l 插值器 以及m i n m a x 插值器的情况下，其计算代价仅接近于传统的双线性插值傅里叶方法，而投影误 差则明显降低，其归一化均方误差和归一化绝对误差均等于或小于标准的重投影方式( 即系统 矩阵方法) ，从主观图像结果上也已可以看到令人满意的表现。因此本文认为采用基于n u f f t 的傅里叶方法进行重投影来实现对迭代中的重投影步骤加速是可行的，能够以小于标准的系统 矩阵重投影方法的时间代价来保障重投影精度，适用于不需要进行分块运算的任何c t 迭代重建 算法。 关键词：抗噪性c t 成像，迭代重建，加速算法，投影估计，高斯统计模型，傅里叶重投影， n u f f t ，插值优化卷积核。 i i 东南大学硕士学位论文 t i t l e ： a b s t r a c t a c c e l e r a t i n g m e t h o d so fc o m p u t e r i z e d t o m o g r a p h i c i t e r a t i v e r e c o n s t r u c t i o n a u t h o r ：l ih u i s u p e r v i s o r ：p r o f b a ox u d o n g s c h o o l ：s o u t h e a s tu n i v e r s i t y af a s tr e s i s t - n o i s er e c o n s t r u c t i o nm e t h o dw h i c hc o m b i n e st h es t a t i s t i c si t e r a t i v ea l g o r i t h mw i t h t h ec l a s s i c a lf i l t e r e db a c kp r o j e c t i o n ( f b p ) a l g o r i t h m ，a c c o m p a n i e db yaf a s tr e p r o j e c t i o nm e t h o d w h i c hc a nb eu s e dt oa c c e l e r a t et h ee s s e n t i a la n dr e p e a t e ds t e p si na l li t e r a t i v ea l g o r i t h m sf o rc t r e c o n s t r u c t i o no n l ye x c e p tb l o c ki t e r a t i v eo n e s ，a r em a i n l yd i s c u s s e di nt h i sd i s s e r t a t i o n t h ee x p e r i m e n t sa r em a d et ot h e s et w om e t h o d sb yu s i n gd e s i g n e ds i m u l a t i o nd a t aa c c o r d i n gt o w h i c ht h ea n a l y s i so nr e c o n s t r u c t i o ns p e e d ，n o r m a li z e ds q u a r ee r r o ra n da b s o l u t ee r r o ra sw e l la st h e a d v a n t a g e sa n dd i s a d v a n t a g e so ft h e mb o t ha r eg i v e n t h es t a t i s t i c a lm e t h o d si nc ti m a g i n ga r ec a p a b l eo fp r o d u c eh i g h e rq u a l i t yi m a g e st h r o u g h m o d e l i n gt h ei m a g i n gs y s t e mw h i l ee x p l i c i t l yt a k i n gi n t oa c c o u n tt h en o i s e ，a n dt h e ya r ca b l et o i n c o r p o r a t eap r i o r ii n f o r m a t i o na b o u tt h es o l u t i o nt or e g u l a r i z et h ep r o b l e m c o n s i d e r i n gb o t ht h e a d v a n t a g e so ft h a ta n dt h ee f f i c i e n c yo ft h ec l a s s i c a lf b pa l g o r i t h m ，t h em e t h o do fe s t i m a t i n gt h ei d e a l s i n o g r a m ( p r o j e c t i o n si na l ld i r e c t i o n s ) f o ri n v e r s er a d o nt r a n s f o r mb yf b pi sd i s c u s s e da n da p p r o v e d i nc h a p t e r3 t h ea p p r o a c hi sb a s e do nag a u s s i a nn o i s em o d e lw h i c hd i f f e r sf r o mt h et r a d i t i o n a lp o i s s o n m o d e l t h ep w l so rm a pc o s tf u n c t i o ni s o p t i m i z e dt h r o u g ht h ei c ma l g o r i t h mw i t hah i g h c o n v e r g e n c es p e e d t h er e s u l t i n ge s t i m a t e ds i n o g r a m sa n dr e c o n s t r u c t e di m a g e si ni t e r a t i o n st o g e t h e r w i t ha c c o r d i n ga n a l y s i so nt h ee r r o r sa r eg i v e n ，f o l l o w e db yad i s c u s s i o no nt h em e r i t sa n dl i m i t so f t h ea p p r o a c h af o u r i e r - b a s e dr e - p r o j e c t i o n ( o rf o r w a r dp r o j e c t i o n ) m e t h o di su s e dt or e d u c ec o m p u t a t i o nf o r t h eb a s i cs t e p si nt h ei t e r a t i v er e c o n s t r u c t i o na l g o r i t h m sw h i l es t i ll p r o v i d i n ga c c u r a c yc o m p a r a b l et o t h et r a d i t i o n a ls p a c eb a s e dm e t h o dw h i c hr e q u i r e si ne a c h i t e r a t i o no n eo rm o r em a t r i x - v e c t o r m u l t i p l i c a t i o ni n v o l v i n gt h es y s t e mm a t r i xt h a ti sh u g ee v e nf o rs m a l lr e c o n s t r u c t i o n s t h ea p p r o a c he s s e n t i a l l yc o n v e y st h ep r o b l e mt oo n eo fs p e c t r a le x t r a p o l a t i o nv i at h ef o u r i e r s l i c et h e o r y t od e a lw i t ht h ed a t an o tf a l l i n go nac a r t e s i a ng r i dt h en u f f t t e c h n i q u ei n v o l v i n g c o n v o l u t i o ni n t e r p o l a t i o ni sa d o p t e dt or e d u c et h er e - p r o j e c t i o ne r r o r t h ef o u r i e r - b a s e dr e p r o j e c t i o n i i i 东南大学硕士学位论文 u s i n go p t i m i z e dk a i s e r - b e s s e lk e r n e la l o n gw i t hm i n - m a xk e r n e li sc o m p a r e dt ot h es p a c em e t h o di n s p e e da n de r r o r s k e yw o r d s ：c ti t e r a t i v er e c o n s t r u c t i o n ，f a s tm e t h o d ，s i n o g r a me s t i m a t i o n ，n o nu n i f o r mf a s t f o u r i e rt r a n s f o r m ( n u f f t ) ，f o u r i e rb a s e dr e p r o j e c t i o n 东南大学学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。 尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他入已经发表或撰写过 的研究成果，也不包含为获得东南大学或其它教育机构的学位或证书而使用过的材料。与我 一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 研究生签名： 东南大学学位论文使用授权声明 东南大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位论文的复印 件和电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸 质论文的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅，可以公布( 包 括刊登) 论文的全部或部分内容。论文的公布( 包括刊登) 授权东南大学研究生院办理。 研究生签名： 第章绪论 第一章绪论 1 1 图像重建与c t 成像技术 由物体的二维截面或断面向该平面内的各个方向作投影，可获得一系列一维投影函数。由 这些一维投影函数来重建该二维截面，称为断层图像重建。这一随着计算机技术的进步而获得 发展和应h j 的图像处理技术，称为计算机断层成像术( c o m p u t e dt o m o g r a p h y 或c o m p u t e r i z e d t o m o g r a p h y ) ，简称c t 技术。“t o m o g r a p h y ”一词来源于希腊语的“t o m o g r a p h i ”， “t o m o ”意 为切片，“g r a p h i ”意为记录。 这一技术最典型的应用是医学上的x 线c r 成像( 用x 射线作为辐射源) ，用于获取人体头 颅、心肺、腹部等内部器官的二维断层图像( 亦称断层摄影术) 。它对x 射线放射诊断是一个重 大突破，冈为c t 扫描检查人体所获得的图像，是真正的横断面或冠状断面的图像，这些图像没 有不同组织器官，病变等的影像相互重叠，又能提供受检切面组织器官和病灶等的解剖细节， 可对病变或组织器官的形态、大小、部位、解剖邻属关系等，做出准确和“立体”的判断。 1 1 1 图像重建问题的提出 1 9 1 7 年，德国数学家r a d o n 【1 】在一篇关丁解决引力场问题的数学论文中，描述了这样一个 问题：设一个具有某种物理性质的二维分布函数f ( x ，y ) 定义在平面域q 内，p ( r ，p ) 表示f ( x ，) 沿 着m r 和0 定义的一条直线的线积分，f ( x ，剪) 的所有线积分的集合为： p ( r ，0 ) = f f ，( z ，y ) 6 ( z c o s o + y s i n o r ) d x d y ( 1 1 ) 其中r ( 一( 2 0 ，。o ) ，0 1 0 ，7 r ) ，6 ( ) 是d i r a cd e l t a 函数。于是，( z ，剪) 可有它的所有线积分的集合 唯一地确定。对于一个固定的0 ，p ( r ，p ) 称为f ( x ，可) 沿0 方向在r 轴上的一维投影。 r a d o n 所提出的问题，其实就是要从包含有目标信息的投影数据得到目标的二维分布函数 ( 图像) ，这一问题被称为从投影重建图像，简称图像重建。在他提出这一问题之后有4 0 年的 沉寂，到1 9 5 7 年才有一位美国物理学家在哈佛大学的电子加速实验室里，开始研究图像重建理 论，并于1 9 6 3 年首次提出了一种多方向投影重建的代数计算方法i z j 。 从那时起，由投影重建图像问题的应用和相关讨论，以不同的形式反复地并且往往是独立 地出现在科学和技术的众多领域，比如射电天文学，诊断医学，电子显微技术，工业无损检测， 光学相干测量，合成孔径雷达和地球物理勘探等 3 1 。这些看起来各不相同的应用有着共同的数学 基本原理，那就是r a d o n 在他的论文中所描述的问题。 1 1 2c t 成像技术及应用 c t 成像技术能将人体中某一薄层中的组织分布情况，通过射线对该薄层的扫描、检测器对 透射信息的采集、计算机对数据的处理，并利用可视化技术在显示器或其他介质上显示出来， 在辅助医生诊断疾病、指导手术计划等方面起到有效的作用。 世界上第一台商用c t 由英国e m l 公司的电子工程师g n h o u n s f i e l d 发明，他与一位神经 放射学专家合作，研制成了可观察头颅内部结构的c t 机，并y - 1 9 7 2 年公开报道了其成果1 4 j 。 h o u n s f i e l d 因此获得了1 9 7 9 年的诺贝尔医学生理学奖；这一奖项同时还授予了c o r m a c k ，正是 东南人学硕i j 学位论文 他在1 9 6 3 年的文章中，不仅对诊断医学中断层图像重建的可行性进行了分析，而且还利用 f o u r i e r 级数理论，给出了一种由无限多个观测角度的投影值米求解目标切片的精确解析公式1 2 】。 经过儿十年的发展，c t 技术在某些方面已相当成熟，设备结构经过了数代变化，扫描速度 和重建分辨率逐渐提高，现代c t 机人多具备弧毫米极的二维分辨率，螺旋式和锥形束超高速 c t 则惊人地做剑了对慢运动组织进行高分辨率的检查，并提供人体组织的三维图像【5 1 。除了 x - c t 以外，还发展出了其他形式的c t ，如单光子发射c t ( s p e c t ) 、正电子发射c t ( p e t ) 、 核磁共振c t ( n m r c t ) 等均已付诸临床应用，透射c t ( t c t ) 、反射c t ( r c t ) 、超声c t 、微 波c t 的研究也取得了极人进展，这些c t 在投影数据的获取上有所不同，重建细节也有一定差 异，而在功能上相互补充。本文的讨论是基于x c t 的。 临床研究表明，c t 在心脏动态成像方面，和在对人体内部器官进行人规模普查( m a s s i v e s c r e e n i n g ) 方面是有效的手段，在肺癌早期检测、乳腺三维检测( 3 dm a m m o g r a p h y ) 、结肠内 息肉检测( 虚拟结肠镜，v i r t u a lc o l o n o s c o p y ) 等方面也具有很高的应用价值 6 1 。美国早期肺癌行 动计划( e l c a p ) 对1 0 0 0 名6 0 岁以上无症状志愿者进行的每年一次c t 普查表明，采用低剂 量c t 肺部扫描与薄层靶扫描相结合的方式，c t 可比常规x 线胸片普查多发现8 倍左右的肺癌， 特别是发现了更多的早期肺癌及非钙化结：宵【7j 。 与此同时，医学可视化技术的发展使得c t 的应用进一步扩展到临床诊断与治疗的各个方 面，其中基于c t 三维重建技术的放疗计划模拟系统、虚拟内窥镜系统、以及三维外科手术计划、 引导系统已在临床上获得了应用。 重建算法是c t 技术中的一个关键问题，它的实质是一个由采集到的数据求解图像矩阵中的 象素值，然后重新构造图像的过程。目前已得到广泛应用的是滤波反投影( f i l t e r e db a c k p r o j e c t i o n ， f b p ) 成像法，它的基本算法容易通过软硬件实现，且在投影数据质量高的情况下，能在较短的 时间内重建出准确清晰的图像，因而倍受青睐。 1 1 3 关注抗噪性重建算法的加速 一项统计预测表明，当人活至7 5 岁时，由于x 线辐射造成的患癌症的累积风险率在 0 6 1 8 之间。世界卫生组织n ，t o ) 和国际放射委员会( i c r p ) 以及国际医学物理组织( i o m p ) 曾制定了医疗照射质量保证和质量控制标准，主张以最小的代价和剂量获得最佳的诊断效果。 对于典型的脑部c t 研究而言，一次扫描产生的辐射剂量约在3r a d s ( 辐射的基本单位) 左右，而 普通人群的辐射剂量安全限制为5r a d s 左右，一年内两次脑部扫描就会超过这一界限i 引。过高 的辐射剂量极大地限制了c t 及上述衍生系统的进一步临床应用，特别是作为常规检查手段及在 妇女、儿童中的应用。 同时，c t 图像的质量又与辐射剂量密切相关，在其他条件等同的情况下，辐射剂量越高， 图像质簧越好。低剂量c t 的噪声比高剂鼙情况下严重。而在x 射线管电压( 或k v p 值) 及系 统其他设置一定的情况下，c t 扫描的辐射剂量主要由通过x 射线管的电流决定，目前所用c t 设备需要较高的电流通过x 射线管( 人部分临床研究所用电流一般超过2 0 0 m a ) ，以获得高质 量的诊断图像p j 。 若在c t 扫描中降低辐射剂量后，成像质量也明显下降，那就增加了对图像进行分析处理的 难度，会给临床诊断带来严重影响。在这种情况下，基于解析变换求闭合解的f b p 算法得到的 图像就难以令人满意了。 迭代成像技术基于与解析重建方法完全不同的数学模型，比如采用最大似然估计和晟小二 乘估计的统计迭代方法。一般地，在迭代中，根据由估计图像计算得到的重投影与实际投影之 间的误差，逐步更新图像。它不要求标准的几何模型，可以在不规则采样和数据缺失的情况下 2 第一章绪论 重建出图像。最有吸引力的是，能将实际系统的多个冈素纳入考虑，建立更精确的模型，来矫 正一些不理想的系统响应，如探测器延时、非线性容积误差、散射和射束硬化等，从而改善图 像质量【10 1 。迭代算法的求解采取估计与逐步逼近的思路解决问题，抗噪声能力强，比较灵活， 因而得到了越来越多的关注。 然而，迭代法中多次的前向( 重) 投影( f o r w a r dp r o j e c t i o n 或r e p r o j e c t i o n ) 和反投影步骤 的计算代价较高，成像速度比解析法慢很多，目前只应h i 于m r 稿lp e t 等信号中噪声很强、需 要长的信号采集时间的成像仪器中。对于c t 来说，它的一个显著特点就是其快速的数据采集， 可以在1 秒钟内采集多幅图像的数据，需要实时处理的数据量庞大，尽管为了提高迭代重建算 法的效率提出了很多的方法，现在已经可以在儿分钟内重建出一幅图像，但仍然显得过_ 丁缓慢， 无法与c t 的扫描速度相匹配，需要进一步的提高计算速度。 1 1 4 重建算法的评价标准 对一种图像重建算法，评价主要在两方面：一是它能达到的成果，即重建图像的质量，二 是它所需要的代价，即算法的可行性。 衡量重建图像的质量一般以三个方面作为标准：图像的空间分辨率；图像对比度；伪影情 况。比如医生对人体c t 图的主观感觉，也大致与这三个因素有关。在实际的医疗和工业应用中， 投影重建图像的完全准确的结果是未知的，没有参照对比，主要是以人的主观意识为评价标准。 但是对算法进行评价时，可以采用计算机数据模型获得模拟的投影数据，这样就可以做一些客 观指标的对比，衡量原始图像和重建结果之间的吻合度。另外还可以用已知的物体进行数据采 集，比如制作简单的实物模型，在实际成像装置中获取观测数据。 采用实物模型的好处是能够反映实际应用环境下各种合理因素的影响，坏处是模型本身的 材质等多余影响有时无法剔除，有时模型的灰度分布、形状和频率范围上的特殊性会使实验结 果具有迷惑性，另外也要考虑成像装置之间的差别，不利于比较在不同硬件环境下进行测试的 不同的算法。 模拟数据测试的好处首先是简单方便，也比较灵活，比如可以任意选取成像装置的参数， 采用广泛使用的机器型号参数等。而且重建结果的比较对象也相当精确，便于客观评价指标的 计算。但是观测数据也就是投影的形成过程必须合理，能反映客观物理过程，否则就没有说服 力了。常见的测试模型如s h e p p - l o g a n 标准头部模型，由许多大小密度不等的椭圆组合而成。 在算法性能测试中，最好采用几种不同的测试模型，避免样本的个体因素影响。 在客观评价指标与主观视觉判别上有时也存在矛盾，因为不同的图像适用于不同的任务， 专科医生看一副图的视觉感受和注意点，与学术研究人员是不可能完全相同的，比如将高频噪 声当作噪声滤除，可能一般人看起来会觉得图像视觉效果变好了，没有那些干扰点了，但是医 生有可能专注的恰恰是某一个边缘、细节或质感，他恰恰能够撇开干扰点把握住他想要的东硝。 而往往在客观指标与主观评价之间往往存在权衡。当然这些矛盾是在一定的可接受范围内的， 大部分时候是以人眼视觉的主观评价为重，图像的基本清晰度是大家都能判断的。在更精确的 评价中，将更多的考虑目标图像的特点和图像的用途，有针对性地进行更具体的评价。 算法的可行性评价主要是两方面，一是对数据的操作是否简单易行，二是计算量火小，也 可以说计算速度。评价和比较计算速度要在同样的计算机配置下进行，而比较基本运算项的计 算量则更清晰。不过很多算法并不能完全分解成基本运算的组合。一般地，如果算法的主要部 分能够分解成简单算术运算，它应该是简单易行，可以朋硬件实现的了，那它的速度也可以随 着硬件速度的提升而不断加快。另外对于迭代算法，“泛的评价指标是收敛速度，迭代多少次 收敛。可以根据重建结果米衡量收敛速度，建立一些误差评价指标，根据误差的卜降速度，再 3 东南人学硕i j 学位论文 综合主观视觉，进行判断。 1 2c t 成像的基本原理 1 2 1 沿直线衰减 c t 成像的物理基础是x 射线的穿透性，当一束x 射线通过人体或物件时，由丁光子的吸 收与散射等效应而衰减，射线被吸收的光子总数取决于它穿过的物质的密度、原子组成及x 射 线的光能频谱。在x 射线能量相同的情况下，密度大的物质衰减更大，在x 照片上显示为更浅 的图像。测量不同组织或物质的衰减系数，就能从反映密度的图像上把它们区分开。 假设一束入射强度为如的单能量x 射线穿过一种均匀介质，则衰减之后接收到的射线强度 为i = o e 一肛，其中p 为该均匀介质对射线的线性衰减系数，z 为射线所穿越的路径长度。对前 式作对数变换可得：p = 昙f n ( 争) 。但通常x 射线路径上的物体密度是不均匀的。 图1 1 假想x 射线传播示意图 假设将物体分隔为很小的段，每段长度为w ，w 足够小以致每段中的密度可设想为均匀( 见 图2 1 ) ，进一步假定小段内的衰减系数相同，而长度为w 的第一段的射线入射强度为而，从第 一段穿出的射线强度为 ，从第二段穿出的射线强度为厶最后一段出来的透射强度为厶，则 = o e p l ( 1 2 ) 如= 1 1 e p 2 w ( 1 3 ) 将式( 2 1 ) 代入( 2 2 ) 得： 如= 如e 一( p l + p 2 ) w ( 1 4 ) 继续这个过程可得最后一段出来的透射强度厶： 厶= e 一( p l + p 2 + p n ) w ( 1 5 ) 上式表明：在x 射线束的一条路径上，如果已知w ， 而及厶，就能计算出物体的衰减系数总 和( p 1 + p 2 + p 。) 。但要分别求得各段的衰减系数，需要多个不同方向路径上的数据，即要获 得与未知数个数相当的方程数。二维的情况类似，见图1 2 。 当w 趋向无穷小时，x 射线在非均匀介质中的传输可以积分形式描述为： i ：o e - - 正p ( z ) a x ( 1 6 ) 其中z 为射线穿过物体的路径，p ( z ) 为路径上点z 处的线性衰减系数。把沿线信息的累加结果， 即线积分。，：p ( z ) 如定义为投影。由投影反求衰减系数函数p ( z ) ，称为反投影。如果我们用二维 密度函数f ( x ，j ) 来描述目标切面( 即断层) ，f ( x ，y ) 的值表示切面上点 ，剪) 的衰减系数，那么 c t 成像的问题就可描述为根据多个方向的投影求得各点的衰减系数，以产生切面的二维密度数 据。 相应的理论工作最早由数学家r a d o n 解决，获得c t 投影的过程可以用r a d o n 变换来描述， 而由投影重建c t 图像的算法本质上都可以归结为r a d o n 反变换在离散域的实现【。 把所有的切片投影数据排列成二维形式，称为正弦图。从函数f ( x ，耖) 到正弦图p ( r ，秒) 的映射 称为r a d o n 变换， 如图1 3 所示，二维函数f ( x ，) 的r a d o n 变换定义为： 4 第一章绪论 q ( r ，0 ) = f ( x ，v ) d z ( 1 7 ) 其中平行的x 射线束与z 轴所成角度为0 ，r 为射线束的法线向量，又： 圈= 瞄c o s0 - 咖s i n p 0 圈 8 ， 由式( 2 6 ) 可得沿s 方向的投影： g ( r ，0 ) = 虑( r c o s 0 一ss i n 0 ，r s i n 0 + s c o s o ) d s ( 1 9 ) 相应的r a d o n 反变换为： ( x ，y ) = 嘉盯d o 盅占掣出 ( 1 1 0 ) 上式中的变换包括了对投影量求差分、h i l b e r t 变换和反投影三个步骤，直接离散化用计算机实 现是很困难的，而且由于有对投影的差分步骤，对投影误差特别敏感，所以一般不直接采用r a d o n 反变换进行图像重建，而是对上述几式进行解析变换后，得剑各种不同形式的逼近表达式。 毛一， if i ?p 1 3 一i r ，一，且甜f 竹却j p 扣 k fp j n一叠 p 扣 n k 一 1 ， - 一五：厶竹9 i _ 一从+ 一j + h 五：k “p i 【h 助 _ - i i - 一厶= , u p l i 肿 i m 号 - y 、 ( x y 焐 淋。f 、心 堞龄、- 图1 2 二维平面衰减示意图 图1 3x 线二维扫描示意图 1 2 2 扫描与成像 x 传统的二维c t 获取投影的扫描方式有两种，平行束和扇形束，最早的c t 机用平行束射线 采样，后来主要用扇形束采样方式，如图1 4 和1 5 。这两种投影数据可以通过重新排列和插值 进行转换。螺旋c t 的数据采集之后则要经过层间插值以获得二维投影。 在实际应用中，c t 探测器接收透过该层面的x 线，转变为可见光后，由光电转换变为电信 号，再经模拟数字转换器( a n a l o g d i g i t a lc o n v e r t e r ) 转为数字，输入计算机处理得到多个投影值。 由投影实现r a d o n 反变换之后得剑图像数据。图像形成就是对二维( 行和列) 排列的方格即数 字矩阵做处理，实际上是一幅纵横二维排列的单位方格，经数模转换器把数字矩阵中的每个数 字转为由黑到白不等灰度的小方块，即像素( p i x e l ) 。像素是构成c t 图像最小的单位。 5 东南人学顺l ：学位论文 o b j e c t s o u r o e ( b ) 图1 4 平行束扫描示意图 图1 5 扇形束扫描示意图 1 2 3 傅里叶切片定理 傅里叶中心切片定理把投影的一维傅里叶变换与图像的二维傅里叶变换的一个切片联系起 来。即：图像，( z ，可) 在角度日得到的投影的一维傅里叶变换等于的二维傅里叶变换f ( 豇，t ，) 的一个 切片，该切片与牡轴夹角为口，如图1 6 所示。 图1 6 投影的一维傅里叶变换和图像二维傅里叶交换的关系 如果我们在不同的角度下取得足够多的投影函数数据，并作它们的傅立叶变换，那么，变 换后的数据就将充满整个( u ，u ) 平面。一旦频域函数f ( u ，口) 或f ( p ，p ) 的全部值都得到后，将其作 一次傅立叶反变换，就能得到原始的密度函数f ( z ，剪) ，也就是我们所要重建的图像。 此方法即为如图1 7 所示的傅立叶变换重建方法，它也是广泛应用的f b p 重建方法的推导 基础。 由图可知，对丁二每次测得的投影数据先作一维傅立叶变换，根据中心切片定理，可将此变 换结果看成二维频率域中同样角度下过原点的直线上的值。在不同投影角下所得的一维变换函 数可在频域中构成完整的二维傅立叶变换函数，将此二维变换函数作一次逆变换就得到了所要 求的空间域中的密度函数。为了在二维逆变换中采朋快速傅立叶变换算法，通常在逆变换前将 极坐标形式的频域函数变换成直角坐标形式的数据。 6 l 1 蝴 第一章绪论 图1 7 傅里叶重建法示意图 要能够实现直接傅里叶重建，笛卡尔坐标与极坐标之间的插值转换精度是个关键问题。由 不精确插值导致的伪影限制了该方法的应用。人们曾考虑过很多插值方法，最近邻点、双线性 和截断s i n c 函数f i r 插值器，理论上还可以证明，使用雅克比加权的二维周期s i n c 核插值函数 的直接傅里叶变换法等价于f b p 方法。换个角度想，这种插值器并不是最理想的，那么使用理 想插值器的傅里叶重建效果应该优于f b p 重建，应该是可以从重建速度和重建效果上改进傅里 叶重建法的。 1 3 相关研究回顾 加快迭代成像的速度有两个最直观的思路，一是减少迭代次数，这就要求加快算法的收敛， 二是减少每次迭代计算所需要的时间。 加快算法收敛的努力一直在进行，g u a n 等人通过选择特定的投影访问顺序来加快代数迭代 重建的收敛速度【1 1 】。h u d s o n 等人提出的有序子集方法，在统计迭代中已经得剑了广泛的应j j 【1 2 1 。 由于该算法本身并不收敛于最人似然解，并且可能在迭代过程中导致局部收敛，因此研究重点 转移到了收敛的有序子集方法上，其中一类方法如行替换、法【1 3 l 和b s r e m 1 4 】等，利用松弛因子 来保证收敛，不过松弛因子的取值变化规律的经验表格的选取有难度，且对结果的影响很大。 另一类是像c o s e m 1 5 】这样的利用所有子集信息的方法，具有比最大期望似然法( m l e m ) 更快的收敛速度。 迭代中反复进行的一个步骤是，从图像空间向投影空间转换的前向( 重) 投影。对于较大 的图像，这需要存取和计算相当庞大而常常是稀疏的系统矩阵，对存储空间和处理器的性能要 求都很高。如果采用实时计算，则其耗时占总计算时间的比例大于5 0 。可以考虑尽量利用成 像系统几何结构的对称性，来实现对系统矩阵的加速计算，比如引入极格模剁1 6 j 。 极端的解决办法是放弃重投影的努力，把迭代逼近步骤完全转移到投影域，这就需要对投 影数据有针对性的、较精确的建模，如w a n g 等人提出了基于新的统计模型的理想投影估计方 法，以使j f b p 方法进行快速重建【1 7 】。另一种投影估计思路则采用传统的泊松模型【1 8 j 。 l i n 等人也考虑到了在投影域和图像域转换的矩阵运算代价，尝试仅在图像域施加正则化约 束，而在投影域进行更新和估计，提高了迭代速度【1 9 1 。 更具吸引力的方法是不采用系统矩阵相关运算米完成重投影，而用1 e 均匀快速傅里叶变换 ( n u f f t ) 实现这一步骤，其中的关键是保证插值精度，采川优化插值方法米减小重建误型驯。 7 东南人学硕l ：学位论义 已有很多基t - n u f f f 思路的尝试应_ h j 丁直接傅里叶变换重建【2 m i i2 1 。 1 4 本文工作及章节安排 本课题主要关注具有抗噪性能的快速c t 成像方法，通过对算法的分析和讨论，结合计算机 仿真实验来得出结论。基本思路是把迭代法与解析法相结合，提高迭代成像的速度。 本文的i ：作分为三个部分：第一章和第二章讨论常见的解析重建算法和迭代重建算法，以 及相关的加速方法。 第三章将迭代法与滤波反投影法结合，采用基丁投影域统计迭代的快速重建方法，即在对 投影数据进行分析并建立合理的统计模型的基础上，估计理想投影数据，再将得剑的估计投影 用具有速度优势的滤波反投影法重建。 第四章将迭代法与傅里叶重建法结合，采用优化k a i s e r - b e s s e l 和m i n m a x 插值实现基于非 均匀快速傅里叶变换的重投影，改变传统迭代方法中耗时严重的重投影方式，以提高迭代速度。 本文的章节安排如下： 第一章：绪论部分。 阐述课题的研究背景和社会意义，c t 基本原理，研究现状。说明论文结构。 第二章：主要的c t 成像重建算法及其加速。 ( 1 ) 滤波反投影算法及其快速算法。 ( 2 ) 迭代重建算法中系统矩阵的快速计算方法。 第三章：基于投影域迭代估计的快速成像。 ( 1 ) 介绍不同于经典泊松模型的高斯统计模型，根据该模型来建立目标函数，寻找合适的 优化方法进行投影估计。 ( 2 ) 对基于投影估计的快速重建思路进行数值实验和结果分析。 第四章：基于非均匀快速傅里叶变换的快速重投影。 ( 1 ) 讨论非均匀快速傅里叶变换在直接傅里叶重建中的应用。 ( 2 ) 采用该思路进行投影算子的数值实验和结果分析。 第五章：总结。 ( 1 ) 本文有哪些内容。 ( 2 ) 需要进一步探讨的内容。 8 第- 二章c r 成像系建方法 第二章c t 成像重建方法 2 1 滤波反投影重建及其加速 2 1 1 滤波反投影( f b p ) 重建法 滤波反投影重建算法的核心问题是：根据数据的采集方法、被重建物体的类型以及重建时 所选择的标准，决定所需要的投影数目，选择适当的滤波和插值函数。 为了重建物体内部的衰减系数分布，f b p 算法把投影值按射线穿越路径反向投射过去，把 该投影值赋予该路线上所有像素，将各个角度的反投影替加，最后重建出断面图像。直接反投 影会使图像发生模糊，为了解决这个问题，通常把投影数据与加权函数卷积之后再投影。 该算法的不足之处在于：第一，它是基于解析公式反求的闭合形式，因而要求投影数据是 精确的、完全的，也就是说，该算法对噪声是敏感的。第二，其计算复杂程度强烈依赖丁- 数据 采集的扫描方式，如果投影数据不是沿直线的简单积分，就得不到该闭合解，算法失效。 标准的f b p 算法已经相当成熟，其推导过程和实现步骤可参考经典文献【2 3 】和【川，下图为算 法流程示意： 图2 1 滤波反投影重建算法框图 算法中滤波器的理想频率响应为：日( j d ) = 洲。这是一个频带无限的滤波器，实际应用中， 通过对理想的频率响应函数加窗来实现。即，滤波器的频域表示为h ( p ) = pw c p ) ，式q 3 w ( p ) 为窗函数，使用不同的窗函数，就可以得到不同的滤波器。 2 1 2f b p 算法的加速 标准f b p 算法对于大型数据量的重建来说，其重建速度还是有待提高。s a h i n e r 等人1 2 5 j 与 h e 等人【2 6 】都在1 9 9 3 年基于像素与射线在不同投影角度上的儿何关系，提出了多方位同时反投 影的快速f b p 算法，将反投影和线性插值操作中的乘法次数减少了一半。他们采用四点对称的 计算思路，每次同时计算4 个投影角：0 ，7 r 2 0 ，7 r 2 + 0 ，7 r p ，而用来做投影插值的像素定 位操作只需要计算一次( 要求投影数为偶数) 。 y j h e 等人推导了h a r t l e y 变换与r a d o n 变换之间的关系，提出了类似于傅里叶中心切片定 理的新的投影定理，在此基础上可以导出基于h a r t l e y 变换的f 1 3 p 算法。由于不涉及复数运算， 在计算餐和存储空间方面有优势唧。 基于逐步求精思路的快速方法如冈数分解反投影法，在不同阶段依次提高整幅图像的分辨 率，直到达到最后的最佳分辨率例。w e iy u c h u a n 等人于2 0 0 5 年推导了反投影与滤波反投影之 间的关系，分别对扇形束情况下的滤波反投影罅法与直接反投影算法的连续形式进行级数展开， 发现了一个有趣的关系，如果把滤波反投影式用无穷级数来表示，推出它的n 阶近似公式，直 接反投影则是它的一阶近似【2 9 1 。这一关系式也可看作逐步求精的快速重建的基础，将可以通过 q 东南人学硕i ：学位论文 不同阶数近似来权衡重建时间与精度。 b a s u 等人于2 0 0 0 年提出了平行束快速分割反投影( f a s th i e r a r c h i c a lb a c k p r o i e c t i o n f h b p ) 重建算法，其基本思想是将一幅图分割成四幅子图。假设现有p 个滤波投影可被川于重建整幅 图像，为了重建其中一个子图像妒，p 个投影首先被截取为该子图像的投影，然后移到以对应 该它的中心，最终投影数被减为尸2 个。然后，使用这p 2 个投影重建出子图像妒；完成四幅子 图像的重建，加以整合便可重构出完整的图像，依次递! j ： 执行上述分割过程，最终总重建计算 量就会减少至：w j o ( n 2l o gn ) i 驯。 该算法中有儿个特定的参数来控制其性能，包括一个径向插值核，投影角上的防混叠滤波 器，还有一个调节参数。这些参数对算法性能的影响，2 0 0 1 年他们进一步进行了讨论，采用相 应的简化算法做了统计误差分析【3 l l 。不同于【2 8 】的是，在该算法的每个阶段，每一个子图像都具 有最终图像的分辨率。 。 g e o r g e 等人把b a s u 递归分解图像的算法称为空域抽取法，他们提出了具有类似优点的快速 算法，按投影角递归分解投影数据，称为投影抽取法1 3 2 j ，对丁从尸个投影角重建的图像， 这类算法的计算复杂度奠j o ( n 2l o gp ) 。通过投影的旋转获得不同的子图像实现分级反投影，子 图像则由更少的投影重建出来的其他i 割像相加得到。之后他们将算法推广到了扇形束情况【3 3 1 。 其思路如- 卜图所示： 点 掣 。 点n 3 lo 。 ( 1 ) 1 n i ( a ) 燮绘 锄萨 io 。 i t 1 ( c ) 芝 。 0 。 ( | ) r t l ( b ) 、 ，一一 ，一，、- - ( | ) 1 7 1 ： l ( d ) 图2 2 投影抽取算法频谱示意图 ( a ) 为单个投影在9 = o 的反投影，( b ) 为旋转后的单个投影的反投影，( c ) 为从所有旋转后的 单个投影得到的图像，( d ) 为该分级算法中典型的过渡图像 1 0 第二章c t 成像重建方法 2 2 迭代重建算法 2 2 1 离散化与迭代改进 成像系统中的已知条什即探测器获得的投影数是有限的，因此要将图像近似用有限个参数 来表示。另外现有的数字化图像显示方式也需要对连续图像进行离散化采样。我们可以利h j 级 数展开法将连续图像近似表示为z = ( x l ，x 2 ，x n p ) 与基函数 b j ( r - 3 ) 的加权和的形式，其中 j = 1 ，2 ，叩为索引值： ，( 力( 刁 ( 2 1 ) 将上式代入正向投影的表达式，得： - s ic 刁，c 刁d 尹= c 力 霎巧如c 而 = n p s i c 而幻c 力d 习巧 。2 2 ， 2 乙o i j x j = 【a x i 其中。巧= ，s i ( 而幻( 力d 尹表示第歹个基函数对第i 个探测器获得的投影值的贡献；s i ( 力为反 映不同系统投影方式的函数。可以看到，n 讨和都取决于基函数 ( 开 。我们定义矩阵a = a i j 为系统矩阵，其中i 和歹分别为投影和图像基函数的索引。它表示的是图像域和投影域之间的映 本文采用像素基函数，即像素方格内部均匀地为1 ，外部为0 ，在边界