摘 要：針對大尺度形變醫學圖像配準速度慢和精度低的特點，提出一種結合薄板樣條(TPS)和B 樣條的彈性配準方法。該方法采用尺度不變特征變換算法(SIFT)進行圖像特征提取與匹配，利用TPS 算法將特征點對作為輸入進行預處理，以降低浮動圖像的形變尺度，從而提高下一步B 樣條配準的速度與精度。然后使用局部區域細化層次B 樣條方法將TPS生成的較稀疏的形變網格作為初始網格，結合有限記憶優化算法(L-BFGS)對控制網格做進一步地處理，此過程只對形變較大的局部區域進行細化，以實現與參考圖像的快速精確配準。實驗結果表明，該方法較層次B 樣條方法有效地提高了配準的速度和精度。

關鍵詞：醫學圖像；彈性配準；薄板樣條；層次B 樣條；大尺度形變

引言：醫學圖像配準是醫學圖像融合、圖像與標準圖譜的匹配、外科手術規劃與設計等研究的基礎，對提高臨床診斷治療、病情監測、外科手術水平等有著非常重要的積極作用[1,2,3]。根據不同的臨床需求，在配準時可采用不同的算法。在某些情況下，使用快速的剛性配準算法就能夠達到要求，但是在浮動圖像與參考圖像之間存在大尺度形變的情況下，就需要通過非線性變換模型來描述復雜的形變過程從而解決問題，如形變復雜多樣的軟組織圖像的配準、不同模態圖像之間的配準、不同個體圖像之間的配準等。針對這類問題，可將目前的方法大體分為基于徑向基函數和基于偽物理模型兩大類。前者采用徑向基函數來表示形變域，如多項式函數[4]、薄板樣條(Thin Plate Spline，TPS)[5,6]、B 樣條[7,8,9]、小波基函數等；后者使用偏微分方程來描述形變的平衡狀態，如彈性模型、粘性流體模型、光流場模型等[10]。

層次B 樣條[11,12,13]是目前常采用的彈性配準方法，它綜合考慮了密集和稀疏的控制網格特點，力求在保證配準精度的同時，提高配準的速度。但其不足在于對每一層都要進行全局細化，不僅造成了大量不必要的計算，而且對于在上一層已經配準好的區域也將會受到整個形變域的影響。TPS 算法[5]由于具有全局變形且沒有冗長迭代的特點，同樣是一種建模變形的有效工具。然而，TPS 雖然可以實現快速計算，但是與層次B 樣條相比會產生較粗糙的配準結果。

針對層次B 樣條配準算法很難對存在大尺度形變的浮動圖像快速地達到理想配準精度的不足，本文結合TPS 的快速和B樣條的精度，采用TPS 算法做預處理，然后使用局部區域細化層次B 樣條方法做進一步精確處理并與參考圖像配準。實驗結果表明，在配準的時間和質量上，該方法較層次B 樣條配準算法都得到了提高。1 算法框架本文將TPS 算法與局部區域細化層次B 樣條方法相結合，進行大尺度醫學圖像彈性配準。算法框架如圖1 所示。用函數) (XR 和 ) (Y F 分別表示兩幅輸入圖像，其中R為參考圖像，F為浮動圖像，X 和Y 分別代表各自圖像的定義域。) (t S 為相似度函數，其中t 為控制參數，用均方差MSD(Mean SquareDifference，MSD)作為相似性測度來衡量兩幅圖像的匹配效果。

此時，配準就是尋找一個形變函數，以使) (t S 取得最小值。

RR((XX)) MMSSDDBB樣樣條條變變換換g(X)S(t)S*F(g(X))SIFT 算法進行特征提取與匹配SIFT 算法進行特征提取與匹配特征點對RR((XX))FF((YY))初步形變網格TTPPSS變變換換f(X)t計算所有控制點領域的MSD計算所有控制點領域的MSD插插值值提高局部區域的控制網格分辨率提高局部區域的插值 控制網格分辨率 插值LL--BBFFGGSS圖1 算法框架首先，輸入待配準圖像 ) (X R 和 ) (YF ，使用尺度不變特征變換算法(Scale Invariant Feature Transform，SIFT)[14,15]進行圖像特征提取與匹配；將上一步得到的特征點對坐標作為TPS 算法的輸入，對 X 進行 TPS 變換，從而得到新的區域 ) (Xf 坐標，通過插值生成初步形變網格；將該形變網格作為局部區域細化層次B 樣條的初始網格，對X 進行B 樣條變換，從而得到區域 ) (Xg 坐標，通過插值得到浮動圖像在區域 ) (Xg 的取值)) ( ( X g F ；利用 MSD 計算參考圖像 ) (XR 和插值圖像 )) ( ( X g F間的相似度 ) (t S ；將 ) (t S 輸入到優化模塊中，利用有限記憶優化算法(Limited memory Broyden-Fletcher-Goldfarb-Shanno ，L-BFGS)[16]進行計算得到最優變換參數t ，此過程需要通過迭代來實現，直到在該層使 ) (t S 取得最小值；然后，對形變較大的局部區域，提高控制網格分辨率進行再一次配準，直到滿足配準要求；最后，輸出配準結果圖像* S 。

2 TPS 算法降低形變尺度為了能夠同時提高配準的速度和精度，本文采用TPS 算法做預處理，而SIFT 算法提取圖像特征點并匹配是進行TPS 預處理的前提。一般情況下，SIFT 算法會產生大量的、健壯的匹配點對，根據經驗，總共選用200 對左右的匹配點對，就能夠使TPS 算法獲得較好的初步形變網格。TPS 算法是將插值問題模擬為一個薄金屬板在點約束下的彎曲變形。它首先確保配準點集的一一對應，然后利用最小化TPS 彎曲能來聯合求解該點集之間的匹配矩陣及映射參數，使圖像能夠在控制點上達到精確吻合，并且在其約束下通過插值以使其他點也盡可能獲得較好的糾正。表示浮動圖像變形能量的代數式可定義為：

2 2 z(x, y)  U(r)  r logr (1)它的解為 U Ux y2 2 ( ) 2222    ，其中 22 y x r  ， )(rU 是構建薄板樣條的基函數。通過SIFT 算法得到的浮動圖像上的特征點坐標記為(x, y) 。要使金屬板在點(xi , yi )處的高度為 iz ，并且讓其具有最小彎曲能量，即在二維空間中找到一個形變函數 ) , ( y x f ，使薄板樣條插值的罰函數 ) ( f Es 最小。該罰函數可定義為：

E f dxdyyfx yfxfs   ( )  (( )  2( )  ( ) ) 2 2 222 222(2)令 ) , ( i i i y x P ， ni , , 2, 1 ，設 ) , ( y x Pr i   ，為控制點對之間的歐幾里德距離。定義3  n 矩陣1 12 2111 n nx yx yx y       P (3)及n n  矩陣12 121 21 20 ( ) ( )( ) 0 ( )( ) ( ) 0nnn nU r U rU r U rU r U r       K (4)從而得到矩陣   TK PL =P 0(5)其中，0 為3 行3 列的零矩陣， T P 代表P 的轉置。構建行向量1 1 1 2 2 2 ( ( , ) , ( , ) , , ( , ) ) n n n V  z x y z x y ，z 列x 向y 量( , 0 , 0 , 0 ) T Y  V 。通過下式：

11 ( | )Tx y a a a  L Y  W (6)可以得到1 2 ( , , , )Tn W  w w w 和系數 1a , x a , y a 。進而構造函數f (x, y) ：

    nif x y a ax x ay y wiU Pi x y1( , ) 1 ( ( , ) ) (7)該函數對于所有i 有f (xi , yi )  zi ，并使罰函數 Es ( f ) 最小。通過 f (x, y) 得到的形變網格，將作為局部區域細化層次 B樣條的初始網格。

3 B 樣條算法細化形變網格并配準3.1 B 樣條初步細化采用B樣條自由形變模型[12]將TPS 算法生成的較稀疏的形變網格覆蓋在浮動圖像上，通過調整控制點的位置，帶動所嵌套的浮動圖像發生形變。B 樣條具有良好的局部控制特性，浮動圖像上某點處的位移值僅與其相鄰的2 4 個控制點有關。因此，圖像中部分控制點移動僅會影響其周圍有限范圍內的空間變化，而不會影響整幅圖像的像素點位置。__由于三次B 樣條具有二階連續性和計算高效性等特性，因此本文采用三次均勻B 樣條基函數。為了能夠根據TPS 算法得到的網格控制點的位移值求得任意位置處的值，在二維空間中構造一個形變函數 ) , ( y x g ，將與控制點相關的信息盡可能光滑地傳播于浮動圖像中的所有點。該形變函數的表達式為：

m p l q ml mg x y Bl u B v    ,3030( , ) ( ) ( ) (8)其中1  xp ， 1  yq ，   x xu ，   y yv ，表示取整運算。 ) 3, , 0( l Bl 代表三次均勻B 樣條的l 次基函數，其形式為： ( ) (1 ) / 6 3B0 u   u ， ( ) (3 6 4) / 6 3 2B1 u  u  u  ，( ) ( 3 3 3 1) / 6 3 2B2 u   u  u  u  ， ( ) / 6 33 u u B  ，其中 10u 。

它們代表控制點對 ) , ( y x g 的權重值，即根據控制點到 ) , ( y x 的距離來加權各個控制點對 ) , ( y x g 的貢獻。

形變函數 ) , ( y x g 生成的曲面分別由 x 方向上的 3r 個節點和y 方向上的3  s 個節點同時控制。記TPS 生成的初步形變網格為 ，每個控制點間的距離為j i,  ，它控制著形變程度。

將浮動圖像中點的運動位移分解為x 和y 兩個方向，并分別求出所有點在這兩個方向上的運動位移。各控制點在x 和y 方向上的位移值即為待優化的參數，網格中的每個節點有2 個參數表示這個節點的位移，因此參數總數為2(r  3)(s  3) ，通過最優化過程求得。網格中的每一個節點都需要迭代計算，迭代時間取決于L-BFGS 優化算法和公差要求。通常要找到最佳的變形參數設置需要較長的計算時間。與TPS 相比它雖然花費了更長的時間，但得到了一個更精確的結果。

3.2 局部區域精確細化當浮動圖像的局部存在較大形變時，如果直接采用較小的控制網格進行配準，很可能會導致發生局部扭曲的現象。其原因在于需要優化的參數多且變動范圍大，尋優結果很容易陷入局部極值。層次B 樣條采用的是由粗到細的控制網格，每一層的控制點間距由上一層控制點的間距二分而得，并利用上一層的配準結果進行更精確的形變配準，可以較好地解決這類問題。

但是，層次B 樣條的每一層控制網格都是針對整幅浮動圖像的，即每層都要進行全局細化。當浮動圖像在初步配準后，大部分區域已經達到了較好的配準結果，只有局部區域還存在一些差異時，層次B 樣條方法不僅會造成不必要的計算，且對于已經配準好的區域也將會受到整個形變域的影響。因此，本文提出局部區域細化層次B 樣條方法來解決上述問題。

在使用TPS 生成的較稀疏的控制網格配準后，對兩幅圖像中所有控制點所對應的像素點的8 領域計算MSD，對大于給定閾值的控制點領域，鎖定該控制點所在的（4  3）（4  3）控制網格，并對此區域內其他點不再進行計算。當與該區域相鄰的控制網格也含有領域的MSD 大于所給閾值的控制點時，將控制網格合并。然后，提高這些需要進行進一步細化的局部區域的控制網格分辨率，以完成精確配準，如圖2 所示。在所有局部區域中使用同樣的方法，直到所有的控制點的領域都滿足要求時，結束配準。

.圖2 局部區域精確細化4 算法流程本文針對大尺度醫學圖像彈性配準問題，提出了以TPS 算法做預處理，降低浮動圖像的形變尺度，然后使用局部區域細化層次B 樣條方法進行精確配準的算法，解決了大尺度形變圖像配準速度慢且精度低的問題。算法流程如圖3 所示：

開開始始SSIIFFTT算算法法進進行行圖圖像像特特征征提提取取與與匹匹配配結結果果是是否否滿滿足足要要求求結結束束YNTTPPSS算算法法得得到到初初步步形形變變網網格格BB樣樣條條算算法法對對形形變變網網格格做做進進一一步步形形變變計計算算所所有有控控制制點點領領域域的的MMSSDD誤誤差差是是否否大大于于給給定定閾閾值值YN得得到到配配準準結結果果圖圖像像提高局部區域控制網格分辨率提高局部區域控制網格分辨率圖3 算法流程圖配準的基本步驟如下：

1) 用SIFT 算法提取R 、F 的特征點并匹配，得到特征點對坐標；2) 將1)中得到的特征點對坐標作為TPS 算法的輸入，根據式（7）求得形變函數 f (x, y) ，生成初步形變網格 ；3) 將2)中得到的形變網格 作為局部區域細化層次B 樣條的初始網格，根據式（8）結合L-BFGS 優化算法計算變換參數；4) 使用相似性測度函數判斷參數是否達到最優，若是，則得到在較稀疏的控制網格下的最優解，否則，返回3)，計算新的變換參數，直到在該層使MSD 取得最小值；5) 計算兩幅圖像中所有控制點的領域的MSD，對于大于給定閾值的領域，提高該控制點所在的（4  3）（4  3）

控制網格的分辨率，返回3)，計算新的變換參數，直到在該層使MSD 取得最小值。在所有局部區域中使用同樣的方法直到所有的控制點領域都滿足要求時，結束配準。

5 實驗結果及分析5.1 圖像配準結果在本實驗中，使用本文提出的配準方法對大量的醫學圖像進行了仿真實驗驗證，圖像大小均為512 512 ，選取控制點領域的MSD 閾值為80。實驗平臺為Intel(R) Core(TM) i3-2120CPU(雙核四線程，3.30GHz)，4G 內存。在 Matlab V8.0 環境下編寫應用程序。

圖4 和圖5 分別顯示了兩組肺部CT 圖像使用層次B 樣條方法和本文方法的初始網格及最后配準結果。層次B 樣條方法是以較粗的規格化控制網格作為初始網格，而本文方法以TPS做預處理得到的形變網格作為局部區域細化層次B 樣條的初始網格。

(a) 層次B 樣條的初始網格 (b) 本文方法的初始網格(c) 參考圖像 (d) 浮動圖像(e) 層次B 樣條配準結果 (f) 本文方法配準結果圖4 實驗圖像與配準結果1(a) 層次B 樣條的初始網格 (b) 本文方法的初始網格(c) 參考圖像 (d) 浮動圖像(e) 層次B 樣條配準結果 (f) 本文方法配準結果圖5 實驗圖像與配準結果25.2 圖像配準效果的評價5.2.1 評價方法本文采用MSD 和互信息（mutual information，MI）來評價兩幅圖像的匹配效果。設R和F 在點(x, y)處的對應像素灰度分別為 ) , ( y x R 和 ) , ( y x F ，MSD 的計算公式為：

( , )2 ( ( , ) ( , ))x yR x y F x yMSD (9)其中，  代表整個圖像區域。MSD 的大小反映了R 和F 對應像素差異的大小，其值越小，兩幅圖像差異越小，配準結果越精確。

另一評價指標為互信息，圖像R 和F 之間的互信息可表示為：

MI (R, F)  H(R)  H(F)  H(R, F) (10)其中 ) (R H 和 ) (FH 分別代表 R和 F 的熵， ) , ( F R H 代表 R和F 的聯合熵。當R 和F 實現最優配準時，互信息值最大。

5.2.2 配準效果評價表1 給出了層次B 樣條方法和本文方法配準結果的MSD和MI 等定量評價指標的比較，可以看出，本文方法較層次B樣條方法能夠得到更精確的配準結果。

表1 配準結果對比項 MSD 互信息 運行時間/s配準前 840.523 1.1797層次B 樣條 57.486 1.5137 81.21本文方法 10.179 1.5341 63.345.3 圖像配準效率的評價圖6 顯示了使用層次B 樣條方法和本文方法配準時在計算時間上的比較。沿x 軸，公差逐漸增大，配準精度隨之降低。

可以看出，當要求較高精度時，使用層次B 樣條配準需要更長的時間，而本文提出的方法采用TPS 算法做預處理且只對局部區域進行細化，可以在花費更少計算時間的同時達到相同的精度。然而，當要求較低精度時，這兩種方法之間的差異就不明顯了。

0204060801001201.00E-04 1.50E-04 2.00E-04 2.50E-04 3.00E-04公差計算時間（秒）

層次B樣條方法本文方法圖6 層次B 樣條方法和本文方法計算時間的比較表2 給出了層次B 樣條方法和本文方法使用相同的相似性測度且公差同為1.00E-4 時，多組圖像的計算時間的比較。從實驗結果可以看出，本文方法較層次B 樣條方法具有更高的配準效率。

表2 計算時間（秒）

CT 圖像 層次B 樣條 本文方法實例1 81.21 63.34實例2 96.49 71.19實例3 75.43 56.57實例4 57.17 45.36實例5 73.39 57.426 結束語本文提出了一種結合TPS 算法和局部區域細化層次B 樣條方法的彈性配準方法，該方法對大尺度形變醫學圖像具有較好的配準結果。當花費相同的時間時，它明顯提高了配準的精度；當精度相同時，它較層次B 樣條配準方法花費更少的時間。通過對大量醫學圖像進行配準實驗及分析，結果表明，本文提出的方法較層次B 樣條配準方法在配準效率和精度上都有一定程度的提高。

參考文獻：

[1] Sotiras A, Davatzikos C, Paragios N. Deformable medical imageregistration: a survey[J]. IEEE Trans on Medical Imaging, 2013, 32(7):1153-1190.

[2] Markelj P, Tomaževič D, Likar B, et al. A review of 3D/2D registrationmethods for image-guided interventions[J]. Medical Image Analysis, 2012,16(3): 642-661.

[3] 趙增軍, 李寶生. 基于模型的醫學圖像形變配準在圖像引導放療中的應用研究進展[J]. 中國圖象圖形學報, 2009, 14(8): 1504-1509.

[4] Pluim J P W, Maintz J B A, Viergever M A. Mutual-information-basedregistration of medical images: a survey[J]. IEEE Trans on MedicalImaging, 2003, 22(8): 986-1004.

[5] Bookstein F L. Principal warps: thin-plate splines and the decomposition ofdeformations[J]. IEEE Trans on Pattern Analysis and MachineIntelligence, 1989, 11(6): 567-585.

[6] Rohr K, Stiehl H S, Sprengel R, et al. Landmark-based elastic registrationusing approximating thin-plate splines[J]. IEEE Trans on MedicalImaging, 2001, 20(6): 526-534.

[7] Chun S Y, Fessler J A. A simple regularizer for B-spline nonrigid imageregistration that encourages local invertibility[J]. IEEE Journal ofSelected Topics in Signal Processing, 2009, 3(1): 159-169.

[8] Jacobson T J, Murphy M J. Optimized knot placement for B-splines indeformable image registration[J]. Joint AAPMCOMP Meeting Program,2011, 38(6): 334-338.

[9] 龐皓文, 孫小楊, 楊波, 等. 基于B 樣條的彈性配準在放療CT 圖像間匹配的應用[J]. 中國醫學物理學雜志, 2011, 28(5): 2887-2889.

[10] 申艷平. 醫學圖像配準技術[J]. 中國醫學物理學雜志, 2013, 30(1):3885-3889.

[11] Zhang Hongying, Zhang Jiawan, Sun Jizhou, et al. Non-rigid imageregistration algorithm based on B-splines approximation[J]. Transactionof Tianjin University, 2007, 13(6): 447-451.

[12] Zhuang Xiahai, Arridge S, Hawkes D J, et al. A nonrigid registrationframework using spatially encoded mutual information and free-formdeformations[J]. IEEE Trans on Medical Imaging, 2011, 30(10):1819-1828.

[13] 李彬, 歐陜興, 田聯房, 等. 基于自適應自由變形法和梯度下降法的胸部多模醫學圖像配準[J]. 計算機應用研究, 2009, 26(10): 3978-3982.

[14] Huo Ju, Yang Ning, Cao Maoyong, et al. A reliable algorithm for imagematching based on SIFT[J]. Journal of Harbin Institute of Technology,2012, 19(4): 90-95.

[15] Goncalves H, Corte-Real L, Goncalves J A. Automatic image registrationthrough image segmentation and SIFT[J]. IEEE Trans on Geoscience andRemote Sensing, 2011, 49(7): 2589-2600.

[16] Morales J L. A numerical study of limited memory BFGS method[J].Journal of Applied Mathematics, 2002, 15(4): 481-487.