白蘭蘭，鄧雙城，朱潤澤，鄧寧升

(北京石油化工學(xué)院，北京 102617)

經(jīng)顱磁刺激是一種能夠無創(chuàng)傷地在大腦中產(chǎn)生局部刺激的方法，通過刺激大腦的相應(yīng)的功能區(qū)達(dá)到治療效果，可用于治療抑郁癥、狂躁癥、帕金森病、癲癇、中風(fēng)等精神及神經(jīng)科疾病。目前，經(jīng)顱磁刺激是由醫(yī)生手動放置控制線圈的位置，依賴于醫(yī)生的經(jīng)驗(yàn)，為提高經(jīng)顱磁刺激的治療效果需要準(zhǔn)確把握刺激部位，因此需要將拍攝的MRI影像三維重建之后與患者真實(shí)頭部進(jìn)行配準(zhǔn)。筆者利用點(diǎn)云配準(zhǔn)技術(shù)實(shí)現(xiàn)頭部模型的配準(zhǔn)。

配準(zhǔn)技術(shù)的應(yīng)用最早出現(xiàn)在醫(yī)學(xué)診斷和圖像處理領(lǐng)域[1]，將功能影像和解剖結(jié)構(gòu)影像結(jié)合起來，在同一幅圖像上表達(dá)多種信息，同時反應(yīng)人體的內(nèi)部結(jié)構(gòu)和功能狀態(tài)，幫助醫(yī)生進(jìn)行醫(yī)學(xué)診斷，目前配準(zhǔn)技術(shù)廣泛應(yīng)用于逆向工程、機(jī)器人和虛擬現(xiàn)實(shí)等新的技術(shù)領(lǐng)域。點(diǎn)云配準(zhǔn)是指將不同坐標(biāo)系下的兩片或多片點(diǎn)云數(shù)據(jù)，經(jīng)過一定的旋轉(zhuǎn)和平移變換使其統(tǒng)一到同一坐標(biāo)系下。目前，幾種常用的算法為ICP系列算法(經(jīng)典ICP算法及改進(jìn)ICP算法)、基于特征的配準(zhǔn)方法和基于RANSAC(random sample consensus)的方法。1992年，Besl和Mckay利用四元數(shù)算法實(shí)現(xiàn)了自由曲面點(diǎn)云配準(zhǔn)的迭代最近點(diǎn)算法，即ICP(Iterative Closest Point)算法[2]。ICP算法的基本原理是對于一片點(diǎn)云中的每一個點(diǎn)，在另一片點(diǎn)云中尋找距離該點(diǎn)最近的點(diǎn)作為其對應(yīng)點(diǎn)，形成了對應(yīng)點(diǎn)對，計算對應(yīng)點(diǎn)對的歐式距離的平方的平均值，利用迭代算法使其平均值最小化，并不斷更新兩片點(diǎn)云的相對位置。因此，ICP算法是在不斷重復(fù)“尋找目標(biāo)點(diǎn)云集與源點(diǎn)集中距離最近的點(diǎn)—計算新的剛性變換矩陣并計算匹配誤差”，直到滿足收斂條件才停止迭代。

經(jīng)典ICP算法存在局限性：(1)收斂速度及配準(zhǔn)結(jié)果受初始位置的影響較大；(2)算法容易收斂到局部最優(yōu)而非全局最優(yōu)；(3)在搜索目標(biāo)點(diǎn)的最近點(diǎn)時，每次都必須計算源點(diǎn)集中的每一點(diǎn)，算法的時間復(fù)雜度非常高；(4)經(jīng)典ICP算法要求兩個點(diǎn)云集必須是嚴(yán)格的包含關(guān)系。盡管原始的ICP算法存在一些局限性，但它為點(diǎn)云配準(zhǔn)提供了一種很好的解決思想，因此，許多學(xué)者在經(jīng)典ICP算法的基礎(chǔ)上進(jìn)行改進(jìn)。Zhang等[3]采用k-dtree算法搜索對應(yīng)點(diǎn)；Dorai C等[4]采取隨機(jī)采樣的方式選取對應(yīng)點(diǎn)對，以加快搜索對應(yīng)點(diǎn)的速度。Pulli等[5]采用距離的百分比衡量誤差，將距離過大的10%的點(diǎn)對刪除；Godin等[6]將紋理一致性和頂點(diǎn)的法向一致性作為約束條件，從而刪除無效的對應(yīng)點(diǎn)對，改進(jìn)算法的迭代條件。Chen等[7]用一種基于點(diǎn)到切平面的歐氏距離代替原始ICP算法中基于點(diǎn)對點(diǎn)的歐氏距離；Turk G等[8]提出一種基于點(diǎn)到投影的評價函數(shù)。

1 初始配準(zhǔn)

由于ICP算法容易受初始位置的影響，因此在利用ICP算法精確配準(zhǔn)之前，先對兩片點(diǎn)云進(jìn)行初始配準(zhǔn)。初始配準(zhǔn)是將原本相距甚遠(yuǎn)的兩點(diǎn)云集通過旋轉(zhuǎn)和平移變換使其處在大致相同的位置。

常見的初始配準(zhǔn)方法有：標(biāo)簽法[9]、重心重合法[10]、特征提取法[11-12]。標(biāo)簽法是指在測量時人為設(shè)置一些標(biāo)志點(diǎn)，根據(jù)這些標(biāo)志點(diǎn)進(jìn)行定位，但此方法對儀器及人工經(jīng)驗(yàn)的依賴度較高；重心重合法是計算兩點(diǎn)云的重心，之后將兩重心重合，此方法無法解決旋轉(zhuǎn)錯位；特征提取法是指提取輪廓曲線等，將表面特征相同的部分重合來實(shí)現(xiàn)初始配準(zhǔn)，此方法只適用于具有突出特征的點(diǎn)云。

采用能考慮全局點(diǎn)云的主成分分析法(Principal Component Analysis，PCA)算法[13-15]，通過識別兩點(diǎn)云數(shù)據(jù)的主軸走向和主要形態(tài)特征得到兩點(diǎn)云數(shù)據(jù)的轉(zhuǎn)換矩陣。PCA算法是一種廣泛應(yīng)用的數(shù)據(jù)降維的算法，通過將點(diǎn)云數(shù)據(jù)精簡，留下點(diǎn)云集中貢獻(xiàn)最大的特征。設(shè)一點(diǎn)集為P={p1,p2,…,pn}，首先求其均值并計算點(diǎn)集的協(xié)方差矩陣cov。協(xié)方差矩陣cov的3個特征向量作為點(diǎn)集P的空間直角坐標(biāo)系的3個坐標(biāo)軸，并以均值為坐標(biāo)系的坐標(biāo)原點(diǎn)。計算源點(diǎn)集點(diǎn)云數(shù)據(jù)對應(yīng)目標(biāo)點(diǎn)集點(diǎn)云數(shù)據(jù)的坐標(biāo)變換矩陣參數(shù)，使用坐標(biāo)變換矩陣參數(shù)將源點(diǎn)集數(shù)據(jù)統(tǒng)一變換到目標(biāo)點(diǎn)集數(shù)據(jù)所在的空間直角坐標(biāo)系上。經(jīng)過初始配準(zhǔn)之后，兩點(diǎn)云基本統(tǒng)一到同一坐標(biāo)系下，為之后的ICP算法提供了良好的初始位置。

基于PCA算法的初始配準(zhǔn)步驟如下：

(1)對源點(diǎn)云集和目標(biāo)點(diǎn)云集進(jìn)行矩陣構(gòu)建：

目標(biāo)點(diǎn)云矩陣

(1)

源點(diǎn)云集

(2)

(2)計算兩點(diǎn)云集的質(zhì)心分別為μT和μS：

(3)

(3)分別對兩點(diǎn)云矩陣求協(xié)方差：

目標(biāo)點(diǎn)集點(diǎn)云協(xié)方差矩陣為

(4)

源點(diǎn)云協(xié)方差矩陣為

(5)

其中:X、Y、Z和X′、Y′、Z′分別為目標(biāo)點(diǎn)集矩陣和源點(diǎn)集矩陣的3個列向量。

(4)求解兩點(diǎn)云協(xié)方差矩陣的特征向量和特征值，并將特征值降序排列，得到兩特征向量矩陣T1和T2：

T1=T2*(T*R)

(6)

其中：T表示矩陣T2到矩陣T1的平移矩陣；R表示矩陣T2到矩陣T1 的旋轉(zhuǎn)矩陣。

(7)

至此，完成源點(diǎn)集與目標(biāo)點(diǎn)集的初始配準(zhǔn)過程。

2 ICP算法的解算步驟

S和T分別表示源點(diǎn)集和目標(biāo)點(diǎn)集，可表示為

S={Si}i=1,2,3,…,n

T={Ti}i=1,2,3,…,m

(8)

旋轉(zhuǎn)變換向量用單位四元數(shù)表示為qrotat=[q0,q1,q2,q3]T[1]，并且q0≥0,q02+q12+q22+q32=1,旋轉(zhuǎn)變換矩陣可表示為R(qrotat)。設(shè)平移變換向量為qtrans=[q4,q5,q6]T，則求對應(yīng)點(diǎn)對間的歐式距離之和最小化問題轉(zhuǎn)化為求解qrotat、qtrans使函數(shù)f(qtrans,qrotat)最小化問題。

(9)

算法計算步驟如下：

(1)計算目標(biāo)點(diǎn)集和源點(diǎn)集的重心分別為μT和μS:

(10)

(2)根據(jù)點(diǎn)集S和T構(gòu)造協(xié)方差矩陣：

(11)

(3)根據(jù)協(xié)方差矩陣構(gòu)造4×4：

(12)

(4)計算Q(ΣT,S)的特征值和特征向量，其最大特征值對應(yīng)的特征向量為最佳旋轉(zhuǎn)向量qrotat，則最佳旋轉(zhuǎn)矩陣為：

(13)

(5)計算最佳平移向量為：

qtrans=μT-R(qrotat)μS

(14)

采用均方根誤差(Root Mean Square Error，RMS)進(jìn)行精度評價，均方根誤差可表示為：

(15)

3 實(shí)驗(yàn)

本文中的算法采用Matlab編程實(shí)現(xiàn)，實(shí)驗(yàn)中最大迭代次數(shù)預(yù)設(shè)最大值為200，針對頭部模型點(diǎn)云數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，并通過采樣精簡點(diǎn)云數(shù)量，從而減小程序的運(yùn)行時間。實(shí)驗(yàn)的硬件環(huán)境為i7-8550U 2.00 GHz處理器，8 G內(nèi)存，Win10 64位系統(tǒng)，基于Matlab平臺進(jìn)行測試實(shí)驗(yàn)。

圖1 原始數(shù)據(jù)Fig.1 Head primitive point cloud data

目標(biāo)點(diǎn)集包含68 421個點(diǎn)，源點(diǎn)集包含40 556個點(diǎn)。進(jìn)行PCA初配準(zhǔn)之前，兩點(diǎn)云初始位置如圖1所示，對兩點(diǎn)云集進(jìn)行PCA初始配準(zhǔn)之后，兩點(diǎn)云之間的相對位置如圖2所示。由此可見，PCA初始配準(zhǔn)使兩點(diǎn)云集的位置初步調(diào)整，使其基本處于相同位置。

圖2 PCA初始配準(zhǔn)后Fig.2 Point cloud data after PCA initial registration

選取約100 000個點(diǎn)作為目標(biāo)點(diǎn)集，源點(diǎn)集分別進(jìn)行隨機(jī)采樣精簡點(diǎn)云數(shù)量依次為90 000、80 000、70 000、60 000個點(diǎn)，如圖3所示。目標(biāo)點(diǎn)集與各源點(diǎn)集的初始位置相同進(jìn)行點(diǎn)云配準(zhǔn)，記錄配準(zhǔn)過程所需時間。實(shí)驗(yàn)結(jié)果如表1所示。

圖3 目標(biāo)點(diǎn)集與精簡的源點(diǎn)集 Fig.3 Target point set and reduced source point set

由表1可知，精簡點(diǎn)云數(shù)量能夠減少程序運(yùn)行的時間，但同時均方根誤差RMS在增加，即ICP算法的精度稍微降低。因此在實(shí)際應(yīng)用中可以在保證配準(zhǔn)精度的前提下，合理地精簡點(diǎn)云數(shù)量，有效提高算法的效率。

表1 不同數(shù)量點(diǎn)云ICP算法耗時對比

Table 1 Time-consuming comparison of different point cloud ICP algorithm

點(diǎn)云數(shù)量/個耗時/sRMS/mm90000221.260.090278980000198.530.090499270000162.350.091023660000129.790.0912947

在經(jīng)顱磁刺激治療中，患者真實(shí)頭部點(diǎn)云的獲取欲采用在其頭部選取點(diǎn)，形成能夠代表真實(shí)頭部的點(diǎn)云。因此實(shí)驗(yàn)中采取在頭部模型上沿著眉毛在頭部表面選取一周的點(diǎn)，以頭頂為中心放射狀向下選取點(diǎn)直到與沿眉毛選取的那一周點(diǎn)重合，共選取約200個點(diǎn)形成1個“帽子”形狀的點(diǎn)云集，同時選取易于識別的左右2個眼角及鼻尖3個點(diǎn)作為特征點(diǎn)，如圖4所示。將此“帽子”形狀的稀疏點(diǎn)云數(shù)據(jù)與一個完整的頭部點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn)。經(jīng)過ICP算法配準(zhǔn)之后，如圖5所示，均方根誤差為0.080 0821 mm。

圖4 “帽子”形的頭部點(diǎn)云Fig.4 Cap-shaped point clouds

圖5 ICP配準(zhǔn)Fig.5 Registration of ICP algorithm

4 結(jié)束語

頭部是人體中最重要的組成部分，隨著科技和醫(yī)學(xué)的快速發(fā)展，越來越多的頭部頑疾能夠得以治療，頭部配準(zhǔn)技術(shù)是很多頭部手術(shù)的關(guān)鍵技術(shù)，頭部配準(zhǔn)的精度對治療效果的好壞有直接影響。利用PCA算法進(jìn)行初始配準(zhǔn)，之后用ICP算法進(jìn)行精確配準(zhǔn)。通過精簡點(diǎn)云數(shù)量提高算法的計算速度。同時采集了稀疏點(diǎn)云數(shù)據(jù)與完整模型進(jìn)行了配準(zhǔn)實(shí)驗(yàn)，并取得良好的效果。然而，算法還可以通過采用k-dtree搜尋最近點(diǎn)，利用三維網(wǎng)格分割技術(shù)將目標(biāo)點(diǎn)集中多余的部分刪除來減少運(yùn)行時間，這也將成為下一步研究的內(nèi)容。