連續(xù)全場數(shù)字圖像相關(guān)的反向組合高斯－牛頓算法*

盧 宇，陳 果，申 遠，陳曉婷

(合肥師范學院電子信息與電氣工程學院，安徽 合肥230601)

數(shù)字圖像相關(guān)(Digital Image Correlation，DIC)[1－3]技術(shù)是一種廣泛應(yīng)用于位移和變形測量的非接觸圖像測量方法，是材料科學領(lǐng)域一種流行而有效的方法[4－6]?；?D或3D-DIC技術(shù)的視頻引伸計使用CCD攝像機拍攝試樣表面散斑的全場位移與變形[7－9]。

DIC算法是由Peter[10]和Yamaguchi[11]分別在20世紀80年代提出的。對DIC算法的研究主要集中在算法的準確性和效率[12－13]。DIC算法按精度可分為整數(shù)像素法和亞像素法。整數(shù)像素算法采用零均值歸一化互相關(guān)(zero-mean nomalized cross-correlation，ZNCC)[14－15]準則。相關(guān)系數(shù)的最大或最小峰值位置即為最佳配準點。為了提高ZNCC算法的計算效率，Lewis[15]在計算過程中引入了求和表和快速傅立葉變換，提出了快速零均值歸一化互相關(guān)(fast zero-mean nomalized cross-correlation，F(xiàn)ZNCC)。

為了提高配準精度，研究人員在不增加相機和鏡頭成本的前提下，提出了各種亞像素配準方法。這些方法包括相關(guān)系數(shù)擬合法和圖像灰度梯度下降法。其中灰度梯度下降法比相關(guān)系數(shù)擬合法更適合于圖像變形的測量[16]，也是DIC的標準方法。前向加性牛頓－拉弗森(forward additive Newton-Raphson，F(xiàn)ANR)[17]算法是第一個提出來的梯度下降算法，目前已成為DIC技術(shù)中應(yīng)用最廣泛的算法。然而，F(xiàn)ANR算法必須計算兩個圖像子集的相似度，并在每個迭代步驟中計算海森矩陣。這樣的過程往往很耗時。Baker等人[18]提出了反向組合高斯－牛頓(inverse composition Gauss-Newton，IC-GN)算法，該算法比FA-NR算法效率更高，并獲得了幾乎一致的精度。IC-GN算法是當前DIC研究的熱點之一[19]。然而，IC-GN算法在提出時只考慮了單點圖像配準。為了進行全場圖像變形測量，Pan等人[20]提出了雙立方插值系數(shù)緩存表和初始值轉(zhuǎn)移策略來加速全場配準計算的效率。Lu等人[21]提出了亞像素擬合初始值策略和子區(qū)域圖像邊緣擴展插值方法。這些方法提高了兩幅圖像全場配準的精度和效率。但是，在材料的力學性能測試過程中，CCD攝像機通常會獲取數(shù)百、上千幅試樣的圖像。如果采用上述全場圖像配準方式，則圖像序列中每幅變形圖像都必須與參考(未變形)圖像進行配準，并且每次配準過程都需要計算參考圖像中每個子區(qū)域的海森矩陣。這個過程消耗了大量的計算時間。因此，本研究提出了一種連續(xù)的全場圖像配準方法，將海森矩陣緩存在緩存表中，以提高實際工程應(yīng)用的效率。

本文首先介紹了傳統(tǒng)的基于ZNCC和IC-GN算法的全場變形測量算法。隨后，在基于雙立方插值系數(shù)緩存表的 IC-GN(IC-GN with Bicubic interpolation coefficient cache table，IC-GN_BICCT)算法[20]的基礎(chǔ)上又將海森矩陣緩存到IC-GN算法中，形成了基于雙立方插值系數(shù)緩存表和海森矩陣緩存表的IC-GN(IC-GN with Bicubic interpolation coefficient cache table and Hessian matrix cache table，ICGN_BIC_HMCT)算法，以縮短后續(xù)圖像的配準時間。通過對計算機生成的人工模擬散斑圖像和CCD攝像機采集的真實散斑圖像進行實驗，驗證了本文提出的基于算法有效性和準確性。結(jié)果表明，該算法可以在不損失精度的情況下提高IC-GN_BICCT算法的效率。

1 數(shù)字圖像相關(guān)方法測量散斑圖像變形的原理

2D-DIC技術(shù)被應(yīng)用于材料力學性能測試實驗中。在被測試樣表面噴涂上黑白散斑并使用萬能材料試驗機拉伸，與此同時CCD攝像機采集散斑的平移和變形圖像，可以通過這些圖像來計算試樣的位移和變形。試樣拉伸前采集的圖像稱為參考圖像，拉伸后的圖像稱為變形圖像。圖1給出了用視頻引伸計測量材料力學性能的實驗裝置示意圖。該系統(tǒng)由高分辨率定焦鏡頭、CCD攝像機、萬能材料試驗機、單色光源和高性能計算機組成。測試前被測試樣固定在萬能材料試驗機的夾具上。CCD攝像機的光軸應(yīng)盡可能垂直于試樣表面。

圖1 材料力學性能測量實驗裝置示意圖

1.1 從粗到細的數(shù)字圖像相關(guān)方法

圖2 為CCD攝像機采集到的真實散斑圖像。通常在圖2(a)參考圖像中選取一塊尺寸為(2n＋1)×(2n＋1)像素尺寸的參考子區(qū)域，該子區(qū)域中心坐標為Pref(xref_c，yref_c)，與之對應(yīng)的來自于變形圖像的變形子區(qū)域的中心坐標為Pcur(xcur_c，ycur_c)，如圖2(b)所示。為了找到配準點Pcur，可根據(jù)經(jīng)驗在變形圖像中劃出一個尺寸為(M×N)像素的搜索區(qū)域，然后使用由粗到細的方法搜索出Pcur。

圖2 圖像配準策略

ZNCC算法用于確定整像素精度級別的最佳配準點。ZNCC計算公式如式(1)所示[15]:

式中:g(x m，y n)表示以(x m，y n)為中心的變形子區(qū)域的圖像灰度分布，(x m，y n)∈(M×N)。ˉf和ˉg分別為參考子區(qū)域和變形子區(qū)域的平均灰度值。S表示遍歷子區(qū)域中所有的像素。ZNCC算法的優(yōu)點是消除了背景光強度變化對配準的影響，而FZNCC算法通常用于提高ZNCC算法的計算效率。通過找到最大值CZNCC來確定變形子集在搜索區(qū)域中的位置(xcur_int，ycur_int)。

然而，ZNCC法只能獲得位移而不能獲得變形。在亞像素圖像配準中，基于零均值歸一化誤差平方和準則(zero-mean normalized sum of squared difference，ZNSSD)的IC-GN算法是一種常用的亞像素配準方法。在說明IC-GN算法前，首先需要引入形函數(shù)W來描述參考子區(qū)域與變形子區(qū)域的對應(yīng)像素位置關(guān)系:

式中:p＝{u，v，u x，u y，vx，v y}T是2D仿射變換向量，其由6個獨立的參數(shù)組成并定義了三種線性變換(即平移、剪切和縮放)。x＝(x，y)表示參考子區(qū)域內(nèi)的像素坐標，(x′，y′)表示變形子區(qū)域內(nèi)的像素坐標。另外，dx＝x－xref_c，dy＝y(tǒng)－yref_c。dx和dy的取值范圍是－n≤dx，dy≤n。CZNCC的最大值坐標為(xcur_int，ycur_int)作為p中(u，v)的初始值。

基于ZNSSD的IC-GN算法可以表示為[12]:

式中:Δp是p的增量。式(3)中

式(3)的解為:

因此，Δp＝{Δu，Δv，Δu x，Δu y，Δv x，Δv y}T可通過式(4)得到。式(4)中的海森矩陣可以描述為:

式中:?f＝((?f/?x)，(?f/?y))是參考子區(qū)域的灰度梯度，?W/?p表示式(2)中對p求偏導:

上述算法過程即為由粗到細，由整像素到亞像素的數(shù)字圖像相關(guān)算法，可以簡單表示為FZNCC＋IC-GN。

1.2 基于雙立方插值系數(shù)緩存表和海森矩陣緩存表的IC-GN算法

在本節(jié)中，我們將IC-GN_BICCT算法引入到連續(xù)幀圖像配準過程中。在式(4)中，參考子區(qū)域的?f·(?W/?p)、H、H－1都需要計算，而且計算過程非常耗時。對于連續(xù)的圖像配準，采取的配準策略往往如圖3所示。第1步使用變形圖像1與參考圖像配準，第2步使用變形圖像2與參考圖像配準，以此類推。如果式(4)中參考圖像的海森矩陣緩存表可以在后續(xù)圖像配準過程中重復使用，則計算時間將大大減少。

圖3 連續(xù)數(shù)字圖像相關(guān)計算的示意圖

全場配準計算時，需要在參考圖像中選取被測區(qū)域，即感興趣區(qū)域(RoI)。為了減少計算量，一般不會計算RoI中的所有像素點的位移與變形情況，而是在圖像上設(shè)置具有固定間隔的網(wǎng)格，只計算網(wǎng)格上的種子點的位移與變形。

根據(jù)IC-GN_BICCT的步驟，使用IC-GN_BIC_HMCT算法計算連續(xù)全場數(shù)字圖像相關(guān)的步驟如下:

(1)整個計算過程的預(yù)計算:

①在參考圖像中選擇出要計算的種子點矩陣的像素坐標；

②建立參考圖像的二維海森矩陣緩存表M H，其行列數(shù)與種子點矩陣行列數(shù)相同；

③計算參考圖像的?f，ˉf，f s；

(2)每幅變形圖像的預(yù)計算:

建立變形圖像的二維插值系數(shù)緩存表M c，其行列數(shù)與種子點矩陣行列數(shù)相同；

(3)每個參考子區(qū)域的預(yù)計算:

①根據(jù)式(5)計算種子點的H及H－1；

②如果是首次計算H及H－1，將H－1存入緩存表M H；否則在M H中查詢；

(4)每幅圖像的迭代:

①如果是首次計算插值系數(shù)，將系數(shù)存入緩存表M c；否則在M c中查詢；

②使用形函數(shù)W(x；p)和M c計算g(W(x；

⑤使用式(4)計算Δp；

⑥更新W(x；p)←W(x；p)W(x；Δp)－1直至‖且ε＝10－3；

在IC-GN_BIC_HMCT算法步驟中(1)②、(3)①、(3)②是較IC-GN_BICCT算法增加的步驟。

2 實驗和分析

本節(jié)利用計算機生成的人工模擬散斑圖像和CCD攝像機采集的真實散斑圖像，對比了IC-GN_BIC_HMCT算法和IC-GN_BICCT算法的計算效率。這兩種算法的初始值均來自式(1)的ZNCC算法。

實驗的硬件條件為4.0 GHz的Intel Core i3－8350K CPU，16 GB內(nèi)存的臺式機。IC-GN_BIC_HMCT、IC-GN_BICCT和IC-GN、FZNCC算法都用C＋＋和OpenCV庫編寫的。

2.1 使用人工模擬散斑圖像的驗證

我們首先使用一組來自實驗力學學會(SEM)的2D-DIC challenge數(shù)據(jù)集[22]的散斑圖像來驗證本文所提出的算法的計算效率。其中來自該數(shù)據(jù)集的Sample 1是由一組具有恒定噪聲數(shù)值的人工模擬散斑圖像組成。圖像的尺寸為512×512像素，灰度階數(shù)為8位，灰度值的標準偏差為1.5。Sample 1使用TexGen[23]方法生成。Sample 1中的圖像在X和Y方向上同時移動，范圍從0像素到1像素，兩幅連續(xù)圖像之間的移動增量為0.05像素。Sample 1中有21幅圖像，并且這些圖像的對比度逐級下降。測試時，Sample 1的第1幅圖像被視為參考圖像。

在參考圖像中選取從點(100，100)到點(400，400)的RoI用于計算，如圖4的矩形區(qū)域所示。為了減少計算量，在X和Y方向上每隔10個像素選取1個種子點。因此RoI中共有961個種子點。

圖4 Sample 1中的參考圖像和RoI

為了定量研究IC-GN_BIC_HMCT和IC-GN_BICCT算法的誤差，我們使用了如下的平均誤差精度標準:

式中:N表示需要計算的種子點數(shù)量，u i，j為計算得到的參考圖像中坐標為(i，j)的種子點的位移量，upreset是預(yù)設(shè)值.

使用包括參考圖像在內(nèi)的21幅連續(xù)圖像與參考圖像來計算全場位移并記錄誤差，結(jié)果如圖5所示。圖中分別使用尺寸為21×21、31×31、41×41像素大小的子區(qū)域。由于這兩種算法的原理是相同的，所以兩種算法的平均誤差曲線是重疊的。圖中可見兩種算法的最大平均誤差接近0.01像素，但對比度的變化不影響測量精度。圖5中只給出了X方向位移的誤差曲線。

圖5 Sample 1的IC-GN_BIC_HMCT和IC-GN_BICCT算法全場位移平均誤差曲線

下一步是比較計算時間。圖6給出了使用不同尺寸子區(qū)域的IC-GN_BIC_HMCT和IC-GN_BICCT的計算時間。在圖6(a)、(b)、(c)中除第1個數(shù)據(jù)點外，IC-GN_BIC_HMCT的計算時間明顯短于ICGN_BICCT。使用21×21尺寸子區(qū)域的IC-GN_BICCT的總時間為7.50 s，平均每幅圖像為0.36 s。而IC-GN_BIC_HMCT的總時間為6.66 s，平均每幅圖像為0.32 s。后者的效率比前者高11.35%。圖6(a)、(b)、(c)中的第1個數(shù)據(jù)點，IC-GN_BIC_HMCT消耗的時間略長于IC-GN_BICCT的算法。這是因為IC-GN_BIC_HMCT算法有數(shù)據(jù)緩存步驟，消耗更多存儲時間。

圖6 使用3種不同尺寸子區(qū)域的IC-GN_BIC_HMCT和IC-GN_BICCT的計算時間比較

對于31×31和41×41的子區(qū)域，IC-GN_BIC_HMCT的效率分別提高了14.95%和18.93%。可以得出結(jié)論，子區(qū)域的尺寸越大，本算法的效率提高越多。

為了更進一步解釋為什么IC-GN_BIC_HMCT算法比IC-GN_BICCT耗時少，我們在表1中列出了IC-GN_BICCT每一步的執(zhí)行時間。表1中記錄了圖4所示的RoI左上方起的種子點1的計算時間。在表1中，第1行的海森矩陣及逆的計算步驟占算法總耗時的20.60%。如果這個計算時間減少，那么計算效率就會提高。

表1 用3種不同尺寸子區(qū)域的IC-GN_BICCT主要步驟的消耗時間 單位:ms

在表1中，當子區(qū)域大小為31×31和41×41像素時，計算海森矩陣及逆的時間成本分別為IC-GN_BICCT總時間成本的22.80%和29.79%。這也說明了為什么子區(qū)域尺寸越大，效率提高越多。

2.2 使用CCD攝像機采集的真實圖像進行位移和變形測量

為了實際驗證本算法對真實散斑圖像的有效性和效率，我們使用萬能材料試驗機進行了試樣拉伸實驗，計算試樣的位移和變形量。通過噴涂油漆，被測試樣上覆蓋了黑白散斑圖案。CCD攝像機被用來采集拉伸前后的一組圖像，圖7為試樣未變形時的參考圖像。試樣受力方向為圖中Y軸正方向。圖像的分辨率為1 624×1 234像素，灰度階數(shù)為8位。攝像機采樣頻率為2幀/s。在試樣斷裂之前，總拉伸時長約為45 s。RoI尺寸設(shè)置為121像素×365像素。在X、Y坐標方向上，每兩個種子點的間隔為10個像素。因此需要計算的種子點數(shù)量為13×37＝481個。為了提高計算效率，計算變形圖時使用21×21像素尺寸子區(qū)域。

圖7 被測試樣與RoI

獲得RoI中所有種子點的位移后，還要使用雙立方插值方法重建全場變形圖，每兩個配準的種子點之間需要插入額外的9個點。圖8顯示了第10幀變形圖像的X軸和Y軸上RoI的位移和變形分布圖。圖中的標尺上的數(shù)字表示位移和變形量的大小。圖中可見試樣在Y軸方向上伸長，而在X軸方向上收縮。

圖8 第10幅變形圖像的X、Y方向上試樣的全場變形圖

下面分析本算法對連續(xù)變形場計算效率的提高。表2給出了分別使用IC-GN_BIC_HMCT和IC-GN_BICCT[20]、IC-GN算法[18－19]計算連續(xù)20幀圖像配準的平均時間。可以看出由于使用了插值系數(shù)緩存方法，使用21×21像素尺寸子區(qū)域時，計算時間從IC-GN算法的0.42 s/幀減少到0.21 s/幀，而使用41×41像素尺寸子區(qū)域時，計算時間更是從IC-GN算法的1.23 s/幀減少到0.55 s/幀，效率提高了123.63%。

表2 使用3種不同尺寸子區(qū)域的IC-GN算法、IC-GN_BIC_HMCT算法、IC-GN_BICCT消耗時間單位:s/frame

為了進一步比較IC-GN_BIC_HMCT和IC-GN_BICCT兩種算法的計算時間，在圖9給出了每幀圖像的計算效率對比曲線。圖9中使用的是21×21尺寸子區(qū)域。當?shù)?幀變形圖像與參考圖像配準時，由于參考圖像的H矩陣沒有被緩存，IC-GN_BIC_HMCT和IC-GN_BICCT[20]兩種算法計算時間基本相等。從第2幀變形圖像開始，IC-GN_BIC_HMCT的計算時間明顯低于IC-GN_BICCT。IC-GN_BIC_HMCT算法計算20幀圖像配準所用時間為4.27 s，而IC-GN_BICCT算法為5.22 s。前者比后者快18.20%。該結(jié)果顯示了新算法在效率方面的優(yōu)越性。除去第1幅變形圖像，后續(xù)的19幅圖像的計算效率提高了19.02%。

圖9 使用真實散斑圖像的計算效率比對

3 結(jié)論

本文采用了一種新的連續(xù)散斑圖像配準方法，稱為IC-GN_BIC_HMCT算法?？紤]到IC-GN中需要計算海森矩陣及逆，本研究致力于減少連續(xù)數(shù)字圖像相關(guān)中的IC-GN計算時間。在已有變形圖像的雙立方插值系數(shù)緩存表的基礎(chǔ)上，又建立了參考圖像海森矩陣緩存表。利用緩存表避免了海森矩陣及逆的重復計算。通過實驗驗證了該算法的正確性和有效性。在人工模擬散斑圖像實驗中，對于3種不同大小的子區(qū)域，本算法的效率分別比沒有海森矩陣緩存表的IC-GN-BICCT算法高11.35%、14.95%和18.93%，但精度相同。在真實散斑圖像實驗中，計算效率提高了18.20%。測試結(jié)果表明了該算法在效率上的優(yōu)越性。