數字圖像相關的多尺度圖像子區匹配算法

林 哲，蔡 恬，王燕鋒

（1.汕頭職業技術學院計算機系，廣東 汕頭 515078；2.汕頭職業技術學院網絡與信息中心，廣東 汕頭 515078；3.中原工學院建筑工程學院，河南 鄭州 450007；4.汕頭大學工學院，廣東 汕頭 515063）

0 引 言

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

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

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

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

上述選擇圖像子區大小的方法有一個共同點，即圖像子區大小是固定的、統一的.這種方法在試件連續變形時是有效的，但是當試件存在損傷區域（如裂紋時），由于應力集中，在后續變形過程中，損傷區域附近各點的位移可能會發生更為復雜的變化.因此，為了更好地近似表示位移，應該更多地關注圖像子區中心附近的點，同時，遠離圖像子區中心的點也不能完全忽略，因為它們也可能包含有用的匹配信息.

為此，本文提出一種多尺度圖像子區匹配算法，利用權重函數調整圖像子區內各個像素點對圖像子區相關性的貢獻，使圖像子區的邊界完全變得模糊，實現動態調整圖像子區大小，更加合適大形變或不連續形變的表面，通過實驗證明了該方法能夠得到足夠精確的測量結果.

1 數字圖像相關方法

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

圖1 圖像子區變形示意圖

在一個圖像子區中，假設變形是均勻分布的，即點的位移在水平和垂直方向上都是線性變化的，由形狀函數進行描述：

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

其中xPQ和yPQ分別是P 點距離其圖像子區中心點Q 的水平距離和垂直距離.

由一階泰勒展開得到

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

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

2 多尺度圖像子區匹配算法

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

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

在DIC 中，為了將變形圖像中的一個點與參考圖像中的另一個點進行匹配，僅使用圖像本身的信息是不夠的，因為任何一張圖像中都可能有大量相同灰度值的點.因此才必須考慮圖像子區中包含的周圍信息，根據一定相關準則來區分各個點.

權重函數在評價兩個圖像子區相似性時確定各個點的重要性，其中Δx 和Δy 表示當前點到圖像子區中心點的x 方向距離和y 方向距離.

設f（x，y）表示參考圖像，g（x，y）表示測試圖像，圖像子區的邊長為2M+1，以圖像子區中心點為坐標原點，則最小平方（Least Square，LS）相關準則[33]定義為：

其中，

將其推廣為加權系數：

其中

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

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

對每個要處理的點搜索其最佳尺度因子a.首先設置a 的初始值，計算第一個點時，a 的初始值設置為0，計算其余點時，a 的初始值從前一個已計算的點繼承而來的.然后，對當前點進行并行計算，找出最大加權系數Cw（a）并更新接下來，使a 變大或變小，直到使Cw（a）得到最優值，此時，最后更新就是所求結果.

3 實驗結果與分析

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

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

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

圖3 模擬生成的數字散斑圖像

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

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

表1 模擬DIC 實驗中迭代計算過程中的相關性系數 C（）

表1 模擬DIC 實驗中迭代計算過程中的相關性系數 C（）

迭代次數 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 模擬實驗中迭代計算過程中的相關性系數 C（）的變化趨勢圖

最后得到的計算結果及誤差如表2 所示.

表2 模擬DIC 實驗的計算結果與誤差

表2 同時給出了文獻[31]方法和文獻[32]方法的計算結果作為對比，可以看出，本文方法的計算精度優于另外兩種方法，其中，本文方法的計算結果與文獻[32]方法的計算結果更接近，說明像素加權策略對提高計算精度起主要作用，同時由于本文方法吸收了文獻[31]方法的動態調整策略，兩者結合進一步提高了計算精度.

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

圖5 模擬DIC 實驗的計算結果

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

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

圖6 真實DIC 實驗的計算結果

4 結論

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