數(shù)字圖像相關(guān)的多尺度圖像子區(qū)匹配算法

林 哲，蔡 恬，王燕鋒

（1.汕頭職業(yè)技術(shù)學(xué)院計(jì)算機(jī)系，廣東 汕頭 515078；2.汕頭職業(yè)技術(shù)學(xué)院網(wǎng)絡(luò)與信息中心，廣東 汕頭 515078；3.中原工學(xué)院建筑工程學(xué)院，河南 鄭州 450007；4.汕頭大學(xué)工學(xué)院，廣東 汕頭 515063）

0 引 言

自1980 年以來，數(shù)字圖像相關(guān)（DIC）被廣泛用于測量固體表面的位移和應(yīng)變[1-4]，它是一種重要的非接觸式光學(xué)實(shí)驗(yàn)方法.通常的做法是用一臺固定的數(shù)碼相機(jī)在試件變形前后的瞬間對其表面拍攝兩次照片，這兩張圖像分別稱為參考圖像和變形圖像，然后用計(jì)算機(jī)算法對它們進(jìn)行處理，確定表面上每個點(diǎn)的相對位移.為了計(jì)算參考圖像與變形圖像之間點(diǎn)與點(diǎn)的關(guān)系，基于圖像子區(qū)（subset）的方法在實(shí)踐中得到了廣泛的應(yīng)用.對于變形圖像中的每一個點(diǎn)，利用其周圍區(qū)域的局部信息（即圖像子區(qū)）對參考圖像中的最佳對應(yīng)點(diǎn)進(jìn)行跟蹤.

在過去30 年里，為了提高基于圖像子區(qū)的DIC 的精度和效率，人們做了大量的工作.例如，Schreier 等利用高階樣條插值函數(shù)實(shí)現(xiàn)了0.01 像素的測量精度[5]，并且分析了因形狀函數(shù)匹配不足而產(chǎn)生的系統(tǒng)誤差，研究了非均勻應(yīng)變場匹配過程中的二階形狀函數(shù)[6].Cheng 等人提出了一個高效、準(zhǔn)確、健壯的框架.采用均勻參數(shù)b 樣條曲面函數(shù).通過迭代優(yōu)化，將整個試件表面積表示為未知的二維變形場[7].此外，關(guān)于搜索過程[8-9]、相關(guān)方法[10-11]、配準(zhǔn)方法[12-18]、測量精度[19-21]、可靠性[22-24]、計(jì)算速度[22-29]等方面研究也有很多報道.

目前，準(zhǔn)確、有效、魯棒地計(jì)算全場測量結(jié)果仍然是一個難題，尤其是在大變形或表面不連續(xù)的情況下.對于基于圖像子區(qū)的方法，如何選擇圖像子區(qū)仍然是一個關(guān)鍵問題.一方面，為了可靠地為變形圖像的每個圖像子區(qū)在參考圖像中找到最佳匹配，圖像子區(qū)的大小必須相對較大，使每個圖像子區(qū)包含足夠的強(qiáng)度信息，便于相互區(qū)分.而另一方面，由于圖像子區(qū)中每個點(diǎn)的位移近似表示為一階或二階形狀函數(shù)，圖像子區(qū)的大小必須相對較小，以盡可能減少近似誤差.一個合適的圖像子區(qū)大小必須足夠大又足夠小，這是一個兩難的問題.如何處理這個問題？

Pan 等人提出了一種選擇圖像子區(qū)大小的方法[30]，利用圖像噪聲方差和圖像子區(qū)強(qiáng)度梯度平方和（Sum of Square of Subset Intensity Gradients，SSSIG），通過理論推導(dǎo)預(yù)測測量精度，提出了一種在計(jì)算位移前選擇合適圖像子區(qū)大小的算法，從11×11 最小圖像子區(qū)開始逐步擴(kuò)大，直到SSSIG 剛好大于閾值.Ghulam Mubashar Hassan 等人提出了一種動態(tài)圖像子區(qū)選擇（Dynamic Subset Selection，DSS）算法.它使用另一個參數(shù)強(qiáng)度變化率在迭代中動態(tài)檢測最佳的圖像子區(qū)大小[31].Liang Z 等人提出了一種新的方法來解決選擇小圖像子區(qū)或大圖像子區(qū)的困境[32]，首先選擇一個相對較大的圖像子區(qū)大小，圖像子區(qū)中的每個像素都有一個系數(shù)，該系數(shù)根據(jù)像素到圖像子區(qū)中心的距離計(jì)算，然后對相關(guān)準(zhǔn)則進(jìn)行了修正，對圖像子區(qū)中的像素做不同的處理，實(shí)際上達(dá)到模糊圖像子區(qū)邊界的效果，從而巧妙地避免了確定最佳圖像子區(qū)大小的問題.

上述選擇圖像子區(qū)大小的方法有一個共同點(diǎn)，即圖像子區(qū)大小是固定的、統(tǒng)一的.這種方法在試件連續(xù)變形時是有效的，但是當(dāng)試件存在損傷區(qū)域（如裂紋時），由于應(yīng)力集中，在后續(xù)變形過程中，損傷區(qū)域附近各點(diǎn)的位移可能會發(fā)生更為復(fù)雜的變化.因此，為了更好地近似表示位移，應(yīng)該更多地關(guān)注圖像子區(qū)中心附近的點(diǎn)，同時，遠(yuǎn)離圖像子區(qū)中心的點(diǎn)也不能完全忽略，因?yàn)樗鼈円部赡馨杏玫钠ヅ湫畔?

為此，本文提出一種多尺度圖像子區(qū)匹配算法，利用權(quán)重函數(shù)調(diào)整圖像子區(qū)內(nèi)各個像素點(diǎn)對圖像子區(qū)相關(guān)性的貢獻(xiàn)，使圖像子區(qū)的邊界完全變得模糊，實(shí)現(xiàn)動態(tài)調(diào)整圖像子區(qū)大小，更加合適大形變或不連續(xù)形變的表面，通過實(shí)驗(yàn)證明了該方法能夠得到足夠精確的測量結(jié)果.

1 數(shù)字圖像相關(guān)方法

為了獲得全場位移，首先必須能夠識別出試件表面的每一點(diǎn)，通常的做法在試件表面噴涂隨機(jī)散斑圖案，使每一個點(diǎn)周圍區(qū)域互不相同.隨著試件的變形，每個區(qū)域也發(fā)生了變形.然后，選取變形圖像中某一點(diǎn)，根據(jù)其周圍圖像子區(qū)的灰度信息在參考圖像中跟蹤對應(yīng)的點(diǎn)，如圖1 所示.

圖1 圖像子區(qū)變形示意圖

在一個圖像子區(qū)中，假設(shè)變形是均勻分布的，即點(diǎn)的位移在水平和垂直方向上都是線性變化的，由形狀函數(shù)進(jìn)行描述：

對照組-男、女占比各為29:21;年齡范圍上限值86歲,下限值62歲,年齡平均值(74.21±1.35)歲。

其中xPQ和yPQ分別是P 點(diǎn)距離其圖像子區(qū)中心點(diǎn)Q 的水平距離和垂直距離.

由一階泰勒展開得到

其中，xrefQ和yrefQ表示參考圖像中Q 點(diǎn)的坐標(biāo).xrefP和yrefP表示參考圖像中P 點(diǎn)的坐標(biāo).xcurP和ycurP表示變形圖像中P 的坐標(biāo).u 和v 分別為Q 在x 方向和y 方向的位移.它們的偏導(dǎo)數(shù)表示為DIC 分析的目的就是找出

為了評價變形圖像子區(qū)與參考圖像子區(qū)之間的相關(guān)性，需要一個相關(guān)準(zhǔn)則.到目前為止，已經(jīng)報道了許多相關(guān)標(biāo)準(zhǔn)，如交叉相關(guān)（cross-correlation，CC）、歸一化交叉相關(guān)（normalized cross-correlation，NCC）、平方和差（sum of squared difference，SSD）等.然后，對每個點(diǎn)應(yīng)用FA-NR 算法[19]或IC-NG 算法[23]求p→的最優(yōu)值，使相關(guān)系數(shù)取得最小值.

2 多尺度圖像子區(qū)匹配算法

從觀測尺度來說，每一個點(diǎn)都是根據(jù)其周圍的灰度信息來識別的，這些灰度信息通過一個窗口即圖像子區(qū)來觀察.然而，通常很難確定觀測窗的大小，因?yàn)橛^測窗必須足夠大，才能包含唯一的灰度信息，但又必須足夠小，以用形狀函數(shù)近似表示內(nèi)部位移分布.這就是在現(xiàn)有方法中難以自動確定最優(yōu)圖像子區(qū)大小的根本原因.為此，本文提出在多尺度下重新考慮圖像子區(qū)大小選擇問題.

正如我們所知，多尺度分析在圖像處理中已經(jīng)被廣泛應(yīng)用多年，已經(jīng)開發(fā)出許多分析工具，例如 Curvlet[34]、Wedgelet[35]、Bandelet[36]、Grouplet[37-38]、Beamlet[39-40].多尺度分析的主要思想是將一個信號（圖像）分解成多個尺度的分量，通過不同的觀測窗口提取出特定的特征.

在DIC 中，為了將變形圖像中的一個點(diǎn)與參考圖像中的另一個點(diǎn)進(jìn)行匹配，僅使用圖像本身的信息是不夠的，因?yàn)槿魏我粡垐D像中都可能有大量相同灰度值的點(diǎn).因此才必須考慮圖像子區(qū)中包含的周圍信息，根據(jù)一定相關(guān)準(zhǔn)則來區(qū)分各個點(diǎn).

權(quán)重函數(shù)在評價兩個圖像子區(qū)相似性時確定各個點(diǎn)的重要性，其中Δx 和Δy 表示當(dāng)前點(diǎn)到圖像子區(qū)中心點(diǎn)的x 方向距離和y 方向距離.

設(shè)f（x，y）表示參考圖像，g（x，y）表示測試圖像，圖像子區(qū)的邊長為2M+1，以圖像子區(qū)中心點(diǎn)為坐標(biāo)原點(diǎn)，則最小平方（Least Square，LS）相關(guān)準(zhǔn)則[33]定義為：

其中，

將其推廣為加權(quán)系數(shù)：

其中

基于多尺度圖像子區(qū)匹配的算法流程圖如圖2 所示.

圖2 基于多尺度圖像子區(qū)匹配的算法流程圖

對每個要處理的點(diǎn)搜索其最佳尺度因子a.首先設(shè)置a 的初始值，計(jì)算第一個點(diǎn)時，a 的初始值設(shè)置為0，計(jì)算其余點(diǎn)時，a 的初始值從前一個已計(jì)算的點(diǎn)繼承而來的.然后，對當(dāng)前點(diǎn)進(jìn)行并行計(jì)算，找出最大加權(quán)系數(shù)Cw（a）并更新接下來，使a 變大或變小，直到使Cw（a）得到最優(yōu)值，此時，最后更新就是所求結(jié)果.

3 實(shí)驗(yàn)結(jié)果與分析

為了準(zhǔn)確評估測量精度，首先采用Matlab 模擬生成數(shù)字散斑圖像.散斑圖像的大小為W×H，散斑密度為ρ，每個散斑具有固定半徑R，中心點(diǎn)坐標(biāo)隨機(jī)均勻分布，由這N 個散斑疊加形成散斑圖像的表達(dá)式為：

其中N=W×H×ρ 表示散斑個數(shù)，（xn，yn）表示第n 個散斑的中心點(diǎn)坐標(biāo).xn和yn分別是均勻分布在[1，W]和[1，H]之間的隨機(jī)整數(shù).但是，由于數(shù)字圖像的像素值的取值范圍是[0，1]，故對I 必須再做規(guī)范化處理：

設(shè)置散斑密度ρ=0.05，散斑半徑R=1，生成的256×256 數(shù)字散斑圖像如圖3.

圖3 模擬生成的數(shù)字散斑圖像

假設(shè)材料泊松比μ=0.2，在計(jì)算機(jī)上模擬試件受y 方向的平行壓力時表面散斑圖像的變化，具體地，將位移參數(shù)設(shè)置為ux=0.02，vy=-0.1，uy=vx=0.

接下來，取圖像中心點(diǎn)（128，128）作為種子點(diǎn)，經(jīng)過IC-NG 算法10 次迭代得到最終的變形向量

表1 模擬DIC 實(shí)驗(yàn)中迭代計(jì)算過程中的相關(guān)性系數(shù) C（）

表1 模擬DIC 實(shí)驗(yàn)中迭代計(jì)算過程中的相關(guān)性系數(shù) C（）

迭代次數(shù) C（Δp→）1 0.078 780 758 2 0.071 938 208 3 0.071 969 003 4 0.071 825 928 5 0.070 442 646 6 0.069 724 213 7 0.069 920 925 8 0.068 876 958 9 0.068 103 011 10 0.068 103 011

圖4 模擬實(shí)驗(yàn)中迭代計(jì)算過程中的相關(guān)性系數(shù) C（）的變化趨勢圖

最后得到的計(jì)算結(jié)果及誤差如表2 所示.

表2 模擬DIC 實(shí)驗(yàn)的計(jì)算結(jié)果與誤差

表2 同時給出了文獻(xiàn)[31]方法和文獻(xiàn)[32]方法的計(jì)算結(jié)果作為對比，可以看出，本文方法的計(jì)算精度優(yōu)于另外兩種方法，其中，本文方法的計(jì)算結(jié)果與文獻(xiàn)[32]方法的計(jì)算結(jié)果更接近，說明像素加權(quán)策略對提高計(jì)算精度起主要作用，同時由于本文方法吸收了文獻(xiàn)[31]方法的動態(tài)調(diào)整策略，兩者結(jié)合進(jìn)一步提高了計(jì)算精度.

本文方法計(jì)算所有點(diǎn)后可以得到x 方向的位移場u 和y 方向的位移場v，如圖5 所示.

圖5 模擬DIC 實(shí)驗(yàn)的計(jì)算結(jié)果

為了進(jìn)一步驗(yàn)證本文方法的性能，在汕頭大學(xué)材料實(shí)驗(yàn)室進(jìn)行一次真實(shí)的DIC 實(shí)驗(yàn)，試件是一塊混凝土塊，利用萬能試驗(yàn)機(jī)在y 軸方向加壓直至破碎.

實(shí)驗(yàn)開始前先拍攝一張參考圖像，如圖6（a）所示，隨著壓力增加，試件開始變形，左上角區(qū)域出現(xiàn)裂紋并逐步擴(kuò)展，如圖6（b）所示，DIC 計(jì)算結(jié)果得到的x 方向的位移場u和應(yīng)變場 Exx 分別如圖 6（c）和圖 6（d）所示.圖 6（c）中藍(lán)綠色表示向負(fù)方向（左）位移，紅黃色表示向正方向（右）位移，可以看出在，在左側(cè)有明顯的黃-藍(lán)分界線，說明試件在此處發(fā)生了不連續(xù)形變，進(jìn)一步從圖6（d）可以看到很明顯的應(yīng)力集中現(xiàn)象（紅色區(qū)域），對比圖6（b）可知此處正是裂紋所在的區(qū)域.本文方法的結(jié)果很好地反映了試件發(fā)生不連續(xù)形變時的位移場和應(yīng)變場.

圖6 真實(shí)DIC 實(shí)驗(yàn)的計(jì)算結(jié)果

4 結(jié)論

針對數(shù)字圖像相關(guān)（DIC）中如何選擇圖像子區(qū)大小的難題，本文提出一種多尺度圖像子區(qū)匹配算法，從多尺度的角度給出一種解決方案.首先定義權(quán)重函數(shù)族，用于調(diào)整圖像子區(qū)內(nèi)各個像素點(diǎn)對圖像子區(qū)相關(guān)性的貢獻(xiàn)，通過迭代過程動態(tài)調(diào)整圖像子區(qū)中各點(diǎn)權(quán)重，使不同位置的圖像子區(qū)具有不同的觀測尺度因子，從而更加合適不均勻形變、大形變或不連續(xù)形變的表面.最后通過實(shí)驗(yàn)驗(yàn)證本文算法的性能，模擬DIC 實(shí)驗(yàn)中測量位移的相對誤差為1.6%和2.8%，在材料實(shí)驗(yàn)室中進(jìn)行的真實(shí)DIC 實(shí)驗(yàn)也驗(yàn)證了本文方法在材料表面發(fā)生不連續(xù)形變時能夠有效監(jiān)測到位移和應(yīng)變的變化.