（计算数学专业论文）基于一种新的扫描方式的滤波反投影算法.pdf

中文摘要 扇形束c t ( e o m p u t e dt o m o g r a p h y ) 扫描结构在实践中易于实现和控制，被广泛 应用于临床。尽管新的快速三维螺旋c t 体积成像方法将来是c t 重建的研究主流， 而扇形束c t 也仍然有研究价值。为此本文定出了将扇形束扫描方式下的重建作为 课题的研究内容。 本文主要介绍探测器等距离排列的扇形束投影，在全局重建与局部重建的数值 模拟中，本文采用结合一种新型滤波函数 1 ，2 】的扇形束滤波反投影( f b p ) 【3 】，通过数 字验证了在全局重建中采用该滤波函数比现有常用的滤波函数效果更好，以及该 特殊函数在局部重建中的有效性。 针对工业c t 中存在感兴趣区域( r o t ) 的旋转半径有时过大，当感兴趣区域 ( r o i ) 绕旋转中心旋转时会落在扫描范围外的问题，本文在张慧韬、张朋 4 1 i 作 的基础上，推导出了在新的扫描方式下的扇形束滤波反投影( f b p ) 的图像重建公 式，并将扇形束重排法用到新的成像模式中，构造了新的成像模式下的局部重建 算法。在数值试验中实现了用局部投影数据进行重建，用原物体的数据与重建后 的数据对比，发现滤波反投影算法在新的扫描方式下的局部重建效果较好，并且容 易实现。 关键词：图像重建；r a d o n 变换；滤波反投影算法；扇形束局部重建；重排法 j 匕塞变通太堂亟堂焦论塞旦s 至! a bs t r a c t f a n - b e a mc t ( c o m p u t e dt o m o g r a p h y ) ，w h i c hi se a s yt ob ei m p l e m e n t e da n db e c o n t r o l l e di np r a c t i c e i sw i d e l yu s e di nc l i n i c a lm e d i c i n e 删l et h ef a s t3 ds p i r a l r e c o n s t r u c t i n gi m a g e sm e t h o dw i l lb et h em a i n s t r e a mi nf e a t u r er e s e a r c h ，i ti s s t i l l v a l u a b l et or e s e a r c ht h ef a n - b c a mc t s ot h er e c o n s t r u c t i n ga l g o r i t h mi n v o l v i n gt h e f a n b e a ms c a n n i n gi st h er e s e a r c hs u b j e c ti nt h i sa r t i c l e w ep r o p o s et h ef a n - b e a mc ti n v o l v i n gc o l l i n e a re q u i s p a c e dd e t e c t o r i nt h eg l o b a l a n dl o c a lr e c o n s t r u c t i o n ，w ep r o p o s et h e2 一df a n b e a mf i l t e r e db a c kp r o j e c t i o n a l g o r i t h m ( f b p ) 3 】w i man e w w i n d o wf u n c t i o n 1 ，2 】i ns i m u l a t i o ne x p e r i e n c e ，b y c o m p a r i n gt h eg l o b a la n dl o c a lr e c o n s t r u c t i o ni m a g e sb yn e ww i n d o wf u n c t i o na n d g e n e r a lw i n d o wf u n c t i o n ，t h en e ww i n d o wf u n c t i o ni se s t i m a t e d s o m e l ：i m e st h er o t a t i o nr a d i u so ft h er e g i o no fi n t e r e s ti st o ol a r g et ob ei n s i d et h e s c a n n i n gr a n g ei ni n d u s t r i a lc t o nt h eb a s eo f t h eh u i t a oz h a n g 、p e n g z h a n g sw o r k , t h ei m a g er e c o n s t r u c t i o nf o r m u l ai n v o l v i n gt h en e ws c a nm o d ei sg i v e ni nt h i sp a p e r a n e wl o c a lr e c o n s t r u c t i o na l g o r i t h mi sc o n s t r u c t i n gu s i n gt h er e c o m p o s i t i o nm e t h o di n t h en e wr e c o n s t r u c t i n gm o d e d u r i n gt h es i m u l a t i o ne x p e r i e n c e ，t h er e c o n s t r u c t i n g u s i n gt h el o c a lp r o j e c t i o nd a t ai si m p l e m e n t e d b yc o m p a r i n gt h ed a t ao fo r i g i n a li m a g e a n dr e c o n s t r u c t i o ni m a g e ，i ti se a s yt of i n dt h a tt h ef i l t e r e db a c kp r o j e c t i o na l g o r i t h mi s e f f e c t i v ea n de a s yt ob ei m p l e m e n t e di nt h en e ws c a n n i n gm o d e k e y w o r d s ：i m a g er e c o n s t r u c t i o n ；r a d o nt r a n s f o r m ；f i l t e r e db a c k p r o j e c t i o n a l g o r i t h m ；f a n b e a ml o c a lr e c o n s t r u c t i o n ；r e c o m p o s i t i o nm e t h o d 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的研 究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表或 撰写过的研究成果，也不包含为获得北京交通大学或其他教育机构的学位或证书 而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作 了明确的说明并表示了谢意。 学位论文作者签名：啊匆 3 5 7 月g - 日 f 学位论文版权使用授权书 本学位论文作者完全了解北京交通大学有关保留、使用学位论文的规定。特 授权北京交通大学可以将学位论文的全部或部分内容编入有关数据库进行检索， 并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校向国 家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 靴敝储戤：枋栊为 签字日期：砂7 年7 腿日 夕 月 匿 d 轭 年 臻刁 名 ： 签 期 爵：- j 师 日 导 字签 致谢 本论文的工作是在我的导师渠刚荣副教授的悉心指导下完成的，渠刚荣副教 授严谨的治学态度和科学的工作方法给了我极大的帮助和影响。在此衷心感谢两 年来渠刚荣老师对我的关心和指导。 渠刚荣副教授悉心指导我们完成了实验室的科研工作，在学习上和生活上都 给予了我很大的关心和帮助，在此向渠刚荣老师表示衷心的谢意。 渠刚荣副教授对于我的科研工作和论文都提出了许多的宝贵意见，在此表示 衷心的感谢。 在实验室工作及撰写论文期间，顾兴勇、肖波等同学对我论文中的图像重 建算法研究工作给予了热情帮助，在此向他们表达我的感激之情。 另外也感谢家人朋友，他们的理解和支持使我能够在学校专心完成我的学 匕。 1 绪论 1 1图像重建问题的发展及展望 计算机断层扫描( c o m p u t e r i z e dt o m o g r a p h y 简写成c t ) 在医学上一个重要 的方式是由多个x 射线投影值获得人体内部的密度分布。1 9 1 7 年奥地利数学家 r a d o n 系统地论证了由积分值确定被积函数的理论方法，为c t 技术的形成和发展 提供了可靠的理论依据。随后，美国物理学家a m c o r m a c k 在此基础上作了进一 步研究，并且完成了仿真与实验研究，这标志着由x 射线投影重建图像的解析数 学方法的确立。1 9 7 1 年，英国e m i 公司工程师g n h o u n s f i e l d 研制出了第一台c t 扫描装置。至此，c t 理论终于取得了实质性的突破。 在二维的情况下，r a d o n 变换是非局部的，即为了恢复图像中的任何一个特定 的点需要待重建物体所有的投影信息。在医学上，即使医生只想了解病人的某个 局部的情况，也要对病人的全身进行扫描，这意味着病人要承受更多的x 射线辐 射，过多的电离辐射可能会引起病人的癌变或者基因变异，对病人的身体带来危 害。而在实际情况中，医生更多时候只是想了解病人的某个特定器官的情况，如： 心脏，大脑等等。为了减少x 射线对病人带来的伤害，同时也为了减少扫描的时 间，提高资源利用率和计算的时间，针对感兴趣区域( r o d 局部重建迫切需要发 展。 由于局部重建算法( l o c a lr e c o n s t r u c t i o na l g o r i t h m ) 的实用性，局部重建算 法在近几十年来是国内外比较关注的话题。1 9 8 4 年s m i t h 首先提出了l a m b d a t o m o g r a p h y ( l t ) 5 6 】。他把c a l d e r o n 在解偏微分方程时使用人引入到c t 重建 算法中，利用算子的数学性质实现了局部重建。l a m b d a t o m o g r a p h y 是最早提出的 局部重建算法，有很多文献对该算法进行了推倒和更进一步地研究。1 9 9 2 年w a l n u t 第一次把小波分析应用到了c t 局部重建中 7 8 。2 0 0 7 年宋沛然在范惠荣、徐茂 林、邱钧工作的基础上，结合重排法，实现了扇形束的滤波反投影全局与局部重 建。关于三维c t 重建算法，虽然早在上世纪8 0 年代就给出了精确重建的充分必 要条件，但在数值上很难实现。2 0 0 2 年k e t s e v i c h 提出螺旋锥束c t 精确重建算法 【9 ，该算法是三维c t 领域的突破性进展。基于他的工作，z o u 和p a n 对投影数据 进行微分、反投影，然后对反投影的结果进行h i l b e r 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 0 1 6 】。该算法在减少数据冗余和区域 ( r e g i o no fi n t e r e s t , r o i ) 重建方面有很大的优势。随后e n o o 等人用经典有限 h i l b e r t 反演公式推倒出了局部重建方法 1 7 】。 c t 设备的不断发展，由原来的第一代的二维扫描发展到了今天的螺旋c t ， 同时，对锥束重建算法的研究也随着活跃起来。虽然理论上、算法上和应用上正 在不断地完善，但是这三个方面还有很多问题等待我们去解决。特别是研究更高 效得长物体、短物体精确三维锥束重建算法 18 1 。与此同时，尽管新的快速三维螺 旋c t 体积成像方法将是未来c t 重建的主流，而扇形束c t 也仍然有其研究价值。 近来，主要研究两方面的问题：一方面是感兴趣区域超短扫描精确重建方法和理 论分析证明，另一方面是如何处理带有部分截断投影数据的精确重建问题 1 9 】。 1 2 本文的研究工作 本文基于r a d o n 变换理论，以滤波反投影为研究对象，对c t 的扇形束扫描模式 和图像重建的算法作了进一步的研究。 论文各章内容安排如下： 第一章介绍c t 局部重建的研究背景、发展历史和未来c t 局部重建的发展 趋势。 第二章介绍c t 图像重建的原理和相关的数学理论基础，本章主要针对解析 重建，为此重点介绍r a d o n 变换、r a d o n 变换的求逆公式、f o u r i e r 分析和f o u r i e r 切片定理。 第三章本文主要针对与探测器等距离排列的扇形束的图像重建算法，通过平 行束的投影数据与扇形束投影数据之间的关系，通过平行束滤波反投影的重建公 式，推导出了探测器等距离排列的扇形束图像重建算法。在计算机实现部分中， 本文采用数据重排的方法，即通过线性插值把扇形束投影数据转化为平行束投影 数据 2 0 2 2 ，再进行平行束滤波反投影重建。在滤波函数的选择上，采用范惠荣、 徐茂林、邱钧等研究的一种新型滤波函数，将其用于全局和局部重建。数值模拟 部分把滤波函数与其它常用滤波函数做了对比，证明其有效性。 第四章首先介绍了新型的c t 旋转台边旋转边平移的扇形束扫描模式，其次 在扇形束滤波反投影重建算法的基础上，推导出了新的扫描模式下的图像重建算 法，最后将采用新型滤波函数的扇形束重排法应用到新的成像模式下，通过数字 实验，将模拟数据与待重建图像原始数据进行比较，说明在新的扫描模式下，扇 形束滤波反投影重建算法的有效性。 2 2c t 图像重建的原理及相关的数学理论基础 在这一章，我们首先对二维图像重建问题进行简单介绍，主要介绍解析成像 重建，给出了解析成像重建的重要的数学理论基础一r a d o n 变换的概念、r a d o n 变 换的求逆公式以及f o u r i e r 切片理论，为下面的研究做准备。 2 1图像重建的原理 对于给定某种公式来估计密度厂的近似值厂，我们称这种方法为解析重建法。 目前，对应不同的x 光的投影结构，投影数据与密度厂之间的解析关系是不同的， 数学计算也不同。在二维的x 光的投影结构中，主要归纳为两种：平行束投影和 扇形束投影。但是，由于扇形束投影与平行束可以相互转化，在本节中，将主要 介绍平行束投影。 图2 1 平行束扫描模式 平行束扫描的机器包括射线源和用于接受投射光的探测器。在每一组数据采 集中，射线源和探测器之间作平行移动，探测器记录该组的投射数据；然后射线 源和探测器同时围绕被探测物体的中心旋转一个固定角度，接着再进行平行移动 扫描，直至旋转角度到1 8 0 0 。再根据这些投影数据进行重建。 当x 射线穿过被探测物体( 比如人体，地球物理表面等) ，由于吸收和微射， x 射线会发生衰减，又由于被探测物体内部组织的不同，衰减系数随之不同，我 们用函数( 五y ) 来表示衰减系数。由于受到密度函数厂的影响，只要知道了的 分布，也就知道了密度分布。 3 设x 射束穿过某个不均匀介质时，射束强度下降可用下面方程表示： i = i o e 一胆 厶是输入射线的强度，i 是测得的数据强度，经过变形上式可以简写为： 互鹏y ) 凼一l n ( 丢) 1 2 - l _ 1 ) 三是x 射线经过的一条路径。c t 图像重建通过探测厶和，得到一系列的投影 数据，利用( 2 1 1 ) 式，确定的分布函数。 2 2r a d o n 变换 本节先给出甩维空间中的r a d o n 变换的定义 2 3 】；再介绍二维图像重建中的 r a d o n 变换的具体形式。 2 2 1r a d o n 变换的定义 r a d o n 变换起源于积分几何，所谓的积分几何就是根据流形的几何性质( 包 括几何形态和集合分布) 研究函数在此流行上的积分。在c t 中的含义是获得投影 数据，即根据流形( 即x 射线束) 的几何性质( 各种平面直线) 研究函数( x 射 线穿过物体的衰减系数分布函数4 x ，y ) ) 在此流形各直线上的积分。 考虑，z 维空间的超平面： l ：p = 石国 其中x = ( _ ，x 2 ，x n ) r ”，s ”1 ，缈是刀维空间中的单位向量，p r 1 是 原点到超平面的距离，p ，缈为超平面的位置参数。 定义 2 2 1 r a d o n 变换是函数f ( x ) 在超平面上的积分， f ( x ) = ( ，x 2 ，) 矽( r ”) r f ( p ，国) = 少 ) d x m x = p 4 r a d o n 已经证明，若函数厂( x ) 连续且具有紧支集，则r f ( p ，缈) 可由沿所有线 积分为一确定。 显然，r f ( p ，秒) 是单位柱面z = s ”1xr 1 上的偶函数，即 r f ( - p ，- c o ) = r f ( p ，c o ) 2 2 2 二维图像重建中r a d o n 变换的具体形式 在二维的图像空间中，设厂 ，y ) 表示某一图像密度的二维空间分布，且该函 数在二维空间中具有紧支集，r a d o n 变换可以表现为： r f ( p ，p ) = i ( x , y ) d s ( 2 - 2 - 1 ) p 为原点到直线的垂线，乡为法线与z 轴的交角，p ，口为直线l 的位置参数。 它实际上是对图像厂o ，y ) 的线积分。 设如图2 - 2 建立直角坐标系，则直线l 上的任意一点y ) 可以表示为： 式子( 2 2 - 1 ) 可以改写成： r f ( p ，9 ) = 王厂( z ，y ) d s = 州p c 。s 秒一s s i n 秒) ，( p s i n p + s c o s o ) ) d s 图2 2 二维r a d o n 变换 5 p 秒 m 宝 s c s s 一 矽秒 s l出 p p = = x y ，f【 l 2 3r a d o n 变换的求逆公式 在上一节中给出了图像重建中的r a d o n 变换的具体形式，在本节中将给出其 相应的求逆公式 2 3 】： f ( x , y ) - - 萨1 胁南半勿 ( 2 - 3 - 1 ) 由式子( 2 3 1 ) 可以看出，已知投影数据r f ( p ，秒) ，要求其逆变换，可以经过 以下几个步骤： ( 1 ) 关于影( p ，p ) 的第一个变量p 求偏导数； ( 2 ) 对其偏导数做关于p 的h i l b e r t 变换； ( 3 ) 对经过h i l b e n 变换后的函数作反投影和归一化运算。 在反投影运算前经过微分、h i l b e r t 变换两个环节，由于存在微分项，公式对 实际投影数据中的微小误差极为敏感 2 4 2 6 。并且h i l b e r t 算子又具有奇异点，在 实际应用中带来许多困难，因此r a d o n 求逆公式目前在商用c t 中未得到应用。 2 4f o u r i e r 分析和f o u r i e r 切片定理 1 8 2 2 年法国数学家傅立叶发表了“热的分析理论”著作，提出并证明了周期 函数都可以展开成正弦技术，奠定了傅立叶级数的理论。其后，该理论被应用到 了处理信号中，傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域 信号( 信号的频谱) ，在频域内对信号进行处理、加工。最后利用傅立叶逆变换将 频域信号转换成时域信号。在计算机图像处理中，也用到了傅立叶变换。 2 4 1 二维f o u r i e r 变换 设f ( x ，y ) 是两个独立变量x 、y 函数，那么它的二维傅立叶变换可以表示为： i u f ( u ，v ) = f ：f ：f ( x , y ) e t 酬卅 ，) d x d y，v ) = il 州卅 p ， ( 2 4 1 ) 对f ( u ，力进行傅立叶逆变换，即把公式( 2 4 1 ) 两边分别乘以 e x p ( j 2 7 r ( u a + 印) ) 后再对变量“1 ，积分，可以得到： 6 f ( ) e j 2 r ( u a + v ，8 ) d u d v = 厂( z ，y ) e - 2 q ( u ( x - a ) + v ( y - p ) ) d x d y d u d v = 厂( 工，y ) 出方艮2 咖( p 咖咖呐) d u d v ( 2 - 4 2 ) 为了更好地说明下面的工作内容，在本节中先引入狄拉克函数的概念： 定义2 4 1 函数8 ( x , y ) 【2 7 若满足： j m ，y ) d x d y = 1 【万( 砖j ，) = o ，z 0 ，y 0 则称函数8 ( x ，y ) 是狄拉克函数。 狄拉克函数8 ( x ，y ) 具有如下性质： 性质1 ：厂( z ，y ) 8 ( x , y ) d x d y = f ( o ，o ) 性质2 ：厂( 石，y ) 艿( x 一而，y - y o ) d x d y = f ( x o ，y o ) 性质3 ：8 ( - x ，- y ) = 8 ( x ，y ) 根据狄拉克函数的定义知函数；0 。；, o e - 2 n j t u x + v y ) 幽咖是一个狄拉克函数 1 9 】，利 用狄拉克函数的性质可以把( 2 4 2 ) 式进一步改写成： f ( 叩) “甜棚d u d v = m ，y ) d x d y ；, 。艮。2 酬吣训w o 邛 d u d v = ；f ( x , y ) 8 ( x - a ，y - f 1 ) d x d y = f ( 口，) 7 ( 2 - 4 - 3 ) 将( 2 4 3 ) 式中的口、分别用x 和少代替，则傅立叶逆变换可以写成： 似，y ) = f ( “，v ) e j 2 r ( n x + v y ) d u d v 2 g 2f o u r i e r 切片理论 ( 2 - 4 4 ) 利用前面介绍的r a d o n 变换计算图像f ( x ，y ) 在秒方向上的投影矽( p ，0 ) 另( f ) = l f ( x , y ) d s p = x c o s o + y s i n 8 ，利用狄拉克函数的性质，上式可以写成： 局( f ) = 厂( 工，y ) 8 ( x c o s o + y s i n 0 一t ) d x d y 再对易( f ) 作关于f 的一维傅立叶变换得： & ( 国) = f 。p o ( t ) e - j 2 m t d t 我们考虑傅立叶变换在，= 0 时的状况，( 2 4 1 ) 式可以简化为： f ( u ，o ) ：f ：f ：f ( x ，y ) e - j 2 ”d x d y，o ) = ii ， j p = m ，y ) d y e - j 2 撇d x ( 2 - 4 - 5 ) ( 2 4 6 ) ( 2 - 4 - 7 ) 以t 为横坐标，与直角坐标系x o y 同原点，建立直角坐标系( 如图2 3 所示) ， 当0 = 0 时，( f ，s ) 坐标系和( x ，y ) 坐标系是重合的，即t = x 、s = y ，故( 2 - 4 - 5 ) 式可 以写成： 易：o ( z ) 2 上。( x ，y ) d y o 从而可得 8 ( 2 - 4 8 ) f 似，0 ) = e p e = o ( t ) e - j 2 d t ( 2 - 4 9 ) ( 2 - 4 9 ) 式说明厂o ，y ) 在垂直于x 轴方向上一维投影的傅立叶变换等于厂j ，) 的傅立叶变换，( “，d 在，= o 处经线上的值，( “，o ) ，即f ( u ，o ) = ( 缈) 。 协一 ， 、 f ，_ 奠而一一 匕 羚 图2 3f o u r i e r 切片定理 可以把上面的结论加以推广，设f ( r o ，乡) 是，( “，) 沿秒角方向经线上的值，而 且s 。佃) 设为易( f ) 投影的一维傅立叶变换，由( 2 4 - 9 ) 式定义的投影的傅立叶变 换和二维傅立叶变换在角度秒所对应的值相等，即f ( r o ，0 ) = 岛( 国) 2 8 】。这也就是 说一个方向投影的傅立叶变换，是二维傅立叶变换的一个切线，且该切面垂直于 投影射线并通过频域中一l , - 零点 2 9 3 1 1 ；反过来，有无限个不同方向的切线叠加在 一起，可以构成一个完整的二维傅立叶变换，这也就是傅立叶切片理论【1 9 】。 9 3 二维滤波反投影重建算法 滤波反投影算法是目前c t 中普遍使用的一种重建算法。本章首先分别介绍了 平行束和扇形束的滤波反投影重建算法，然后介绍扇形束的局部重建算法，为第 四章的新的扫描模式的扇形束局部重建的研究打下基础。 3 1平行滤波反投影重建算法 从前面介绍的傅立叶切片理论可以看出，一个方向投影的一维傅立叶变换， 是二维傅立叶变换的一个切线；反过来有无限多个不同方向的切线叠加在一起， 可以构成一个近似的二维傅立叶变换。再对近似的傅立叶变换求逆变换可以得到 横截面影像。这一原理给了断层图像重建的基本方法，实际应用中还存在一些其 它的问题要解决。 3 1 1 理论推导 滤波反投影算法可以通过对傅立叶变换进行坐标变换后得到设f ( x ，y ) 是需要 重建的图像，它由傅立叶逆变换获得的表达式为： m ，j ，) = e f ( 州) e j 2 x ( u x + v y ) d u d v 把上式的直角坐标变换为极坐标，令 u = 缈c o s 臼 1 ，= c o s i n 秒 则有 d u d v = a l c l c o c l o ( 3 i i ) 式可以改写成： 0 ) e2 彻“瞄“州相) ( a d a 2 d 0，y ) 5j ：土，( 缈， 2 彻瞄芦 l o ( 3 1 1 ) ( 3 - 1 2 ) ( 3 - i - 3 ) ( 3 - i - 4 ) ( 3 1 - 5 ) 利用f ( r o ，0 + 1 8 0 。) = f ( - c o ，0 ) 的性质得 f ( x ，y ) ：f f 一：f ( o j ，p ) 怫j 2 j r m t d t o d o ，y ) 5j ：，p ) 怫 ( 3 - 1 - 6 ) 上式是简单式，令t = x c o s o + y s i n o ，如果利用投影切片定理，用0 角度投影的傅 立叶变换s o ( 国) 代替二维f o u r i e r 变换： d t o d o ，y ) 2j - 【l & ( 缈) 怫口删 将( 3 1 - 7 ) 式改写成 f ( x ，y ) = f q o ( t ) d t 其中 q o ( t ) - e & ( 缈) 彬2 删d a o = 一 d t = 一j z 删 ( 3 - 1 - 7 ) ( 3 - 1 8 ) ( 3 - 1 - 9 ) ( 3 - 1 - 1 0 ) q o ( t ) 表示一个滤波过程，h 表示理想滤波函数的频率。( 3 1 8 ) 就是利用每个 角度的滤波投影幺( f ) 进行反投影。 3 1 2 计算机实现 在实际情况中，t 的取值是有限值。设 e 0 ，p o ( t ) = 0 。用d 表示相邻两 平行线之间的距离，用表示相邻两个方位图之问的夹角( a = 州m ) ，p ( n d ，m a ) 已 知，其中一n 咒n ，0 m m 一1 。理想滤波函数缈有上限，设为形，根据采 样定理，采样间隔为1 ( 2 w ) 。 第一步，用f f t ( 快速傅立叶变换) 计算p ( n d ，m a ) 的傅立叶变换 跗万w 舢州蝥n ( n d , m a 矿2 删2 ) 第二步，乘以理想滤波函数和进行傅立叶逆变换 ( 3 - 1 - 1 1 ) 掰妒号言跗矿w 帕卜伽一，4 埘 上式说明g ( f ) 在投影函数离散点的值，可以通过计算s ( 后州，m ) 陋叫卅的 离散傅立叶逆变换得到。然而， 川是理想的滤波函数，根据佩利- 维纳准则 3 2 】， 这一理想滤波器是不可实现的。但是如果结合具体的成像过程，则不但能够实现， 而且能够达到足够的精度。为此，需要对理想滤波函数h 进行窗处理。 常见的滤波函数的窗函数有f 3 3 - 3 4 1 ，其中b = 西1 ： r a m - l a k 滤波函数的窗函数一w l ( 缈) ： d 掣= $ j 搿 ( 窑) s h e p p l o g a n 滤波函数的窗函数- - s i n c 函数： s ；n c c 詈，= 等 h a m m i n g 滤波函数的窗函数一j l ( 缈) ： h ( c o ) = 州h ) c o s ( 警) 其中口- o 5 4 相应的滤波函数只要用理想滤波函数俐乘窗函数即可。在本文中，窗函数用 h ( c o ) 表示，则将( 3 1 1 2 ) 式改为： q ( n d , m a ) = 熹言跗矿w 蝴础矽w 蒯伽( 3 - 1 - 1 3 ) 第三步计算f ( x ，y ) j l l ，一l f ( x ，y ) = q ( x c o s ( m a ) + y s i n ( m a ) ，m a ) m = o 1 2 ( 3 1 - 1 4 ) j e 立窑煎厶堂蚯堂位监塞三缍逮蘧厦控墅重建簋莹 上式中jc o s ( 胍) + ys i n ( m a ) 7 f 定和t 完全相等，可以用线性插值的方法， 利用口( n d ，m a ) 近似q ( xc o s ( m a ) + ys i n ( m ) ，埘) 。 在本章节中，对s h e p p l o g a n ( 5 1 2 x5 1 2 ) 进行全局重建，f b p 的滤波函数选 月j 了r a m - l a k 滤波函数、s h e p p - l o g a n 滤波函数和h a m m i n g 滤波函数。 结果如f ： ( c ) s h e p p 一o g a r t 滤波函数重建图像( d ) h a m m i n g 滤波荫数重建图像 图3 一l 平行柬滤波反投影重建 由圈3 - l ( b ) 、图3 - 1 ( c ) 和图3 - 1 ( d ) 可以看出，与原始模l o 相比，采用三种滤波 函数均获得了很好的重建质量，崩肉眼无法区分。下一章节中将引用评价函数对 这二种滤波函数的重建精度作进一步客观的评价。 3 2扇形束投影重建算法 上一节介绍的平行束滤波反投影重建要求x 射线是平行束，但是在实际情况 中，由于探测器和焦距的距离很小，很难的到平行光，只能获得扇形束投影。在 这一节中，将介绍扇形束的滤波反投影重建。扇形束扫描方式采集数据时，射线 源和接受器同步绕着受检物体旋转。x - 射线源对应一排足够多检测器，使所有检 测器和射线源所对应的角包含整个待重建区域。根据探测器类型的不同，扇形束 投影算法分为等距( 线阵探测器) ( 如图3 - 2 ( a ) ) ，等角( 弧形探测器) ( 如图3 - 2 ( b ) ) 两种。 l j 。紧西过_ 7 潞曩 ( a ) 等距扇形束扫描模式 一沁一， 一 ( b ) 等角扇形束扫描模式 图3 2 扇形束扫描模式 1 4 3 2 1 探测器等距离排列扇形束投影重建算法 蝴淼一娥瓤探测器等距椭懒 竺竺道，表示旋转的角度。探测器；忑乏芜嚣i 蕃笺篓置二萎竺粤竺体运动 叠嚣耋雷徽嚣爱蒹蓬揣? ：与射线船与研q 的交点彳的坐标，二二菇7 ：? 习兰竺鬯交点亏的坐标 拿! - 个投影吩 ，) 与q 皿上的投影值足。三j 磊萎：j 掌竺? 。，鼍此在研呸上 投影? 瓣篱二0 芝雾麓蛐黼 蕊汝 j 图3 4 、，和秒、f 的关系 1 5 、，和秒、t 的关系可用以下公式表示： f = ，c o s y = x ! 而雨d 矛 秒= + 7 = + 删觚五x i ( 3 - 2 - 1 ) 根据平行束的滤波反投影公式 2 3 】有 厂( ，) = re 另( f ) ( ，c o s ( 乡一矽) 一t ) d t d o ( 3 - 2 2 ) 由于reb ( f ) j i l ( ，c o s ( p 一) 一f ) d t d o = f ”易( f ) ( ，c o s ( 乡一) 一f ) d t d 8 ，故有 m ，) = 互1f 。e 删j i l ( ，c o s ( 9 一矽) 一f ) d t d o ( 3 - 2 - 3 ) 利用( 3 2 1 ) ，把秒、f 变换为关于、x 的式子，先计算捌p = i j i d z u x 。 怍i i 淼别a d a 筇, 8 i 计算得 西苏，：d ( d 2 + x 2 ) - - 矿d x 2 ( d 2 + ，2 ) 恐 ：旦可；a u a p ：o 一( d 2 + 2 ) _ u 玳2 币1 ；郴_ l 怍i 丽d 3 将( 3 2 1 ) 代入( 3 2 - 3 ) ，整理得： m 脚= 三p 仨c 考每， ( 3 - 2 - 4 ) ( 3 - 2 - 5 )邯 志番 考虑到积分函数对自变量是以2 万为周期的函数，则上式积分可以改为： 竹聊= 互1f 仨( 番) m 州+ 删a n c 争沪面d 甭x ，瓦匆出妒 ( 3 - 2 6 ) 现在再来研究函数h ，设m ( r ，) 为待重建区域中的任意一点，射线s m 与z 7 轴 的交点坐标为石 令7 = 删a n ( 丢) ， d 湖厂2 历鼋矛、d + x ，c o s ( + a r c t a n ( - 丢) 一矽) 一了荔告 = r c o s ( p 4 - y 一矽) - - x c o s 7 = r c o s ( ( f l 一) + y ) 一x c o s y ， 龇n 72 历鼋矛 = r c o s ( p 一) c o s y r s i n ( f l 一矽) s i n 7 一x c o s y = ( 厂c o s ( p 一矽) 一x ) c o s y r s i n ( f l 一) s i n 厂 ( 3 - 2 7 ) 式可以改写成 ( 3 2 7 ) ( 爿) c 唧+ ，s i n ( 卅吾c 。s 7 一心n ( 一加i n 厂( 3 - 2 - 8 ) 把c o s y = 了翥子和s i n7 - 了荔丢i 矛分别代入式( 3 2 8 ) 得到 1 7 (工：一工)南+r、s。in。(fl+-x，矽。)(x：一x) (x；-x)(d+rsin(f1-矽) 一矗瓦了 ( 3 - 2 9 ) 令u = d + ，s i n ( f l 一扔，三= 以丽，则( 3 2 9 ) 可简写为兰( 一曲。 所以( 3 2 - 6 ) 可以改写为 现在我们具体研究一t h ( u ( x ；一工) ) 。前面已知办( f ) 是频率函数h 的傅立叶逆 变换。j i l ( f ) = e h p 2 删d 国于是可以得到 办( 一u = 肼p 胁 a 拿c o u l 缈，可以把( 3 - 2 - 1 1 ) 改写为 c i u 俨啪十2 耐口弘= 等肼耐7 = 等坳， 这样( 3 2 1 0 ) 可以写为 ( 3 - 2 1 1 ) 删) 2 互1 膀1 仨即坝扣) 赤蝴( 3 - 2 - 1 2 ) 3 2 2 计算机实现 本文将采用数据重排法进行计算机实现。扇形束投影数据重排法就是把扇形 束投影数据进行重排，使其接近平行束投影数据重排法不限制投影的角度和数据 采集的角度。因此它非常实用于射线源与探测器之间距离较短的情况。 本文采用插值逼近的方法，根据扇形束的扫描数据进行内插估计，近似得出 平行束扫描的投影数据集。再应用平行束扫描数据的重建方法实现图像重建。 假设，只在等间隔的y ( 0 y v 一1 ) 处不同的角度获取投影数据，每个 相邻的角问的变化量为1 1 ，并且有 v i - = 2 万 对于每个投影，仅在长度变化量为r 的2 u + 1 个等间隔长度进行采样。并有 聊x 一 又因为密度函数厂在每一条射线上的线积分尺已知，即r 在点( 刁，订) ， 一u 2 u ，0 y v l 上的值已知，根据公式( 3 2 1 ) 知 尸( 刁丽d ，订+ a r c t a n ( 鲁) ) = r ( 肌哪 p ( 刁了手丽d ，订+ a r c t a n ( - 警) ) ，通过给定的插值方法，估计得出平行束扫描 方式中的数据p ( n d ，m a ) ，d 平行束射线间距。具体操作过程分为以下两个部分( 以 线性插值为例) ： 假设固定不变，一【，妣根据舷靠，订+ a r c t a n 铮) 的值， 对第二变量进行线性插值估计础赢，删，0 鲕州- l 1 9 令= r ，o y 矿一1 则y r + a r c t a n ( 詈) 肌( 山1 ) r + 删a i l ( 詈) 成立 则有 础宫籍，删 叭1 ) r + 舭蟹) _ 砸 。一 r + 竺竿塑鼬刁南 降2 舶， 此时得到近似的不等间隔平行束射线扫描的投影数据。 假设m 固定不变，0 m m l ，对舷印亏等，删进行线性插值，估 、+ 功2 计p ( n d , m a ) ，其中- n _ n n 钳= l 【- 刁扫n ：d d 丽j ，- u u u 贝o 7 刁了 乒丽d刀d ( 7 + ) 刁了 乒南 则有 p ( n d ，m a ) + 刀d 一刁了芝丽d 只+ 1 ) 刁万丽d ，删 删咿丽d ，删( 3 - 2 - 1 4 ) 这就是用于重建的近似的等间距平行束射线扫描的投影数据。 3 2 3 局部重建及误差分析 上一小节介绍的重建算法是传统意义上的c t 重建算法( g l o b a lc t ) ，它是一 种全局的算法，它应用于整个物体的二维的断层重建。也就是说要得到断层中某 一点的数据，需要采集该断层所有的点的投影值，需要花费很多时间和资源，在 医学里，x 射线还会对人的身体产生副作用，不利于身体的健康。在现实的情况 中，我们往往只对一部分区域感兴趣( r e g i o no f i n t e r e s t ，r o i ) ，为了减少x 射线 对人体产生的副作用，节约资源和提高计算的速度，迫切需要发展局部重建算法。 因此在最近十几年的时间里一些针对感兴趣区域( r i o ) 的局部重建算法得到了飞 速发展，这些局部算法仅需要围绕感兴趣极其邻域采集投影数据，即可重建出该 区的图像。 具体地说，对任一感兴趣的局部圆形区域d c y d ，其半径为r ，设d 7 的圆心 为d ，圆心区域b ( o ，r + f ) 是d 的同心圆，利用区域b ( 0 7 ，r + f ) 的投影数据重 建区域d 7 内的任一点的值。 图3 5 局部重建数据区域 在本章节中，在采用滤波反投影( f b p ) 基础上，将采用范惠荣、徐茂林、邱 钧等研究的一种新型滤波函数【1 】，用于局部重建。 局部重建将用到的滤波函数是： 式中，厶是归一化常数，r 为感兴趣区域内的一个点到中心的距离，a = 与， 口 d 为平行束的距离。 在全局重建中，为了让新型滤波函数与其它滤波函数相比较，在这里引进常 用的评价函数【3 6 】 ( 1 ) 归一化均方根距离测量值d ，它强调少量大误差，少数点大的偏差使得 d 较大； 叫一# ( 2 ) 归一化平均绝对距离测量值r ，它强调大量小误差： r = 黜 ( 3 ) 最坏情况距离测量值e ，表示原图像和重建结果之间的最大密度偏差： e = m a x q o 、 ， 、 ， ：i - ， 、 、， 、， ： ，7 、 ， j ， 、 - ， ，r - 一， ，-， 图4 1 扫描模式示意图 在图4 1 中，坐标系 x ，y 以d 为坐标原点，该坐标系是固定的，o 代表被 测物体旋转中心，该点可在x 轴上来回移动，以d 为圆心建立一个运动的坐标系 x ，) ，) ，表示探测物体旋转的角度( 0 2 刀) 。设射线源s 在y 正半轴上，s 到x 的距离为d ，a b 为线阵列探测器，t l b 与y 垂直，并且a b 的中点d ”在y 负 半轴上，d ，表示s 到a b 的距离，2 h 表示a b 的长度。，点表示感兴趣区域的中 心，为了确保感兴趣区域始终在探测的区域内，每次移动旋转中心，都使，在y 轴 上。风表示优的距离，则0 0 的距离可用d =