CN102262781A - 一种基于单元分解光流场的喷墨印花纹理图像配准方法 - Google Patents

Abstract

一种基于单元分解光流场的喷墨印花纹理图像配准方法，包括以下步骤：①输入参考图像和变形图像，设置迭代实施参数；②对配准位移向量的系数矩阵进行不完全Cholesky分解；③置基函数为线性基函数，根据迭代误差值结果，采用预优共轭梯度算法求解单元分解系数列向量；④计算全局误差估测值，进行尺度空间调整，得到粗尺度空间；⑤置基函数为二阶基函数，在粗尺度空间上，根据局部误差估测结果，求解单元分解系数列向量，得到细尺度空间；⑥在细尺度空间上，将单元分解系数列向量堆叠展开，得到配准位移向量，完成配准。本发明能有效提高特征纹理曲线配准表征能力、提高在噪声环境下的配准精度和速度，适用于喷墨印花纹理图像的配准处理。

Description

一种基于单元分解光流场的喷墨印花纹理图像配准方法

技术领域

本发明涉及喷墨印花纹理图像配准领域，尤其是一种喷墨印花纹理图像配准方法。

背景技术

喷墨印花织物是一种织物图案成型较为复杂的织物类型，它采用高密度数码喷墨印花技术，可以极大提高织物画面的精细度和色彩丰富程度。喷墨印花织物的纹理层次丰富、质地细腻逼真，足可表现油画、山水画等逼真画面和纹理效果，已被广泛应用于云锦、丝绸等高端纺织产品中。由于喷墨印花纹理是由喷嘴将染料在织物上喷绘而成，通常具有复杂的细微油墨组织结构，组织结构的好坏是获得精细花纹的先决条件，它将直接影响印花面料的质量，因此需要通过跟踪喷墨印花纹理的表现行为来控制喷嘴的绘制效果。

为了确保喷墨印花纹理的绘制效果，喷墨印花纹理CAD系统需要对不同时刻视频采集的喷墨印花纹理进行图像配准，并根据配准的情况对喷墨纹理的表现行为和特征进行控制。喷墨印花纹理图像配准是寻找从源图像变形到目标图像的最佳空间变换过程。图像配准的效果好坏也将直接影响后续的纹理分割和三维纹理重建。

目前，图像配准按照变换的性质可将其分为刚性配准和弹性配准两大类。由于喷墨纹理组织结构复杂，形态精细且不规则，单纯依靠旋转和平移的刚性变换往往达不到喷墨纹理结构的精确配准要求，无法表示局部复杂结构的变形细节，特别是对于喷墨液滴形成复杂图案纹理的非线性变形行为，因此必须采用具有更高自由度的弹性配准方法。

基于光流场的弹性配准模型是一类用于提高可变形物体配准变化能力的方法，并在医学图像、遥感图像等领域中得到了广泛的应用。基于光流场的弹性图像配准模型可以用下面的方程表示：

$E={\∫}_{\mathrm{\Ω}}{|T(x+U)-R|}^2\mathrm{dx}+\mathrm{\α}{\∫}_{\mathrm{\Ω}}{|\▿U|}^2\mathrm{dx}---\left(1\right)$

其中，区域Ω是满足Lipschitz边界条件的有界开集，函数R是参考图像，T是变形图像，R和T是视频序列中的相邻图像对，U是R和T之间的映射变换位移场。

光流场图像配准模型式(1)的2个组成项分别对应于位移场U的逼近度约束和光滑度约束，其中逼近度约束确保相邻待配准图像R和T在Ω上偏差较小，光滑度约束确保位移场U在Ω上正则光滑，α是控制两个约束项作用的调节参数。

然而，直接利用现有光流场配准模型式(1)对喷墨印花纹理图案进行配准将存在一定不足，对喷墨印花图案的精确配准存在较大困难，具体表现为：①当待配准喷墨印花图像的局部区域具有较大变形位移时，光流场配准方法将会产生较大偏差，甚至失效。②光流场配准模型受噪声影响很大，仅添加全局光滑度约束不能较好地保持图像的局部不连续性，从而容易导致喷墨印花纹理图像在配准演化过程中，因局部细节模糊而丢失重要的特征信息。

如何提高对大变形位移的处理能力，以及确保含噪环境下的稳健性是光流场弹性配准模型应用于喷墨印花纹理图像配准中亟需解决的问题。文献[1]M.Reuter，“Hierarchical shape segmentation and registration via topological features oflaplace-beltrami eigenfunctions，”International Journal of Computer Vision，vol.89，no.2-3，pp.287-308，2010中提出一种基于阶谱分层的光流场模型，通过阶谱分层调整对大位移场的估计求解效果很好。文献[2]T.Brox，J.Malik，“LargeDisplacement Optical Flow：Descriptor Matching in Variational Motion Estimation，”IEEE Transactions on Pattern Analysis and Machine Intelligence，vol.33，no.3，pp.500-513，2011中采用扩展相位技术，能够在亚像素级别上对大尺度位移场做精确计算。文献[3]A.Doshia，A.G.Bors，“Smoothing of optical flow using robustifieddiffusion kernels，”Image and Vision Computing，vol.28，no.12，pp.1575-1589，2010中采用各项异性扩散方程对光流场模型的光滑正则项进行改进，同时结合运动模糊图像复原方法，改善了光流场模型造成的图像模糊。文献[4]P.M.Jodoin，M.Mignotte，“Optical-flow based on an edge-avoidance procedure，”Computer Vision andImage Understanding，vol.113，no.4，pp.511-531，2009中提出了一种基于目标的局部窗口级Lucas-Kanade光流估计方法，该方法能够较好地保持纹理图像的特征，并且对配准图像的模糊赝像也具有很好的抑制作用。

光流场弹性配准模型是典型的非线性二次泛函，具有很强的非凸性，从而导致其在数值计算时存在一定的难度。目前，国内外很多学者采用有限差分法对光流场弹性配准模型进行数值求解。然而，由于有限差分法是基于网格的数值方法，它在处理喷墨印花纹理中复杂拓扑结构的不连续区域(如纹理曲线的分裂与合并)和小尺度变形方面遇到较大困难。为了达到较高的变形精度，需要对计算网格实施精细单元的剖分，虽然能够取得较好的计算结果，但计算开销比较大，不适合处理具有大容量特征的喷墨印花视频图像序列的配准工作。此外，喷墨印花纹理图像在视频获取过程中，将不可避免地包含噪声，噪声将恶化图像的质量，有限差分配准方法在计算过程中容易出现网格畸变的情况，从而影响配准解的精度，降低配准的效果甚至导致配准失效。

发明内容

为了克服已有喷墨印花纹理图像配准方法对特征纹理曲线的配准表征能力弱、在噪声环境下配准精度低，以及计算开销大、运行时间较长的不足，本发明提供一种有效提高特征纹理曲线配准表征能力、提高在噪声环境下配准精度，以及减少计算开销、缩短运行时间的基于单元分解光流场的喷墨印花纹理图像配准方法。

本发明解决其技术问题所采用的技术方案是：

一种基于单元分解光流场的喷墨印花纹理图像配准方法，所述图像配准方法包括以下步骤：

步骤①：首先输入参考图像R和变形图像T，置迭代序标k＝0，迭代误差阀值η_ε，单元分解系数列向量M_k＝0，初始尺度参数h，全局表征误差阀值η_τ和局部表征误差阀值η_ρ，位移向量U_k＝0；

步骤②：计算U_k的系数矩阵COEF，并对系数矩阵COEF进行不完全Cholesky分解，得到预优矩阵；

步骤③：置k＝k+1，置基函数为线性基函数ψ＝(1 x y)，其中，x和y是二维坐标系，采用预优共轭梯度算法对k次迭代时的单元分解系数列向量M_k＝(GM_k-1+K)(COFF)^-1进行迭代求解，并计算迭代误差值err_P，其中，G是2r×2r对称正定矩阵，K是长度为2r的列向量，r是基函数的项数，M_k-1是k-1次迭代时的单元分解系数列向量。

如果迭代误差值err_P大于η_ε，则转至步骤③；否则转至步骤④；

步骤④：计算全局误差估测值err_G，如果该值大于η_τ时，则需进一步进行尺度空间调整，置尺度参数h＝h/2，调整η_ε取值，并转至步骤③；否则得到粗尺度空间h_coarse，转至步骤⑤进行局部自适应调整；

步骤⑤：置基函数为二阶基函数ψ＝(1 x y x² xy y²)，其中，x和y是二维坐标系，并在粗尺度空间h_coarse上调用预优共轭梯度算法对M_k进行迭代求解，得到局部自适应调整的M_k，并计算局部误差估测值err_L；

如果局部误差估测值err_L大于η_ρ时，则转至步骤④调整尺度空间因子，否则输出最终的细尺度空间h_fine，并转至步骤⑥；

步骤⑥：输出系数向量M_k，并将M_k堆叠展开，得到位移向量U_k，完成配准。

所述的步骤②中，U_k是配准能量表达式的约束变量，配准能量表达式定义为：

其中，Ω是矩形有界开区域，p和q是任意两个相邻面片，{Ω_p|p∈(1，...，N)}是一个有交叠的方形面片的覆盖集合，覆盖了计算域Ω，N是面片集合的势；

E(U_k)的三个组成分量E₁(U_k)、E₂(U_k)和E₃(U_k)关于单元分解系数列向量

和

的定义分别为：

$E_1\left(U_k\right)={[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]}^T{G}^{\left(p\right)}[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]-2{\left[K^{\left(p\right)}\right]}^T[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]---\left(10\right)$

$E_2\left(U_k\right)={\left[M_k^{\left(p\right)}\right]}^T{D}^{\left(p\right)}\left[M_k^{\left(p\right)}\right]---\left(11\right)$

$E_3\left(U_k\right)=2{\mathrm{\δ}}_{\mathrm{pq}}{\left[M_k^{\left(p\right)}\right]}^T{S}^{(p,q)}\left[M_k^{\left(p\right)}\right]-2{\left[M_k^{\left(p\right)}\right]}^T{W}^{(p,q)}\left[M_k^{\left(q\right)}\right]---\left(12\right)$

其中，δ_pq是Kronecker函数，当p＝q时，δ_pq＝1，否则δ_pq＝0，G^(p)是2r×2r对称正定矩阵，K^(p)是长度为2r向量，D^(p)、S^(p，q)和W^(p，q)是2r×2r对称正定矩阵，它们的各自组成分量为：

其中，μ，γ∈(1，...，r)是U分量中的基函数指标，r是基函数的项数。将K^(p)、G^(p)、D^(p)、S^(p，q)和W^(p，q)分量堆叠得到矩阵K、G、D、S和W。

进一步，所述的步骤②中，系数矩阵COEF定义为：COEF＝G+D+S-W，COEF是大型稀疏矩阵，对系数矩阵COEF按LL^T+Q作不完全Cholesky分解，其中，L是下三角矩阵，Q是剩余矩阵，然后将矩阵LL^T作为共轭梯度法的预优矩阵。

本发明的技术构思为：

1)采用具有紧支撑特性的单元分解法，即通过微分流形中的“有限覆盖技术”，首先构造喷墨印花纹理图像中有界区域合适的有限覆盖面片集合，并在每个独立覆盖面片上逼近局部变形函数，然后再将各覆盖面片“粘合”，从而最终形成对变形函数的全局逼近。采用单元分解法不仅解决了光流场配准模型在处理不连续曲线特征上的困难，而且还有效避免局部拓扑纹理特征被光滑而模糊的配准结果。

2)为了实现图像配准模型对特征纹理曲线的整体光滑表征，采用阶谱分层策略来控制单元分解网格的全局配准效果。阶谱分层策略可以在不同尺度空间上，对大位移场进行估计求解，能够在亚像素级别上对大尺度位移场做精确计算，实现“由粗到精”的配准计算求解，从而取得较好的配准效果。

3)采用一种新的单元分解结构特征能量项，提高阶谱分层和基函数自适应调整对单元分解配准模型的有效控制，从而提高配准方法在噪声环境下的细节配准效果。结构特征能量项分为局部能量项和全局能量项。局部特征能量项主要用于减小配准方法在变形区域上产生的局部配准偏差，全局特征能量项用于提高配准方法在噪声环境下的配准精度，从而确保模型对含噪图像和局部变形位移图像的配准处理具有较好的一致性，克服原有弹性配准模型存在的不足。

4)为了解决光流场模型采用有限差分网格数值求解方法时，在处理网格划分问题上存在的困难，采用基于单元分解的光流场数值求解方法。该方法只需要节点的信息，不需要节点与节点之间相互联系的信息，这样很容易在复杂计算区域内布置节点，从而提高其对单元分解网格生成过程的整体和局部的灵活度控制，有效提高网格单元对纹理图案的配准表征效果。

本发明不仅能够克服含噪环境下纹理配准模糊的缺陷，而且能够增强单元分解网格对纹理图案的表征效果，提高单元分解光流场模型对精细纹理边缘的配准精度，从而可以获得稳定、高精度的配准数值结果，非常适合含噪环境下喷墨印花纹理中复杂变形结构的配准操作。

本发明有益的效果是：

1)利用单元分解光流场技术对喷墨印花纹理图像实施配准建模，可以有效解决喷墨印花纹理图像在噪声环境下配准精确度低，并减少其对局部特征纹理区域的变形配准所产生的偏差缺陷。

2)对单元分解光流场采用阶谱分层和基函数自适应调整策略，通过对单元分解网格生成过程实施局部和整体的灵活度控制，可以有效解决配准模型在处理不连续特征问题上的困难，有效避免局部拓扑纹理特征被光滑而模糊的配准结果。

3)通过引入控制光滑性效果的结构特征能量项，可以增强单元分解网格对纹理图案的表征效果，提高单元分解光流场模型对精细纹理边缘的配准精度。

4)采用基于单元分解的光流场数值求解方法，便于在复杂计算区域内布置节点，减少了因网格细分而产生的计算开销，有效缩短了方法的运行时间。

附图说明

图1为本发明方法的执行流程图；

图2为权函数

的3个特征函数；其中图2(a)是函数

(λ)；图2(b)是函数图2(c)是函数

图3为待配准图像对和变形区域光流场矢量图；其中图3(a)是参考图像；图3(b)是印花喷嘴操作失误导致图3(a)的变形图像，图3(a)中的方框区域标记了图3(b)变形的两个主要区域；图3(c)和图3(d)分别给出了图3(a)中变形方框A区域和B区域上的光流场矢量流图；

图4为本发明方法在不同迭代次数和基函数下的配准效果比较示意图；其中图4(a)-(c)和图4(d)-(f)分别是迭代100次和200次时采用一阶基函数的配准结果，图4(g)-(i)是迭代200次时采用二阶基函数的配准结果；

图5为本发明方法与两种典型的光流场配准方法(文献[1]中具有阶谱分层特性的Reuter光流场方法和文献[4]中具有噪声抑制特性的Jodoin光流场方法)对含噪喷墨印花图像的配准结果，以及形变区域A和B的局部放大效果图；其中图5(a)和图5(b)分别是基准图像和含噪待配准图像，图5(c)-(e)分别是Reuter方法、Jodoin方法和本发明方法的配准结果对比。

图6为图5中的三种方法对图5(a)中的强变形区域A配准后，在误差精度、信噪比和运算时间上的性能比较示意图。

图7为图5中的三种方法对图5(a)中的弱变形区域B配准后，在误差精度、信噪比和运算时间上的性能比较示意图。

图8为含噪环境下(噪声方差σ＝40)，本发明方法与两种典型的非光流场配准方法(文献[5]J.P.KERN，M.S.PATTICHIS，“Robust multispectral imageregistration using mutual-information models，”IEEE Transactions on Geoscience andRemote Sensing，vol.45，no.5，pp.1494-1505，2007中基于区域互信息配准的Kern方法和文献[6]Y.YANG，X.GAO，“Remote sensing image registration via activecontour model，”International Journal of Electronics and Communications，vol.63，no.4，pp.227-234，2009中基于特征轮廓配准的Yang方法)对图5(b)进行配准后的三维特征结构比较示意图；其中，图8(a)是基准图像；图8(b)-(d)分别为Kern方法、Yang方法和本发明方法的配准结果对比；

图9为图8中的三种方法对图5(b)中的变形区域配准后，在信噪比和运算时间上的性能比较示意图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图9，一种基于单元分解光流场的喷墨印花纹理图像配准方法，所述图像配准方法包括以下步骤：

步骤①：首先输入参考图像R和变形图像T，置迭代序标k＝0，迭代误差阀值η_ε，单元分解系数列向量M_k＝0，初始尺度参数h，全局表征误差阀值η_τ和局部表征误差阀值η_ρ，位移向量U_k＝0；

步骤②：计算U_k的系数矩阵COEF，并对系数矩阵COEF进行不完全Cholesky分解，得到预优矩阵；

步骤③：置k＝k+1，置基函数为线性基函数ψ＝(1 x y)，其中，x和y是二维坐标系，采用预优共轭梯度算法对k次迭代时的单元分解系数列向量M_k＝(GM_k-1+K)(COFF)^-1进行迭代求解，并计算迭代误差值err_P，其中，G是2r×2r对称正定矩阵，K是长度为2r的列向量，r是基函数的项数，M_k-1是k-1次迭代时的单元分解系数列向量。

如果迭代误差值err_P大于η_ε，则转至步骤③；否则转至步骤④；

步骤④：计算全局误差估测值err_G，如果该值大于η_τ时，则需进一步进行尺度空间调整，置尺度参数h＝h/2，调整η_ε取值，并转至步骤③；否则得到粗尺度空间h_coarse，转至步骤⑤进行局部自适应调整；

步骤⑤：置基函数为二阶基函数ψ＝(1 x y x² xy y²)，其中，x和y是二维坐标系，并在粗尺度空间h_coarse上调用预优共轭梯度算法对M_k进行迭代求解，得到局部自适应调整的M_k，并计算局部误差估测值err_L；

如果局部误差估测值err_L大于η_ρ时，则转至步骤④调整尺度空间因子，否则输出最终的细尺度空间h_fine，并转至步骤⑥。

步骤⑥：输出系数向量M_k，并将M_k堆叠展开，得到位移向量U_k，完成配准。

为了克服光流场配准模型(1)在处理不连续曲线特征上的困难，本发明采用具有紧支撑特性的单元分解法。单元分解法的核心思想是采用具有紧支撑特性的有限覆盖技术，构造有界区域合适的有限覆盖面片集合，并在每个独立覆盖面片上逼近局部函数，然后再将各覆盖面片“粘合”，从而最终形成对函数的全局逼近。

采用单元分解法后，光滑度约束将主要对各个独立覆盖域(同质连续区域)有效，而对于覆盖域交叠区域(异质不连续区域)的影响较小，这样就可以克服模型式(1)在处理不连续特征问题上的困难，有效避免局部拓扑纹理特征被光滑而模糊的配准结果。

定义1设Ω是矩形有界开区域，{Ω_p|p∈(1，...，N)}是一个有交叠的方形面片的覆盖集合，覆盖了计算域Ω，p是面片集合的指标，N是面片集合的势。定义2设{Ω_p}是覆盖域集合，

是属于{Ω_p}的一个Lipschitz单元分解，即对于

满足则称

是附属于{Ω_p}的一个单元分解，

被称作单元分解的权函数。

为了使逼近函数U具有足够的分片光滑效果，权函数

应确保在各个面片区域Ω_p上连续，通常权函数是x和y两个分量的乘积，即有：

对于二维喷墨印花纹理图像的配准模型，权函数分量

应在节点中心处取最大值，并由近及远逐渐衰减，且在支撑域半径外衰减为零或可忽略不计。

针对权函数的上述要求，本发明设计采用一种新的权函数，该权函数的分量表达式如下：

图2给出了权函数的3个特征函数的三维表示。

定义3设是属于{Ω_p}的一个单元分解，

是给定Ω_p上的函数空间，则称

为单元分解空间，V^(p)称作局部逼近空间，通常V^(p)被选为多项式空间。

本发明采用完备多项式作为基函数，当阶数越高时，收敛率越好，精度越高。对于二维喷墨印花纹理图像的配准模型，本发明采用的线性和二阶基函数ψ^(p)为：

●线性基函数：ψ^(p)＝(1，x，y)            (3)

●二阶基函数：ψ^(p)＝(1，x，y，x²，xy，y²)(4)

定义4设

是给定Ω_p上的一组基函数，r是基函数的项数，

是一组标量系数，则称局部逼近空间V^(p)(Ω_p)中的函数

是Ω_p上的局部逼近函数。

由定义1-4，可得单元分解特性如下：设{v^(p)|p∈(1，...，N)}是{Ω_p}上的局部逼近函数集，

是附属于{Ω_p}的一个单元分解，则全局逼近函数U是局部逼近函数集{v^(p)}的加权和，即有

利用单元分解特性，图像配准模型式(1)可以分解为覆盖面片区域集{Ω_p|p∈(1，...，N)}上的能量值E^(P)之和，即有：

$E={\mathrm{\Σ}}_{p=1}^N{E}^{\left(P\right)}\left(U\right)={\mathrm{\Σ}}_{p=1}^N{{\∫}_{{\mathrm{\Ω}}_p}|T(x+U)-R|}^2\mathrm{dx}+\mathrm{\α}{\mathrm{\Σ}}_{p=1}^N{{\∫}_{{\mathrm{\Ω}}_p}|\▿U^{\left(p\right)}|}^2\mathrm{dx}---\left(5\right)$

图像配准模型式(5)实现了位移场U在覆盖面片区域集上的正则光滑度约束，从而保证其在各个独立面片上的分片光滑性。

进一步地，为了实现图像配准模型式(5)对局部不连续结构的表征，本发明在相邻覆盖面片的重叠区域上引入一种新的结构特征能量项，以实现其对位移向量场U不连续表示的正则化约束，从而提高位移向量场U对复杂纹理结构特征的表征能力。结构特征能量项分为局部能量项和全局能量项。局部特征能量项主要用于减小配准方法在变形区域上产生的局部配准偏差，全局特征能量项用于提高配准方法在噪声环境下的配准精度，从而确保模型对含噪图像和局部变形位移图像的配准处理具有较好的一致性，克服原有弹性配准模型存在的不足。

定义5设p和q是任意两个相邻面片，则面片交叠区域Ω_(p，q)上的局部结构特征能量

其中，

是偏微分算子，

是全偏差阶数，k是最高Sobolev微分平滑阶数。

局部结构特征能量可以对位移向量场U在局部交叠区域Ω_(p，q)上做局部不连续正则约束，即在区域Ω_(p，q)上将累积到k阶偏导数的偏差和作为约束惩罚项，在本发明我们选取k＝2，即累计到2阶偏导数。

定义6设面片p与相邻的4个面片存在交叠区域，q∈(1，...，4)是相邻面片的指标，则全局结构特征能量

$E_s=\frac{1}{2}{\mathrm{\Σ}}_{p=1}^N{E}_s^{\left(p\right)}\left(U\right)=\frac{1}{2}{\mathrm{\Σ}}_{p=1}^N{\mathrm{\Σ}}_{q=1}^4{E}_s^{(p,q)}\left(U\right)---\left(7\right)$

通过对全局结构特征能量的约束求解，可以实现位移向量场U在所有交叠区域上的不连续正则约束，从而提高位移场U配准所得图像中纹理特征细节的全局表征效果。

综上所述，基于单元分解的光流场图像配准模型表示如下：

$E={\mathrm{\Σ}}_{p=1}^N{{\∫}_{{\mathrm{\Ω}}_p}|T(x+U)-R|}^2\mathrm{dx}+\mathrm{\α}{\mathrm{\Σ}}_{p=1}^N{{\∫}_{{\mathrm{\Ω}}_p}|\▿U^{\left(p\right)}|}^2\mathrm{dx}+\frac{1}{2}\mathrm{\β}{\mathrm{\Σ}}_{p=1}^N{\mathrm{\Σ}}_{q=1}^4{E}_s^{(p,q)}\left(U\right)---\left(8\right)$

由单元分解原理可知，在给定基函数和权函数后，式(8)中U是由单元分解的标量系数组决定的，为了矩阵计算简便，本发明对标量系数进行向量堆叠。

对于任意给定面片p∈(1，...，N)和U的独立分量U_i(i＝1，2)，首先将U_i的单元分解标量系数

堆叠成列向量接着将U所有分量的单元分解列向量堆叠成列向量

最后将所有面片上的U全部分解列向量堆叠成列向量M。

式(8)中，项是关于U的非独立项，无法直接用单元分解系数向量来表示。本发明利用泰勒级数对该式进行一阶逼近展开，使其成为U的独立项，以便于进行模型的二次型优化求解。

设U_k和U_k-1是U的本次迭代k和前次迭代k-1时的值，M_k和M_k-1是对应的单元分解系数列向量，u_k是位移场U_k和U_k-1在两次迭代间的微小增量。由光流法可知，T(x+U_k)≈T(x+U_k-1+u_k)，对式(8)中的T(x+U)在点x附近的一阶泰勒级数展开，则位移场U_k的能量表达式为：

略去和u_k无关的常量项，式(9)中E₁(U_k)、E₂(U_k)和E₃(U_k)关于单元分解系数列向量

和

的定义分别为：

$E_1\left(U_k\right)={[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]}^T{G}^{\left(p\right)}[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]-2{\left[K^{\left(p\right)}\right]}^T[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]---\left(10\right)$

$E_2\left(U_k\right)={\left[M_k^{\left(p\right)}\right]}^T{D}^{\left(p\right)}\left[M_k^{\left(p\right)}\right]---\left(11\right)$

$E_3\left(U_k\right)=2{\mathrm{\δ}}_{\mathrm{pq}}{\left[M_k^{\left(p\right)}\right]}^T{S}^{(p,q)}\left[M_k^{\left(p\right)}\right]-2{\left[M_k^{\left(p\right)}\right]}^T{W}^{(p,q)}\left[M_k^{\left(q\right)}\right]---\left(12\right)$

其中，δ_pq是Kronecker函数，当p＝q时，δ_pq＝1，否则δ_pq＝0，K^(p)是长度为2r向量，G^(p)、D^(p)、S^(p，q)和W^(p，q)是2r×2r对称正定矩阵，它们的各自组成分量为：

其中，δ_pq是Kronecker函数(当p＝q时，δ_pq＝1，否则δ_pq＝0)，其中，μ，γ∈(1，...，r)是U分量中的基函数指标。

式(10，11，12)都是

的二次型，将所有面片上的U全部分解列向量堆叠成列向量M_k。由变分最优化原理可知，式(9)的最优化求解等价于如下方程的求解，即：

G[M_k-M_k-1]+D[M_k]+S[M_k]-W[M_k]＝K     (18)

整理上式可得：

(G+D+S-W)M_k＝GM_k-1+K                (19)

令系数矩阵COEF＝G+D+S-W，COEF是大型稀疏矩阵，且条件数很大，收敛速度很慢。为此，本发明采用基于不完全Cholesky分解的预优共轭梯度法对线性方程组(19)进行求解，即首先对系数矩阵COEF按LL^T+Q作不完全Cholesky分解，其中，L是下三角矩阵，Q是剩余矩阵，然后再将矩阵LL^T作为共轭梯度法的预优矩阵，并按照共轭梯度法进行迭代求解。

由于喷墨印花纹理图像中包含大量拓扑结构复杂的精细纹理曲线，为了实现图像配准模型对特征曲线的全局和局部区域上的有效表征，需要对单元分解网格生成的灵活度和基函数阶数的高低控制提出较高的要求，单元分解网格的细密程度和基函数阶数的高低差异将直接影响图像配准模型方法对纹理图案的配准精度。

为此，本发明采用阶谱分层策略来全局控制单元分解网格的生成效果，阶谱分层策略可以在不同尺度空间上，实现“由粗到精”的计算求解，即将在粗尺度空间上的迭代结果作为下一次精尺度空间上的迭代初始值，阶谱分层策略具有快速收敛的优点。

此外，相同尺度上纹理曲线配准精度的局部调整则可以通过调节基函数的阶数来实现。当纹理曲线的局部精细度较高时，本发明采用二阶基函数式(4)，反之则采用线性基函数式(3)。缺省情况下，采用线性基函数以确保本发明方法具有较快的计算速度。设η_ρ是局部表征误差阀值，当同一尺度上相邻两次迭代值之间的误差值大于η_ρ时，将采用二阶基函数式(4)。

图1给出了本发明方法的执行流程图；

图2给出权函数

的3个特征函数(即

和

)的三维表示。由权函数的特征函数可知，本发明使用的权函数式(2)满足如下特性：

(1)当λ＜-1或λ＞1时，

即(λ)在支撑域半径外衰减为0；

(2)当λ＝0时，

即

在节点中心处的作用效果最大；

(3)当0≤λ≤1时，

即

关于节点中心处是对称的；

(4)即

的作用效果与距离节点中心处的λ取值无关，是恒定的。

本发明提出的权函数不仅计算简单(仅涉及多项式计算)，而且能自动满足正规化条件，由权函数的上述特性可知，即对

(x，y)所属的4个非零权函数所支撑的区域满足

图3给出了本发明的一个配准实施例，其中图3(a)是参考图像，图3(b)是印花喷嘴操作失误导致图3(a)的变形图像，图3(a)中的方框区域标记了图3(b)变形的两个主要区域。图3(c)和图3(d)分别给出了图3(a)中变形方框A区域和B区域上的光流场矢量流图。光流场矢量演化的方向和强度反映了喷墨液滴的晕染效果和变形趋势，由图3(c)和图3(d)可知，图3(a)中变形方框A区域中喷墨纹理的整体变形程度要高于B区域。

图4给出了采用本发明对图3的配准结果，其中图4(a)-(c)和图4(d)-(f)分别是迭代100次和200次时采用一阶基函数的配准结果，图4(g)-(i)是迭代200次时采用二阶基函数的配准结果。从图中情况可以看出：①单元分解光流场模型在迭代100次和200次后形成的图4(b)和图4(e)，在单元分解网格的有效夹逼作用下，都取得了与参考图像(图4(a)和图4(d))在相关区域近似匹配的配准图像(图4(c)和图4(f))。这主要是因为图3(a)中强变形区域A和弱变形区域B中的拓扑纹理结构通常都处于单元分解的交叠区域上，由于本发明方法添加了配准模型对交叠区域上的结构特征能量约束项，使得本发明方法对不连续特征纹线的配准演化具有较好的控制能力，从而兼顾了正则光滑约束和特征保留，有效提高了配准模型对不连续特征纹线的表征能力，配准效果令人比较满意。②随着迭代次数的增加，阶谱分层尺度因子h不断减小，单元分解网格对特征纹理的配准效果将不断提高，图4(d)-(e)中单元分解网格对特征纹理的逼近程度要高于图4(a)-(b)中的单元分解网格，因此图4(f)在变形区域A和B中的配准效果明显优于图4(c)的配准效果。③阶谱分层尺度因子h缩减到一定程度后，改变基函数的阶数可以有效提高特征纹理的局部配准效果。图4(g)-(i)中采用二阶基函数对局部网格中的特征纹理配准逼近效果要优于图4(d)-(f)，配准后图像中的特征细节表现情况更平顺。

表1给出了图4中的配准实例在不同迭代次数下的尺度因子h、预优共轭梯度算法的迭代误差值err_P、全局误差估测值err_G和局部误差估测值err_L。分析表1数据可得：①当迭代次数增加时，err_P的取值单调减小，表明采用预优共轭梯度算法在计算配准系数矩阵时的线搜索方向正确，取得较好的配准收敛效果；②阶谱分层尺度因子h随迭代次数的增加将不断减小，实现了变形图像“由粗到精”的分层配准；③err_G的取值随阶谱分层因子h的减小而不断减小，但在相同阶谱层次上的变化幅度较小。当迭代次数在200次左右时，err_G随h变化的幅度非常小，这表明h调正到一定程度后，单纯改变h而形成的全局配准效果并不理想，此时就需要通过调整单元分解模型的基函数阶数，进行局部误差估测驱动下的局部配准调整；④err_L在h较大时，变化幅度较小，但当h不断减小后，err_L的变化幅度将变化很大，这表明在细尺度上通过调整局部误差估测可以获得更好的配准效果。

表1不同迭代次数下的评价项比较

图5给出了在含噪环境下(噪声方差σ＝40)，采用本发明方法与两种典型的光流场配准方法(文献[1]中具有阶谱分层特性的Reuter光流场方法和文献[4]中具有噪声抑制特性的Jodoin光流场方法)对含噪喷墨印花图像的配准结果，以及形变区域A和B的局部放大效果图。图5(a)和4(b)分别是基准图像和含噪待配准图像(噪声方差σ＝40)，图5(c)-(e)分别为Reuter方法、Jodoin方法和本发明方法的配准结果对比。从图中情况可以看出：Reuter方法对大位移形变区域A的配准效果明显优于Jodoin方法，Jodoin方法对区域A的配准操作存在较大偏差。但Reuter方法对噪声较敏感，导致变形区域A和B中许多局部边缘细节模糊而错误配准，而Jodoin方法对噪声具有较好的稳健性。本发明方法由于同时结合Reuter方法和Jodoin方法的各自优势，对含噪环境下大位移形变区域的配准取得较好的效果。

图6和图7分别给出了图5中的三种方法对图5(a)中强变形区域A和弱变形区域B配准后，在误差精度、信噪比和运算时间上的性能比较。本发明选取作为衡量变形区域与配准区域之间的误差精度，其中Ψ表示变形区域，|Ψ|表示区域Ψ所含像素个数，X_i，j和Y_i，j分别表示变形区域和配准区域在像素点(i，j)上的灰度值。从图6和图7可以看出，本发明方法在误差精度上具有和Reuter方法相近的精度值，在信噪比上具有和Jodoin方法相近的比值，克服了其他两种方法的各自不足，取得了两种方法相折衷的效果。此外，本发明方法所需运算时间较其他两种方法要少15％-25％。这主要是由于其他两种方法采用差分网格数值计算，在对变形区域的配准过程中，需要实施精细单元的剖分，因此计算开销较大，而本发明方法采用单元分解的数值求解方法，减少了因网格细分而产生的计算开销。

图8给出了在含噪环境下(噪声方差σ＝40)，采用本发明方法与两种典型的非光流场配准方法(文献[5]中基于区域互信息配准的Kern方法和文献[6]中基于特征轮廓配准的Yang方法)对图5(b)进行配准后的三维特征结构比较。从图中情况可以看出：本发明方法(图6(d))要优于Kern方法(图6(b))和Yang方法(图6(c))，它能够更好地保留变形区域中配准图案细节的拓扑形状。

图9给出了图8中的三种方法配准后在信噪比和运算时间上的比较。从图中情况可以看出：本发明方法所得配准图像中含有的噪声较少，特征信号较多，因而图像的细节特征保持情况较好。此外，由于本发明方法不需要在配准过程中进行区域信息或特征提取，以及后续的匹配操作，因此减小了配准方法的复杂度，缩短了运行时间。

Claims (3)

1.一种基于单元分解光流场的喷墨印花纹理图像配准方法，其特征在于：所述图像配准方法包括以下步骤：

步骤①：首先输入参考图像R和变形图像T，置迭代序标k＝0，迭代误差阀值η_ε，单元分解系数列向量M_k＝0，初始尺度参数h，全局表征误差阀值η_τ和局部表征误差阀值η_ρ，位移向量U_k＝0；

步骤②：计算U_k的系数矩阵COEF，并对系数矩阵COEF进行不完全Cholesky分解，得到预优矩阵；

步骤③：置k＝k+1，置基函数为线性基函数ψ＝(1 x y)，其中，x和y是二维坐标系，采用预优共轭梯度算法对k次迭代时的单元分解系数列向量M_k＝(GM_k-1+K)(COFF)^-1进行迭代求解，并计算迭代误差值err_P，其中，G是2r×2r对称正定矩阵，K是长度为2r的列向量，r是基函数的项数，M_k-1是k-1次迭代时的单元分解系数列向量。

如果迭代误差值err_P大于η_ε，则转至步骤③；否则转至步骤④；

步骤④：计算全局误差估测值err_G，如果该值大于η_τ时，则需进一步进行尺度空间调整，置尺度参数h＝h/2，调整η_ε取值，并转至步骤③；否则得到粗尺度空间h_coarse，转至步骤⑤进行局部自适应调整；

步骤⑤：置基函数为二阶基函数ψ＝(1 x y x² xy y²)，其中，x和y是二维坐标系，并在粗尺度空间h_coarse上调用预优共轭梯度算法对M_k进行迭代求解，得到局部自适应调整的M_k，并计算局部误差估测值err_L；

如果局部误差估测值err_L大于η_ρ时，则转至步骤④调整尺度空间因子，否则输出最终的细尺度空间h_fine，并转至步骤⑥；

步骤⑥：输出系数向量M_k，并将M_k堆叠展开，得到位移向量U_k，完成配准。

2.根据权利要求1所述的基于单元分解光流场的喷墨印花纹理图像配准方法，其特征是，步骤②中所述U_k是配准能量表达式的约束变量，配准能量表达式定义为：

其中，Ω是矩形有界开区域，p和q是任意两个相邻面片，{Ω_p|p∈(1，...，N)}是一个有交叠的方形面片的覆盖集合，覆盖了计算域Ω，N是面片集合的势；

E(U_k)的三个组成分量E₁(U_k)、E₂(U_k)和E₃(U_k)关于单元分解系数列向量

和

的定义分别为：

$E_1\left(U_k\right)={[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]}^T{G}^{\left(p\right)}[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]-2{\left[K^{\left(p\right)}\right]}^T[M_k^{\left(p\right)}-M_{k-1}^{\left(p\right)}]---\left(10\right)$

$E_2\left(U_k\right)={\left[M_k^{\left(p\right)}\right]}^T{D}^{\left(p\right)}\left[M_k^{\left(p\right)}\right]---\left(11\right)$

$E_3\left(U_k\right)=2{\mathrm{\δ}}_{\mathrm{pq}}{\left[M_k^{\left(p\right)}\right]}^T{S}^{(p,q)}\left[M_k^{\left(p\right)}\right]-2{\left[M_k^{\left(p\right)}\right]}^T{W}^{(p,q)}\left[M_k^{\left(q\right)}\right]---\left(12\right)$

其中，δ_pq是Kronecker函数，当p＝q时，δ_pq＝1，否则δ_pq＝0，K^(p)是长度为2r向量，G^(p)、D^(p)、S^(p，q)和W^(p，q)是2r×2r对称正定矩阵，它们的各自组成分量为：

其中，μ，γ∈(1，...，r)是U分量中的基函数指标，r是基函数的项数。将K^(p)、G^(p)、D^(p)、S^(p，q)和W^(p，q)分量堆叠得到矩阵K、G、D、S和W。

3.根据权利要求2所述的基于单元分解光流场的喷墨印花纹理图像配准方法，其特征是，步骤②中所述的系数矩阵COEF定义为：COEF＝G+D+S-W，COEF是大型稀疏矩阵，对系数矩阵COEF按LL^T+Q作不完全Cholesky分解，其中，L是下三角矩阵，Q是剩余矩阵，然后将矩阵LL^T作为共轭梯度法的预优矩阵。

Images