Finding planar regions in a terrain

Finding Planar Regions in a Terrain Stefan Funke 1 , Theo haris Malamatos, Rahul Ray Max-Plank-Institut fur Informatik, Stuhlsatzenhausweg 85, 66123 Saarbr uken, Germany Abstrat We onsider the problem of omputing large onneted regions in a triangulated terrain of size n for whih the normals of the triangles deviate by at most some small xed angle. In previous work an exat near-quadrati algorithm was presented, but only a heuristi implementation with no guarantee was pratiable. We present a new approximation algorithm for the problem whih runs in O ( n= 2 ) time and|apart from giving a guarantee on the quality of the pro dued solution|has b een implemented and shows go o d p erformane on real data sets representing frature surfaes with around half a million triangles. 1. Intro dution A terrain is a surfae in R 3 dened by a fun- tion f : R  R ! R . If f is pieewise linear and the surfae onsists of a olletion of triangles, the terrain is alled a triangulated irregular network (TIN) . Given a TIN T , the goal is to nd large, nearly planar regions in T . We dene `near planarity' as follows: for a real parameter Æ > 0, we say that a subset of triangles T in T is Æ -planar if (i) the adjaeny graph of the triangles in T is onneted, and (ii) there is a vetor ! r suh that for eah t 2 T , \ ( n t ; r )  Æ , where n t denotes the normal of triangle t . We all ! r the referene normal . (Throughout the pap er \ ( v ; u ) denotes the angle b etween two vetors ! v and ! u .) Researhers in the material sienes examine surfae top ographies of materials as they provide useful information ab out the generation pro ess and the internal struture of the material. Surfaes generated by frature, wear, orrosion and mahining are of interest. Among many other riteria, they want to examine feature-related parameters like faets in brittle frature surfaes. This appliation motivated our work. Email addresses: funkempi-sb.mpg.de (Stefan Funke), tmalamatmpi-sb.mpg.de (Theo haris Malamatos), rahulmpi-sb.mpg.de (Rahul Ray). 1 Part of this researh was onduted while the author was visiting the University of Illinois at Urbana-Champaign, USA. Related work. Assume that eah triangle in T is assigned a weight and that any set of triangles has weight equal to the sum of the weights of its triangles. Smid, Ray, Wendt and Lange [2℄ onsidered the problem of nding a Æ -planar subset of T of maximum weight. They presented an O ( n 2 log n (log log n ) 3 ) algorithm, where n is the number of triangles in T . Sine this is impratial they prop osed a heuristi whih omputes nearly planar regions quite quikly but with no guarantee. Also they suggested that nding a Æ -planar set with weight at least a onstant fration of the optimum may be equally diÆult as solving the problem exatly, see [2℄. Our results. In this pap er we adopt the following notion of approximation. Given a real parameter  > 0, we say that a subset of triangles T of T is  - approximate Æ -planar if it is Æ (1 +  )-planar and has weight at least as large as an optimal Æ -planar set. Under this notion, we present a new approximation algorithm that runs in O ( n= 2 ) time whih is independent of Æ . For n suÆiently large, the algorithm uses optimal O ( n ) spae. Reently HarPeled and Mazumdar [1℄ have used a similar notion of approximation for the problem of omputing the smallest k -enlosing disk. We have implemented and empirially evaluated the algorithm on real test data from the appliation domain onsisting of frature surfae terrains with more than 500,000 triangles. It provides very go o d quality results in reasonable time. Submitted paper to 20th EWCG Seville, Spain (2004) 20th Europ ean Workshop on Computational Geometry 2. Preliminaries Let T b e a TIN. We asso iate with T an undireted weighted graph G T ( V ; E ) as follows. Eah triangle t in T has an asso iated weight w ( t ) and orresp onds to a vertex v t in V . An edge onnets two verties of V if and only if the orresp onding triangles in T are adjaent. Note that G T is the dual graph of T , is planar and has degree three. Eah vertex v t 2 V is assigned the weight of its asso iated triangle w ( t ) whih an b e, for example, equal to the area of the triangle t or one (when we want to maximize the area of the deteted region or the number of triangles, resp etively). The weight w ( V 0 ) of any subset V 0 of V is dened as the sum of the weights of the verties in V 0 . For a p oint u 2 R 3 we denote by ! u the vetor ! O u , where O is the origin. Let S 2 denote the unit sphere, i.e., the boundary of the three-dimensional ball of radius one entered at the origin. For eah triangle t 2 T , we an asso iate a p oint p t 2 S 2 whih represents the normalized normal of triangle t . Speially, ! p t = ! n t = j ! n t j . Our goal is to approximate the spae S 2 of all p ossible normals by a - nite set of points V  S 2 suh that for any p 2 S 2 , there is a point v 2 V nearby. Denition 1 A set of point V  S 2 is al led a Æ  - disretization of S 2 if 8 p 2 S 2 : 9 v 2 V with \ ( v ; p )  Æ   . Lemma 2 There exists a Æ  -disretization of S 2 of size O (1 = ( Æ  ) 2 ) whih an be omputed in the same time. The following onstrution yields a Æ  -disretization of S 2 as needed in Lemma 2. (We omit its pro of in this abstrat.) Consider a ub e L with sidelength two entered at the origin. Note that S 2  L . Plae a 2-dimensional grid of size k  k with k = d p 2 = ( Æ  ) e over eah of the six faes of L . This generates k 2 equally sized square grid ells on eah fae of L , where eah ell has sidelength at most ( Æ  p 2), and 6 k 2 + 2 grid p oints overall. See Figure 1. Let A b e the set onsisting of these grid p oints. Our Æ  -disretization V of S 2 is dened as V = n ! z j ! z j : z 2 A o ; that is, for eah gridpoint z we sho ot a ray from the origin through z and we inlude the p oint where the ray leaves S 2 into set V . k k Fig. 1. Cub e with sidelength two ontaining S 2 and with a k  k grid on eah of its faes. 3. Finding Large Planar Regions 3.1. The Basi Algorithm We rst desribe a simple algorithm for the problem that omputes an  -approximate solution in O ( n= ( Æ  ) 2 ) time. The algorithm proeeds as follows: 1. Compute a Æ  -disretization V of S 2 . 2. For eah p 2 V , (a) Compute the set V p of verties v t with \ ( p; n t )  (1 +  )  Æ . (b) Consider the subgraph of G T indued by set V p and determine its heaviest onneted omp onent C p . 3. Rep ort the set of triangles T orresp onding to the heaviest omp onent C p found in Step 2 and the asso iated referene normal ! p . We now disuss the running time and orretness of the algorithm. Computing the Æ  -disretization takes O (1 = ( Æ  ) 2 ) by Lemma 2. For eah element p 2 V we determine the subgraph indued by V p and ompute its onneted omp onents, whih an b e done in O ( n ) time. So the total running time of Step 2 is O ( n= ( Æ  ) 2 ) whih also dominates the overall running time. For orretness, observe rst that learly the omputed set T is Æ (1 +  )-planar. It remains to show that the weight of T is at least that of an optimal Æ -planar set T  . For set T  there exists a vetor ! r  suh that for all triangles t 2 T  , \ ( r  ; n t )  Æ . Let p b e a point in V for whih \ ( p; r  ) = min u 2 V \ ( u; r opt ). By the denition of V , the angle \ ( p; r  ) must be at most Æ  . Then, for any triangle t 2 T  the angle b etween ! n t and ! p is at most Æ + ( Æ  ) = Æ (1 +  ). Therefore for all t 2 T  , Marh 25-26, 2004 Seville (Spain) v t 2 V p and hene, the algorithm will nd a onneted omp onent with at least the same weight. Lemma 3 Given a triangulated irregular network T , and two real parameters Æ > 0 and  > 0 we an ompute in O ( n= ( Æ  ) 2 ) time an  -approximate Æ -planar set of triangles in T . The running time of the basi algorithm is optimal in terms of n . But one may ask whether the dep endene on  or Æ an b e improved. In the following we will rene the algorithm to ahieve O ( n= 2 ) running time, whih is indep endent of Æ . 3.2. The Rened Algorithm The improvement in the running time of the algorithm omes from two fators. First we determine a set of referene normals V 0 of size at most O ( n= 2 ) that ontains all relevant referene normals, avoiding the insp etion of (1 = ( Æ  ) 2 ) p otential referene normals. Seond by a buketing sheme, we redue signiantly the numb er of times eah triangle is onsidered. The rened algorithm pro eeds as follows: 1. For eah triangle t 2 T with normal n t , let p t b e a p oint in V for whih \ ( p t ; n t ) = min u 2 V \ ( u; n t ); store t in the buket asso iated with p t . 2. Determine a set V 0  V of p otential referene normals as V 0 = f p 2 V : 9 p t with non-empty buket and \ ( p; p t )  (1 + 2  )  Æ g . 3. For eah v 2 V 0 , (a) Collet the set of triangles N v whih are ontained in bukets of referene normals r 0 2 V 0 with \ ( r 0 ; v )  (1 + 2  )  Æ . (b) Prune N v keeping only those triangles t with \ ( n t ; v )  (1 +  )  Æ . Let N 0 v b e the pruned set. () Consider the subgraph of G T indued by the verties orresponding to triangles in N 0 v and determine its heaviest onneted omp onent C v . 4. Output the heaviest onneted omp onent C v found in Step 3. Before analysing the algorithm, we state a small lemma whih informally says that in the Æ  -disretization, for any p oint p there is a small number of other p oints nearby. Lemma 4 Let p be a point in the Æ  -disretization V . Then the number of points p 0 2 V with \ ( p; p 0 ) < (1 + 2  )  Æ is at most O (1 = 2 ) . The pro of of this lemma follows easily from our onstrution of V and it is omitted. We note also that given a triangle normal n t we an ompute the point p t with \ ( n t ; p t ) = min u 2 V \ ( n t ; u ) in onstant time by rst determining whih fae of the ub e L is hit by the ray ! n t and then lo ating the p osition of the intersetion p oint within the grid on that fae. We now state the main result of this setion: Theorem 5 Given a triangulated irregular network T , and two real parameters Æ > 0 and  > 0 we an ompute in O ( n= 2 ) time an  -approximate Æ -planar set of triangles in T . PROOF. We laim that the output of the rened algorithm give us suh a set. To show orretness, it suÆes to prove that the algorithm omputes the same solution as the basi algorithm. We leave this pro of to the reader. We fo us now on the running time. Step 1 of the algorithm takes O ( n ) sine for eah t we an in onstant time determine p t and a- ess the asso iated buket using a hashing sheme. Step 2, where we form the set V 0 , takes O ( n= 2 ) sine there are at most n non-empty bukets. For eah of the non-empty bukets, we explore O (1 = 2 ) p oints in the neighb orho o d by Lemma 4. Finally, for the overall running time of Step 3, observe that again by Lemma 4, eah triangle an be olleted by at most O (1 = 2 ) referene normals and that the running time of one iteration of Step 3 is O ( j N v j ). Therefore Step 3 takes O ( n= 2 ) time in total. 2 4. Exp erimental Results We have implemented the rened algorithm of Setion 3.2 in C++ using the LEDA library [3℄. Exp eriments were arried out on a 1.8 GHz Pentium 4 mahine with 256 MB of RAM. We used several data sets representing frature surfaes from the appliation domain in material sienes. Input data were given as 512  512 raster images with the intensity of eah pixel orresponding to its height value. To obtain the TIN, we triangulated the p oint set by reating triangles ( i; j ) ; ( i + 1 ; j ) ; ( i + 1 ; j + 1) and ( i; j ) ; ( i; j + 1) ; ( i + 1 ; j + 1). See Figure 2. To sp eed-up our program (while preserving the same guarantee), we have ome up with the following three heuristis: 20th Europ ean Workshop on Computational Geometry (i,j+1) (i+1,j+1) (i+1,j)(i,j) Fig. 2. Triangulation sheme for the array of height values. ( i ) With eah potential referene normal we asso- iate the weight of all triangles in bukets at distane at most (1 + 2  )  Æ and examine the referene normals in dereasing order of weight. If the urrent best solution exeeds the asso iated weight of the next referene normal, we an stop examining further. ( ii ) If the normal of a triangle t forms an angle larger than 2(1 +  ) Æ with the normal of eah of the three adjaent triangles, we an prune t out (at Step 1). It an only form a singleton solution. ( iii ) In Step 3(b) where we ompute the set of relevant triangles N 0 v we only hek triangles in bukets at angle distane more than Æ , thus saving some exp ensive oating-p oint op erations. In our exp eriments, the three heuristis ombined redue the running time by nearly a fator of two. For our test instanes onsisting of roughly 500 ; 000 triangles and a hoie of Æ = 0 : 2 (ab out 11.5 degrees) and  = 0 : 2, the running time varied b etween 70 and 110 se. The best solution was in all ases found within the rst 20 seonds due to the prioritization sheme. Thus most of the running time was sp ent on heking that no b etter solution exists. (A heuristi ould just rep ort the solution omputed, for example, after 30 seonds.) We also ran exp eriments to determine the variations of the running time as a funtion of n and  . The dep endene graph on  for a xed data set is shown in Figure 3. The upp er urve denotes the total running time and the lower urve denotes the time when the reported solution was deteted. In Figure 4 we examine the dep endene on n (number of triangles). Note that again the rep orted solution was found after only few seonds. For the heuristi implementation in [2℄, running times rep orted are in the range of 30{40 seonds but as it is mentioned the solution omputed may b e far from optimal. 1 10 100 1000 10000 0 0.3 0.6 0.9 1.2 1.5 CPU Seconds Epsilon total time detection time Fig. 3. CPU time versus  for n = 522K, Æ = 0 : 2. 0 10 20 30 40 50 60 70 80 90 100 0 100 200 300 400 500 CPU Seconds Data Size (K) total time detection time Fig. 4. CPU time versus n for Æ = 0 : 2,  = 0 : 2. Referenes [1℄ S. Har-Peled and S. Mazumdar. Fast algorithms for omputing the smallest k -enlosing dis. In Pro. 11th Annu. European Sympos. Algorithms , Leture Notes Comput. Si. Springer-Verlag, 2003. [2℄ K. Lange, R. Ray, M. Smid, and U. Wendt. Computing large planar regions in terrains. In Proeedings of the 8th International Workshop on Combinatorial Image Analysis (IWCIA) , Eletronis Notes in Computer Siene, Volume 46, 2002. [3℄ K. Mehlhorn and S. Naher. LEDA: A Platform for Combinatorial and Geometri Computing . Cambridge University Press, Cambridge, U.K., 1999. [4℄ U. Wendt, K. Lange, M. Smid, R. Ray, and K. Tonnies. Surfae top ography quantiation by integral and feature-related parameters. Materialwissenshaft und Werkstotehnik , 33(10):621{627, Otob er 2002.