改進的最小包圍球隨機增量算法

李世林， 李紅軍

(北京林業大學理學院，北京 100083)

改進的最小包圍球隨機增量算法

李世林， 李紅軍

(北京林業大學理學院，北京 100083)

三維空間中離散點集的最小包圍球，在碰撞檢測、計算幾何和模式識別等領域都有廣泛應用。為了更好地理解和構造最小包圍球算法，首先對最小包圍球的性質進行分析。然后，基于對隨機增量算法的分析，提出了構造較大初始包圍球和減少迭代過程中最小包圍球更新次數兩種策略。依據后一種策略提出的方法稱為隨機點組-重算最遠點算法。計算機隨機生成數據和現實三維模型采樣數據的多組實驗結果表明，隨機點組-重算最遠點算法相比于之前的經典算法能夠有效地提高時間效率。

最小包圍球；隨機增量算法；隨機點組-重算最遠點算法

對于三維空間中的離散點集P，求解其最小包圍球(minimum enclosing ball，MEB)以尋找一個半徑最小的球，其包含P中所有的點，球體有許多良好的性質[1]，因此在計算機圖形學和計算幾何中，與離散點集的凸包圍多面體[2]相比，包圍球常作為邊界，可以更快地進行近似的碰撞檢測、范圍界定和形狀分析。在空間數據庫中，MEB的求解可用于建立空間數據索引，以便提高查詢速度[3]。同時MEB問題在模式識別[4]、計算幾何[5]、機器學習[6]等領域也有廣泛應用。模型局部的包圍球可以用于模型裁剪[7]或者重采樣，因這些領域處理的數據集通常都是大規模的，所以提高MEB算法運行效率是很有必要的。

MEB問題可以追溯到19世紀中期，Sylvester[8]在1857年提出了最小包圍圓問題，并在1860年，給出了求解該問題的一種近似線性算法[9]。之后人們開始系統地研究該問題，比如最近的基于α-殼的最小包圍圓求解算法[10]，并研究從二維平面到三維空間再到更高維空間的算法；1982年，Megiddo[11]首次給出了一種理論上的線性時間算法，其時間復雜度是最小的，但算法本身的實現卻很復雜；1991年，Welzl[12]提出了一種簡單快捷的隨機算法，其時間復雜度是線性級的，并指出該算法易于推廣到高維空間的 MEB和最小包圍橢球(minimum enclosing ellipsoid)中。后來，Berg等[5]對該方法進行了進一步的改進，稱之為隨機增量算法(randomized incremental algorithm，RIA)，該算法是求解MEB問題有代表性的經典算法。分析該算法可看出，影響求解的時間效率主要因素是MEB更新的次數，因此提高算法效率的策略有2種：①構造較大的初始 MEB；②修改迭代過程中的選點規則。在給出這兩種策略對應的具體算法之前，為便于理解和構造算法，先概述MEB的性質。

1　離散點集的最小包圍球及其性質

定義1. 點集P的MEB定義為：

其中， o*為最優球心，r*為最優半徑，表示Euclidean范數。

2012年，文獻[13]給出一系列有關最小包圍圓的性質，這里對應給出三維空間中MEB的有關性質，這將有助于RIA的理解和實現。

性質1. 三維空間中，有限離散點集的MEB是唯一的。

此性質的證明見文獻[12]。

若P中只有一個點p1，則MEB是退化的，半徑為0，球心為點p1。

性質2. 非退化的MEB可以由2～4個邊界點定義。

證明. 點集P的MEB的構成有3種情況：

(1) 點集 P中距離最遠的兩個點(p1,p2)作為邊界點，p1,p2兩點之間距離的一半作為球半徑，兩點連線段的中點作為球心。

(2) 點集P中任意3個點求其最小包圍圓，其中最大的圓半徑為r，圓心為o，確定該圓的那 3個點(p1,p2,p3)，若點集P中的任何一點與o之間的距離都小于r，則點集 P的 MEB的球心為o，半徑為r，邊界點是p1,p2,p3。

(3) 若有4個邊界點，那么點集P的MEB為這4個點所構成的空間四面體的外接球。由空間幾何的知識可知，空間四面體的包圍球是唯一的，若邊界點多于4個，那么邊界點集中任意4個點所確定的MEB是等同的，都可以歸納為4個邊界點。

性質3. 若點集P的MEB是由p1,p2,p3,p44個邊界點確定，那么這4個點的MEB是其所構成的空間四面體的外接球。

引理1. 設點集P的子集的 MEB為Di：

(1) 若點pi+ 1∈ Di，則點集的 MEB為Di+1= Di；

(2) 若點pi+ 1? Di，則點pi+1必定是在點集 P′的MEB即Di+1的邊界上。

引理1是RIA的遞增思想，其證明可參考文獻[5]。

證明. 由性質 3的條件可知，任意 3個點所確定的MEB都不包含第4個點。設p1,p2,p3所確定的MEB是球o，p4?o，由引理1可知，點p4是在點集的MEB的邊界上。同理p1,p2,p3也在其邊界上，又因為四面體的外接球是唯一的，所以該4個點的MEB是這4個點所構成空間四面體的外接球 。

2　較大初始最小包圍球的隨機增量算法

Berg等[5]給出一種易于理解的RIA，該算法的思想與動態規劃問題中動態維護最優解有相似之處。根據文獻[5]中關于平面點集最小包圍圓 RIA的描述，容易給出三維空間MEB的RIA。

RIA還有很大的改進空間，初始點對的選取直接影響到整個算法的迭代次數。傳統RIA只是隨機選取了兩個點作為初始MEB的直徑，盡管可以從平均意義上獲得比較好的時間效率且能避免最差的情況，但是迭代次數卻很大。本文采用較遠兩點構造較大初始 MEB，以便減少后續迭代過程中MEB的更新次數。

步驟1. 構造較大初始MEB。在點集P中任取點p1，找出距離p1最遠的點pi，再找出距離pi最遠的點pj，以點pi,pj的連線作為直徑，構造初始MEB，記為D2。此時的初始MEB大小相對于RIA更接近最優的 MEB D的大小，迭代中重新構造MEB的次數會變得更少，所以理論上改進后的RIA效率比原始算法更高。

后續步驟采用經典RIA，迭代更新當前MEB。

步驟2. 設MEB D2的兩個點分別為p1,p2, 令i=3轉到步驟3。

步驟3. 檢驗點pi，如果pi在包圍球Di-1中，那么MEB不變，即Di=Di-1；否則需根據性質3重新構造MEB Di。

步驟4. 令i=i+1，若i≤n，則轉到步驟3，否則算法終止，輸出Dn的半徑和球心。

顯然，步驟1的時間復雜度為O(n)，步驟2為O(1)，步驟 3與步驟 4構成循環，時間復雜度為O(n)，因此該方法確定初始點對的時間復雜度也是O(n)。由于新的改進算法在步驟1選取較大的包圍球，所以步驟3中的MEB重新構造的次數會大大減少，因此新的算法能夠有效地提高時間效率，稱該算法為較大初始MEB的RIA，記為iRIA。

3　隨機點組-重算最遠點算法

通過修改RIA迭代過程中的選點規則，提出了一種新算法，稱為隨機點組-重算最遠點算法(random point group - recalculation the farthest point，R-RCF)。

3.1算法介紹

步驟1. 將P中各點打亂成一個隨機排列p1,p2,···,pn。

步驟2. 作前4個點p1, p2,p3,p4確定的MEB。根據性質2和性質3，球的邊界可能通過這4個點(4點的外接球)；也可能只通過其中的3個點(3點的外接圓構成球的最大截圓)，但包含第4點；還可能只通過其中的兩個點，后一種情況球邊界上的兩點一定是位于球的一條直徑的兩端。令i=5，轉到步驟3。

步驟3. 考察點pi，如果該點位于球外部，轉到步驟4。否則，令i=i+1，若i≤n，執行步驟3；若i＞n，執行步驟5。

步驟4. 根據引理1，在p1,p2,p3,p4中任意選擇3個點，加上pi，構造包含這5個點的MEB，且半徑是所有滿足條件所構造的MEB中最小的，此時確定該MEB的4個點記為新的 p1,p2,p3,p4，令i=i+1，若i≤n，執行步驟3；若i＞n，執行步驟5。

步驟5. 在點集P中找出距離MEB的球心最遠點F。若F點已在球內或球邊界上，則該球即為所求的球，算法結束。否則，執行步驟6。

步驟6. 以F為邊界點，在點集｛p1, p2, p3, p4｝中選擇3個點和點F一起構造新的MEB，使得這5個點都在球內，且半徑是所構造的MEB中最小的，此時選擇的3個點和點F記作新的 p1, p2,p3,p4，轉到步驟5。

R-RCF將求解MEB問題分為2部分進行：①動態更新最優解利用了RIA的核心思想，該算法的優點是每次迭代所需的計算量遠遠小于RIA，但缺點是每次更新最優解會遺漏少部分已經位于上一步的MEB內部的點。②重算遺漏部分則是利用了最遠點依次選取算法的思想[14]，且相對于最遠點依次選取算法，該部分迭代次數更少，同時又重算了動態更新中所遺漏的部分，彌補了動態更新最優解的缺點。本文稱該算法為R-RCF。

3.2算法分析

3.2.1算法重算部分的必要性

由于R-RCF算法在步驟4動態更新最優解的過程中會遺漏少部分已經包含在球內部的點，所以該算法加入了步驟5和步驟6，對步驟4動態更新最優解的過程中所遺漏的區域進行重算。為了更形象地說明該問題，以二維平面最小包圍圓為例進行分析。

如圖1所示，ci是點ABC所確定的最小包圍圓，當遍歷到點D時，由于點D不在ci里，所以需要更新最小包圍圓，ci+1是滿足算法步驟 4的圓，此時s區域內遍歷過的點已經不包含于更新后的最小包圍圓中，要想使其重新位于最小包圍圓中，就需要對該區域進行重算。該算法根據參考文獻[14]中算法的步驟3，采用重算最遠點的方法進行重算，同時該算法動態更新最優解的思想恰是參考RIA，且利用邊界點動態更新最優解。

圖1　R-RCF更新最優解時遺漏s區域的點

3.2.2算法收斂性

本節要證明上述 R-RCF算法一定能終止，且最后一次求得的MEB就是包含點集P中所有點的MEB。

引理2. 算法步驟4和步驟6所生成的MEB的半徑是遞增的。

證明. 設步驟4生成的包含A, B, C, D 4點的MEB為B1，半徑為 r1；當考慮球外的點E后，生成的包含A, B, C, D, E 5點的MEB記為B2，半徑為r2，顯然B2包含A, B, C, D 4點。考慮到B2的構造，會有B2≠B1。若r1=r2，與空間離散點集的 MEB的唯一性存在矛盾；若r1＞r2，與點A, B, C, D的MEB為 B1也存在矛盾；所以r1＜r2，即更新后的MEB半徑比原MEB的半徑大。

定理1. 上述R-RCF算法是收斂的，且最后一次得到的MEB是包含點集P中所有點的MEB。

證明. 因為在點集P中任取4個點、3個點或2個點所生成的MEB的個數是有限的。由引理2可知，算法進行過程中所生成的MEB半徑是遞增的，經過有限次迭代后，可求得這些最小球中半徑最大的一個。從算法步驟5可知，只有當點集P中所有點都在由4個點、3個點或2個點生成的MEB內部時，算法結束，因而最后得到的MEB會是包含點集P中所有點的MEB。

4　實驗驗證

上述算法的實驗是在筆記本電腦上進行的，計算機配置是 Intel(R) Core(TM) i3-2310M CPU @ 2.10 GHz處理器和 6 GB內存，程序運行環境MATLAB R2012a。實驗針對RIA，iRIA和R-RCF的時間效率進行對比，實驗數據分為兩類：計算機隨機生成點集(簡稱為隨機點集)和現實對象采集的三維點集(簡稱為采樣點集)。

4.1實驗比較

為了對比在不同規模下的數據求取MEB的時間效率，設計了10組點集，點數n分別是100，200，300，400，500，1 000，2 000，3 000，4 000，5 000，每組數據分別進行500次實驗。隨機點集的生成范圍有2種，矩形域隨機點集和球形域隨機點集。

實驗 1. 是矩形域隨機點集，隨機點按照均勻分布產生于一個單位正方體，如下式：

表 1是在矩形域上采樣所得出的實驗數據，Red1表示iRIA相對于RIA的時間壓縮比，Red2表示R-RCF相對于RIA的時間壓縮比。從表1可看出，本文提出的 iRIA算法在時間代價方面能壓縮 8%以上，且隨著點數的增加，時間效率的提升越來越高。而本文提出的 R-RCF算法在時間代價方面能壓縮66%以上，這將大幅度地提高求解MEB的時間運行效率。MEB的實驗效果如圖2(a)所示。

表1　矩形域上不同點數的平均運行時間比較

實驗2. 針對球形域隨機點集進行，隨機點按照均勻分布產生于一個平移后的單位球。

球形域上MEB的求取效果如圖2(b)所示。實驗結果見表2。從時間效率上看，本文所提出的iRIA算法由于第一步計算較遠的兩個點增加了不少時間，同時，由于球形域點集形狀對稱的特點，本文改進后的 iRIA算法并不能有效地減少迭代次數，所以總時間代價上表現的略差些。但是本文提出的R-RCF算法在計算球形域點集的MEB在時間代價上能壓縮71%以上。

圖2　離散點集最小包圍球計算結果

表2　球形域上不同點數的平均運行時間比較

在表2中，由于RIA和iRIA的運行時間相差并不大，故進行了方差未知均值相等(H0:μ1=μ2)的假設驗證。當取顯著性水平為0.05時，接受H0，即認為此兩種算法針對球形區域的時間效率無顯著差異。

4.2采樣點集的最小包圍球

在三維游戲制作或者三維激光測量中，需獲得各種形狀的點云數據。為了驗證RIA、iRIA和R-RCF針對不同復雜數據的時間效率，將進行MEB計算，針對6個現實對象采集的三維點集，分別為網絡公開數據庫(http://visionair.ge.imati.cnr.it/ontologies/shapes/)人手、花瓶和佛像模型的采樣點集，以及實驗室掃描數據：樹1、樹2和場景模型，如圖3所示。

圖3　現實三維模型采集數據的最小包圍球

6組數據的點數和3種算法的時間及其效率比較見表3，其中運行時間是實驗重復30次的平均運行時間。從表3看出，本文提出的改進算法iRIA的平均運行時間比RIA快10%～33%，R-RCF比RIA 快51%～93%。具體來說，iRIA與RIA之間的比較效果：對于比較“修長”的點集數據模型，如人手模型、樹1模型和樹2模型，iRIA的程序運行效率比RIA提高了22%以上，這是因為對于該類的點集數據模型，iRIA算法改進的初始值能大大減少求解MEB的迭代次數，這與矩形域隨機點集得出的結論一致；而對于一些比較“飽滿”的點集數據模型，如花瓶模型、佛像模型和場景模型，iRIA的程序運行效率比RIA提高的程度明顯低于“修長”的點集數據模型，這與球形域隨機點集得出的結論一致。從表3中的實驗結果可以看出，針對不同形狀的采樣點集，iRIA算法的運行效率普遍高于RIA。R-RCF與RIA之間的比較效果：對于比較“飽滿”的點集數據模型，R-RCF的程序運行效率比RIA提高了89%以上。可見，針對不同形狀的大規模點集，本文提出的新算法R-RCF的運行效率遠高于原始算法，同時，對于“飽滿”的點集數據模型，效果尤其明顯，這將為計算離散點集的MEB的工程應用節省大量的時間。

表3　3種算法針對激光掃描點云數據的平均運行時間比較

5　結 束 語

針對三維空間中求取離散點集MEB問題，為便于構造求解MEB的算法，本文先概述了三維空間中離散點集的MEB的幾條主要性質。在分析求解MEB的RIA的基礎上，提出一種稱為較大初始MEB的 RIA，這種算法通過在第一階段選取較遠點對構造初始 MEB，有效地減少了第二階段的MEB更新次數，提高了算法的時間效率。同時，本文又提出了一種新算法，稱之為R-RCF，該算法通過減少迭代過程中MEB的更新次數，大幅度提高了時間效率。通過隨機數據和現實三維模型采樣數據等多組實驗驗證，結果表明，與原始RIA相比，本文提出的R-RCF算法能大幅度的提高時間效率。本文算法分析詳細，有助于學生或者工程技術人員在算法實現過程中較好地處理相關細節。

此外，因為最小包絡長方體比MEB更加接近物體表面，所以參考本文的研究思路和二維空間中最狹長包絡矩形的求解原理及方法[15]，研究三維空間的最小包絡長方體的快速計算是未來的一個研究方向。

[1] Hu P, Xu D, Li H. A normalization method of moment invariants for 3D objects on different manifolds [J]. Computer Aided Drafting, Design and Manufacturing, 2014, 24(2): 15-22.

[2] 唐磊, 施侃樂, 雍俊海, 等. 模型適應的凸包圍多面體并行生成算法[J]. 中國科學: 信息科學, 2014, 44(12): 1515-1526.

[3] Samet H. The design and analysis of spatial data structures [M]. New Jersey: Addison Wesley Press, 1990: 86-89.

[4] 來疆亮, 王守覺. 最小球覆蓋幾何算法及其在模式識別中的應用[J]. 模式識別與人工智能, 2006, 19(2): 271-276.

[5] Berg M D, Jreveld M V. Overmars M, et al. 計算幾何:算法與應用[M]. 鄧俊輝, 譯. 北京: 清華大學出版社, 2005: 99-103.

[6] Ben-Hur A, Horn D, Siegelmann H T, et al. Support vector clustering [J]. Journal of Machine Learning Research, 2001, 2(2): 125-137.

[7] 耿博, 張慧娟, 王衡, 等. 離散曲面的近似 Poisson盤采樣[J]. 中國科學: 信息科學, 2012, 42(6): 703-716.

[8] Sylvester J J. A question on the geometry of situation [J]. The Quarterly Journal of Pure and Applied Mathematics, 1857, (1): 79.

[9] Sylvester J J. On Poncelet’s approximate linear valuation of surd forms [J]. Philosophical Magazine, 1860, 20(132): 203-222.

[10] 張勇, 陳強. 一種基于計算幾何方法的最小包容圓求解算法[J]. 工程圖學學報, 2007, 28(3): 97-101.

[11] Megiddo N. Linear-time algorithms for linear programming in R3 and related problems [J]. Siam Journal of Computing. IEEE, 1983, 12(4): 329-338.

[12] Welzl E. Smallest enclosing disks (balls and ellipsoids) [J]. Lecture Notes in Computer Science, 1991, 555(1): 359-370.

[13] 李紅軍, 張曉鵬. 離散點集最小包圍圓算法分析與改進[J].圖學學報, 2012, 33(2): 34-38.

[14] 汪衛, 王文平, 汪嘉業. 求一個包含點集所有點的最小圓的算法[J]. 軟件學報, 2000, 11(9): 1237-1240.

[15] 周敏, 鄭國磊, 陳樹林. 二維圖形最狹長包絡矩形的求解原理及方法[J]. 圖學學報, 2013, 34(4): 46-53.

An Im proved Random ized Incremental A lgorithm for Generating M inimum Enclosing Ball of Discrete Point Set

Li Shilin,Li Hongjun

(College of Science, Beijing Forestry University, Beijing 100083, China)

The minimum enclosing ball (MEB) of discrete point set in the 3D space has been w idely used in many fields, such as collision detection, computational geometry, pattern recognition and so on. In order to better understand and construct MEB, its character is analyzed firstly. A fter analyzing the classical randomized incremental algorithm, we propose two strategies to improve time efficiency. One strategy is to construct a bigger initial enclosing ball and another strategy is to decrease the number of updating MEB. Based on the second strategy, a new algorithm called random point group-recalculation farthest point is proposed. Multi-group experimental results, whose data are both random ly generated by computer and realistic 3D model sampling, show that the random point group-recalculation farthest point algorithm can effectively improve the time efficiency compared with the classical algorithms.

minimum enclosing ball; randomized incremental algorithm; random point grouprecalculation farthest point

TP 391

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

A

2095-302X(2016)02-0166-06

2015-09-24；定稿日期：2015-10-03

國家自然科學基金項目(61372190)

李世林(1992–)，男，山東青島人，碩士研究生。主要研究方向為計算幾何。E-mail：lishilin_2015@bjfu.edu.cn

李紅軍(1969–)，男，湖南郴州人，副教授，博士。主要研究方向為計算機圖形學和計算幾何等。E-mail：lihongjun69@bjfu.edu.cn