基于B樣條的醫(yī)學圖像彈性配準*

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

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

1 引 言

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

本研究采用后者。該方法通過操縱套用在圖像區(qū)域上的控制點網格來模擬彈性形變，然后通過優(yōu)化算法優(yōu)化網格控制點坐標位置的方法來完成彈性配準。

2 配準流程

配準流程圖見圖1。

圖1 配準流程圖Fig 1 Registration flow diagram

2.1 待配準圖像網格化

2.1.1 建立控制網格 在待配準圖像上建立控制網格，見圖2（以11×11控制網格為例），相對控制網格，定義像素點坐標位置為（i，j），則像素空間（m×n控制網格）表示為：

控制空間（（m＋3）×（n＋3）控制網格）表示為：

圖211 ×11控制網格Fig 2 Congtrolmesh of 11×11

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

2.1.2 建立像素網格 設圖像在X和Y方向的像素點數量分別為δxm和δyn，其中，δx和δy分別表示單位控制網格在X和Y方向的間距，相對像素網格，像素點坐標位置為（δxi，δyj），則：

像素空間表示為：

控制空間表示為：

基于B樣條的空間變換需要用到控制點相對像素網格的坐標位置作為初始參數。

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

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

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

對于任意一像素點（x，y），經三次B樣條彈性變形后的坐標位置（相對像素網格）可表示為：

其中：φi＋l，j＋m表示最近鄰 4×4個控制點的坐標位置，i、j分別表示最近鄰4×4個控制點的單元索引號，i＝［x／δx］－1，j＝［y／δy］－1；［］表示取下整數；u、v分別表示x、y的相對單元控制網格的位置，u＝x／δx－［x／δx］，v＝y(tǒng)／δy－［y／δy］；

B0（u）～B3（u）表示三次B樣條基函數，其表達式為：

其中，0≤u＜1，這些函數起到權函數的作用，它們根據控制點到（x，y）的距離大小來加權每個控制點對 T（x，y）的貢獻［7］。

其模型構造見圖2（見圖中區(qū)域4）。

2.3 雙線性插值

雙線性插值是將待配準圖像的像素點（x，y）經B樣條變換后的坐標位置（Tx，Ty）作為插值點，利用該插值點最近鄰4個像素點的灰度值，計算該坐標位置的灰度值，并賦給配準后圖像的對應位置（x，y），見圖 3。具體插值過程［8］如下：

圖3 雙線性插值示意圖Fig 3 The sketch map of bilinear interpolation

圖4 雙線性插值原理圖Fig 4 The illustrative diagram of bilinear interpolation

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

則插值點（Tx，Ty）的灰度值 g（Tx，Ty）為：

2.4 SSD測度

SSD測度公式［9］見式（11）：

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

2.5 梯度下降法［10］

根據式（6），求解最優(yōu)Φ為最主要的步驟，其主要思想是尋求使相似性測度Essd最小的空間變換位置。求解最優(yōu)Φ是通過先求解Essd的梯度▽Essd，然后通過梯度下降法對Φ更新來實現的，Φ的更新公式見式（12）：

其中，stepsize為步長，▽Essd是對控制點的每個網格位置處的變形的梯度。這里步長也要在速度和精度之間取得折中，通過變步長法實現由粗到細的快速精確配準，變形模型越精細，步長越小。

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

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

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

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

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

（1）Essd的求解過程

Essd的求解過程見式（11）。

（2）▽Essd的求解過程

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

▽Essd的求解過程見式（13）：

（3）step＿init的求解過程

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

（4）step＿best的求解過程

初始化迭代次數w＝0，找到標志found＝0，停止標志 stop＝0，step＿best＝0；設置 while循環(huán)條件（～found），滿足，進入循環(huán)，將當前 Φ、當前▽Essd以及當前stepsize代入式（12），更新Φ，求解Essd，比較當前Essd與之前的Essd。若Essd變小，得到最佳步長step＿best（found設 1），退出循環(huán)；若 Essd變大，步長減半，進入下一循環(huán)再次求解Essd比較。w增1，當w大于100，說明線性搜索沿著當前的搜索方向不能找到一個可接受的 Φ，stop設1（尋找 step＿best失敗），退出循環(huán)。

3 實驗結果與分析

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

從每組實驗的配準前差分圖像和配準后差分圖像對比看出，三組實驗的配準效果均令人滿意；從三組實驗的配準時間對比看出，控制點間距為64×64像素的配準結果最佳，控制點間距為128×128像素和32×32像素次之。

圖5 控制點間距為32×32像素的配準效果Fig 5 Registration results of 32×32 pixels as controlmesh

圖6 控制點間距為64×64像素的配準效果Fig 6 Registration results of 64×64 pixels as controlmesh

圖7 控制點間距為128×128像素的配準效果Fig 7 Registration results of 128×128 pixels as controlmesh

表1 三組控制點間距的配準時間Table 1 Registration time of three sets of controlmesh

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

實驗結果分析表明，該算法可得到令人滿意的配準結果（綜合配準效果和配準時間），從而驗證了本研究算法的有效性。但若選取合適的控制點間距，則可大大縮短配準時間，得到更為滿意的配準結果。

4 結論

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