基于點(diǎn)云的膝關(guān)節(jié)脛骨三維CT與MRI圖像配準(zhǔn)

何巍，魏國棟，師為禮，苗語，何飛，楊華民

（長春理工大學(xué)　計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長春　130022）

基于點(diǎn)云的膝關(guān)節(jié)脛骨三維CT與MRI圖像配準(zhǔn)

何巍，魏國棟，師為禮，苗語，何飛，楊華民

（長春理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長春130022）

同時(shí)利用CT和MRI圖像對膝關(guān)節(jié)及脛骨疾病進(jìn)行診斷和治療，是目前臨床上的常規(guī)方法。精確配準(zhǔn)膝關(guān)節(jié)脛骨部分三維CT和MRI圖像可以更有效地利用這兩種圖像中的組織信息。目前，基于圖像特征或者圖像灰度值的方法可以實(shí)現(xiàn)圖像間的配準(zhǔn)，但是都存在一些配準(zhǔn)精度低或者時(shí)間復(fù)雜度高的缺陷。針對膝關(guān)節(jié)脛骨組織的剛性特征，采用一種基于骨組織表面重建的點(diǎn)云配準(zhǔn)方法。首先對數(shù)據(jù)集進(jìn)行骨表面重建，通過重建后的結(jié)果提取點(diǎn)云數(shù)據(jù)，對點(diǎn)云數(shù)據(jù)進(jìn)行下采樣并利用重心進(jìn)行初始變換，最后利用ICP算法對點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn)。對該方法進(jìn)行配準(zhǔn)結(jié)果分析及配準(zhǔn)精度實(shí)驗(yàn)，表明這種方法簡單可靠并且配準(zhǔn)精度較高。

配準(zhǔn)；三維重建；ICP；點(diǎn)云

隨著科學(xué)技術(shù)的不斷創(chuàng)新和發(fā)展，利用不同成像原理的醫(yī)學(xué)成像系統(tǒng)被應(yīng)用于醫(yī)學(xué)輔助診斷和治療當(dāng)中。由于不同的成像方式都有其各自不可替代的特點(diǎn)，以至于在醫(yī)療診斷或治療中通常使用多種成像設(shè)備的影像共同協(xié)作來提供綜合信息。目前，廣泛使用的成像模式有計(jì)算機(jī)斷層掃描成像（CT）、核磁共振成像（MRI）、超聲成像（US）、正電子發(fā)射斷層成像（PET）、單光子發(fā)射斷層成像（SPECT）和數(shù)字剪影血管造影成像（DSA）等［1］。它們可以從各自的角度觀察人體內(nèi)的結(jié)構(gòu)或功能代謝情況。其中CT對密度大、非流動(dòng)性結(jié)構(gòu)表現(xiàn)力強(qiáng)，如：腫瘤、骨等；MRI則對軟組織比較敏感，如：膝關(guān)節(jié)軟骨。為此在膝關(guān)節(jié)部位診斷和治療當(dāng)中，為了更準(zhǔn)確的判斷病情或進(jìn)行治療，需要進(jìn)行CT圖像和MRI圖像的融合，使醫(yī)生可以同時(shí)獲得多方面信息［2］。在對兩種模式的圖像進(jìn)行融合時(shí)，關(guān)鍵的一步就是要將這兩種模式的圖像進(jìn)行精確配準(zhǔn)。

三維圖像配準(zhǔn)就是求使兩個(gè)三維圖像中具有相同特征點(diǎn)的位置在某一坐標(biāo)空間上達(dá)到一致時(shí)的空間變換。目前主要有基于灰度［3］和基于特征［4］的兩種方法。基于特征配準(zhǔn)是通過提取點(diǎn)、線、輪廓、邊緣等圖像特征來計(jì)算變換參數(shù)。這種算法完全依賴于圖像特征的提取精度，所以配準(zhǔn)結(jié)果往往不理想。基于灰度的配準(zhǔn)則是直接利用圖像的灰度信息進(jìn)行圖像間的匹配，而不需要提取任何特征信息。但這種算法需要對三維圖像中所有的點(diǎn)進(jìn)行計(jì)算，所以較基于特征匹配，計(jì)算代價(jià)龐大。

在脛骨CT與MRI圖像進(jìn)行配準(zhǔn)中，骨組織在這兩種圖像上組織區(qū)域明顯且具有無形變的特性。本文就是利用骨組織的這種剛性特征，通過對兩個(gè)圖像中骨組織表面進(jìn)行剛性空間變換，來獲得兩種圖像的空間位置對準(zhǔn)。首先利用面繪制方法對脛骨CT和MRI圖像進(jìn)行三維可視化重建，通過重建后的三維模型提取待配準(zhǔn)點(diǎn)云。之后，通過改進(jìn)的迭代最近點(diǎn)（Iterative Closest Points，ICP）算法對得到的兩個(gè)點(diǎn)云進(jìn)行配準(zhǔn)。最后對配準(zhǔn)過程進(jìn)行分析并對配準(zhǔn)結(jié)果進(jìn)行精度評(píng)價(jià)。

1　點(diǎn)云重建

1.1面繪制

目前醫(yī)學(xué)圖像的存儲(chǔ)格式通常為DICOM （Digital Imaging and Communication in Medicine）格式。CT與MRI數(shù)據(jù)集為一系有序的二維斷層圖像構(gòu)成的三維空間體素（voxel：立體像素）陣列。為了得到兩種圖像中脛骨組織表面的點(diǎn)云進(jìn)行空間位置配準(zhǔn)，首先對CT和MRI數(shù)據(jù)采用面繪制進(jìn)行脛骨三維表面重建，提取重建表面模型的點(diǎn)作為待配準(zhǔn)的點(diǎn)云。

面繪制方法基本思想是通過求解三維數(shù)據(jù)中的等值面，經(jīng)過表面的特征分割處理，再通各個(gè)特征面片的拼接組合，形成三維物體表面。面繪制方法有兩類，一種是基于斷層輪廓線方法，該方法首先構(gòu)建二維圖像輪廓，然后連接相鄰兩二維圖像的輪廓線形成物體表面，這種方法對與具有復(fù)雜輪廓的圖形處理不理想，并且相鄰二維圖像之間的灰度信息沒有得到利用，重建表面往往容易產(chǎn)生跳躍性形變；另一種是基于體素的表面重建方法，基本思想是在體素中提取等值面，按照一定規(guī)則將這些等值面連接形成物體面繪制模型，比如Marching Cube（移動(dòng)立方體）方法［5］就是典型的代表。面繪制方法通常包括三個(gè)過程。分別是斷層圖像預(yù)處理、圖像邊緣分割提取和多邊形擬合，形成物體表面［5］。

1.2Marching Cubes（MC）

Marching Cubes算法由Lorensen等人于1987年提出［6］，該算法是用于提取具有相同灰度特征的等值面（即物體表面）的一種算法。算法的核心思想是利用三維空間圖像中每兩層相鄰的8個(gè)頂點(diǎn)的灰度值與所給的閾值，來計(jì)算這8個(gè)頂點(diǎn)所構(gòu)成的立方體中的等值面位置，然后通過連接這些等值面形成的曲面構(gòu)成物體表面。

MC方法主要分為三步：

（1）根據(jù)閾值提取含有等值面體素，確定分割模式。

圖1為以相鄰兩層斷層圖像k與k+1分別取4個(gè)點(diǎn)構(gòu)成的體素模型，整個(gè)三維數(shù)圖像可看作由多個(gè)體素構(gòu)成。

圖1　體素模型

圖2　Marching cubes的15種拓?fù)浣Y(jié)構(gòu)

對與給定的閾值c，若頂點(diǎn)值大于c，標(biāo)記當(dāng)前頂點(diǎn)為1，則當(dāng)前頂點(diǎn)在形成的等值面的外側(cè)，若頂點(diǎn)灰度值小于c，將該頂點(diǎn)標(biāo)記為0，則該點(diǎn)的位置為等值面的內(nèi)側(cè)。若立方體內(nèi)一條棱的兩頂點(diǎn)分布于等值面的兩邊，則這條棱與等值面有交點(diǎn)。通過8個(gè)頂點(diǎn)不同狀態(tài)可產(chǎn)生256種不同的體素狀態(tài)，根據(jù)旋轉(zhuǎn)對稱性和正反對稱性可精簡得到最終15種拓?fù)浣Y(jié)構(gòu)，如圖2所示［5］。

（2）求等值面與體素邊緣交點(diǎn)。

為方便計(jì)算機(jī)查找，建立256個(gè)全狀態(tài)的查找表，查找索引表共8位，每一位記錄一個(gè)頂點(diǎn)在等值面內(nèi)側(cè)還是外側(cè)的標(biāo)志，用0和1表示，如圖3所示。這樣通過八個(gè)頂點(diǎn)標(biāo)志得到索引值，一個(gè)8位的索引值對應(yīng)已有查找表中的一種狀態(tài)，這樣就得到了交點(diǎn)位置和等值面的構(gòu)成方式。交點(diǎn)位置可以通過對兩個(gè)點(diǎn)的插值獲得，最后結(jié)合交點(diǎn)位置和交點(diǎn)連接方式形成體素內(nèi)等值面。

圖3　狀態(tài)索引表

（3）計(jì)算等值面交點(diǎn)法向量。

等值面交點(diǎn)法向量在三維渲染時(shí)用于計(jì)算光照，取得更好的視覺效果。因?yàn)榈戎得嫔系狞c(diǎn)的在面內(nèi)的梯度分量為零，所以交點(diǎn)處的梯度矢量即為其法向量，可采用中心差分和線性插值的方法來計(jì)算各交點(diǎn)梯度。

1.3生成點(diǎn)云

利用MC算法分別對CT與MRI三維圖像進(jìn)行脛骨和脛骨關(guān)節(jié)的等值面提取，利用每個(gè)體素邊緣與等值面的交點(diǎn)的集合作為待配準(zhǔn)的CT與MRI重建點(diǎn)云。

2　ICP算法

2.1基本原理

ICP（迭代最近點(diǎn)）算法是目前三維點(diǎn)云配準(zhǔn)中使用最多的一種配準(zhǔn)方法。最初是1992年由Besl 和Mckay［7］提出的一種基于特征點(diǎn)的配準(zhǔn)方法。該方法通過迭代優(yōu)化變換參數(shù)矩陣，每次尋找浮動(dòng)點(diǎn)云中每個(gè)點(diǎn)在目標(biāo)點(diǎn)云中的最近點(diǎn)，構(gòu)成對應(yīng)點(diǎn)對，計(jì)算使對應(yīng)點(diǎn)對配準(zhǔn)的變換參數(shù)，若根據(jù)此參數(shù)變換沒有達(dá)到配準(zhǔn)精度則繼續(xù)尋找對近點(diǎn)對進(jìn)行變換參數(shù)求解，直到滿足配準(zhǔn)精度，得到最終的變換參數(shù)。該方法提出較早，但憑借其邏輯簡單、速度快且精度高的特點(diǎn)仍廣泛應(yīng)用在各領(lǐng)域的三維點(diǎn)云配準(zhǔn)當(dāng)中［8，9］。基本原理如圖4所示，在兩個(gè)點(diǎn)云中尋找對應(yīng)點(diǎn)對，計(jì)算變換參數(shù)（圖4（a）），通過變換得到新點(diǎn)云（圖4（b）），如不滿足收斂精度則繼續(xù)迭代直到精確配準(zhǔn)（圖4（c））。

圖4　ICP點(diǎn)云配準(zhǔn)

2.2具體步驟

假設(shè)目標(biāo)點(diǎn)云P，浮動(dòng)點(diǎn)云為Q，Pi、Qi分別為點(diǎn)云P和Q在點(diǎn)i上的坐標(biāo)，迭代k次Q在P中的對應(yīng)點(diǎn)設(shè)為，利用對應(yīng)點(diǎn)對計(jì)算變換矩陣，并通過變換矩陣對浮動(dòng)點(diǎn)云進(jìn)行更新，判斷點(diǎn)云間平均距離是否小于給定值。

對Qk中的每個(gè)點(diǎn)在Pk中搜索相應(yīng)最近點(diǎn)，判斷最近點(diǎn)的依據(jù)是使式（1）中描述的對應(yīng)點(diǎn)均方誤差最小；

采用一種閉合求解方法［9］計(jì)算旋轉(zhuǎn)矩陣R；

旋轉(zhuǎn)矩陣確定后，最優(yōu)平移變換T為經(jīng)過旋轉(zhuǎn)變換后的重心差異；

通過計(jì)算得到的變換得到變換后的新點(diǎn)云Qk′為：

如果dk-dk+1不小于給定的閾值τ（τ＞0），返回第1步繼續(xù)迭代，直到dk-dk+1＜τ或迭代達(dá)到所給最大迭代次數(shù)則計(jì)算結(jié)束［11］。

2.3ICP算法改進(jìn)

傳統(tǒng)的ICP算法可以得到精確的配準(zhǔn)效果，在初值較好的情況下收斂性也較好。但這種方法也存在不少問題。首先，算法要求一個(gè)點(diǎn)集必須包含在另一個(gè)點(diǎn)集當(dāng)中，這一條件在很多情況下無法滿足；其次，要求有相近的初始位置，否則會(huì)陷入局部最優(yōu)得到錯(cuò)誤結(jié)果；而后，算法在尋找最近點(diǎn)對時(shí)計(jì)算量非常大；最后，采用距離最近的點(diǎn)作為對應(yīng)點(diǎn)對有可能產(chǎn)生錯(cuò)誤的對應(yīng)點(diǎn)［11］。對于傳統(tǒng)ICP算法存在的問題，通常采用分階段［12-14］的方法對其中某階段進(jìn)行改進(jìn)，例如可分為控制點(diǎn)的選取、特征度量、空間搜索、點(diǎn)對加權(quán)和變換求解［15］。

配準(zhǔn)過程中我們針對不同的問題采用不同階段的改進(jìn)方法，在保證脛骨關(guān)節(jié)部位配準(zhǔn)精度的同時(shí)減少配準(zhǔn)時(shí)間。由于待配準(zhǔn)的點(diǎn)云位置相差較大，首先利用兩點(diǎn)云重心對點(diǎn)云進(jìn)行初配準(zhǔn)，使兩點(diǎn)云具有相近的初始位置。為了減少計(jì)算量、加快運(yùn)算速度，對點(diǎn)云進(jìn)行一致采樣方法［16］，每隔幾點(diǎn)選取一點(diǎn)作為采樣點(diǎn)云。在特征度量上采用傳統(tǒng)ICP算法的點(diǎn)到點(diǎn)（point-to-point）距離。空間搜索上采用多維二元搜索樹（K-D Tree）［17］加快查找速度。

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

實(shí)驗(yàn)數(shù)據(jù)采用512*512*627的CT數(shù)據(jù)集與512*512*62 MRI的數(shù)據(jù)集。圖5中左側(cè)為CT數(shù)據(jù)脛骨面繪制三維重建渲染結(jié)果，圖5右側(cè)為MRI數(shù)據(jù)脛骨關(guān)節(jié)處面繪制三維重建渲染結(jié)果。

圖5　CT（左）和MRI（右）的三維點(diǎn)云重建渲染結(jié)果

三維重建后得到CT點(diǎn)云中含點(diǎn)210192個(gè)，MRI重建后點(diǎn)的個(gè)數(shù)為106450個(gè)。重建后的點(diǎn)云模型輪廓清晰特征明顯，為之后的配準(zhǔn)提供良好的數(shù)據(jù)基礎(chǔ)。

對于位置相差很大的點(diǎn)云利用兩個(gè)點(diǎn)云的重心進(jìn)行初配準(zhǔn)是必要的，距離較遠(yuǎn)會(huì)使ICP算法陷入局部最優(yōu)配準(zhǔn)當(dāng)中，導(dǎo)致配準(zhǔn)后的結(jié)果錯(cuò)誤，如圖6 （a）為目標(biāo)點(diǎn)云與浮動(dòng)點(diǎn)云的初始位置在z方向相差300mm，此時(shí)，若不對點(diǎn)云進(jìn)行初配準(zhǔn)則配準(zhǔn)的結(jié)果如圖6（b）所示，浮動(dòng)圖像在其初始位置陷入局部最優(yōu)，圖6（c）為經(jīng)過重心初始定位后再進(jìn)行配準(zhǔn)的結(jié)果，浮動(dòng)圖像變換到目標(biāo)圖像的對應(yīng)位置。經(jīng)過初配準(zhǔn)有效減少了局部最優(yōu)的情況，但在某些極端情況下局部最優(yōu)還是會(huì)出現(xiàn)，這時(shí)需要人工介入，確定大致的待配準(zhǔn)區(qū)域，然后再利用ICP算法進(jìn)行配準(zhǔn)。

為了減少計(jì)算量，加快配準(zhǔn)速度，對點(diǎn)云進(jìn)行重采樣。選擇合適的采樣點(diǎn)個(gè)數(shù)對于配準(zhǔn)的精度和時(shí)間影響較大。采樣點(diǎn)個(gè)數(shù)少將缺少原始點(diǎn)云的特征，使得經(jīng)過采樣后的點(diǎn)云與目標(biāo)點(diǎn)云配準(zhǔn)時(shí)產(chǎn)生大量錯(cuò)誤對應(yīng)點(diǎn)對，從而造成配準(zhǔn)得到假的最優(yōu)值，導(dǎo)致配準(zhǔn)失敗。采樣點(diǎn)多會(huì)提高配準(zhǔn)的準(zhǔn)確性，但準(zhǔn)確性的提高代價(jià)是配準(zhǔn)時(shí)間大幅增加。如圖7所示，在相同最大迭代次數(shù)下（1000次），50和100采樣點(diǎn)時(shí)不能正確配準(zhǔn)（紅色），如圖8所示，200個(gè)及以上采樣點(diǎn)時(shí)，隨著采樣點(diǎn)的增加配準(zhǔn)時(shí)間同步增長，而配準(zhǔn)精度不再提高。

圖6　初始相差較大點(diǎn)云配準(zhǔn)（深色為浮動(dòng)圖像）

圖7　采樣點(diǎn)個(gè)數(shù)不同時(shí)配準(zhǔn)結(jié)果

圖8　采樣點(diǎn)個(gè)數(shù)不同時(shí)配準(zhǔn)結(jié)果

對配準(zhǔn)結(jié)果分別采用對應(yīng)點(diǎn)距離、浮動(dòng)點(diǎn)到目標(biāo)點(diǎn)中最近三點(diǎn)構(gòu)成的平面距離和CT與MRI切片圖像中人為標(biāo)記對應(yīng)點(diǎn)距離對配準(zhǔn)精度進(jìn)行評(píng)價(jià)。經(jīng)過多次試驗(yàn)，對應(yīng)點(diǎn)距離平均誤差為1.6mm，點(diǎn)到平面平均誤差為0.9mm，人為標(biāo)記點(diǎn)平均誤差為1.8mm。

4　結(jié)論

為了使用CT圖像與MR圖像融合來提高膝關(guān)節(jié)部分的診斷和治療的精確性。提出了一種基于點(diǎn)云的三維CT圖像與MR圖像膝關(guān)節(jié)脛骨配準(zhǔn)方法。這種方法采用MC算法進(jìn)行點(diǎn)云重建，并對ICP算法進(jìn)行改進(jìn)來進(jìn)行點(diǎn)云配準(zhǔn)。利用多種方法對配準(zhǔn)精度進(jìn)行評(píng)價(jià)表明該方法配準(zhǔn)精度高，計(jì)算量較小，為之后進(jìn)行的融合提供有力的支持，保證膝關(guān)節(jié)部分診斷和治療的精確性。但該方法也存在一個(gè)缺點(diǎn)，在浮動(dòng)圖像表面特征不明顯時(shí)，若浮動(dòng)圖像相對與目標(biāo)圖像為倒置的，可能會(huì)導(dǎo)致配準(zhǔn)后的浮動(dòng)圖像也為倒置的，達(dá)不到配準(zhǔn)的要求，這也是本方法今后的改進(jìn)方向。

［1］Petra A，den Elsen，et al.Medical image matching-a review with classification［J］.IEEE Engineering in Medical and Biology，1993（3）：26-39.

［2］ 楊鎮(zhèn)源.在體三維光聲圖像的腫瘤微血管信息提取方法研究［D］.長春：長春理工大學(xué)，2014.

［3］ Maes F，Collignon A，Vandermeulen D，et al.Multimodality image registration by maximization of mutual information［J］.IEEE Transactions on Medical Imaging，1997（16）：187-198.

［4］Staring M，Vander-Heide U A，Klein S.Registration of Cervical MRI using multi-feature mutual information［J］.IEEE Transactions on Medical Imaging，2009，28（9）：1414-1421.

［5］ 葉青.三維重建技術(shù)在醫(yī)學(xué)圖像中的研究與應(yīng)用［D］.西安：西安電子科技大學(xué)，2009.

［6］Lorensen WE，Cline HE.Marching cubes：a high resolution3Dsurfaceconstructionalgorithm［J］. Computer Graphics，1987，21（4）：163-169.

［7］Besl PJ，Mckay ND.Amethod for registration of 3-d shapes［J］.IEEE Transactionson Pattern Analysisand Machine Intelligence，1992，14（2）：239-256.

［8］ Furong Peng.Street view cross-sourced point cloud matching and registration［J］.Image Processing（ICIP），IEEE International Conference on，2014，2026-2030.

［9］NieBner M，Dai A，F(xiàn)isher M.Combining inertial navigation and ICP for real-time 3D surface reconstruction［J］.Eurographics 2014-Short Papers，2014：13-16.

［10］Berthold K P Horn.Closed-form solution of absolute orientation using unit quaternions［J］.Journal of the Optical Society of America A.1987，4（4）：629-642.

［11］ 劉承香，阮雙琛，劉繁明，等.基于迭代最近點(diǎn)算法的地形匹配算法可靠性分析［J］.深圳大學(xué)學(xué)報(bào)理工版，2015（01）：22-26.

［12］Rusinkiewicz S，Levoy M.Efficient Variants of the ICP Algorithm［C］.The Third International Conference on 3D Digital Image and Modeling.Quebee city，Canada，2001：145-152.

［13］ 伍毅.三維掃描信息獲取的深度圖像配準(zhǔn)算法設(shè)計(jì)及開發(fā)［D］.杭州：浙江大學(xué)，2005.

［14］Nishino K，Ikeuchi K.Robust simultaneous registrationofmultiplerangeimagescomprisinga large number of points［J］.Electroics and Communications in Japan（PartⅡ：Electronics），2004，87（8）：61-74.

［15］ 李世飛，王平，沈振康.迭代最近點(diǎn)算法研究進(jìn)展［J］.信號(hào)處理，2009，25（10）：1582-1588.

［16］Greg T，Marc L.Zippered polygon meshes form range image［C］.Proceedings of the13th InternationalConferenceonPatternRecognition1996：879-883.

［17］Zhang Z.Iterative point matching for registration of freeform curves and surfaces［J］.International Journal of Computer Vision.1994，13（2）：119-15.

Point Cloud Based Registration of 3D CT and MRI Tibia Image

HE Wei，WEI Guodong，SHI Weili，MIAO Yu，HE Fei，YANG Huamin
（School of Computer Science and Technology，Changchun University of Science and Technology，Changchun 130022）

In order to improve the accuracy of diagnosis or operation on knee and tibia，the registration between CT and MRI 3D images is required in clinical application.Accurate registration can help doctors understand anatomical structures correctly.In order to improve the algorithm for registering CT and MRI images of tibia，this paper proposes a registration method based on surface point cloud of bone.First，surface reconstruction of bone is performed by using the Marching Cubes method，and point clouds are extracted as the reconstructed result.Second，a down sample method is utilized for the point clouds to reduce computational consumption，and the transformation of displacement is estimated by using center of gravity.Finally，ICP method is used for point clouds registration.The experimental results show that the proposed method has high registration accuracy and low computational consumption.

registration；3D reconstruction；ICP；point cloud

TP317.4

A

1672-9870（2015）05-0131-05

2015-08-08

何巍（1978-），女，博士，講師，E-mail：hw@cust.edu.cn

楊華民（1963-），男，教授，博士生導(dǎo)師，E-mail：yhm@cust.edu.cn