---

前言

由于本人近期正在展開數字影像相關技術用于測量材料形變方向的研究，其中需要對別人現有演算法的復現和調研，盡管其中很多演算法都已經非常成熟，但對于初學者而言即使明白其中的原理，無法上手實踐和操作的話，依然無法能夠將其完全的應用起來或者在上面進行創新，我希望能將自己作為一個初學者復現他人代碼和學習該原理的程序記錄下來，方便每一個涉足該領域的人能更快應用這些知識，

本文的論述基礎建立在我的前一篇文章Matlab實作二維數字影像相關（2D Digital Image Correlation, 2D-DIC）【ADIC2D代碼復現及原理介紹】，推薦先通過這篇帖子了解DIC整體模型及運算后，再閱讀本文，

本文的算術推導程序主要借鑒：http://www.ncorr.com/index.php/dic-algorithms

---

內容回顧

為保證閱讀的通暢，我會將上一篇帖子中的模型在這一節中進行簡單的回顧與展示，便于在后面的公式推導中方便隨時回頭查看，在上一文中，我們可以將數字影像相關總結歸納為：通過在參考影像上設定好子區 f i = F ( x o + Δ x i ) f_{i}=F\left(x^{o}+\Delta x_{i}\right) fi?=F(xo+Δxi?) ，基于給定的一組形函式引數初值 P 0 P_{0} P0?，對目標函式（相關標準）不斷迭代優化從而得到一組最優解 P ? P^{*} P?，從而得到一組在形變影像上與參考子區最佳匹配的形變子區 g i ? = G ( x o + W ( Δ x i , P ? ) ) g_{i}^{*}=G\left(x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}_{i}, \boldsymbol{P^{*}}\right)\right) gi??=G(xo+W(Δxi?,P?))，最終實作參考影像像素點與形變影像像素點的匹配， 具體模型示意如圖所示，

本文會用到的各引數說明：

1. 參考子區中心點： x 0 = [ x 0 y 0 ] T \boldsymbol{x}^{0}=\left[\begin{array}{ll}x_{0} & y_{0}\end{array}\right]^{T} x0=[x0??y0??]T

2. 參考子區像素點： x i = Δ x i x 0 + x o = [ x i y i ] = [ Δ x i Δ y i ] + [ x 0 y 0 ] \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=[xi?yi??]=[Δxi?Δyi??]+[x0y0?]

3. 形變子區像素點： x i ′ = Δ x i ′ x 0 + x o = [ x i ′ y i ′ ] = [ Δ x i ′ Δ y i ′ ] + [ x 0 y 0 ] \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?]

4. 形變子區像素點偏移量： Δ x i ′ x 0 = W ( Δ x i x 0 , P ) \Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}=\boldsymbol W(\Delta \boldsymbol{x}^{i}\boldsymbol{x}^{0},P) Δxi′x0=W(Δxix0,P)

5. 參考子區： f i = F ( x o + Δ x i x 0 ) = F ( x o + W ( Δ x i x 0 , 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))

6. 形變子區： g i = G ( x o + W ( Δ x i x 0 , 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))

形函式相關引數說明：

1. 0階形函式： W S F 0 ( Δ x i x 0 , P S F 0 ) = [ 1 0 u 0 1 v ] [ Δ x i Δ y i 1 ] \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)=[10?01?uv?]???Δxi?Δyi?1????
2. 1階形函式： W S F 1 ( Δ x i x 0 , P S F 1 ) = [ 1 + u x u y u v x 1 + v y v ] [ Δ x i Δ y i 1 ] \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+ux?vx??uy?1+vy??uv?]???Δxi?Δyi?1????
3. 2階形函式： W S F 2 ( Δ x i x 0 , P S F 2 ) = [ 1 2 u x x u x y 1 2 u y y 1 + u x u y u 1 2 v x x v x y 1 2 v y y v x 1 + v y v ] [ Δ x i 2 Δ x i Δ y i Δ y i 2 Δ x i Δ y i 1 ] \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)=[21?uxx?21?vxx??uxy?vxy??21?uyy?21?vyy??1+ux?vx??uy?1+vy??uv?]?????????Δxi2?Δxi?Δyi?Δyi2?Δxi?Δyi?1??????????

一. 非線性優化數學模型

數字影像相關中的非線性優化示意如圖所示：¹

根據上面的內容回顧，我們可以知道，數字影像相關最重要的一點就是基于一個函式（相關標準）迭代求解其最優解（ P ? P^{*} P?），而相關標準本身是一個非線性函式，其最優解的位置就位于其極大值或極小值處【由選擇的相關標準決定】，而要想找到這個極值點位置，就需要用數值分析中的非線性方程組的數值解法，這些方法中，最為常用的即牛頓迭代法（Newton’s method，也叫做牛頓-拉夫遜法，Newton-Raphson method），下面，我們利用這種方法來思考和推導我們這個問題：

注：采用一階形函式 W S F 1 \boldsymbol{W}^{S F 1} WSF1及 零均值歸一化最小距離平方標準（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD） 來構建我們的數學模型

零均值歸一化最小距離平方標準（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD）
C ZNSSD = ∑ i = 1 I [ f i ? f ˉ f ~ ? g i ? 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?[f ?fi??fˉ???g ?gi??gˉ??]2∈[0,4]其中 f ˉ = ∑ i I f i I \bar{f}=\frac{\sum_{i}^{I} f_{i}}{I} fˉ?=I∑iI?fi??、 g ˉ = ∑ i I f g I \bar{g}=\frac{\sum_{i}^{I} f_{g}}{I} gˉ?=I∑iI?fg??分別代表參考影像和形變影像的灰度值均值， f ~ = ∑ i = 1 I ( f i ? f ˉ ) 2 \widetilde{f}=\sqrt{\sum_{i=1}^{I}\left(f_{i}-\bar{f}\right)^{2}} f ?=∑i=1I?(fi??fˉ?)2 ?、 g ~ = ∑ i = 1 I ( g i ? g ˉ ) 2 \widetilde{g}=\sqrt{\sum_{i=1}^{I}\left(g_{i}-\bar{g}\right)^{2}} g ?=∑i=1I?(gi??gˉ?)2 ?為子區的歸一化函式【相當于沒有除以總數的方差，表征了子區灰度值的離散程度】， C ZNSSD C_{\text {ZNSSD }} CZNSSD ?數值越小，相關性越好，

已知：參考子區 f i = F ( x o + Δ x i x 0 ) f_{i}=F\left(\boldsymbol x^{o}+\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}\right) fi?=F(xo+Δxix0)、形函式引數初值 P 0 \boldsymbol P^{0} P0
目標：求解最優形函式引數值 P ? \boldsymbol{P^{*}} P?

根據已知條件，我們可以推導得到此時的形變子區為 g i = G ( x o + W ( Δ x i x 0 , P 0 ) ) 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的公式，我們可以計算出這一數值，記為 C ZNSSD ( P 0 ) C_{\text {ZNSSD }}(\boldsymbol{P^{0}}) CZNSSD ?(P0)，由于我們目的是要找到 P ? \boldsymbol{P^{*}} P?，在這個數學模型下即 C ZNSSD ( P ? ) ? C ZNSSD ( P 0 ) C_{\text {ZNSSD }}(\boldsymbol{P^{*}})\leqslant C_{\text {ZNSSD }}(\boldsymbol{P^{0}}) CZNSSD ?(P?)?CZNSSD ?(P0)，因此我們要給 P 0 \boldsymbol P^{0} P0一個增量 Δ P \Delta \boldsymbol P ΔP，使得 C ZNSSD ( P 0 + Δ P ) ? C ZNSSD ( P 0 ) 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，首先對 C ZNSSD ( P 0 + Δ P ) C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P}) CZNSSD ?(P0+ΔP)進行二階泰勒展開有
C ZNSSD ( P 0 + Δ P ) = C ZNSSD ( P 0 ) + ? C ZNSSD ( P 0 ) T Δ P + 1 2 Δ P T ? ? C ZNSSD ( P 0 ) Δ P C_{\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由于是求解極大（極小值），即 C ZNSSD ( P 0 + Δ P ) C_{\text {ZNSSD }}(\boldsymbol{P^{0}+\Delta \boldsymbol P}) CZNSSD ?(P0+ΔP)關于 Δ P \Delta \boldsymbol P ΔP的導數為0，我們對上式求導有
d C ZNSSD ( P 0 + Δ P ) d Δ P ≈ ? C ZNSSD ( P 0 ) + ? ? C ZNSSD ( P 0 ) Δ 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=0 dΔPdCZNSSD ?(P0+ΔP)?≈?CZNSSD ?(P0)+??CZNSSD ?(P0)ΔP=0

書本里展示的牛頓迭代法是在求解非線性方程組的零點，而現在要求極值點，即要求的是方程組導數（梯度）的零點，所以需要變通下思路，【網路上那些什么二階牛頓拉夫遜法blabla的其實都是扯犢子，根本沒這種說法】

此時，我們即可得到 Δ P \Delta \boldsymbol P ΔP的運算式為
Δ P = ? ? C ZNSSD ( P 0 ) ? ? ? C ZNSSD ( P 0 ) ? 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?的形函式引數 P n e w = P 0 + Δ P \boldsymbol{P^{new}}=\boldsymbol{P^{0}+\Delta \boldsymbol P} Pnew=P0+ΔP，不斷像這樣如法炮制，迭代優化，直到計算得到的 Δ P \Delta \boldsymbol P ΔP小于停機準則【人為設定】，即可認為此時結果收斂， P ? ≈ P n e w \boldsymbol{P^{*}}\approx\boldsymbol{P^{new}} P?≈Pnew，這種方法和思路就是標準的牛頓-拉夫遜法，而接下來涉及到的高斯牛頓法就是在這一基礎上對其運算or思路進行了優化后的方法，

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

前向累加高斯-牛頓法（FA-GN，Forward Additive Gauss-Newton method） 其整體思路與上文的牛頓-拉夫遜法基本一致，只是做了一個運算上的簡化，為了體現出這一簡化的程序，這里將上面的公式細化，設上一次迭代的到的形函式引數為 P o l d \boldsymbol{P^{old}} Pold，我們有本次迭代的的目標函式 C ZNSSD ( P o l d + Δ P ) C_{\text {ZNSSD }}(\boldsymbol{P^{old}+\Delta \boldsymbol P}) CZNSSD ?(Pold+ΔP)，即
C ZNSSD ( P o l d + Δ P ) = ∑ i = 1 I [ F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d + Δ P ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d + Δ P ) ) ? g ˉ ) 2 ] 2 C_{\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ˉ?)2 ?F(xo+W(Δxix0,0))?fˉ???∑i=1I?(G(xo+W(Δxix0,Pold+ΔP))?gˉ?)2 ?G(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 0 dΔPd(gˉ?)?≈0 d ( ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? 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 = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? 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 ?就可以看做兩個常數項來進行處理了，

通過上述假設，求解目標函式在 P o l d \boldsymbol{P^{old}} Pold的梯度為：
? C ZNSSD ( P o l d ) = d C ZNSSD ( P o l d ) d Δ P ≈ ? 2 ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ∑ [ [ F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] [ d G ( x o + W ( Δ x i x 0 , P o l d ) ) 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ˉ?)2 ?F(xo+W(Δxix0,0))?fˉ???∑i=1I?(G(xo+W(Δxix0,Pold))?gˉ?)2 ?G(xo+W(Δxix0,Pold))?gˉ?????[dΔPdG(xo+W(Δxix0,Pold))?]????
再求解目標函式在 P o l d \boldsymbol{P^{old}} Pold的Hessian為：
? ? C ZNSSD ( P o l d ) = d 2 C ZNSSD ( P o l d ) d ( Δ P ) 2 ≈ ? 2 ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 { ∑ [ ? d d Δ P G ( x o + W ( Δ x i x 0 , P o l d ) ) ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] [ d G ( x o + W ( Δ x i x 0 , P o l d ) ) d Δ P ] T + ∑ [ F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] [ d 2 G ( x o + W ( Δ x i x 0 , P o l d ) ) 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ˉ?)2 ?dΔPd?G(xo+W(Δxix0,Pold))????[dΔPdG(xo+W(Δxix0,Pold))?]T+∑???∑i=1I?(F(xo+W(Δxix0,0))?fˉ?)2 ?F(xo+W(Δxix0,0))?fˉ???∑i=1I?(G(xo+W(Δxix0,Pold))?gˉ?)2 ?G(xo+W(Δxix0,Pold))?gˉ?????[d(ΔP)2d2G(xo+W(Δxix0,Pold))?]?????如果直接使用這個Hessian來計算增量 Δ P \Delta \boldsymbol P ΔP的話就是傳統的牛頓拉夫遜的做法了，但我們從式子上可以看到，這個Hessian的計算量是很大的，而高斯牛頓法的精髓就是將這個公式的進行了簡化，將下面加號后的式子看做是0來簡化計算，當然要滿足這一假設需要滿足兩個條件：

1. P o l d \boldsymbol{P^{old}} Pold接近于最優解 P ? \boldsymbol{P^{*}} P?【初值估計 P 0 \boldsymbol{P^{0}} P0很重要】
2. 灰度值的二階導接近于0（這一要求主要取決于影像本身和插值平滑性，但一般來說影像灰度過渡都很均勻）【插值演算法的平滑程度很重要】

在這兩項條件滿足的情況下，后面的式子相當于是兩個小量相乘，等于一個更小的量可以近似為0，大大簡化了計算目標函式Hessian的時間，

這種假設是有代價的，通常來說高斯牛頓法會比傳統牛頓拉夫遜法有更多的迭代次數，但從單次迭代上來看，運算代價減少了很多，因此實際使用時需要進行權衡，通常來說子區越小，則用高斯牛頓法的效果會更好，當然，除了計算代價的優化外，高斯牛頓法還擁有更好的收斂性且更容易實作¹

這樣一來我們便可以得到在FA-GN中,所使用的目標函式在 P o l d \boldsymbol{P^{old}} Pold的Hessian為：
? ? C ZNSSD ( P o l d ) ≈ 2 ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ∑ [ d G ( x o + W ( Δ x i x 0 , P o l d ) ) d Δ P ] [ d G ( x o + W ( Δ x i x 0 , P o l d ) ) 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 = ? ? C ZNSSD ( P o l d ) ? ? ? C ZNSSD ( P o l d ) ? 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后按照下面的迭代公式更新即可
P n e w = P o l d + Δ 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的核心思想，那我建議你直接跳到下面的數學推導部分，如果你還是很蒙，那試試這一節能不能讓你弄通順一些，在本文的第一章中，我們在已知形函式引數初值的情況下，可以得到一對參考子區和形變子區的運算式，即

1. f i = F ( x o + Δ x i x 0 ) = F ( x o + W ( Δ x i x 0 , 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))
2. g i = G ( x o + W ( Δ x i x 0 , P 0 ) ) 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去更新形變子區，即在形變影像上找到一個與原參考子區更加匹配的形變子區，那我們是否能反過來，用同樣的方法在參考影像上找到一個與形變子區更加匹配的參考子區呢？其答案是可以的，從下圖¹可以看到這三者的一個關系

圖中綠色區域即為原始參考子區， P r c P_{rc} Prc?代表由原參考子區變換到形變子區的形函式引數，即上面公式中的 P 0 \boldsymbol{P^{0}} P0； P r r P_{rr} Prr?代表由原參考子區變換到新參考子區的形函式引數，而這個塊新參考子區與形變子區的匹配程度無疑是更高的，這個新的參考子區可以表示為：
f i ^ = F ( x o + W ( Δ x i x 0 , P r r ) ) \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?))當 P r r = 0 P_{rr}=0 Prr?=0時，即代表著新參考子區與原參考子區完全重合，那這樣一來，我們就有了一個新的目標，與FA-GN或者傳統牛頓拉夫遜不同，他們是在單純的找滿足最大相關標準時候的 P r c ? P_{rc}^{*} Prc??，而直接求解這個值是很麻煩的，而對于 P r r P_{rr} Prr?其最優解就是 P r r ? = 0 P_{rr}^{*}=0 Prr??=0，只要想辦法讓這個形函式引數逼近0即可，

這樣一來我們就可以更新下我們對于數字影像相關的歸納：通過在參考影像上設定好子區 f i = F ( x o + Δ x i ) f_{i}=F\left(x^{o}+\Delta x_{i}\right) fi?=F(xo+Δxi?) ，基于給定的一組形函式引數初值 P r c 0 P_{rc}^{0} Prc0?及 P r r 0 P_{rr}^{0} Prr0?，對目標函式（相關標準）不斷迭代優化從而得到一個在 P r r = 0 P_{rr}=0 Prr?=0的情況下的最優解 P r c ? P_{rc}^{*} Prc??，從而得到一組在形變影像上與參考子區最佳匹配的形變子區 g i ? = G ( x o + W ( Δ x i , P r c ? ) ) g_{i}^{*}=G\left(x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}_{i}, \boldsymbol{P_{rc}^{*}}\right)\right) gi??=G(xo+W(Δxi?,Prc??))，最終實作參考影像像素點與形變影像像素點的匹配， 而我們針對優化 P r c P_{rc} Prc?去使用的高斯牛頓法就被叫做FA-GN，而針對優化 P r r P_{rr} Prr?去使用的高斯牛頓法就被叫做IC-GN，

2.數學推導

同樣的，為了展示IC-GN與其他演算法的不同性以及其優點，這里將其具體運算展開來寫，設上一次迭代的到的形函式引數為 P o l d \boldsymbol{P^{old}} Pold，由于這一次是針對 P r r P_{rr} Prr?來做，故求解的增量 Δ P \Delta \boldsymbol P ΔP即為 P r r P_{rr} Prr?，其每次的起點都是原參考子區，即形函式引數零點，如此我們的目標函式為 C ZNSSD ( Δ P ) C_{\text {ZNSSD }}(\boldsymbol{\Delta \boldsymbol P}) CZNSSD ?(ΔP)，即
C ZNSSD ( Δ P ) = ∑ i = 1 I [ F ( x o + W ( Δ x i x 0 , Δ P ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , Δ P ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] 2 C_{\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ˉ?)2 ?F(xo+W(Δxix0,ΔP))?fˉ???∑i=1I?(G(xo+W(Δxix0,Pold))?gˉ?)2 ?G(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 0 dΔPd(fˉ?)?≈0 d ( ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 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處的梯度為：
? C ZNSSD ( 0 ) = d C ZNSSD ( 0 ) d Δ P ≈ 2 ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ∑ [ [ F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] [ d F ( x o + W ( Δ x i x 0 , 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ˉ?)2 ?2?∑??????∑i=1I?(F(xo+W(Δxix0,0))?fˉ?)2 ?F(xo+W(Δxix0,0))?fˉ???∑i=1I?(G(xo+W(Δxix0,Pold))?gˉ?)2 ?G(xo+W(Δxix0,Pold))?gˉ?????[dΔPdF(xo+W(Δxix0,0))?]????
再求解目標函式在0處的Hessian為：
? ? C ZNSSD ( 0 ) = d 2 C ZNSSD ( 0 ) d ( Δ P ) 2 ≈ 2 ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 { ∑ [ d d Δ P F ( x o + W ( Δ x i x 0 , 0 ) ) ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ] [ d F ( x o + W ( Δ x i x 0 , 0 ) ) d Δ P ] T + ∑ [ F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ∑ i = 1 I ( F ( x o + W ( Δ x i x 0 , 0 ) ) ? f ˉ ) 2 ? G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ∑ i = 1 I ( G ( x o + W ( Δ x i x 0 , P o l d ) ) ? g ˉ ) 2 ] [ d 2 F ( x o + W ( Δ x i x 0 , 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ˉ?)2 ?2?????∑???∑i=1I?(F(xo+W(Δxix0,0))?fˉ?)2 ?dΔPd?F(xo+W(Δxix0,0))????[dΔPdF(xo+W(Δxix0,0))?]T+∑???∑

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/380864.html

標籤：AI

上一篇：實作Strcpy函式 - 通過函式發現 “程式之美” | 不斷優化、優化、再優化~

下一篇：機器學習（一）：回歸