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

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

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

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

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

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

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

配準；三維重建；ICP；點云

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

三維圖像配準就是求使兩個三維圖像中具有相同特征點的位置在某一坐標空間上達到一致時的空間變換。目前主要有基于灰度［3］和基于特征［4］的兩種方法?；谔卣髋錅适峭ㄟ^提取點、線、輪廓、邊緣等圖像特征來計算變換參數(shù)。這種算法完全依賴于圖像特征的提取精度，所以配準結(jié)果往往不理想?；诨叶鹊呐錅蕜t是直接利用圖像的灰度信息進行圖像間的匹配，而不需要提取任何特征信息。但這種算法需要對三維圖像中所有的點進行計算，所以較基于特征匹配，計算代價龐大。

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

1　點云重建

1.1面繪制

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

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

1.2Marching Cubes（MC）

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

MC方法主要分為三步：

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

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

圖1　體素模型

圖2　Marching cubes的15種拓撲結(jié)構(gòu)

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

（2）求等值面與體素邊緣交點。

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

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

（3）計算等值面交點法向量。

等值面交點法向量在三維渲染時用于計算光照，取得更好的視覺效果。因為等值面上的點的在面內(nèi)的梯度分量為零，所以交點處的梯度矢量即為其法向量，可采用中心差分和線性插值的方法來計算各交點梯度。

1.3生成點云

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

2　ICP算法

2.1基本原理

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

圖4　ICP點云配準

2.2具體步驟

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

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

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

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

通過計算得到的變換得到變換后的新點云Qk′為：

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

2.3ICP算法改進

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

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

3　實驗結(jié)果及分析

實驗數(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（右）的三維點云重建渲染結(jié)果

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

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

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

圖6　初始相差較大點云配準（深色為浮動圖像）

圖7　采樣點個數(shù)不同時配準結(jié)果

圖8　采樣點個數(shù)不同時配準結(jié)果

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

4　結(jié)論

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

［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］ 劉承香，阮雙琛，劉繁明，等.基于迭代最近點算法的地形匹配算法可靠性分析［J］.深圳大學(xué)學(xué)報理工版，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］ 伍毅.三維掃描信息獲取的深度圖像配準算法設(shè)計及開發(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］ 李世飛，王平，沈振康.迭代最近點算法研究進展［J］.信號處理，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