基于特征分布的三維流線相似性研究

張伯雄， 謝伙生

(福州大學數學與計算機科學學院， 福建 福州 350116)

基于特征分布的三維流線相似性研究

張伯雄， 謝伙生

(福州大學數學與計算機科學學院， 福建 福州 350116)

針對利用流線的形狀特征對流線進行分類和選取， 可以方便用戶洞察三維流場的特征, 提出一種高效的基于特征分布的三維流線相似性比較方法. 該方法在Lu Kewei等流線相似性比較方法的基礎上， 引入曲折度和速度方向熵兩個全局屬性. 首先將流線均勻分段， 然后計算每段的曲率直方圖、 扭率直方圖、 曲折度直方圖、 速度方向熵直方圖， 構建相應的二維直方圖， 最后利用堆土機距離(EMD)及k-means聚類方法進行流線相似度計算和分類. 實驗結果表明， 該方法在引入全局幾何屬性后能夠產生更魯棒性的查詢和分類結果.

流線相似性； 特征分布； 二維直方圖； 流線聚類

0 引言

流場可視化是科學計算可視化的一個重要課題, 在許多科學和工程領域， 可視化矢量場在理解潛在的流場特征和模式方面起著重要的作用. 流場可視化研究方法包括基于圖標的方法[1]、 基于紋理的方法[2]、 基于幾何圖形的方法[3]和基于特征結構的可視化方法[4]. 其中因為流線便于計算， 可在不同分辨率下以交互的速率進行渲染， 所以流線可視化方法是最普遍使用的方法之一. 但當成千上百條流線被渲染描繪一個流場時， 由于雜亂和遮擋問題， 用戶很難發(fā)覺潛在的流場特征， 而且用戶往往只想查看具有特定形狀特征的流線. 因此， 將流線基于形狀特征進行分類， 展示和用戶目標相關的流線是非常重要的.

流線相似性度量和分類工作中， 有許多是基于點和點之間的距離. Chen等[5]把流線間局部相似性度量的方法引入到流線生成中， 這種度量方法涉及流線形狀和方向的統計信息， 以及每對流線的歐式距離； Yu等[6]計算流線間點對的歐式距離， 使用平均最短點距離度量流線相似性， 最后結合單鏈接方法進行層次聚類. 但這些相似性度量方法不具有旋轉、 平移、 縮放不變性， 兩條相似的流線可能因為方向不同、 距離過遠、 長短不一而導致度量距離過大不能劃分為同一類. 魯大營等[7]提出了一種基于迭代最鄰近點(ICP)算法與k-means聚類方法的流線生成算法， 首先利用ICP算法實現流線間輪廓特征上的配準并根據幾何相似性進行排序， 然后利用k-means聚類算法對流線分組， 這種方法雖然具有旋轉、 平移不變性， 但計算量過大.

近年來， 基于流線特征分布的相似性度量方法越來越受到重視. Mc Loughlin等[8]計算流線每點的各種特征值， 構建1D直方圖， 并用χ2測試比較相似性， 但由于形狀不相近的流線也可能具有相似的1D直方圖， 文獻[8]的方法將不相似的流線歸為一類. Lu等[9]考慮特征的空間分布信息， 將流線分段， 構建2D直方圖， 解決了文獻[8]方法中由于空間信息丟失導致的流線相似性度量錯誤， 但他的相似度量方法只考慮流線的曲率和扭率兩種局部屬性， 當流線選取的比例增大時， 新增流線的整體形狀可能與目標流線有較大的偏差. Li等[10]受特征包想法的啟發(fā)， 選擇k個具有代表性的特征矢量作為詞匯， 在詞匯組合中考慮兩點的弧長信息從而構建空間敏感的特征包， 通過加權的曼哈頓距離進行流線間的相似性計算， 并考慮流線的曲折度和速度方向熵兩個全局屬性， 使得查詢的流線結果與目標流線在整體上更加貼近. 但由于文獻[10]中加權曼哈頓距離不一定能正確表達人對流線相似性的理解， 在流線聚類實驗中， 可能使一條螺旋上升的流線被錯誤地劃到與其形狀相差甚遠的流線組中. 相比之下， 文獻[9]中基于2D直方圖的方法仍能正確地聚類， 具有更好的魯棒性.

因此， 采用2D直方圖比較流線相似性， 保證流線相似性比較和分類的魯棒性， 并引入文獻[10]中所用的曲折度和速度方向熵兩個全局屬性， 使得流線相似性比較在全局屬性上也能更加接近目標流線， 產生更準確的相似性度量結果.

1 三維流線特征的選擇

選取合適的特征描述子是構建直方圖重要的一步. 為了選取形狀相類似的流線, 本研究選取能夠描述流線形狀的相關特征描述子. 首先選擇曲率和扭率， 因為它們分別表示了曲線的彎曲程度和脫離密切面的變化程度. 曲率和扭率屬于流線的局部幾何屬性， 為了能結合流線的全局幾何屬性進行比較， 接著引入曲折度和速度方向熵. 曲折度反映了曲線整體的彎曲程度； 速度方向熵反映了速度方向的均勻程度， 亦即流線的平坦程度.

1.1 曲率

曲率是正切矢量關于弧長的變化率. 曲線上某點p的曲率， 可以通過下式進行計算：

(1)

1.2 扭率

扭率的計算公式如下：

(2)

其中:det表示矩陣的行列式.

1.3 曲折度

曲折度是曲線長度和始末兩點最短距離之比. 曲線的長度通過組成流線的直線段的累計長度來近似， 始末兩點的最短距離采用歐幾里得距離進行計算.

對于流線上點p的曲折度， 先算出起點和該點間的弧長s(p)， 然后除以兩點間的歐幾里得距離， 其計算公式如下所示：

(3)

其中:p0是起始點位置.

1.4 速度方向熵

對于流線上每一點， 使用該點小范圍領域的N點來計算熵值. 為了計算熵值， 使用文獻[11]的算法把單位球體分成50份相同面積的區(qū)域， 并判斷每個矢量指向哪一塊. 對于流線上點p的速度方向熵可通過下面的公式進行計算：

(4)

其中:di是領域中樣本點落入i方向塊的頻率.

2 基于特征分布的流線相似性查詢與分類

2.1 流線特征的2D直方圖

流線特征分布可選1D直方圖與2D直方圖. 對于1D直方圖， 存在兩條形狀不相同的流線并可能有著相似的特征分布的問題. 例如， 圖1(d)曲線由圖1(a)曲線經過分割后拼接生成. 顯然， 它們有著相同的一維特征分布， 但是兩者的形狀并不相同. 兩條曲線的1D曲率直方圖分別如圖1(b)和圖1(e)所示.

本研究采用2D直方圖， 對于給定的流線， 將其分成M段， 每段具有相同數目的樣本點數， 接著對每一段建立一個1D直方圖， 從而構建2D直方圖. 由于攜帶了空間信息， 2D直方圖可彌補由于丟失空間信息而導致的流線相似性度量錯誤的不足.

圖1(c)和圖1(f)分別為曲線1(a)和曲線1(d)的2D直方圖； 在把曲線分成4等分后， 曲線各段的特征分布情況便體現在2D直方圖中.

2.2 基于k-means的三維流線分類

本文中k-means新的聚類中心計算公式為：

(5)

其中:C(j， t)為第j簇新中心的第t個特征直方圖；H(i， t)為流線i(或聚類中心i)的第t個特征直方圖；nj為第j簇所包含的流線下標集合.

流線(或聚類中心)i，j的EMD相似性距離度量公式為：

(6)

k-means聚類的目標函數為：

(7)

計算具體步驟如下：

Step1: 選擇k條流線作為初始聚類中心；

Step2: 使用公式(6)計算出流線到各簇類域間的距離， 并把流線劃分到距離最短的一類中；

Step3: 使用公式(5)更新各簇類中心；

Step4: 重復Step2和Step3直到目標函數(7)小于一定的誤差或者達到一定的迭代次數， 終止計算.

2.3 交互式流線相似性查詢與分類方法

首先通過密集播種， 生成大量的流線； 然后根據每條流線的弧長， 均勻地插值采樣相同數量的N個點； 接著計算流線上每點的特征值并建立2D直方圖； 在最后的用戶交互里， 根據用戶選擇的流線， 利用EMD計算每條流線與目標流線的直方圖距離， 根據指定的數量比例顯示最相似的流線組； 然后使用k-means聚類方法進行流線分類， 按用戶要求顯示感興趣的聚類結果. 在交互期間用戶可以調整流線段數和特征值等分數； 也可以采用不同屬性組合或者賦予不同權值進行相似性比較試驗.

3 實驗結果與分析

本研究使用兩組數據進行實驗. 第一組數據由Roger Crawfis的臺風模擬算法生成， 第二組數據是來自文獻[12]提供的五臨界點數據集. 下面對這兩組數據分別進行流線相似性查詢和流線聚類實驗.

3.1 運行時間測試

實驗采用CPU-GPU混合解決方案. 硬件配置為： Intel Core i5-3470 CPU， 頻率 3.2 GHz； 16 G內存； GT630顯卡. 使用CPU進行流線的生成， 采用CUDA技術對特征以及直方圖進行并行計算. 表1列出了特征值估計、 直方圖構建及相似性比較所用的時間.

表1 實現時間Tab.1 Implementation time

從表1中可以看出， 特征值估計和直方圖構建所用時間都非常短， 用戶可以調整流線劃分的段數M和特征值的劃分數N. 由于調整是對單個特征的直方圖進行， 所以在1～2 s內便可獲得調整效果反饋.

3.2 流線相似性查詢實驗

試驗中臺風數據集和五臨界點數據集分別生成600條和500條流線， 每條流線采樣1 000點. 流線均分成10段， 四種特征值均劃分為30塊. 速度方向熵中領域的取點數左右各取20點. 在選取完目標流線后， 加入各種屬性進行相似性比較. 實驗對比效果如圖2～4所示.

實驗顯示結果可以明顯地看出：

1) 圖2(b)中, 當臺風流線最相似流線增加到50%時， 新增流線雖然也是漩渦狀， 但已開始向外伸展， 而圖2(c)中在加入全局屬性后新增流線依然貼近目標流線.

2) 圖3～4中, 五臨界點數據流線在考慮全局屬性后整體也更向目標流線靠攏.

3) 圖4(d)中, 不加入扭率， 選取的流線更為平坦.

由此可見, 不同屬性體現流線的不同形狀特征， 簡單地將各特征直方圖通過EMD計算出距離后相加可能不足以突顯流線的形狀特征. 不同屬性間的組合， 賦予各特征直方圖相似性距離以不同的權值可以更好地體現流線形狀特征.

3.3 流線分類實驗

實驗中兩組數據均用k-means聚類方法， 在相同初始聚類中心的條件下聚成3類， 不同顏色的流線組表示不同的類簇.

圖5～6為聚類效果對比圖. 從結果中可以看出， 圖5(c)中紅色簇下端流線部分更接近于圖5(a)中藍色的流線簇. 考慮了全局屬性的聚類結果將紅色部分歸為藍色簇， 聚類后整體在形狀上更加貼合.

圖6為五臨界點數據流線聚類效果. 從結果可以看出， 加入全局屬性后黃色簇流線都較順直， 紅色簇流線中都是兩端螺旋. 從兩組數據的聚類結果可以看出， 加入全局屬性后簇內流線更加相似， 簇間差別更加明顯.

4 結論

本研究通過特征2D直方圖距離進行流線相似性比較， 克服了1D直方圖空間信息丟失的缺陷， 同時引入曲折度和速度方向熵兩個全局屬性， 使得流線相似性查詢和聚類結果都更具魯棒性， 整體形狀更加緊湊. 該方法采用CUDA技術對特征值和直方圖進行并行計算， 消耗的時間都非常短， 為用戶提供了實時的交互查詢體驗. 該方法中直方圖如何量化， EMD距離如何計算， 是影響相似性比較結果的關鍵. 在流線特征值分布極不均勻時， 如何以較小的塊數進行劃分， 如何設計EMD距離， 不僅影響相似性比較的計算量， 也影響其比較效果. 這些問題將是下一步的研究方向.

[1] KLASSEN R V， HARRINGTON S J. Shadowed hedgehogs: technique for visualizing 2D slices of 3D vector fields[C]//Proceedings of the 2nd Conference on Visualization′91. [S.l.]: IEEE Computer Society Press， 1991: 148-153.

[2] SUNDQUIST A. Dynamic line integral convolution for visualizing streamline evolution[J]. IEEE Transactions on Visualization and Computer Graphics， 2003， 9(3): 273-282.

[3] MA J， WANG C， SHENE C K. Coherent view-dependent streamline selection for importance-driven flow visualization[C]//Proceedings of SPIE - The International Society for Optical Engineering. San Francisco: SPIE， 2013， 8 654(5 755): 1-15.

[4] OZEN C A ， ROCKWELL D. Flow structure on a rotating plate[J]. Experiments in fluids， 2012， 52(1): 207-223.

[5] CHEN Y， COHEN J， KROLIK J. Similarity-guided streamline placement with error evaluation[J]. IEEE Transactions on Visualization and Computer Graphics， 2007， 13(6): 1 448-1 455 .

[6] YU H F， WANG C L， SHENE C K，etal. Hierarchical streamline bundles[J]. IEEE Transactions on Visualization and Computer Graphics， 2012， 18(8): 1 353-1 367.

[7] 魯大營， 朱登明， 王兆其. 三維流場的流線提取算法[J]. 計算機輔助設計與圖形學學報， 2013(5): 666-673.

[8] MCLOUGHLIN T， JONES M W， LARAMEE R S，etal. Similarity measures for enhancing interactive streamline seeding[J]. IEEE Transactions on Visualization and Computer Graphics， 2013， 19(8): 1 342-1 353.

[9] LU K， CHAUDHURI A， LEE T Y，etal. Exploring vector fields with distribution-based streamline analysis[C]//2013 IEEE Pacific Visualization Symposium. Sydney: IEEE， 2013: 257-264.

[10] LI Y F， WANG C L， SHENE C K. Streamline similarity analysis using bag-of-features[C]//Proceedings of SPIE - The International Society for Optical Engineering. Francisco: SPIE， 2013， 9017: 628-637.

[11] LEOPARDI P. A partition of the unit sphere into regions of equal area and small diameter[J]. Electronic Transactions on Numerical Analysis， 2006， 25(12): 309-327.

[12] YE X， KAO D， PANG A. Strategy for seeding 3D streamlines[C]//16th IEEE Visualization Conference VIS 2005. Minneapolis: IEEE Computer Society， 2005: 471-478.

(責任編輯： 林曉)

3D stremlines similarity analysis based on distribution of measurements

ZHANG Boxiong， XIE Huosheng

(College of Mathematics and Computer Science， Fuzhou University， Fuzhou， Fujian 350116， China)

Considering the feature of shape， classification and selection of streamlines advance the understanding of the features in 3D flow field.In this paper， we propose an efficient 3D streamlines similarity comparison method based on distribution of features. This method is based on the streamline similarity comparison method proposed by Lu Kewei et al. and introduces two global geometric properties: tortuosity and velocity direction entropy. At first， we divide each streamline into segments evenly， then we construct curvature histogram、 torsion histogram、 tortuosity histogram and velocity direction entropy histogram for each segment， and form corresponding 2D histograms.In the end， we use the Earth Mover’s Distance(EMD) and k-means clustering method to measure the similarity between streamlines and classify these streamlines. The final experiment showed that after combining the two global features， this method can produce more robust query and clustering results.

streamline similarity; feature distribution; 2D histogram; streamline clustering

10.7631/issn.1000-2243.2016.05.0633

1000-2243(2016)05-0633-06

2014-09-24

謝伙生(1964-)， 副教授， 主要從事智能圖形圖像處理、 數據挖掘、 大規(guī)模數據可視化、 機器學習等方面研究， xiehs@sina.com

福建省自然科學基金資助項目(2014J01229)

TP391.41

A