CN103942752B - 快速一致性图像变换方法及变换系统 - Google Patents

Abstract

本发明涉及一种快速一致性图像变换方法及变换系统。所述方法包括如下步骤：设定一致性误差预设值及最大迭代次数，并根据源控制点与目标控制点的对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，据此计算相应的正向一致性误差及反向一致性误差，然后判断是否正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，如果是，则结束，否则根据所述正向一致性误差及反向一致性误差调整正向变换函数及反向变换函数，直到达到结束条件。本发明实现了快速图像一致性变换，可以解决图像小形变弹性配准中的一致性变换问题，有效提高求解一致性变换的效率。

Description

快速一致性图像变换方法及变换系统

技术领域

本发明属于图像处理领域，主要涉及一种快速一致性图像变换方法及变换系统。

背景技术

图像配准是图像处理、图像分析、图像融合和图像识别与检测等技术的重要步骤，其已经被广泛运用于医学图像处理、计算机视觉、模式识别等领域。基于控制点对的图像配准方法是图像配准的研究内容，而一致性变换函数的求解是其重要研究内容。基于控制点对的一致性图像配准方法同时求解正反形变函数，保持正向形变函数与反向形变函数的一致性，可以同时得到更准确的正向形变函数及反向形变函数，具有重要的应用价值。其中，H.J.Johnson和G.E.Christensen于2002年提出了一种针对小形变的一致性变换求解方法，通过迭代求解一致性变换，使得正反变换具有最小的扭曲能量和一致性误差。假设正向形变为h，形变位移场为u(x)；反向变换为g，形变位移场为w(x)，则h(x)＝x+u(x)，g(x)＝x+w(x)。定义正向变换的反函数为h^-1，其位移场为反向变换的反函数为g^-1，其位移场为则

假设控制点对应关系(p_i,q_i)已知，其中p_i是第i个源控制点，q_i是第i个目标控制点，H.J.Johnson和G.E.Christensen给出的求解过程如下：

步骤1，r_i＝p_i，s_i＝q_i；u(x)＝0；w(x)＝0，设定优化步长α和β，最大控制点偏移误差迭代次数iter，最大迭代次数miter等；

步骤2，基于控制点的对应关系利用薄板样条插值方法求解正向形变函数f₁(x)，满足f₁(r_i)＝q_i；基于控制点的对应关系利用薄板样条插值方法求解反向形变函数f₂(x)，满足f₂(s_i)＝p_i；

步骤3，u(x)＝u(x)+α[f₁(x)-x]，w(x)＝w(x)+α[f₂(x)-x]；

步骤4，r_i＝p_i+u(r_i)，s_i＝q_i+w(s_i)，iter＝iter+1；

步骤5，求取正向形变的逆函数h^-1(x)，反向形变的逆函数g^-1(x)；

步骤6，更新正、反形变的位移场。u(x)＝u(x)-β[u(x)-g^-1(x)+x]，w(x)＝w(x)-β[w(x)-h^-1(x)+x]；

步骤7，检查是否满足终止准则，用avgerr_q→p表示控制点偏移误差，如果iter＞miter，或者或者时，迭代结束；否则返回步骤2。

H.J.Johnson和G.E.Christensen提出的上述数值求解方法非常耗时，在每次迭代中，其步骤2中需要进行两次薄板样条插值运算，其时间复杂度为O(nN)(N为图像中所有像素点个数，n为控制点个数)；步骤3调整u(x)、w(x)的运算复杂度是O(N)；步骤5的求逆计算需要对每个像素点进行多次迭代，时间复杂度为O(NT)(T为求逆过程中的迭代次数)；步骤6对每个点进行一次误差计算，时间复杂度为O(N)。该算法迭代一共需要进行两次形变、四次调整(包括：两次调整正、反变换的位移场u(x)、w(x))、两次求逆、一次判定，每次迭代需要花费较长时间。而算法需要多次迭代才能收敛，其计算时间是非常惊人的。

发明内容

本发明所要解决的技术问题是：提供一种快速一致性图像变换方法及变换系统，旨在解决现有的基于控制点对应关系的小形变一致性变换方法中耗时较长的问题。本发明是这样实现的：

一种快速一致性图像变换方法，包括如下步骤：

步骤A，设定一致性误差预设值及最大迭代次数，并根据源控制点与目标控制点的对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，然后执行步骤B；

步骤B，根据所述正向变换函数及反向变换函数计算相应的正向一致性误差及反向一致性误差，然后执行步骤C；

步骤C，判断所述正向一致性误差及反向一致性误差，并统计迭代次数；如果正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，则结束；否则执行步骤D；

步骤D，根据所述正向一致性误差及反向一致性误差调整正向变换函数及反向变换函数，然后返回步骤B。

进一步地，所述步骤B具体包括如下步骤：

用户输入n对控制点对，利用薄板样条基函数构造变换函数，得到正向变换函数f₁(x)，以及反向变换函数f₂(x)；

所述正向一致性误差及反向一致性误差通过如下公式计算：

正向一致性误差：

反向一致性误差：

其中，n₁,n₂分别为正向及反向变换初始的均匀离散网格点坐标。

进一步地，所述步骤C包括如下步骤：

计算最大正向一致性误差max_δ₁及最大反向一致性误差max_δ₂：

并统计迭代次数；

判断是否max(max_δ₁,max_δ₂)≤ε，ε为设定的一致性误差预设值，或者迭代达到设定的最大迭代次数；如果是，则结束，否则执行步骤D。

进一步地，正向变换函数及反向变换函数按如下公式调整：

f₁(x)＝f₁(x)+α*δ₁；

f₂(x)＝f₂(x)+β*δ₂；

其中α、β分别是正向变换函数及反向变换函数调整的步长，第k次调整步长设置为：

α_k+1＝ηα_k

β_k+1＝ηβ_k

其中，η是预设的调整因子，0＜η＜1。

一种快速一致性图像变换系统，包括：函数构造模块、误差计算模块、判断模块、函数调整模块；其中：

函数构造模块用于设定一致性误差预设值及最大迭代次数，并根据源控制点与目标控制点对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，然后跳转至误差计算模块；

误差计算模块用于根据所述正向变换函数及反向变换函数计算相应的正向一致性误差及反向一致性误差，然后跳转至判断模块；

判断模块用于判断所述正向一致性误差及反向一致性误差，并统计迭代次数；如果正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，则结束；否则跳转至函数调整模块；

函数调整模块用于根据所述正向一致性误差及反向一致性误差调整所述正向变换函数及反向变换函数，然后跳转至误差计算模块。

进一步地，所述误差计算模块具体用于利用薄板样条基函数构造变换函数，得到正向变换函数f₁(x)，以及反向变换函数f₂(x)；

所述正向一致性误差及反向一致性误差通过如下公式计算：

正向一致性误差：

反向一致性误差：

其中，n₁,n₂分别为正向及反向变换初始的均匀离散网格点坐标。

进一步地，所述判断模块用于：

计算最大正向一致性误差max_δ₁及最大反向一致性误差max_δ₂：

并统计迭代次数；

判断是否max(max_δ₁,max_δ₂)≤ε，ε为设定的一致性误差预设值，或者迭代达到设定的最大迭代次数；如果是，则结束，否则跳转至函数调整模块。

进一步地，正向变换函数及反向变换函数按如下公式调整：

f₁(x)＝f₁(x)+α*δ₁；

f₂(x)＝f₂(x)+β*δ₂；

其中α、β分别是正向变换函数及反向变换函数调整的步长，第k次调整步长设置为：

α_k+1＝ηα_k

β_k+1＝ηβ_k

其中，η是预设的调整因子，0＜η＜1。

与现有技术相比，本发明出的快速一致性图像变换方法具有简单、快速的特点，算法运行一次相当于H.J.Johnson和G.E.Christensen算法的两次求逆的计算量，所需要花费的时间很少，有利于处理运算量较大的小形变一致性图像变换问题。

附图说明

图1：本发明实施例提供的快速一致性图像变换方法的流程示意图；

图2：本发明实施例提供的快速一致性图像变换系统的组成示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。

图1是本发明实施例提供的快速一致性图像变换方法的流程示意图，图2是本发明实施例提供的快速一致性图像变换系统的组成示意图。结合图1及图2，在快速一致性图像变换系统的基础上，对快速一致性图像变换方法各步骤详述如下：

在步骤S101中，函数构造模块1设定一致性误差预设值ε及最大迭代次数,根据源控制点与目标控制点的对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，然后执行步骤S102。正向变换函数和反向变换函数按如下方法构造：

首先，用户输入n对控制点对，其中，源控制点p_i及目标控制点q_i表示所述n对控制点对中的第i对控制点对，对用户输入的源控制点p_i，目标控制点q_i，利用薄板样条基函数构造变换函数，得到从源控制点映射到目标控制点的正向变换函数f₁(x)，以及从目标控制点映射到源控制点的反向变换函数f₂(x)。

在步骤S102中，误差计算模块2根据正向变换函数及反向变换函数求解相应的正向一致性误差及反向一致性误差，然后执行步骤S103。假设n₁,n₂分别为正向、反向变换初始的均匀离散网格点坐标，正向一致性误差及反向一致性误差通过如下公式计算：

正向一致性误差：

反向一致性误差：

在步骤S103中，判断模块3判断是否正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，如果是，则结束，否则执行步骤S104。具体方法如下：

判断模块3首先计算最大正向一致性误差max_δ₁及最大反向一致性误差max_δ₂；

并统计迭代次数；

然后判断是否max(max_δ₁,max_δ₂)≤ε(ε为设定的一致性误差预设值)，或者迭代达到设定的最大迭代次数，如果是，则结束，否则执行步骤S104。

在步骤S104中，函数调整模块4根据正向一致性误差及反向一致性误差调整正向变换函数及反向变换函数，然后返回步骤S102。正向变换函数及反向变换函数按如下公式调整：

f₁(x)＝f₁(x)+α*δ₁；

f₂(x)＝f₂(x)+β*δ₂；

其中α、β分别是正向变换函数及反向变换函数调整的步长，随着迭代次数的增加，正向一致性误差及反向一致性误差越来越小，需要调整的步长就越小。第k次调整步长设置为：

α_k+1＝ηα_k

β_k+1＝ηβ_k

其中，η为调整因子，0＜η＜1，可以取0.99。

对比本发明与H.J.Johnson和G.E.Christensen算法可以看出，本发明只有一次正向变换函数及反向变换函数求解过程，没有函数求逆过程，正向变换函数及反向变换函数的调整步骤也要少一半，从而减小了大量的计算时间。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种快速一致性图像变换方法，其特征在于，包括如下步骤：

步骤A，设定一致性误差预设值及最大迭代次数，并根据源控制点与目标控制点的对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，然后执行步骤B；

步骤B，根据所述正向变换函数及反向变换函数计算相应的正向一致性误差及反向一致性误差，然后执行步骤C；

步骤C，判断所述正向一致性误差及反向一致性误差，并统计迭代次数；如果正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，则结束；否则执行步骤D；

步骤D，根据所述正向一致性误差及反向一致性误差调整正向变换函数及反向变换函数，然后返回步骤B；

所述步骤B具体包括如下步骤：

用户输入n对控制点对，利用薄板样条基函数构造变换函数，得到正向变换函数f₁(x)，以及反向变换函数f₂(x)；

所述正向一致性误差及反向一致性误差通过如下公式计算：

正向一致性误差：

反向一致性误差：

其中，n₁,n₂分别为正向及反向变换初始的均匀离散网格点坐标。

2.如权利要求1所述的快速一致性图像变换方法，其特征在于，所述步骤C包括如下步骤：

计算最大正向一致性误差max_δ₁及最大反向一致性误差max_δ₂：

$max\_{\mathrm{\δ}}_1=\underset{i}{max}\left({\mathrm{\δ}}_1^i\right);$

并统计迭代次数；

判断是否max(max_δ₁,max_δ₂)≤ε，ε为设定的一致性误差预设值，或者迭代达到设定的最大迭代次数；如果是，则结束，否则执行步骤D。

3.如权利要求1所述的快速一致性图像变换方法，其特征在于，正向变换函数及反向变换函数按如下公式调整：

f₁(x)＝f₁(x)+α*δ₁；

f₂(x)＝f₂(x)+β*δ₂；

其中α、β分别是正向变换函数及反向变换函数调整的步长，第k次调整步长设置为：

α_k+1＝ηα_k

β_k+1＝ηβ_k

其中，η是预设的调整因子，0＜η＜1。

4.一种快速一致性图像变换系统，其特征在于，包括：函数构造模块、误差计算模块、判断模块、函数调整模块；其中：

函数构造模块用于设定一致性误差预设值及最大迭代次数，并根据源控制点与目标控制点对应关系，构造从源控制点映射到目标控制点的正向变换函数和从目标控制点映射到源控制点的反向变换函数，然后跳转至误差计算模块；

误差计算模块用于根据所述正向变换函数及反向变换函数计算相应的正向一致性误差及反向一致性误差，然后跳转至判断模块；

判断模块用于判断所述正向一致性误差及反向一致性误差，并统计迭代次数；如果正向一致性误差及反向一致性误差都不大于设定的一致性误差预设值，或者迭代达到设定的最大迭代次数，则结束；否则跳转至函数调整模块；

函数调整模块用于根据所述正向一致性误差及反向一致性误差调整所述正向变换函数及反向变换函数，然后跳转至误差计算模块；

所述误差计算模块具体用于利用薄板样条基函数构造变换函数，得到正向变换函数f₁(x)，以及反向变换函数f₂(x)；

所述正向一致性误差及反向一致性误差通过如下公式计算：

正向一致性误差：

反向一致性误差：

其中，n₁,n₂分别为正向及反向变换初始的均匀离散网格点坐标。

5.如权利要求4所述的快速一致性图像变换系统，其特征在于，所述判断模块用于：

计算最大正向一致性误差max_δ₁及最大反向一致性误差max_δ₂：

$max\_{\mathrm{\δ}}_1=\underset{i}{max}\left({\mathrm{\δ}}_1^i\right);$

并统计迭代次数；

判断是否max(max_δ₁,max_δ₂)≤ε，ε为设定的一致性误差预设值，或者迭代达到设定的最大迭代次数；如果是，则结束，否则跳转至函数调整模块。

6.如权利要求4所述的快速一致性图像变换系统，其特征在于，正向变换函数及反向变换函数按如下公式调整：

f₁(x)＝f₁(x)+α*δ₁；

f₂(x)＝f₂(x)+β*δ₂；

其中α、β分别是正向变换函数及反向变换函数调整的步长，第k次调整步长设置为：

α_k+1＝ηα_k

β_k+1＝ηβ_k

其中，η是预设的调整因子，0＜η＜1。