湯春明,牛百鈴,李光旭
(天津工業(yè)大學(xué)電子與信息工程學(xué)院,天津 300387)
?
基于球面保角映射理論的表面特征標(biāo)記點(diǎn)匹配法
湯春明,牛百鈴,李光旭
(天津工業(yè)大學(xué)電子與信息工程學(xué)院,天津300387)
摘要:提出一種基于曲面參數(shù)化理論的半自動(dòng)匹配方法.首先,選擇任一樣本數(shù)據(jù)作為參照模型,根據(jù)其表面幾何特征手動(dòng)標(biāo)定特征標(biāo)記點(diǎn);其次,利用保角映射方法將全部樣本數(shù)據(jù)映射到同一球形域,根據(jù)球面坐標(biāo)對(duì)應(yīng)關(guān)系實(shí)現(xiàn)在其他樣本表面上特征標(biāo)記點(diǎn)的自動(dòng)標(biāo)定;最后,利用迭代逼近點(diǎn)方法歸一化特征標(biāo)記點(diǎn)在空間的位置.實(shí)驗(yàn)考察了統(tǒng)計(jì)形狀模型的通用性原則和專一性原則,并且比較了球面保角映射中不同約束條件對(duì)獲得模型質(zhì)量的影響.實(shí)驗(yàn)結(jié)果表明:利用球面保角映射3個(gè)基準(zhǔn)點(diǎn)約束所得模型的專一性指標(biāo)明顯好于零質(zhì)心約束方法,并且實(shí)驗(yàn)所需時(shí)間減少近50%.
關(guān)鍵詞:球面保角映射理論;特征標(biāo)記點(diǎn)匹配;統(tǒng)計(jì)形狀模型構(gòu)建;活動(dòng)形狀模型;圖像分割
model;image segmentation
由于將對(duì)象形狀信息作為先驗(yàn)知識(shí),統(tǒng)計(jì)形狀模型(statistical shape model,SSM)已廣泛應(yīng)用于計(jì)算機(jī)醫(yī)學(xué)輔助診斷(computer aided diagnosis,CAD)領(lǐng)域.常見的統(tǒng)計(jì)形狀模型包括點(diǎn)分布模型(point distribution module,PDM)[1]和概率地圖集PA[2].其中,點(diǎn)分布模型由于表達(dá)簡(jiǎn)潔、魯棒性強(qiáng)、豐富的顯示算法支持而常被用于醫(yī)學(xué)圖像分割、結(jié)構(gòu)分析、形狀再現(xiàn)等.
點(diǎn)分布模型是一種用來描述目標(biāo)輪廓的可變形數(shù)學(xué)模型.首先在樣本圖像的目標(biāo)輪廓上選取特征標(biāo)記點(diǎn)構(gòu)成訓(xùn)練集,然后利用主成份分析法對(duì)其進(jìn)行統(tǒng)計(jì)分析,得到點(diǎn)分布統(tǒng)計(jì)模型.該模型主要包括2部分:第1部分是特征點(diǎn)組成的平均值,即平均輪廓;第2部分是特征點(diǎn)的形變方式,即特征點(diǎn)相對(duì)于平均輪廓的整體變化趨勢(shì).形變范圍通過模型參數(shù)的變化范圍加以限制,避免任意變形的可能性.
在不同的樣本圖像的目標(biāo)輪廓之間選擇對(duì)應(yīng)的,具有相同局部幾何特征的特征標(biāo)記點(diǎn)是關(guān)系點(diǎn)分布模型質(zhì)量的關(guān)鍵要素之一[3].這種對(duì)應(yīng)點(diǎn)的匹配,最簡(jiǎn)單的方法是在每個(gè)例子中選擇一個(gè)起點(diǎn),以同等的距離放置并且使每個(gè)邊界或表面的點(diǎn)數(shù)目相等.該方法在文獻(xiàn)[4]中提出,然而,這種等距離的放置方法并不能給出一個(gè)合理的分組對(duì)應(yīng),在訓(xùn)練集中等效點(diǎn)的相對(duì)位置可能有很大的不同.用這些點(diǎn)訓(xùn)練得到的統(tǒng)計(jì)形狀模型質(zhì)量較差.另一種方法是提取表面的某些幾何特性,如曲率或脊線形狀[5-6].這些功能之間的對(duì)應(yīng)關(guān)系可以用一個(gè)通用的數(shù)值優(yōu)化算法建立[7].但是這些參考局部特征方法不能達(dá)到全局最優(yōu).另一類是考慮全局優(yōu)化策略,如球面諧波(SPHARM)映射[8],重參數(shù)化方法[5,9].本文以臟器圖像作為研究對(duì)象,使用重參數(shù)化方法,引入球面保角映射作為零虧格表面模型的映射域,并比較了在球面保角映射中的2個(gè)約束條件對(duì)獲得模型質(zhì)量的影響.
1.1匹配過程
本文方法為半自動(dòng)方法.首先選擇一個(gè)表面樣本作為參照模型,在其表面手動(dòng)標(biāo)定特征標(biāo)記點(diǎn).然后將所有表面樣本映射到同一球面域中,根據(jù)球面坐標(biāo)位置關(guān)系,自動(dòng)在其他表面上找出特征標(biāo)記點(diǎn)位置,如圖1所示.“映射(mapping)”是指三維肺部表面特征標(biāo)記點(diǎn)映射到球形參數(shù)域,“放置(placing)”表示圖1中特征標(biāo)記點(diǎn)放置的步驟順序.首先計(jì)算球面形映射,把三維肺部表面特征標(biāo)記點(diǎn)映射到單位球形結(jié)構(gòu)域中.假設(shè)形狀的旋轉(zhuǎn)變換可以獨(dú)立完成.因此,根據(jù)曲面曲率的參數(shù)域分布,可以匹配設(shè)置參考模型的所有訓(xùn)練集.同時(shí),參考映射到球形參數(shù)域特征標(biāo)記點(diǎn)的模型,完成三維肺部訓(xùn)練集的對(duì)應(yīng)關(guān)系.
圖1 特征標(biāo)志點(diǎn)匹配過程Fig.1 Outline of landmarks matching
1.2球面保角映射
許多臟器的形狀是屬零拓?fù)浣Y(jié)構(gòu),即封閉無孔的表面或自相交.由于具有相同的拓?fù)浣Y(jié)構(gòu),可以把一個(gè)封閉的零虧格表面調(diào)和映射到單位球面[10].如圖2所示的三角網(wǎng)格曲面,假設(shè)K為簡(jiǎn)單復(fù)合體,u、v為頂點(diǎn),euv為跨越u、v的邊緣線.在三維空間表面的位置矢量r(u,v)通常被表示為
使用f來表示K上定義的分段線性函數(shù).能量函數(shù)被定義為
圖2 三角網(wǎng)表面中相鄰2個(gè)單元Fig.2 Triangular mesh surface adjacent two units
字符串常量kuv描述了表面和參數(shù)域之間的變化因素.改變它可以定義不同段的能量.如果kuv= 1,該段能量稱為Tuette能量.字符串常量定義的諧波能量計(jì)算為
式中:αij和βij分別為相對(duì)于原三角曲面邊緣的2個(gè)相對(duì)角度.
文獻(xiàn)中有許多不同的方法來處理封閉零虧格表面,如諧波能量最小化、球面參數(shù)化和拉普拉斯等[11-12].球形映射是從一個(gè)封閉的零規(guī)格表面映射到單位球面的調(diào)和映射. Gu等[13]提出了一個(gè)快速下降算法使得到的諧波能量最小化和采用零質(zhì)量中心約束指定映射結(jié)果.然而,通過使用快速下降法的共形映射的結(jié)構(gòu)不是唯一的,會(huì)形成莫比烏斯帶群.莫比烏斯帶變換不同的約束將導(dǎo)致不同的映射結(jié)果.本文比較分別采用零質(zhì)心約束和3個(gè)基準(zhǔn)點(diǎn)約束的配準(zhǔn)結(jié)果.
1.3球面保角映射到約束條件
本文將把三維肺部表面特征標(biāo)記點(diǎn)映射到二維球面域.由于映射不是唯一的,球形映射有6個(gè)自由度并且所有球形的映射形式構(gòu)成了莫比烏斯變換,稱為莫比烏斯映射.它在復(fù)平面C中定義為
莫比烏斯變換是保形變換.確保共形映射是唯一的及算法的收斂性,需要添加額外的約束.在文獻(xiàn)[13]中,Gu等使球形映射圖的質(zhì)心為原點(diǎn):
式中:M2為網(wǎng)格M1的圖;δM1為M1的區(qū)域元素.此約束確定訓(xùn)練樣本是唯一的旋轉(zhuǎn),稱此為零質(zhì)心約束.
同時(shí),另一種約束形狀變化的方法是在表面固定3個(gè)基準(zhǔn)點(diǎn).給出3個(gè)不同的點(diǎn)Z1,Z2,Z3∈C,則莫比烏斯變換分別映射到(0,1,∞)
為了保證圖形的對(duì)稱性,選擇P0、Pinit作為中心軸和表面的交叉點(diǎn),Z1、Z3分別為它們的映射點(diǎn). Z2是P1的復(fù)數(shù)坐標(biāo)點(diǎn),P1被選擇為x軸和表面的交點(diǎn),顯示在圖3(a)中.點(diǎn)P0,P1,Pinit分別映射到南極,CIO是北極.該算法如下:
輸入:網(wǎng)格M1,能量誤差閾值δE,在頂點(diǎn)上的3個(gè)基準(zhǔn)點(diǎn)P0,P1,Pinit;輸出:3個(gè)點(diǎn)P0,P1,Pinit分別被映射到南極,CIO和北極的球面共性圖M2.
(a)計(jì)算高斯圖.
(b)最小化Tuette能量得到Tuette映射.
(c)最小化共性能量得到球面共性映射.
(d)由公式(16)計(jì)算球面共形圖的立體投影.
(e)Z1,Z2,Z3分別是P0,P1,Pinit的投影,計(jì)算全部頂點(diǎn)的高斯變換.
比較這2種約束下的匹配結(jié)果.當(dāng)繪制球面上的經(jīng)度線和緯度線時(shí),可以在原始表面上進(jìn)一步匹配,如圖3(a)和(c)采用零質(zhì)量中心約束的肺部,圖3(b)和(d)是采用3個(gè)基準(zhǔn)點(diǎn)約束的肺部.
圖3 莫比烏斯變換下不同的約束條件相同的肺部樣本映射圖Fig.3 Mobius transformation under same constraints of different lung sample map
1.4特征標(biāo)記點(diǎn)的配準(zhǔn)
圖像配準(zhǔn)實(shí)際就是通過分析待配準(zhǔn)圖像中存在的幾何物體的形狀,選擇其中一種最合適的空間變換,使2幅圖像中相對(duì)應(yīng)部位可以一一對(duì)應(yīng)起來[14].所以說,圖像配準(zhǔn)最基本的問題就是如何找到一種圖像變換方法,這是一個(gè)尋找最優(yōu)解的問題.
在前面的小節(jié)中,保角映射保證了原有的表面點(diǎn)的相對(duì)位置.找到坐標(biāo)點(diǎn)的參考模型,首先,用參考模型映射的相同方式將訓(xùn)練樣本映射到球形域,然后比較參考,根據(jù)球面坐標(biāo)(θ,準(zhǔn))獲得每一個(gè)點(diǎn)集相應(yīng)的位置.然而,相應(yīng)的位置通常不位于頂點(diǎn),所以必須在物體表面周圍的參考標(biāo)記的位置找到近似的頂點(diǎn).一種有序的頂點(diǎn)智能搜索方法是近鄰搜索(NNS),這個(gè)搜索可以用樹的特性迅速消除大部分的搜索空間,從而有效地進(jìn)行.
1.5訓(xùn)練模型一致化
獲得每個(gè)樣本的特征標(biāo)記點(diǎn)后,需要將這些樣本的位置、大小在三維空間上進(jìn)行統(tǒng)一.本文采用迭代逼近點(diǎn)法(Iterative Closest Point,ICP)實(shí)現(xiàn).假設(shè)用表示空間第1個(gè)點(diǎn)集,第2個(gè)點(diǎn)集的對(duì)齊匹配轉(zhuǎn)換為使下式的目標(biāo)函數(shù)最小.
ICP算法本質(zhì)上是基于最佳匹配方法的最小二乘法,通過重復(fù)點(diǎn)的集合,以計(jì)算出最佳的剛性變換過程,直到最后達(dá)到最優(yōu)收斂準(zhǔn)則. ICP算法主要是找到目標(biāo)點(diǎn)集與參考點(diǎn)之間的旋轉(zhuǎn)R和平移Τ變換,使得兩匹配數(shù)據(jù)滿足某種程度上的最優(yōu)匹配.假設(shè)目標(biāo)點(diǎn)集P的坐標(biāo)及參考點(diǎn)集Q的坐標(biāo)為,在第k次迭代中計(jì)算與點(diǎn)集P相對(duì)應(yīng)的坐標(biāo)為{QkQk∈R,i = 1,2,…,N},計(jì)算P與Qk之間的變換矩陣并對(duì)原變換進(jìn)行更新,直到數(shù)據(jù)間平均距離小于給定閾值子,即滿足式(7)最小.步驟如下:
(1)在目標(biāo)點(diǎn)集P中取點(diǎn)集Pki∈Pk;
(2)在參考點(diǎn)集Q中取對(duì)應(yīng)點(diǎn)集Qki∈Qk,并計(jì)算使
(3)求得旋轉(zhuǎn)矩陣Rk與平移向量Tk,使得
(4)計(jì)算
(5)如果dk+1大于或等于給定的子則返回到(2),直到dk + 1<子或迭代次數(shù)大于預(yù)設(shè)的最大迭代次數(shù)為止.
對(duì)于ICP每一次的迭代,最小化對(duì)應(yīng)點(diǎn)的均方差均使得點(diǎn)集Pki離Qki更近,而Qki則是Pki在Qi的最近點(diǎn).因此,每次的迭代會(huì)使得Pi離Qi更近.
2.1實(shí)驗(yàn)數(shù)據(jù)
本研究采用人體肺部的CT圖像作為研究對(duì)象.除生理器官本身的差異以外,受呼吸作用影響肺部的形變量較大.本研究采用的86例肺部圖像數(shù)據(jù)由The Lung Image Database Consortium(LIDC)提供.實(shí)驗(yàn)比較了在物體表面均勻分布特征標(biāo)記點(diǎn)方法和提案方法在模型質(zhì)量上的差異.實(shí)驗(yàn)過程及結(jié)果如圖4所示.
圖4 實(shí)驗(yàn)過程及結(jié)果Fig.4 Experimental procedure and results
2.2評(píng)價(jià)指標(biāo)
模型的通用性能力用來檢測(cè)所生成的形狀模型對(duì)合理的新形狀實(shí)例表達(dá)的能力,即模型能夠描述任一對(duì)象目標(biāo)實(shí)例,不僅限于那些在訓(xùn)練集中出現(xiàn)的.它通常用leave-one-out算法執(zhí)行.首先在N例訓(xùn)練中選出N-1例用于訓(xùn)練統(tǒng)計(jì)形狀模型.然后用所得到的統(tǒng)計(jì)模型去和剩余一例進(jìn)行匹配.此過程重復(fù)執(zhí)行,其匹配誤差的平均值作為此統(tǒng)計(jì)形狀模型形狀表達(dá)能力的評(píng)價(jià)指標(biāo).通用性能力被定義為一個(gè)關(guān)于形狀特征矢量M的方程.例如,用2種方法獲得的統(tǒng)計(jì)形狀模型A和B,對(duì)于大多數(shù)的形狀變量M有GA(M)≤GB(M),則表示A模型優(yōu)于B模型.
模型的專一性能力是用來估計(jì)產(chǎn)生形狀有效性的能力,即利用模型描述的形狀皆為對(duì)象物體的合法實(shí)例.首先,利用統(tǒng)計(jì)形狀模型產(chǎn)生一系列形狀實(shí)例.然后,將產(chǎn)生的形狀實(shí)例與訓(xùn)練集中的訓(xùn)練樣本進(jìn)行比較.則模型的專一性定義為
式中:xj代表訓(xùn)練樣本產(chǎn)生的形狀實(shí)例.同樣,若SA(M)≤SB(M)則說明模型A更具專一性.關(guān)于模型通用性能力和專一性能力的實(shí)現(xiàn)方法參見文獻(xiàn)[15].
2.3實(shí)驗(yàn)結(jié)果與分析
本文首先對(duì)86個(gè)肺部模型的通用性能力和專一性能力進(jìn)行評(píng)價(jià),結(jié)果如圖5所示.由圖5可知,通用性模型能夠描述任意的目標(biāo)實(shí)例,而不僅僅是那些訓(xùn)練集中出現(xiàn)的.專一性模型僅僅能夠描繪有效的目標(biāo)實(shí)例,有簡(jiǎn)潔性.
圖5 模型質(zhì)量的評(píng)價(jià)Fig.5 Quality evaluation model
再對(duì)20個(gè)肺部模型的通用性能力和專一性能力進(jìn)行衡量,結(jié)果如圖6所示.
圖6的比較結(jié)果表明:雖然在模型通用性上,利用球面保角映射零質(zhì)心約束得到的模型并沒有明顯的優(yōu)越性,但是,利用球面保角映射3個(gè)基準(zhǔn)點(diǎn)約束所得模型的專一性指標(biāo)明顯好于零質(zhì)心約束方法,并且實(shí)驗(yàn)所需時(shí)間減少近50%.一般來說,特征標(biāo)記點(diǎn)的匹配程度越高,特征點(diǎn)分布越“集中”,所獲得的統(tǒng)計(jì)形狀模型有效變化范圍越小,即產(chǎn)生的樣本實(shí)例越“緊湊”.因此,由于采用3個(gè)基準(zhǔn)點(diǎn)約束方法獲得的特征點(diǎn)匹配程度高于零質(zhì)心約束方法,“模型專一性能力”也隨之得到增強(qiáng).
在計(jì)算效率上3個(gè)基準(zhǔn)點(diǎn)約束方法也具有很好的可推廣性.實(shí)踐中,利用通用計(jì)算機(jī)(CPU:Xeon E5-1607 v2,內(nèi)存:8 G)完成一例表面數(shù)據(jù)配準(zhǔn)的執(zhí)行時(shí)間通常少于30 min.
圖6 模型質(zhì)量的比較結(jié)果Fig.6 Comparison of results of model quality
自動(dòng)進(jìn)行表面模型的特征標(biāo)記點(diǎn)匹配是建立統(tǒng)計(jì)形狀模型的第一步.本文提出了利用表面變量化方式和球面保角映射來實(shí)現(xiàn)三維三角網(wǎng)表面的特征標(biāo)記點(diǎn)匹配問題.本研究利用表面數(shù)據(jù)的局部平均曲率和形狀模型的緊致性特性作為對(duì)應(yīng)特征的選取原則.由于3個(gè)基準(zhǔn)點(diǎn)約束方法綜合考慮了形狀的局部特性和整體特性,優(yōu)化過程具有一定的穩(wěn)定性來抑制優(yōu)化算法的局部極小問題.
此外,標(biāo)記點(diǎn)的數(shù)目會(huì)影響模型表達(dá)能力,本文用1 069個(gè)標(biāo)記點(diǎn)來表示一個(gè)肺部圖形,如果采樣點(diǎn)數(shù)目進(jìn)一步增加,模型會(huì)有表達(dá)能力上的提高,后續(xù)工作是增加標(biāo)記點(diǎn)個(gè)數(shù),提高模型的表達(dá)能力.
參考文獻(xiàn):
[1] COOTES T F,TAYLOR C J. Statistical models of appearance
for medical image analysis and computer vision [C]//Medical Imaging 2001. [s.n.]:International Society for Optics and Photonics,2001:236-248.
[2] PARK H,MEYER C R,HYUNLIN PARK,et al. Construction of an abdominal probabilistic atlas and its application in segmentation [J]. IEEE Trans on Medical Imaging,2003,22 (4):483-492.
[3] SHEFFER A,PRAUN E,ROSE K. Mesh parameterization methods and their applications [J]. Foundations and Trends in Computer Graphics and Vision,2006,2(2):105-171.
[4] BAUMBERG A,HOGG D. Learning flexible models from image sequences[M]. Berlin:Springer Berlin Heidelberg,1994.
[5] WANG Y,PETERSON B S,STAIB L H. Shape-based 3D surface correspondence using geodesics and local geometry[C]// Computer Vision and Pattern Recognition,2000. Proceedings. [s.n.]:IEEE Conference on. IEEE,2000,2:644-651.
[6] BELONGIE S,MALIK J,PUZICHAJ. Shape matching and object recognition using shape contexts[J]. IEEE Trans on Pattern Analysis and Machine Intelligence,2002,42(4):509-522.
[7] SCOTT G L,LONGUET-HIGGINS H C. An algorithm for associating the features of two images[J]. Proc of the Royal Society of London,1991,244(1309):21-26.
[8] KELEMEN A,SZEKELY G,GERIG G. Elastic model-based segmentation of 3-Dneuroradiological data sets[J]. IEEE Trans. on Medical Imaging,1999,18(10):828-839.
[9] MEIER D,F(xiàn)ISHER E. Parameter space warping:Shape -based correspondence between morphologically different objects[J]. IEEE Trans on Medical Imaging,2002,21(1):31-47. [10] GU X,WANG Y,CHAN T F,et al. Genus zero surface conformal mapping and its application to brain surface mapping[J]. IEEE Trans on Medical Imaging,2004,23(8):949-958.
[11] GU X,WANG Y,YAU S T. Geometric compression using riemann surface structure [J]. Communications in Information and Systems,2004,3(3):171-182.
[12] FLOATER M S,HORMANN K. Surface Parameterization:A tutorial and survey [C]//Advances in Multiresolution for Geometric Modelling Mathematics and Visualization. 2005:157-186.
[13]夏述高.三角曲面參數(shù)化若干問題研究[D].大連:大連理工大學(xué),2011. XIA S G. Some issues triangle parametric surfaces[D]. Dalian:Dalian University of Technology,2011(in Chinese).
[14] POWELL M J D. An efficient method for finding the minimum of a function of several variables without calculating derivatives [J]. The Computer Journal,1964,7(2):155-162.
[15] DAVIES R H. Learning shape:Optimal models for analysing natural variability[D]. Manchester:Dissertation University of Manchester,2002.
Signature point matching method based on spherical conformal mapping theory
TANG Chun-ming,NIU Bai-ling,LI Guang-xu
(School of Electronics and Information Engineering,Tianjin Polytechnic University,Tianjin 300387,China)
Abstract:A semi-automatic matching method based on surface parameterization theory is presented. Firstly,one of surface sample data as a reference model is selected. According to the geometric characteristics of surface,some landmark points are choosed manually. Then,all the sample data are mapped to the same sphere by the conformal mapping method. The landmark points on the other surface samples are completed automatically according to the positions corresponding to their spherical coordinates. Finally,to normalize the positions of landmarks in threedimensional space,the Iterative Closed Points method is utilized. The experiments measure the general principle and specificity principle of the model. Moreover,the influence on the point distribution model due to the different constraints of spherical conformal mapping method is compared. The results show that the specificity of the model based on spherical confomcal mapping theory is much letter than that of zero centroid constraint method,the time for the experiment is decreased by 50%.
Key words:spherical conformal mapping theory;signature point matching;statistical shape model building;active shape
通信作者:湯春明(1971—),女,教授,主要研究方向?yàn)閳D像處理與模式識(shí)別. E-mail:tangchunminga@hotmail.com
收稿日期:2015-09-01基金項(xiàng)目:天津市應(yīng)用基礎(chǔ)與前沿技術(shù)研究計(jì)劃項(xiàng)目(14JCYBJC42300)
DOI:10.3969/j.issn.1671-024x.2016.01.011
中圖分類號(hào):TP311
文獻(xiàn)標(biāo)志碼:A
文章編號(hào):1671-024X(2016)01-0054-05
天津工業(yè)大學(xué)學(xué)報(bào)2016年1期