（生物医学工程专业论文）医学图像三维体积重建中插值算法研究.pdf

东南大学硕士学位论文 a b s t r a c t a so n eo ft h em o s ti m p o r t a n ts t e p si nt h r e e - d i m e n s i o n a lr e c o n s t r u c t i o n , i m a g ei n t e r p o l a t i o n c a l lo f f e rr e a l3 ds c e n ef o rc l i n i c a ld i a g n o s i sa n dt h e r a p y b ya p p l y i n gi n t e r p o l a t i o nb e t w e e ns l i c e s ， a n i s o t r o p i cd a t ac a n b ec h a l l g e di n t oi s o t r o p i cd a t a , w h i c hc o u l dh i g h l yi m p r o v e3 dv o l u m e r e c o n s t r u c t i o n sq u a l i t y t h i sp a p e ri n t r o d u c e ds e v e r a li n t e r p o l a t i o nm e t h o d s ，w h i c hw e r ef u r t h e r c l a s s i f i e da ss c e n e - b a s e di n t e r p o l a t i o na n do b j e c t - b a s e di n t e r p o l a t i o n b a s e do ns i n ei n t e r p o l a t i o n ，s c e n e - b a s e di n t e r p o l a t i o nm e t h o dw a si n s ti n t r o d u c e dw h i c h i n c l u d e sp o l y n o m i a l ，l a g r a n g c ，b s p l i n ea n dw i n d o w e ds i n ci n t e r p o l a t i o n t h en e a r e s tn e i g h b o r a n dl i n e a ri n t e r p o l a t i o na l s ob e l o n g st ot h i sc a t e g o r y n e x t ，o b j e c t - b a s e di n t e r p o l a t i o nw a si n t r o d u c e dw h i c hi n c l u d e sm a t c hi n t e r p o l a t i o n ，f r e ef o r m n o n r i g i dr i e 酉s t r a t i o nb a s e di n t e r p o l a t i o na n dl e v e ls e tm o t i o nr e g i s t r a t i o nb a s e di n t e r p o l a t i o n d u r i n gt h i ss e c t i o n ，s o m es c a t t e r e dd a t ai n t e r p o l a t i o nm e t h o d sw e r ea l s os h o w e du p t h ef n s ta n d t h i r di n t e r p o l a t i o n sw e r em o d i f i e dt oe x p e c t i n gm o r er e a s o n a b l er e s u l t s f i n a l l y , c l a s s i c a lr a yt r a c i n ga l g o r i t h m sw a su s e dt or e n d e r i n g t h e3 dv o l u m ed a t ap r o d u c e d b ya b o v em e n t i o n e di m a g ei n t e r p o l a t i o nm e t h o d s k e yw o r d s ：i n t e r p o l a t i o nb e t w e e ns l i c e s ；t h r e e d i m e n s i o nr e c o n s t r u c t i o n ；s c a t t e r e dd a t a i n t e r p o l a t i o n ；m a t c hi n t e r p o l a t i o n ；n o nr i g i dr e g i s t r a t i o n 东南大学硕士学位论文 东南大学学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。 尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人己经发表或撰写 过的研究成果，也不包含为获得东南大学或其它教育机构的学位或证书而使用过的材料。 与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢 意。 研究生签名： 东南大学学位论文使用授权声明 东南大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位论文的复印 件和电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸 质论文的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅，可以公布( 包 括刊登) 论文的全部或部分内容。论文的公布( 包括刊登) 授权东南大学研究生院办理。 研究生签名：丕叠望酗导师签名： 期： 砰j 哆 j 第一章绪论 第一章绪论 1 1 课题背景和研究意义 在现代临床医学中，作为诊断和治疗规划工具的放射医学图像已经得到广泛应用： x 线透视、x c t 、超声图像等提供了人体组织器官的形态信息，而核医学图像( s p e c t 和p e t ) 以及核磁共振图像不但提供人体组织器官的形态信息，还可以提供有关化学元 素的分布及其生物功能活动信息。医学影像在临床诊断、教学科研等方面正发挥着极其 重要的作用。传统的影像技术只是获得人体某一断层的影像数据，然后医生通过胶片进 行诊断或者通过显示屏幕进行观察。但是无论是胶片还是屏幕显示，医务人员所观察到 的仍然是一组二维断层图像，且只能以固定方式对图像进行观察，所得到的诊断结果带 有医生的主观判断，这在很大程度上取决于医生的临床经验。这种方法能够达到一般的 诊断目的，如发现肿瘤或者病变等，但对于制定手术计划等需要理解复杂甚至变形的三 维结构的临床任务，仅通过在脑中重建断层图像是非常难以完成的，而且这种方式很大 程度上依赖于医生丰富的临床经验和空间想象力。三维重建技术的应用可以改变这种状 况。通过对一组二维断层图像进行处理，从而恢复物体三维形状，建立虚拟的器官和组 织，以便使医生从多方位、多层次的观察角度对影像数据进行详细的医学分析，从而对 病变体及其它感兴趣的区域进行定性直至准确的定量分析，比如放射治疗计划中三维剂 量场的计算等。这无疑对提高影像数据的利用价值有深远的意义，而且对提高临床诊断 的准确性和正确性有很大益处，同时大大减轻了医生的负担。因此三维重建技术具有很 现实的意义。 三维重建工作包括三维图像数据的获取、处理及显示。其中三维医学图像数据的获 取一般通过对序列切片图像进行转换和插值得到。通常的医学成像系统( c t m 砌等) 都是以层的方式进行3 d 扫描的，受采集时间的限制和为减少病人受照剂量和诊断花费 等原因，一般断层图像相邻两层间的距离要大于同一层图像中相邻像素间的距离，这种 图像数据如果直接用于显示，将会在边界处形成明显的台阶状结构，对于一些细微结构， 如冠状动脉等，还会造成结构的断裂。图像插值正是解决这一问题的关键。通过对断层 图像进行某种形式的层间插值，可将非等分辨率的图像数据转变为等分辨率的数据，从 而大大提高最终三维重建结果的质量。因此，在图像的显示、处理和分析中，需要使用 插值来获得所需分辨率的图像。寻找和设计快速、有效、适用性强的插值算法，是三维 重建工作的一个重要内容，也是本课题的研究方向。 东南大学硕上学位论文 1 2 国内外研究现状 所谓图像插值就是一个图像数据再生的过程它由原始图像数据再生出具有更 高分辨率的图像数据。若将一幅较低分辨率的图像再生出一幅较高分辨率的图像，这种 插值可看作“图像内的插值”。若在若干幅图像之间再生出几幅新的图像，这种插值可看 作“图像间的插值”。重点研究的是图像间的插值方法，对于这类问题引入如下约束条件： ( 1 ) 插值图像要与原始断层图像相似； ( 2 ) 插值图像与两个原始断层图像的相似度应该分别和它与这两个断层图像的距离 成反比关系； ( 3 ) 插值图像序列应该呈现出从一幅原始断层图像到另一幅原始断层图像的渐变过 程。 插值只是改变了断层图像的空间分辨率，使三维数据的处理、分析和显示更加方便， 并没有产生新信息。图像插值方法已有大量的文献涉及，从总体上讲，现有的研究方法 大体上可分为两类：基于场景的( s c e n eb 勰c d ) 及基于对象的( o b j e c tb a s e d ) 。 这里有必要说明一点，基于对象的插值方法中有很大一类是需要进行分割预处理 的，这类算法选择一组特征点，一般是目标物体边界区域上的点，然后确定这些特征点 之间的对应关系，用这些特征对应关系确定一幅图像映射到另一幅图像所用变换函数以 引导后续的插值过程。在两幅图像中自动选出相同特征的点是一件困难的事，这一步骤 通常需要人的援助。而本文主要研究无须分割的自动化插值方法，所以对涉及分割的这 一大类方法不作讨论。 基于场景的插值方法是用已有的切片体素数据直接决定被插值体素数据 1 l 。这类方 法的研究进行得比较早，主要有最邻近插值( n e a 佗s tn e i g l l b o r ) 、线性插值( l m e 神、样条 插值( s p l i n e ) 和克里金插值( 蹦n g i n g ) 口】。最邻近插值方法是以欧氏距离最近的己知体素值 作为待插值体素值。线性插值是根据待插体素到相邻断层的距离对上下两个相邻的断层 图像进行加权平均。样条插值中使用的最多的是三次样条和改进的三次样条，它假定灰 度在三个正交方向上都按样条曲线的趋势变化，用三元三次样条张量积的形式来估计插 值点的灰度。 最早同时也是应用最广泛的插值算法是线性插值，因为它计算量小，插值速度快。 随后h c r m 衄等人3 1 将三次样条用于医学图像插值上。由于样条插值预先假定了数据按 某种样条曲线变化，当数据的实际变化规律与此有较大偏差时，则会产生不可忽略的错 误。近年来部分学者提出用统计学的克里金方法进行灰度插值，成为医学图像插值的研 究热点。本质上说，它们都是计算采样点灰度的加权平均，并且都根据插值点和采样点 第一章绪论 之间的相对距离来调整权值的大小。这些插值方法插值精度不高，产生的新断面通常会 出现边缘模糊，由此重建出的三维真实感图像表面会产生伪像，当断层间距较大时这一 点尤其明显。造成这种情况的主要原因是这些方法没有考虑到物体几何形状的变化。 无须分割的基于对象的插值方法( 以下简称基于对象的插值方法) 则从已有的切片 体素中获取对象信息来寻找被插值体素数据。大体可分为以下两步：首先寻找待插值点 在相邻断层内的对应点，然后使用基于场景的插值方法计算该点的像素值，不同的算法 可能另有不同的步骤。这类方法主要有基于匹配的插值方法、基于二进制体素的插值方 法和基于非刚性形变配准的插值方法。 最早提出的是动态弹性插值( d y n a m i ce l a s t i ci n t e r p o l a t i o n ) 。它的思想起源于b u r r 在研究图像配准时提出的弹性匹配概念，并由l e v i n 【4 】首先将该思想运用于三维重建。 9 0 年代初，对于灰度断层图像插值问题的研究主要集中在如何有效的寻找断层图像的 匹配点上，产生了一系列基于匹配的图像插值方法，它是在灰度线性插值和形状插值基 础上提出的新的改进方法。它充分考虑到生物组织在三维空间中的几何形状的不规则 性，可产生较好的中间切片图像。特征匹配的效果是插值成败的关键，边界轮廓信息是 重要的匹配信息，但不是唯一的，距离、对象结构信息、灰度数的相关信息、以及结构 走向的先验知识均可用来改善匹配效果。比较典型的方法是g o s h t a s b y 提出的基于梯度 的算法【5 】，在相邻断层的一个小窗内根据一个包含梯度和距离的评价函数找出最匹配的 点对进行插值。 9 6 年r a y a 和u d u p a 提出了一种新的基于二进制体素形状的插值方法【6 j ，这种算法 的思想来源于二维二值图像插值算法。首先将二维灰度图变为三维二值图，再计算每一 体素到器官边界的( 带符号) 距离，形成一“距离图像”，然后插值此“距离图像”到指定的 分辨率( 如用线性插值、改进的三次样条方法等) ，最后将插值后的“距离图像”转换成最 初的二进制形式，用相反的过程将目标二值图降维为二维灰度图。 9 0 年代后期，基于非刚性形变配准的插值方法逐渐兴起。这种方法的思想是相邻 切片的解剖结构在形状和灰度上的变化存在连续性，利用非刚性配准能准确模拟这种形 变，然后利用得到的形变关系用灰度插值方法得到目标体素的灰度值。这类方法使用的 非刚性形变有弹性模型( e l a s t i c ) 、流体模型( f l u i d ) 、有限元模型( f i n i t ee l e m e n tm o d e l s ， f e m ) 和基于b 样条的自由曲面形变( f r e e - f o r md e f o r m a t i o n ，f f d ) p j 。其中以p e n n e y 和 r u e c k e r t 提出的基于f f d 的插值方法哺1 为典型。 g r e v e r a 和u d u p a 等人对很多插值算法都做过很好的比较。总的说来，基于场景的 插值方法比基于对象的插值方法在计算上方便、直观，但精度较差；而后者由于它是从 几何体的外表面出发来决定未知部分，在许多场合比前者有更高的精度，因此可看作“整 体插值算法”，但此类方法在计算复杂度和所占用的内存空间一般方面都要远超前者。 东南大学硕士学位论文 1 3 论文组织结构 本文的主要工作是研究三维图像重建问题的实用化方法，在对几种常用的基于场景 的插值方法，和在基于匹配的插值方法、基于二进制体素的插值方法和基于非刚性配准 的插值方法等方法上，均进行了系统的研究。论文共分为五章，各章节安排如下： 第一章为绪论。主要介绍了三维图像重建方法的研究意义，图像重建问题的两大类 解决方案及其发展近况，以及论文的内容安排。 基于匹配的插值方法是本文的主要介绍内容。分为三章。主要介绍了无须分割的灰 度图像层间插值算法，并介绍了离散点插值算法。这类算法分别是：基于点匹配的插值 方法、基于自由曲面形变配准的插值方法。 第六章总结了本论文的研究工作，并提出了今后的研究方向。 第二章插值方法 第二章插值方法 将插值广义地分为两类一一基于场景( s c e n e - b a s e d ) 的图像插值和基于对象 ( o b j e c t - b a s e d ) 的图像插值。前者用已有的切片体素数据直接决定被插值体素数据。而 后者则从已有的像素中获取对象信息来决定被插值体素数据，完成插值。 2 1 基于场景( s c e n e - b a s e d ) 的插值方法 2 1 1 最近邻域法 该方法为最简单、最快速、但最粗糙的插值方法，它以欧氏距离最近的已知体素 值作为待插值体素值。 2 1 2 灰度级线性插值 将未知像素灰度值用相邻切片上对应位置的两已知像素灰度值线性组合表示。假设 图像在v “及y “处的灰度值为( v “) 和f ( v “) ，则图像在未知的v 处( 三点在同 一直线上) 的灰度值可由下式近似计算 f ( v ) = a f “) + b “1 ) ( 2 1 ) 其中 4 = 瑚扣尚 而此处的模表示欧氏距离。上式常用于两相邻切片间的插值。 2 1 3 双线性插值 ( 2 2 ) 这是对最近领域法的一种改进，即待插点处的灰度值用离待插点最近的四个点像素 的灰度值加权求得。这种方法在图像内插中使用较多。 2 1 4 三次多项式 根据连续信号采样定理，若对采样值用插值函数 东南大学硕士学位论文 c ( 工) ：s i n ( 五x ) 瓜 ( 2 3 ) 则可准确恢复原函数，当然也就可以准确得到采样点间的任意点的值。三次插值法 将c ( x ) 近似为三次多项式，利用待插点周围的1 6 个点的灰度值加权求得待插点处的灰 度值。 其它多项式形式的插值还有三次样条插值、改进的三次样条插值、双三次样条插值、 三次b 一样条插值等。 如三次样条插值是将“，f ( v ”) ) ) ，k - i ，f ( v “1 ) ) ，k + l ，f ( v “) ) ， “3 ，f ( v “) ) 作为节点，构造三次样条插值插值多项式，) ，求得f ( v ) 。这种方 法常用于序列切片间的插值。 2 2 基于对象( o b j e c t b a s e d ) 的插值方法 2 2 1 二进制形式的基于形状的插值 该方法分以下几步： ( 1 ) 对二值图像施行距离变换，将每一像素赋以其到物体边界的带符号距离，得 到距离图( d i s t a n c em a p ) 。 ( 2 ) 插值此距离图到所需的分辨率， 得到目的图像的距离图。 ( 3 ) 将目的图像的距离图逆变换成二值形式的目的图像。 2 2 2 灰度形式的基于形状的插值 该方法分以下几步： ( 1 ) 将相邻两己知切片转化成三维二值图像；转化的方法用。提升”( 1 i f t i n g ) 操作， 即在每一像素处拉出一个三维坐标，其坐标值为该像素的灰度值； ( 2 ) 将距离变换应用到这两幅二值图像，即将每一像素赋以其到边界的距离； ( 3 ) 插值这两幅距离图像得到未知切片的距离图像，其中插值方法通常用线性插 值； ( 4 ) 将未知切片的距离图像转化为未知切片的三维二值图像； 第二章插值方法 ( 5 ) 将上述三维二值图像“崩塌”( c o l l a p s e ) 成二维灰度图像 2 2 3g o s h t a s b y 方法 g o s h t a s b y 等提出了先找对应点，再进行线性插值的方法。这是切片间的插值方法， 在每次插值时仅利用对相邻切片的像素信息。g o s h t a s b y 将这两片称作参考图像 ( r e f e r e n c ei m a g e ) 和目标图像( t a r g e ti m a g e ) ，并将具有大梯度值的像素称作特征点 ( f e a t u r ep o i n t ) 。算法将位于参考图像和目标图像上的特征点建立对应关系，这种对应关 系允许参考图像上的一个特征点可以和目标图像上n 个特征点建立对应( n 为非负整 数) 。然后，按每一组对应点建立线性插值， 求得未知切片上的相应像素的灰度值。 而对于梯度较小的像素应用常规的灰度形式的线性插值，求得未知切片上其它像素的 灰度值。 2 2 4 其它的基于对象的插值方法 h i g g i i l s 等1 9 9 6 年提出了基于非线性过滤( n o n l i n e a rf i l t e r i n g ) 的序列切片灰度插 值方法，比传统的插值算法( 如最近邻域法、线性、三次卷积及基于样条的插值算法) 更好地保持了区域边界。t a m 等1 9 9 8 年将二维图像每一像素的灰度值看成第三维坐 标，建立一种称为联合球( u n i o no f s p h e r e s ) 的三维模型，再根据任两幅切片图像的联 合球产生任一新切片的联合球，从而反求出此新切片这种方法较好地解决了数值不 稳定性问题一一即许多处理技术所体现的算法依赖于输入数据的情况，也无需用户过多 干预。2 0 0 0 年，c h a t z i s 等引入了形态学骨架( m o r p h o l o g i c a ls k e l e t o n ) 方法。这是一种 有效的基于形状的插值方法，l e e 等利用形态学的膨胀与腐蚀算法进行切片间的插 值2 0 0 1 年，b a l l e s t e r 等在中研究了将图像插值归结为二阶偏微分方程的变分问题的 求解方法 所有以上方法均可用作图像间插值。 关于图像内的插值，除了前述的常规插值方法( 如线性插值) 外，还有以下一些技 术1 9 9 3 年，m a h m o o d i 等给出了自适应的插值算法，对于文本区域使用5 阶插值多 项式，对于图形区域利用三次样条多项式，具有较好的处理效果。1 9 9 6 年，a p p l e d o m 在 中提出了利用高斯函数及其偶数阶导数的插值新方法g a o 等1 9 9 7 年将一幅二维图像 看成是对某二维连续函数( 称作图像场) 的采样结果将图像场( i m a g ef i e l d ) 看成能 量场，则插值问题转化成确定任一点处的势能，最后归结为沿等势能线的插值该算法 同时可用于图像放大。1 9 9 9 年，p l a z i a e 利用神经网路的b p ( b a c kp r o p a g a t i o n ) 模型， 通过自学习，改善图像插值的质量。 东南大学硕士学位论文 2 3 配准过程中的插值方法 目标图像像素经过变换后所得到的点集坐标不一定是整数，所以需要通过插值的方 法来获得变换后图像网格上的灰度值。如果是变换是线性的，那么可以直接求得其反变 换的解析式，并且变换后的点也是均匀分布的，对于这种情况采用部分体积插值方法。 如果不存在直接的逆变换且插值后的点是非均匀分布的，采用离散点插值方法。 2 3 1 部分体积插值( p a r t i a lv o l u m e - - p v ) 部分体积插值方法最早由m a t s 等人提出，主要是为了克服双线性插值方法在图像中 产生新的灰度而引起图像灰度分布发生变化的缺点，以便得到较为光滑的目标函数，有 利于优化搜索。 这种方法是对双线性插值的一种改进，它并不直接计算插值点的灰度，而是根据线 性插值的权重分配原则，将每对像素对联合直方图的贡献分散到联合直方图上与之相邻 的各个像素对上，这样联合直方图上各个像素对的频度值以非整数的幅度增加，不会像 三线性插值那样因为出现新的灰度值而破坏目标函数值分布的光滑性。其具体计算公式 为： | i l ( ，( h ) ，f ( v i ) ) = _ i l ( ，( “) ，f ( v ；) ) + h 其中f ( v ，) 为其灰度值，m 为各临近点权重，与它们到目标点的距离成反比： w lf f i ( 1 一d x ) ( 1 一d y ) w 2 = d x ( 1 一d y 、 w 3 = ( 1 一d x ) 匆 w 4 = 螂 其中d x ，d y 分别为u 与v i 之间的沿x ，y 方向的距离。 ( 2 4 1 第二章插值方法 n d y卜 晒 ， 卜 呱 图2 1p v 插值示意图 p v 插值方法严格来说并不是一种真正的插值方法，因为它并没有求得插值点的灰 度，它只是一种对灰度进行统计时用来代替插值方法的有效辅助手段。该算法对于互信 息的计算有很好的减少误差的作用。 2 3 2 邻域均值插值 对于待插值的图像上的每一点p ，在变换后的离散点集中寻找离它最近的像素点集 ( 以p 为中心的一个像素宽度范围内) ，将它们的平均值赋给p ，并将p 点标记为1 。若其 附近没有变换后的像素点，则标记为0 。第一遍遍历之后图像上可能存在一些没有灰度 值的空洞，然后再次遍历，对于被标记为0 的点，判断其四邻域是否标有记为l ，若至少 有一个，则将标记为1 的邻居灰度值作平均该点的灰度，并且将标记改为1 。如此重复扫 描几次，便可将图像中的空洞全部填满。 2 4 点对应关系引导下的插值方法 本方法隶属于基于对象的插值方法，由于篇幅关系，故另起一节。此类方法主要的 想法就是将序列切片先予以配准，再进行线性插值，算法自动根据两片之间的组织信 息进行两幅图像上像素之间的对应( 包括误配的检测及校正) 。前面的g o s h t a s b y 方 法就是属于这一类方法。 层间距过大时，采用基于场景的插值方法进行重建不能得到令人满意的效果。因为 原始切片数量较少，不能包含精确重建原始体积数据所需要的全部信息，所以采用基于 s i n c 函数发展出的插值方法会产生混叠和边缘模糊。从另外一个角度来说，这类方法最 大的不足之处在于简单地采用z 方向的若干邻近点进行加权和，这里隐含一个前提，那 就是这些点都属于同一组织器官，具有类似的灰度，也就是说它们是天然匹配好的。在 东南大学硕士学位论文 层间距较小的情况下，这个前提假设基本能够成立，当层间距大到超出了同一组织轮廓 的渐变区域，那么这类方法将不再适用。必须为每一点在临近的切片中寻找其对应点， 确保插值过程在属于同一组织的灰度点之间进行，这是提高插值效果的关键。在本章中 可以看到，当断层图像之间的距离远大于像素间的距离时，使用图像匹配方法能够大大 地提高图像插值的质量。 图2 2 两种插值方法的示意图 下面所介绍的方法都是基于点匹配的自动插值算法，这类算法大致可以分解为两个 过程：首先，根据某种算法寻找点的对应关系；然后根据对应关系引导插值过程最终得 到目标图像。寻找断层间点的匹配关系在十多年前就已经有人在研究了，其 o o s h t a s b y ， h i g g i n s 1 l 】等人提出的利用梯度信息的匹配点算法最具有代表性。最近几年随着非刚性配 准的发展，p e n n e y 提出可以利用非刚性配准来描述图像上点的对应关系。还有一大类 算法是b c v e m u d 等人0 2 提出的基于水平集运动配准的插值方法，思路比较新颖， 但仔细分析最终也可以归结为上面提到的匹配一插值思想。 既然不能直接利用z 方向的邻近点进行插值，必然会导致非整数坐标的像素点出 现，要么是在原始切片中，要么是在待插值的切片中。采用已知切片的原始对应点连线 进行插值，不可避免会在待插值层产生分布不均匀的离散数据，所以在插值后的处理过 程中会涉及到离散数据的插值或者拟合问题，本节对这一类问题也作了初步的讨论。 有了两幅灰度图像中每个点的对应关系，就可以根据对应关系对每个点进行较为准 确的插值。本章所讲的插值方法主要有两种思路，一种是主动式的，从待插值断层出发， 遍历每个像素点，寻找该点在上下两层中的对应点，然后根据层间距权重进行灰度插值， 此类算法的典型代表有p e n n e y 和h i g g i n s ；另一种是被动式的，即根据已知图像的对应点 第二章插值方法 连线进行插值，插值线与目标断层相交于一点，一般来说，这点不一定落在像素网格点 上，遍历每组像素对，最后得到一幅离散的图像数据，用某种重采样方法对此图进行插 值，使得图像点都在像素网格上，这类算法就是离散点插值算法。 2 4 1 主动插值 这类方法中以p e n n e y 的插值线算法最为典型。下面简单介绍一下他的这种算法： 。匕兰二二二石歹乏a k 二i 圣每；= e 了。鼍妒，k 二：粤唑粤些j e 。，一7 ”霉嚣铲” k 之三二蔓= 歹乏b 图2 3 计算中间层的插值点x ( 1 ) 对于待插值层上的每个像素u ，计算其世界坐标x = 五雎。 ( 2 ) 在u 周围取一搜索窗，窗口半径为最大形变r ，在该窗口内所有匹配点对的连 线记为t 。遍历。中每条直线，计算x 到直线的距离，得到最小距离插值线 k ( 3 ) 计算平行于z 。且通过x 点的直线k 。 ( 4 ) 计算k 与上下两层的交点x 和x 。( x 和x 。可能为非整数值) 。 ( 5 ) 用双线性插值计算x 和邑对应的灰度值( x ) 和( x 。) 。 根据( x ) 、l ( x 。) 和插值层到上下两层的距离，用线性插值计算x ( 即u ) 点的灰度 值。其中l ，瓦和。分别是切片a 、i 、b 上的二维点到世界坐标系的转换矩阵，上恤是 距离x 最近的插值线。插值线算法的主要工作量在于为断层上每一点寻找最小距离插值 线恤，而其实由于图像的空间连续性，在某个较小的邻域内工。是不变的，这样就带 1 1 ， 一 - k。|。，：、。_，。，泌，!=! 东南大学硕士学位论文 来了很多的冗余计算。在后面的章节里会对这种算法的性能作全面的评价。 2 4 2 被动插值 根据己知的对应关系在待插值图像上插值出一系列中间点，对于待插值图像来说， 这个过程是一个被动接受的过程，因为它真正需要的是从规格化的像素网格点出发的插 值对应关系。为了转移矛盾，需要后继算法来弥补这一点。而已知图像上的若干采样点， 如何求得规格化的图像数据，这便是离散点插值算法的作用。 离散数据插值是根据离散的或者说是非均匀分布的数据采样拟合出光滑的表面。许 多科学和工程领域都要用到这种方法。文献【”1 对这类问题的不同算法作了比较，它们在 数据连贯性、计算复杂性等方面都有较大的不同。这里仅仅介绍两种有代表性的算法： 2 4 2 1 基于d e l a u n a y 三角化的三角形内插值 这种方法适用于平面离散点较少，需要精确保证现有点灰度值的应用中。它的思想 是分而治之，首先用已知点形成的多边形( 这里是三角形) 将目标点分割成互相独立的 若干个插值单元，然后在每个插值单元内利用顶点的已知数据用某种方法进行插值。由 于它的这种思路，其插值结果在数据连贯性上不尽如人意，而且随着项点( 已知离散点) 的增加，多边形的数目将快速增加，从而导致插值时间的大幅上升。但是这种方法提出 较早，目前已经比较成熟，且能很好的保留目前已知点的灰度值，目前应用得比较多。 知道，d c l a u n a y 三角剖分具有最小内角最大以及平均形态比最大的性质，因此它是给定 区域的点集的最佳三角剖分。 s t e p1 d e l a u n a y 三角化 由于d c l a u n a y 三角网中每个三角形的外接圆都不包含其他的已知点。逐点插入算 法就利用这个性质，其核心思想埽1 是：当添加一个新点时，找出外接圆包含此点的三角 形。将这些三角形原来的边按一定规则重新与新加入的顶点组合三角形： ( 1 ) 创建一个包围所有已知点的超级三角形，作为初始三角形。 ( 2 ) 从己知点集中每次加入一个点，调用三角形分割算法。直到所有点添加完毕。 ( 3 ) 对于每一个三角形，如果与初始的超级三角形共享顶点，则将其从三角形队列 中删除。 这样便得到了平面点集的d c l a u n a y = 角网，可以用来在下面引导插值过程。 s t e p2 三角形内面积插值 第= 章插值方法 这里采用l a f a n g e 二次三角形单元内插值，具体算法如下： ( 1 ) 遍历图像上的每一点，找出包围它的三角形，以便于下面的插值；对于不属于 任何三角形的点，认为它是背景点，不作处理或简单地进行灰度插值处理。下 面的算法描述了如何判断点在三角形内： 设三角形为a b c ，目标点为p 。计算向量叉积p a x a b 、p b x b c 、p c x c a ，若这三者 的方向相一致，则点在三角形内，否则点便在三角形之外。 ( 2 ) 计算p 点与其相对边形成的三个子三角形，根据其面积确定三个顶点的灰度对p 点的贡献。把三角形三个顶点编号为1 、2 、3 ，其三条边的中点为4 、5 ，6 。根 据下面的式子可以计算出三角形单元三个顶点对应的面积比厶： 厶= a i a ，i = 1 , 2 ，3 其中， ，如， 分别为三角形p 2 3 、p 1 3 、p 2 1 的面积，a 为整个三角形的面积。 由( i _ 1 , 2 ，3 ) 计算六个节点的插值函数m ： 吼t = ( 2 一1 ) ，k = 1 , 2 ，3 仍，4 = 4 厶如，竹5 = 4 岛厶，仍6 = 4 厶厶， 这里采用的是二阶面积插值，公式如下： 6 ，( p ) = f ( x 牲) 砚，f ( p ) k = l 2 4 2 2 基于多层b 样条函数的曲面拟合算法 ( 2 7 ) ( 2 8 ) ( 2 9 ) 多层b 样条函数拟合( m u l t i l e v e lb s p l i n ea l g o d t h m m b a ) 2 0 1 ，并不是一种真正 意义上的图像插值算法，它只是尽可能光滑地拟合已知点所张成的曲面，所以它在已知 点处的取值不一定与原来的值相等，但由于采用了三次b 样条作为插值基函数，其生成 的结果图像达到了二阶连续。而且其计算复杂度较小，计算时间随离散点数的增多并无 显著增加。 令q = ( 工，) ，) l o j m ，o ) ， n ) 为x y 平面上的矩形域。3 d 空间中有一组离散点 东南大学硕士学位论文 p = ( t ，y c ，z 。) ) ，点( t ，y 。) 在q 中。为了拟合离散数据p ，用覆盖在q 上的控制网格m 所 定义的均匀三次b 样条函数建立拟合函数。 令疵为网格m 上第( i ，j ) 个控制点的值，其中i _ 一1 0 i l l + 1 ；j - 一l 0 n + 1 。 由这些控制点定义的拟合函数为： ( 2 1 0 ) 其中f _ b ， ，j 一1 ，j = l y n ，j 一1 。“= 工i n , - l x l n , l v = y l n ，- l y l n y j ，e 表示i 阶b 样条基 函数： 玩( “) = ( 1 - “) 3 6 b i ( u ) = f 3 3 - 6 u 2 + 4 ) 6 b 2 ( “) = ( - 孙3 + 3 u 2 + 3 u + 1 ) 6 b 3 ( u ) = u 3 6 其中o s t t 的点执行以下操作： 2 1 在目标图像中以p 点坐标为中心开一个半径为w 的搜索窗口，对窗i = 1 中所有 梯度大于t 的点根据下面的公式计算相异度，并选择相异度最小的点作为p 的对应点。 i c ( p k ，p k + - ) l2 = “- 2 【，( ，i ) j ：( f t + 1 ) 2 + , 2 2 【( 占( p k ) g ( p k + - ) 】2( 3 2 ) + “3 【耿p t ) 一a ( p t + 】) 】2 + “4 【( 毛一x k + l 】+ ( y i y k + 1 ) 】2 其中x l 【、y k 、x k + l 、y k + 1 分别是点p k 、p k + l 的平面坐标 2 2 根据已经得到的对应关系计算每点的水平偏移和垂直偏移： 然后分别对水平偏移和垂直偏移进行中值滤波。 s t e p3 交换参考图像和目标图像，重新执行s t e p 2 ，保留两次匹配一致的对应关系。 s t e p 4 对梯度值小于t 的点和s t e p 3 q b 丢弃的误匹配点，根据已知的高梯度点偏移用基于 d e l a u n a y - - - 角化的三角形内插值得到整幅图像的偏移关系，进而确定每点的对 应关系。 3 2 匹配参数的选择 由于不同图像的性质会对匹配的准确性产生影响，所以匹配过程需要根据匹配图像 的类型改变参数。算法中使用了梯度阈值t 来屏蔽对均质区域的匹配，文献【l 习中提到在 灰度级为2 5 5 的图像中，t 取值在5 以上比较合适。对于匹配中的另一个参数搜索窗 口的半径w ，发现窗口太小则不能在较大的偏移点之间建立对应关系，窗口太大则不仅 仅是减慢了搜寻的过程，还增加了误匹配的可能性。【1 5 】指出搜寻窗口的半径应当与层间 像素距离相等或是层间像素距离的1 2 倍，令r 为断层的像素间距，d 为断层间距，则窗口 半径可通过下式确定： 脚：1 2 一d ， ( 3 4 ) 权重因子u l , u ：，吻，决定了灰度、梯度、梯度方向和偏移在距离相异度评价函数中 的相对重要性。当被匹配图像的灰度差异很小时，灰度可以可靠地用来匹配图像，此时 l i l 应当取较大的值。如果图像具有相似的灰度变化，那么与“：和鸭就应取高的值。一般 o h 一 一 k 儿 = = 出方 东南丈学硕士学位论文 来说梯度和梯度方向均被归一化到和灰度值相同的取值范围，以便使前三项权重因子在 同一数量级上变化。由于图像的相似性主要取决于灰度和梯度信息，所以对于第四项， 只是将它作为灰度梯度都相同时的辅助判断信息。 应该说，在较大的权重因子取值范围内，这个匹配算法的表现是稳定的。在下一小 节，将使用多参数优化算法，研究不同参量对应的算法表现，并给出了对不同断层图像 进行匹配算法的结果。 3 3 实验结果及讨论 采用一组头部m r 断层来测试匹配点算法的效果，下面两节实验所用的数据也是这 组图像。进行实验的数据是2 5 6 x 2 5 6 头部m r 断层，其层内分辨率为0 9 3 7 5 m m ，层间 距为4 m m 。为了与原始图进行对比，取间隔的两幅图像进行实验，这样进行匹配的两幅 图像层间距为8 棚，其搜索窗口半径为生：8 5 3 ，梯度阈值取5 0 。 ， 图3 2 其中a , b 两幅图像为原始图像，c 为线性插值结果，d 为匹配点插值结果 第三章基于梯度信息匹配点的插值方法 为了找到使得两幅图像尽可能匹配的四项权重因子，采用p o w e i l 算法进行最优化搜 索，出发点为 0 2 ，0 4 ，0 1 ，1 6 ) ，采用匹配点插值算法结果与原始图的均方误差为代价函 数，进行最优化搜索。最后得到一组最优的权重因子： 1 8 3 ，1 2 2 ，4 2 7 ，1 6 0 ，此时总的 匹配点数目为3 0 5 ，均方误差为2 3 8 2 4 9 。 图3 2 是利用匹配信息得到的插值实验结果。 东南大学硕士学位论文 第四章基于自由曲面形变配准的插值方法 4 1 引言 相比于前面提出的点匹配算法，基于非刚体配准的方法由于利用了图像的相似性信 息，无疑要比只依赖于梯度或是距离的单点匹配算法能更好的体现断层间的对应关系。 当然基于配准的插值算法也必须满足如下两个前提：相邻断层包含相似的解剖特征，并 且配准算法能够找出这些相似特征的映射关系。r u e c k e r t 提出的算法利用f f d 使依附 于图像上的一组网格控制点发生形变，从而体现解剖特征的映射关系。根据这组形变网 格点可以得到图像的形变数据场，然后可以利用离散点插值方法进行精确的层间插值。 下图描述了图像配准的一般步骤： 4 2 图像变换 图4 - l 图像配准的一般步骤 图像配准的目的就是要找出从一幅图像到另一幅图像每一点的映射关系。相邻两层 断层的解剖组织特征基本是连续变化的，但不能完全用刚体变换来描述，需要采用非刚 体形变。r u e c k e t t 提出的算法包含一个全局刚体变换和一个局部变换： t ( 】【，”= 巳( x ，y ) + l 叫( 】【，” ( 4 1 ) 全局刚体变换描述图像的整体运动，局部变换描述图像的局部形变。局部变换也使用基 于b 样条的f f d 模型。 第网章基于自由曲面形变配准的插值方法 4 2 1 刚体变换 刚体变换是指在变换前后的两个平面中，任意两点间的距离保持不变的一种坐标变 换方法。刚体变换又可以分解为平移( t r a n s l a t i o n ) 和旋转( r o t a t i o n ) 。在二维刚体变 换中，一般包含三个变换参数：两个平移量t x 、t y 和一个旋转量0 ，其变换公式如下： m 川 痂c o s 8 口c 咖o s 口o j h l y j 嘲 4 2 2 自由曲面形变( f r e e f o r md e f o r m a t i o n ) ( 4 2 ) 图4 2 图像配准中的好形变和差形变 局部变换使用基于b 样条的f f d 模型，与m b a 方法使用的均匀三次b 样条函数 类似，所不同的是那节里主要是以求得网格控制点的值为主，通过调整网格控制点的值 来拟合己知的插值点，控制点的值代表的是图像的灰度；这里的网格控制点都已经初始 化图像上点的坐标，通过某种方法移动网格控制点以使图像上的点发生期望的形变。这 东南大学硕士学位论文 里使用的是r u c k e r t 算法的二维形式，定义力为目标图像上的点，上面附着均匀分布 的n i x n ，个网格控制点噍j 。3 = ：是f f d 可以写成1 维b 样条的2 维张量积： 33 瓦。( j ，) ，) = 且( h ) b ( v ) 丸。，州 i t = o l = 0 ( 4 3 ) 其中f - b ，j 一1 ，j = l y l n ，j 一1 ，u = x n x l x n ， v = y l n ，一l y l n ，j ，b i 表示i 阶b 样 条基函数网格控制点( c p g ) 尺寸：j | v ，= c x 8 ，+ 2 ，n ，= c y a , + 2 ( 图像尺寸“、c y ， 网格间距正，氏) 妒o = ( f 疋。j 正) ( 4 4 ) 其中一1 n x 一1 ，一1 匀 n y 一1 。控制点破，作为b 样条f f d 的参数控制着非刚体形变 的程度。点距比较大时描述了全局的非刚体形变，点距比较小时则刻画局部形变。 按照这样的定义，变换t 与这些控制点就成为了线性关系，这样初始的刚体变换a 可以叠加到变换t 的控制点上： 33 ( a o t ) ( x ，y ) =