張若男， 王仁芳

（1. 上海海洋大學(xué)信息學(xué)院，上海 201306；2. 浙江萬里學(xué)院計算機與信息學(xué)院，浙江 寧波 315100）

基于Snake的點模型谷脊線提取與優(yōu)化

張若男1,2， 王仁芳2

（1. 上海海洋大學(xué)信息學(xué)院，上海 201306；2. 浙江萬里學(xué)院計算機與信息學(xué)院，浙江 寧波 315100）

基于主動輪廓模型（Snake模型），提出一種點模型的谷脊線提取與優(yōu)化方法。首先構(gòu)建點模型的局部隱式曲面，并求出采樣點的曲率值；然后通過求解主曲率極值點得到潛在谷脊點，依據(jù)特征點的主方向連接谷脊點生成谷脊折線段；最后利用主動輪廓模型對谷脊折線段進行優(yōu)化。實驗結(jié)果表明，算法是魯棒的且能夠生成光滑的谷脊線。

點模型；特征提??；谷；脊；Snake模型

隨著三維掃描設(shè)備和獲取技術(shù)的快速發(fā)展，點模型的數(shù)字幾何處理已經(jīng)成為數(shù)字幾何處理領(lǐng)域的研究熱點[1-3]。特征線是三維幾何模型的重要組成部分，對其外觀及準確表達具有重要的作用，因此特征線的提取在幾何分析、曲面重建與編輯以及數(shù)據(jù)分塊等方面有著非常廣泛地應(yīng)用。在數(shù)字幾何處理領(lǐng)域，網(wǎng)格模型特征提取的研究工作已取得一定的成果[4]，然而點模型特征提取的研究工作相對較少，主要原因是點模型缺乏拓撲連接關(guān)系，且在獲取過程中易受噪聲影響。因此，點模型特征線的提取仍是一項挑戰(zhàn)性工作。

現(xiàn)有的特征線提取方法大致分為3類：①基于曲面擬合的方法[5-11]。文獻[5]提出了一種基于魯棒的移動最小二乘法(robust moving least square, RMLS)的特征線提取方法，該方法首先通過RMLS算子收集候選特征點并對候選點分別進行 RMLS局部擬合，然后將候選特征點投影到最近的曲面交線上，對投影點進行光順處理，最后通過特征折線生長獲得特征線。該方法對特征點的識別主要取決于尖銳邊上的投影點，因此特征線的提取效果比較精確，但該方法基于RMLS擬合，時間代價較大。文獻[8]將 RMLS引入到點云谷脊特征線提取的研究中，實現(xiàn)了點云模型谷脊線的較準確提取，但由于沒有進行優(yōu)化，最終效果不夠光滑。針對網(wǎng)格模型，文獻[11]提出了一種基于局部擬合的谷脊線提取方法，該方法是采用緊支撐徑向函數(shù)(compactly supported radial basis functions, CS-RBFs)擬合三角網(wǎng)格曲面，并將網(wǎng)格點投影到擬合曲面上，計算曲率梯度和曲率張量，然后計算網(wǎng)格邊端點的最大最小主曲率，跟蹤檢測出所有谷脊特征點。該方法在網(wǎng)格模型上取得了良好的效果，但仍然有谷脊線連接性較弱的問題，并且該方法效率較低。②基于主成分分析的方法[12-15]。文獻[12]通過主成分分析法(principal component analysis, PCA)來計算采樣點的法向，并根據(jù)采樣點局部鄰域內(nèi)法向的變化對點模型進行分類，然后在不同分類的邊界點構(gòu)造最小生成樹，最后生成特征線。文獻[13]通過構(gòu)造黎曼樹建立點云的拓撲連接信息，利用局部鄰域的協(xié)方差矩陣的特征值計算某一點屬于特征點的可能，構(gòu)造最小生成樹來確定特征線。文獻[14]提出了點模型的特征檢測算法，該算法首先基于多尺度的協(xié)方差分析得到點模型的特征區(qū)域，然后利用最小生成樹算法將特征區(qū)域連接，最后對特征線進行非真實感繪制。盡管這些方法能夠提取點模型的特征線，但是對噪聲較為敏感，同時對細微特征的提取效果有限。③基于統(tǒng)計的方法[16-18]。文獻[16]提出一種基于高斯法向聚類的特征提取方法，首先在一點局部鄰域內(nèi)構(gòu)建當前點及其鄰域點組成的所有三角形集合，利用高斯法向聚類算法對三角形法向進行聚類，然后依據(jù)法向聚類的個數(shù)判別當前點是否為特征點，最后用特征折線生長算法連接得到特征線。由于該算法需要對采樣點鄰域進行三角網(wǎng)格化，從而增加了運算的復(fù)雜度，且存在跨越特征邊的三角形，由此降低了特征識別的準確性。

為了能魯棒地提取出光滑的特征線，本文提出一種基于主動輪廓模型的谷脊特征線提取方法。通過構(gòu)建點模型的局部隱式曲面，求出采樣點的曲率值；通過求解主曲率極值得到潛在谷脊點，并依據(jù)谷脊點的主方向連接谷脊點生成谷脊折線；利用主動輪廓模型對谷脊線進行優(yōu)化。

1 算法概述

本文算法思想是通過三次多項式擬合點模型的局部隱式曲面，根據(jù)最大、最小主曲率標記出潛在的谷脊特征點，然后采用啟發(fā)式搜索法分別對谷脊點進行連接，并采用Snake模型對谷脊線進行平滑優(yōu)化。設(shè)定的點模型為P={p1, p2,…, pN}, pi∈R3。算法步驟如下（圖1）：

步驟1. 計算主曲率。采用三次多項式擬合采樣點的局部鄰域，得到采樣點的局部隱式曲面表達式，從而計算出采樣點的主曲率。

步驟2. 提取谷脊特征點。根據(jù)谷脊點定義標記出谷脊特征點。

步驟3. 生成谷脊折線。根據(jù)最小主曲率對應(yīng)的主方向，采用啟發(fā)式搜索的方法對谷脊點進行連接，生成初始谷脊折線。

步驟4. 谷脊折線優(yōu)化。對于生成的每一條折線，將其作為初始的主動輪廓線，定義能量函數(shù)并求其最小值，以優(yōu)化谷脊折線得到光滑的谷脊線。

圖1 谷脊線提取與優(yōu)化算法流程

2 谷脊點的提取

根據(jù)擬合的采樣點的局部隱式曲面，計算出采樣點的主曲率值，并通過設(shè)置曲率閾值篩選出潛在的谷脊特征點。

2.1 谷脊點的定義

谷點和脊點是曲面局部鄰域內(nèi)主曲率沿主方向變化的極值點，能夠很好地表征三維物體的幾何形狀。谷點為待重建曲面凹部的極值點，脊點為待重建曲面凸部的極值點。在三維曲面中，曲面上某一點處的主曲率衡量了曲面在該點的不同方向的彎曲程度。因此，待重建曲面上的離散點在不同方向有多個不同的曲率值，最大的曲率值稱為最大主曲率記為kmax，取得最大主曲率的方向是最大主方向，記為 tmax；最小的曲率值稱為最小主曲率記為kmin，取得最小主曲率的方向稱為最小主方向，記為 tmin；數(shù)學(xué)上可以證明，最大主方向與最小主方向是相互垂直的。谷脊點的定義[19]為：谷點是在tmin方向上局部主曲率變化取得極小值且 ki<0的點，即：

脊點是在tmax方向上局部主曲率變化取得極大值且ki>0的點，即：

2.2 谷脊點提取

根據(jù)谷脊點的定義可知，要提取谷脊點需要求得采樣點的主曲率，為了更準確地計算采樣點主曲率ki，對采樣點進行局部擬合。通過已有研究工作對點模型進行基于Mean Shift的幾何特征性聚類[20]，得到多個聚類單元，然后采用三次多項式對每個聚類單元進行局部曲面擬合。設(shè)c為聚類單元中心，在c處建立局部坐標系(c; u, v, w)，c為原點，(u, v)平面是與c點處法向正交的平面，w為c點處法向的方向。局部待擬合曲面g(x)設(shè)為：

擬合曲面與原始曲面的誤差度量為：

式(1)、(2)中，x=(u, v, w)為點pi在局部坐標系中的坐標表示。其中的未知系數(shù)通過對下式求最小值得出。

其中，Ps為當前聚類單元，ωi是類高斯函數(shù)，用以控制鄰域點到中心點的權(quán)重。聚類單元中的點離中心點越近，對擬合精度的影響越大，其權(quán)重較大；相反，聚類單元中離中心點較遠的點，對擬合精度的影響較小，其權(quán)重較小。

根據(jù)局部曲面的隱式表達式g(x)，計算曲面的第一、第二基本量進而求得pi點的高斯曲率、平均曲率、最大主曲率和最小主曲率。圖 2給出了Fandisk點模型各曲率的可視化效果，最大、最小主曲率有效地識別出谷脊特征屬性。設(shè)定曲率閾值ε，當kmax>ε時，將pi標記為脊點，加入脊點集 ；當kmin＜-ε時，將pi標記為谷點，加入谷點集 。如圖2所示，粉色表示脊點區(qū)域，深藍色表示谷點區(qū)域。

圖2 使用平均曲率、高斯曲率、最大、最小主曲率標記Fandisk模型

3 谷脊線生成與優(yōu)化

3.1 谷脊線的生成

最大主曲率對應(yīng)的主方向是采樣點曲率變化最大的方向，即凹凸變化最大的方向；最小主曲率對應(yīng)的主方向是采樣點曲率變化最小的方向，即相對平緩的方向（即現(xiàn)實模型中山脊或者山谷對應(yīng)的方向），因此，本文通過在采樣點pi的k-近鄰域Nk( pi)內(nèi)根據(jù)最小主曲率對應(yīng)的主方向tmin將谷脊點進行連接，以避免采用最小生成樹構(gòu)造特征線的復(fù)雜過程。由于谷點和脊點可以通過變換z軸的方向相互轉(zhuǎn)換，因此，以脊點為例描述連接方法，谷點連接不再贅述。

采用啟發(fā)式搜索[21]的思想對脊點進行連接，步驟如下：

步驟2. 若Nk(ri)內(nèi)不存在脊點，那么特征線連接結(jié)束；返回步驟1重新選取未連接過的脊點進行連接，直到所有脊點都被連接。

步驟3. 遍歷谷脊特征線的每個節(jié)點，如果一個節(jié)點連接的兩條短線段之間的夾角，小于某一角度，那么就消除這個節(jié)點，將其他兩個點直接連接成線。

圖3 脊點的連接

3.2 谷脊線的優(yōu)化

在3.1節(jié)所得到的初始谷脊線是谷脊點直接首尾連接的結(jié)果，呈現(xiàn)為鋸齒狀的折線。在三維模型分割應(yīng)用中，折線是可以滿足分割要求的，但在特征線提取中，折線會造成后期渲染質(zhì)量低劣，可視化效果不好，因此需要一種機制來對初始谷脊特征線進行平滑優(yōu)化。常用的特征線優(yōu)化方法有B樣條擬合法[22]，然而該方法缺乏靈活性，且對于特征線平滑度和精確度地控制不好。Kass等[23]首次提出了主動輪廓模型(active contour model, ACM)，又稱為Snake模型，用于檢測圖像的特征。

Snake模型是一種目標輪廓提取模型，其本質(zhì)是一條二維可形變的彈性曲線。Snake模型的基本思想是在圖像目標特征邊界附近給出初始輪廓，該輪廓在自身內(nèi)力、圖像力等外力的共同作用下作變形運動，通過能量最小化使Snake收斂，最終逼近圖像目標特征邊界。

原始的Snake模型定義為在s∈[0, 1]上的參數(shù)曲線，即v(s)=[x(s), y(s)]，其中x(s)和y(s)分別表示每個控制點在圖像中的坐標位置，s是以傅立葉變換形式描述邊界的自變量。與模型相關(guān)的Snake能量函數(shù)定義如下：

內(nèi)部能量函數(shù)定義為：

外部能量函數(shù)定義為：

Snake模型廣泛應(yīng)用于二維圖像的輪廓提取和圖像分割，文獻[14]將其引入到三維模型特征線的提取中，取得了較好的效果，本文將Snake模型引入到谷脊線的優(yōu)化中。

根據(jù)Snake模型的相關(guān)知識，結(jié)合本文谷脊特征線優(yōu)化的目的，對 Snake模型進行分析。Snake模型有許多其他優(yōu)化方法無法企及的特性：

(1) 將特征、初始輪廓、目標輪廓和基于知識的約束都融合在同一過程中；

(2) 初始輪廓確定后，可以自主地移動直至收斂至能量最小狀態(tài)；

(3) 能量極小化可以極大地擴展捕獲區(qū)域，降低計算復(fù)雜度。

這些好的特性在三維模型谷脊特征線優(yōu)化的應(yīng)用中發(fā)揮著重要的作用，然而Snake模型也有其自身的缺點：

(1) 對初始位置要求很高，需要確保初始輪廓線在特征區(qū)域附近；

(2) 由于主動輪廓線的非凸性，收斂的過程可能會存在發(fā)散的可能。

基于以上分析，3.1節(jié)中所得初始谷脊線在特征區(qū)域附近，因此將谷脊折線作為初始輪廓線完全符合Snake模型對初始位置的要求。

將脊（或谷）折線作為初始輪廓線，它在內(nèi)力與外力的共同作用下存在變形能，使其偏向于能量最小的位置移動，經(jīng)過多次移動收斂得到光滑的谷脊線。內(nèi)部能量函數(shù)由拉伸能和彎曲能組成，在pi點的內(nèi)部能量可以表示為：

其中，Et(pi)是在pi點處的拉伸能，Eb(pi)是在pi點處的彎曲能，分別為：

其中，θi為向量 pi-1pi和向量pipi+1之間的夾角，α和β為控制參數(shù)，用于將能量分量歸一化到[0,1]范圍內(nèi)。

當初始輪廓線位于目標輪廓線附近時，外部能量通常取為點云在該點處的特征能。本文中連接得到的谷脊折線顯然位于目標谷脊線附近，因此，外部能Eext(pi)為該點處的特征能 Ef(pi)，其主要相關(guān)特征屬性為曲率值，因此定義如下：

其中，k( pi)表示點pi處的最大主曲率，γ用于將外部能量歸一化。

在Snake迭代過程中，谷脊折線的節(jié)點都會移動到局部鄰域內(nèi)能量最小的位置。單次迭代的實現(xiàn)算法如下，其中 resample( pi)是重采樣[24]函數(shù)，m是重采樣點的個數(shù)。

圖 4給出了脊線平滑的局部過程和局部效果圖，圖4(a)對相鄰3個脊點進行分析，將中間脊點調(diào)整至局部重采樣點中能量最小的點處，圖4(b)顯示了一段脊線平滑前后的效果比較。

圖4 脊線的平滑

4 實驗結(jié)果與分析

在VS2010環(huán)境下采用C++和OpenGL實現(xiàn)了本文的谷脊特征線提取算法，硬件平臺為1.8 GHz奔騰處理器、2 GB內(nèi)存的PC機。

圖5為本文算法得到的一系列點模型的谷脊線效果圖，可以看出，本文方法能很好地計算出點模型上的關(guān)鍵谷脊特征線條，可以很好地把握模型的主要幾何形狀。

為測試本文算法的魯棒性，本文對點模型加入不同程度的噪聲進行實驗結(jié)果比較。加入噪聲的方法是使點云中的點偏離原始位置隨機值，隨機值分布服從高斯分布，通過改變μ的值控制加入噪聲的程度。圖6中μ分別取0.00、0.01、0.05時，本文算法對Vase點模型的谷脊線提取效果，從圖6可以看出，本文算法是魯棒的，在有噪聲的情況下仍然可以較好地提取出谷脊特征線，其主要原因是本文采用抗噪的Mean Shift進行聚類，并進行三次隱式曲面擬合，在處理含有大量噪聲的模型時仍然可以保持良好的穩(wěn)定性。

圖7為文獻[9]和本文方法的Dolphin點模型和RockArm點模型的谷脊線提取效果對比，文獻[9]中首先對點模型進行空間分割，然后對點模型進行全局擬合，求得主曲率值進行谷脊特征的提取；本文方法首先對點模型進行Mean Shift聚類，對聚類單元進行局部擬合，求得主曲率值進行谷脊特征提取。由圖7中Dolphin點模型的身體部分(圖7(c)、(d)中放大區(qū)域)和RockArm點模型的型號烙印部分（圖7(g)、(h)中放大部分）可以看出，本文方法的谷脊線提取效果更加光滑，其主要原因是本文采用表面特征相似性的Mean Shift聚類作為局部曲面擬合域，從而提高了擬合質(zhì)量，所求曲率值更加準確。

圖8出示了Venus點模型和Dragon點模型的文獻[8]方法與本文方法的效果對比，由圖8(b)、(c)中放大區(qū)域可以看出，本文算法對特征更加敏感，谷脊線更加光滑；同時通過對比圖8(d)、(e)，很顯然本文算法得到的 Dragon點模型頭部的谷脊線更加光滑，連續(xù)性更好。其主要原因是本文對點模型直接進行三次曲面擬合，并采用Snake模型對谷脊線進行優(yōu)化；而文獻[8]是對點模型進行局部平面擬合后投影到曲面上，該過程會產(chǎn)生較大誤差，影響曲率值的準確性，進而影響谷脊特征判斷。

圖5 Buddha、Armadillo、Lasso、Dragon點模型的谷脊線效果

圖6 不同噪聲影響下Vase點模型谷脊線效果

圖7 文獻[9]方法與本文方法比較

圖8 文獻[8]方法與本文方法效果比較

5 結(jié) 束 語

本文提出了一種基于Snake模型的點模型谷脊線提取與優(yōu)化方法。首先對點模型進行了表面特征相似性的Mean Shift聚類，并在聚類單元上進行曲面擬合，以求得準確的曲率值。為了使得提取的谷脊線更加光滑，采用基于Snake模型的迭代方法對谷脊線進行優(yōu)化處理，通過定義能量函數(shù)并求解能量函數(shù)的最小值，使谷脊線收斂至局部能量最小處，得到光滑谷脊線。實驗結(jié)果表明，本文算法是魯棒的，對特征有較強的敏感性。

從實驗結(jié)果可知，本文算法能提取出較為光滑的谷脊特征線，但仍然存在谷脊線之間有裂縫的情況。為解決這個問題，未來工作需要對谷脊線的連接進行更加深入的研究，同時Snake模型迭代優(yōu)化較為耗時，尋求一種高效穩(wěn)定的優(yōu)化算法也是研究方向之一。

[1] 王仁芳, 李繼芳, 楊 慶, 等. 點模型的快速高質(zhì)量繪制[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2010, 22(2): 191-197.

[2] 王仁芳, 張三元. 數(shù)字幾何處理的若干問題研究進展[M].北京: 清華大學(xué)出版社, 2012: 33-34.

[3] 朱 建, 楊 龍, 劉維星, 等. 顯著性特征保持的點云模型縮放[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2014, 26(10): 1624-1632.

[4] 王愛霖, 劉 弘, 張桂娟. 基于谷脊線特征的三維網(wǎng)格模型簡化方法[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2014, 26(5): 788-793.

[5] Daniels II J, Ha L K, Ochotta T, et al. Robust smooth feature extraction from point clouds [C]//Proceedings of IEEE International Conference on Shape Modeling and Applications Los Alamitos. Lyon, France, 2007: 1123-1136.

[6] Kim S K. Extraction of ridge and valley lines from unorganized points [J]. Multimedia Tools and Applications, 2012, 63(1): 265-279.

[7] Weber C, Hahmann S, Hagen H, et al. Sharp feature preserving MLS surface reconstruction based on local feature line approximations [J]. Graphical Models, 2012, 74(6): 335-345.

[8] 龐旭芳, 龐明勇, 肖春霞. 點云模型谷脊特征的提取與增強算法[J]. 自動化學(xué)報, 2010, 36(8): 1073-1083.

[9] Ohtake Y, Belyaev A, Alexa M. Sparse low-degree implicit surfaces with application to high quality rendering, feature extraction, and smoothing [C]// Proceedings of Eurographics Symposium on Geometry. Vienna, Austria, 2005: 148-158.

[10] Ohtake Y, Belyaev A, Seidel H P. 3D scattered data approximation with adaptive compactly supported radial basis functions [C]//Shaping Modeling International. Riken, Japan, 2004: 31-39.

[11] Ohtake Y, Belyaev A, Seidel H P. Ridge-valley lines on meshes via implicit surface fitting [C]//Proceedings of ACM SIGGRAPH. Los Angeles, Califonia, USA, 2004: 6, 8. [12] Demarsin K, Vanderstraeten D, Volodine T, et al. Detection of closed sharp feature lines in point clouds for reverse engineering applications [C]//Report TW458 Department of Computer Science. Pittsburgh, PA, USA, 2006: 571-577.

[13] Gumhold S, Wang Xinlong, MacLeod R. Feature extraction from point cloud [C]//Proceedings of 10th International Meshing Roundtable. Berlin, Germany, 2001: 293-305.

[14] Pauly M, Keiser R, Gross M. Multi-scale feature extraction on point-sampled surfaces [J]. Computer Graphics Forum, 2003, 22(3): 281-289.

[15] Park M K, Lee S J, Lee K H. Multi-scale tensor voting for feature extraction from unstructured point clouds [J]. Graphical Models, 2012, 74(4): 197-208.

[16] Weber C, Hahmann S, Hagen H. Sharp feature detection in point clouds [C]//Proceedings of the Shape Modeling International Conference. Washington DC, USA, 2010: 175-186.

[17] Kalogerakis E, Simari P, Nowrouzezahrai D, et al. Robust statistical estimation of curvature on discretized surfaces [C]// Proceedings of Symposium on Geometry. Switzerland, 2007: 13-22.

[18] 王小超, 劉秀平, 李寶軍, 等. 基于局部重建的點云特征點提取[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2013, 25(5): 659-665.

[19] Belyaev A G, Pasko A A, Kunii T L. Ridges and ravines on implicit surfaces [C]//Computer Graphics International. Washington, DC, USA, 1998: 530-535.

[20] 王仁芳, 徐慧霞, 陳仲委, 等. 點模型微分屬性的估算及其應(yīng)用[J]. 自動化學(xué)報, 2011, 37 (12): 1474-1482.

[21] 張長東, 戴 寧, 廖文和, 等. 基于啟發(fā)式搜索策略的牙齒生物特征線提取技術(shù)[J]. 中國機械工程學(xué)報, 2012, 23(13): 1567-1586.

[22] 徐 進. 帶約束的 B樣條曲線曲面延伸技術(shù)[J]. 圖學(xué)學(xué)報, 2013, 34(3): 36-42.

[23] Kass M, Witkin A, Terzopoulos D. Snakes: active contour models [J]. International Journal of Computer Vision, 1987, 1(4): 321-331.

[24] 左軍毅, 張怡哲, 梁 彥. 自適應(yīng)不完全重采樣粒子濾波器[J]. 自動化學(xué)報, 2012, 38(4): 647-652.

Extraction of Valley-Ridge Lines from Point-Sampled Geometry with Snake

Zhang Ruonan1,2, Wang Renfang2
(1. College of Information Science, Shanghai Ocean University, Shanghai 201306, China; 2. College of Computer Science and Information Technology, Zhejiang Wanli University, Ningbo Zhejiang 315100, China)

Based on the active contour model, a robust algorithm is proposed for extracting and optimizing valley-ridge lines from point-sampled geometry (PSG). First, local implicit surface of the PSG is reconstructed, and the curvature of every sampling point is computed. Potential valley-ridge points are then obtained by solving the extreme of the main curvature, and potential valley-ridge points is connected into valley-ridge polylines according to the main direction of the feature points. Finally, snake model is used to smooth the valley-ridge polyline. Experimental results show that the algorithm is robust and can generate smooth valley-ridge lines.

point-sampled geometry; feature extraction; valley; ridge; Snake model

TP 391

A

2095-302X(2015)05-0662-09

2015-04-17；定稿日期：2015-05-18

浙江省科技計劃資助項目(2012C21004)；浙江省大學(xué)生科技創(chuàng)新(新苗人才計劃)(2014R419028)；寧波市自然科學(xué)基金資助項目(2013A610111)

張若男(1990-)，女，山東博興人，碩士研究生。主要研究方向為數(shù)字幾何處理。E-mail：ruonanzh@sina.com

王仁芳(1974-)，男，河南南陽人，教授，博士。主要研究方向為數(shù)字幾何處理。E-mail：renfang_wang@126.com