宋湘芬，林仕令，宋磊，李文

（廣西醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院，廣西南寧530021）

1 引 言

隨著醫(yī)學(xué)圖像學(xué)的不斷發(fā)展，醫(yī)學(xué)圖像配準(zhǔn)已經(jīng)成為醫(yī)學(xué)圖像處理技術(shù)的研究熱點(diǎn)和難點(diǎn)。醫(yī)學(xué)圖像配準(zhǔn)方法可以分為剛性配準(zhǔn)和彈性配準(zhǔn)，當(dāng)前國(guó)內(nèi)外對(duì)剛性配準(zhǔn)的研究已經(jīng)成熟，而彈性配準(zhǔn)仍在發(fā)展階段。實(shí)際上，臨床上醫(yī)學(xué)圖像的許多形變和偏移需要彈性配準(zhǔn)來(lái)實(shí)現(xiàn)。例如在研究腹部、胸部臟器的圖像配準(zhǔn)時(shí)，不自主的生理運(yùn)動(dòng)或患者移動(dòng)等使其內(nèi)部的器官和組織的位置、尺寸和形狀發(fā)生改變，或者由于手術(shù)干涉產(chǎn)生的組織或器官的變形，這些都是剛性配準(zhǔn)所不能解決的，需要彈性配準(zhǔn)來(lái)實(shí)現(xiàn)，因此，彈性配準(zhǔn)有很高的臨床實(shí)用價(jià)值［1］。而B(niǎo)樣條變換因?yàn)槠浜芎玫木植靠刂菩?，可以?duì)局部區(qū)域的變形進(jìn)行更好的逼近，避免了進(jìn)行全局計(jì)算的重復(fù)進(jìn)行，節(jié)省了時(shí)間，因而，基于B樣條的彈性配準(zhǔn)方法受到了廣泛的關(guān)注?；贐樣條的彈性配準(zhǔn)方法主要有兩種，一種是以求得網(wǎng)格控制點(diǎn)灰度值為主，通過(guò)調(diào)整網(wǎng)格控制點(diǎn)值來(lái)擬合已知的插值點(diǎn)；另一種是由Lee等最先提出的基于B樣條FFD模型［2－3］，該模型由于其優(yōu)越的局部形變的建模能力、不受限制的自由度的特點(diǎn)能夠描述復(fù)雜的幾何形體［4］。

本研究采用后者。該方法通過(guò)操縱套用在圖像區(qū)域上的控制點(diǎn)網(wǎng)格來(lái)模擬彈性形變，然后通過(guò)優(yōu)化算法優(yōu)化網(wǎng)格控制點(diǎn)坐標(biāo)位置的方法來(lái)完成彈性配準(zhǔn)。

2 配準(zhǔn)流程

配準(zhǔn)流程圖見(jiàn)圖1。

圖1 配準(zhǔn)流程圖Fig 1 Registration flow diagram

2.1 待配準(zhǔn)圖像網(wǎng)格化

2.1.1 建立控制網(wǎng)格 在待配準(zhǔn)圖像上建立控制網(wǎng)格，見(jiàn)圖2（以11×11控制網(wǎng)格為例），相對(duì)控制網(wǎng)格，定義像素點(diǎn)坐標(biāo)位置為（i，j），則像素空間（m×n控制網(wǎng)格）表示為：

控制空間（（m＋3）×（n＋3）控制網(wǎng)格）表示為：

圖211 ×11控制網(wǎng)格Fig 2 Congtrolmesh of 11×11

記單個(gè)控制點(diǎn)為 φi，j，共（m＋3）×（n＋3）個(gè)控制點(diǎn)，每個(gè)控制點(diǎn)具有水平和垂直方向上的自由度，所以 Φ的維數(shù)（自由度）是N＝2（m＋3）（n＋3），為了簡(jiǎn)單起見(jiàn)，可將Φ二維矩陣位置轉(zhuǎn)換成一維矩陣位置，即：

2.1.2 建立像素網(wǎng)格 設(shè)圖像在X和Y方向的像素點(diǎn)數(shù)量分別為δxm和δyn，其中，δx和δy分別表示單位控制網(wǎng)格在X和Y方向的間距，相對(duì)像素網(wǎng)格，像素點(diǎn)坐標(biāo)位置為（δxi，δyj），則：

像素空間表示為：

控制空間表示為：

基于B樣條的空間變換需要用到控制點(diǎn)相對(duì)像素網(wǎng)格的坐標(biāo)位置作為初始參數(shù)。

2.1.3 控制網(wǎng)格坐標(biāo)范圍的確定原因 由B樣條理論可知，控制點(diǎn)移動(dòng)范圍由選擇的B樣條次數(shù)決定，不同的B樣條函數(shù)的支撐區(qū)域不一樣。本研究選擇三次B樣條，因此，每一個(gè)控制點(diǎn)在其最近鄰4×4網(wǎng)格點(diǎn)上移動(dòng)，見(jiàn)圖2（圖中區(qū)域3），一個(gè)控制點(diǎn)移動(dòng)將引起最近鄰4×4網(wǎng)格內(nèi)所有的像素點(diǎn)發(fā)生移動(dòng)［5］，見(jiàn)圖2（圖中區(qū)域4），一個(gè)像素點(diǎn)移動(dòng)后的坐標(biāo)位置由最近鄰4×4個(gè)控制點(diǎn)的坐標(biāo)位置和該像素點(diǎn)在單元控制網(wǎng)格內(nèi)的相對(duì)位置共同決定。為了確保每個(gè)像素點(diǎn)都有對(duì)應(yīng)的4×4個(gè)控制點(diǎn)，即像素空間左邊界的單元控制網(wǎng)格內(nèi)的每個(gè)像素點(diǎn)必須要有X坐標(biāo)為－1的4個(gè)控制點(diǎn)，下邊界的單元控制網(wǎng)格內(nèi)的每個(gè)像素點(diǎn)必須要有Y坐標(biāo)為－1的4個(gè)控制點(diǎn)，右邊界的單元控制網(wǎng)格內(nèi)的每個(gè)像素點(diǎn)必須要有X坐標(biāo)為（m＋1）的4個(gè)控制點(diǎn)，上邊界的單元控制網(wǎng)格內(nèi)的每個(gè)像素點(diǎn)必須要有Y坐標(biāo)為（n＋1）的4個(gè)控制點(diǎn)。究其控制網(wǎng)格在像素空間外向X、Y負(fù)方向各擴(kuò)充一個(gè)坐標(biāo)，卻向X、Y正方向各擴(kuò)充兩個(gè)坐標(biāo)的原因，是因?yàn)閄坐標(biāo)為（m＋1）、Y坐標(biāo)為（n＋1）的控制點(diǎn)都屬于鄰近右上方向的單元控制網(wǎng)格內(nèi)，所以必須要在向X、Y正方向各擴(kuò)充一個(gè)坐標(biāo)的基礎(chǔ)上再擴(kuò)充一個(gè)坐標(biāo)。

2.2 基于B樣條的自由形變模型

基于B樣條的FFD模型的參數(shù)變量是各控制點(diǎn)的坐標(biāo)位置，該模型的主要思想是移動(dòng)控制點(diǎn)來(lái)實(shí)現(xiàn)局部變形，而控制點(diǎn)以外的像素點(diǎn)移動(dòng)后的坐標(biāo)位置通過(guò)B樣條擬合來(lái)得到［6］。除待配準(zhǔn)圖像網(wǎng)格化后得到控制點(diǎn)的初始坐標(biāo)位置外，控制點(diǎn)移動(dòng)后的坐標(biāo)位置通過(guò)優(yōu)化算法更新。了解這些控制點(diǎn)的移動(dòng)后的坐標(biāo)位置后，基于B樣條的FFD模型，求出圖像中每個(gè)像素點(diǎn)的移動(dòng)后的坐標(biāo)位置。在這里把圖像中像素點(diǎn)的移動(dòng)分解為X方向和Y方向，分別求出移動(dòng)后的X和Y坐標(biāo)位置。

對(duì)于任意一像素點(diǎn)（x，y），經(jīng)三次B樣條彈性變形后的坐標(biāo)位置（相對(duì)像素網(wǎng)格）可表示為：

其中：φi＋l，j＋m表示最近鄰 4×4個(gè)控制點(diǎn)的坐標(biāo)位置，i、j分別表示最近鄰4×4個(gè)控制點(diǎn)的單元索引號(hào)，i＝［x／δx］－1，j＝［y／δy］－1；［］表示取下整數(shù)；u、v分別表示x、y的相對(duì)單元控制網(wǎng)格的位置，u＝x／δx－［x／δx］，v＝y(tǒng)／δy－［y／δy］；

B0（u）～B3（u）表示三次B樣條基函數(shù)，其表達(dá)式為：

其中，0≤u＜1，這些函數(shù)起到權(quán)函數(shù)的作用，它們根據(jù)控制點(diǎn)到（x，y）的距離大小來(lái)加權(quán)每個(gè)控制點(diǎn)對(duì) T（x，y）的貢獻(xiàn)［7］。

其模型構(gòu)造見(jiàn)圖2（見(jiàn)圖中區(qū)域4）。

2.3 雙線(xiàn)性插值

雙線(xiàn)性插值是將待配準(zhǔn)圖像的像素點(diǎn)（x，y）經(jīng)B樣條變換后的坐標(biāo)位置（Tx，Ty）作為插值點(diǎn)，利用該插值點(diǎn)最近鄰4個(gè)像素點(diǎn)的灰度值，計(jì)算該坐標(biāo)位置的灰度值，并賦給配準(zhǔn)后圖像的對(duì)應(yīng)位置（x，y），見(jiàn)圖 3。具體插值過(guò)程［8］如下：

圖3 雙線(xiàn)性插值示意圖Fig 3 The sketch map of bilinear interpolation

圖4 雙線(xiàn)性插值原理圖Fig 4 The illustrative diagram of bilinear interpolation

如圖4所示，設(shè)插值點(diǎn)（Tx，Ty）的最近鄰4個(gè)像素點(diǎn)為 A、B、C、D，其坐標(biāo)分別為（i，j）、（i＋1，j）、（i，j＋1）、（i＋1，j＋1），其灰度值分別是 g（A）、g（B）、g（C）、g（D）。首先計(jì)算 E和 F兩點(diǎn)的灰度值g（E）、g（F）：

則插值點(diǎn)（Tx，Ty）的灰度值 g（Tx，Ty）為：

2.4 SSD測(cè)度

SSD測(cè)度公式［9］見(jiàn)式（11）：

其中，N為圖像配準(zhǔn)區(qū)域內(nèi)像素總數(shù)，I2（x，y）和I1（x，y）分別表示二維空間中的參考圖像和待配準(zhǔn)圖像的灰度函數(shù)，當(dāng)兩幅圖像達(dá)到配準(zhǔn)狀態(tài)時(shí)，它們之間Essd最小。

2.5 梯度下降法［10］

根據(jù)式（6），求解最優(yōu)Φ為最主要的步驟，其主要思想是尋求使相似性測(cè)度Essd最小的空間變換位置。求解最優(yōu)Φ是通過(guò)先求解Essd的梯度▽Essd，然后通過(guò)梯度下降法對(duì)Φ更新來(lái)實(shí)現(xiàn)的，Φ的更新公式見(jiàn)式（12）：

其中，stepsize為步長(zhǎng)，▽Essd是對(duì)控制點(diǎn)的每個(gè)網(wǎng)格位置處的變形的梯度。這里步長(zhǎng)也要在速度和精度之間取得折中，通過(guò)變步長(zhǎng)法實(shí)現(xiàn)由粗到細(xì)的快速精確配準(zhǔn)，變形模型越精細(xì)，步長(zhǎng)越小。

最優(yōu)Φ的求解過(guò)程如下：

（1）取步長(zhǎng)stepsize為0.01，初始化控制參數(shù)Φ，初始化迭代次數(shù)iter，求解初始Essd和初始▽Essd；

（2）根據(jù)初始Φ和初始▽Essd，找到最佳初始步長(zhǎng) step＿init（1e－4～1e＋3），更新初始 Essd，保存當(dāng)前參數(shù)，以備下一次迭代用；

（3）設(shè)置 while循環(huán)條件（‖▽Essd‖≥le－4），滿(mǎn)足，則進(jìn)入循環(huán)，根據(jù)當(dāng)前參數(shù)，尋找最佳步長(zhǎng)step＿best，尋找 step＿best失敗則退出循環(huán)，否則將step＿best代入式（12），更新 Φ，求解 Essd、▽Essd，迭代次數(shù)iter增1，保存當(dāng)前參數(shù)，進(jìn)入下一循環(huán)，直到迭代次數(shù)達(dá)到20或者不滿(mǎn)足循環(huán)條件（即找到最優(yōu)Φ），退出循環(huán)。

其中，Essd、▽Essd、step＿init、step＿best的求解過(guò)程如下：

（1）Essd的求解過(guò)程

Essd的求解過(guò)程見(jiàn)式（11）。

（2）▽Essd的求解過(guò)程

定義 topx、topy、bottomx、bottomy，分別代表當(dāng)前中所有控制點(diǎn)（x，y）向右、上、左、下移動(dòng)當(dāng)前 stepsize后的控制網(wǎng)格坐標(biāo)位置；定義 I＿topx、I＿topy、I＿bottomx、I＿bottomy，分別代表 topx、topy、bottomx、bottomy經(jīng)三次 B樣條彈性配準(zhǔn)后的圖像；定義 E＿topx、E＿topy、E＿bottomx、E＿bottomy，分別代表每個(gè)控制點(diǎn)影響的像素區(qū)域內(nèi) I＿topx、I＿topy、I＿bottomx、I＿bottomy與參考圖像的配準(zhǔn)誤差Essd。

▽Essd的求解過(guò)程見(jiàn)式（13）：

（3）step＿init的求解過(guò)程

初始化Essd為無(wú)窮大，設(shè)置for循環(huán)條件（stepsize＝le－4：le＋3），滿(mǎn)足，則進(jìn)入循環(huán)，將初始 Φ、初始▽Essd以及當(dāng)前stepsize代入式（12），更新Φ，求解Essd，比較當(dāng)前Essd與之前的Essd。若Essd變小，得到較佳初始步長(zhǎng)step＿init，保存當(dāng)前Essd，以備下一次循環(huán)比較；若Essd變大，則直接進(jìn)入下一循環(huán)。當(dāng)不滿(mǎn)足循環(huán)條件，退出循環(huán)，得到最佳初始步長(zhǎng)step＿init，更新當(dāng)前 Essd為初始 Essd。

（4）step＿best的求解過(guò)程

初始化迭代次數(shù)w＝0，找到標(biāo)志found＝0，停止標(biāo)志 stop＝0，step＿best＝0；設(shè)置 while循環(huán)條件（～found），滿(mǎn)足，進(jìn)入循環(huán)，將當(dāng)前 Φ、當(dāng)前▽Essd以及當(dāng)前stepsize代入式（12），更新Φ，求解Essd，比較當(dāng)前Essd與之前的Essd。若Essd變小，得到最佳步長(zhǎng)step＿best（found設(shè) 1），退出循環(huán)；若 Essd變大，步長(zhǎng)減半，進(jìn)入下一循環(huán)再次求解Essd比較。w增1，當(dāng)w大于100，說(shuō)明線(xiàn)性搜索沿著當(dāng)前的搜索方向不能找到一個(gè)可接受的 Φ，stop設(shè)1（尋找 step＿best失?。顺鲅h(huán)。

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

為了驗(yàn)證本研究算法的有效性，選取同一個(gè)人不同時(shí)間的兩幅腦部圖像（一幅作為待配準(zhǔn)圖像，一幅作為參考圖像），大小均為217×181像素，做控制點(diǎn)間距分別為32×32像素、64×64像素、128×128像素（即控制網(wǎng)格分別為9×8、6×5、4×4）的三組對(duì)比實(shí)驗(yàn)。圖5、圖6、圖7分別為其三組實(shí)驗(yàn)的配準(zhǔn)效果，表1為三組實(shí)驗(yàn)的配準(zhǔn)時(shí)間。實(shí)驗(yàn)平臺(tái)為 Intel（R）Core（TM）i5－3210 M CPU，2.50 GHz，4 GB內(nèi)存（2.74 GB可用）。利用 MATLAB R2014a（8.3.0.532）和 Microsoft Visual Studio 2012（C／C＋＋）混合編程完成算法的實(shí)現(xiàn)工作。

從每組實(shí)驗(yàn)的配準(zhǔn)前差分圖像和配準(zhǔn)后差分圖像對(duì)比看出，三組實(shí)驗(yàn)的配準(zhǔn)效果均令人滿(mǎn)意；從三組實(shí)驗(yàn)的配準(zhǔn)時(shí)間對(duì)比看出，控制點(diǎn)間距為64×64像素的配準(zhǔn)結(jié)果最佳，控制點(diǎn)間距為128×128像素和32×32像素次之。

圖5 控制點(diǎn)間距為32×32像素的配準(zhǔn)效果Fig 5 Registration results of 32×32 pixels as controlmesh

圖6 控制點(diǎn)間距為64×64像素的配準(zhǔn)效果Fig 6 Registration results of 64×64 pixels as controlmesh

圖7 控制點(diǎn)間距為128×128像素的配準(zhǔn)效果Fig 7 Registration results of 128×128 pixels as controlmesh

表1 三組控制點(diǎn)間距的配準(zhǔn)時(shí)間Table 1 Registration time of three sets of controlmesh

究其控制點(diǎn)間距為32×32像素和128×128像素的配準(zhǔn)時(shí)間相對(duì)64×64像素較長(zhǎng)的原因，是因?yàn)楫?dāng)控制點(diǎn)間距為64×64像素時(shí)，迭代次數(shù)只達(dá)到12，就得到最優(yōu)控制參數(shù)；而當(dāng)控制點(diǎn)間距為32×32像素或128×128像素時(shí)，迭代次數(shù)均需達(dá)到20才退出循環(huán)，得到最優(yōu)控制參數(shù)。迭代次數(shù)越多，配準(zhǔn)時(shí)間越長(zhǎng)。

實(shí)驗(yàn)結(jié)果分析表明，該算法可得到令人滿(mǎn)意的配準(zhǔn)結(jié)果（綜合配準(zhǔn)效果和配準(zhǔn)時(shí)間），從而驗(yàn)證了本研究算法的有效性。但若選取合適的控制點(diǎn)間距，則可大大縮短配準(zhǔn)時(shí)間，得到更為滿(mǎn)意的配準(zhǔn)結(jié)果。

4 結(jié)論

本研究算法以基于B樣條的FFD模型為彈性形變模型，其中，B樣條網(wǎng)格控制點(diǎn)定義了B樣條自由形變參數(shù)，控制點(diǎn)的移動(dòng)（即參數(shù)的改變）將引起待配準(zhǔn)圖像的變形。以雙線(xiàn)性插值為插值算法，計(jì)算控制點(diǎn)移動(dòng)后位置的灰度，以SSD測(cè)度為配準(zhǔn)依據(jù)，運(yùn)用梯度下降法對(duì)配準(zhǔn)過(guò)程進(jìn)行優(yōu)化，使待配準(zhǔn)圖像逐漸向參考圖像靠攏。實(shí)驗(yàn)證明，本研究算法有效。