當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

Matlab实现二维数字图像相关（2D Digital Image Correlation, 2D-DIC）【ADIC2D代码复现及原理介绍】

發布時間：2023/12/31 编程问答 36 豆豆

生活随笔收集整理的這篇文章主要介紹了 Matlab实现二维数字图像相关（2D Digital Image Correlation, 2D-DIC）【ADIC2D代码复现及原理介绍】小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

前言

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

本文所復現的論文——Atkinson D, Becker T. A 117 Line 2D Digital Image Correlation Code Written in MATLAB[J]. Remote Sensing, 2020, 12(18): 2906.¹
GitHub源碼地址：https://github.com/SUMatEng/ADIC2D.git

由于官方數據集可能存在變動或者被墻，這里提供一下我已經上傳好的版本
數據集地址：

百度云：鏈接：https://pan.baidu.com/s/1In4ZAFRIzfWIxoBSMAYGYQ
提取碼：8xmi
谷歌云盤：https://drive.google.com/file/d/1JUTGmZirxXHYg4fPno-wE3YUmyTqWKzI/view?usp=sharing

數字圖像相關專欄目錄：

Matlab實現二維數字圖像相關（2D Digital Image Correlation, 2D-DIC）【ADIC2D代碼復現及原理介紹】

數字圖像相關（Digital Image Correlation, DIC）中的非線性優化方法（FA-GN與IC-GN）

數字圖像相關（Digital Image Correlation, DIC）中的非線性優化方法IC-GN的數值解計算

用MATLAB繪制隨機散斑圖案【源碼+正確的橢圓旋轉公式】

這篇論文是讓我真正想要將自己正在做的這些事進行記錄的源動力，因為這篇文章的作者也是考慮到盡管數字圖像相關技術很成熟，但是網絡上缺乏可以實操的源碼便于理解，從而將簡化的代碼以及對于原理的解釋共享在網絡上便于他人學習。非常感謝Devan Atkinson，Thorsten Becker兩位作者對于他們知識成果的共享。我將在這里將我學到的東西以中文的形式以方便國內的學者進行進一步的研究。

This paper aims to bridge the gap between the theory and implementation of DIC. It does this by firstly presenting the theory for a 2D, subset based DiC framework that is predominantly consistent with current state-of-the-art practices. ¹

一.數字圖像相關（Digital Image Correlation）

數字圖像相關技術(Digital Image Correlation, DIC) 是廣泛應用于剛體結構形變檢測的視覺測量技術。它通過在被測物表面投影或繪制隨機散斑圖案，定義圖像的相似度函數，對物體形變前后的兩幅圖像進行分析即可得到采樣點的位移場，從而得到形變²。數字圖像相關技術按照其測量的維數，可以被分為二維數字圖像相關（Two-Dimensional Digital Image Correlation, 2D-DIC）和三維數字圖像相關（Three-Dimensional Digital Image Correlation, 3D-DIC）³。如圖所示⁴就是2D-DIC獲取形變后子區位移場的過程。

對于數字圖像相關而言，其流程可以分為標定、相關運算、位移變換以及最后根據實際需要計算應變或者其他信息四個步驟，標定就是通過一系列拍攝的特定圖案的圖片，利用標定算法求解出相機圖像坐標系與真實世界坐標系下的映射關系，之后我單獨寫一篇來分享好了，這里主要分享一下相關運算的流程與推導。

二.相關運算

數字圖像相關的本質就是在參考圖像上設置好子區后，找到形變后圖片上對應的已經產生過位移的子區，從而計算出形變子區相對于參考子區在圖像坐標系上的位移，最終利用標定好的映射關系求解出該區域在實際物理空間中的位移，最后利用得到的位移計算得到所拍攝表面所發生的形變。在這個過程中參考子區是已知的，但形變后的圖片上原本的子區必然發生了變化，其具體位置是未知的，因此如何才能找到和原先參考子區對應的那塊區域（求解這一未知量）正是關鍵所在，這正是相關運算所做的事情。

1.數學模型

我強烈建議大家能一起推導一次這一部分，對于理解相關運算有很大的作用！！

首先設參考圖像為 $F$ ，它很多時候代表了 $t = 0$ 時刻所拍攝到的散斑圖案，其上的參考子區設為 $f$ 。

子區的排布有很多種，可以加步進那樣整列排布，也可以隨機的點選；子區的形狀也可以各種各樣，矩形、圓型，只要保證相關運算所處理的所有子區為同一種形狀、大小的即可。

接下來設形變圖像為 $G$ ，其上的形變子區為 $g$ ，如下為一個在相機圖像坐標系下的示意圖¹

圖中左上角的矩形框即為參考子區，右下角的平行四邊形即為形變子區，其中參考子區 $f$ 中心點為 $x0=[x0y0]T\boldsymbol{x}^{0}=\left[\begin{array}{ll}x_{0} & y_{0}\end{array}\right]^{T}$ （后面的公式表達請注意上下標），參考子區上的其他像素點為 $xi=[xiyi]T\boldsymbol{x}^{i}=\left[\begin{array}{ll}x_{i} & y_{i}\end{array}\right]^{T}$ 。以中心點 $x0\boldsymbol{x}^{0}$ 作為參考，我們可以得到
$xi=Δxix0+xo=[xiyi]=[ΔxiΔyi]+[x0y0]\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]$ 其中 $Δxix0\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}$ 表示 $x0\boldsymbol{x}^{0}$ 到 $xi\boldsymbol{x}^{i}$ 的坐標偏移。現在我們設置形變子區中與之前 $xi\boldsymbol{x}^{i}$ 相對應的當前坐標點為 $xi′=[xi′yi′]T\boldsymbol{x}^{i'} =\left[\begin{array}{ll}x_{i}' & y_{i}'\end{array}\right]^{T}$ ，與之前中心點 $x0\boldsymbol{x}^{0}$ 相對應的形變子區中心點為 $xd\boldsymbol{x}^ze8trgl8bvbq$ 。與上式同理，我們同樣可以得到 $xi′\boldsymbol{x}^{i'}$ 的表達式為
$xi′=Δxi′x0+xo=[xi′yi′]=[Δxi′Δyi′]+[x0y0]\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′x0\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}$ 表示 $x0\boldsymbol{x}^{0}$ 到 $xi′\boldsymbol{x}^{i'}$ 的坐標偏移。注意，沒有使用形變子區中心點作為參考點，是有兩個原因的，第一個原因是由于本身 $xd\boldsymbol{x}^ze8trgl8bvbq$ 對于我們來說依然是未知點，它只不過是一個 $xi′\boldsymbol{x}^{i'}$ 特例罷了，這里需要用一個已知點作為參考，還有一個原因會在形函數推導那一部分提到，這里先不引入。

在這個式子當中，最重要的就是計算這里的 $Δxi′x0\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}$ ，因為在獲得這個值后就找到形變后與之前子區一一對應的像素點坐標，從而可以計算整個子區圖像的位移場，最后即可計算該區域的形變或者受力了。而這個未知量的求解，需要引入一個有限元分析中的概念——形函數（Shape Function）。對于一個區域或者一個單元，其形變是可以通過數學模型來進行表示的，而這個用于描述的數學模型就是形函數（圖像處理中常見的仿射變換實質就是一個形函數）。形函數的階數越高，則可以描述的形變就越復雜。因此對于 $Δxi′x0\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}$ 是可以用形函數來表達的，這里設形函數關系為 $W$ ，形函數參數（Shape Function Parameters, SFPs）為 $P$ ，我們有
$Δxi′x0=W(Δxix0,P)\Delta \boldsymbol{x}^{i'} \boldsymbol{x}^{0}=W(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0},P)$
那么現在關系式就全部建立完成了，接下來，如何建立這個形函數以及如何知道最合適的形函數參數成為了下一步的目標，首先是形函數的建立，這需要根據實際情況選擇合適階數的形函數關系，留到下一章進行探討。那么如何知道最合適的形函數參數呢？這就需要用到數字圖像相關最核心的內容——相關運算。在信號與系統中，評判兩個信號的相似程度就是對他們進行相關運算，二維圖像實質也是一種信號，因此評判參考子區和形變子區是否相似也是可以用相關運算的。

這樣的話我們就可以將兩者相關運算作為一個目標函數（其結果為相關標準，用于衡量兩者的相似性， Correlation Criterion），通過建立好的形函數關系式以及一組形函數參數初值 $P_{0}$ ，通過對相關標準的梯度下降法迭代求解，即可以找到一組最合適的形函數參數最優解 $P^{*}$ ，此時所得到的形變子區就是與真實的形變情況最為接近的，從而實現了參考子區與形變子區之間的匹配！

二維數字圖像本身是對于光強的一種采樣，即本身是離散點，考慮到匹配的精度，需要對形變圖像進行插值，以實現亞像素級別的匹配。經過插值運算后，參考子區和形變子區都可以寫成函數的形式，即
$fi=F(xo+Δxix0)f_{i}=F\left(\boldsymbol x^{o}+\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}\right)$ $gi=G(xo+W(Δxix0,P))g_{i}=G\left(\boldsymbol x^{o}+\boldsymbol{W}\left(\Delta \boldsymbol{x}^{i} \boldsymbol{x}^{0}, \boldsymbol{P}\right)\right)$ 將上面的推導進行歸納，整個運算的實質就是通過給定一組形函數參數初值 $P_{0}$ ，對目標函數（相關標準）不斷迭代優化從而得到一組最優解 $P^{*}$ 。

2.形函數

在二維數字圖像相關中，最常采用的形函數有0階形函數（ $W^{SF0}$ ）、1階形函數（ $W^{SF1}$ ）和2階形函數（ $W^{SF2}$ ）。階數越高，能描述的形變就越復雜，但待求解的的位置量也會越多。如圖所示的就是這三種形函數的示意圖¹

0階形函數：只能描述位移變化，無法描述拉伸或者扭曲（速度較快，適合用于小形變系統的實時檢測）
$WSF0(Δxix0,PSF0)=[10u01v][ΔxiΔyi1]\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]$
1階形函數：能描述子區位移和拉伸（仿射變換），但不能描述扭曲（適合用于小形變系統）
$WSF1(Δxix0,PSF1)=[1+uxuyuvx1+vyv][ΔxiΔyi1]\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]$
2階形函數：能描述子區位移、拉伸和扭曲（適合用于大形變系統）
$WSF2(Δxix0,PSF2)=[12uxxuxy12uyy1+uxuyu12vxxvxy12vyyvx1+vyv][Δxi2ΔxiΔyiΔyi2ΔxiΔyi1]\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]$ 由上述式子我們可以得到不同階數形函數需要求解的參數為：
$PSF0=[uv]T\boldsymbol{P}^{S F 0}=\left[\begin{array}{llllllllllll} u & v\end{array}\right]^{T}$ $PSF1=[uuxuyvvxvy]T\boldsymbol{P}^{S F 1}=\left[\begin{array}{llllllllllll} u & u_{x} & u_{y} & v & v_{x} & v_{y} \end{array}\right]^{T}$ $PSF2=[uuxuyuxxuxyuyyvvxvyvxxvxyvyy]T\boldsymbol{P}^{S F 2}=\left[\begin{array}{llllllllllll} u & u_{x} & u_{y} & u_{x x} & u_{x y} & u_{y y} & v & v_{x} & v_{y} & v_{x x} & v_{x y} & v_{y y} \end{array}\right]^{T}$ 那這里的每一項參數代表什么呢？從上一節的內容，我們可以知道，對于參考子區中任意一個像素點，其對應的形變子區上的映射點坐標都可以寫成參考子區對應點的坐標值+ $X / Y$ 方向的偏移量的形式。同時，由于兩個點本身相關，因此這個偏移量和參考子區對應點坐標是有關系的。因此，我們可以得到以下關系式，設參考子區上的某一點 $f (x, y)$ ，設其在形變子區上的映射點為 $g(x~,y~)g(\widetilde{x},\widetilde{y})$ ，我們有：
$x~=x+u(x,y)\widetilde{x}=x+u(x,y)$ $y~=y+v(x,y)\widetilde{y}=y+v(x,y)$ 對兩式進行二階泰勒展開⁵，我們可以獲得
$x~=x0+u+uxΔx+uyΔy+12uxxΔx2+12uyyΔy2+uxyΔxΔy\tilde{x}= x_{0}+u+u_{x} \Delta x+u_{y} \Delta y+\frac{1}{2} u_{x x} \Delta x^{2} +\frac{1}{2} u_{y y} \Delta y^{2}+u_{x y} \Delta x \Delta y$ $y~=y0+v+UxΔx+vyΔy+12vxxΔx2+12vyyΔy2+vxyΔxΔy\tilde{y}= y_{0}+v+U_{x} \Delta x+v_{y} \Delta y+\frac{1}{2} v_{x x} \Delta x^{2} +\frac{1}{2} v_{y y} \Delta y^{2}+v_{x y} \Delta x \Delta y$ 而形函數的參數就是展開式中各項的系數，值得注意的是，展開式中的常數項 $x_{0},y_{0})$ ，考慮到每一個子區只會計算一個形函數，因此單一子區內各點的計算均要使用相同的常數項，這也正是上一節中每個點均以參考子區中心點坐標作為參考的原因了。

在實際進行相關運算的過程中，需要根據實際情況選擇合適的形函數，形函數選擇的越好則計算得到的形變子區越接近真實情況，計算的值也會更加準確！

3.相關標準

相關標準在相關運算中相當于最后的評價函數，選擇合適的相關標準對于迭代優化的速度有很大作用。一個評價函數的好壞我個人認為是要滿足快收斂、足夠“凸”、“抗噪聲”，快收斂使得整體優化能盡快逼近最優解。足夠“凸”使得每次迭代優化步進都會很大，同時最優解點確實夠好。而抗噪聲則確保了所逼近的最優解點和實際情況相符，而不是由于噪聲導致的異常值，減小誤差。這一塊我自己還不太能吃透，就直接搬結果吧。

在Bing Pan等人的文章Equivalence of digital image correlation criteria for pattern matching⁶中介紹了各種相關標準的算術推導以及相互之間的轉換，其中最好的三種相關標準為零均值歸一化互相關標準（Zero-Normalized Cross-Correlation Criterion, ZNCC）、零均值歸一化最小距離平方標準（Zero-Normalized Sum of Squared Differences Criterion, ZNSSD）和最優化參數的距離平方標準（Parametric Sum of Squared Difference, PSSDab）。⁷

這三種標準相互之間可以互相轉換，在對最終結果的推導上是一致的。Bing Pan等人還在文章中測試分別使用PSSD與ZNSSD作為相關標準時候的計算效率，PSSD會更好一些。當然選用哪種標準還是要看實際情況來選擇。

其他知識

考慮篇幅（我太懶了而且沒搞太懂），相關運算的其他知識點我就在這里提一嘴，之后有空再單獨寫帖子補充吧。

由于需要對形變圖像進行精確插值，目前最常采用的插值運算為雙三次b樣條插值（bicubic b-spline interpolation）

上述的三種相關標準，對于高頻噪聲都比較敏感，因此需要對參考圖像和形變圖像都進行高斯濾波，濾除高頻噪聲（具體的參數需要研究，由于低通濾波器會將很多紋理模糊）

基于相關標準來優化形函數參數實質是個非線性優化問題，目前最常采用的計算方式為逆合成高斯-牛頓法——IC-GN算法（Inverse compositional Gauss-Newton method），推薦閱讀數字圖像相關（Digital Image Correlation, DIC）中的非線性優化方法（FA-GN與IC-GN）

優化的開始需要提供一組形函數參數初值，這組初值越接近最優解，迭代速度提升也就越快，同時其合理性也決定了所得最優解與真實值之間的誤差（因為優化等于一個尋找峰頂的過程，我們想要的是最高的峰頂，但如果初值位置不位于最高峰，則找到的峰頂并不是最理想的）。

三.ADIC2D代碼解釋

在github上下載好源碼，根據readme下載好需要的數據集，我們即可完成對于ADIC2D項目的部署，注：我是在Matlab 2018a上進行的調試。如圖，是部署好后的各文件含義

整個項目的調用關系及框架如圖所示

部署完成之后直接運行runme.m文件即可，數據集的選擇、ADIC2D算法輸入參數均可在文件開頭進行修改。具體可修改參數如下：

% runme.m參數輸入部分解釋 Sample=3; % EN:which sample to analyse CH:選擇使用的數據集序號1——3 FileNames=ImageNamesLoad(Sample); % load appropriate images for chosen sample Mask=ones(size(im2double(imread(FileNames{1})))); % EN:define a mask that includes all pixels of the image CH:設定一個蒙版決定使用圖片中的那些區域，即ROI% Function inputs (which set up the DIC analysis) GaussFilter=[0.4,5];% 高斯濾波參數[sigma, kernal size]，用于濾除高頻噪聲 StepSize=5; % 子區步進 SubSize=41; % 子區大小 SubShape='Circle'; % 子區形狀 SFOrder=1; % 形函數階數（0,1,2） RefStrat=0; % 子區開始編號 StopCritVal=1e-4; % 停止迭代的標準 WorldCTs=0; % 世界坐標系下的標定點位置 ImageCTs=0; % 圖像坐標系下的標定點位置 rho=0; % 標定板厚度（mm）

之后runme.m會調用ADIC2D算法函數將讀入的數據集進行處理，之后會輸出處理好的數據集結構體。runme.m的后半部分是利用輸出的數據結果進行顯示和計算誤差，包括某副圖像子區沿 $X / Y$ 方向的位移情況、位移結果的均方根誤差等。這部分代碼都比較易讀，直接說一下ADIC2D及其子函數的算法邏輯吧

1.ImgCorr

ImgCorr函數就是按照上一個大章節中的思路來書寫的，下面是它的算法概要

2.SubCorr

SubCorr函數是利用IC-GN的方法來優化形函數參數，直到得到最優解，從而輸出最佳匹配的形變子區信息，具體代碼邏輯請參照論文¹及其源碼，這里僅展示其算法框架圖

四.寫在最后

這是我在自己做課題過程中讀到的一篇寫的很通透的論文，數字圖像相關其實是很成熟的技術了，但是就目前看來在國內，這一塊的研究還很少，即使是有相關的應用也大都是照搬理論，并沒有對于這種技術詳細的描述和深入的解釋。本人也是個小白，希望能通過這篇帖子將這種技術說的稍微清楚些，其中如果有任何理解上的偏差和錯誤還望業界的大佬能多多指正，確保輸出的知識足夠準確和客觀。之后我把自己的課題理清楚后，應該會寫一個小的demo作為分享，方便相關學者們能更快的理解和使用這種技術。最后非常感謝Atkinson D, Becker T對于這份知識的分享。

推薦個應用該技術的網站，方便大家進一步深入
http://www.ncorr.com/index.php/dic-algorithms

參考引用

Atkinson D, Becker T. A 117 Line 2D Digital Image Correlation Code Written in MATLAB[J]. Remote Sensing, 2020, 12(18): 2906. ?? ?? ?? ?? ??

Chu T C, Ranson W F, Sutton M A. Applications of digital-image-correlation techniques to experimental mechanics[J]. Experimental mechanics, 1985, 25(3): 232-244. ??

周軼昊. 基于雙目視覺的物體表面三維復雜運動重建及其應用[D].復旦大學,2012. ??

Khoo S W, Karuppanan S, Tan C S. A review of surface deformation and strain measurement using two-dimensional digital image correlation[J]. Metrology and Measurement Systems, 2016, 23(3). ??

Lu H, Cary P D. Deformation measurements by digital image correlation: implementation of a second-order displacement gradient[J]. Experimental mechanics, 2000, 40(4): 393-400. ??

Pan B, Xie H, Wang Z. Equivalence of digital image correlation criteria for pattern matching[J]. Applied optics, 2010, 49(28): 5501-5509. ??

Sutton M A, Orteu J J, Schreier H. Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[M]. Springer Science & Business Media, 2009. ??

總結

以上是生活随笔為你收集整理的Matlab实现二维数字图像相关（2D Digital Image Correlation, 2D-DIC）【ADIC2D代码复现及原理介绍】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：常见封装及参数
下一篇： html图像设计代码,html——图像设