前言
内容回顾
一. 非线性优化数学模型
二. 前向累加高斯-牛顿法——FA-GN（Forward Additive Gauss-Newton method）
三. 逆合成高斯-牛顿法——IC-GN（Inverse compositional Gauss-Newton method）
- 1.非线性优化数学模型变形
- 2.数学推导
四.写在最后
参考引用

前言

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

本文的论述基础建立在我的前一篇文章Matlab实现二维数字图像相关（2D Digital Image Correlation, 2D-DIC）【ADIC2D代码复现及原理介绍】。推荐先通过这篇帖子了解DIC整体模型及运算后，再阅读本文。

数字图像相关专栏目录：

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

本文的算术推导过程主要借鉴：http://www.ncorr.com/index.php/dic-algorithms

内容回顾

为保证阅读的通畅，我会将上一篇帖子中的模型在这一节中进行简单的回顾与展示，便于在后面的公式推导中方便随时回头查看。在上一文中，我们可以将数字图像相关总结归纳为：通过在参考图像上设置好子区fi=F(xo+Δxi)f_{i}=F\left(x^{o}+\Delta x_{i}\right)fi=F(xo+Δxi) ，基于给定的一组形函数参数初值P0P_{0}P0，对目标函数（相关标准）不断迭代优化从而得到一组最优解P∗P^{*}P∗，从而得到一组在形变图像上与参考子区最佳匹配的形变子区gi∗=G(xo+W(Δxi,P∗))g_{i}^{*}=G\left(x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}_{i}, \boldsymbol{P^{*}}\right)\right)gi∗=G(xo+W(Δxi,P∗))，最终实现参考图像像素点与形变图像像素点的匹配。具体模型示意如图所示。

本文会用到的各参数说明：

参考子区中心点：x0=[x0y0]T\boldsymbol{x}^{0}=\left[\begin{array}{ll}x_{0} & y_{0}\end{array}\right]^{T}x0=[x0y0]T
参考子区像素点：xi=Δxix0+xo=[xiyi]=[ΔxiΔyi]+[x0y0]\boldsymbol{x}^{i}=\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}+\boldsymbol{x}^{o}=\left[\begin{array}{l} x_{i} \\ y_{i} \end{array}\right]=\left[\begin{array}{l} \Delta x_{i} \\ \Delta y_{i} \end{array}\right]+\left[\begin{array}{l} x^{0} \\ y^{0} \end{array}\right]xi=Δxix0+xo=[xiyi]=[ΔxiΔyi]+[x0y0]
形变子区像素点：xi′=Δxi′x0+xo=[xi′yi′]=[Δxi′Δyi′]+[x0y0]\boldsymbol{x}^{i'}=\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}+\boldsymbol{x}^{o}=\left[\begin{array}{l} x_{i}' \\ y_{i}' \end{array}\right]=\left[\begin{array}{l} \Delta x_{i}' \\ \Delta y_{i}' \end{array}\right]+\left[\begin{array}{l} x^{0} \\ y^{0} \end{array}\right]xi′=Δxi′x0+xo=[xi′yi′]=[Δxi′Δyi′]+[x0y0]
形变子区像素点偏移量：Δxi′x0=W(Δxix0,P)\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}=\boldsymbol W(\Delta \boldsymbol{x}^{i}\boldsymbol{x}^{0},P)Δxi′x0=W(Δxix0,P)
参考子区：fi=F(xo+Δxix0)=F(xo+W(Δxix0,0))f_{i}=F\left(\boldsymbol x^{o}+\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}\right)=F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)fi=F(xo+Δxix0)=F(xo+W(Δxix0,0))
形变子区：gi=G(xo+W(Δxix0,P))g_{i}=G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P}\right)\right)gi=G(xo+W(Δxix0,P))

形函数相关参数说明：

0阶形函数：WSF0(Δxix0,PSF0)=[10u01v][ΔxiΔyi1]\boldsymbol{W}^{S F 0}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P}^{S F 0}\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]WSF0(Δxix0,PSF0)=[1001uv]⎣⎡ΔxiΔyi1⎦⎤
1阶形函数：WSF1(Δxix0,PSF1)=[1+uxuyuvx1+vyv][ΔxiΔyi1]\boldsymbol{W}^{S F 1}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P}^{S F 1}\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]WSF1(Δxix0,PSF1)=[1+uxvxuy1+vyuv]⎣⎡ΔxiΔyi1⎦⎤
2阶形函数：WSF2(Δxix0,PSF2)=[12uxxuxy12uyy1+uxuyu12vxxvxy12vyyvx1+vyv][Δxi2ΔxiΔyiΔyi2ΔxiΔyi1]\boldsymbol{W}^{S F 2}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P}^{S F 2}\right)=\left[\begin{array}{cccccc} \frac{1}{2} u_{x x} & u_{x y} & \frac{1}{2} u_{y y} & 1+u_{x} & u_{y} & u \\ \frac{1}{2} v_{x x} & v_{x y} & \frac{1}{2} v_{y y} & 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]WSF2(Δxix0,PSF2)=[21uxx21vxxuxyvxy21uyy21vyy1+uxvxuy1+vyuv]⎣⎢⎢⎢⎢⎢⎢⎡Δxi2ΔxiΔyiΔyi2ΔxiΔyi1⎦⎥⎥⎥⎥⎥⎥⎤

一. 非线性优化数学模型

数字图像相关中的非线性优化示意如图所示：¹

根据上面的内容回顾，我们可以知道，数字图像相关最重要的一点就是基于一个函数（相关标准）迭代求解其最优解（P∗P^{*}P∗）。而相关标准本身是一个非线性函数，其最优解的位置就位于其极大值或极小值处【由选择的相关标准决定】，而要想找到这个极值点位置，就需要用数值分析中的非线性方程组的数值解法。这些方法中，最为常用的即牛顿迭代法（Newton’s method，也叫做牛顿-拉夫逊法，Newton-Raphson method）。下面，我们利用这种方法来思考和推导我们这个问题：

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

零均值归一化最小距离平方标准（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD）
CZNSSD =∑i=1I[fi−fˉf~−gi−gˉg~]2∈[0,4]C_{\text {ZNSSD }}=\sum_{i=1}^{I}\left[\frac{f_{i}-\bar{f}}{\widetilde{f}}-\frac{g_{i}-\bar{g}}{\widetilde{g}}\right]^{2}\in [0,4]CZNSSD =i=1∑I[ffi−fˉ−ggi−gˉ]2∈[0,4]其中fˉ=∑iIfiI\bar{f}=\frac{\sum_{i}^{I} f_{i}}{I}fˉ=I∑iIfi、gˉ=∑iIfgI\bar{g}=\frac{\sum_{i}^{I} f_{g}}{I}gˉ=I∑iIfg分别代表参考图像和形变图像的灰度值均值，f~=∑i=1I(fi−fˉ)2\widetilde{f}=\sqrt{\sum_{i=1}^{I}\left(f_{i}-\bar{f}\right)^{2}}f=∑i=1I(fi−fˉ)2、g~=∑i=1I(gi−gˉ)2\widetilde{g}=\sqrt{\sum_{i=1}^{I}\left(g_{i}-\bar{g}\right)^{2}}g=∑i=1I(gi−gˉ)2为子区的归一化函数【相当于没有除以总数的方差，表征了子区灰度值的离散程度】。CZNSSD C_{\text {ZNSSD }}CZNSSD 数值越小，相关性越好。

已知：参考子区fi=F(xo+Δxix0)f_{i}=F\left(\boldsymbol x^{o}+\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}\right)fi=F(xo+Δxix0)、形函数参数初值P0\boldsymbol P^{0}P0
目标：求解最优形函数参数值P∗\boldsymbol{P^{*}}P∗

根据已知条件，我们可以推导得到此时的形变子区为gi=G(xo+W(Δxix0,P0))g_{i}=G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P^{0}}\right)\right)gi=G(xo+W(Δxix0,P0))，根据ZNSSD的公式，我们可以计算出这一数值，记为CZNSSD (P0)C_{\text {ZNSSD }}(\boldsymbol{P^{0}})CZNSSD (P0)。由于我们目的是要找到P∗\boldsymbol{P^{*}}P∗，在这个数学模型下即CZNSSD (P∗)⩽CZNSSD (P0)C_{\text {ZNSSD }}(\boldsymbol{P^{*}})\leqslant C_{\text {ZNSSD }}(\boldsymbol{P^{0}})CZNSSD (P∗)⩽CZNSSD (P0)。因此我们要给P0\boldsymbol P^{0}P0一个增量ΔP\Delta \boldsymbol PΔP，使得CZNSSD (P0+ΔP)⩽CZNSSD (P0)C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P})\leqslant C_{\text {ZNSSD }}(\boldsymbol{P^{0}})CZNSSD (P0+ΔP)⩽CZNSSD (P0)。

那现在的问题就是求解出这个ΔP\Delta \boldsymbol PΔP。首先对CZNSSD (P0+ΔP)C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P})CZNSSD (P0+ΔP)进行二阶泰勒展开有
CZNSSD (P0+ΔP)=CZNSSD (P0)+∇CZNSSD (P0)TΔP+12ΔPT∇∇CZNSSD (P0)ΔPC_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P})=C_{\text {ZNSSD }}(\boldsymbol{P^{0}})+\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})^{T}\Delta \boldsymbol P+\frac{1}{2}\Delta \boldsymbol P^{T}\nabla \nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})\Delta \boldsymbol P CZNSSD (P0+ΔP)=CZNSSD (P0)+∇CZNSSD (P0)TΔP+21ΔPT∇∇CZNSSD (P0)ΔP由于是求解极大（极小值），即CZNSSD (P0+ΔP)C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P})CZNSSD (P0+ΔP)关于ΔP\Delta \boldsymbol PΔP的导数为0，我们对上式求导有
dCZNSSD (P0+ΔP)dΔP≈∇CZNSSD (P0)+∇∇CZNSSD (P0)ΔP=0\frac{\mathrm{d} C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P})}{\mathrm{d}\Delta \boldsymbol P}\approx\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})+\nabla \nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})\Delta \boldsymbol P=0dΔPdCZNSSD (P0+ΔP)≈∇CZNSSD (P0)+∇∇CZNSSD (P0)ΔP=0

书本里展示的牛顿迭代法是在求解非线性方程组的零点，而现在要求极值点，即要求的是方程组导数（梯度）的零点，所以需要变通下思路。【网络上那些什么二阶牛顿拉夫逊法blabla的其实都是扯犊子，根本没这种说法】

此时，我们即可得到ΔP\Delta \boldsymbol PΔP的表达式为
ΔP=−∇CZNSSD (P0)∗∇∇CZNSSD (P0)−1\Delta \boldsymbol P = -\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})*\nabla \nabla C_{\text {ZNSSD }}(\boldsymbol{P^{0}})^{-1} ΔP=−∇CZNSSD (P0)∗∇∇CZNSSD (P0)−1这样一来，我们便可以更新得到一个逼近于P∗\boldsymbol{P^{*}}P∗的形函数参数Pnew=P0+ΔP\boldsymbol{P^{new}}=\boldsymbol{P^{0}+\Delta \boldsymbol P}Pnew=P0+ΔP，不断像这样如法炮制，迭代优化，直到计算得到的ΔP\Delta \boldsymbol PΔP小于停机准则【人为设置】，即可认为此时结果收敛，P∗≈Pnew\boldsymbol{P^{*}}\approx\boldsymbol{P^{new}}P∗≈Pnew。这种方法和思路就是标准的牛顿-拉夫逊法。而接下来涉及到的高斯牛顿法就是在这一基础上对其运算or思路进行了优化后的方法。

二. 前向累加高斯-牛顿法——FA-GN（Forward Additive Gauss-Newton method）

前向累加高斯-牛顿法（FA-GN，Forward Additive Gauss-Newton method） 其整体思路与上文的牛顿-拉夫逊法基本一致，只是做了一个运算上的简化，为了体现出这一简化的过程，这里将上面的公式细化，设上一次迭代的到的形函数参数为Pold\boldsymbol{P^{old}}Pold，我们有本次迭代的的目标函数CZNSSD (Pold+ΔP)C_{\text {ZNSSD }}(\boldsymbol{P^{old}+\Delta \boldsymbol P})CZNSSD (Pold+ΔP)，即
CZNSSD (Pold+ΔP)=∑i=1I[F(xo+W(Δxix0,0))−fˉ∑i=1I(F(xo+W(Δxix0,0))−fˉ)2−G(xo+W(Δxix0,Pold+ΔP))−gˉ∑i=1I(G(xo+W(Δxix0,Pold+ΔP))−gˉ)2]2C_{\text {ZNSSD }}(\boldsymbol{P^{old}+\Delta \boldsymbol P})=\sum_{i=1}^{I}\left[\frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}+\Delta \boldsymbol P\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}+\Delta \boldsymbol P\right)\right)-\bar{g}\right)^{2}}}\right]^{2} CZNSSD (Pold+ΔP)=i=1∑I⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2F(xo+W(Δxix0,0))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold+ΔP))−gˉ)2G(xo+W(Δxix0,Pold+ΔP))−gˉ⎦⎤2直接对他进行求导是非常复杂的，所以需要通过两个假设将其进行简化，即
d(gˉ)dΔP≈0\frac{\mathrm{d}\left ( \bar{g} \right )}{\mathrm{d} \Delta \boldsymbol P}\approx 0dΔPd(gˉ)≈0d(∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2)dΔP≈0\frac{\mathrm{d}\left ( \sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}} \right )}{\mathrm{d} \Delta \boldsymbol P}\approx 0 dΔPd(∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2)≈0数学公式看起来比较抽象，这里做下解释，第一个假设的实际意义可以理解为子区图像的灰度值均值不会因为子区形状变化产生改变，即图像上各位值的平均灰度值变化不剧烈【要求图像光照均匀，亮度一致】；第二个假设的实际意义可以理解为，从统计特性上来看，子区的灰度值离散程度不会因为子区形状变化产生改变，即图像上的散斑分布尽可能随机【要求图像上的散斑喷涂不要的出现一大片白或者一大片黑的情况】。这两个假设在满足拍摄要求的情况下是合理的，这样一来的话，进行求导时gˉ\bar{g}gˉ和∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2就可以看做两个常数项来进行处理了。

通过上述假设，求解目标函数在Pold\boldsymbol{P^{old}}Pold的梯度为：
∇CZNSSD (Pold)=dCZNSSD (Pold)dΔP≈−2∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2∑[[F(xo+W(Δxix0,0))−fˉ∑i=1I(F(xo+W(Δxix0,0))−fˉ)2−G(xo+W(Δxix0,Pold))−gˉ∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2][dG(xo+W(Δxix0,Pold))dΔP]]\begin{aligned} \nabla C_{\text {ZNSSD }}(\boldsymbol{P^{old}})&=\frac{\mathrm{d}C_{\text {ZNSSD }}(\boldsymbol{P^{old}})}{\mathrm{d} \Delta \boldsymbol P} \\ &\approx \frac{-2}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}}\sum \left [\left [ \frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}} \right ] \left [ \frac{\mathrm{d}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ]\right ] \end{aligned}∇CZNSSD (Pold)=dΔPdCZNSSD (Pold)≈∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2−2∑⎣⎡⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2F(xo+W(Δxix0,0))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2G(xo+W(Δxix0,Pold))−gˉ⎦⎤[dΔPdG(xo+W(Δxix0,Pold))]⎦⎤
再求解目标函数在Pold\boldsymbol{P^{old}}Pold的Hessian为：
∇∇CZNSSD (Pold)=d2CZNSSD (Pold)d(ΔP)2≈−2∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2{∑[−ddΔPG(xo+W(Δxix0,Pold))∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2][dG(xo+W(Δxix0,Pold))dΔP]T+∑[F(xo+W(Δxix0,0))−fˉ∑i=1I(F(xo+W(Δxix0,0))−fˉ)2−G(xo+W(Δxix0,Pold))−gˉ∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2][d2G(xo+W(Δxix0,Pold))d(ΔP)2]}\begin{aligned} \nabla\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{old}})&=\frac{\mathrm{d^{2}}C_{\text {ZNSSD }}(\boldsymbol{P^{old}})}{\mathrm{d} (\Delta \boldsymbol P)^{2}} \\ &\approx \frac{-2}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}}\left \{ \sum \left [ -\frac{\frac{\mathrm{d}}{\mathrm{d}\Delta \boldsymbol P}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}} \right ] \left [ \frac{\mathrm{d}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ]^{T} \right. \\\\ &\left. {+\sum \left [ \frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}} \right ] \left [ \frac{\mathrm{d^{2}}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\mathrm{d} (\Delta \boldsymbol P)^{2}} \right ]} \right \} \end{aligned} ∇∇CZNSSD (Pold)=d(ΔP)2d2CZNSSD (Pold)≈∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2−2⎩⎨⎧∑⎣⎡−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2dΔPdG(xo+W(Δxix0,Pold))⎦⎤[dΔPdG(xo+W(Δxix0,Pold))]T+∑⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2F(xo+W(Δxix0,0))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2G(xo+W(Δxix0,Pold))−gˉ⎦⎤[d(ΔP)2d2G(xo+W(Δxix0,Pold))]⎭⎬⎫如果直接使用这个Hessian来计算增量ΔP\Delta \boldsymbol PΔP的话就是传统的牛顿拉夫逊的做法了，但我们从式子上可以看到，这个Hessian的计算量是很大的，而高斯牛顿法的精髓就是将这个公式的进行了简化，将下面加号后的式子看做是0来简化计算。当然要满足这一假设需要满足两个条件：

Pold\boldsymbol{P^{old}}Pold接近于最优解P∗\boldsymbol{P^{*}}P∗【初值估计P0\boldsymbol{P^{0}}P0很重要】
灰度值的二阶导接近于0（这一要求主要取决于图像本身和插值平滑性，但一般来说图像灰度过渡都很均匀）【插值算法的平滑程度很重要】

在这两项条件满足的情况下，后面的式子相当于是两个小量相乘，等于一个更小的量可以近似为0，大大简化了计算目标函数Hessian的时间。

这种假设是有代价的，通常来说高斯牛顿法会比传统牛顿拉夫逊法有更多的迭代次数，但从单次迭代上来看，运算代价减少了很多。因此实际使用时需要进行权衡。通常来说子区越小，则用高斯牛顿法的效果会更好。当然，除了计算代价的优化外，高斯牛顿法还拥有更好的收敛性且更容易实现¹

这样一来我们便可以得到在FA-GN中,所使用的目标函数在Pold\boldsymbol{P^{old}}Pold的Hessian为：
∇∇CZNSSD (Pold)≈2∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2∑[dG(xo+W(Δxix0,Pold))dΔP][dG(xo+W(Δxix0,Pold))dΔP]T\nabla\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{old}}) \approx \frac{2}{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}} \sum \left [ \frac{\mathrm{d}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ] \left [ \frac{\mathrm{d}G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ]^{T} ∇∇CZNSSD (Pold)≈∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)22∑[dΔPdG(xo+W(Δxix0,Pold))][dΔPdG(xo+W(Δxix0,Pold))]T之后利用上文提到的计算公式ΔP=−∇CZNSSD (Pold)∗∇∇CZNSSD (Pold)−1\Delta \boldsymbol P = -\nabla C_{\text {ZNSSD }}(\boldsymbol{P^{old}})*\nabla \nabla C_{\text {ZNSSD }}(\boldsymbol{P^{old}})^{-1}ΔP=−∇CZNSSD (Pold)∗∇∇CZNSSD (Pold)−1得出ΔP\Delta \boldsymbol PΔP后按照下面的迭代公式更新即可
Pnew=Pold+ΔP\boldsymbol{P^{new}}=\boldsymbol{P^{old}+\Delta \boldsymbol P}Pnew=Pold+ΔP
FA-GN的整体流程图如下¹

这里总结一下FA-GN的整体算法逻辑，通过设定好的参考子区和原始形函数参数，得到与之对应的原始形变子区，利用高斯牛顿法计算得到一组相较于原始形变子区的增量，更新得到新的形函数参数，然后将这一组新的形函数参数作为下一次迭代的原始形函数参数，不断进行迭代，直到逼近最优形函数参数。总而言之，这种方法始终都是在改变形变子区的形状来找到与参考子区最佳匹配的时候。

三. 逆合成高斯-牛顿法——IC-GN（Inverse compositional Gauss-Newton method）

与FA-GN不同，逆合成高斯-牛顿法——IC-GN（Inverse compositional Gauss-Newton method）在整体思路上有很大的变化。从上文的结论中可以看到，FA-GN始终都是在计算相较于原始形变子区的增量，然后找到与参考子区匹配程度最好形变子区。而IC-GN不同，它得到原始形变子区后，计算的是相较于参考子区的增量，然后通过增量变形参考子区，找到与原始形变子区匹配最好的变形参考子区，从而更新得到新的形函数参数，将其作为下一次迭代的原始形函数参数不断迭代优化。为方便理解， IC-GN的整体流程图如下¹

1.非线性优化数学模型变形

如果你通过上面的文字和流程图已经基本理解了IC-GN的核心思想，那我建议你直接跳到下面的数学推导部分，如果你还是很蒙，那试试这一节能不能让你弄通顺一些。在本文的第一章中，我们在已知形函数参数初值的情况下，可以得到一对参考子区和形变子区的表达式，即

fi=F(xo+Δxix0)=F(xo+W(Δxix0,0))f_{i}=F\left(\boldsymbol x^{o}+\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}\right)=F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)fi=F(xo+Δxix0)=F(xo+W(Δxix0,0))
gi=G(xo+W(Δxix0,P0))g_{i}=G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P^{0}}\right)\right)gi=G(xo+W(Δxix0,P0))

然后利用牛顿拉夫逊法的思路可以计算得到增量ΔP\Delta \boldsymbol PΔP去更新形变子区，即在形变图像上找到一个与原参考子区更加匹配的形变子区。那我们是否能反过来，用同样的方法在参考图像上找到一个与形变子区更加匹配的参考子区呢？其答案是可以的，从下图¹可以看到这三者的一个关系

图中绿色区域即为原始参考子区，PrcP_{rc}Prc代表由原参考子区变换到形变子区的形函数参数，即上面公式中的P0\boldsymbol{P^{0}}P0；PrrP_{rr}Prr代表由原参考子区变换到新参考子区的形函数参数。而这个块新参考子区与形变子区的匹配程度无疑是更高的，这个新的参考子区可以表示为：
fi^=F(xo+W(Δxix0,Prr))\widehat{f_{i}}=F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P_{rr}\right)\right) fi=F(xo+W(Δxix0,Prr))当Prr=0P_{rr}=0Prr=0时，即代表着新参考子区与原参考子区完全重合。那这样一来，我们就有了一个新的目标，与FA-GN或者传统牛顿拉夫逊不同，他们是在单纯的找满足最大相关标准时候的Prc∗P_{rc}^{*}Prc∗，而直接求解这个值是很麻烦的，而对于PrrP_{rr}Prr其最优解就是Prr∗=0P_{rr}^{*}=0Prr∗=0，只要想办法让这个形函数参数逼近0即可。

这样一来我们就可以更新下我们对于数字图像相关的归纳：通过在参考图像上设置好子区fi=F(xo+Δxi)f_{i}=F\left(x^{o}+\Delta x_{i}\right)fi=F(xo+Δxi) ，基于给定的一组形函数参数初值Prc0P_{rc}^{0}Prc0及Prr0P_{rr}^{0}Prr0，对目标函数（相关标准）不断迭代优化从而得到一个在Prr=0P_{rr}=0Prr=0的情况下的最优解Prc∗P_{rc}^{*}Prc∗，从而得到一组在形变图像上与参考子区最佳匹配的形变子区gi∗=G(xo+W(Δxi,Prc∗))g_{i}^{*}=G\left(x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}_{i}, \boldsymbol{P_{rc}^{*}}\right)\right)gi∗=G(xo+W(Δxi,Prc∗))，最终实现参考图像像素点与形变图像像素点的匹配。而我们针对优化PrcP_{rc}Prc去使用的高斯牛顿法就被叫做FA-GN，而针对优化PrrP_{rr}Prr去使用的高斯牛顿法就被叫做IC-GN。

2.数学推导

同样的，为了展示IC-GN与其他算法的不同性以及其优点，这里将其具体运算展开来写，设上一次迭代的到的形函数参数为Pold\boldsymbol{P^{old}}Pold，由于这一次是针对PrrP_{rr}Prr来做，故求解的增量ΔP\Delta \boldsymbol PΔP即为PrrP_{rr}Prr，其每次的起点都是原参考子区，即形函数参数零点，如此我们的目标函数为CZNSSD (ΔP)C_{\text {ZNSSD }}(\boldsymbol{\Delta \boldsymbol P})CZNSSD (ΔP)，即
CZNSSD (ΔP)=∑i=1I[F(xo+W(Δxix0,ΔP))−fˉ∑i=1I(F(xo+W(Δxix0,ΔP))−fˉ)2−G(xo+W(Δxix0,Pold))−gˉ∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2]2C_{\text {ZNSSD }}(\boldsymbol{\Delta \boldsymbol P})=\sum_{i=1}^{I}\left[\frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \Delta \boldsymbol P\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \Delta \boldsymbol P\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}}\right]^{2} CZNSSD (ΔP)=i=1∑I⎣⎡∑i=1I(F(xo+W(Δxix0,ΔP))−fˉ)2F(xo+W(Δxix0,ΔP))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2G(xo+W(Δxix0,Pold))−gˉ⎦⎤2同样的，由于参考图像和形变图像基本是在同样场景下拍摄的，我们可以和FA-GN一样提出以下假设
d(fˉ)dΔP≈0\frac{\mathrm{d}\left ( \bar{f} \right )}{\mathrm{d} \Delta \boldsymbol P}\approx 0dΔPd(fˉ)≈0d(∑i=1I(F(xo+W(Δxix0,0))−fˉ)2)dΔP≈0\frac{\mathrm{d}\left ( \sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}} \right )}{\mathrm{d} \Delta \boldsymbol P}\approx 0 dΔPd(∑i=1I(F(xo+W(Δxix0,0))−fˉ)2)≈0通过上述假设，可以计算得到目标函数在0处的梯度为：
∇CZNSSD (0)=dCZNSSD (0)dΔP≈2∑i=1I(F(xo+W(Δxix0,0))−fˉ)2∑[[F(xo+W(Δxix0,0))−fˉ∑i=1I(F(xo+W(Δxix0,0))−fˉ)2−G(xo+W(Δxix0,Pold))−gˉ∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2][dF(xo+W(Δxix0,0))dΔP]]\begin{aligned} \nabla C_{\text {ZNSSD }}(\boldsymbol{0})&=\frac{\mathrm{d}C_{\text {ZNSSD }}(0)}{\mathrm{d} \Delta \boldsymbol P} \\ &\approx \frac{2}{ \sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}\sum \left [\left [ \frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}} \right ] \left [ \frac{\mathrm{d}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) }{\mathrm{d}\Delta \boldsymbol P} \right ]\right ] \end{aligned}∇CZNSSD (0)=dΔPdCZNSSD (0)≈∑i=1I(F(xo+W(Δxix0,0))−fˉ)22∑⎣⎡⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2F(xo+W(Δxix0,0))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2G(xo+W(Δxix0,Pold))−gˉ⎦⎤[dΔPdF(xo+W(Δxix0,0))]⎦⎤
再求解目标函数在0处的Hessian为：
∇∇CZNSSD (0)=d2CZNSSD (0)d(ΔP)2≈2∑i=1I(F(xo+W(Δxix0,0))−fˉ)2{∑[ddΔPF(xo+W(Δxix0,0))∑i=1I(F(xo+W(Δxix0,0))−fˉ)2][dF(xo+W(Δxix0,0))dΔP]T+∑[F(xo+W(Δxix0,0))−fˉ∑i=1I(F(xo+W(Δxix0,0))−fˉ)2−G(xo+W(Δxix0,Pold))−gˉ∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2][d2F(xo+W(Δxix0,0))d(ΔP)2]}\begin{aligned} \nabla\nabla C_{\text {ZNSSD }}(\boldsymbol{0})&=\frac{\mathrm{d^{2}}C_{\text {ZNSSD }}(0)}{\mathrm{d} (\Delta \boldsymbol P)^{2}} \\ &\approx \frac{2}{ \sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}\left \{ \sum \left [ \frac{\frac{\mathrm{d}}{\mathrm{d}\Delta \boldsymbol P}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}} \right ] \left [ \frac{\mathrm{d}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ]^{T} \right. \\\\ &\left. {+\sum \left [ \frac{F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right) -\bar{f}}{\sqrt{\sum_{i=1}^{I}\left(F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}}}-\frac{G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}}{\sqrt{\sum_{i=1}^{I}\left(G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, P^{old}\right)\right)-\bar{g}\right)^{2}}} \right ] \left [ \frac{\mathrm{d^{2}}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)}{\mathrm{d} (\Delta \boldsymbol P)^{2}} \right ]} \right \} \end{aligned} ∇∇CZNSSD (0)=d(ΔP)2d2CZNSSD (0)≈∑i=1I(F(xo+W(Δxix0,0))−fˉ)22⎩⎨⎧∑⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2dΔPdF(xo+W(Δxix0,0))⎦⎤[dΔPdF(xo+W(Δxix0,0))]T+∑⎣⎡∑i=1I(F(xo+W(Δxix0,0))−fˉ)2F(xo+W(Δxix0,0))−fˉ−∑i=1I(G(xo+W(Δxix0,Pold))−gˉ)2G(xo+W(Δxix0,Pold))−gˉ⎦⎤[d(ΔP)2d2F(xo+W(Δxix0,0))]⎭⎬⎫同样的，利用和之前FA-GN的方法，上式中下半部分可以当做是0来处理，这样的话就可以得到在在IC-GN中,所使用的目标函数在0处的Hessian为：
∇∇CZNSSD (0)≈2∑i=1I(F((xo+W(Δxix0,0))−fˉ)2∑[dF(xo+W(Δxix0,0))dΔP][dF(xo+W(Δxix0,0))dΔP]T\nabla\nabla C_{\text {ZNSSD }}(0) \approx \frac{2}{\sum_{i=1}^{I}\left(F\left((\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)-\bar{f}\right)^{2}} \sum \left [ \frac{\mathrm{d}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ] \left [ \frac{\mathrm{d}F\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, 0\right)\right)}{\mathrm{d}\Delta \boldsymbol P} \right ]^{T} ∇∇CZNSSD (0)≈∑i=1I(F((xo+W(Δxix0,0))−fˉ)22∑[dΔPdF(xo+W(Δxix0,0))][dΔPdF(xo+W(Δxix0,0))]T从这个式子可以看到，在IC-GN中迭代所需要使用的Hessian，参与其计算的都是一些已知量，不会因为迭代次数而改变，意味着可以在迭代前就将Hessian这一数值计算出来，大大加快了每次迭代的运算速度。在我的实测中，对一副512*512的图片在同等子区大小的情况下进行非线性优化计算，IC-GN的计算时间仅需要FA-GN的一半左右，而且最后收敛的结果是基本一致的！！

根据得到的梯度和Hessian，我们便可以利用ΔP=−∇CZNSSD (0)∗∇∇CZNSSD (0)−1\Delta \boldsymbol P = -\nabla C_{\text {ZNSSD }}(0)*\nabla \nabla C_{\text {ZNSSD }}(0)^{-1}ΔP=−∇CZNSSD (0)∗∇∇CZNSSD (0)−1得出增量ΔP\Delta \boldsymbol PΔP，接下来就是利用这个增量去更新参考子区与形变子区之间对应的形函数参数（Prc\boldsymbol{P_{rc}}Prc）。但要注意的是，此时我们得到的增量是相较于参考子区的（Prr\boldsymbol{P_{rr}}Prr），因此更新式并不是直接相加的关系！两者的关系为：
W(Δxix0,Pnew)=W(W(Δxix0,ΔP)−1,Pold)\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P^{new}}\right)=\boldsymbol{W}\left(\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \Delta\boldsymbol{P}\right)^{-1}, \boldsymbol{P^{old}}\right) W(Δxix0,Pnew)=W(W(Δxix0,ΔP)−1,Pold)并且为了计算上式，通常会将形函数中的系数矩阵増广为一个“方阵”（对角占优矩阵，因为ΔP\Delta \boldsymbol PΔP的运算需要使用Cholesky分解），设由形函数参数P\boldsymbol PP所形成的系数矩阵的增广矩阵为w(P)\boldsymbol w(\boldsymbol P)w(P)，则有更新式为
w(Pnew)=w(Pold)w(ΔP)−1\boldsymbol w(\boldsymbol {P^{new}})=\boldsymbol w(\boldsymbol {P^{old}})\boldsymbol w(\Delta \boldsymbol P)^{-1}w(Pnew)=w(Pold)w(Δ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.

四.写在最后

非线性优化这一块我是纯手算推导了四页纸才大致搞明白了一些，但依然因为自己的线性代数水平较差很多地方卡了很久，很多论文中都有不同的推导计算公式，但Ncorr中的这个推导是最为完整的同时表述上不会存在二义性，建议作为参考。上一篇帖子中所用到的那篇论文²，在非线性优化这里推导的时候，计算的实际是ΔP\Delta \boldsymbol PΔP的转置（PPP实际是一个列向量，但是他算法实现和推导中其实都是用的行向量），这样的二义性在很多论文中都存在，同样的参数名却有不同的表达式。最后作为小白，很感谢Ncorr作者Justin Blaber对于知识的分享，也希望各位业界中的大佬对于我文中理解上的偏差和错误进行指正。

参考引用

http://www.ncorr.com/index.php/dic-algorithms ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
Atkinson D, Becker T. A 117 Line 2D Digital Image Correlation Code Written in MATLAB[J]. Remote Sensing, 2020, 12(18): 2906. ↩︎

数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法（FA-GN与IC-GN）相关推荐

数字图像相关（Digital Image Correlation, DIC）中的非线性优化方法IC-GN的数值解计算
目录前言内容回顾一.IC-GN中增量 Δ P \Delta \boldsymbol P ΔP的数值解二.写在最后参考引用前言由于本人近期正在展开数字图像相关技术用于测量材料形变方向的研究 ...
java mysbatis select_java相关：详解Mybatis中的select方法
java相关:详解Mybatis中的select方法发布于 2020-7-3| 复制链接摘记: selectById方法根据id,查询记录 ```java public void updateRe ...
JDBC中的setObject方法时干什么的
JDBC中的setObject方法时干什么的 2013-05-27 14:45zyfysukhhj | 分类:JAVA相关 | 浏览1146次 JDBC中的setObject方法时干什么的分享到: ...
#DIC#数字图像相关
1.1DIC基本原理在实验中DIC特指一种种光学测量技术,⽤于在整个⼒学试验过程中测量试样表⾯上不断变化的全场⼆维或三维坐标.测量出的坐标场可⽤于进⼀步导出位移.应变.应变率.速度和曲率等感兴趣量( ...
数字图像相关DIC算法，Ubuntu16.04，Ncorr项目C++版本开源环境配置
数字图像相关DIC算法,Ubuntu16.04,NcorrC++版本开源代码环境配置流程. 本文介绍C++版本数字图像相关法DIC环境配置过程,配置了好几天,痛苦踩坑经历. DIC算法.资料.源码.实 ...
应变测量的传统方法与DIC（数字图像相关）方法
应变测量的传统方法与DIC(数字图像相关)方法 1. 传统应变片测量应变原理: 1.1. 应变的定义: 1.2 应变片的工作原理: 1.3 传统应变片测量方法与注意事项 1.4 应变片的布置与组桥 1 ...
二维数字图像相关算法软件Ncorr的使用心得
二维数字图像相关(2D Digital Image Correlation)是一种非接触式的光学测量方法,常应用于图像分析处理上,它可以根据变形前后的2张或多张图像,求解出规定区域近似的位移与应变情 ...
中科大何向南团队+快手App联合出品 KuaiRec | 快手首个稠密为99.6%的数据集 | 相关介绍、下载、处理、使用方法
文章目录 1. 数据集介绍 1.1 相关链接: 1.2 构建方法 1.3 代表性验证 1.4 相关实验 2. 数据集下载 2.1 big matrix 2.1 small matrix 2.3 ite ...
计算机综合课设交通运输相关,计算机在道路运输管理中的应用课程设计.doc
计算机在道路运输管理中的应用课程设计课程设计论文题目:课程名称:计算机在道路运输管理中的应用学院: 交通运输专业: 交通运输班级: 学生姓名: 学号: 指导教师: ====2010 ...
基于激光雷达的里程计及3D点云地图中的定位方法
本文转载自公众号@点云PCL,基于激光雷达的里程计及3D点云地图中的定位方法 :https://mp.weixin.qq.com/s/laA1YAPBCpqlzdGi0yb2cQ 论文:LOL: Li ...

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

目录

前言