數(shù)字體圖像相關(guān)方法的亞體素算法性能研究

代延輝，許 蔚，李 坤

（昆明理工大學(xué) 建筑工程學(xué)院，云南 昆明 650500）

0 引 言

如何精準(zhǔn)獲得材料的力學(xué)性能始終是當(dāng)今研究的熱點(diǎn)問(wèn)題之一，根本原因在于無(wú)法準(zhǔn)確測(cè)量其內(nèi)部變形，傳統(tǒng)的測(cè)量方法（如應(yīng)變片法、散斑干涉法、光彈性法等）具有一定的局限性，導(dǎo)致測(cè)量結(jié)果與實(shí)際變形相差甚大。近年來(lái)，隨著三維層析設(shè)備與計(jì)算機(jī)圖像處理技術(shù)的成熟，許多新型測(cè)量技術(shù)噴涌而現(xiàn)，最具代表性的是數(shù)字體圖像相關(guān)方法（Digital Volume Correlation，DVC）。數(shù)字體圖像相關(guān)方法是一種有效的光學(xué)測(cè)量技術(shù)，具有無(wú)損性、抗干擾能力強(qiáng)、測(cè)量范圍廣等多種優(yōu)勢(shì)，廣泛地應(yīng)用于生物醫(yī)學(xué)、實(shí)驗(yàn)力學(xué)和材料科學(xué)等多學(xué)科交叉領(lǐng)域。

數(shù)字體圖像相關(guān)方法通過(guò)對(duì)比變形前后物體圖像的位置變化，以獲取其內(nèi)部的變形場(chǎng)信息。該方法在計(jì)算過(guò)程中通常分為如下兩步，即整體素位移搜索和亞體素位移測(cè)量。整體素位移搜索可采用空域或頻域搜索以獲得物體內(nèi)部的位移，然而該方式測(cè)得的位移精度不足以達(dá)到預(yù)期理想，因此需要采用亞體素位移算法進(jìn)一步提高測(cè)量精度。亞體素位移測(cè)量算法是將整體素搜索的位移初值不斷進(jìn)行優(yōu)化直至達(dá)到亞體素精度，是提高測(cè)量精度的關(guān)鍵性技術(shù)，對(duì)此研究人員先后提出擬合算法、灰度梯度算法、前向累加牛頓-拉弗森（FA-NR）算法以及三維反向高斯-牛頓（IC-GN）算法等多種算法用于實(shí)現(xiàn)亞體素位移的測(cè)量?；诓煌愋偷乃惴?，其測(cè)量精度和計(jì)算效率各不相同，如何選取最佳算法尤為重要，因此定量對(duì)比多種算法的計(jì)算性能具有十分重要的指導(dǎo)意義。

本文介紹了在數(shù)字體圖像相關(guān)方法中最為常見(jiàn)的三種亞體素位移測(cè)量算法，通過(guò)計(jì)算機(jī)模擬三維散斑體圖像的剛體平移實(shí)驗(yàn)，對(duì)比不同算法的測(cè)量精度和計(jì)算效率，隨后利用椎骨零應(yīng)變?cè)粧呙鑼?shí)驗(yàn)進(jìn)一步對(duì)比三種算法的測(cè)量誤差，實(shí)驗(yàn)均表明IC-GN算法的性能最優(yōu)。

1 數(shù)字體圖像相關(guān)方法的基本原理

采用數(shù)字體圖像相關(guān)方法測(cè)量物體內(nèi)部變形時(shí)，需要結(jié)合三維層析設(shè)備記錄測(cè)試對(duì)象在不同狀態(tài)下的體圖像，通過(guò)匹配技術(shù)尋找這些圖像中灰度相同或相似的體素點(diǎn)，對(duì)比該點(diǎn)在變形前后位置的變化以獲取物體內(nèi)部變形場(chǎng)信息。具體步驟如下：首先在變形前的體圖像中選取一個(gè)尺寸為(2+1)voxel×(2+1)voxel×(2+1)voxel的正方體為參考子體塊，通過(guò)搜索尋找變形圖像中灰度最為相似的目標(biāo)子體塊，并且采用相關(guān)系數(shù)（如最小平方距離函數(shù)、歸一化互相關(guān)函數(shù)等）評(píng)價(jià)兩者的相似程度。由于搜索時(shí)只能以體素為單位，該方式測(cè)得位移無(wú)法達(dá)到亞體素級(jí)別的精度，因此需要借用亞體素位移算法提高測(cè)量精度。

1.1 三元二次多項(xiàng)式擬合算法

數(shù)字體圖像相關(guān)法的擬合算法與二維算法類似，如圖1所示。首先通過(guò)整體素搜索獲得位移初值，選取一個(gè)以整體素點(diǎn)為中心、大小為3×3×3的子體塊為研究對(duì)象，采用三元二次多項(xiàng)式對(duì)該點(diǎn)周圍的幾個(gè)離散點(diǎn)進(jìn)行插值擬合，該多項(xiàng)式極值點(diǎn)對(duì)應(yīng)的坐標(biāo)便是所求點(diǎn)的位置，進(jìn)而計(jì)算出亞體素位移。

圖1 擬合點(diǎn)示意圖

三元二次多項(xiàng)式擬合函數(shù)為：

式中共有10個(gè)未知參數(shù)（～），為了求取其對(duì)應(yīng)的數(shù)值，通常選取等距的12個(gè)離散點(diǎn)依次代入三元二次多項(xiàng)式中，可計(jì)算出所有未知量及多項(xiàng)式表達(dá)式。擬合獲取多項(xiàng)式的極值點(diǎn)應(yīng)滿足下列關(guān)系：

最終計(jì)算得出的亞體素位移為：

1.2 前向累加牛頓-拉弗森（FA-NR）算法

FA-NR算法是通過(guò)逐步非線性迭代優(yōu)化整體素位移的初值直至收斂到最佳值，其實(shí)質(zhì)是將相關(guān)系數(shù)求導(dǎo)問(wèn)題轉(zhuǎn)化為參數(shù)求最優(yōu)值問(wèn)題。FA-NR算法的基本原理與擬合算法截然不同，需要用一階形函數(shù)描述子體塊的形狀變化，其定義如下：

式中：(,,)和(′,′,′)分別表示子體塊變形前后的中心點(diǎn)；，，為中心點(diǎn)在，，方向上的位移變化量；Δ，Δ，Δ為子體塊內(nèi)一點(diǎn)到中心點(diǎn)的距離；u，u，u，v，v，v，w，w，w分別表示子體塊在，，方向上的一階位移梯度。

一階形函數(shù)=[,,,u,u,u,v,v,v,w,w,w]共有12個(gè)參數(shù)，將整體素位移作為迭代的初始值，采用抗噪性能強(qiáng)的零均值歸一化互相關(guān)函數(shù)（ZNCC）作為相關(guān)系數(shù)對(duì)離散數(shù)據(jù)進(jìn)行處理，當(dāng)變形前后的圖像子體塊灰度分布最為相似時(shí)，相關(guān)系數(shù)達(dá)到極值。

零均值歸一化互相關(guān)函數(shù)（ZNCC）如下：

相關(guān)系數(shù)極值點(diǎn)應(yīng)滿足下列關(guān)系：

采用正向迭代的表達(dá)式為：

式中：p為第次迭代值；?(p)為相關(guān)函數(shù)的一階導(dǎo)數(shù) 向 量，即Jacobian矩陣；??(p)為二階導(dǎo)數(shù)，即Hessian矩陣?？刂频氖諗織l件，當(dāng)滿足 |p-p|≤0.001或達(dá)到規(guī)定迭代次數(shù)時(shí)，迭代終止。由于迭代過(guò)程中需要亞體素位置的灰度及灰度梯度信息，通常采用高精度的三元三次插值法對(duì)子體塊的灰度進(jìn)行插值處理以獲得所需信息。

1.3 三維反向高斯-牛頓（IC-GN）算法

IC-GN算法的原理與FA-NR算法類似，但是迭代方向相反，IC-GN每次迭代過(guò)程中以目標(biāo)子體塊為不變量，參考子體塊按照扭曲函數(shù)發(fā)生改變，其迭代參數(shù)以相乘的形式進(jìn)行更新，如圖2所示。IC-GN算法常采用零均值歸一化最小平方距離相關(guān)函數(shù)（ZNSSD）用于評(píng)價(jià)子體塊變形前后相似程度的衡量標(biāo)準(zhǔn)。

圖2 IC-GN算法

零均值歸一化最小平方距離相關(guān)函數(shù)（ZNSSD）為：

與正向算法類似，IC-GN算法達(dá)到收斂條件時(shí)停止迭代。除此之外，IC-GN算法同樣也需要進(jìn)行三元三次插值處理獲得灰度梯度。

2 數(shù)值模擬平移實(shí)驗(yàn)

在實(shí)際應(yīng)用數(shù)字體圖像相關(guān)法測(cè)量位移時(shí)會(huì)受到三維層析圖像采集系統(tǒng)、加載條件、室外環(huán)境等眾多因素的影響，因而采用數(shù)值模擬實(shí)驗(yàn)的效果會(huì)比較理想。文獻(xiàn)[12]將二維計(jì)算機(jī)模擬散斑圖制作的方法擴(kuò)展至三維可有效地模擬出三維散斑體圖包含的灰度信息，該方法分為以下三步：

1）在一定范圍區(qū)域內(nèi)隨機(jī)分布給定數(shù)量的點(diǎn)作為散斑的中心。

2）利用下式計(jì)算出該區(qū)域任意一點(diǎn)灰度：

式中：為散斑顆粒的總數(shù)；是散斑顆粒中心的光強(qiáng)強(qiáng)度且隨機(jī)分布；為散斑顆粒的大??；(x,y,z)表示散斑顆粒中心的位置。

3）將以上所有信息保存為數(shù)字圖像，便可以得到相應(yīng)的散斑體圖像。

通過(guò)設(shè)置以上參數(shù)可獲得散斑數(shù)量為12 000個(gè)、半徑為3 voxel、尺寸為150 voxel×150 voxel×150 voxel的三維散斑體圖像，將該圖像連續(xù)平移10次，每次移動(dòng)0.1 voxel，直到最大位移為1 voxel時(shí)停止，得到10組剛性平移的體圖像，并分別給平移前后的圖像添加均值為0、方差為2的高斯噪聲，如圖3所示。將平移后的圖像作為變形體圖像，利用三種不同的亞體素算法計(jì)算圖像的變形。

圖3 計(jì)算機(jī)制作的散斑體圖像

2.1 三種亞體素算法的精度比較

為了精準(zhǔn)評(píng)價(jià)數(shù)字體圖像相關(guān)方法在平移實(shí)驗(yàn)的測(cè)量誤差，通常將計(jì)算結(jié)果與真實(shí)施加的位移對(duì)比。

均值誤差（系統(tǒng)誤差）為：

標(biāo)準(zhǔn)差（隨機(jī)誤差）為：

式中：u為數(shù)字體圖像相關(guān)方法計(jì)算出第點(diǎn)的位移；為預(yù)加載的真實(shí)位移值；為計(jì)算點(diǎn)的個(gè)數(shù)。

使用三種方法對(duì)變形前后的體圖像進(jìn)行運(yùn)算時(shí)，計(jì)算子體塊的尺寸設(shè)置為21 voxel×21 voxel×21 voxel和41 voxel×41 voxel×41 voxel，計(jì)算步長(zhǎng)均為10 voxel。不同算法對(duì)應(yīng)的測(cè)量誤差如圖4所示，三種算法的均值誤差都呈現(xiàn)出周期為1 voxel、近似于正弦曲線的分布規(guī)律，這與文獻(xiàn)[13]中的誤差規(guī)律基本一致。位移位于0、0.5 voxel和1 voxel時(shí)，測(cè)得的位移平均值最接近于真實(shí)值，均值誤差接近于0；當(dāng)施加的位移為0.2 voxel和0.8 voxel時(shí)，F(xiàn)A-NR算法和IC-GN算法測(cè)量位移的均值與真實(shí)位移相差最大，即均值誤差為最大，而擬合算法在平移0.3 voxel，0.7 voxel對(duì)應(yīng)的均值誤差為最大。其次對(duì)比兩個(gè)不同的計(jì)算子體塊，位移的均值誤差變化甚小，而標(biāo)準(zhǔn)差卻成倍的減少，即計(jì)算子體塊大小的選擇對(duì)數(shù)字體圖像相關(guān)法測(cè)量位移的均值誤差影響不大，對(duì)標(biāo)準(zhǔn)差的影響卻十分顯著。

圖4 不同計(jì)算子體塊對(duì)應(yīng)的測(cè)量誤差

比較三種算法的測(cè)量誤差，均值誤差方面，IC-GN算法在不同尺寸的計(jì)算子體塊中測(cè)量誤差最大為0.006 0 voxel，約為擬合算法的1 5，F(xiàn)A-NR算法測(cè)量的最大值約為其2倍，而擬合算法的測(cè)量誤差最大達(dá)到0.028 2 voxel，這說(shuō)明IC-GN算法的測(cè)量精度要優(yōu)于其他兩種算法。標(biāo)準(zhǔn)差方面，IC-GN算法和FA-NR算法相近，且均小于擬合算法。在子體塊為21 voxel×21 voxel×21 voxel時(shí)，IC-GN算法和FA-NR算法測(cè)得標(biāo)準(zhǔn)差均小于0.002 voxel，而擬合算法最大為0.013 8 voxel，由此可見(jiàn)擬合算法的穩(wěn)定性明顯更低。當(dāng)子體塊增大1倍時(shí)，三種算法的標(biāo)準(zhǔn)差減小了約4倍。

2.2 三種亞體素算法的計(jì)算效率比較

在數(shù)字體圖像相關(guān)方法中，除了要考慮計(jì)算的精度外，計(jì)算效率也是值得考慮的因素之一，隨著計(jì)算子體塊的增大，計(jì)算效率均有不同程度的降低。三種算法均在計(jì)算子體塊為21 voxel×21 voxel×21 voxel時(shí)計(jì)算效率最高，在子體塊為41 voxel×41 voxel×41 voxel時(shí)計(jì)算速度最低。其次擬合算法的計(jì)算效率最高，IC-GN算法次之，F(xiàn)A-NR算法效率最低，如圖5所示，以FA-NR的計(jì)算效率為標(biāo)準(zhǔn)。

圖5 計(jì)算效率比

在子體塊為41 voxel×41 voxel×41 voxel時(shí)，擬合算法的計(jì)算效率是IC-GN的6.41倍，是FA-NR算法的20.34倍。探究其原因可能在于：擬合算法在計(jì)算過(guò)程中，不需要進(jìn)行Hessian矩陣和子體塊灰度插值的計(jì)算，且每次多項(xiàng)式只需計(jì)算一次即可，無(wú)需迭代，因此擬合算法具有較高的計(jì)算效率。FA-NR算法在每次迭代過(guò)程中會(huì)反復(fù)計(jì)算Hessian矩陣，且亞體素位置的灰度和灰度梯度信息需要反復(fù)利用三元三次插值獲得，這將極大消耗計(jì)算時(shí)間和內(nèi)存。而在IC-GN算法的迭代過(guò)程中，由于采用反向迭代的策略，Hessian矩陣只需要計(jì)算一次即可，迭代次數(shù)顯著減少，提高了計(jì)算效率。綜上分析，IC-GN既保證了計(jì)算精度的同時(shí)也提高了計(jì)算效率，F(xiàn)A-NR算法則是以犧牲時(shí)間而獲得較高的測(cè)量精度，而擬合算法的計(jì)算效率最高但是計(jì)算精度最低。

3 零應(yīng)變?cè)粧呙鑼?shí)驗(yàn)

零應(yīng)變?cè)粧呙璧膶?shí)驗(yàn)對(duì)象為正常成年男性腰椎標(biāo)本，由昆明醫(yī)科大學(xué)捐贈(zèng)，如圖6所示。測(cè)試時(shí)，將實(shí)驗(yàn)樣品放置在Micro CT-225上進(jìn)行零應(yīng)變?cè)粧呙?，掃描期間樣品沒(méi)有任何移動(dòng)，掃描時(shí)長(zhǎng)約為15 min，并將掃描圖像以39μm的各向同性體素尺寸重建。CT設(shè)備掃描試件時(shí)并沒(méi)有發(fā)生移動(dòng)，理論上物體沒(méi)有發(fā)生變形，然而由于多種因素的影響，如CT射線熱誤差、算法誤差、體圖像質(zhì)量等，必然導(dǎo)致測(cè)量結(jié)果有少許偏差，且這些誤差并不能完全消除。由于運(yùn)算時(shí)采用相同的體圖像，測(cè)量結(jié)果可直接用于對(duì)比三種算法的測(cè)量誤差大小。

圖6 椎骨體圖像及其灰度分布圖

利用數(shù)字體圖像相關(guān)方法測(cè)量椎骨變形時(shí)，計(jì)算子體塊大小選擇為21 voxel×21 voxel×21 voxel，計(jì)算步長(zhǎng)設(shè)置為5 voxel。三種算法的測(cè)量誤差如圖7所示，其中擬合算法的誤差最大，F(xiàn)A-NR算法和IC-GN算法的誤差相對(duì)較小且基本相同。

圖7 三種算法的測(cè)量均值絕對(duì)誤差和標(biāo)準(zhǔn)差

雖然算法不同，但是誤差最值出現(xiàn)的方向一致，即均值絕對(duì)誤差最大值在軸，標(biāo)準(zhǔn)差則出現(xiàn)在軸方向。在計(jì)算精度方面，IC-GN算法精度最高，其均值絕對(duì)誤差在0.057 5～0.146 9 voxel范圍內(nèi)，F(xiàn)A-NR算法的誤差為0.058 9～0.148 5 voxel，而擬合算法精度最低，其誤差范圍為0.063 1～0.159 7 voxel。在計(jì)算穩(wěn)定性方面，IC-GN算法的表現(xiàn)最優(yōu)，其標(biāo)準(zhǔn)差為0.057 5～0.146 9 voxel，F(xiàn)A-NR算法的標(biāo)準(zhǔn)差為0.059 6～0.014 85 voxel，擬合算法穩(wěn)定性最差，其標(biāo)準(zhǔn)差范圍為0.062 1～0.155 2 voxel。需要特別注意的是，三種算法在零應(yīng)變實(shí)驗(yàn)中的誤差最值要明顯大于平移實(shí)驗(yàn)中的最值，由此可見(jiàn)零應(yīng)變實(shí)驗(yàn)的誤差主要是由外部因素引起，而算法誤差所占比例相對(duì)較小。綜合分析，F(xiàn)A-NR算法和IC-GN算法在椎骨原位掃描下測(cè)量的精度及穩(wěn)定性也均優(yōu)于擬合算法。

4 結(jié) 論

本文簡(jiǎn)要地介紹了三種常見(jiàn)的亞體素位移測(cè)量算法：擬合算法、FA-NR算法和IC-GN算法，通過(guò)兩組不同的實(shí)驗(yàn)對(duì)比三種算法的計(jì)算精度及其穩(wěn)定性，結(jié)果顯示IC-GN算法的性能最優(yōu)。在散斑體圖像的平移實(shí)驗(yàn)中，三種算法的計(jì)算精度均隨著子體塊的增大而提高，但是計(jì)算效率卻成倍的遞減。其中IC-GN算法的測(cè)量精度最高，而擬合算法的測(cè)量誤差和穩(wěn)定性均為最差，而其計(jì)算效率最高。在椎骨零應(yīng)變測(cè)試中三種算法表現(xiàn)的性能也不盡相同，IC-GN算法和FA-NR算法的測(cè)量精度均優(yōu)于擬合算法。其次對(duì)比兩組實(shí)驗(yàn)可發(fā)現(xiàn)平移實(shí)驗(yàn)的測(cè)量精度要高于零應(yīng)變測(cè)試的精度，主要原因是外部因素引起的誤差要大于算法誤差。綜合對(duì)比計(jì)算精度和計(jì)算效率，可確定IC-GN算法為數(shù)字體圖像相關(guān)方法中的最佳亞體素位移測(cè)量算法。

注：本文通訊作者為許蔚。