CN110517300A - 基于局部结构算子的弹性图像配准算法 - Google Patents

Abstract

一种基于局部结构算子的弹性图像配准算法，通过分析局部边界熵算子存在的不足，对其进行有效改进从而得到两个新的特征描述算子，即局部结构熵和局部变异系数；然后将上述特征算子和图像梯度方向信息整合起来，设计一个新的相似性测度执行图像的弹性配准，能够较为准确地获取所需的配准性能，为医学图像的分析与理解，多模态信息的融合提供有效的辅助。

Description

基于局部结构算子的弹性图像配准算法

技术领域

本发明具体涉及图像处理技术领域，特别是针对图像背景复杂、组织对比度弱、灰度均匀性差的图像间的弹性配准，可用于临床医学图像中兴趣目标的空间定位和信息融合中。

背景技术

图像配准是一个以空间变换参数为参数，以相似性测度为目标函数的优化问题，通过寻找合适的空间变换，使待配准的两幅图像之间信息差异最小化，从而实现图像信息的整合。这种配准对图像理解与分析，多模态信息融合，目标定位等领域具有重要的价值而得到广泛的研究。目前，已有的图像配准，尤其是多模态弹性配准，通常可分为基于信息论的配准和基于特征信息的配准。基于信息论的配准方法主要通过开发具有模态独立特性的相似性测度，实现图像的配准。这类方法存在以下不足：

(1)灰度概率分布的评估不够准确；

(2)概率分布的计算中通常会忽略重要的局部纹理信息，从而使得这类方法在某些局部区域具有较大的配准误差。基于特征信息的配准方法利用预先定义的特征描述算子实现图像模态独立和特征纹理的探测。这类方法能够使用简单的相似性测度完成复杂的图像配准任务而成为一个重要的研究方向，而如何设计具有独特探测能力的局部结构描述算子就成为研究的关键。这不仅关系到图像中不同特征的提取，而且决定了配准的精确性和时效性。

为了提取图像中各种不同的纹理特征，就需要对局部图像区域内相邻像素间的灰度和位置信息进行对比分析，寻找他们间潜在的信息差异，并将这种差异转化为一个与局部灰度信息对应的特征物理量。这种特征量的分布要在一个相对大的取值范围内，以便凸显局部区域内细微的灰度波动，实现弱组织对比度和灰度不均匀的抑制。现有特征算子中，局部边界熵和图像梯度能够在一定程度上实现图像特征信息的探测，但是它们都存在明显的不足。具体表现为，局部边界熵由于取模运算符的存在而无法分辨灰度值小于局部均值的像素，而图像梯度对噪声较为明显且易受灰度不均匀的影响。因此，合理改进这些特征描述算子是图像配准的关键。

除了开发合适的特征描述算子，还需要设计一个能够有效量化特征描述算子之间信息差异的相似性测度，用于评估待配准图像间的不一致性(dissimilarities)。这种不一致性刻画的越全面准确，配准算法就会具有越高的准确性，因此相似性测度的设计在一定程度上关系到配准算法的性能。基于特征的配准方法中，较为常用的相似性测度为差的平方和(sum of squared differences,SSD)。这种测度可以作用在像素灰度上实现单模态图像的配准，也可以作用在特征描述算子上执行多模态的配准任务。此外，该测度的计算较为简洁，能够通过并行处理实现快速的配准操作而得到广泛的应用。但是该测度也存在一个问题：当特征描述算子的取值范围相差很大的时候，会导致图像信息量化存在一定的偏差而导致较低的配准性能。

发明内容

为了解决现有技术中针对局部边界熵算子存在的不足，本发明提出新的特征描述算子降低弱组织对比度和灰度不均匀现象对配准的干扰，并基于这些特征算子实施图像的弹性配准，实现兼顾全局和局部的高精度配准。

本发明的主要计算思路为：通过分析局部边界熵算子存在的不足，对其进行有效改进从而得到两个新的特征描述算子(即局部结构熵和局部变异系数)；然后将上述特征算子和图像梯度方向信息整合起来，设计一个新的相似性测度执行弹性图像的配准。具体步骤如下：一种基于局部结构算子的弹性图像配准算法，包括以下步骤：

(1)局部结构算子的开发：将像素灰度及其对应的取模操作结果整合到局部边界熵的计算中，从而引出两种新的局部结构描述算子，即局部结构熵和局部变异系数，实现局部区域内不同图像特征的提取；

(2)相似性测度的设计：使用局部边界熵、局部结构熵、以及局部变异系数对图像进行处理后可以得到三幅不同的特征图，将它们组合起来构成一幅图像对应的特征矢量图，即完成灰度图像到特征矢量图的转换，将待配准的参考图像和浮动图像通过上述方法分别转换为对应的特征矢量图，即参考矢量图和浮动矢量图，计算待配准图像对应特征矢量图之间差的平方和，辅助待配准图像间信息差异的评估，将待配准图像梯度方向角之间的差异整合到差的平方和计算中，从而得到一种新的相似性测度，该测度能同时从不同的角度准确量化待配准图像间的信息差异，极大地降低弱组织对比度和灰度不均匀现象在配准中的影响。

(3)自由形变配准框架：使用相似性测度量化待配准图像间的信息差异并将这种差异最小化，便需要通过合适的空间变换将浮动图像变换到参考图像空间使二者尽可能地重叠起来，在弹性图像配准中，常使用的空间变换是基于三次样条函数的自由形变变换，在该变换下，浮动图像中每个像素点的位置变动通过样条函数计算得到，由于样条函数具有二阶导连续的特性，这使得空间变换能够有效模拟全局与局部区域内像素的位置变动，浮动图像在像素位置变动后通过局部结构描述算子更新特征矢量图和相似性测度值，利用优化算法将相似性测度值最小化处理，即可得到所需的配准结果。

所述的局部边界熵可表示为：

其中，Π_X是一个中心在X处边长为r的三维正方形区域，该区域内的灰度平均值为m，I(y)为坐标位置y处的像素灰度，p(y)为灰度I(y)对应的概率分布，log(·)为自然对数运算符，mod(·)为取模运算符。取模运算符的存在使得概率分布具有如下特性：

其中，K为不小于0的正整数，它能够在一定程度上反映局部区域内像素灰度相对于灰度均值的波动程度性；

此外，取模运算符无法有效探测灰度值小于局部均值(即I<m)的像素，导致局部边界熵算子在某些情况下具有较低的边界探测能力,为了缓解取模运算符存在的不足，我们对p(y)进行改进得到一种新的概率p^*(y)，具体表示如下：

基于这个新的概率分布可以得到一个局部结构熵，公式表示为：

此外，概率p(y)之和也可以作为一个重要的图像纹理特征反映局部区域内像素灰度相对于灰度均值的波动程度，因此将其作为局部变异系数，用于图像配准中。该变异系数公式表示为：

根据上面提到的局部结构描述算子LEE(X)，LSE(X)和LVF(X)可以看出：同一个像素位置X对应三种不同的特征纹理信息，利用这些算子可以将待配准的参考图像和浮动图像分别转化为参考特征图和浮动特征图。

所述的步骤(2)相似性测度的设计中所述的相似性测度构建步骤如下：

(2a)利用三种局部结构描述算子LEE(X)，LSE(X)和LVF(X)，将待配准图像转化为特征矢量图，即V(X)＝(LEE(X),LSE(X),LVF(X))，从而导致两个不同的特征矢量图，分别表示为参考特征矢量图V_R和浮动特征矢量图V_F；计算这两个特征矢量图之间差的平方和，可辅助待配准图像之间的信息差异的量化，基于参考与浮动特征矢量图对应的差的平方和可公式表示为：

其中，Ω为全局图像区域，T(X,Φ)为一个以Φ为参数的空间变换；(2b)为进一步量化待配准图像间的信息差异，我们计算图像对应灰度梯度方向角之间的差异，并将这种差异整合到差的平方和计算中，得到所需的相似性测度，具体表示为：

在参考图像中坐标X处和浮动图像中坐标Y处对应的梯度方向角计算公式为：

其中，arccos(·)为反余弦函数，为梯度运算符，I_R和I_F分别表示参考图像和浮动图像的灰度，|·|表示梯度模值运算符。

所述的步骤(3)自由形变配准框架中空间变换的求解步骤如下：

(3a)在自由形变配准框架下，空间变换参数Φ表示一个N_x×N_y×N_z的网格，在每个坐标轴方向上的网格间距(spacing)为s＝(s_x,s_y,s_z)，每个网格结点在坐标(i,j,k)处的位置偏移为d_i,j,k，对浮动图像中任一像素点X＝(x,y,z)，其对应的空间变换T(X,Φ)为：

其中，u＝x/s_x-i-1,v＝y/s_y-i-1,w＝z/s_z-i-1,为向下取整运算符，L_l表示样条函数中第l个基函数，它可以表示为L₀(t)＝(1-t)³/6，L₁(t)＝(3t³-6t²+4)/6，L₂(t)＝(-3t³+3t²+3t+1)/6和L₃(t)＝t³/6，其中0≤t<1；

(3b)在图像弹性配准中，为了保证变形场的平滑且符合图像目标的实际形变，我们需要对变形场进行约束处理，约束处理的计算公式为：

其中，B(·)为形变约束项，N为参与运算的像素数目。将这个约束项整合到相似性测度中，即可得到最终的配准代价函数为C＝(1-λ)S(Φ)+λB(Φ)，其中λ为权重因子。

本发明的有益效果是：本发明提供了一种基于局部结构算子的弹性图像配准算法，通过分析局部边界熵算子存在的不足，对其进行有效改进从而得到两个新的特征描述算子，即局部结构熵和局部变异系数；然后将上述特征算子和图像梯度方向信息整合起来，设计一个新的相似性测度执行弹性图像的配准，能够较为准确的弹性配准性能，为医学图像的分析与理解，多模态信息的融合提供有效的辅助。

附图说明

图1是本发明的流程图。

图2是本发明针对公开的人工合成磁共振图像进行的仿真实验，第一列为参考与浮动图像，第二列为参考与浮动手工标注图，第三列为配准前后手工标注图之间的差图。

图3是本发明针对临床磁共振图像进行的仿真实验，第一列为参考与浮动图像，第二列为参考与浮动手工标注图，第三列为配准前后手工标注图之间的差图。

图4是本发明针对临床磁共振图像进行的仿真实验，第一列为参考与浮动图像，第二列为参考与浮动手工标注图，第三列为配准前后手工标注图之间的差图。

图5是本发明算法与现有的几种图像配准算法的对比结果。

具体实施方式

本发明的总体思路为：(1)局部边界熵算子对应的概率分布计算中，取模运算符虽然能够在一定程度上探测局部区域内灰度的相对差异，但是它无法探测灰度值小于局部均值的像素。为探测灰度值小于局部均值的像素，可将像素灰度、局部灰度均值、以及灰度取模操作结果三个不同的物理量整合起来，得到一种新的灰度特征。最后，用取模运算符对这种灰度特征进行处理，即可有效探测局部区域内相邻像素间的灰度波动。

(2)局部边界熵和它的两个改进算子在特征空间的取值范围是不同的，直接将它们用于差的平方和测度中可能会导致一定程度的计算偏差。为此，我们将待配准图像对应的梯度方向信息整合到由特征描述算子构建相似性测度中，用以降低特征算子取值范围不同造成的影响，提升待配准图像间信息差异的量化效果。

(3)上述的特征描述算子和相似性测度可用于刚性配准和弹性配准中，用于弹性配准时，便需要一个能够模拟图像局部变形的空间变换。这种变换通常基于三次样条函数构建而成，被称为自由形变(free form deformation,FFD)。在该形变框架下，每个像素点的位置偏移量通过样条函数的插值计算得到。

下面结合附图对基于局部结构算子的弹性图像配准算法做进一步描述；

参考图1，本发明基于局部结构算子的弹性图像配准算法，包括如下步骤：

步骤1，构建能够有效探测细微灰度波动的局部结构算子

(1a)根据局部边界熵算子对应概率分布计算存在的不足，在图像中取一个大小为r×r×r的矩形区域，将区域内的像素灰度正则化处理，通过公式和p(y)＝mod(I(y),m)/m，以及p^*(y)＝p(y)/∑p(y)计算两种不同的概率分布，将它们应用到图像熵的计算公式中得到两种不同的熵信息(即局部边界熵和局部结构熵)。此外，由于取模运算符(modulus operator)mod(·)具有数值选择性截断的特性，上述概率分布使用这种运算符会存在一部分灰度被截断而另一部分灰度被保存的现象，进而导致取模操作处理后的局部灰度之和发生变化。这种变化与局部灰度分布密切相关，因此可以作为一种新的特征信息用于图像配准中，即

(1b)将上述三种不同的图像特征纹理整合起来，可以构成一个特征纹理矢量，它们从不同的方面刻画同一幅图像，因此能够有效降低弱组织对比度和灰度不均匀等现象的影响，并且比像素灰度更适于图像弹性配准任务中。

步骤2，基于局部结构算子的相似性测度

使用上述三种局部描述算子处理一副图像中将产生三幅不同的特征纹理图(feature map)，整合这些特征图可构成一个特征图矢量。基于此，可将待配准的参考与浮动图像转化为参考与浮动特征图矢量。利用差的平方和(sum of squared differences)测度可评估两组特征图矢量之间的信息差异，辅助原始图像的配准。

使用差的平方和测度计算特征图矢量之间的信息差异具有操作简单且可并行处理的特点，但该计算过程存在一定的偏差，因为三种特征算子的取值范围是不同的，它们在测度的计算中起的作用也是不一样的。为了降低计算偏差，需要引入图像梯度方向角信息对差的平方和测度进行约束，从而得到本发明所需的相似性测度。

步骤3，自由形变配准框架

图像配准过程中，浮动图像需要逐渐向参考图像靠拢，以实现它们间信息差异的最小化，即相似性测度的最小化。浮动图像向参考图像靠拢的过程会发生局部和全局区域的形变。这种形变可通过自由形变框架进行模拟，使用优化算法搜索自由形变框架中样条函数对应的最优空间变换参数即可完成图像的弹性配准。

1、仿真条件：

本发明在Windows 10 64bit Intel(R)Core(TM)i7-6700HQ CPU@2.60GHz2.60GHz RAM 16GB平台上的MATLAB 2013a软件上进行仿真模拟的，仿真数据选用公开的磁共振图像数据进行弹性配准实验，实验数据来源网址分别为https://brainweb.bic.mni.mcgill.ca/brainweb/和http://www.loni.ucla.edu/Atlases/LPBA40/。

2、仿真内容与结果

1)仿真实验1：

本仿真实验使用人工合成的多模态磁共振图像进行配准实验，验证算法的有效性，实验结果呈现在图2中：

从图像配准前后的对比中可以看出：配准前待配准的多模态两幅图像存在较大的灰度、位置以及形态学上的差异；但是配准后两幅图像间的不一致性明显降低。

图2中，第一列为参考和浮动磁共振图像，第二列分别为参考和浮动图像对应的手工标注图，分别标注出了脑部白质(white matter,WM)、灰质(gray matter,GM)以及脑脊液(cerebrospinal fluid,CSF)三种组织，第三列为配准前和配准后两幅手工标注图之间的差图(difference image)。从图中可以看出，本算法能够有效降低手工标注图像间的形态学差异。

2)仿真实验2：

本仿真实验使用相同模态的临床磁共振图像进行弹性配准实验，验证算法对临床实际图像的配准性能，配准结果呈现在图3中：

图3中，第一列为相同模态的参考和浮动磁共振图像，第二列分别为参考和浮动图像对应的手工标注图，标注兴趣区域分别为白质、灰质和脑脊液，第三列为配准前和配准后两幅手工标注图之间的差图。从图中可以看出，配准前两幅图像存在较大的位置误差，配准后图像间的位置和形态学差异显著降低了。

3)仿真实验3：

本仿真实验使用相同模态的临床磁共振图像进行弹性配准实验，验证算法对临床实际图像的配准性能，配准结果呈现在图4中：

图4中，第一列为相同模态的参考和浮动磁共振图像，第二列分别为参考和浮动图像中54个兴趣区域(region of interest,ROI)对应的手工标注，第三列为配准前和配准后两幅手工标注图之间的差图。从图中可以看出，配准前两幅图像存在较大的位置误差，配准后图像间的位置和形态学差异显著降低了。

3)仿真实验4：

在本仿真实验中，通过对比本算法与已有配准算法(即eSSD算法，NMI算法以及MIND算法)间的性能差异，验证本文发明的图像配准有效性。

对比实验使用相同模态的临床磁共振图像进行的配准性能对比，实验结果通过重叠精度(Dice Similarity Coefficient,DSC)量化三个兴趣区域(即脑白质、灰质、脑脊液)间的重合比例，该指标可表示如下：

其中，A表示参考图像的手工标注结果，B为算法配准后的结果，∩表示交集操作符，N(·)为给定区域内像素数目求解函数.DSC的取值范围在0到1之间，且DSC越大，表示算法配准越准确。

仿真实验的对比结果参考图5，从图中能够较为直接地观察到MIND算法在三个兴趣区域上都获得了最高的DSC值，因此具有最好的配准精度；本文的算法稍微弱于MIND算法，但是明显比eSSD和NMI算法具有更好的配准性能。

对比结果表明，本文发明的算法能够与现有的配准算法(如MIND)在配准性能上进行比较，而且优于某些配准算法(eSSD和NMI)。

各位技术人员须知：虽然本发明已按照上述具体实施方式做了描述，但是本发明的发明思想并不仅限于此发明，任何运用本发明思想的改装，都将纳入本专利专利权保护范围内。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种基于局部结构算子的弹性图像配准算法，其特征在于，包括以下步骤：

(1)局部结构算子的开发：将像素灰度及其对应的取模操作结果整合到局部边界熵的计算中，从而引出两种新的局部结构描述算子，即局部结构熵和局部变异系数，实现局部区域内不同图像特征的提取；

(2)相似性测度的设计：使用局部边界熵、局部结构熵、以及局部变异系数对图像进行处理后可以得到三幅不同的特征图，将它们组合起来构成一幅图像对应的特征矢量图，即完成灰度图像到特征矢量图的转换，将待配准的参考图像和浮动图像通过上述方法分别转换为对应的特征矢量图，即参考矢量图和浮动矢量图，计算待配准图像对应特征矢量图之间差的平方和，辅助待配准图像间信息差异的评估，将待配准图像梯度方向角之间的差异整合到差的平方和计算中，从而得到一种新的相似性测度，该测度能同时从不同的角度准确量化待配准图像间的信息差异，极大地降低弱组织对比度和灰度不均匀现象在配准中的影响。

(3)自由形变配准框架：使用相似性测度量化待配准图像间的信息差异并将这种差异最小化，便需要通过合适的空间变换将浮动图像变换到参考图像空间使二者尽可能地重叠起来，在弹性图像配准中，常使用的空间变换是基于三次样条函数的自由形变变换，在该变换下，浮动图像中每个像素点的位置变动通过样条函数计算得到，由于样条函数具有二阶导连续的特性，这使得空间变换能够有效模拟全局与局部区域内像素的位置变动，浮动图像在像素位置变动后通过局部结构描述算子更新特征矢量图和相似性测度值，利用优化算法将相似性测度值最小化处理，即可得到所需的配准结果。

2.根据权利要求1所述的基于局部结构算子的弹性图像配准算法，其特征在于，所述的局部边界熵可表示为：

其中，Π_X是一个中心在X处边长为r的三维正方形区域，该区域内的灰度平均值为m，I(y)为坐标位置y处的像素灰度，p(y)为灰度I(y)对应的概率分布，log(·)为自然对数运算符，mod(·)为取模运算符。取模运算符的存在使得概率分布具有如下特性：

其中，K为不小于0的正整数，它能够在一定程度上反映局部区域内像素灰度相对于灰度均值的波动程度性；

此外，取模运算符无法有效探测灰度值小于局部均值(即I<m)的像素，导致局部边界熵算子在某些情况下具有较低的边界探测能力,为了缓解取模运算符存在的不足，我们对p(y)进行改进得到一种新的概率p^*(y)，具体表示如下：

基于这个新的概率分布可以得到一个局部结构熵，公式表示为：

此外，概率p(y)之和也可以作为一个重要的图像纹理特征反映局部区域内像素灰度相对于灰度均值的波动程度，因此将其作为局部变异系数，用于图像配准中。该变异系数公式表示为：

根据上面提到的局部结构描述算子LEE(X)，LSE(X)和LVF(X)可以看出：同一个像素位置X对应三种不同的特征纹理信息，利用这些算子可以将待配准的参考图像和浮动图像分别转化为参考特征图和浮动特征图。

3.根据权利要求2所述的基于局部结构算子的弹性图像配准算法，其特征在于，所述的步骤(2)相似性测度的设计中所述的相似性测度构建步骤如下：

(2a)利用三种局部结构描述算子LEE(X)，LSE(X)和LVF(X)，将待配准图像转化为特征矢量图，即V(X)＝(LEE(X),LSE(X),LVF(X))，从而导致两个不同的特征矢量图，分别表示为参考特征矢量图V_R和浮动特征矢量图V_F；计算这两个特征矢量图之间差的平方和，可辅助待配准图像之间的信息差异的量化，基于参考与浮动特征矢量图对应的差的平方和可公式表示为：

其中，Ω为全局图像区域，T(X,Φ)为一个以Φ为参数的空间变换；

(2b)为进一步量化待配准图像间的信息差异，我们计算图像对应灰度梯度方向角之间的差异，并将这种差异整合到差的平方和计算中，得到所需的相似性测度，具体表示为：

在参考图像中坐标X处和浮动图像中坐标Y处对应的梯度方向角计算公式为：

其中，arccos(·)为反余弦函数，为梯度运算符，I_R和I_F分别表示参考图像和浮动图像的灰度，|·|表示梯度模值运算符。

4.根据权利要求1所述的基于局部结构算子的弹性图像配准算法，其特征在于，所述的步骤(3)自由形变配准框架中空间变换的求解步骤如下：

(3a)在自由形变配准框架下，空间变换参数Φ表示一个N_x×N_y×N_z的网格，在每个坐标轴方向上的网格间距(spacing)为s＝(s_x,s_y,s_z)，每个网格结点在坐标(i,j,k)处的位置偏移为d_i,j,k，对浮动图像中任一像素点X＝(x,y,z)，其对应的空间变换T(X,Φ)为：

其中，u＝x/s_x-i-1,v＝y/s_y-i-1,w＝z/s_z-i-1，为向下取整运算符，L_l表示样条函数中第l个基函数，它可以表示为L₀(t)＝(1-t)³/6，L₁(t)＝(3t³-6t²+4)/6，L₂(t)＝(-3t³+3t²+3t+1)/6和L₃(t)＝t³/6，其中0≤t<1；

(3b)在图像弹性配准中，为了保证变形场的平滑且符合图像目标的实际形变，我们需要对变形场进行约束处理，约束处理的计算公式为：

其中，B(·)为形变约束项，N为参与运算的像素数目。将这个约束项整合到相似性测度中，即可得到最终的配准代价函数为C＝(1-λ)S(Φ)+λB(Φ)，其中λ为权重因子。