歸一化互信息與多分辨率融合的2D-3D配準方法

陳龍， 張峰峰,2， 于凌濤， 孫立寧,2， 干旻峰， 詹蔚

(1.蘇州大學 機電工程學院，江蘇 蘇州 215006；2.蘇州大學 蘇州納米科技協同創新中心，江蘇 蘇州 215123；3.哈爾濱工程大學 機電工程學院，黑龍江 哈爾濱 150001；4.蘇州大學 附屬第一醫院，江蘇 蘇州 215000)

隨著計算機輔助手術技術的成熟發展，借助計算機來幫助醫生進行術前規劃、術中導航，已經逐漸成為一種普通的方式。一般計算機手術導航主要有2類：第一種方法試將術前CT圖像與術中病人暴露的脊骨的三維空間進行3D-3D配準。此方法雖然可以術中較精確的實行導航，但由于需要暴露患者脊柱來放標志點，給患者造成了巨大的傷害，不滿足當前對微創手術的定義。另一種方法主要通過患者術前CT三維圖像與術中X光圖像進行3D-2D圖像的實時配準。雖然精度相對3D-3D配準略低，但可以減小病人的解剖創口，降低術中病人的痛苦。除此以外，X光機成像設備相對較為緊湊，操作比較方便，相比于CT成像設備，它具有更小的輻射，減少了對醫生和病人的傷害[1]。

從配準的實質出發，目前研究中所涉及到的2D-3D配準方法大致有以下3種[2-4]，即基于灰度、梯度和特征的配準方法。基于特征的配準方法一般需要進行邊緣檢測，將邊緣檢測到的輪廓圖像作為配準的框架或基準。但是其需要對影像進行分割，這在某種程度上也展示了配準的準確性往往取決于分割的準確性。與此同時此方法對圖像的要求較高，對邊緣檢測、圖像分割等算法也具有較大的依賴性，所以實施起來有一定的難度，而且配準精度不高。基于梯度的方法一般需要從術前CT圖像圖像提取3D梯度矢量場(表面法線)，需要通過X射線源發出穿過當前3D位置的光線來定義的位置處提取X射線圖像的2D梯度矢量。這種方法非常的復雜，可實施性比較低，而且配準失敗的概率比較大。在目前的研究中發現基于圖像灰度的配準方法相對于上述的兩種方法較為穩定，但是由于使用了大量的灰度信息，導致在配準效率上要略低。

針對上述的情況，本文提出一種基于區域貢獻的歸一化互信息與多分辨率策略相結合的2D-3D配準方法，首先對參考圖像進行輪廓提取，然后對所提取的圖像輪廓進行區域分類，再通過每塊區域的貢獻計算歸一化互信息的值，以減少在使用互信息測度方法時由于灰度級數差距較大導致的誤配準，多分辨率策略在2D-3D配準中的應用有效減少了由于插值與測度方法計算量過大導致的效率較低的問題。

1 基于區域貢獻的歸一化互信息與多分辨率的配準算法

1.1 基于灰度信息的2D-3D配準框架

如圖1所示，一般基于灰度信息的2D-3D配準流程大多數要經過幾何變換，優化算法，圖像插值，和相似性測度這4個部分。正常在做醫學圖像研究分析的時候，大多將患者的幾個圖像放在一起進行對比分析結果。醫學圖像配準主要通過保持一幅圖像固定不變，對另外一幅圖像進行空間變換，使其上標志點與第1幅圖像上的標記點對應重合。這樣可以人體上的解剖點能夠在2幅圖像上對應空間位置一致，能夠幫助醫生在術中實時找到病灶點位置。

圖1 基于灰度信息的圖像配準系統框架Fig.1 Image registration system framework based on gray information

在2D-3D配準過程中，首先保持參考圖像不變，然后開始對生成浮動圖像設置一個大致合理的初始幾何變換參數。在幾何變換的過程中，難免出現變換后的圖像的像素點坐標不是整數的問題。要解決此問題，此時需要引入插值法，通過插值法來找出原來圖像控制點上的像素值；隨后設定一個合適的相似性測度值，對參考圖像和浮動圖像的相似性進行評判，是否到達預先設定的最優值；如果未達到，則重復上述步驟，直至能夠找到一個合適的幾何變換參數來滿足之前設定的精度要求[5-7]。

1.2 基于CT與X光圖的2D-3D配準方法

配準過程中的浮動圖像主要將從虛擬射線源發射出來的X射線穿過CT切片，和每個CT的切片的交點在使用插值法后得到每個交點的CT值。將CT累加值通過映射成灰度圖像就能得到DRR圖像。參考圖像一般選擇脊柱模型的X光圖。為了能使浮動圖像在空間中能更加自由的變換，故此配準框架采用了歐拉幾何變換，能夠使浮動圖像更加方便的實現平移、旋轉。為了能夠最快找到幾何變換參數，本文還引用了梯度下降法，沿梯度下降方向尋求極小值[8]。除上述外，配準框架還包含了基于區域貢獻的歸一化互信息法、光線投射法以及多分辨金子塔配準法。這些方法都能有效提高配準速度。

1.2.1 光線投射插值法

手術過程中獲取的X光圖像為二維圖像，術前CT圖像為三維圖像。因此需要將3D體素投射到二維空間中，采用數字影像重建技術[9-12]方法來進行插值。通過光線投射法可以模擬該方法模擬 X光射線穿過CT體數據的整個過程，同時可以計算出計算在不同介質中X光經衰減、吸收后，最終在投影平面上的像素值。具體的數字影像重建原理如圖2所示。

圖2 數字影像重建原理Fig.2 Principle of digital image reconstruction

一般情況下，當X射線穿過某均勻介質時，其的衰減大致表示為：

I=I0e-μd

(1)

式中：I0為初始時刻輸入射線的強度；μ為不同介質中的衰減系數；d為X射線穿過的距離。人體脊柱受血肉的包圍，因此各部分的衰減系數一般不同。在拍攝X光圖的過程中，式(1)一般可改寫為：

(2)

式中：μi和di表示人體不同組織i的大致線性衰減系數和穿過人體血肉組織i的距離。X射線的衰減模型如圖3所示。

圖3 X射線衰減模型Fig.3 X-ray attenuation model

一般通過CT值可以計算出線性衰減系數，推廣之后得到一般計算公式可以表示為：

(3)

X光射線穿過體數據的衰減可以通過查表獲得，式(3)改寫為：

(4)

式中：F為轉換因子。

1.2.2 基于區域貢獻的歸一化互信息測度

互信息在在早期的信息論中，經常用它來作為一種信息的度量，以表達兩件事情之間的關聯性[13]。目前通常將2幅醫學圖像之間的相似程度通過互信息來衡量。在對互信息進行定義的過程中，必須要對熵的定義有一個清晰的了解。在醫學圖像配準環節中，經常一個系統的不確定性和它的復雜性可以通過熵來描述[14]。假設通過以下定義來定義圖像A的熵：

(5)

此時針對圖像A和B的相關性進行統計，一般與之相對就需要用聯合熵表示，大致定義如下所示：

(6)

圖像A和B的互信息可以由上的定義得到:

I(A,B)=H(B)-H(B|A)=H(A)+H(B)-

(7)

互信息測度方法的缺點是一般情況下在圖像重疊區域較少時，由于圖像采樣點的減少，參與互信息統計的像素點數也隨之減少，因此在圖像灰度值差距較小的情況下，容易造成誤配準。根據聯合熵與個體熵之間的關系，得到歸一化互信息：

(8)

本文提出了一種基于區域貢獻的歸一化互信息測度方法，在進行區域分割后，對圖像中包含各種不同信息的區域進行貢獻評價，以減少背景等無關區域對圖像配準的影響。如圖所示，將配準所用的圖像分割為一系列的區域，對參考圖像進行輪廓提取后，根據各部分區域對互信息測度的貢獻度，將能夠表現出脊柱特征的輪廓所在區域設為一類區域，將脊柱輪廓內部包含的區域設為二類區域，將背景等與配準無關的區域設為三類區域，由此可以實現對配準圖像區域貢獻的劃分。主要劃分區域如圖4所示。

圖4 參考圖像的輪廓提取與區域劃分Fig.4 Outline extraction and regional division of reference images

對三類區域之間的歸一化互信息分別進行求解，根據其貢獻度進行加權，可以得到：

(9)

式中：α+β+θ=1，α>β>θ。由此，將基于區域貢獻的互信息測度方法用于對參考圖像以及經過光線投射插值后的浮動圖像進行相似性測度對比，從而進一步判斷2D-3D配準的結果是否滿足預先設定值的要求。

1.2.3 多分辨率策略

由于在圖像配準框架中采用了光線投射與互信息測度算法，這2種算法具有較高的復雜度，因此可能會導致配準的效率較低。在保證配準精度的前提下，為了提高配準效率，本文在配準環節中引入了多分辨率策略。可以用類似金字塔的形式來闡釋這一算法，通過將一開始具有高分辨率的圖像細分到n層具有不同分辨率的子圖像，按從下往上的順序依次排列。采取從底部到頂端的方式進行配準，將最初的高分辨率圖像放在底部，將第n級低分辨率的圖像放在最頂部[15-16]。

引入多分辨率配準方法需要首先將參考圖像和浮動圖像進行對應n層的劃分。配準過程遵循由上而下，從最頂端最低分辨率的圖像開始，將其配準的結果作為配準初始解。由于一開始配準計算量相對較小，時間較短，這樣就可以依次將上一層的配準結果作為下一層配準的初始解，從而縮短大量的配準時間。整個實驗過程從粗配準到精配準，配準效率得到很大的提高。

圖5 多分辨率策略下的金字塔Fig.5 Multi-resolution pyramid

2 實驗和結果

2.1 實驗平臺硬件組成

整個操作導航的實驗平臺大概通過C形臂X光機，電腦工作站，帶標記的脊柱模型，NDI光學跟蹤器等組成。具體如圖6所示。

圖6 實驗平臺硬件圖Fig.6 Experimental platform hardware diagram

2.2 脊柱正位圖像2D-3D配準

本次2D-3D配準實驗基于ITK與MFC開發包聯合編程的基礎之上完成。一開始需要進行脊柱正位圖像的初始參數選取，一般包括旋轉和平移，選取的初始配準參數:rx=270°，ry=3°，rz=180°，tx=30 mm,ty=60 mm,tz=-120 mm。

為了控制合理的配準時間和配準精度，在配準過程中，將最大的配準步長設為0.1，最小的配準步長設為0.001，通過50次迭代。在CT切片當中選取80張分辨率為768×768的DICOM格式的切片作為浮動圖像。將如圖7所示的X光圖像作為本次實驗的參考圖像，選取干擾性較小的脊柱中間區域進行配準實驗。

圖7 脊柱正位參考圖像Fig.7 Front reference image of the spine

分別采用互信息測度法與基于區域貢獻的歸一化互信息測度法得到的最終配準參數如表1所示，圖8為用配準后參數生成的DRR圖像。

表1 脊柱正位圖像配準后參數Table 1 Front image registration parameters of the spine

圖8 脊柱正位配準后的DRR圖像Fig.8 DRR image of the front registration of the spine

從圖8中可以看出，采用基于區域貢獻的歸一化互信息測度方法進行配準，具有較高的配準精度。為了對其進行定量分析，對其相似性測度值與配準時間進行了對比，通過相似性測度來展現配準的精度。如表2所示，對2種測度方法得到的最終脊柱正位DRR圖像與脊柱正位參考圖像的互信息與相關系數進行了計算，并給出了使用多分辨率策略前后的配準時間。

2.3 脊柱側位圖像2D-3D配準

與此同時為了驗證基于區域貢獻歸一化互信息測度的多分辨率2D-3D配準方法的可行性以及優越性，嘗試進行多角度配準。選擇同段脊柱側位圖像進行配準，選取的同段脊柱側位初始配準參數：選取的初始配準參數:rx=270°，ry=3°，rz=-90°，tx=-200 mm,ty=-140 mm,tz=-120 mm。

表2改進前后算法脊柱正位圖像相似性測度值與配準時間對比

Table2Comparison of image similarity measure and registration time on the front of the spine before and after the improved algorithm

衡量指標配準圖像相似性配準時間/s互信息相關系數未使用多分辨率使用多分辨率互信息測度0.697 80.457 91 281679本文方法0.998 00.689 01 325701

配準中保持迭代次數，最大、最小步長不變。浮動圖像依然選用80張分辨率為768×768的DICOM格式CT切片。采取實驗的步驟，脊柱側位參考圖像為如圖9所示的X光圖像。

圖9 脊柱側位參考圖像Fig.9 Spine lateral reference image

重復實驗步驟，分別采用用互信息測度法與基于區域貢獻的歸一化互信息測度法進行配準實驗，得到的最終配準參數如表3所示，其中圖10為用配準后參數生成的脊柱側位DRR圖像。

表3脊柱側位圖像配準后參數

Table3Parameters after registration of the lateral image of the spine

參數rxryrztxtytz互信息測度271.6-0.6-92.9-201.3-138.8-119.3本文方法270.7-1.2-91.1-199.6-137.2-118.9

從圖10中依然可以看出，采用基于區域貢獻的歸一化互信息測度方法在脊柱側位圖像進行配準，依然適用。對2種測度方法得到的最終脊柱側位DRR圖像與脊柱側位參考圖像的互信息與相關系數進行計算，并給出了使用多分辨率策略前后的配準時間，具體實驗數據如表4所示。

從表2和表4中可以看出，用基于區域貢獻的歸一化互信息測度的2D-3D配準方法得到的配準參數，不管是在脊柱正位還是脊柱側位生成的DRR圖像與參考圖像的相似性測度值都要高于僅用互信息測度的配準方法，其相關系數提高了約58%，其整體配準精度提高了40%左右，證明該相似性測度方法具有較高的配準精度。同時可以適用于多個角度的脊柱圖像的2D-3D配準，具有一定的可操作性以及較廣的適用性。另外，從配準時間可以看出，將多分辨率配準策略應用在2D-3D配準中，可以將配準時間縮短為原來的53%，較為有效地提高了實際配準的效率。

圖10 脊柱側位配準后生成的DRR圖像Fig.10 DRR image generated after lateral registration of the spine

表4改進前后算法脊柱側位圖像配準相似性測度值與配準時間對比

Table4Comparison of image registration similarity measure value and registration time before and after the improved algorithm

衡量指標配準圖像相似性配準時間/s互信息相關系數未使用多分辨率使用多分辨率互信息測度0.685 20.432 91 306682本文方法0.968 30.692 11 340714

2.4 改變步長后的脊柱正位圖像2D-3D圖像配準

在實驗過程中將最大步長和最小步長分別改為0.05和0.000 1，其余參數保持不變，研究此時步長的改變對配準精度以及時間是否有直接的影響。采用實驗1中的脊柱正位初始參數，表5所示為改變配準步長后，分別用互信息測度法與基于區域貢獻的歸一化互信息測度法得到的最終配準參數。圖11為用改變步長后的配準后參數生成的脊柱正位DRR圖像。

表5改變步長后脊柱正位圖像配準后參數

Table5Parameters after registration of the lateral image of the spine

參數rxryrztxtytz互信息測度268.11.5181.3 33.561.8-120.9本文方法269.05.0 180.230.357.9-120.5

圖11 改變步長后的脊柱正位配準后生成的DRR圖像Fig.11 DRR image generated after registration of the frontal image of the spine after changing the step size

從圖11中可以看出，改變配準步長后采用基于區域貢獻的歸一化互信息測度方法在脊柱側位圖像進行配準，配準精度較高。通過相似性測度值與配準時間進行了對比，其所要的配準時間相對于之前的步長更久，精度略有提升，有一定的提升效果。具體實驗數據如表6所示。

表6步長改變后改進前后算法脊柱正位圖像相似性測度值與配準時間對比

Table6Comparison of the image similarity measure and registration time of the frontal image of the algorithm before and after the step change

衡量指標配準圖像相似性配準時間/s互信息相關系數未使用多分辨率使用多分辨率互信息測度0.701 80.469 21 562723本文方法0.999 40.724 01 808786

從表6可以看出，當配準步長變為0.05，最小步長設為0.000 1之后，脊柱正位圖像的相似性測度值在改進算法前后都有了一定程度的提高，即配準的精度相對有了一定的提高。脊柱正位圖像的相似性測度值相對于改變步長之前提升了大概8%，其使用多分辨率算法花費的配準時間相較于改變步長之前的分別增加了6%和12%。從實驗數據可以看出，在步長變小的時候，配準精度相對提高8%，但配準時間增加12%。由此得出雖然精度有所提高，但在手術過程中時間的延遲將會對醫生的手術造成一定的影響，也會直接導致患者的手術風險進一步的提高。所以在滿足精度要求的同時，盡量縮短手術時間，可以極大的提升手術的成功的概率。

3 結論

1)通過實驗可以看出，在不同的步長情況下，同段脊柱的同位置配準時間和精度有所不同，當步長變小時，配準精度相對提高，但取而代之的是配準時間相對變長。需要尋求一個合適的步長，能夠在配準時間和精度之間找到一個平衡。

2)采用基于區域貢獻的歸一化互信息與多分辨率算法后，所需要的配準時間仍然較長，需要進一步縮短配準時間，使醫生能夠提高手術成功率。

3)所采取的實驗對象為脊柱模型，與實際動物或人體實驗仍然存在一定的差距，在人體或動物脊柱外圍存在一定的肌肉組織和皮膚，在實際過程中可能會導致可操作性降低。

4)在手術過程中，把脊柱圖像配準看成類似“剛體”的配準，其實質上脊柱在手術過程中仍然存在一定的變形，解決變形問題，有利于進一步提高配準精度，降低手術風險。