数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法IC-GN的数值解计算

---

前言

由于本人近期正在展开数字图像相关技术用于测量材料形变方向的研究，其中需要对别人现有算法的复现和调研，尽管其中很多算法都已经非常成熟，但对于初学者而言即使明白其中的原理，无法上手实践和操作的话，依然无法能够将其完全的应用起来或者在上面进行创新，我希望能将自己作为一个初学者复现他人代码和学习该原理的过程记录下来，方便每一个涉足该领域的人能更快应用这些知识。

本文的论述基础建立在我的前一篇文章数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法（FA-GN与IC-GN）。推荐先通过这篇帖子了解DIC非线性迭代优化方法后，再阅读本文。

数字图像相关专栏目录：

1. Matlab实现二维数字图像相关（2D Digital Image Correlation, 2D-DIC）【ADIC2D代码复现及原理介绍】
2. 数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法（FA-GN与IC-GN）
3. 数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法IC-GN的数值解计算
4. 用MATLAB绘制随机散斑图案【源码+正确的椭圆旋转公式】

由于本人的论文已经完工，如果本文对你有所帮助，欢迎届时关注或引用我的论文。

---

内容回顾

为保证阅读的通畅，我会将上一篇帖子中IC-GN的算法逻辑在本章节进行简要的回顾与展示，便于在后面的公式推导中方便随时回头查看。

形函数相关参数说明：

1. 0阶形函数：${\mathbf{W}}^{SF0}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,{\mathbf{P}}^{SF0}\right)=\left[\begin{array}{ccc}1& 0& u\\ 0& 1& v\end{array}\right]\left[\begin{array}{c}\Delta x_i\\ \Delta y_i\\ 1\end{array}\right]$
2. 1阶形函数：${\mathbf{W}}^{SF1}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,{\mathbf{P}}^{SF1}\right)=\left[\begin{array}{ccc}1+u_x& u_y& u\\ v_x& 1+v_y& v\end{array}\right]\left[\begin{array}{c}\Delta x_i\\ \Delta y_i\\ 1\end{array}\right]$
3. 2阶形函数：${\mathbf{W}}^{SF2}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,{\mathbf{P}}^{SF2}\right)=\left[\begin{array}{cccccc}\frac{1}{2}{u}_{xx}& u_{xy}& \frac{1}{2}{u}_{yy}& 1+u_x& u_y& u\\ \frac{1}{2}{v}_{xx}& v_{xy}& \frac{1}{2}{v}_{yy}& v_x& 1+v_y& v\end{array}\right]\left[\begin{array}{c}\Delta x_i^2\\ \Delta x_i\Delta y_i\\ \Delta y_i^2\\ \Delta x_i\\ \Delta y_i\\ 1\end{array}\right]$

在上一文中我们详细推导了数字图像相关中的非线性迭代优化的算术逻辑：通过在参考图像上设置好子区$f_i=F\left(x^o+\Delta x_i\right)$ ，基于给定的一组形函数参数初值$P_{rc}^0$及$P_{rr}^0$，对目标函数（相关标准）不断迭代优化从而得到一个在$P_{rr}=0$的情况下的最优解$P_{rc}^{\ast }$，从而得到一组在形变图像上与参考子区最佳匹配的形变子区$g_i^{\ast }=G\left(x^o+\mathbf{W}\left(\Delta {\mathbf{x}}_i,{\mathbf{P}}_{\mathbf{r}\mathbf{c}}^{\ast }\right)\right)$，最终实现参考图像像素点与形变图像像素点的匹配。而我们针对优化$P_{rc}$去使用的高斯牛顿法就被叫做FA-GN，而针对优化$P_{rr}$去使用的高斯牛顿法就被叫做IC-GN。

采用一阶形函数${\mathbf{W}}^{SF1}$及 零均值归一化最小距离平方标准（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD） 来构建我们的数学模型

零均值归一化最小距离平方标准（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD）
$$C_{\text{ZNSSD }}=\sum _{i=1}^I{\left[\frac{f_i-\bar{f}}{\tilde{f}}-\frac{g_i-\bar{g}}{\tilde{g}}\right]}^2\in [0,4]$$其中$\bar{f}=\frac{{\sum }_i^I{f}_i}{I}$、$\bar{g}=\frac{{\sum }_i^I{f}_g}{I}$分别代表参考图像和形变图像的灰度值均值，$\tilde{f}=\sqrt{{\sum }_{i=1}^I{\left(f_i-\bar{f}\right)}^2}$、$\tilde{g}=\sqrt{{\sum }_{i=1}^I{\left(g_i-\bar{g}\right)}^2}$为子区的归一化函数【相当于没有除以总数的方差，表征了子区灰度值的离散程度】。$C_{\text{ZNSSD }}$数值越小，相关性越好。

设上一次迭代的到的形函数参数为${\mathbf{P}}^{\mathbf{o}\mathbf{l}\mathbf{d}}$，由于这一次是针对$P_{rr}$来做，故求解的增量$\Delta \mathbf{P}$即为$P_{rr}$，其每次的起点都是原参考子区，即形函数参数零点，如此我们的目标函数为$C_{\text{ZNSSD }}(\mathbf{\Delta }\mathbf{P})$，即
$$C_{\text{ZNSSD }}(\mathbf{\Delta }\mathbf{P})=\sum _{i=1}^I{\left[\frac{F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,\Delta \mathbf{P}\right)\right)-\bar{f}}{\sqrt{{\sum }_{i=1}^I{\left(F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,\Delta \mathbf{P}\right)\right)-\bar{f}\right)}^2}}-\frac{G\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,P^{old}\right)\right)-\bar{g}}{\sqrt{{\sum }_{i=1}^I{\left(G\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,P^{old}\right)\right)-\bar{g}\right)}^2}}\right]}^2$$根据IC-GN算法可以得到在该算法中目标函数在0处的梯度为：
$$\begin{array}{rl}\nabla C_{\text{ZNSSD }}(0)& =\frac{\mathrm{d}{C}_{\text{ZNSSD }}(0)}{\mathrm{d}\Delta \mathbf{P}}\\ & \approx \frac{2}{\sqrt{{\sum }_{i=1}^I{\left(F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)-\bar{f}\right)}^2}}\sum \left[\left[\frac{F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)-\bar{f}}{\sqrt{{\sum }_{i=1}^I{\left(F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)-\bar{f}\right)}^2}}-\frac{G\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,P^{old}\right)\right)-\bar{g}}{\sqrt{{\sum }_{i=1}^I{\left(G\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,P^{old}\right)\right)-\bar{g}\right)}^2}}\right]\left[\frac{\mathrm{d}F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)}{\mathrm{d}\Delta \mathbf{P}}\right]\right]\end{array}$$
在IC-GN中，所使用的目标函数在0处的Hessian为：
$$\nabla \nabla C_{\text{ZNSSD }}(0)\approx \frac{2}{{\sum }_{i=1}^I{\left(F\left(({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)-\bar{f}\right)}^2}\sum \left[\frac{\mathrm{d}F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)}{\mathrm{d}\Delta \mathbf{P}}\right]{\left[\frac{\mathrm{d}F\left({\mathbf{x}}^o+\mathbf{W}\left(\Delta {\mathbf{x}}^i{\mathbf{x}}^0,0\right)\right)}{\mathrm{d}\Delta \mathbf{P}}\right]}^T$$
根据得到的梯度和Hessian，我们便可以利用$\Delta \mathbf{P}=-\nabla C_{\text{ZNSSD }}(0)\ast \nabla \nabla C_{\text{ZNSSD }}(0{)}^{-1}$得出增量$\Delta \mathbf{P}$之后利用如下的更新式完成对形函数参数的迭代更新，设由形函数参数$\mathbf{P}$所形成的系数矩阵的增广矩阵为$\mathbf{w}(\mathbf{P})$，则有更新式为
$$\mathbf{w}({\mathbf{P}}^{\mathbf{n}\mathbf{e}\mathbf{w}})=\mathbf{w}({\mathbf{P}}^{\mathbf{o}\mathbf{l}\mathbf{d}})\mathbf{w}(\Delta \mathbf{P}{)}^{-1}$$

关于不同阶数形函数参数对应的増广系数矩阵，建议参考Gao, Y.; Cheng, T.; Su, Y.; Xu, X.; Zhang, Y.; Zhang, Q. High-efficiency and high-accuracy digital image correlation for three-dimensional measurement. Opt. Lasers Eng. 2015, 65, 73–80.

一.IC-GN中增量$\Delta \mathbf{P}$的数值解

从内容回顾中，可以知道只要得到增量$\Delta \mathbf{P}$便可以完成对于迭代的计算，但这里的$\Delta \mathbf{P}=-\nabla C_{\text{ZNSSD }}(0)\ast \nabla \nabla C_{\text{ZNSSD }}(0{)}^{-1}$该怎么去书写其数值解，成为了本文的目标。

利用ZNSSD描述中的符号$\tilde{f}，\tilde{g}$对IC-GN方法中的梯度和Hessian进行简化，可以得到
$$\begin{array}{rl}\nabla C_{\text{ZNSSD }}(0)& =\frac{2}{\tilde{f}}\sum \left[\left[\frac{f_i-\bar{f}}{\tilde{f}}-\frac{g_i-\bar{g}}{\tilde{g}}\right]\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}\right]\\ & =\frac{2}{{\tilde{f}}^2}\sum \left[\left[f_i-\bar{f}-\frac{\tilde{f}}{\tilde{g}}\left(g_i-\bar{g}\right)\right]\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}\right]\end{array}$$

$$\nabla \nabla C_{\text{ZNSSD }}(0)=\frac{2}{{\tilde{f}}^2}\sum \left[\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}\cdot {\left(\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}\right)}^T\right]$$根据$\Delta \mathbf{P}$的公式，则Hessian和梯度中的系数$\frac{2}{{\tilde{f}}^2}$被消去，仅需额外计算$\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}$即可得到增量$\Delta \mathbf{P}$的数值结果。利用链式法则有：
$$\frac{\mathrm{d}{f}_i}{\mathrm{d}\Delta \mathbf{P}}=\frac{\mathrm{d}{f}_i}{\mathrm{d}{\mathbf{x}}^i}\cdot \frac{\mathrm{d}{\mathbf{x}}^i}{\mathrm{d}\Delta \mathbf{P}}$$由于${\mathbf{x}}^i$代表了子区中其他像素点的坐标，因此第一项即为图像的梯度，在MATLAB中可以用imgradientxy()函数计算得到。第二项则是${\mathbf{x}}^i$对$\mathbf{P}$中各项求偏导的结果，根据形函数表达式，以一阶形函数为例，可以得到如下的表达式
$$\{\begin{array}{c}{x}_i=x_0+\Delta u+\left(1+\Delta u_x\right)\cdot \Delta x_i+\Delta u_y\cdot \Delta y_i\\ \\ y_i=y_0+\Delta v+\left(1+\Delta v_y\right)\cdot \Delta y_i+\Delta v_x\cdot \Delta x_i\end{array}$$对其各项求偏导即可得到第二项，即$\frac{\mathrm{d}{\mathbf{x}}^i}{\mathrm{d}\Delta \mathbf{P}}$的数值结果，经过整理有其关于不同阶数形函数的数值表达¹为：

1. 0阶形函数：$$\frac{\partial {\mathbf{x}}^i}{\partial \Delta {\mathbf{P}}^{SF0}}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

2. 1阶形函数：$$\frac{\partial {\mathbf{x}}^i}{\partial \Delta {\mathbf{P}}^{SF1}}=\left[\begin{array}{cccccc}1& \Delta x_i& \Delta y_i& 0& 0& 0\\ 0& 0& 0& 1& \Delta x_i& \Delta y_i\end{array}\right]$$

3. 2阶形函数：
   $$\frac{\partial {\mathbf{x}}^i}{\partial \Delta {\mathbf{P}}^{SF2}}=\left[\begin{array}{cccccccccccc}1& \Delta x_i& \Delta y_i& \Delta x_i^2/2& \Delta x_i\Delta y_i& \Delta y_i^2/2& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& \Delta x_i& \Delta y_i& \Delta x_i^2/2& \Delta x_i\Delta y_i& \Delta y_i^2/2\end{array}\right]$$

因此要完成IC-GN非线性迭代优化数值计算，需要引入的数据有

• 图像的梯度
• 子区像素点相对于中心点的位移量

利用这些数据信息配合上算法，即可完成非线性迭代的数值计算了。同理FA-GN也是如此，利用类似的方法也可以计算得到FA-GN的数值解结果

二.写在最后

这部分内容主要是作为辅助，方便大家能更快的通过MATLAB来搭建完整的数字图像相关匹配算法，里面的公式也是我亲自推导验证的，如果存在错误，还希望您可以及时指正，以确保所分享知识的准确性。很多知识的数学原理都是看上去很复杂，但是细细推导也就是高数里很普通的内容，希望能通过这样的方式让自己的知识活学活用起来。

参考引用

---

1. Lu H, Cary P D. Deformation measurements by digital image correlation: implementation of a second-order displacement gradient[J]. Experimental mechanics, 2000, 40(4): 393-400. ↩︎