Elastic registration of medical images using radial basis functions with compact support
- 格式:pdf
- 大小:236.83 KB
- 文档页数:6
基于结构图像表示和微分同胚 Demons 算法的多模态医学图像配准李碧草;张俊峰;杨冠羽;舒华忠【摘要】为了提高多模态医学图像的配准精度,利用参考图像和浮动图像的结构信息,提出了一种基于结构图像表示和微分同胚 Demons 算法的多模态医学图像配准算法。
该算法由结构图像表示和图像配准组成。
在结构图像表示阶段,采用Arimoto 熵图像来描述参考图像和浮动图像的结构信息,首先计算2幅图像中所有像素点周围指定大小邻域的熵值,然后把计算的熵值作为熵图像中对应点的灰度值以生成2幅熵图像。
在图像配准阶段,使用微分同胚 Demons 配准算法对2幅熵图像进行配准。
基于 Brainweb 数据库中磁共振数据的测试结果表明:与微分同胚 Demons 算法和基于香农熵的 Demons 算法相比,利用 Arimoto 熵表示图像的结构信息可以进一步提高配准精度。
%In order to improve the registration accuracy of multi-modal medical images,a multi-mo-dal medical image registration algorithm based on structural image representation and diffeomorphic Demons algorithm is proposed by using the structural information of the reference image and the float image.This algorithm contains structural image representation and image registration.In the process of structural image representation,the Arimoto entropy is used to describe the structural information of the reference image and the float image.First,the entropic values of the neighbourhoods with the designed size of all pixel points in the two images are calculated.Then,these entropic values are re-garded as the intensity values of the corresponding points in the entropy image and two entropy ima-ges can be obtained.In theprocess of image registration,the diffeomorphic Demons registration ap-proach is used to register these two entropy images.The experimental results of the MRI (magnetic resonance imaging)data in Brainweb database show that applying Arimoto entropy to represent the structural information can further improve the registration accuracy compared with the diffeomorphic Demons algorithm and the Demons method based on Shannon entropy.【期刊名称】《东南大学学报(自然科学版)》【年(卷),期】2015(000)005【总页数】5页(P851-855)【关键词】Arimoto 熵;Demons 算法;结构图像表示;图像配准【作者】李碧草;张俊峰;杨冠羽;舒华忠【作者单位】东南大学影像科学与技术实验室,南京 210096; 东南大学计算机网络和信息集成教育部重点实验室,南京 211189; 东南大学中法生物医学信息研究中心,南京 210096;东南大学影像科学与技术实验室,南京 210096; 东南大学计算机网络和信息集成教育部重点实验室,南京 211189; 东南大学中法生物医学信息研究中心,南京 210096;东南大学影像科学与技术实验室,南京 210096; 东南大学计算机网络和信息集成教育部重点实验室,南京 211189; 东南大学中法生物医学信息研究中心,南京 210096;东南大学影像科学与技术实验室,南京 210096; 东南大学计算机网络和信息集成教育部重点实验室,南京 211189; 东南大学中法生物医学信息研究中心,南京 210096【正文语种】中文【中图分类】TP391Key words: Arimoto entropy;Demons algorithm;structural image represent ation;image registration引用本文: 李碧草,张俊峰,杨冠羽,等.基于结构图像表示和微分同胚Demons算法的多模态医学图像配准[J].东南大学学报:自然科学版,2015,45(5):851-855. [doi:10.3969/j.issn.1001-0505.2015.05.007]在医学成像和计算机视觉领域,图像配准引起了广泛关注.目前,图像配准算法通常可以分为基于特征和基于密度两大类[1].基于特征的配准算法包括特征提取和特征匹配,其配准精度很大程度上取决于特征点的提取;鉴于人体组织和器官及医学成像设备自身的特点,利用常见的特征提取方法提取特征点相对较难而且精度不高.基于密度的配准算法常应用于医学图像处理中,其中最典型的是Demons算法[2].部分学者对Demons算法进行了改进,提出了组合的Demons算法[3]、对称的Demons 算法[4-5]及微分同胚Demons算法[6]等;然而,这些算法对配准图像之间的灰度值变化比较敏感,在处理来自不同成像设备的图像和亮度差异较大的图像时,配准效果不佳甚至可能配准失败.文献[7-8]利用结构图像表示,将多模态配准转化为单模态配准问题,把不同模态的图像都转化成第3种模态的图像,从而降低了灰度值的差异对配准结果带来的影响.本文利用图像的结构信息,提出了一种基于图像结构表示和Demons算法的多模态医学图像配准算法.该算法首先获得参考图像和浮动图像的结构信息,然后采用Demons算法对结构图像进行配准,从而提高多模态图像配准的精度.给定一幅图像f(x, y),其Arimoto熵定义为[9]式中,α为不等于1的正数;B为计算图像直方图的箱子数;pi为图像中灰度值等于i 的概率.文献[10]指出,当α→1时,根据罗比达法则可以得出Arimoto熵的极限等于香农熵.因此,Arimoto熵是香农熵的一种推广形式.对于图像f(x, y),计算每个像素点(x, y)周围指定大小图像块的Arimoto熵;然后,将此熵值作为熵图像A(x, y)对应点的像素值,从而得到f(x, y)的Arimoto熵图像为式中,R(i, j)表示以(x, y)为中心选取的尺寸为n×n像素的图像块,其中,n为奇数;Aα(R(i, j))表示R(i, j)的Arimoto熵.T1加权的磁共振 (MR-T1) 图像及其相应的Arimoto图像见图1.图1(a)中的2个小方块表示所选取的图像块,分别计算可得Arimoto熵,即为图1(b)中2个圆点的像素值.式(2)中参数α的取值和局部图像块R的大小决定了生成的熵图像.当α=0.25,0.75,1.00,1.50,图像块R(s, t)尺寸N=3×3,5×5,7×7,9×9,11×11像素时,图1(a)生成的Arimoto熵图像见图2.由图2可以看出,随着图像块R(s, t)尺寸的增大,熵图像的边缘逐渐变得模糊.当α<1(即α=0.25,0.75)时,熵图像的对比度较低;当α>1(即α=1.50)时,对比度则较高.因此,需要在实验中选择合适的α和参考图像S.给定参考图像和浮动图像,按照第1节中的方法得到其熵图像;然后,利用微分同胚Demons配准算法对2幅熵图像进行配准.根据式(1),熵的计算需要估计pi.常用的概率密度估计方法包括图像直方图和帕曾窗估计,且后者通常能够得到更为鲁棒的结果,尤其是在样本数较少的情况下.为了使熵图像满足结构等价的必要条件,在更新直方图时增加一个空间加权.在局部邻域块R内,直方图更新公式为式中,h(R(s, t))表示图像块R(s, t)的灰度直方图,其中,s,t分别为参考图像S和浮动图像T的灰度值;Gσ(s, t)表示标准差为σ的高斯核函数,并且将其作为图像块内 (s, t) 点处的加权值.然后,将帕曾窗估计与式(3)计算得到的直方图进行卷积,得到鲁棒的直方图,并对直方图进行归一化处理以得到概率密度.文献[2]利用气体扩散模型来执行图像配准过程,把图像的弹性配准看作一个扩散过程.Demons算法通过迭代产生位移场,即式中,d(k)表示第k次迭代的形变场;Gσ表示标准差为σ的高斯函数;v(k-1)表示第k-1次迭代的形变场更新,且式中,s为图像S的梯度.参考图像S和浮动图像T间的配准可以看成是求解二者间最优空间变换的一个优化过程.微分同胚Demons算法通常将均方误差作为相似性度量.文献[11]指出,如果仅考虑相似度项,优化过程是一个病态问题.为了避免这一问题,在能量函数中引入一个正则化项来约束形变场[11],改进后的能量函数为式中,M为S和T重叠部分的图像像素;σ1,σ2为权重,以平衡相似度项D和正则项G;d为形变场;T∘d表示采用形变场d对浮动图像T进行形变.采用形变场d的雅克比行列式作为式(6)中的正则化项,即配准算法可以看作是求如下迭代方程的最优解:迭代结束时,获得的形变场d即为2幅待配准图像的最优形变.算法1为基于Arimoto熵图像的Demons非刚体配准算法.算法1 基于Arimoto熵图像的Demons配准算法输入:S,T.输出:S,T的Arimoto熵图像AS和AM.给定形变场的初始值d(0).for k=1:Nmax利用式(6)计算当前迭代下的形变场更新v(k)利用式(5)计算当前迭代下的形变场d(k)如果E(AS,AT∘(d(k)+v(k)))满足停止条件,跳出循环得到最终形变场d;否则继续循环,直至达到最大迭代次数Nmaxend基于形变场d完成S和T的配准为了验证所提配准算法的性能,选择T1加权的MRI图像、T2加权的MRI图像和PD加权的MRI图像进行实验.为简便起见,这3种图像分别记为MR-T1,MR-T2和MR-PD.测试图像来源于蒙特利尔神经医学部的Brainweb数据库[12],3幅图像的尺寸均为181×217×181像素,在物理坐标下每个体素尺寸为1 mm×1mm×1 mm.此外,这3种图像的所有层都是相互对应的.将这3种MRI图像的中间层作为测试图像(见图3(a)~(c)).当α=0.50,N=5×5像素时,其对应的Arimoto熵图像见图4(d)~(f).为了检验参数α对配准算法的影响,利用3种MRI图像设计了3组配准实验:① MR-T1和MR-T2;② MR-T1和MR-PD;③ MR-T2和MR-PD.在每组实验中,前者为参考图像,通过如下计算对后者进行形变:式中,(x′,y′)为形变后坐标;m为形变量.在每组配准实验中,m取3个不同的值以产生3幅浮动图像,由此共产生9对测试图像.对于每对测试图像,计算α=0.01,0.25,0.50,0.99,1.25,1.50时的配准结果.为了衡量配准结果,将采用配准算法得到的形变与真实形变量间的差值作为配准误差.表1给出了N=5×5像素时不同α取值下的配准误差.由表可知,α=0.50时,利用所提算法能够得到相对较低的配准误差.因此,在以下实验中,选取α=0.50.然后,验证图像块尺寸对配准算法的影响.选取MR-T1作为参考图像,形变后的MR-PD为浮动图像.同上,选取m=2,3,4以生成3幅浮动图像,选择α=0.50以生成Arimoto熵图像,图像块尺寸从3×3像素增大到17×17像素.不同图像块尺寸下的配准结果见表2.由表可知,当m=2,3时,配准误差随图像块尺寸的增大而增大;而当m=4时,N=9×9像素时获得的配准误差最小.综上可知,N=5×5像素时的图像块能获得折中的效果.为了进一步测试所提算法在处理多模态图像配准中的优势,分别采用微分同胚Demons算法[11]、基于香农熵的Demons算法[8]进行非刚体配准实验.实验中选择α=0.50,N=5×5像素.3种算法的配准误差结果见表3.由表可知,所提算法的配准误差最小;微分同胚Demons算法的配准效果最差,在配准MR-T1和MR-T2数据时甚至出现了失败.因此,所提算法在多模态图像配准方面明显优于其他2种算法.为了更直观地比较3种算法的配准结果,将图4(a)和(b)分别作为参考图像和浮动图像进行配准实验,配准前的棋盘格图像见图4(c),利用3种算法配准后的棋盘格图像分别见图4(d)~(f). 针对多模态医学图像的配准问题,提出了一种基于结构图像表示的微分同胚Demons算法.首先,利用Arimoto熵计算参考图像和浮动图像的熵图像,用熵图像表示2幅图像的结构信息;然后采用微分同胚Demons算法配准熵图像.与常规的微分同胚Demons相比,所提算法采用图像的结构信息,而非灰度信息计算图像间的相似度.实验结果表明,与微分同胚Demons算法和基于香农熵的Demons算法相比,所提算法能够获得更高的配准精度.下一步的研究工作将考虑三维医学图像的弹性配准.【相关文献】[1]Khader M, Hamza A B. Nonrigid image registration using an entropic similarity [J]. IEEE Transactions on Information Technology in Biomedicine, 2011, 15(5): 681-690.[2]Thirion J P. Image matching as a diffusion process: an analogy with Maxwell’s demons [J]. Medical Image Analysis, 1998, 2(3): 243-260.[3]Cachier P, Bardinet E, Dormont D, et al. Iconic feature based nonrigid registration: the P ASHA algorithm [J]. Computer Vision and Image Understanding, 2003, 89(2/3): 272-298.[4]Wang H, Dong L, O’Daniel J, et al. Validation of an accelerated ‘demons’ algorithm for deformable image registration in radiation therapy [J]. Physics in Medicine and Biology , 2005, 50(12): 2887-2905.[5]Rogelj P, Kovacic S. Symmetric image registration[J].Medical Image Analysis, 2006, 10(3 ): 484-493.[6]Vercauteren T, Pennec X, Perchant A, et al. Diffeomorphic demons: efficient non-parametric image registration [J]. NeuroImage, 2009, 45(1): S61-S72.[7]Wachinger C, Navab N. Structural image representation for image registration [C]//Proc eedings of 2010 IEEE Computer Society Conference on Computer Vision and Pattern Reco gnition Workshops. San Francisco, CA, USA, 2010: 23-30.[8]Wachinger C, Navab N. Entropy and Laplacian images: structural representations for m ulti-modal registration [J]. Medical Image Analysis, 2012, 16(1): 1-17.[9]Arimoto S. Information-theoretical considerations on estimation problems [J]. Information and Control, 1971, 19(3 ): 181-194.[10]Boekee D E, van der Lubbe J C A. The R-norm information measure [J]. Information and Control, 1980, 45(2): 136-155.[11]Vercauteren T, Pennec X, Perchant A, et al. Non-parametric diffeomorphic image registration with the demons algorithm [C]//2007 Medic al Imaging and Computer-Assisted Intervention. Brisbane, Australia, 2007: 319-326. [12]Kwan R K S, Evans A C, Pike G B. MRI simulation-based evaluation of image-processing and classification methods [J]. IEEE Transactions on Medical Imaging, 1999, 18( 11): 1085-1097.。
elastix_registration_method -回复elastix_registration_method——基于弹性变形的图像配准方法图像配准是计算机视觉领域的关键问题之一,它旨在将两幅或多幅图像的位置、尺度、旋转或其他形变参数进行匹配,以便在后续的图像处理或分析任务中实现更准确的结果。
elastix_registration_method(弹性变形配准方法)是一种常用的图像配准算法,特点是能够提供更精确的配准结果。
一、弹性变形配准的背景和意义图像配准在医学影像处理、地理信息系统、计算机视觉、机器人导航等领域扮演着重要角色。
然而,由于图像间的形变和变形,所得到的结果往往不够精确,而且传统的配准方法更多地依赖于刚性变换模型,无法适应非刚性变换。
因此,弹性变形配准应运而生。
弹性变形配准算法可以通过对图像进行弹性变形,实现更准确的配准。
它模拟了物体在外部作用下的弹性变形过程,将刚性变形模型扩展为非刚性变形模型。
弹性变形配准方法的基本思路是通过对图像进行形变,使得相应的图像特征能够达到最佳匹配,从而实现图像的准确配准。
弹性变形配准方法在医学影像处理中得到广泛应用。
例如,对于脑部MRI 图像的配准,在手术前后或多个时间点的比较中,可以帮助医生更准确地确定病变区域的位置和变化。
此外,弹性变形配准还在图像分析、形状建模以及计算机视觉等领域发挥着重要作用。
二、elastix_registration_method的基本原理elastix_registration_method是一个开源的弹性变形图像配准框架,基于C++实现。
它结合了B-spline弹性变形模型和优化算法,能够实现高效准确的图像配准。
1. 弹性变形模型elastix_registration_method采用了B-spline弹性变形模型。
B-spline 是一种常用的曲线和曲面拟合技术。
它的基本思想是通过一组控制点来表示整个图像的形状,而变形过程中的每个点都可以根据控制点的位置得到相应的弹性变形。
一种基于轮廓的医学图像弹性配准方法
舒小华;沈振康
【期刊名称】《计算机工程与应用》
【年(卷),期】2007(043)014
【摘要】提出了一种基于轮廓的医学图像弹性配准方法.该方法先对医学图像进行轮廓提取,计算轮廓所围部分(目标)的质心.围绕质心每隔一定角度抽得轮廓点作为特征点.通过计算相似度确定两组特征点的对应关系.采用Wendland提出的ψ函数作为弹性变换的基函数.通过两组对应点确定变换系数.最后根据配准参数对形变图像进行变换,从而得到配准图像.实验表明,本方法对具有形变的医学图像的配准是快速的、有效的.
【总页数】3页(P190-191,248)
【作者】舒小华;沈振康
【作者单位】国防科技大学,ATR国家重点实验室,长沙,410073;湖南工业大学,电气工程系,湖南,株洲,412008;国防科技大学,ATR国家重点实验室,长沙,410073
【正文语种】中文
【中图分类】TP309
【相关文献】
1.一种医学图像弹性配准方法 [J], 于佳;陆丹
2.一种基于轮廓形状的胸部医学图像插值方法 [J], 屈景怡;史浩山
3.基于层次B样条的医学图像弹性配准方法 [J], 张红颖;张加万;孙济洲;杨甲东
4.基于活动轮廓模型的脑部医学图像弹性配准 [J], 唐祚;闫德勤
5.基于互信息和薄板样条的医学图像弹性配准方法研究 [J], 曹国刚;罗立民
因版权原因,仅展示原文概要,查看原文内容请购买。
Intensity gradient based registration and fusionof multi-modal imagesEldad Haber1and Jan Modersitzki21Mathematics and Computer Science,Emory University,Atlanta,GA,USA,haber@mathcs.emory,edu2Institute of Mathematics,L¨u beck,Germany,modersitzki@math.uni-luebeck.de Abstract.A particular problem in image registration arises for multi-modal images taken from different imaging devices and/or modalities.Starting in1995,mutual information has shown to be a very successfuldistance measure for multi-modal image registration.However,mutualinformation has also a number of well-known drawbacks.Its main dis-advantage is that it is known to be highly non-convex and has typicallymany local maxima.This observation motivate us to seek a different image similarity mea-sure which is better suited for optimization but as well capable to handlemulti-modal images.In this work we investigate an alternative distancemeasure which is based on normalized gradients and compare its perfor-mance to Mutual Information.We call the new distance measure Nor-malized Gradient Fields(NGF).1IntroductionImage registration is one of today’s challenging medical image processing prob-lems.The objective is tofind a geometrical transformation that aligns points in one view of an object with corresponding points in another view of the same object or a similar one.An important area is the need for combining informa-tion from multiple images acquired using different modalities,sometimes called fusion.Typical examples include the fusion of computer tomography(CT)and magnetic resonance(MRI)images or of CT and positron emission tomography (PET).In the past two decades computerized image registration has played an increasingly important role in medical imaging(see,e.g.,[2,10]and references therein).One of the challenges in image registration arises for multi-modal images taken from different imaging devices and/or modalities.In many applications, the relation between the gray values of multi-modal images is complex and a functional dependency is generally missing.However,for the images under con-sideration,the gray value patterns are typically not completely arbitrary or random.This observation motivated the usage of mutual information(MI)as a distance measure between two images cf.[4,18].Starting in1995,mutual in-formation has shown to be a successful distance measure for multi-modal imageregistration.Therefore,it is considered to be the state-of-the-art approach to multi-modal image registration.However,mutual information has a number of well-known drawbacks;cf.e.g.,[14,13,12,17,20].Firstly,mutual information is known to be highly non-convex and has typically many local maxima;see for example the discussion in [10],[20],and Section3.Therefore,the non-convexity and hence non-linearity of the registration problem is enhanced by the usage of mutual information. Secondly,mutual information is defined via the joint density of the gray value distribution,and therefore,approximations of the density are required.These approximations are nontrivial to compute and typically involve some very sen-sitive smoothing parameters(e.g.a binning size or a Parzen window width,see [16]).Finally,mutual information decouples the gray value from the location information.Therefore judging the output of the registration process is difficult. These difficulties had stem a vast amount of research into mutual information registration,introducing many nuisance parameters to help and bypass at least some of the difficulties;see,e.g,[14].These observations motivate us to seek a different image similarity measure which is capable to handle multi-modal images but better suited for optimization and interpretation.In this paper we investigate an alternative distance measure which is based on normalized gradients.As we show,the alternative approach is deterministic,simpler,easier to interpret,fast and straightforward to implement, faster to compute,and also much more suitable to optimization.The idea of using derivatives to characterize similarity between images is based on the observation that image structure can be defined by intensity changes. The idea is not new.In inverse problems arising in geophysics,previous work on joint inversion[9,21,7]discussed the use of gradients in order to solve/fuse inverse problems of different modalities.In image registration,a more general framework similar to the one previously suggested for joint inversion was given in[6]and[14].Our approach is similar to[14]however,we have some main differences.Unlike the work in[14]we use only gradient information and avoid using mutual information.Furthermore,our computation of the gradient is less sensitive and allow us to deal with noisy images.This allows us to use simple optimization routines(e.g.Gauss Newton).The rest of the paper is organized as follows:In Section2we shortly lay the mathematical foundation of image registration.Section3presents an illustrative example showing some of the drawbacks of mutual information.In Section4we discuss the proposed alternative image distance measure.Finally,in Section5 we demonstrate the effectiveness of our method.2The mathematical settingGiven a reference image R and a template image T,the goal of image registration is tofind a“reasonable”transformation,ϕ,such that the“distance”between the reference image and a deformed template image is small.In this paper we use both affine linear and nonparametric registration.In the linear registration,we set the transformationϕtoϕ(γ,x)= γ1γ2γ3γ4x1x2+γ5γ6.(1)Given a distance measure D,the registration problem is tofind a minimizerγoff(γ):=DR(·),T(ϕ(γ,·)).(2)For nonparametric registration we setϕ=x+u(x)and use regularizationminimizingf(γ):=DR(x),T(x+u))+αS(u).(3)whereαS(u)is a regularization term such as elastic(see[10]).3An illustrative exampleTo emphasize the difficulty explained above,we present an illustrative example. Figure1shows a T1and a T2weighted magnetic resonance image(MRI)of a brain.Since the image modalities are different,a direct comparison of gray values is not advisable and we hence study a mutual information based approach.Figure1c)displays our approximation to the joint density which is based on a kernel estimator,where the kernel is a compactly supported smooth func-tion.Note that the joint density is completely unrelated to the spatial image content(though there is some interpretation[11]).We now slide the template along the horizontal axis.In the language of equation(1),wefixγ1,...,5and obtain the transformed image by changingγ6.Figure1e)shows the negative mutual information versus the shift ranging from-2to2pixels.Thisfigure clearly demonstrates that mutual information is a highly non-convex function with respect to the shift parameter.We have experimented with different interpolation methods and all have resulted in similar behavior. In particular,the curve suggests that there are many pronounced local minima which are closed in value to the global minimum.Our particular example is not by any means pathologic.Similar,and even less convex curves of Mutual Information appear also in[20]pp293who used different interpolation.Figure1d)displays a typical visualization of our alternative NGF distance between R and T(discussed in the next section).Note that for the alternative distance measure,image differences are related to spatial positions.Figure1f) shows the alternative distance measure versus the shift parameter.For this par-ticular example,it is obvious that the alternative measure is capable for multi-modal registration and it is much better suited for optimization.4A simple and robust alternativeThe alternative multi-modal distance measure is based on the following simple though general interpretation of similarity:(a)(b)(c)(d)(e)(f)Fig.1.Distance measures versus shifts;(a)Original BrainWeb [3]T1(b)Original Brain-Web [3]T2(c)the joint density approximation for the images,(d)the normalized gradient field d d (T,R )for the images,(e)negative mutual information versus shift,(f)normalized gradient field versus shift.two image are considered similar,if intensity changes occur at the same locations.An image intensity change can be detected via the image gradient.Since the mag-nitude of the gradient is dependent upon the modality of the image,it would be unwise to base an image similarity measure on gradient magnitude.We therefore consider the normalized gradient field,i.e.,the local orientation,which is purely geometric information.n (I,x ):=∇I (x )∇I (x ) ,∇I (x )=0,0,otherwise.(4)For two related points x in R and ϕ(x )in T or,equivalently,x in T ◦ϕ,we look at the vectors n (R,x )and n (T ◦ϕ,x ).These two vectors form an angle θ(x ).Since the gradient fields are normalized,the inner product (dot-product)of the vectors is related to the cosine of this angle,while the norm of the outer product (cross-product)is related to the sine.In order to align the two images,we can either minimize the square of the sine or,equivalently,maximize thesquare of the cosine.Thus we define the following distance measuresd c(T,R)= n(R,x)×n(T,x) 2;D c(T,R)=12Ωd c(T,R)d x,(5)d d(T,R)= n(R,x),n(T,x) 2;D d(T,R)=−12Ωd d(T,R)d x,(6)Note that from an optimization point of view,the distances D c or D d are equiv-alent.The definition of the normalized gradientfield(4)is not differentiable in areas where the image is constant and highly sensitive to small values of the gra-dientfield.To avoid this problem we define the following regularized normalized gradientfieldsn E(I,x):=∇I(x)∇I(x) E, ∇I(x) E:=∇I(x) ∇I(x)+E2.(7)In regions where E is much larger than the gradients the maps n E(I,x)are almost zero and therefore do not have a significant effect on the measures D c or D d,respectively.Indeed,the choice of the parameter E in(7)answers the question,“what can be interpreted as a jump?”.As suggested in[1],we propose the following automatic choice:E=ηVΩ|∇I(x)|d x,(8)whereηis the estimated noise level in the image and V is the volume of the domainΩ.The measures(5)and(6)are based on local quantities and are easy to com-pute.Another advantage of these measures is that they are directly related to the resolution of the images.This property enables a straightforward multi-resolution approach.In addition,we can also provide plots of the distancefields d c and d d,which enables a further analysis of image similarity;see,e.g.,Fig-ure1d).Note that in particular d c=0everywhere if the images match perfectly. Therefore,if in some areas the function d c takes large values,we know that these areas did not register well.5Numerical ExperimentsIn our numerical experience,we have use various examples with results along the same lines.Since it is impossible to show every result we restrict ourselves to three illustrative and representative examples.In thefirst example we use the images in Figure1.We take the T1image and generate transformed versions of the image by using an affine linear trans-formation(1).We then use the transformed images to recoverγand the original image.The advantage of this synthetic experiment is that it is controlled,i.e., we know the exact answer to the registration problem and therefore we are able(a)(b)(c)(d)(e)(f)Fig.2.Experiments with Viola’s example;(a)reference R ,(b)template T ,(c)reg-istered T ,(d)overlay of T and R (202pixels checkerboard presentation),(e)cross product n T ×n R ,(f)joint density at the minimum.to test our algorithm under a perfectly controlled environment.We randomly choose γand then run our code.Repeating the experiment 100times we find that we are able to converge to the correct γwith accuracy of less than 0.25pixels.In the second example we use the images from Viola’s Ph.D thesis [18].We compare the results of both MI and our NGF registration.The difference between the MI registration and the NGF registration was less than 0.25of a pixel,thus we conclude that the methods give virtually identical minima.However,to obtain the minima using MI we needed to use a random search technique to probe the space which is rather slow.The global MI minima has the value of about −9.250×10−2while the guess γ=0has the value of about −9.115×10−2.Thus the “landscape”of the MI function for this example is similar to the one plotted in Figure 1.In comparison,our NGF algorithm used only 5iteration on the finest grid.The registration was achieved in a matter of seconds and no special space probing was needed to obtain the minima.The value of the NGF function at γ=0was −4.63×101while at the minima its value was −2.16×102thus our minima is much deeper compared with the MI minima.The results of our experiments are also presented in Figure 2.(a)(b)(c)(d)Fig.3.Experiments with PET/CT images:(a)reference R,(b)template T,(c)MI registered T and and deformed grid,(d)NGF registered T and deformed grid;Another advantage of our method is the ability to quickly evaluate the reg-istration result by looking at the absolute value of the cross-product|n T×n R|. This product has the same dimension as the images and therefore can be viewed in the same scale.If the match is perfect then the cross product should vanish and therefore any deviation from zero implies a imperfectfit.Such deviations are expected to be present due to two different imaging processes,the noise in the images and possible additional nonlinear features.Figure2e)shows the cross product for the image matched above.It is evident that the matching is very good besides a small number of locations where we believe the image to be noisy.The previous two examples demonstrate the ability of our new algorithm to successfully perform parametric registration.Nevertheless,when the number of parameters is very large then stochastic optimization techniques stall and differentiable smooth distance measure functions are much more important.We therefore test our algorithm on a PET/CT registration,see Figure3.We perform an elastic registration of CT and a PET images(transmission)of a human thorax.For this example,the differences between the MI and NGF approach are very small;cf.Figure3.For the NGF no special effort was needed in order to register the images.For MI we needed a good starting guess to converge to a local minima.References1.U.Ascher,E.Haber and H.Haung,On effective methods for implicit piecewisesmooth surface recovery,SISC28(2006),no.1,339–358.2.Lisa Gottesfeld Brown,A survey of image registration techniques,ACM ComputingSurveys24(1992),no.4,325–376.3. C.A.Cocosco,V.Kollokian,R.K.-S.Kwan,and A.C.Evans,BrainWeb MR sim-ulator,Available at http://www.bic.mni.mcgill.ca/brainweb/.4. A.Collignon,A.Vandermeulen,P.Suetens,and G.Marchal,multi-modality med-ical image registration based on information theory,Kluwer Academic Publishers: Computational Imaging and Vision3(1995),263–274.5.J.E.Dennis and R.B.Schnabel,Numerical methods for unconstrained optimizationand nonlinear equations,SIAM,Philadelphia,1996.6.M.Droske and M.Rumpf,A variational approach to non-rigid morphological reg-istration,SIAM Appl.Math.64(2004),no.2,668–687.7.L.A.Gallardo and M.A.Meju,Characterization of heterogeneous near-surface ma-terials by joint2d inversion of dc resistivity and seismic data,Geophys.Res.Lett.30(2003),no.13,1658–1664.8.G.Golub,M.Heath,and G.Wahba,Generalized cross-validation as a method forchoosing a good ridge parameter,Technometrics21(1979),215–223.9. E.Haber and D.Oldenburg,Joint inversion a structural approach,Inverse Prob-lems13(1997),63–67.10.J.Modersitzki,Numerical methods for image registration,Oxford,2004.11.H.Park,P.H.Bland,K.K.Brock,and C.R.Meyer,Adaptive registration usinglocal information measures,Medical Image Analysis8(2004),465–473.12.Josien PW Pluim,J.B.Antoine Maintz,and Max A.Viergever,Image registrationby maximization of combined mutual information and gradient information,IEEE TMI19(2000),no.8,809–814.13.J.P.W.Pluim,J.B.A.Maintz,and M.A.Viergever,Interpolation artefacts in mu-tual information based image registration,Proceedings of the SPIE2004,Medical Imaging,1999(K.M.Hanson,ed.),vol.3661,SPIE,1999,pp.56–65.14.,Mutual-information-based registration of medical images:a survey,IEEETransactions on Medical Imaging22,1999,986–1004.15.L.Rudin,S.Osher and E.Fatemi,Nonlinear total variation based noise removalalgorithms,Proceedings of the eleventh annual international conference of the Cen-ter for Nonlinear Studies on Experimental mathematics:computational issues in nonlinear science,1992,259–268.16.R.Silverman,Density estimation for statistics and data analysis,Chapman&Hall,1992.17.M.Unser and P.Th´e venaz,Stochastic sampling for computing the mutual infor-mation of two images,Proceedings of the Fifth International Workshop on Sam-pling Theory and Applications(SampTA’03)(Strobl,Austria),May26-30,2003, pp.102–109.18.Paul A.Viola,Alignment by maximization of mutual information,Ph.D.thesis,Massachusetts Institute of Technology,1995.19.G.Wahba,Spline models for observational data,SIAM,Philadelphia,1990.20.Terry S.Yoo,Insight into images:Principles and practice for segmentation,regis-tration,and image analysis,AK Peters Ltd,2004.21.J.Zhang and F.D.Morgan,Joint seismic and electrical tomography,Proceedingsof the EEGS Symposium on Applications of Geophysics to Engineering and Envi-ronmental Problems,1996,pp.391–396.。
Proc.ComputerVisonandPatternRecognition(CVPR’99),FortCollins,CO,USA,23-25June1999,IEEEComputerSocietyPR00149,pp.402–407
ElasticRegistrationofMedicalImagesUsingRadialBasisFunctionswithCompactSupport
MikeFornefett,KarlRohr,andH.SiegfriedStiehlUniversit¨atHamburg,FachbereichInformatik,ArbeitsbereichKognitiveSystemeVogt-K¨olln-Str.30,D-22527Hamburg,GermanyEmail:fornefett@informatik.uni-hamburg.de
AbstractWeintroduceradialbasisfunctionswithcompactsup-portforelasticregistrationofmedicalimages.Withthesebasisfunctionstheinfluenceofalandmarkontheregistra-tionresultislimitedtoacirclein2Dand,respectively,toaspherein3D.Therefore,theregistrationcanbelocallyconstrainedwhichespeciallyallowstodealwithratherlo-calchangesinmedicalimagesdueto,e.g.,tumorresec-tion.AnimportantpropertyoftheusedRBFsisthattheyarepositivedefinite.Thus,thesolvabilityoftheresultingsystemofequationsisalwaysguaranteed.Wedemonstrateourapproachforsyntheticaswellasfor2Dand3Dtomo-graphicimages.
1.IntroductionRegistrationisanimportanttechniqueinmedicalimageanalysis.Rigidandaffineregistrationmethodscanonlycopewithglobaldifferences,forexample,translation,ro-tation,andscaling.Inmanycases,however,elasticornon-rigidmethodsarerequiredtocopewithlocaldifferencesbetweentheimages.Suchdifferencesaredueto,forexam-ple,scanner-induceddeformations,movementsofthepa-tient,surgicalinterventions,ordifferentanatomy(e.g.im-ageatlas-registration).Inthispaper,weconsiderapoint-basedelasticregistra-tionapproachbasedontheradialbasisfunction(RBF)in-terpolationmethod.Withthisapproachthetransformationiscomposedofradiallysymmetricfunctionsthatserveasbasisfunctions.ThechoiceofthetypeoftheRBFiscru-cialfortheoverallcharacteristicssuchasthesmoothnessorthelocalityofthetransformationfunction.Bookstein[2]hasintroducedthin-platesplinesformed-icalimageregistration.Thisapproachyieldsminimalbendingenergypropertiesmeasuredoverthewholeim-age,butthedeformationisnotlimitedtoregionswherethepointlandmarksareplaced.Thisbehaviourisadvan-tageousforyieldinganoverallsmoothdeformation,butitisproblematicwhenratherlocaldeformationslimitedtoimagepartsaredesired.Tocopewithlocaldeformations,thelandmarkshavetobewell-distributedovertheimagestopreventdeformationsinregionswherenochangesaredesired[1].OthershaveinvestigatedmultiquadricsasRBFsforreg-istration,e.g.[5],andforimagedeformations[7].TheseRBFshaveaparameterwhichcontrolstheirlocality.How-ever,thefunctionvaluesofmultiquadricsareincreasingwithgrowingdistancefromthelandmarkpositionandthustheregistrationresultatlocationsfaroffthecenteroftheRBFislargelyinfluenced.OtherRBFsdecreasewithgrowingdistancefromthelandmarkpositionsuchasin-versemultiquadrics,e.g.[7],andtheGaussian,e.g.[1].SincetheseRBFsasymptoticallyapproachzero,theglobalinfluenceisreduced,butitisnotspatiallylimited,i.e.theseRBFshavenocompactsupport.Inthispaper,weintroduceRBFswithcompactsupportfortheregistrationofmedicalimages.ThebasisfunctionsweemployhaveasimilarshapeastheGaussian,buttheyhavetheadvantagethattheirinfluenceislimitedaroundalandmark(in2Dand3Dimagesonacircleorasphere,resp.).Thispropertyallowstheregistrationofmedicalimageswherechangesoccuronlylocally.Theapplica-tionscenariowehaveinmindistheregistrationoflocalchangesinmedicalimagesduetotheresectionofatu-mororduetoothersurgicalinterventions.Thisapproachhasalsoverynicetheoreticalproperties.Actually,forthebasisfunctionsweuse,itcanbeshownthattheresultingsystemofequationsisalwayssolvable.Thus,weprovideananswertoapreviouslyposedquestionin[1],whereGaussian-shapedRBFswithcompactsupportaresoughtwhilesolvabilityisalwaysensured.Below,wefirstgiveanoverviewofthegeneralschemeforregistrationbasedonRBFs(Sect.2).InSect.3,weintroduceanelasticregistrationschemeusingRBFswithcompactsupportanddiscussitsproperties.Finally,inSect.4,wepresentexperimentalresultsfor2Dand3Dimages.2.ImageregistrationwithRBFInthissection,webrieflydescribetheradialbasisfunc-tioninterpolationschemeanddiscussitspropertiesde-pendingonthechoiceofthebasisfunction.
2.1.GeneralschemeGenerally,inregistrationapplicationsonehastodeter-mineatransformationfunctionwhereistheimagedimension,e.g.for2Dand3Dimages,resp.Aninterpolationtransformationfunctionbasedonpoint-landmarksmustfulfillthefollowingconstraints:
(1)whereconstituteagivensetofpoint-landmarks
inthesourceimageandarethecorrespondinglandmarksinthetargetimage.Often,eachcoordinateofthetransformationfunctioniscalculatedseparately,i.e.theinterpolationproblemissolvedforeachcoordinatewiththecorrespondingconstraints.Inthefollowing,wewriteinsteadof.In2D,iscalculatedseparatelyforandandin3Dfor,,and.Ifweapplyaradialbasisfunctionapproach,thentheinterpolationfunctiongenerallyconsistsoftwoparts: