基于ICP算法的頭部模型配準方法研究

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

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

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

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

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

1 初始配準

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

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

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

基于PCA算法的初始配準步驟如下：

(1)對源點云集和目標點云集進行矩陣構建：

目標點云矩陣

(1)

源點云集

(2)

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

(3)

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

目標點集點云協(xié)方差矩陣為

(4)

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

(5)

其中:X、Y、Z和X′、Y′、Z′分別為目標點集矩陣和源點集矩陣的3個列向量。

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

T1=T2*(T*R)

(6)

其中：T表示矩陣T2到矩陣T1的平移矩陣；R表示矩陣T2到矩陣T1 的旋轉矩陣。

(7)

至此，完成源點集與目標點集的初始配準過程。

2 ICP算法的解算步驟

S和T分別表示源點集和目標點集，可表示為

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

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

(8)

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

(9)

算法計算步驟如下：

(1)計算目標點集和源點集的重心分別為μT和μS:

(10)

(2)根據(jù)點集S和T構造協(xié)方差矩陣：

(11)

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

(12)

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

(13)

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

qtrans=μT-R(qrotat)μS

(14)

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

(15)

3 實驗

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

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

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

圖2 PCA初始配準后Fig.2 Point cloud data after PCA initial registration

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

圖3 目標點集與精簡的源點集 Fig.3 Target point set and reduced source point set

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

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

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

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

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

圖4 “帽子”形的頭部點云Fig.4 Cap-shaped point clouds

圖5 ICP配準Fig.5 Registration of ICP algorithm

4 結束語

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