趙晶潔，趙瑞斌，龐明勇

?

保持尖銳特征的隱式曲線繪制算法

趙晶潔，趙瑞斌，龐明勇

(南京師范大學(xué)教育信息工程研究所，江蘇 南京 210097)

隱式曲線在生物、醫(yī)學(xué)、氣象、地學(xué)、石油勘探及物探等領(lǐng)域有著廣泛的應(yīng)用。提出一種繪制帶有尖銳特征的平面隱式曲線的算法，能有效地提取隱式曲線的尖銳特征。該算法首先確定曲線的繪制區(qū)域，采用自上而下的方式生成繪制區(qū)域的四叉樹表示，并在四叉樹節(jié)點表示的每個單元格內(nèi)生成一個數(shù)值場特征點；然后連接特征點生成對偶網(wǎng)格；最后，利用Marching Squares算法生成曲線。實驗結(jié)果表明，該算法能在網(wǎng)格較稀松的情況下繪制出隱式曲線，并且可以實現(xiàn)曲線的尖銳特征。

隱式曲線；圖形繪制；可視化；移動四邊形

隱式曲線的繪制工作一般包括“采樣”和“構(gòu)造”兩個階段[10]。采樣階段的主要任務(wù)是運用特定策略或技術(shù)采集位于隱式曲線上的點；構(gòu)造階段的目標(biāo)則是將采集到的點以恰當(dāng)方式連接起來，構(gòu)造出分段線性的折線逼近原來的隱式曲線。目前，隱式曲線的繪制方法主要有：基于空間剖分的方法[11]、基于連續(xù)性的方法[12]、以及基于物理模擬的方法[13]等幾類。在基于空間剖分的方法中，最為經(jīng)典流行的是被稱為Marching Squares (MS)的方法[14]，其是可視化領(lǐng)域中著名的Marching Cube方法的二維特例。該方法利用了隱式曲線的一個良好性質(zhì)，即的隱式曲線是平面2上分別滿足()>0和()<0的2個點集的邊界。對于平面上的任一線段，若其2端點處的隱函數(shù)值異號，則該線段上至少有的一個零點，可采用二分迭代法或插值法逼近該零點。在繪制隱式曲線的過程中，①在繪制區(qū)域上構(gòu)造一個四邊形格網(wǎng)；②計算各個格網(wǎng)頂點處的隱函數(shù)值，并確定每個格網(wǎng)邊上可能存在的零點；③用線段連接格網(wǎng)中每個四邊形的不同邊上的零點，且逼近隱式曲線在該四邊形內(nèi)的部分，從而在繪制區(qū)域中繪制出隱式曲線。為了實現(xiàn)隱式曲線的精細(xì)化繪制，往往需要構(gòu)造高密度的格網(wǎng)，這會帶來巨大計算量[15]。為了降低計算量，學(xué)者們提出了空間自適應(yīng)剖分方法，即先在繪制區(qū)域中構(gòu)造低密度格網(wǎng)，再對各個包含曲線的四邊形進(jìn)行自適應(yīng)細(xì)分，并用四叉樹數(shù)據(jù)結(jié)構(gòu)存儲平面區(qū)域的細(xì)分?jǐn)?shù)據(jù)，從而消除部分冗余計算，從而降低了繪制過程中的總計算量?；谶B續(xù)性的方法一般利用空間和曲線的連續(xù)性特性，通過分析隱式曲線的切向、曲率等局部幾何性質(zhì)來實現(xiàn)曲線的繪制[16]。如文獻(xiàn)[17]的算法首先找到隱式曲線上的一個點，然后以該點為中心構(gòu)造一個正方形，再根據(jù)隱式曲線與該正方形四條邊的相交情況，不斷向四周擴展出新的四條邊，最后通過對隱式曲線在各四邊形內(nèi)的部分進(jìn)行線性化處理，構(gòu)造出整條隱式曲線(圖1(a))。文獻(xiàn)[18]則從曲線上的一點出發(fā)，根據(jù)曲線的局部曲率大小沿曲線在該點處的切向前進(jìn)一定的距離，此過程稱為預(yù)測；再根據(jù)函數(shù)定義的梯度方向，將預(yù)測的新點以投射方式“更正”到曲線上，然后從新點出發(fā)，重復(fù)“預(yù)測-更正”過程，從而繪制出整條曲線(圖1(b))?；谖锢砟M的方法將隱函數(shù)視為二維空間的數(shù)值場函數(shù)，然后采用模擬某種物理過程來驅(qū)動采樣點運動，最終實現(xiàn)曲線的優(yōu)化繪制。如文獻(xiàn)[19]用一條封閉折線包圍隱式曲線，然后通過場函數(shù)移動折線上的點，使之迭代“收縮”到隱式曲線上。文獻(xiàn)[20]提出采用粒子系統(tǒng)在曲線上采樣，所得到的采樣點在曲線上的分布是優(yōu)化的。

圖1 2種基于連續(xù)性的隱式曲線繪制方法

當(dāng)隱式曲線是連續(xù)光滑的，或隱函數(shù)數(shù)值場平滑變化時，現(xiàn)有方法均能很好地工作。但當(dāng)隱式曲線上具有尖銳特征或數(shù)值場局部變化劇烈時，基于空間剖分的方法往往不能捕捉到相應(yīng)的尖銳特征，而基于連續(xù)性與基于物理模擬的方法通常不能魯棒地工作。另一方面，在醫(yī)學(xué)、地學(xué)等領(lǐng)域中，有大量的應(yīng)用涉及到帶有尖銳特征的隱式曲線的繪制需求。

針對上述問題，本文提出一種保持隱式曲線尖銳特征的繪制算法。算法以MS技術(shù)為基礎(chǔ)，采用自適應(yīng)剖分方法來減少冗余數(shù)據(jù)，其特別之處在于，通過引入數(shù)值最優(yōu)化技術(shù)等實現(xiàn)曲線上尖銳特征的捕捉，最終繪制出帶有尖銳特征的隱式曲線。與現(xiàn)有方法相比，本文算法不但能夠有效地保持隱式曲線上的尖銳特征，而且其在處理不同隱式曲線時，具有較好的魯棒性。

1 算法梗概

圖2 繪制區(qū)域上兩種不同四邊形格網(wǎng)的構(gòu)造方法

圖3 單元格內(nèi)特征點與采樣點的關(guān)系

步驟4. 生成對偶網(wǎng)格。對于提取到的特征點，根據(jù)其所在四叉樹單元與鄰近其他單元之間的毗鄰關(guān)系，使用文獻(xiàn)[21]提出的遍歷算法構(gòu)造自適應(yīng)細(xì)分四叉樹網(wǎng)格的對偶網(wǎng)格。其中，對偶網(wǎng)格的頂點為數(shù)值場的特征點。上述毗鄰關(guān)系決定了對偶網(wǎng)格頂點間的拓?fù)溥B接。

步驟5. 隱式曲線的逼近與繪制。對于對偶網(wǎng)格中的每個單元，基于MS方法對其內(nèi)部的隱式曲線進(jìn)行分段線性離散逼近，從而實現(xiàn)隱式曲線的繪制。

2 算法步驟

本文算法主要包含4個步驟，本節(jié)將分別對其進(jìn)行詳細(xì)討論。

2.1 構(gòu)造自適應(yīng)細(xì)分四叉樹平面剖分結(jié)構(gòu)

為了進(jìn)一步逼近隱式曲線并降低計算復(fù)雜度，在均勻四叉樹的基礎(chǔ)上，定義了如下細(xì)分規(guī)則并進(jìn)一步構(gòu)造自適應(yīng)細(xì)分的四叉樹結(jié)構(gòu)。

(1) 對于均勻四叉樹中的每個節(jié)點，將其4個頂點分別代入，得到4個隱函數(shù)值。如果這4個值均為正或均為負(fù)，則認(rèn)為該節(jié)點具有一致的頂點符號，同時隱式曲線不經(jīng)過該節(jié)點且無需對其繼續(xù)細(xì)分。

2.2 提取隱函數(shù)數(shù)值場特征點

為了逼近隱式曲線并突出其尖銳特征，本文算法將視為二元數(shù)值場函數(shù)，并借助計算數(shù)值場的特征來探測隱式曲線的特征。算法在已構(gòu)造自適應(yīng)細(xì)分四叉樹的基礎(chǔ)上，通過對四叉樹的每個葉子節(jié)點進(jìn)行數(shù)值場特征探測生成一個對應(yīng)的特征點。

本文隱式曲線特征點的提取，基于如下發(fā)現(xiàn)：給定一個四叉樹葉子節(jié)點，對于經(jīng)過該節(jié)點中的采樣點且以在該采樣點處的梯度方向為法向的每條直線，其在該節(jié)點內(nèi)特征點處的值與在特征點處的值之間存在誤差，把每個誤差疊加，所得的誤差和最小。因此，可根據(jù)上述發(fā)現(xiàn)構(gòu)造一個二次誤差函數(shù)(quadratic error functions，QEF)來求取每個葉子節(jié)點中的特征點。QEF定義如下

其中，是采樣點數(shù)量；T()是經(jīng)過采樣點且以在該采樣點處的梯度方向為法向的直線在特征點處的值，其可表示為

其中，為采樣點；T()為在采樣點p處的梯度；(,)為采樣點與特征點之間的距離。

為了更直觀地說明上述原理，圖3給出了一個一維情況下的特例。其中，曲線為數(shù)值場；1和2為2個采樣點；虛線1，2分別為經(jīng)過1和2，且以隱函數(shù)在1和2處的梯度方向為法線的直線；如果在軸上選擇2個點3和4，可以發(fā)現(xiàn)，在3點處1(3)與(3)之間的差值為0，而在4點處的差值為(4)–(4)。如果在軸上選取除3外任意一點，該點處的差值都會大于0，在3點處，差值最小。由圖3可知，場曲線在3點處呈現(xiàn)出尖角特征，即3是該數(shù)值場的特征點。因此，可通過最小化QEF找到場特征點。

本文算法通過如下步驟提取特征點：

步驟1. 對于自適應(yīng)細(xì)分四叉樹中的每個葉子節(jié)點，將其4個頂點看作采樣點；

步驟2. 將采樣點代入QEF函數(shù)，定義每個采樣點的T()與()之間的誤差平方和；

步驟3. 通過最小化值確定特征點的位置。

2.3 創(chuàng)建對偶網(wǎng)格

在構(gòu)建隱式曲線時，傳統(tǒng)的MS算法通過連接在網(wǎng)格中的零點來逼近隱式曲線在該網(wǎng)格單元內(nèi)的部分。然而，若曲線在該單元中有尖銳特征，則該方法往往會丟失這些特征。為此，本文算法通過連接自適應(yīng)細(xì)分四叉樹中的特征點來創(chuàng)建對偶網(wǎng)格。其能夠逼近隱式曲線的形狀并有效保留曲線的尖銳特征?；谝韵?個規(guī)則來創(chuàng)建對偶網(wǎng)格：

(1) 若自適應(yīng)細(xì)分四叉樹中，一個節(jié)點的4個孩子節(jié)點均為葉子節(jié)點，則將該4個孩子節(jié)點中的特征點依次相連，生成一個單元格，參見圖4(a)；

(2) 對于兩個相鄰節(jié)點，如果一個是葉子節(jié)點，一個是非葉子節(jié)點，算法首先查找到非葉子節(jié)點中的所有葉子節(jié)點，然后將最初的葉子節(jié)點與相鄰的新查找到的葉子節(jié)點的特征點依次相連組成一個單元格，參見圖4(b)；

(3) 如圖4(c)所示，對于2個具有相同邊界的非葉子節(jié)點，算法首先分別遍歷并找到其所對應(yīng)的所有葉子節(jié)點，然后將相同邊界周圍的葉子節(jié)點依次相連組成新單元格。

根據(jù)上述3個規(guī)則，通過對自適應(yīng)細(xì)分四叉樹每個節(jié)點內(nèi)的特征點依次相連，構(gòu)造出一個對偶網(wǎng)格，如圖5所示。

圖4 創(chuàng)建對偶網(wǎng)格的3個規(guī)則

圖5 對偶網(wǎng)格示意圖

2.4 繪制曲線

在構(gòu)造對偶網(wǎng)格的基礎(chǔ)上，本文算法借助經(jīng)典的MS算法實現(xiàn)隱式曲線的繪制。由于在繪制曲線時MS算法通常在形狀規(guī)則的四邊形上操作，而本文算法所創(chuàng)建的對偶網(wǎng)格中，既有不規(guī)則的四邊形又有三角形。因此，算法將對偶網(wǎng)格中的每個單元格在拓?fù)渖系韧谒倪呅?，再運用MS算法繪制隱式曲線。本文算法首先判斷對偶網(wǎng)格單元格的每條邊上是否存在隱函數(shù)的零點，然后對于有零點的邊則通過線性插值確定零點位置，最后依次連接單元格上的零點逼近隱式曲線在該單元格內(nèi)的部分，從而繪制出隱式曲線。

對于MS中可能出現(xiàn)的二義性問題，本文算法對有二義性的單元格進(jìn)行三角分割，求出單元格對角線上的中點坐標(biāo)以及其在該中點的函數(shù)值，該點與其他4個頂點構(gòu)成4個三角形。因為在三角形中總有1個頂點的符號與其他2個不同，所以可以避免二義性，如圖6所示。

圖6 MS中二義性情況處理示意

由于在對偶網(wǎng)格的創(chuàng)建過程中，對偶網(wǎng)格的頂點位于函數(shù)的數(shù)值場特征上，所以本文算法所繪制的隱式曲線能有效地保留尖銳特征。

3 實驗與結(jié)果

為了驗證本文算法的魯棒性和準(zhǔn)確性，本文使用Dev C++和OpenGL圖形庫在Windows平臺上實現(xiàn)了本文算法。并采用如下5個隱式曲線進(jìn)行曲線繪制的測試，并將繪制結(jié)果與傳統(tǒng)的MS算法以及文獻(xiàn)[22]中的算法在二維圖形上的運行結(jié)果進(jìn)行了比較。

對比圖7中給出的3種方法的繪制結(jié)果，可以發(fā)現(xiàn)：①傳統(tǒng)的MS方法(下文統(tǒng)稱方法1)明顯會丟失了原始隱式曲線上的尖銳特征。如圖7所示，第1行中所繪制出的曲線均未能有效地捕捉到相應(yīng)的尖銳特征，曲線上的尖銳特征被鈍化處理；②文獻(xiàn)[22]中的方法(下文統(tǒng)稱方法2)只保留了原曲線的部分尖銳特征，且所保留的尖銳特征精確度不夠。如圖7所示，第2行中所繪制的曲線沒有準(zhǔn)確捕捉到相應(yīng)的尖銳特征，如該行第1列中的曲線理論上應(yīng)繪制出“V”形局部，但該曲線明顯丟失了尖角特征；③相比于方法2，本文算法能準(zhǔn)確捕捉到曲線的尖銳特征。圖7第3行中所繪制的曲線所示，所有曲線均保留了尖銳特征，例如該行第一列中的曲線呈現(xiàn)出尖角特征，逼近于“V”形，比較準(zhǔn)確地繪制出了隱式曲線。

表1 3種算法性能效果對比

曲線網(wǎng)格行列數(shù)/四叉樹深度線段數(shù)計算時間(s) MS方法文獻(xiàn)[22]方法本文方法MS方法文獻(xiàn)[22]方法本文方法MS方法文獻(xiàn)[22]方法本文方法 式(3)23881915150.018 8590.016 0930.016 465 式(4)26661 1545185400.492 5310.087 5550.089 490 式(5)26661 4195656000.544 8910.282 2330.314 777 式(6)255521388880.384 6640.345 3090.382 034 式(7)26661 0962673290.453 8320.374 1450.382 160

圖7 不同算法對各函數(shù)曲線的繪制結(jié)果

在性能方面，對比表1數(shù)據(jù)可知: 在網(wǎng)格深度相同的情況下，方法2與本文算法繪制曲線所用的線段數(shù)少于方法1，并且這2種方法運行速度均快于方法1。當(dāng)隱式曲線比較復(fù)雜時，本文算法需要耗費額外的時間來精確計算數(shù)值場的特征點，除此之外，利用這些特征點生成對偶網(wǎng)格也需要額外的時間，所以相對于方法2，本文算法生成曲線的速度較慢。

4 結(jié)束語

本文提出一種能保持尖銳特征的隱式曲線繪制算法。算法通過構(gòu)造QEF定位、捕捉隱函數(shù)數(shù)值場的特征點來捕捉隱式曲線的特征，通過連接特征點生成對偶網(wǎng)格來逼近了隱式曲線的形狀。實驗表明，本文算法能在較稀松的網(wǎng)格中捕捉隱函數(shù)數(shù)值場的特征點，并且能在由這些特征點構(gòu)建的對偶網(wǎng)格中提取出較為精確的曲線，提升了曲線的繪制精確性。

為了逼近隱式曲線并且保留曲線的尖銳特征，本文算法需要生成自適應(yīng)細(xì)分四叉樹以及對偶網(wǎng)格，因此算法的運行時間較長。若要解決此問題，可開發(fā)只在包含曲線的單元格內(nèi)提取特征點或者創(chuàng)建對偶網(wǎng)格的算法，預(yù)計該算法會提高運算速度。

[1] 徐國良. CAGD中的隱式曲線與曲面[J]. 數(shù)值計算與計算機應(yīng)用, 1997, 18(2): 114-124.

[2] BLOOMENTHAL J. An implicit surface polygonizer [M]. San Diego: Academic Press, 1994: 324-349.

[3] SHELBERG M C, MOELLERING H, LAM N. Measuring the fractal dimensions of empirical cartographic curves [M]. Rome: Food and Agriculture Organization of the United Nation, 1982: 481-490.

[4] 趙偉, 趙卓寧, 李五生. 一種有效的離散數(shù)據(jù)場等值線生成方法[J]. 成都信息工程學(xué)院學(xué)報, 2007, 22(1): 116-121.

[5] 溫維亮, 郭新宇, 陸聲鏈, 等. 隱式曲面在三維植物建模中的應(yīng)用研究綜述[J]. 圖學(xué)學(xué)報, 2010, 31(3): 1-10.

[6] SUNDARAMOORTHI G, YEZZI A, MENNUCCI A. Coarse-to-fine segmentation and tracking using Sobolev active contours [J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2008, 30(5): 851-864.

[7] 壽華好, 何蘋, 繆永偉. 自動微分在隱式曲線繪制中的應(yīng)用[C]//第四屆全國幾何設(shè)計與計算學(xué)術(shù)會議. 北京: 中國工業(yè)與應(yīng)用數(shù)學(xué)學(xué)會, 2009: 150-153.

[8] 李華, 蒙培生, 王乘. 醫(yī)學(xué)圖像重建MC算法三角片的合并與實現(xiàn)[J]. 計算機應(yīng)用, 2003, 23(6): 104-106.

[9] 李彩云, 朱春鋼, 王仁宏. 參數(shù)曲線的分段近似隱式化[J]. 高校應(yīng)用數(shù)學(xué)學(xué)報A輯, 2010, 25(2): 202-210.

[10] 駱沛, 吳壯志, 夏春和, 等. 基于e1范數(shù)最小化的非流形曲線族重構(gòu)[J]. 計算機學(xué)報, 2013, 36(9): 1917- 1928.

[11] LAGUARDIA J J, CUETO E, DOBLARE M. A natural neighbour Galerkin method with quadtree structure [J]. International Journal for Numerical Methods in Engineering, 2006, 15(5): 529-548.

[12] ZHU X J, JI L X, JIN X B. Fitting and reconstruction of three-dimensional curve based on orthogonal curvature [C]//2009 9th International Conference on Electronic Measurement and Instruments. New York: IEEE Press, 2009: 323-328.

[13] 李云夕, 馮結(jié)青, 金小剛. 基于有向距離場的代數(shù)B-樣條曲線重建[J]. 軟件學(xué)報, 2007, 18(9): 2306-2317.

[14] MAPLE C. Geometric designing and space planning using the marching squares and marching cube algorithms [C]//2003 International Conference on Geometric Modeling and Graphics. New York: IEEE Press, 2003: 90-95.

[15] MANTZ H, JACOBS K, MECKE K. Utilizing Minkowski functionals for image analysis: a marching square algorithm [J]. Journal of Statistical Mechanics Theory and Experiment, 2008, 12 (12): 12015.

[16] 肖海, 章亞男, 沈林勇, 等. 光纖光柵曲線重建算法中的曲率連續(xù)化研究[J]. 儀器儀表學(xué)報, 2016, 17(5): 993-999.

[17] CIPOLLETTI M P, DELRIEUX C A, PERILLO G M E, et al. Superresolution border segmentation and measurement in remote sensing images [J]. Computers and Geosciences, 2012, 40(3): 87-96.

[18] AKKOUCHE S, GALIN E. Adaptive implicit surface polygonization using marching triangles [J]. Computer Graphics Forum, 2001, 20(2): 67-80.

[19] 劉穎. 復(fù)雜的隱式函數(shù)曲線繪制算法的研究[D]. 長春: 吉林大學(xué), 2006.

[20] 范麗鵬, 王麗英, 龐明勇. 平面簡單閉合曲線離散采樣與重建算法[J]. 圖學(xué)學(xué)報, 2015, 36(4): 511-515.

[21] JU T, LOSASSO F, SCHAEFER S, et al. Dual contouring of Hermite data [J]. ACM Transactions on Graphics, 2002, 21(3): 339-346.

[22] RAMAN S, WENGER R. Quality isosurface mesh generation using an extended marching cubes look up table [J]. Computer Graphics Forum, 2010, 27(3): 791-798.

An Algorithm for Visualizing Implicit Curves with Sharp Features

ZHAO Jing-jie, ZHAO Rui-bin, PANG Ming-yong

(Institute of EduInfo Science and Engineering, Nanjing Normal University, Nanjing Jiangsu 210097, China)

Implicit curve plays an essential role in the fields of medicine, meteorology, geology, petroleum exploration, geophysics and so on. In this paper, we propose an algorithm to visualize implicit curves with sharp features, which can effectively extract the sharp features of such curves. The algorithm first defines the visualizing area of the curve and then adopts a quadtree that generates visualizing area by a top-down method. In each cell, the method produces a feature point of the numerical field, and connects different feature points to generate the dual mesh. Finally, the algorithm employs the Marching Squares algorithm to generate the curves. Experiments show that our method can efficiently extract the sharp features of implicit curves, and it can work with various implicit curves with or without sharp features robustly.

implicit curve; graphic plotting; visualization; marching squares

TP 391

10.11996/JG.j.2095-302X.2019020373

A

2095-302X(2019)02-0373-06

2018-06-21；

2018-09-12

全國教育科學(xué)“十三五”規(guī)劃教育部重點課題(DCA170302)；江蘇省社會科學(xué)基金項目(15TQB005)

趙晶潔(1994-)，女，山西晉城人，碩士研究生。主要研究方向為數(shù)字媒體技術(shù)。E-mail：Ginjer@yeah.net

龐明勇(1968-)，男，安徽淮南人，教授，博士，博士生導(dǎo)師。主要研究方向為幾何建模與數(shù)字幾何處理。E-mail：panion@netease.com