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

盧 宇，陳 果，申 遠(yuǎn)，陳曉婷

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

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

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

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

本文首先介紹了傳統(tǒng)的基于ZNCC和IC-GN算法的全場(chǎng)變形測(cè)量算法。隨后，在基于雙立方插值系數(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ù)圖像的配準(zhǔn)時(shí)間。通過(guò)對(duì)計(jì)算機(jī)生成的人工模擬散斑圖像和CCD攝像機(jī)采集的真實(shí)散斑圖像進(jìn)行實(shí)驗(yàn)，驗(yàn)證了本文提出的基于算法有效性和準(zhǔn)確性。結(jié)果表明，該算法可以在不損失精度的情況下提高IC-GN_BICCT算法的效率。

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

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

圖1 材料力學(xué)性能測(cè)量實(shí)驗(yàn)裝置示意圖

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

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

圖2 圖像配準(zhǔn)策略

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

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

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

式中:p＝{u，v，u x，u y，vx，v y}T是2D仿射變換向量，其由6個(gè)獨(dú)立的參數(shù)組成并定義了三種線性變換(即平移、剪切和縮放)。x＝(x，y)表示參考子區(qū)域內(nèi)的像素坐標(biāo)，(x′，y′)表示變形子區(qū)域內(nèi)的像素坐標(biāo)。另外，dx＝x－xref_c，dy＝y(tǒng)－yref_c。dx和dy的取值范圍是－n≤dx，dy≤n。CZNCC的最大值坐標(biāo)為(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可通過(guò)式(4)得到。式(4)中的海森矩陣可以描述為:

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

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

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

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

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

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

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

(1)整個(gè)計(jì)算過(guò)程的預(yù)計(jì)算:

①在參考圖像中選擇出要計(jì)算的種子點(diǎn)矩陣的像素坐標(biāo)；

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

③計(jì)算參考圖像的?f，ˉf，f s；

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

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

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

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

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

(4)每幅圖像的迭代:

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

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

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

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

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

2 實(shí)驗(yàn)和分析

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

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

2.1 使用人工模擬散斑圖像的驗(yàn)證

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

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

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

為了定量研究IC-GN_BIC_HMCT和IC-GN_BICCT算法的誤差，我們使用了如下的平均誤差精度標(biāo)準(zhǔn):

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

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

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

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

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

對(duì)于31×31和41×41的子區(qū)域，IC-GN_BIC_HMCT的效率分別提高了14.95%和18.93%?？梢缘贸鼋Y(jié)論，子區(qū)域的尺寸越大，本算法的效率提高越多。

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

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

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

2.2 使用CCD攝像機(jī)采集的真實(shí)圖像進(jìn)行位移和變形測(cè)量

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

圖7 被測(cè)試樣與RoI

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

圖8 第10幅變形圖像的X、Y方向上試樣的全場(chǎng)變形圖

下面分析本算法對(duì)連續(xù)變形場(chǎng)計(jì)算效率的提高。表2給出了分別使用IC-GN_BIC_HMCT和IC-GN_BICCT[20]、IC-GN算法[18－19]計(jì)算連續(xù)20幀圖像配準(zhǔn)的平均時(shí)間?？梢钥闯鲇捎谑褂昧瞬逯迪禂?shù)緩存方法，使用21×21像素尺寸子區(qū)域時(shí)，計(jì)算時(shí)間從IC-GN算法的0.42 s/幀減少到0.21 s/幀，而使用41×41像素尺寸子區(qū)域時(shí)，計(jì)算時(shí)間更是從IC-GN算法的1.23 s/幀減少到0.55 s/幀，效率提高了123.63%。

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

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

圖9 使用真實(shí)散斑圖像的計(jì)算效率比對(duì)

3 結(jié)論

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