一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法技术

本发明专利技术涉及一种基于B样条和Levenberg‑Marquardt优化的快速图像配准方法，本发明专利技术将图像配准问题转化为目标函数最小值求解问题，采用三次B样条来模拟图像中的非刚性形变，采用Levenberg‑Marquardt优化方法对问题进行迭代求解，并在Levenberg‑Marquardt方法的基础上，重复利用计算量最大的Jacobian矩阵，引入近似迭代步长，并通过线性搜索寻找近似迭代步长的最佳权重。为了解决迭代过程中正则项系数对收敛速度的影响，本发明专利技术还引入了一种自适应的方法来确定系数，通过建立系数的变化与目标函数值的变化之间的相关关系，可以有效加速算法的收敛。

A fast image registration method based on B spline and Levenberg-Marquardt optimization

The invention relates to a fast image registration method based on B spline and Levenberg Marquardt optimization. The invention transforms the image registration problem into the minimum value solving problem of the target function. The three B spline is used to simulate the non rigid deformation in the image, and the Levenberg Marquardt optimization method is used to solve the problem iteratively. On the basis of the Levenberg Marquardt method, the maximum Jacobian matrix is reused, the approximate iterative step is introduced, and the optimal weight of the approximate iteration step is searched by linear search. In order to solve the influence of the regular term coefficient on the convergence speed in the iterative process, an adaptive method is introduced to determine the coefficient. By establishing the correlation between the change of the coefficient and the change of the target function, the convergence of the algorithm can be speeded up effectively.

【技术实现步骤摘要】
一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法
本专利技术涉及医学图像分析与处理
，具体涉及一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法。
技术介绍
图像配准是医学图像处理中最为关键的技术之一。借助医学图像配准技术，可以将源于不同的成像设备的两种或多种图像结合起来，取其各自的信息优势，在一幅图像上表达更加丰富的信息。同时，在临床上有时需要将患者的图像与正常人相同部分的图像对比，对辅助医生的诊断和治疗具有重要意义，图像配准也是该类应用的关键技术。近几十年以来，随着计算机及医疗设备的快速发展，图像配准技术还应用在图像引导手术、图像引导放疗、器官运动模型的建立等方面。随着计算机技术的不断发展、普及和应用，医学图像处理技术也有了飞速的发展，并且获得的图像数据也从二维扩展到三维甚至是四维，数据量随之产生爆炸性的增长。图像配准技术得到了更为广泛的应用，找到一种快速的医学图像配准算法也成为越来越关键的研究问题。图像配准分为刚性配准和弹性配准。在医学图像中，由于刚性形变并不能很好的模拟脏器的不规则运动和形变，因此往往采用弹性配准。基于B样条的弹性配准是应用最为广泛的一种配准方法。相比于其他方法，B样条弹性配准具有局部特性，即改变某一控制点只会影响该控制点周围邻域的像素值，因此节省了大量的计算，使得快速配准成为了可能。然而，随着医学影像数据量的迅速增长，现有的配准方法的速度与精度已经不能满足需求。
技术实现思路
本专利技术的目的在于提供一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法，用以解决现有配准方法速度慢、精度低的问题，采用基于B样条的自由形变模型作为图像的形变模型，并利用快速Levenberg-Marquardt优化方法对自由形变模型中的控制点最佳位置进行求解，该方法能够适应数据量较大的二维医学图像甚至三维医学图像。该基于B样条和Levenberg-Marquardt优化的快速图像配准方法包括以下步骤：(1)读入浮动图像T和参考图像R，在浮动图像R上增加控制点网格，并定义控制点网格的间距，控制点网格中x方向相邻两个控制点的间距为δx，控制点网格中y方向相邻两个控制点的间距为δy，此时任意一个控制点的位移向量均为0向量；(2)利用B样条基函数，根据控制点网格中控制点的位置计算得到浮动图像中每个像素点的位移向量，根据位移向量计算得到变形后的浮动图像T'；(3)用Levenberg-Marquardt优化方法对变形后的浮动图像T'的控制点网格位移向量进行更新；(4)重复上述步骤(2)、(3)，直至达到迭代终止条件；迭代终止条件为设定相邻两次迭代相似性测度变化的阈值，若相邻两次迭代相似性测度变化小于该阈值，则判定迭代终止。步骤(4)中相邻两次迭代的相似性测度为灰度均方差；所述步骤(2)具体包括以下步骤：201)浮动图像T中任意像素点p(x,y)位移向量的计算：像素p(x,y)位置的位移向量d(p,Sk)计算方法如下：其中，Sk为整个控制点网格第k次迭代的位移向量，为任意一个控制点的位移向量；由于B样条具有良好的紧支性和模拟非刚性形变的能力，因此采用三次B样条来模拟图像中的非刚性形变，Bl为三次B样条基函数，表示如下：202)计算变形后的浮动图像T'：变形后的浮动图像T'上像素点p(x,y)的灰度值由以下公式计算得到：T'(p)＝T(p-d(p,Sk))。图像配准的目的是对浮动图像施加一定的形变，使得变形后的浮动图像T'与参考图像尽可能地接近。描述相邻两次迭代的两幅图像的相似性可以有多种方法，以灰度均方差测度为例，设F＝T'-R为变形后的浮动图像与参考图像的残差向量，则F＝0时达到最佳配准，但实际上往往无法直接得到上述方程的解，因此利用Levenberg-Marquardt优化方法，转而求解min||F||2，也即将配准问题转化为优化问题。所述求解min||F||2所采用的迭代公式如下：Sk+1＝Sk+Dk+αDk'其中Sk+1为整个控制点网格第k+1次迭代的位移向量，Dk为迭代步长，Dk'为补充步长；迭代步长Dk与补充步长Dk'计算方法如下：通过对F进行Taylor展开，并引入正则项λI，λ为正则项系数，计算公式如下：(J(Sk)TJ(Sk)+λI)Dk＝-J(Sk)TF(Sk)其中，J(Sk)为在第k迭代控制点网格位移向量为Sk时残差向量F的Jacobian矩阵，J(Sk)T为J(Sk)的转置矩阵，迭代步长Dk如下：Dk＝-(J(Sk)TJ(Sk)+λI)-1J(Sk)TF(Sk)中间变量Sk'＝Sk+Dk，变形后的浮动图像T'与参考图像R的残差向量F(Sk')，补充步长Dk'如下：Dk'＝-(J(Sk)TJ(Sk)+λI)-1J(Sk)TF(Sk')引入线性搜索，在指定区间内寻找参数α，使得更新后相似性测度下降最大，也即max||F(Sk)||2-||F(Sk+1)||2由于迭代过程中的Jacobian矩阵非常庞大，为了减少计算量，在原Levenberg-Marquardt优化方法的基础上，本专利技术引入了补充步长Dk'，这样就可以重复利用计算量庞大的Jacobian矩阵，也正是为了计算补充步长Dk'，本专利技术引入了中间变量Sk'。由于迭代公式(J(Sk)TJ(Sk)+λI)Dk＝-J(Sk)TF(Sk)中的正则项系数λ对算法收敛速度具有显著影响，因此，正则项系数λ采用以下自适应更新方法：其中参数μ根据置信域半径rk的大小进行更新，更新公式如下：其中：0＜p1＜p2＜1，置信域半径rk的计算公式如下：本专利技术具有如下优点：本专利技术采用三次B样条来模拟图像中的非刚性形变，采用Levenberg-Marquardt优化方法对问题进行迭代求解，能够适应数据量较大的二维医学图像甚至三维医学图像。本专利技术引入近似迭代步长，在Levenberg-Marquardt方法的基础上，重复利用计算量最大的Jacobian矩阵，可以在保证配准精度的基础上有效提高算法的运行速度。本专利技术还通过线性搜索寻找近似迭代步长的最佳权重，为了解决迭代过程中正则项系数对收敛速度的影响，本专利技术引入了一种自适应的方法来确定系数，建立了系数的变化与目标函数值的变化之间的相关关系。附图说明图1为本专利技术配准方法的流程图。图2为基于三次B样条计算变形图像的示意图。图3为参考图像。图4为浮动图像。图5为配准前浮动图像与参考图像的残差图像。图6为采用本专利技术方法变形后的浮动图像。图7为采用本专利技术方法配准后的残差图像。图8为采用本专利技术方法最终匹配时控制点网格的状态。图9为本专利技术方法同其他方法在迭代次数上的对比结果。图10为本专利技术方法同其他方法在配准效果上的对比。具体实施方式以下实施例用于说明本专利技术，但不用来限制本专利技术的范围。实施例1参见图1～10，本专利技术基于B样条和快速Levenberg-Marquardt优化的快速图像配准方法，采用基于B样条的自由形变模型作为图像的形变模型，并利用快速Levenberg-Marquardt优化方法对自由形变模型中的控制点最佳位置进行求解，该方法能够适应数据量较大的二维医学图像甚至三维医学图像。所述配准方法涉及到的变量命名如下：R←参考图像，T←浮动图像本文档来自技高网...

【技术保护点】
1.一种基于B样条和快速Levenberg‑Marquardt优化的快速图像配准方法，其特征在于，所述配准方法包括以下步骤：(1)读入浮动图像T和参考图像R，在浮动图像R上增加控制点网格，并定义控制点网格的间距，控制点网格中x方向相邻两个控制点的间距为δx，控制点网格中y方向相邻两个控制点的间距为δy，此时任意一个控制点的位移向量均为0向量；(2)利用B样条基函数，根据控制点网格中控制点的位置计算得到浮动图像中每个像素点的位移向量，根据位移向量计算得到变形后的浮动图像T'；(3)用Levenberg‑Marquardt优化方法对变形后的浮动图像T'的控制点网格位移向量进行更新；(4)重复上述步骤(2)、(3)，直至达到迭代终止条件；迭代终止条件为设定相邻两次迭代相似性测度变化的阈值，若相邻两次迭代相似性测度变化小于该阈值，则判定迭代终止。

【技术特征摘要】
1.一种基于B样条和快速Levenberg-Marquardt优化的快速图像配准方法，其特征在于，所述配准方法包括以下步骤：(1)读入浮动图像T和参考图像R，在浮动图像R上增加控制点网格，并定义控制点网格的间距，控制点网格中x方向相邻两个控制点的间距为δx，控制点网格中y方向相邻两个控制点的间距为δy，此时任意一个控制点的位移向量均为0向量；(2)利用B样条基函数，根据控制点网格中控制点的位置计算得到浮动图像中每个像素点的位移向量，根据位移向量计算得到变形后的浮动图像T'；(3)用Levenberg-Marquardt优化方法对变形后的浮动图像T'的控制点网格位移向量进行更新；(4)重复上述步骤(2)、(3)，直至达到迭代终止条件；迭代终止条件为设定相邻两次迭代相似性测度变化的阈值，若相邻两次迭代相似性测度变化小于该阈值，则判定迭代终止。2.根据权利要求1所述的基于B样条和Levenberg-Marquardt优化的快速图像配准方法，其特征在于：步骤(4)中相邻两次迭代的相似性测度为灰度均方差；所述步骤(2)具体包括以下步骤：201)浮动图像T中任意像素点p(x,y)位移向量的计算：像素p(x,y)位置的位移向量d(p,Sk)计算方法如下：其中，Sk为整个控制点网格第k次迭代的位移向量，为任意一个控制点的位移向量；Bl为三次B样条基函数，表示如下：202)计算变形后的浮动图像T'：变形后的浮动图像T'上像素点p(x,y)的灰度值由以下公式计算得到：T'(p)＝T(p-d(p,Sk))。3.根据权利要求2所述的基于B样条和Levenberg-Marquardt优化的快速图像配准方法，其特征...

【专利技术属性】
技术研发人员：吕科，董继阳，
申请(专利权)人：北京中科嘉宁科技有限公司，
类型：发明
国别省市：北京,11