專利名稱：基于分級的正反互逆的三維稠密點(diǎn)集快速配準(zhǔn)方法
技術(shù)領(lǐng)域：
本發(fā)明屬于計(jì)算機(jī)應(yīng)用領(lǐng)域，涉及計(jì)算機(jī)視覺和模式識別的一種三維稠密點(diǎn)集對 點(diǎn)集的配準(zhǔn)方法。
背景技術(shù)：
配準(zhǔn)是一個基礎(chǔ)問題，它源自多個領(lǐng)域的很多實(shí)際問題，如不同傳感器獲得的信息融合；不同時(shí)間、條件獲得圖像或模型的差異檢測；成像系統(tǒng)和物體場景變化情況下獲 得的圖像的三維信息獲取；圖像中的模式或目標(biāo)識別等。雖然參與配準(zhǔn)的物體可以有多種 形式，如三維物體的幾何模型、二維或三維圖像，但其實(shí)質(zhì)是兩個二維或三維點(diǎn)集的配準(zhǔn)問 題，即找到兩個點(diǎn)集之間的對應(yīng)關(guān)系。迭代最近點(diǎn)算法(Iterated Closest Point, I CP)是最常用的點(diǎn)集配準(zhǔn)算法，它通 過迭代優(yōu)化矩陣，在每次迭代過程中，對參考點(diǎn)集上的每個點(diǎn)，在目標(biāo)點(diǎn)集中尋找最近點(diǎn)， 并利用這樣的對應(yīng)點(diǎn)，計(jì)算相應(yīng)的旋轉(zhuǎn)矩陣和平移向量，將其用于參考點(diǎn)集上，得到新的參 考點(diǎn)集并進(jìn)入下次迭代過程，最終得到最優(yōu)的轉(zhuǎn)換矩陣，實(shí)現(xiàn)兩點(diǎn)集的配準(zhǔn)。但最初提出的 ICP算法在尋找最近點(diǎn)時(shí)很耗時(shí)，因此很多學(xué)者提出了改進(jìn)的ICP算法以提高配準(zhǔn)效率。ICP算法雖然使用廣泛，但它屬于剛性配準(zhǔn)，在很多應(yīng)用中不能滿足需要，因?yàn)楹?多形變的性質(zhì)是非線性的。因此必須采用非剛性配準(zhǔn)算法，包括基于B樣條或薄板樣條 (Thin PlateSpline,TPS)的彈性配準(zhǔn)算法、基于互信息的彈性配準(zhǔn)方法等。通常一種配準(zhǔn)算法得到的從參考點(diǎn)集到目標(biāo)點(diǎn)集的變形函數(shù)(稱之為正變換)和 從目標(biāo)點(diǎn)集到參考點(diǎn)集的變形函數(shù)(稱之為反變換)是不互逆的。以基于薄板樣條(TPS) 的配準(zhǔn)方法為例，雖然它能通過全局能量最低的限制得到從一個點(diǎn)集到另一個點(diǎn)集的平滑 的變形函數(shù)，但它只在特征點(diǎn)上滿足互逆的條件，在其他頂點(diǎn)不滿足，因此會在一些局部區(qū) 域出現(xiàn)大的配準(zhǔn)誤差。而正反變換互逆是建立兩個點(diǎn)集元素之間一一對應(yīng)關(guān)系的必要條 件。為此一些學(xué)者展開了對正反互逆的配準(zhǔn)算法的研究。提出了一些有效的方法，可實(shí)現(xiàn) 兩個點(diǎn)集的精確配準(zhǔn)。但現(xiàn)有方法幾乎都涉及逆函數(shù)求解，由于處理對象的離散性，求解過 程耗時(shí)長，在很多具有實(shí)時(shí)性要求的領(lǐng)域和行業(yè)中得不到應(yīng)用。雖然已有一些關(guān)于如何提 高正反互逆配準(zhǔn)效率的研究，但目前結(jié)果還不能滿足實(shí)際應(yīng)用要求。

發(fā)明內(nèi)容
為克服上述現(xiàn)有技術(shù)存在的不足，本發(fā)明提供一種快速而又準(zhǔn)確度高的三維點(diǎn)集 對點(diǎn)集的配準(zhǔn)方法，因?yàn)槎S是三維的特例，所以也支持二維點(diǎn)集間的配準(zhǔn)，主要用于三維 人臉識別、CAD產(chǎn)品質(zhì)量檢測、醫(yī)學(xué)圖像融合和分析、三維模型統(tǒng)計(jì)等。本發(fā)明欲解決的技術(shù)問題是給定兩個點(diǎn)集，其中一個點(diǎn)集為參考點(diǎn)集，表示成 P，包含Np個頂點(diǎn)，S卩{Pi，i = 1，...，Np}；另一個點(diǎn)集為目標(biāo)點(diǎn)集，表示成Q，包含Nq個頂點(diǎn)， 即{qj; j = 1，. . .，NtJ。現(xiàn)需快速求解兩個變換函數(shù)；正變換函數(shù)h 從參考點(diǎn)集到目標(biāo)點(diǎn) 集的變換函數(shù)和反變換函數(shù)g 從目標(biāo)點(diǎn)集到參考點(diǎn)集的變換函數(shù)，并要求h g—1，g h—1，其中h—1和g—1分別表示h和g的逆函數(shù)。最終通過這兩個變換函數(shù)，得到兩個點(diǎn)集的精確 的對應(yīng)關(guān)系。為實(shí)現(xiàn)上述目的，本發(fā)明采用了基于全局和局部配準(zhǔn)的分級的方法。全局配準(zhǔn) 為預(yù)配準(zhǔn)，能使參考點(diǎn)集和目標(biāo)點(diǎn)集大體匹配上，但在一些局部區(qū)域會因存在較大的逆 一致性誤差而出現(xiàn)錯配現(xiàn)象，局部配準(zhǔn)即對這些局部區(qū)域進(jìn)行調(diào)整，逐步減小逆一致性 誤差，實(shí)現(xiàn)精確配準(zhǔn)。根據(jù)這一思想。函數(shù)h和g被分解成一系列函數(shù)的復(fù)合，即h = hln ο . . . hn ο hg, g = glno . . . gll ο gg，其中 \和 全局的正反變換函數(shù) ^1Ji = 1，...， η)和gli(i = 1，...，η)為局部的正反變換函數(shù)。此外，為了提高配準(zhǔn)的效率，采用了簡化 的公式求取逆一致性誤差，避免了耗時(shí)的函數(shù)求逆運(yùn)算以及大矩陣運(yùn)算。本發(fā)明算法的主要步驟包括1)輸入兩個點(diǎn)集參考點(diǎn)集P和目標(biāo)點(diǎn)集Q ；2)全局配準(zhǔn)；全局配準(zhǔn)可采用兩種方法薄板樣條函數(shù)(TPS)算法或迭代最近點(diǎn)(ICP)算法。如采用TPS方法，包含兩個子步驟a)從兩個點(diǎn)集中采用人機(jī)交互的方式選取兩組對應(yīng)特征點(diǎn)；b)利用TPS算法和特征點(diǎn)，計(jì)算從參考點(diǎn)集到目標(biāo)點(diǎn)集的變換函數(shù)(正變換hg)， 以及從目標(biāo)點(diǎn)集到參考點(diǎn)集的變換函數(shù)(反變換gg)，并把正反變換分別作用于參考點(diǎn)集和 目標(biāo)點(diǎn)集，得到變形后的兩個點(diǎn)集P’和Q’。如采用ICP方法，也包含兩個子步驟a)判斷待配準(zhǔn)的兩個點(diǎn)集中的頂點(diǎn)數(shù)量是否超過用戶設(shè)定的閾值，如果是從兩個 點(diǎn)集中隨機(jī)選取兩組頂點(diǎn)，組成兩個新的點(diǎn)集，新點(diǎn)集中各頂點(diǎn)間的距離不小于用戶設(shè)定 的距離閾值，且頂點(diǎn)數(shù)量滿足數(shù)量閾值要求；b)利用ICP算法計(jì)算兩個新點(diǎn)集的正變換\和反變換gg，并把正反變換分別作用 于原始參考點(diǎn)集和目標(biāo)點(diǎn)集，得到變形后的兩個點(diǎn)集P’和Q’。3)局部配準(zhǔn)；包括局部特征點(diǎn)的自動定義和局部變形兩個步驟。3.1)局部特征點(diǎn)的自動定義，包含5個子步驟a)求變形后的參考點(diǎn)集上的每一點(diǎn)到目標(biāo)點(diǎn)集上的最近點(diǎn)，如果兩點(diǎn)之間的距離 小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，最終得到一個對應(yīng)點(diǎn)對列表；b)求變形后的目標(biāo)點(diǎn)集上的每一點(diǎn)到參考點(diǎn)集上的最近點(diǎn)，如果兩點(diǎn)之間的距離 小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，最終得到另一個對應(yīng)點(diǎn)對列表；c)根據(jù)前兩個步驟得到的結(jié)果，計(jì)算參考點(diǎn)集和目標(biāo)點(diǎn)集中存在對應(yīng)點(diǎn)的每一個 頂點(diǎn)的逆一致性誤差，并對其對應(yīng)點(diǎn)進(jìn)行調(diào)整；d)把兩個對應(yīng)點(diǎn)對列表中的經(jīng)過調(diào)整操作的對應(yīng)點(diǎn)對挑選出來，存成一個新列 表，并根據(jù)逆一致性誤差從大到小的順序?qū)π铝斜磉M(jìn)行排序；e)從排序后的列表中提取用于局部配準(zhǔn)的兩組特征點(diǎn)。3. 2)局部變形采用具有緊支撐的徑向基函數(shù)(CSRBF)，并根據(jù)3. 1)步驟中得到 的特征點(diǎn)計(jì)算正變換和反變換gli，并把正反變換分別作用于已變形的參考點(diǎn)集P’和目 標(biāo)點(diǎn)集Q’，得到兩個新的變形點(diǎn)集，依然采用P’和Q’表示。4)判斷局部配準(zhǔn)的次數(shù)是否超過用戶設(shè)定的閾值，如果是，則配準(zhǔn)結(jié)束，保存兩個點(diǎn)集中頂點(diǎn)的對應(yīng)關(guān)系；否則返回步驟3)繼續(xù)執(zhí)行局部配準(zhǔn)。與現(xiàn)有技術(shù)相比，本方法的特點(diǎn)在于采用基于全局配準(zhǔn)和局部配準(zhǔn)的分級方式進(jìn) 行配準(zhǔn)，并利用正方向和反方向上對應(yīng)關(guān)系的差異來定義逆一致性誤差和用于局部配準(zhǔn)的 特征點(diǎn)，避免了現(xiàn)有正反互逆配準(zhǔn)算法中復(fù)雜耗時(shí)的函數(shù)求逆運(yùn)算。實(shí)現(xiàn)簡單、效率高，同 時(shí)配準(zhǔn)精度也高。


圖1為本發(fā)明所述一種基于分級的正反互逆的三維稠密點(diǎn)集快速配準(zhǔn)方法的算 法流程圖；圖2為本發(fā)明實(shí)施例1的兩套顱骨點(diǎn)云模型的初始位置；圖3為本發(fā)明實(shí)施例1基于TPS的正方向和反方向的配準(zhǔn)結(jié)果；圖4為本發(fā)明實(shí)施例1基于TPS和一次CSRBF的正方向和反方向的配準(zhǔn)結(jié)果；圖5為本發(fā)明實(shí)施例1基于TPS和五次CSRBF的正方向和反方向的配準(zhǔn)結(jié)果；圖6為本發(fā)明實(shí)施例1經(jīng)TPS和若干次局部配準(zhǔn)后目標(biāo)點(diǎn)集的逆一致性誤差；圖7為本發(fā)明實(shí)施例1經(jīng)TPS和若干次局部配準(zhǔn)后參考點(diǎn)集的逆一致性誤差；圖8為本發(fā)明實(shí)施例1基于ICP和若干次CSRBF的正方向和反方向的配準(zhǔn)結(jié)果；圖9為本發(fā)明實(shí)施例1經(jīng)ICP和若干次局部配準(zhǔn)后目標(biāo)點(diǎn)集的逆一致性誤差；圖10為本發(fā)明實(shí)施例1經(jīng)ICP和若干次局部配準(zhǔn)后參考點(diǎn)集的逆一致性誤差；圖11為本發(fā)明實(shí)施例1的兩套人臉點(diǎn)云模型的初始位置；圖12為本發(fā)明實(shí)施例1人臉在正反兩個方向上的全局配準(zhǔn)結(jié)果和最終結(jié)果。
具體實(shí)施例方式下面結(jié)合具體實(shí)施方式
和優(yōu)選的實(shí)施例作為具體的說明，但不作為對發(fā)明的限 制。實(shí)施例1 一種基于分級的正反互逆的三維稠密點(diǎn)集快速配準(zhǔn)方法，其方法步驟 如圖1所示步驟ISl 給定兩個點(diǎn)集；給定兩個點(diǎn)集，其中一個點(diǎn)集為參考點(diǎn)集，表示成P，另 一個點(diǎn)集為目標(biāo)點(diǎn)集，表示成Q ；步驟2S2 進(jìn)行全局配準(zhǔn)；全局配準(zhǔn)可采用薄板樣條算法(Thin-Plate-Spline，TPS) S200或者迭代最近點(diǎn) 算法(IteratedClosest Point，ICP) S20。其中前者在配準(zhǔn)過程中需要人機(jī)交互，后者是自 動的。1)基于TPS的全局配S200薄板樣條配準(zhǔn)算法(TPS)屬于非剛性配準(zhǔn)算法，其函數(shù)f由兩部分組成，第一部分 由徑向基函數(shù)表示的彈性變換構(gòu)成；第二部分為全局仿射變換200。具體公式如下f(x, γ, ζ) = Rs (χ, y, ζ) + φ5 (χ, γ, ζ)= Σ"=ι α·υ^Pi ~ ^ y^ 2)||) + A + βιχ + P^ + PAz ⑴其中η為特征點(diǎn)的個數(shù)，= 丨為特征點(diǎn)乂=(^,，幻和頂點(diǎn)(x，y，z)之間的歐式距離，aji = 1，...η)，^jU = 1,2,3,4)為待求的權(quán)重。對于彈性變換部分，還有四個附加的邊界條件，分別表示如下
(2)TPS在配準(zhǔn)過程中能使變形模型的全局彎曲能量最低(公式(3))，所以TPS被認(rèn)
為是光滑性最好的一種配準(zhǔn)算法之一。
(3)為求解TPS函數(shù)f中的未知變量aji = 1，...n)和3」(j = 1，2，3，4)，需要以 人機(jī)交互的方式在參考點(diǎn)集和目標(biāo)點(diǎn)集中確定兩組對應(yīng)的特征點(diǎn)乂和《G=l，...《)，其中η 由用戶設(shè)定，一般不超過100 ；然后把特征點(diǎn)乂一一映射到對應(yīng)特征點(diǎn)W，如下公式所示
(4)這一求解過程通常表示成矩陣形式


(5)其中K為ηΧη的矩陣，矩陣中的元
P為nX4的矩陣，矩陣中的
元素 一旦確定了權(quán)重aji = 1，...η)和^(j = 1,2,3,4)的值，就得到了全局的正 變換函數(shù)\201，并可根據(jù)公式(1)對參考點(diǎn)集進(jìn)行變形，使之匹配到目標(biāo)點(diǎn)集上202。反 變換函數(shù)gg的求解過程類似。2)基于ICP的全局配準(zhǔn)S20ICP算法的目標(biāo)是求解兩組點(diǎn)集P和Q的坐標(biāo)系之間的3 X 3旋轉(zhuǎn)矩陣R和平移向 量S，使其滿足以下最優(yōu)化條件20 當(dāng)P和Q的頂點(diǎn)數(shù)量相同，并且一一對應(yīng)時(shí)，可采用四元數(shù)奇異值分解(Singular ValueDecomposition, SVD)得到R和S，并可根據(jù)R和S構(gòu)造4X4的齊次變換矩陣T。但 在實(shí)際應(yīng)用中無法保證P和Q滿足上述條件。因此，ICP方法常用迭代方式尋找最優(yōu)配準(zhǔn) 矩陣。具體步驟如下(1)粗配準(zhǔn)通常采用主元分析法(Principal Component Analysis, PCA)求得 P 和Q兩組點(diǎn)集的3個特征向量，結(jié)合各自的重心，得到4個對應(yīng)的特征點(diǎn)對，按照前述對應(yīng) 點(diǎn)集的配準(zhǔn)方法求得初始配準(zhǔn)矩陣Tc^(2)精配準(zhǔn)，包含4個步驟設(shè)k為迭代次數(shù)，Np和Nq分別為點(diǎn)集P和Q的元素個數(shù)，ε為給定精度(例如ε =ΙΟ"10)。①令k := 0,P° = P,T0 = T。，對P。中的每一點(diǎn)
使用kd樹在Q中尋 找與乂歐氏距離最近的點(diǎn)從而構(gòu)成一個點(diǎn)集

滿足條件\器^凡~qj ，更新平均殘差③如果|rk+1-rk|彡ε或者k達(dá)到了迭代次數(shù)上限則算法結(jié)束，Tk · Tk^Ttl即為 所求的最終配準(zhǔn)矩陣；否則繼續(xù)④。④對Qk+1和Dk+1進(jìn)行對應(yīng)點(diǎn)集的SVD法求得齊次坐標(biāo)下的配準(zhǔn)矩陣T，令Tk+1 = T， 令k : = k+Ι，轉(zhuǎn)步驟②。ICP的配準(zhǔn)效率直接與點(diǎn)集中的頂點(diǎn)數(shù)量相關(guān)，當(dāng)頂點(diǎn)數(shù)量超過1萬時(shí)，配準(zhǔn)效率 會很低。為了提高全局配準(zhǔn)的速度，本發(fā)明設(shè)置了頂點(diǎn)數(shù)量上限4000。當(dāng)點(diǎn)集P和Q中 頂點(diǎn)數(shù)量超過上限時(shí)，采用隨機(jī)采樣的方式從兩個點(diǎn)集中選取4000個頂點(diǎn)，組成新的點(diǎn)集 21。其中點(diǎn)集中的每兩個頂點(diǎn)之間的距離不小于距離閾值，該閾值的初始值由用戶設(shè)定，但 隨采樣情況改變，利用新的點(diǎn)集執(zhí)行ICP配準(zhǔn)，得到旋轉(zhuǎn)矩陣R和平移向量S后，即得到全 局的正變換函數(shù)\22。以此函數(shù)作用于原始參考點(diǎn)集，完成全局配準(zhǔn)23。反變換函數(shù)％的 求解過程類似。步驟3S3 經(jīng)過全局配準(zhǔn)后，兩個點(diǎn)集在大部分區(qū)域能匹配上，但在一些局部區(qū) 域存在較大配準(zhǔn)誤差。本發(fā)明采用基于緊支撐的徑向基函數(shù)配準(zhǔn)算法(Compact Support Radial BasisFunctions, CSRBF)對這些局部區(qū)域進(jìn)行局部配準(zhǔn)，最終達(dá)到精確配準(zhǔn)的目 的?；诰o支撐的徑向基函數(shù)配準(zhǔn)算法CSRBF是一種基于特征點(diǎn)的配準(zhǔn)算法，其特征 點(diǎn)根據(jù)簡化的逆一致性誤差自動生成。1)局部特征點(diǎn)的自動確定。局部配準(zhǔn)中的特征點(diǎn)；^和< 0'= 1,...W)根據(jù)逆一致性 誤差自動獲得。通常，點(diǎn)集P中的頂點(diǎn)PiW逆一致性誤差定義為Ih(Pi)I1(Pi) ι |，同理， 點(diǎn)集Q中的頂點(diǎn)QiW逆一致性誤差定義為I |g (Qi)-IT1(Qi) 11。但這樣的定義方式涉及到逆 函數(shù)h—1和g-1的求解。由于點(diǎn)集的離散性，求解過程非常耗時(shí)。為了加快配準(zhǔn)速度，本發(fā)明采用了一種簡單而又直觀的方式求解頂點(diǎn)的逆一致性 誤差。理想狀態(tài)下，點(diǎn)集P中的每一個頂點(diǎn)Pi通過正變換函數(shù)在映射到點(diǎn)集Q中的頂點(diǎn)
同時(shí)頂點(diǎn)Qi也通過反變換函數(shù)映射到頂點(diǎn)Pi,反之亦然。這樣就得到兩個點(diǎn)集的一一對應(yīng) 關(guān)系，并且點(diǎn)集中任何一個頂點(diǎn)的逆一致性誤差為0。但在實(shí)際應(yīng)用中，兩個點(diǎn)集的頂點(diǎn)數(shù) 量不一定相同，頂點(diǎn)之間也不一定能一一對應(yīng)。但逆一致性誤差的本質(zhì)不變。以點(diǎn)集P上 的頂點(diǎn)為例。對于點(diǎn)集P中的一個頂點(diǎn)Pi,可通過查找最近點(diǎn)找到它在點(diǎn)集Q中的對應(yīng)點(diǎn) 同時(shí)點(diǎn)集Q中可能沒有，也可能有幾個頂點(diǎn)，記作點(diǎn)集Q。 ，在點(diǎn)集P上的對應(yīng)點(diǎn)為Pi。這 時(shí)Qi和Qcot之間的差異就可用于衡量點(diǎn)Pi的逆一致性誤差。它們之間的差異越大，則點(diǎn)Pi 的逆一致性誤差就越大。根據(jù)逆一致性誤差自動獲取局部特征點(diǎn)的具體步驟如下(1)求變形后的參考點(diǎn)集上的每一點(diǎn)到目標(biāo)點(diǎn)集上的最近點(diǎn)，如果兩點(diǎn)之間的距 離小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，得到一個對應(yīng)點(diǎn)對列表30。
(2)求變形后的目標(biāo)點(diǎn)集上的每一點(diǎn)到參考點(diǎn)集上的最近點(diǎn)，如果兩點(diǎn)之間的距 離小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，得到另一個對應(yīng)點(diǎn)對列表31。(3)求有對應(yīng)點(diǎn)的頂點(diǎn)的逆一致性誤差并對其對應(yīng)點(diǎn)進(jìn)行調(diào)整32。以點(diǎn)集P中的頂點(diǎn)Pi為例(點(diǎn)集Q中的頂點(diǎn)Qi類似)，設(shè)頂點(diǎn)Pi的在目標(biāo)點(diǎn)集中 的對應(yīng)點(diǎn)為Qi,同時(shí)在目標(biāo)點(diǎn)集中尋找對應(yīng)點(diǎn)為Pi的頂點(diǎn)，組成點(diǎn)集Q。 ，如果Qcot為空，則 Pi的逆一致性誤差為無效值，同時(shí)也無需對Pi的對應(yīng)點(diǎn)進(jìn)行調(diào)整；否則，Api的逆一致性誤 差定義為Qi和Q。mtCT之間的歐式距離，其中Q。mtCT為點(diǎn)集Qcot的重心。同時(shí)在這種情況下， 把頂點(diǎn)Pi的對應(yīng)點(diǎn)從點(diǎn)Qi調(diào)整為Qi和Qcot組成的點(diǎn)集的重心。(4)選出經(jīng)過調(diào)整操作的對應(yīng)點(diǎn)對，并按照逆一致性誤差從大到小的順序?qū)ζ溥M(jìn) 行排序33。(5)從排序后的點(diǎn)對列表中選出特征點(diǎn)片和<(34)0_ = 1，..^)，特征點(diǎn)對在點(diǎn)集？ 中的分布滿足下列條件Pci - Pj > 0.5a, i 本 j其中a為緊支撐的徑向基函數(shù)配準(zhǔn)算法CSRBF的支撐集大小，定義為a = 3. 66 Δ， Δ為步驟(4)中所有對應(yīng)點(diǎn)對在x、y、ζ軸上的坐標(biāo)的最大差值。如在全局配準(zhǔn)中采用的 是TPS算法，則還有附加條件Jp； - P1j J > Q.Sal k, ρ' e Pcand ρ) e P1其中;/7是TPS配準(zhǔn)時(shí)在點(diǎn)集P中手工定義的特征點(diǎn)，k是局部配準(zhǔn)執(zhí)行的次數(shù)。根據(jù)步驟(2)定義的頂點(diǎn)Pi的逆一致性誤差可近視等價(jià)于標(biāo)準(zhǔn)定義 Ih(Pi)I-1 (Pi) I I，其中 h(Pi) qi;因 g(qp Ppqj e Qc。r，所以 g—Hpi) Qcenter。通過這
種簡化定義，避免了復(fù)雜的函數(shù)求逆運(yùn)算，因此提高了配置效率。從步驟(5)得知，局部特征點(diǎn)的數(shù)量m取決于Δ的值。Δ值越大，則CSRBF的支 撐集a就越大，特征點(diǎn)的數(shù)量m就越小。在實(shí)驗(yàn)中，為避免大矩陣運(yùn)算，m的最大值被限定 為 400。此外，值得注意的是，雖然局部正變換和反變換中的特征點(diǎn)的序號和對應(yīng)關(guān)系是 相同的，但兩種變換中的特征點(diǎn)的坐標(biāo)是不同的。用于局部正變換中的特征點(diǎn)#位于變形 后的參考點(diǎn)集P’中，<位于原始的目標(biāo)點(diǎn)集Q中；而用于局部反變換中的特征點(diǎn)片位于原 始的參考點(diǎn)集P中，仏'位于變形后的目標(biāo)點(diǎn)集Q’中。2)本發(fā)明具體采用了 Wendland緊支撐徑向基函數(shù)(Wendland Compact Support RadialBasis Functions, Wendland CSRBF)進(jìn)行局部變形。其中徑向基函數(shù)定義為 Wendland函數(shù)。采用這一函數(shù)，每一個特征點(diǎn)在三維空間中的作用域?yàn)橐粋€半徑可調(diào)的球 體，而不再是整個點(diǎn)集范圍。因此在配準(zhǔn)過程中可以只調(diào)整沒有匹配上的區(qū)域，同時(shí)保持已 經(jīng)配好的區(qū)域。當(dāng)給定空間維數(shù)d，平滑度C2k(R)以及歐式距離r，Wendland函數(shù)Ψd,k(r)
可表示為
=中+V) (6)其中
9
為截?cái)喽囗?xiàng)式，而 為積分運(yùn)算，在公式(6)中需執(zhí)行k次。從公式(6)可以看到，Wendland函數(shù)￥d,k(r)只在r彡1時(shí)有效。其有效范圍可 以縮放至a，縮放后函數(shù)的數(shù)學(xué)屬性保持不變。 對于3維空間和k = 0，1，2的情況，Wendland函數(shù)￥d,k(r)可分別表示如下 在實(shí)驗(yàn)中，我們采用Ψ^ω作為CSRBF中的徑向基函數(shù)，具體公式如下 其中J = (x，少，ζ)是一個頂點(diǎn)，||#-〒|為頂點(diǎn)無到特征點(diǎn)片,乃，勸之間的歐式距
離，m是特征點(diǎn)的個數(shù)，Cii (i = l，...m)是未知權(quán)重。同TPS類似，這些權(quán)重可通過把參考 點(diǎn)集上的特征點(diǎn)片一一映射到目標(biāo)點(diǎn)集上的對應(yīng)特征點(diǎn)<求解，如下所示 但由于局部配準(zhǔn)中的特征點(diǎn)和<都是估計(jì)得到，可能并不完全對應(yīng)。因此，本發(fā) 明增加了一個誤差閾值ξ，允許特征點(diǎn)在對應(yīng)時(shí)存在一定的誤差 寫成矩陣形式 其中K為mXm的矩陣，矩陣中的元素~
I為單位矩陣，
；I 二 Vf/肌，在本發(fā)明中直接設(shè)置為
, a是Ψ^ΟΟ的支撐集大小
權(quán)重a i(i = 1，. . .m)的值一旦確定，就得到了局部的正變換函數(shù)，并可根據(jù)公 式(7)對參考點(diǎn)集進(jìn)行局部變形，使之更緊密地匹配到目標(biāo)點(diǎn)集上35。局部的反變換函數(shù) gli的求解過程類似。步驟4S4 局部配準(zhǔn)可執(zhí)行多次，直到執(zhí)行次數(shù)超過用戶設(shè)定的閾值。通常對于頂 點(diǎn)數(shù)量較小的點(diǎn)集，執(zhí)行2次局部配準(zhǔn)就能得到很好的結(jié)果，而對于復(fù)雜的點(diǎn)集，一般不超 過5次。步驟5S5 如果大于設(shè)定的閥值，配準(zhǔn)結(jié)束，最近點(diǎn)即為對應(yīng)點(diǎn)。將此方法應(yīng)用于兩個顱骨點(diǎn)云模型和兩個人臉點(diǎn)云模型的配準(zhǔn)，參見圖2和圖11。這些模型都是從CT圖像經(jīng)分割得到。其中目標(biāo)顱骨的頂點(diǎn)數(shù)為166，757，參考顱骨的 頂點(diǎn)數(shù)為166，543;目標(biāo)人臉的頂點(diǎn)數(shù)為62，589，參考人臉的定點(diǎn)數(shù)為59，901。所有實(shí)驗(yàn) 都在一臺HPxw4400工作站上進(jìn)行，該工作站具有頻率為1. 86GHz的Intel (R)Xeon (R) CPU, 3G內(nèi)存。圖3是兩個顱骨點(diǎn)云模型在正反兩個方向經(jīng)TPS全局配準(zhǔn)后的結(jié)果。其中第一列 為TPS中采用的兩組特征點(diǎn)，第二列為正方向上的配準(zhǔn)的結(jié)果，第三列是反方向上的配準(zhǔn) 結(jié)果。由于顱骨頂點(diǎn)數(shù)過大，如采用點(diǎn)繪制很多細(xì)節(jié)看不清楚，因此我們把點(diǎn)云模型進(jìn)行了 三角化操作，并采用面繪制的方法繪制其中一個顱骨，然后采用點(diǎn)繪制或面繪制的方法繪 制另一個顱骨。從圖3可以看出，經(jīng)全局配準(zhǔn)后，在正反兩個方向上兩個顱骨都能基本上匹 配上，但在一些局部區(qū)域，如頭頂和邊緣區(qū)域，兩個顱骨還存在一些差異。圖4顯示了在TPS全局配準(zhǔn)的基礎(chǔ)上再執(zhí)行一次CSRBF局部配準(zhǔn)后的結(jié)果。其中 第一列為125對自動生成的用于CSRBF配準(zhǔn)的特征點(diǎn)，第二列為正方向上的配準(zhǔn)的結(jié)果，第 三列是反方向上的配準(zhǔn)結(jié)果。經(jīng)一次局部調(diào)整后，兩個顱骨在頂部和背部區(qū)域的差異明顯 減小。圖5顯示了最終的配準(zhǔn)結(jié)果。可以看出，經(jīng)過五次局部配準(zhǔn)后，兩個顱骨無論在正 方向上，還是反方向上幾乎都能完全匹配上。與圖3相比，邊緣區(qū)域部分的配準(zhǔn)準(zhǔn)確度有很 大提高。圖6和7分別顯示了經(jīng)TPS和0次/I次/3次/5次局部配準(zhǔn)后目標(biāo)點(diǎn)集和參考點(diǎn) 集的逆一致性誤差。如按照3. 3. 2的逆一致性誤差定義，有些頂點(diǎn)的誤差會為無效值。因 此，我們采用了另一種定義來求參考點(diǎn)集上頂點(diǎn)Pi的誤差| |Pi-goh(Pi) | | ；同理，目標(biāo)點(diǎn) 集上頂點(diǎn)屮的誤差為Ik-hog (qi) ||。在實(shí)驗(yàn)中，誤差的單位為毫米。從結(jié)果可以看出， 隨著局部配準(zhǔn)的不斷執(zhí)行，誤差不斷減小，最后誤差基本控制在1毫米范圍內(nèi)。圖8是全局配準(zhǔn)改用ICP算法后，在正反兩個方向，分別經(jīng)0次、1次和5次CSRBF 的配準(zhǔn)結(jié)果。圖9和10是對應(yīng)的逆一致性誤差結(jié)果。同樣可以看出，誤差隨著局部配準(zhǔn)的 次數(shù)增加而減低。圖12是圖11中顯示的兩套人臉點(diǎn)云模型經(jīng)ICP配準(zhǔn)后的結(jié)果，以及經(jīng)ICP全局 配準(zhǔn)和3次CSRBF局部調(diào)整后的配準(zhǔn)結(jié)果。后者的準(zhǔn)確度明顯比前者高。此外，對配準(zhǔn)耗時(shí)情況也有統(tǒng)計(jì)。其中對于顱骨模型的配準(zhǔn)，采用TPS全局配準(zhǔn)和 5次CSRBF局部配準(zhǔn)共耗時(shí)32. 5秒；采用ICP全局配準(zhǔn)和5次CSRBF局部配準(zhǔn)共耗時(shí)34. 9 秒。人臉模型的配準(zhǔn)采用ICP全局配準(zhǔn)和3次CSRBF局部配準(zhǔn)共耗時(shí)8. 7秒。相對于現(xiàn)有 算法由于涉及函數(shù)求逆操作，處理兩幅分辨率為181X217需近半小時(shí)時(shí)間，我們算法的效 率有很大提高。本方法的特色在于采用基于全局配準(zhǔn)和局部配準(zhǔn)的分級方式進(jìn)行配準(zhǔn)，并利用正 方向和反方向上對應(yīng)關(guān)系的差異來定義逆一致性誤差和用于局部配準(zhǔn)的特征點(diǎn)，避免了現(xiàn) 有正反互逆配準(zhǔn)算法中復(fù)雜耗時(shí)的函數(shù)求逆運(yùn)算。實(shí)現(xiàn)簡單、效率高，同時(shí)配準(zhǔn)精度也高。上述實(shí)驗(yàn)結(jié)果表明本方法對拓?fù)浣Y(jié)構(gòu)復(fù)雜、頂點(diǎn)數(shù)量繁多的點(diǎn)云模型都能快速得 到準(zhǔn)確度很高的配準(zhǔn)結(jié)果，可以用于計(jì)算機(jī)視覺和模式識別等各應(yīng)用領(lǐng)域，具有快速準(zhǔn)確、 操作簡單、應(yīng)用前景廣的特點(diǎn)。最后所應(yīng)說明的是，以上實(shí)施方式僅用以說明本發(fā)明的技術(shù)方案而非限制。盡管
11參照實(shí)施方式對本發(fā)明進(jìn)行了詳細(xì)說明，本領(lǐng)域的普通技術(shù)人員應(yīng)當(dāng)理解，對本發(fā)明的技 術(shù)方案進(jìn)行修改或者等同替換，都不脫離本發(fā)明技術(shù)方案的精神和范圍，其均應(yīng)涵蓋在本 發(fā)明的權(quán)利要求范圍當(dāng)中。
權(quán)利要求
一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征在于包括下列步驟步驟一(S1)輸入兩個點(diǎn)集，一個參考集，一個目標(biāo)集；步驟二(S2)全局配準(zhǔn)；步驟三(S3)局部配準(zhǔn)；步驟四(S4)判斷局部配準(zhǔn)次數(shù)；步驟五(S5)配準(zhǔn)結(jié)束，最近點(diǎn)即為對應(yīng)點(diǎn)；其中，步驟二(S2)所述的全局配準(zhǔn)可采用兩種方法，薄板樣條函數(shù)算法(TPS)(S200)或迭代最近點(diǎn)算法(ICP)(S20)；其中，步驟三(S3)所述的局部配準(zhǔn)方法包含六個子步驟1)求變形后的參考點(diǎn)集上的每一點(diǎn)到目標(biāo)點(diǎn)集上的最近點(diǎn)，兩點(diǎn)之間的距離小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，最終得到一個對應(yīng)點(diǎn)對列表(30)；2)求變形后的目標(biāo)點(diǎn)集上的每一點(diǎn)到參考點(diǎn)集上的最近點(diǎn)，兩點(diǎn)之間的距離小于用戶設(shè)定的閾值，則以該最近點(diǎn)作為它的對應(yīng)點(diǎn)，最終得到另一個對應(yīng)點(diǎn)對列表(31)；3)根據(jù)前兩個步驟得到的結(jié)果，計(jì)算參考點(diǎn)集和目標(biāo)點(diǎn)集中存在對應(yīng)點(diǎn)的每一個頂點(diǎn)的逆一致性誤差，并對其對應(yīng)點(diǎn)進(jìn)行調(diào)整(32)；4)把兩個對應(yīng)點(diǎn)對列表中的經(jīng)過調(diào)整操作的對應(yīng)點(diǎn)對挑選出來，存成一個新列表，并根據(jù)逆一致性誤差從大到小的順序?qū)π铝斜磉M(jìn)行排序(33)；5)從排序后的列表中提取用于局部配準(zhǔn)的兩組特征點(diǎn)(34)；6)采用具有緊支撐的徑向基函數(shù)CSRBF，并根據(jù)上一個步驟得到的特征點(diǎn)計(jì)算正變換hli和反變換gli，并把正反變換分別作用于已變形的參考點(diǎn)集P’和目標(biāo)點(diǎn)集Q’，得到兩個新的變形點(diǎn)集，依然采用P’和Q’表示(35)；其中，步驟四(S4)所述的判斷局部配準(zhǔn)次數(shù)小于預(yù)先設(shè)定的閥值的，重復(fù)執(zhí)行步驟三(S3)。
2.根據(jù)權(quán)利要求1所述一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征在 于步驟二(S2)所述的薄板樣條函數(shù)算法(TPS) (S200)全局配準(zhǔn)法包含如下兩個子步驟1)從兩個點(diǎn)集中采用人機(jī)交互的方式選取兩組對應(yīng)特征點(diǎn)(200)；2)利用TPS算法和特征點(diǎn)，計(jì)算從參考點(diǎn)集到目標(biāo)點(diǎn)集的變換，以及從目標(biāo)點(diǎn)集到 參考點(diǎn)集的變換，并把正反變換分別作用于參考點(diǎn)集和目標(biāo)點(diǎn)集，得到變形后的兩個點(diǎn)集 (201-202)。
3.根據(jù)權(quán)利要求1所述一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征在 于步驟二(S2)所述的迭代最近點(diǎn)算法(ICP) (S20)也包含如下兩個子步驟1)判斷參考點(diǎn)集和目標(biāo)點(diǎn)集中的頂點(diǎn)數(shù)量，頂點(diǎn)數(shù)量大于設(shè)定的閥值，則從兩個點(diǎn)集 中隨機(jī)選取兩組頂點(diǎn)，組成兩個新的點(diǎn)集，新的點(diǎn)集中各頂點(diǎn)間的距離大于用戶設(shè)定的距 離閾值，且頂點(diǎn)數(shù)量滿足數(shù)量閾值要求(20-22)；2)利用ICP算法計(jì)算兩個新點(diǎn)集的正變換和反變換，并把正反變換分別作用于參考點(diǎn) 集和目標(biāo)點(diǎn)集，得到變形后的兩個點(diǎn)集(23)。
4.根據(jù)權(quán)利要求1所述的一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征 在于所述步驟三(S3)中的子步驟1) (30)和2) (31)中采用kd樹結(jié)構(gòu)加快最近點(diǎn)求解。
5.根據(jù)權(quán)利要求1所述的一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征在于所述步驟三(S3)中的子步驟3) (32)中頂點(diǎn)的逆一致性誤差的定義分兩種情況情況1 參考點(diǎn)集中頂點(diǎn)Pi在目標(biāo)點(diǎn)集中的對應(yīng)點(diǎn)為頂點(diǎn)同時(shí)目標(biāo)點(diǎn)集中有一個或 幾個頂點(diǎn)，記作點(diǎn)集Q。 ，其在參考點(diǎn)集中的對應(yīng)點(diǎn)為頂點(diǎn)Pi,此時(shí)頂Api的逆一致性誤差 定義為Qi和Qrente之間的歐式距離，其中Qcenter為點(diǎn)集Qcot的重心；情況2 參考點(diǎn)集中頂點(diǎn)Pi在目標(biāo)點(diǎn)集中的對應(yīng)點(diǎn)為頂點(diǎn)但目標(biāo)點(diǎn)集中沒有頂點(diǎn)的 對應(yīng)點(diǎn)為頂點(diǎn)Pi，此時(shí)頂點(diǎn)Pi的逆一致性誤差定義為無效值； 上述定義也適應(yīng)目標(biāo)點(diǎn)集上的頂點(diǎn)。
6.根據(jù)權(quán)利要求1所述的一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征 在于所述步驟三(S3)中的子步驟3) (32)步驟中對應(yīng)點(diǎn)的調(diào)整只在一種情況下執(zhí)行即當(dāng)參考點(diǎn)集中頂點(diǎn)Pi在目標(biāo)點(diǎn)集中的對應(yīng)點(diǎn)為頂點(diǎn)Qi,同時(shí)目標(biāo)點(diǎn)集中有一個或幾 個頂點(diǎn)，記作點(diǎn)集Q。 ，其在參考點(diǎn)集中的對應(yīng)點(diǎn)為頂點(diǎn)Pi,此時(shí)頂點(diǎn)Pi的對應(yīng)點(diǎn)調(diào)整為^ 和Qcot組成的點(diǎn)集的重心；上述對 應(yīng)點(diǎn)的調(diào)整情況也適應(yīng)目標(biāo)點(diǎn)集。
7.根據(jù)權(quán)利要求1所述的一種基于分級的正反互逆的三維稠密點(diǎn)集配準(zhǔn)方法，其特征 在于所述步驟三(S3)中的子步驟5) (34)中從參考點(diǎn)集中提取出來的特征點(diǎn)之間的距離 大于0.5a，其中a為緊支撐徑向基函數(shù)CSRBF的緊支撐集大小，定義為a = 3.66 Δ，Δ為 列表中對應(yīng)點(diǎn)對在χ、y、ζ軸上的坐標(biāo)的最大差值。
全文摘要
本發(fā)明提供一種基于分級的正反互逆的三維稠密點(diǎn)集快速配準(zhǔn)方法，解決的技術(shù)問題是給定兩個點(diǎn)集，其中一個為參考點(diǎn)集，另一個為目標(biāo)點(diǎn)集，現(xiàn)需求解兩個變換函數(shù)；正變換函數(shù)h從參考點(diǎn)集到目標(biāo)點(diǎn)集的變換函數(shù)，反變換函數(shù)g從目標(biāo)點(diǎn)集到參考點(diǎn)集的變換函數(shù)，并要求h≈g-1，g≈h-1，其中h-1和g-1分別表示h和g的逆函數(shù)。最終通過這兩個變換函數(shù)，得到兩個點(diǎn)集的精確的對應(yīng)關(guān)系。與現(xiàn)有技術(shù)相比，本方法的優(yōu)點(diǎn)在于基于全局配準(zhǔn)和局部配準(zhǔn)的分級方式進(jìn)行配準(zhǔn)，并利用正方向和反方向上對應(yīng)關(guān)系的差異來定義逆一致性誤差和用于局部配準(zhǔn)的特征點(diǎn)，避免了現(xiàn)有正反互逆配準(zhǔn)算法中復(fù)雜耗時(shí)的函數(shù)求逆運(yùn)算。實(shí)現(xiàn)簡單、效率高，同時(shí)配準(zhǔn)精度也高。
文檔編號G06K9/64GK101887525SQ20101022525
公開日2010年11月17日 申請日期2010年7月9日 優(yōu)先權(quán)日2010年7月9日
發(fā)明者周明全, 耿國華, 鄧擎瓊 申請人:北京師范大學(xué)