CN103729849A - 一种数字图像变形初值的计算方法 - Google Patents

Abstract

本发明公开了一种数字图像变形初值的计算方法，针对一般变形初值估计方法在被测物体表面存在较大刚体转动或大变形的情况下失效的问题，提出了一种自动的变形初值估计方法，即通过SIFT算法匹配得到变形前后图像计算子区的对应点，并利用最小二乘法获得变形初值估计作为Newton-Rapshon方法的迭代初值。对大变形下的模拟散斑图进行了仿真计算，结果显示该方法在30°大转角的情况下仍能够有效测量位移。

Description

一种数字图像变形初值的计算方法

技术领域

本发明属于光测力学技术领域，具体涉及一种数字图像变形初值的计算方法。

背景技术

数字图像相关（Digital Image Correlation,DIC）最早出现在上世纪80年代初，由美国南卡罗来纳州大学的Peters和Ranson最早提出，该方法日趋成熟，已经成为一种广为人知并有众多应用的非接触全场光学测量方法。

为计算变形前图像中某点的位移，通常选择以该点为中心的一个正方形子区。由于参考图像正方形子区在变形后图像中形状会发生变化，因此可用形函数来近似描述变形后图像子区的形状。此时，变形前后图像子区相似程度的相关函数则可以表示为变形参数的非线性方程。此方程可以用Newton-Rapshon（后面简称为N-R）方法或其它方法求解。其它方法如遗传算法、差分进化算法虽可以用于求取变形参数，但相比N-R方法极为耗时，现以少见于应用。近30年来，对N-R方法的研究主要涉及以下几个关键问题：①相关函数选择问题；②变形初值估计问题；③形函数的选择问题；④Hessian矩阵的优化问题。由于N-R方法收敛范围在几个像素内（Knauss认为在7个像素范围内收敛），因此就要求提供较准确的初值估计。当变形较小时，简单的整像素相关搜索可以得到较准确的变形初值（此时应变参数设为0）。而当变形较大时，此方法收索到的结果相关系数值小而变得不可靠，甚至于是错误的。基于人机交互的变形初值估计方法可为任何变形状态迅速提供准确可靠的变形初值估计，避免计算量较大的搜索过程，缺点是不能实现计算过程的完全自动化。

现有技术参考文献如下：

[1]潘兵，谢惠民，续伯钦等.数字图像相关中的亚像素位移定位算法进展[J].力学进展,2005,35(3):345～351.

Pan Bing,Xie Huimin,Xu Boqin et al..Development of sub-pixel displacementsregistration algorithms in digital image correlation[J].Advances in Mechanics,2005,34(3):345～351.

[2]Peters W H,Ranson WF.Digital Imaging Techniques in Experimental StressAnalysis[J].Optical Engineering,1981,21(3):427～431.

[3]Pan B,Xie H M,Guo Z Q et al..Full-field strain measurement usingtwo-dimensional Savitzky-Golay digital differentiator in digital imagecorrelation[J].Optical Engineering,2007,46(3):033601-033601-10.

[4]潘兵，谢惠明.数字图像相关中基于位移场局部最小二乘拟合的全场应变测量[J].光学学报,2007,27(11):1980～1986.

Pan Bing,Xie Huiming.Full-field strain measurement based on local least-squarefitting for digital image correlation method[J].Acta Optica Sinica,2007,27(11):1980～1986.

[5]Schreier H M,Sutton M A.Systematic errors in digital image correlation dueto undermatched subset shape functions[J].Experimental Mechanics,1989,29(3):303～310.

[6]Ma S P,Jin G C.Digital speckle correlation method improved by geneticalgorithm[J].Acta Mechanica Solida Sinica,2003,16(4):366～373.

[7]Pan B,Xie H M.Digital image correlation method with differentialevolution[J].Optoelectronics & Laser,2007,18(1):100～103.

[8]潘兵，谢惠明，夏勇等.数字图像相关中基于可靠变形初值估计的大变形测量[J].光学学报,2009,2(29):400～406.

Pan Bing,Xie Huimin,Xia Yong et al..Large-Deformation Measurement Based onReliable Initial Guess in Digital Image Correlation Method[J].Acta Optica Sinica,2009,2(29):400～406.

[9]G.Vendroux,W.G.Knauss.Submicron deformation field measurements:Part2.Improved digital image correlation[J].Experimental Mechanics,1998,38(2):86～92.

[10]Bruck H A,McNeil S R,Sutton M A et al..Digital image correlation usingnewton-rapshon method of partial differential correlation[J].ExperimentalMechanics,1989,29(3):261～267.

[11]LOWE D G.Object Recognition from local scale-invariant features[C].TheInternational Conference on Computer Vision,1999(2):1150-1157.

[12]Peng Zhou,Kenneth E,Goodson.Subpixel displacement and deformationgradient measurement using digital image/speckle correlation[J].OpticalEngineering,2001,40(8):1613～1620.

[13]Hung PC,Voloshin PS.In-plain Strain Measurement by Digital ImageCorrelation[J].Journal of the Brazilian Society of Mechanical Sciences andEngineering,2003,25(3):215-221.

发明内容

本发明所要解决的技术问题是：提供一种数字图像变形初值的计算方法，解决了现有技术中数字图像变形初值计算方法不能实现计算过程的完全自动化以及计算精度低的问题。

本发明为解决上述技术问题采用以下技术方案：

一种数字图像变形初值的计算方法，包括如下步骤：

（1）获取被测物体表面变形前的数字图像及表面变形后的数字图像，将表面变形前的数字图像作为“参考图像”，表面变形后的图像作为“目标图像”；首先，将参考图像及目标图像均置于同一个平面坐标系内，在参考图像上，任意选取一个像素点作为网格节点，并以该网格节点为中心点，以2M+1为边长划分正方形区域，作为参考图像子区，其中M为像素点数，其次，通过特征点提取及匹配算法从目标图像中自动提取出与参考图像子区对应的像素点，以此得到一一对应的多组像素点，所述参考图像中与目标图像中任一对应的像素点之间的变换关系用“形函数”表示；

（2）根据一一对应像素点及参考图像子区中心点坐标计算确定中间计算量，所述中间计算量具体包括：参考图像子区中与目标图像子区中在同一坐标系下所有像素点的坐标值及参考图像子区中所有像素点的坐标相对于中心的变化量；

（3）根据步骤（2）中得到的中间计算量及步骤（1）中得到的一一对应的多组像素点利用最小二乘法求解所述“形函数”的参数值，将得到的“形函数”的参数值作为所述数字图像的变形初值。

所述“形函数”为：

x′=x+u+u_xΔx+u_yΔy

y′=y+v+v_xΔx+v_yΔy

其中，x′,y′为目标图像子区中除中心点外任一像素点的坐标，x,y为参考图像子区中与x′,y′对应的像素点的坐标，Δx,Δy分别为为参考图像子区中的像素点x,y到参考图像子区中心在横纵坐标轴方向的距离，u,v是参考图像子区中心到目标图像子区中心在横纵坐标轴方向的位移，u_x，u_y，v_x，v_y分别为参考图像子区中到目标图像子区中对应点的位移梯度。

所述M为大于等于20且小于等于50的整数。

所述“形函数”的参数值是通过将中间计算量代入一阶或二阶“形函数”得到的。

与现有技术相比，本发明具有以下有益效果：

1、本发明操作简单，自动化程度高。

2、不受结构表面变形大小的限制，在存在较大面内位移或扭转变形的情况下都可以使用。

附图说明

图1为本发明数字图像变形前后图像子区示意图。

图2为本发明小变形全场相关系数分布。

图3为本发明大变形全场相关系数分布。

图4为本发明SIFT特征点匹配。

图5为本发明变形前图像。

图6为本发明变形后图像。

具体实施方式

下面结合附图对本发明的结构及工作过程作进一步说明。

1、数字图像相关方法

数字图像相关方法是通过拍摄得到被测物体表面变形前后数字图像，匹配出两幅图像的指定对应点，计算出全场位移及应变的光学测量方法。变形前图像成为“参考图像”，变形后图像称为“目标图像”，其原理图如图1所示。数字图像相关方法首先预先在参考图像划分一虚拟网格用于求取离散的全场位移。图中p点为一网格节点，以其为中心(2M+1)×(2M+1)像素大小的正方形区域为参考图像子区（虚线框），其中M为像素点数，优选取值为20至50内的整数。在变形后图像中通过一定的搜索方法，按定义好的相关函数来进行相关计算，寻找出使相关函数取极值的目标图像子区（实线框）。结合目标图像子区中心点P′(x₀′,y′₀)坐标即可计算出点P(x₀, y₀)在x,y方向的位移量u,v。

变形后的目标图像子区的形状与位置相比之前都发生了变化。在采用一阶形函数的情况下，此时目标图像子区中的G′(x′,y′)与参考图像中的G(x,y)存在如下线性关系：

x′=x+u+u_xΔx+u_yΔy

                                   ⑴

y′=y+v+v_xΔx+v_yΔy

在⑴式中，Δx,Δy为点(x,y)到参考图像子区中心(x₀,y₀)的距离，u,v是参考图像子区中心在x,y方向的位移，u_x，u_y，v_x，v_y为图像子区的位移梯度。

根据文献0的研究，本文采用如下的相关函数来评判变形前后图像子区的相似程度：

$C\left(P\right)=\frac{\underset{G_s\∈S}{\mathrm{\Σ}}{\{f\left(G_s\right)-[g\left({\mathrm{\χ}}_1\left(G_s\right)\right)-w_0\left]\right\}}^2}{\underset{G_s\∈S}{\mathrm{\Σ}}{f}^2\left(G_s\right)}\left(2\right)$

式中P=(u,u_x,u_y,v,v_x,v_y,w₀)^T是变形参数矢量，w₀是亮度补偿参数。S是参考图像子区的点集，G_s代表其中的点。χ₁为公式⑴所代表的线性映射。f和g分别表示变形前后的图像。

由式⑵可知，相关函数C(P)是关于变形参数P的函数。当在目标图像中匹配到对应的图像子区时，相关函数C(P)取极小值，C(P)一阶倒数为0，亦即有：

$\▿C={\left(\frac{\∂C}{{\∂P}_i}\right)}_{i=1..7}=\frac{-2}{\underset{G_s\∈S}{\mathrm{\Σ}}{f}^2\left(G_s\right)}{\left\{\underset{G_s\∈S}{\mathrm{\Σ}}\right\{f\left(G_s\right)-\stackrel{\‾}{g}(G_s,P)\left\}\frac{\∂\tilde{g}(G_s,P)}{{\∂P}_i}\right\}}_{i=1..7}=0---\left(3\right)$

式中

上式可用N-R迭代法求解，其方程为：

$\▿\▿C\left(P_0\right)(P-P_0)=-\▿C\left(P_0\right)---\left(4\right)$

因此：

$P=P_0-\frac{\▿C\left(P_0\right)}{\▿\▿C\left(P_0\right)}---\left(5\right)$

这里P₀为变形初值估计，P为迭代后的近似值，为相关函数的梯度向量，是相关函数的二次偏导，也被称作Hessian矩阵。对Hessian矩阵的近似处理及具体过程可见文献0。由数值分析知识可知，对于非线性方程组，只有接近真值的初值估计才能在给定条件下收敛。因此可靠的变形初值估计对于N-R方法尤为重要。

2、基于SIFT的变形初值估计

通常情况下，可以通过简单的搜索来获得变形初值。下面首先介绍整像素搜索及其不足之处，随后介绍基于SIFT算法的变形初值估计。

2.1整像素位移搜索

基于目标图像子区的形状不变，仅其位置发生了变化的假设，可以在变形后图像中的搜索区域内逐点移动图像子区并计算相关系数，标准协方差相关函数值最大的点即为目标位置点。整像素位移相关搜索是以整像素为单位来进行的，因此所获得的位移是像素的整数倍。这样N-R方法的迭代初始值就可以写成P₀=(u,v,0,0,0,0,0)^T。

大多情况下，参考图像子区的变形都较小，因此以像素为单位的相关搜索可以提供准确的像素级位移估计。如图2所示，图中相关系数的最大值接近1，表明搜索可靠。问题是，某些情况下变形后图像相对于参考图像可能会有较大的刚体转动或发生大变形。图3显示了参考图像子区在发生了30°刚体转动后的全场相关系数分布，可以看出此时没有一个如图2中出现的相关高峰。在相关系数最大值为0.3829的情况下，虽然找到了相关系数最大值，但这样的目标位置显然是错误的。因此在整像素搜索不起作用的大变形情况下，就需要其它的方法来获取变形初值估计。

2.2基于SIFT算法的变形初值估计

SIFT（Scale Invariant Feature Transform）算法由Lowe0于1999年提出，并于2004年做了完善和总结。SIFT算法实现主要包含4个步骤：①建立尺度空间，检测尺度空间极值；②去除不稳定的特征点；③计算特征点的描述信息，确定特征点的方向；④生成本地特征描述符，组合建立特征向量。SIFT算法的具体实现可见文献0。

变形初值的估计就是求取公式⑴中的6个参数。待定参数可以通过解下面两个方程组直接求出：

$\mathrm{XU}=D_x\⇒\left[\begin{array}{ccc}1& {\mathrm{\Δx}}_1& {\mathrm{\Δy}}_1\\ 1& {\mathrm{\Δx}}_2& {\mathrm{\Δy}}_1\\ & \·\·\·\·\·\·& \\ 1& {\mathrm{\Δx}}_n& {\mathrm{\Δy}}_n\end{array}\right]\left(\begin{array}{c}u\\ u_x\\ u_y\end{array}\right)=\left(\begin{array}{c}{x}_1^{\′}-x_1\\ x_2^{\′}-x_2\\ \·\·\·\\ x_n^{\′}-x_n\end{array}\right)---\left(6\right)$

$\mathrm{XV}=D_y\⇒\left[\begin{array}{ccc}1& {\mathrm{\Δx}}_1& {\mathrm{\Δy}}_1\\ 1& {\mathrm{\Δx}}_2& {\mathrm{\Δy}}_1\\ & \·\·\·\·\·\·& \\ 1& {\mathrm{\Δx}}_n& {\mathrm{\Δy}}_n\end{array}\right]\left(\begin{array}{c}v\\ v_x\\ v_y\end{array}\right)=\left(\begin{array}{c}{y}_1^{\′}-y_1\\ y_2^{\′}-y_2\\ \·\·\·\\ y_n^{\′}-y_n\end{array}\right)---\left(\text{7}\right)$

上式中Δx=x-x₀,Δy=y-y₀，x₀,y₀为参考图像子区的中心坐标。式⑹可通过最小二乘法求解：

$\left(\begin{array}{c}u\\ u_x\\ u_y\end{array}\right)={\left(X^{T}X\right)}^{-1}{X}^T{D}_x---\left(8\right)$

这里(X^TX)^-1X^T为X的伪逆矩阵，对于式⑺可按同样方法得到 $\left(\begin{array}{c}v\\ v_x\\ v_y\end{array}\right)={\left(X^{T}X\right)}^{-1}{X}^T{D}_y.$ 这样，问题就变成了求取3对或更多对的参考图像与目标图像的对应点。

求取对应点的问题是一个图像特征提取及匹配的问题。由于SIFT算子对图像缩放、旋转甚至仿射变换保持不变性，因此可以利用SIFT算法找到多组对应点。一般情况下，SIFT算法会找到大于3组的对应点（如图4所示，19组对应点用直线连接），此时式⑹、⑺就可以用最小二乘法求解迭代初始值，并且有P₀=(u,v,u_x,v_x,u_y,v_y,0)^T（亮度补偿参数w₀设为0）。

3、实验

为验证本文方法的有效性，用刚体转动实验来验证方法的准确性。首先用Peng Zhou等0提出的散斑模拟方法生成一张散斑图作为变形前图像。模拟散斑图如图5所示，具体参量为：图像大小为256pixel×256pixel，散斑颗粒大小为4pixel，散斑颗粒数目为4000。将图5绕中心顺时针旋转30°后的图像作为变形后图像，如图6所示。此时由于存在较大刚体转动角度，转动后的图像子区与原来的参考图像子区相似程度较低，通过整像素相关搜索来搜索变形后图像子区位置的方法失效，文献0认为转角大于7°时整像素相关搜索就完全失效。

由于是转动实验，为便于观察实验结果，在图5中只选取5个点（图5中的红色点）来进行计算。图像计算子区大小为61pixel×61pixel，SIFT算法阈值为0.6。对于第一个计算点，基于SIFT算法的变形初值估计和最终收敛的变形真实值分别如表1所示，表1显示两者非常接近。利用SIFT算法提供的变形初值估计，N-R方法经过几次迭代（实验用的5个点都用了4次迭代）即可迅速收敛。图6中绿色十字叉为变形前点所在位置，红色十字叉为变形后所在位置，从位置的变化可以看出，实验所测结果与预期吻合，表明了本文方法的有效性。

表1第一个计算点的变形初值与真值

Table 1 Initial guess and actual value of deformation for the first calculation point

4、结论

（1）利用SIFT算子对图像缩放、旋转甚至仿射变换保持不变的特性，通过SIFT算法匹配得到变形前后图像计算子区的对应点，并利用最小二乘法可获得变形初值估计作为N-R方法的迭代初值。

（2）该方法在被测物面出现包括较大刚体转动或大变形时仍能进行准确可靠的变形测量。

（3）该方法较之于先前已有的方法具备自动求取迭代初值的作用。

Claims (4)

1.一种数字图像变形初值的计算方法，其特征在于：包括如下步骤：

（1）获取被测物体表面变形前的数字图像及表面变形后的数字图像，将表面变形前的数字图像作为“参考图像”，表面变形后的图像作为“目标图像”；首先，将参考图像及目标图像均置于同一个平面坐标系内，在参考图像上，任意选取一个像素点作为网格节点，并以该网格节点为中心点，以2M+1为边长划分正方形区域，作为参考图像子区，其中M为像素点数，其次，通过特征点提取及匹配算法从目标图像中自动提取出与参考图像子区对应的像素点，以此得到一一对应的多组像素点，所述参考图像中与目标图像中任一对应的像素点之间的变换关系用“形函数”表示；

（2）根据一一对应像素点及参考图像子区中心点坐标计算确定中间计算量，所述中间计算量具体包括：参考图像子区中与目标图像子区中在同一坐标系下所有像素点的坐标值及参考图像子区中所有像素点的坐标相对于中心的变化量；

（3）根据步骤（2）中得到的中间计算量及步骤（1）中得到的一一对应的多组像素点利用最小二乘法求解所述“形函数”的参数值，将得到的“形函数”的参数值作为所述数字图像的变形初值。

2.根据权利要求1所述的数字图像变形初值的计算方法，其特征在于：所述“形函数”为：

x′=x+u+u_xΔx+u_yΔy

y′=y+v+v_xΔx+v_yΔy

其中，x′,y′为目标图像子区中除中心点外任一像素点的坐标，x,y为参考图像子区中与x′,y′对应的像素点的坐标，Δx,Δy分别为为参考图像子区中的像素点x,y到参考图像子区中心在横纵坐标轴方向的距离，u,v是参考图像子区中心到目标图像子区中心在横纵坐标轴方向的位移，u_x，u_y，v_x，v_y分别为参考图像子区中到目标图像子区中对应点的位移梯度。

3.根据权利要求1所述的数字图像变形初值的计算方法，其特征在于：所述M为大于等于20且小于等于50的整数。

4.根据权利要求1所述的数字图像变形初值的计算方法，其特征在于：所述“形函数”的参数值是通过将中间计算量代入一阶或二阶“形函数”得到的。

Images