基于ZMLD与GC离散优化的非刚性多模态医学图像配准方法与流程

本发明涉及医学图像处理技术领域，特别是涉及基于zmld与gc离散优化的非刚性多模态医学图像配准方法。

背景技术：

在临床医学中，不同的成像模式可以提供不同的生理信息。单模态医学图像提供的信息往往是有限的。多模态医学图像配准有利于将不同模态图像之间的信息互补，信息互补的图像提供病变组织或器官的多种信息，为医生做出准确的诊断提供有力的理论依据。

目前针对多模态医学图像配准中的相似性测度的计算方法主要分为两类。一类是使用信息论度量作为相似性测度，互信息(mutualinformation，mi)是目前广泛使用的信息论度量，它利用图像的灰度信息来计算两幅图像的相似度。但是mi忽略了图像的局部特征和结构信息，导致多模态图像配准精度降低。另一类方法是用局部描述符提取不同模态的结构信息，从而把多模态配准转化为单模态配准，使用简单相似性测度进行配准。wachinger等提出了两种用于多模态图像配准的结构表示方法。一种方法是在图像中取每个像素的邻域块，计算邻域块的熵(即该点的邻域结构信息)，将不同模态的图像转化为相同模态的熵图，并使用基于熵图像的误差平方和(sumofsquareddifferencesonentropyimages，essd)作为相似性测度。该方法计算速度快但是熵图像较模糊。另一种方法使用拉普拉斯特征映射，高阶流形通过构建邻域图进行降维，然后计算拉普拉斯图的l2距离。该方法配准精度高但是计算成本非常大，并且特征降维也损失了图像信息。heinrich等提出模态独立邻域描述符(modalityindependentneighborhooddescriptor，mind)用于多模态图像配准，根据相邻图像块之间的相似性计算mind，对非功能强度关系和图像噪声具有较好的鲁棒性。但mind不具有旋转不变性，在图像边缘和复杂纹理区域图像特征存在旋转时，mind会出现配准误差。相比mind，zernike矩具有旋转和平移不变性，以及对噪声的鲁棒性。因此，zernike矩已被广泛应用于图像处理，计算机视觉和模式识别领域。

图像配准问题可看成是能量函数的极值求解问题。通过构造能量函数，采用优化方法求解最小值，则最小值对应最优的配准效果。优化方法分为两类：连续优化和离散优化。连续优化常用算法有梯度下降法(gd)，共轭梯度下降法(cg)，拟牛顿方法(qn)等。连续优化算法大部分依赖于目标函数梯度的计算，导数的计算量较大，并且易陷入局部最小值。基于马尔可夫随机场(mrf)的离散优化策略用来克服连续优化的缺点。离散优化本质上是无梯度的，计算复杂度相对较低，并且可通过较大的邻域搜索空间进行优化，有效避免陷入局部最小值。jiansun等使用置信传播法(bp)来解决立体匹配问题。bp是一种高效的算法，但是计算复杂度较高。wainwright等提出树重加权消息传递法(trw)用于能量最小化。与bp相比，trw可用于更多的能量函数，但trw不能保证其完全收敛。glocker等使用mrf和线性规划(lp)的优化算法，把图像配准问题看做一个离散的三维标签问题。但是lp算法需要较大空间容量，这将限制lp对复杂变形的图像进行精确的配准。boykov等提出了一种交互式的图割法(graphcuts，gc)。gc是一种基于图论的组合优化方法，采用最大流/最小割(max-flow/min-cut)理论来求mrf能量的全局最优解。kolmogorov和rother等比较了常用的离散优化算法，并且得出gc法要优于其他优化算法。

针对非刚性图像存在噪声和强度失真时，现有方法无法同时准确提取图像混合信息，连续优化计算复杂度相对较高且易陷入局部最小值的问题。

技术实现要素：

本发明实施例提供了基于zmld与gc离散优化的非刚性多模态医学图像配准方法，可以解决现有技术中的问题。

本发明提供了基于zmld与gc离散优化的非刚性多模态医学图像配准方法，该方法包括以下步骤：

读取输入的待配准的参考图像i和浮动图像j，两幅图像的分辨率相同；

分别计算图像i和j两幅图像基于zernike矩的局部描述符zmld；

使用图像i和j的zmld之间的绝对误差和sad作为能量函数的数据项，采用位移矢量场的一阶导数作为平滑项，构造能量函数，并对能量函数进行离散化；

利用图割法gc的α扩展优化算法求解离散化后能量函数的最小值；

输出能量函数最小值对应的最佳位移矢量场，即配准后的图像。

本发明实施例中的基于zmld与gc离散优化的非刚性多模态医学图像配准方法，解决了非刚性图像存在噪声和强度失真时，现有方法无法同时准确提取图像强度和边缘、纹理特征，连续优化计算复杂度相对较高且易陷入局部最小值的问题。通过多模态医学图像数据集的实验表明，本发明方法提高了非刚性多模态医学图像配准的精度和效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为brainweb数据库中mr-t1图像，及mr-t1图像的a00、a11和a22图像信息；

图2是本发明实施例中基于zmld与gc离散优化的非刚性多模态医学图像配准方法的流程图；

图3为图的简单构造示意图；

图4为rire数据库中一组脑部mr的t1和t2、pd加权图像，其中(a)为参考t1图像，(b)为浮动t2图像，(c)为浮动pd图像，(d)为essd方法(t2-t1)，(e)为mind方法(t2-t1)，(f)为zmld方法(t2-t1)，(g)为essd方法(pd-t1)，(h)为mind方法(pd-t1)，(i)为zmld方法(pd-t1)；

图5为atlas数据集中的两组待配准图像,配准结果以及配准后图像差，其中(a)为参考t2图像，(b)为浮动ct图像，(c)为参考spect图像，(d)为浮动t2图像，(e)为ffd-lbfgs算法(ct-t2)，(f)为ffd-lbfgs算法(ct-t2)图像差，(g)为本发明方法(ct-t2)，(h)为本发明方法(ct-t2)图像差，(i)为ffd-lbfgs算法(t2-spect)，(j)为ffd-lbfgs算法(t2-spect)图像差，(k)为本发明方法(t2-spect)，(l)为本发明方法(t2-spect)图像差；

图6为使用基于zmld的局部描述符计算相似性测度，分别在gc和lp离散优化算法下的配准结果，其中(a)为参考ct图像，(b)为浮动mr1图像，(c)为浮动mr2图像，(d)为gc法(mr1-ct)，(e)为gc法(mr1-ct)图像差，(f)为lp法(mr1-ct)，(g)为lp法(mr1-ct)图像差，(h)为gc法(mr2-ct)，(i)为gc法(mr2-ct)图像差，(j)为lp法(mr2-ct)，(k)为lp法(mr2-ct)图像差。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在介绍本发明方法的具体内容前，首先对本发明中使用的一些基本知识做简单说明：

(1)zernike矩原理

zernike矩利用基函数vnm(ρ,θ)表示单位圆内的完全正交基，其定义为:

vnm(ρ,θ)＝rnm(ρ)e^jmθ(1)

其中：n表示阶数，m表示重数，0≤|m|≤n，n-|m|为偶数，ρ和θ分别表示极坐标下像素点的半径和角度，并且zernike距径向多项式如式(2)所示：

其中:

在p×p大小的图像块中，对于连续函数f(x,y)，将zernike矩离散归一化表示为：

其中*表示复共轭，λp表示离散归一化到单位圆后圆内像素点的个数，0≤ρxy≤1。对于旋转α角度后的浮动图像，由式(3)可推出和分别表示参考图像和浮动图像的zernike矩。等式两边同时取模，可得则表示目标图像旋转前后其幅值保持不变，所以zernike矩具有旋转不变性。

(2)gc法

gc法是解决能量最小化问题的一个较好的组合优化工具，在特定条件下,它会产生全局或局部最小值。gc法通过以下形式最小化能量函数来解决离散标签问题：

其中n是邻域系统，f是标签函数。为数据项，用来反映图像配准的精确程度。为平滑项，用来惩罚相邻像素间分配标签的差异性。

(3)mrf配准框架

i和j分别表示维度为d的参考图像和浮动图像，x表示图像的连续空间域，对于空间上任意一点x＝(x1,x2,...,xd)∈x，则i(x)和j(x)分别表示参考和浮动图像在x处的灰度值，d为变形矢量场，则

d^*表示配准成功时的最佳位移矢量场，λ是调节参数。ronald等把绝对值差的积分作为相似性测度c，变形矢量场d的一阶导数项作为平滑函数s，则有：

(4)基于zernike矩计算图像自相似性

zmld基本思想是借鉴了基于非局部均值的图像去噪算法中自相似性的概念。本发明根据预定义邻域中图像块的zernike矩来计算图像自相似性。

首先，用anm来表示znm的大小，则有

不同阶数和重数的zernike矩的大小代表不同的图像信息，如图1所示，(a)显示了brainweb数据库中mr-t1图像的a00、a11、a22的图像信息，由图(b)可以看出，a00非常接近原始图像的强度分布，结合式(7)可知a00可以表示图像的强度信息。由式(8)-(9)中a11和a22可以表示图像中局部结构产生的高频信息，如边缘和复杂纹理区域。但阶数越高，对图像噪声变得越敏感，如图(d)所示。基于对噪声鲁棒性和计算复杂度的综合考虑，用a00与a11来计算zmld，可以同时提取图像的强度信息和边缘、纹理特征。

参照图2，本发明实施例提供了基于zmld与gc离散优化的非刚性多模态医学图像配准方法，该方法包括以下步骤：

步骤一，读取输入的待配准的参考图像i和浮动图像j，两幅图像的分辨率相同；

步骤二，分别计算图像i和j两幅图像基于zernike矩的局部描述符zmld。将图像i和j在搜索邻域r中的图像分别划分为5×5的图像块，用anm表示不同的图像信息，a00和a11分别表示图像强度和边缘纹理特征，使用a00和a11计算zmld可以同时提取图像混合信息，图像i和图像j的zmld计算分别为：

其中，zmld(i,x,r)和zmld(j,x,r)分别表示图像i和j的zmld，a为正则化常数，表示图像i中以x为中心点和搜索邻域r中以r为中心点的图像块之间的距离，vnm(i,x)表示图像i中像素点x处的局部方差估计，r满足要求式(7)和(8)中和vnm(i,x)，(n＝m＝1或n＝m＝0)的计算公式分别如下：

和vnm(j,x)的计算和式(9)、(10)相同，p为以x为中心的图像块，式(9)需要对图像i中的所有像素点x和搜索位置r进行计算。由公式(7)和(9)得出，当搜索邻域r是一个自相似性较高的同质区域时，r中的图像块将具有相同的a00和a11，此时zmld(i,x,r)将取最大值。当邻域r是自相似性较低的边缘、纹理区域时，此时zmld(i,x,r)将小于最大值。由此可得，当处于不同的搜索邻域时，zmld将随着邻域的改变而发生变化，从而有效的表示了局部图像的自相似性。

v(i,x)是局部方差的估计，v较小会产生明显的衰减函数，较大则表示广泛的响应，所以v与图像中的噪声量有关。由式(7)和(10)可得，对于以x为中心的图像块p来说，使用zmld可以捕获与p相似的图像块，并且产生高响应；对不相似的图像块产生低响应，由此具有对图像噪声更好的鲁棒性。

距离有效表示了局部图像的自相似性，vnm则表示zmld对图像噪声具有良好的鲁棒性。所以，基于zmld在有噪声和强度失真的情况下，可以准确提取图像的结构信息。

步骤三，构造能量函数，并对能量函数进行离散化。本发明将图像配准问题看做是mrf的标签问题，给浮动图像j中每个图像块的中心点分配位移标签，标签即位移矢量，用来判断参考图像和浮动图像在空间上的位置是否达到一致，所有位移矢量的集合组成位移矢量场d。最佳标签集即为最优的配准效果，使能量函数最小的标签集就是最佳标签集。能量函数分为两部分：

ef＝edata+esmmoth

本发明用i和j的zmld之间的绝对误差和sad作为能量函数的数据项，来判断两幅图像的相似度；采用位移矢量场的一阶导数作为平滑项，用来惩罚相邻像素间有较大变化的位移标签。

i和j的zmld之间的sad为：

由式(11)可得，当i和j具有相同的zmld时，sadzmld将达到最小，此时两幅图像的相似度最高。

然后将能量函数离散化，离散能量方程如下所示：

其中：

||·||是l2范数，x1和x2分别表示两个相邻图像块的中心像素点。

步骤四，利用图割法gc的α扩展优化算法求解离散化后能量函数的最小值。为每个像素分配标签，通过α扩展形成网络图，图中结点表示图像块的中心点，源点以及汇点，边表示标签转换过程中所需能量。通过寻找网络图中的最小割即可得到能量函数的最小值。

α扩展移动的方法是：对于任意x∈x，x的标签变化f^*有两种可能：保持当前标签fx或者改变为α，所以可以把α扩展移动算法看做是一个两标签问题。标签为0则说明标签为1则说明根据三角不等式，可以证明：

由不等式(13)可得，任意相邻图像块的中心点像素x1，x2处的变形矢量场通过α扩展可达到全局最优。boykov等进一步证明在这种条件下，α扩展算法将最终收敛至局部最小值。

当通过每步α扩展形成图时，把浮动图像j中每个图像块的中心点作为图中的结点x，除了相邻结点相互连接外，每个结点还将增加两条边连接到源点和汇点，其中源点s和汇点t分别表示当前分配的标签f和α标签。两个结点之间的边的权值是平滑度，结点和源点之间的边的权值是当前标签f所需的能量，而将结点连接到汇点的边的权值是改变为α标签所需的能量。图的一维图像如图3所示。

在每个α扩展步骤中，图的最小割是mrf框架的最小能量。通过对每个标签进行α扩展，可以保证达到局部最小值，最终得到最佳位移矢量场d^*。gc的α扩展移动算法如下所示：

输入：任意位移矢量场d和能量函数ef,位移标签集l

1:将任意位移矢量场d作为初始的最佳位移矢量场d^*；

2:α标签遍历l中标签；

3:通过gc中d和α来最小化ef(d)；

4:如果ef(d)<ef(d^*)，则此时的d即为最佳矢量场d^*；否则返回第2步；

5:输出最佳位移矢量场d^*。

步骤五，输出最佳位移矢量场，即配准后图像。

实验说明

为验证本发明方法的有效性，分别从两个方面予以验证：

(1)基于结构表示的相似性测度对配准图像的影响

图4给出了rire数据库中一组脑部mr的t1和t2、pd加权图像，显示了在mrf配准框架下分别使用essd，mind和zmld三种不同的基于结构表示的相似性测度，并且使用gc优化的配准结果。为了得出实验结论，将mr的t1加权图像作为参考图像，对t2和pd图像施加人工形变，作为浮动图像，本组实验的图像大小都为256×256。图6(d)-(i)分别显示了t2-t1、pd-t1分别使用essd，mind和zmld的配准效果。

表1列出了分别使用三种方法所得tre的mean和std：

表1使用不同相似性测度配准的tre

相比图4所示的配准结果，mind和zmld要比essd得到的结果更接近参考图像。图中红色箭头所指的边缘和纹理区域中，图像块之间存在旋转。由于zmld具有旋转不变性，则zmld比mind能更有效的提取图像特征。mr的t1、t2加权图像具有相似的结构特征，所以表2中t2-t1的tre值要比pd-t1的值要小。所有这些方法的tre结果表明，在mrf离散的配准框架下，zmld的配准精度更高。综上所述，与essd和mind方法相比，zmld方法配准效果更好，精度更高。

图4(a)参考t1图像；(b)浮动t2图像；(c)浮动pd图像；(d)essd方法(t2-t1)；(e)mind方法(t2-t1)；(f)zmld方法(t2-t1)；(g)essd方法(pd-t1)；(h)mind方法(pd-t1)；(i)zmld方法(pd-t1)。

(2)优化算法对配准图像的影响

本组实验采用zmld来计算相似性测度，使用不同的优化算法来对比图像配准的精度和时间。图像配准的时间主要分为两部分，一部分是构造zmld来计算相似性测度所用的时间，另一部分是使用转换模型和优化算法实现配准的时间。本组实验的图像大小都为256×256。

(a)连续和离散优化算法对配准图像的影响

分别采用了基于b样条的自由形变模型(ffd)和限域拟牛顿算法(l-bfgs)的连续优化算法(简称ffd-lbfgs算法)进行配准及本发明采用mrf配准模型和gc的离散优化来实现配准过程。为验证本发明算法的优越性，图5给出了atlas数据集中的两组待配准图像，配准结果以及配准后图像差如图5的(e)-(l)所示。

图5(a)参考t2图像；(b)浮动ct图像；(c)参考spect图像；(d)浮动t2图像；(e)ffd-lbfgs算法(ct-t2)；(f)ffd-lbfgs算法(ct-t2)图像差；(g)本发明方法(ct-t2)；(h)本发明方法(ct-t2)图像差；(i)ffd-lbfgs算法(t2-spect)；(j)ffd-lbfgs算法(t2-spect)图像差；(k)本发明方法(t2-spect)；(l)本发明方法(t2-spect)图像差。

两种方法的配准后的tre和配准时间如表2所示:

表2不同优化算法下的tre和配准时间

由图5的(f)(h)、(j)(l)配准后的图像差以及表2中的tre的值可以得出，使用本发明的离散优化算法和ffd-lbfgs算法的连续优化得到图像差几乎一致，说明两种方法的配准精度相差不大。但对于配准时间来说，ffd-lbfgs算法中的总共配准运行时间大概在80秒左右，而本发明算法约60秒。这是因为连续优化算法中导数的计算量较大，并且易陷入局部最小值，而使用mrf的gc离散优化算法可以有效避免这些缺点，缩短配准时间。所以，本发明采用的离散优化算法有效提高了配准效率。

(b)不同离散优化算法对配准的影响

本组使用atlas数据集中急性中风的ct图像作为参考图像，把高频纹波变形的mr图像作为浮动图像1，把低频形变较大的mr图像作为浮动图像2。图6显示了使用基于zmld的局部描述符计算相似性测度，分别在gc和lp离散优化算法下的配准结果。图6的(d)-(k)分别显示了在gc和lp优化下的配准结果以及配准后的图像差。

图6(a)参考ct图像；(b)浮动mr1图像；(c)浮动mr2图像；(d)gc法(mr1-ct)；(e)gc法(mr1-ct)图像差；(f)lp法(mr1-ct)；(g)lp法(mr1-ct)图像差；(h)gc法(mr2-ct)；(i)gc法(mr2-ct)图像差；(j)lp法(mr2-ct)；(k)lp法(mr2-ct)图像差。

两种算法下配准图像的精度和运行时间如表3所示：

表3不同离散优化算法下配准图像的tre和时间

由图6的(e)和(g)的图像差可以看出，使用gc法优化的结果要比lp法得到的配准图像和参考图像差异性要小。对于浮动图像1来说，图像进行了复杂的波形变化，由于lp算法需要巨大的空间容量，导致lp不能使用大量的标签，所以会使lp无法对复杂变形的浮动图像1进行精确配准。(i)和(k)的两幅图的图相差几乎一致。浮动图像2虽然产生了大形变，但是由于lp和gc都可以在全局中搜索最小值，不易陷入局部最小值，所以两种方法都可以达到较好的配准效果。对于配准时间来说，lp和gc方法的计算时间都是一分钟左右。所以对比不同的离散优化方法，gc法具有更高的配准精度。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。