Rapid Delaunay triangulation for randomly distributed point cloud data using adaptive Hilbert curve

Author's Accepted Manuscript

Rapid Delaunay Triangulation for Randomly Distributed Point Cloud Data Using Adaptive Hilbert Curve Tianyun Su, Wen Wang, Zhihan Lv, Wei Wu, Xinfang Li

www.elsevier.com/locate/cag

PII: DOI: Reference:

S0097-8493(15)00122-3 http://dx.doi.org/10.1016/j.cag.2015.07.019 CAG2631

To appear in:

Computers & Graphics

Received date: 12 April 2015 Revised date: 13 July 2015 Accepted date: 14 July 2015 Cite this article as: Tianyun Su, Wen Wang, Zhihan Lv, Wei Wu, Xinfang Li, Rapid Delaunay Triangulation for Randomly Distributed Point Cloud Data Using Adaptive Hilbert Curve, Computers & Graphics, http://dx.doi.org/10.1016/j. cag.2015.07.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Rapid Delaunay Triangulation for Randomly Distributed Point Cloud Data Using Adaptive Hilbert Curve Tianyun Sua,1 , Wen Wangb,∗, Zhihan Lvc , Wei Wub , Xinfang Lia a Marine

Information and Computation Center, First Institute of Oceanography, State Oceanic Administration, Qingdao, China, 266061 b College of Information Science and Engineering, Ocean University of China, Qingdao, China, 266100 c Shenzhen Institutes of Advanced Technology (SIAT), Chinese Academy of Science, Shenzhen, China, 518055

Abstract Given the enormous scale and diverse distribution of 2D point cloud data, an adaptive Hilbert curve insertion algorithm which has quasi-linear time complexity is proposed to improve the efficiency of Delaunay triangulation. First of all, a large number of conflicting elongated triangles, which have been created and deleted many times, can be reduced by adopting Hilbert curve traversing multi-grids. In addition, searching steps for point location can be reduced by adjusting Hilbert curve’s opening direction in adjacent grids to avoid the “jumping” phenomenon. Lastly, the number of conflicting elongated triangles can be further decreased by adding control points during traversing grids. The experimental results show that the efficiency of Delaunay triangulation by the adaptive Hilbert curve insertion algorithm can be improved significantly for both uniformly and non-uniformly distributed point cloud data, compared with CGAL, regular grid insertion and multi-grid insertion algorithms. Keywords: Delaunay triangulation, adaptive Hilbert curve, grid division, multi-grid, point cloud data

1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1. Introduction As an important researching content in computational geometry, Delaunay triangulation is applied widely in many areas, including geographical information system, digital elevation model, finite element analysis, surface reconstruction [1], etc. Nowadays, algorithms in Delaunay triangulation field can be divided into several categories, namely, divide-and-conquer algorithm [2], incremental insertion algorithm [3], triangle expanding algorithm [4], sweep line algorithm [5], gift wrapping algorithm [6] and convex hull based algorithm [7]. According to the comparison of different Delaunay algorithms by Su et al. [8], these algorithms have their own characteristics, which can meet the requirements in different areas. With the rapid development of computer hardware, parallel computation based Delaunay algorithm [9, 10, 11, 12, 13] is attracting more and more attention in recent years. But the requirement for multi-

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

∗ Corresponding

author. Tel:+86 15726269702 Email addresses: [email protected] (Tianyun Su), [email protected] (Wen Wang), [email protected] (Zhihan Lv), [email protected] (Wei Wu), [email protected] (Xinfang Li) 1 Tel: +86 18669789597 Preprint submitted to Computers & Graphics

39 40 41 42 43

core CPU computers is a limitation for wide use of parallel Delaunay triangulation. Among these algorithms, the incremental insertion algorithm has become popular for its specific advantages—simplicity, small space occupancy and convenience for dynamic update. With the successive developments by Lee and Schachter [14], Bowyer [15], Watson [16], Sloan [17], Macedonio and Pareschi [18], Floriani [19] and Tsai [20], the incremental insertion algorithm is formed with the major processes of the bounding box creating, point locating, cavity expanding and triangle mesh updating [21]. On this basis, several approaches were proposed to improve the incremental insertion algorithm. The grid index scheme [22] can increase point location speed, but the overall efficiency will be reduced in updating grid index every time when one point is inserted. If the grid index scheme is not used, the algorithm efficiency will greatly depend on the order of the inserted points, which attracts much attention for optimization. Amenta et al. [23] presented a biased randomized insertion order (BRIO) in 2003, which can decrease the frequency of circumcircle judgement. But it wastes large amount of time on point location, namely, searching the triangle which contains the inserted point. Sloan [17] put forward a method to insert points by regular grids, July 17, 2015

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74

75 76

77 78 79

80 81 82 83 84

85 86 87 88 89 90 91

which can nearly lower down the time complexity to linear. Based on such regular grid insertion algorithm, Liu & Snoeyink [24], Zhou & Jones [25], Hornus & Boissonnat [26] and Buchin [27, 28] proposed the scheme of traversing grids according to various space-filling curve (e.g. Hilbert curve and Peano curve) sequence, which further enhances the efficiency of regular grid insertion algorithm. But the regular grid insertion scheme only has good effects on uniformly distributed point set. Devroye et al. [29, 30] presented k-d tree (short form for k dimensional tree) grid division scheme. Compared with the regular grid insertion scheme, the k-d tree scheme is more efficient for non-uniformly distributed point set. However, as the result of high complication of the k-d tree, constructing k-d tree grids would cost much more time and lead to low efficiency ultimately. On that basis, S.H. Lo [31] proposed multi-grid insertion scheme that is not only more simple than k-d tree insertion scheme but also can extremely improve the efficiency of Delaunay triangulation for non-uniformly distributed points. Nonetheless, because of the line-by-line grid traversing manner, the multi-grid insertion will generate many conflicting elongated triangles, which cost extra time on constructing and deleting, and thus decrease the efficiency finally. Based on the analysis of the factors influencing the efficiency of incremental insertion algorithm, this paper proposes an algorithm with adaptive Hilbert curve insertion to improve the efficiency of the incremental Delaunay triangulation for both uniformly and non-uniformly distributed point cloud data.

each other. In this way, however, much more elongated triangles will be generated and deleted while each point 94 is inserted. This will cause extreme rise for the time of 95 Part 3. 96 For Part 3, the time can be minimized when points are 97 inserted at random. But the random insertion scheme 98 will make inserted point far away from the last updated 99 triangle, which results in a sharp increase for the time 100 of Part 2. 101 Based on the analysis of the influencing factors for 102 incremental insertion algorithm, the regular grid inser103 tion scheme can well balance the time consumption be104 tween Part 2 and Part 3 for uniformly distributed point 105 set [17]. However, for the reason of highly different 106 point number in each grid, the regular grid insertion 107 scheme is less efficient for non-uniformly distributed 108 point set. A multi-grid insertion scheme was proposed 109 by S.H. Lo in 2013 [31], which can enhance the effi110 ciency of Delaunay triangulation for non-uniformly dis111 tributed point set by dividing grids iteratively according 112 to the point number of each grid cell. With the compari113 son of S.H. Lo, the regular grid insertion and multi-grid 114 insertion scheme have the same efficiency for uniform 115 point distribution. While for non-uniform point distri116 bution, the efficiency of multi-grid insertion scheme can 117 have nearly 10 fold improvement compared with the 118 regular grid insertion scheme. 92 93

2. The incremental insertion algorithm based on point sorting The running time of the incremental insertion algorithm based on point sorting can be mainly divided into 3 parts. Part 1: preprocessing the order of points to be inserted; Part 2: seeking out triangles where inserted points are located; Part 3: updating the triangulation mesh by BowyerWatson scheme. Among them, the time of Part 1 can be neglected compared with that of Part 2 and Part 3 when point sorting process is not very complex. Therefore, efficiency improvement for Part 2 and Part 3 are mainly concerned in this paper. During the process of Part 2, the running time will decrease remarkably when the inserted points are close to

Figure 1: Grid traversing order of multi-grid

As shown in Figure 1, the grid traversing order of the multi-grid insertion scheme is line-by-line. How121 ever, there are some disadvantages for this kind of grid 122 traversing order. Firstly, a large number of elongated 123 triangles connecting with the boundary points of pre124 vious grid line (red triangles shown in Figure 2) will 125 be generated, which can increase the time of Part 3 in 126 constructing and deleting these triangles. Furthermore, 127 since these elongated triangles are very dense, the time 128 of Part 2 will increase as well when seeking out the 129 triangle where the inserted point is located. Secondly, 130 during construction of Delaunay triangulation, the in119 120

2

serted points are usually linked with four vertices of the bounding box to generate many elongated triangles 133 which need to be deleted (Figure 3). Thirdly, the first in134 serted point of a grid may have relatively long distance 135 with the last inserted point in adjacent grid and has to 136 traverse a lot of triangles (red triangles shown in Figure 137 4) to find its target triangle. Thus, it can cut down the 138 efficiency of Part 2. 131 132

146

3.1. Rule of grid division

In order to average the number of points in every grid, higher-level grids will be divided when the number of 0 149 points in a grid exceeds Nmax . Let G be the grid formed 150 in the area where original point set is located. Then, 151 an iteration process will be carried out for each divided 0 152 grid unit of G during grid division and the relationship 153 among the generated grids in the iteration process is as 154 follows (m is the current grid level when grid division is 0 1 m−1 155 stopped): G ⊃ G ⊃ · · · G ⊃ Gm . 156 In order to satisfy traversing requirements of the 157 Hilbert curve, the total grid number Ngrid and the grid 158 number of each column and row Nside at each level must 159 be guaranteed in accordance with Formula (1) and For160 mula (2) respectively, where r is the order of divided 161 grid with Hilbert curve, r ∈ N. 147 148

Ngrid = 4r Nside = 2r

Figure 2: Elongated triangles generated with the boundary points of previous grid line

(1) (2)

For one grid unit Gm−1 , it should be divided into Gm 163 while its point number Np satisfies the grid division con164 dition as shown in Formula (3). 162

Np > Nmax 165 166

And the average point number Nave of each grid unit of Gm can be calculated according to Formula (4) Nave =

Figure 3: Elongated triangles generated with the vertices of bounding box 167

(3)

Np Np = r Ngrid 4

So, r can be obtained by Formula (5).   Np r = log4 Nave

(4)

(5)

According to the research of S. H. Lo [31], the average point number Nave in every grid unit Gm should be 170 set from 4 to 16. So, after r is calculated from Formula 171 (5), a judgment should be performed as follows: 168 169

If

Figure 4: Points inserted in adjacent grids (a and b) far from each other

Np > 16, then r = r + 1. 4r

After grid division, points need to be assigned to the grid units. 174 Let Pm−1 = {pm−1 = (xi , yi ) | i = 1, 2, · · · Np } be the i m−1 175 points which are located at G and xmax , xmin , ymax , m−1 176 ymin be the maximum and minimum values of G in 177 x-axis and y-axis respectively. The span of x-axis Xspan 178 and the span of y-axis Yspan can be defined as follows: 172

139

3. Algorithm improvement

In view of the disadvantages of the incremental insertion algorithm and advantages of the locality-preserving 142 behavior of Hilbert curve, an adaptive Hilbert curve in143 sertion algorithm is proposed in this paper, which di144 vides grids with multi-grid rules and then traverses grids 145 in terms of adaptive Hilbert curve. 140 141

173

Xspan = xmax − xmin , Yspan = ymax − ymin 3

179 180

And the row width Wr and the column width Wc can be calculated with the following formulas: Wr =

Yspan Xspan , Wc = . Nside Nside

in Pm−1 can be assigned to grid set So, the point pm−1 i 182 G = {(row, col) | row, col = 1, 2, · · · Nside } in accor183 dance with the following rules: 181

m

row =

bottom-right corner. As shown in Figure 6, for two continuous Gm units, the last point p1 in prior unit may be 214 far from the first point p2 in the next unit with normal m+1 215 traversing sequence of G . In this case, the calcula216 tion time will increase during seeking out the triangle 217 where the point p2 is located. Thus, a rotation rule for 218 the opening direction of Hilbert curve is proposed as be219 low to solve this problem. 212 213

yi − ymin xi − xmin + 1, col = +1 Wr Wc

Based on this grid division rule, all points can be iteratively assigned to grid unit Gm from m = 1 until the m 186 point number of each G is not satisfied with Formula 187 (3). 184 185

188

3.2. Adaptive Hilbert curve

Hilbert curve is a continuous fractal space-filling curve which was first described by the German math191 ematician David Hilbert in 1891 [32]. It has locality192 preserving behavior and can fill in the whole space [32, 193 33, 34]. Nowadays Hilbert curve has been widely used 194 in many areas such as surface sampling [35], image 195 management, spatial database indexing and so on. 196 In view of improving the efficiency of Delaunay trian197 gulation, the two-dimensional adaptive Hilbert curve is 198 applied in this paper. Two-dimensional n-order Hilbert n n 199 curve Hn (n = 1, 2, 3 · · · ) can fill in a 2 × 2 plane re200 gion. With this curve, the square-shaped lattice area of 201 each grid unit can be filled and traversed linearly only 202 once without intersection. An example of Hilbert curve 203 with n = 1, 2 and 3 is shown in Figure 5. Hilbert curve 204 with n = 1 is also referred to as Hilbert cell. Due to 205 the locality-preserving behavior of Hilbert curve, it can 206 make continuous inserted points not very far from each 207 other and reduce the number of elongated triangles to 208 improve the efficiency of Delaunay triangulation. 189 190

(a)

(b)

(c)

Figure 5: First-third order Hilbert curve. (a) Hilbert cell/first-order Hilbert curve; (b) Second-order Hilbert curve; (c) Third-order Hilbert curve

However, the “jump” phenomenon will happen if Hilbert curve is applied in multi-grid and the traversing 211 sequence of Hilbert curve is from bottom-left corner to 209 210

Figure 6: The “jump” phenomenon. (Blue dash square means grid unit Gm ; red dash square means grid unit Gm+1 )

Definition 1 (two-dimensional multiple Hilbert gene). A two-dimensional multiple Hilbert gene is a series of 222 transformation commands, which controls the traversm+1 223 ing sequence of G when it is divided from area S mq m 224 where the present G unit is located. 220 221

Let S mq−1 and S mq+1 be the areas prior and next to Gm 226 unit respectively in terms of the traversing sequence, the q−1 q q+1 227 location relationship of S m , S m and S m can be rep228 resented with 20 kinds of permutations. Based on these 229 permutations, the multiple Hilbert gene list (as shown in 230 Table 1) is proposed to define the transformation of the 231 Hilbert curve’s opening direction for each permutation 232 in order to avoid the “jump” phenomenon. 233 In Table 1, the column “The transformation for Gm+1 q 234 in S m ” means the rotation transformation based on the 235 sequence of the Hilbert curve generated by Hilbert cell 236 in Figure 5a. The grid traversing sequence of trans237 formed Hilbert curve is defined in the column “Traversm+1 238 ing direction for G in S mq ”, where “reverse the di239 rection” means the traversing sequence of Hilbert cell m+1 240 in G is opposite to the normal sequence as shown in q−1 241 Figure 5a. The position relationship between S m and q 242 S m is shown in Figure 7, and the character “-” in the first q 243 column of Table 1 means the grid in S m is the first unit q q+1 m 244 of G . The position relationship between S m and S m q−1 q 245 is similar to that of S m and S m , and the character “-” q 246 in the second column of Table 1 means the grid in S m is m 247 the last unit of G . 225

4

Table 1: Multiple Hilbert gene list (The position relation is shown in Figure 7)

Position relation of S mq−1 and S mq Left to right Left to right Right to left Right to left Top to bottom Top to bottom Bottom to top Bottom to top Left to right Right to left Top to bottom Bottom to top Left to right Top to bottom Right to left Bottom to top

(a) Left to right

Position relation of S mq and S mq+1 Bottom to top Top to bottom Bottom to top Top to bottom Left to right Right to left Left to right Right to left Left to right Right to left Top to bottom Bottom to top Left to right Top to bottom Right to left Bottom to top -

The transformation for Gm+1 in S mq Rotate 90◦ clockwise Without change Rotate 180◦ clockwise Rotate 90◦ counterclockwise Without change Rotate 90◦ clockwise Rotate 90◦ counterclockwise Rotate 180◦ clockwise Without change Rotate 180◦ clockwise Rotate 90◦ counterclockwise Rotate 90◦ clockwise Without change Without change Rotate 90◦ counterclockwise Rotate 90◦ counterclockwise Rotate 180◦ clockwise Rotate 180◦ clockwise Rotate 90◦ clockwise Rotate 90◦ clockwise

Traversing direction for Gm+1 in S mq Reverse the direction Without change Without change Reverse the direction Without change Reverse the direction Reverse the direction Without change Without change Without change Reverse the direction Reverse the direction Without change Without change Reverse the direction Reverse the direction Without change Without change Reverse the direction Reverse the direction

(b) Bottom to top

Figure 8: The grid traversing sequence transformed by the multiple Hilbert gene list

(c) Right to left

(d) Top to bottom q−1

Figure 7: Position relationship between S m

q

and S m

Figure 8 shows the traversing sequence transformed by the multiple Hilbert gene list.The grids in Figure 8 250 are identical with those in Figure 6. As illustrated in 251 Figure 8, “jump” phenomenon is avoided successfully, 252 which results in decrease of the time consumption of 253 point location. 248 249

254

3.3. Control points

One major reason for generating elongated triangles 256 is that a vertex of the triangle is far away from the other 257 two. In order to reduce the generated elongated trian255

(a) Not adding control points

(b) Adding control points

Figure 9: The Delaunay triangulation when the 2000th point is inserted

gles, the control points are selected from all data point sets which are inserted into the triangle mesh in ad260 vance. In this way, the non-control points will be linked 261 with adjacent control points while inserted, which can 258 259

5

cut down the number of elongated triangles effectively. For every G1 unit, if there are one or more points in 264 it, the first point is selected as the control point. If not, 265 it needs to skip the unit and check the next one until all 1 266 G units have been judged. 267 It needs to point out that the selection sequence of 268 control points is opposite to the normal sequence of 269 Hilbert curve as illustrated in Figure 5a. By this means, 270 the last inserted point in control point set and the first 271 inserted non-control point can be in the same grid unit, 272 so as to quickly find the base triangle where the lat273 ter is located. Figure 9a shows the Delaunay triangu274 lation of 10000 uniformly distributed points when the th 275 2000 point is inserted without control point. Figure 276 9b shows the Delaunay triangulation of the same point 277 set when inserting the same point as illustrated in Fig278 ure 9a, while control points are added in advance. It 279 can be seen that adding control points can dramatically 280 reduce the number of points connecting to a large num281 ber of distant points through narrow elongated triangles 282 compared with the situation where control points are not 283 added.

point set; delete control point from the grid; 313 end if 314 if the point number in the ith Gm satisfies formula 315 (3), then 316 griddivide(m + 1); 317 else 318 put the points in the ith Gm unit into uncontrolled 319 point set; 320 end if 321 i + +; 322 end while 323 end

262

311

263

312

284

3.4. The algorithm flow

Combined with the improved algorithm above, the pseudo code of the overall sorting process of the algo287 rithm is shown as follows: 285 286

begin if the whole point number satisfies formula (3), 290 then 291 m ← 1; 292 griddivide(m); 293 else 294 put the whole points into uncontrolled point set; 295 end if 296 insert points into the triangle mesh from control 297 point set; 298 insert points into the triangle mesh from uncon299 trolled point set; 300 end 288 289

function griddivide(m) begin 303 divide Gm by the grid division rule; 304 assign points into Gm ; 305 determine the traversing order of Gm by the adaptive 306 Hilbert curve; 307 i ← 0; 308 while Gm has not been traversed over do 309 if m == 1, then 310 put the first point of the ith G1 unit into control 301 302

324

3.5. The time complexity analysis

The algorithm in this paper is mainly divided into two steps, i.e. point sorting and point inserting. 327 The point sorting step is a recursive procedure. Each 328 iteration includes four sub-steps: subdividing grids, as329 signing points into subdivided grids, sorting grids with 330 adaptive Hilbert curve and sorting points by the order of 331 the sorted grids. Among these sub-steps, the time com332 plexity of sub-dividing grid is O(1) and the time com333 plexity of the other three sub-steps is O(m), where m 334 is the point number in the grid to be divided. So, each 335 iteration of point sorting has a linear time complexity. 336 Due to the number of iterations is related to point distri337 butions rather than the point number, it can be inferred 338 that the time complexity of point sorting is O(n), where 339 n is the total point number of a point set. 340 Time complexity of the point inserting is relatively 341 complicated to analyze because it is hard to find the rela342 tionship between the number of triangles created and re343 moved repeatedly and the number of inserted points for 344 various kinds of point distributions. According to the 345 reference 31, the time complexity is quasi-linear only 346 for the distributions that don’t create many elongated 347 triangles. Based on the above analysis, the algorithm in 348 this paper can obviously reduce the number of elongated 349 triangles, which doesn’t vary evidently for each inserted 350 point with control points added in advance. Thus, the 351 point inserting step can be inferred to have quasi-linear 352 complexity no matter what the distribution is. 353 Therefore, the time complexity of the whole algo354 rithm proposed by this paper is quasi-linear. 325 326

355

356 357

6

4. Experiment and result analysis In order to test the efficiency of the algorithm proposed in this paper, an experiment is carried out for

different point distributions, including uniform distribution, non-uniform distribution (line, cross, cir360 cle, spiral) and mixed distribution (as shown in Fig361 ure 10a–10f) with CGAL, regular grid insertion al362 gorithm(RG), multi-grid insertion algorithm(MG) and 363 adaptive Hilbert curve insertion algorithm(HG) via a PC 364 with Intel(R) Xeon(R) CPU E5-1620 @3.60 GHz and 8 365 GB RAM running on Windows 7 64bit using VS2008. 366 According to the experiment, it is more efficient when 367 the value of Nmax is between 50 and 200, and 50 is used 368 in this paper. The results are shown in Table 2–6 and 369 Figure 11–15 for various point number from 2 to 10 milnd 370 lion. As mentioned in the 2 section, the time of incre371 mental insertion algorithm is mainly composed of Part 372 2 and Part 3. The time of Part 2 is proportional to the 373 average number (NP) of traversed triangles during find374 ing the base triangle where the inserted point is located. 375 The time of Part 3 is proportional to the average num376 ber (NR) of removed triangles for every inserted point. 377 In the case where the point sorting process is not very 378 complex, the smaller the NP and NR are, the higher the 379 efficiency is. Because there is no output about average 380 number of the traversed triangles and removed triangles, 381 the value of NP and NR of CGAL is not illustrated in 382 Table 2-6. 358 359

(a) Uniform

(b) Line

(c) Cross

(d) Circle

(e) Spiral

(f) Mixed

Figure 12: The efficiency of 4 algorithms for various distributions of 4 M points

Figure 13: The efficiency of 4 algorithms for various distributions of 6 M points

Figure 14: The efficiency of 4 algorithms for various distributions of 8 M points

Figure 10: Different point distributions (the point number is rarefied to 10000) Figure 15: The efficiency of 4 algorithms for various distributions of 10 M points

For uniform point distribution, as the point number in each grid is roughly equal, the multi-grid insertion is 385 not invoked because the threshold of 50 points is never 386 exceeded. Thus, NP and NR of RG and MG are iden387 tical for the insertion of 2–10 million points. For most 388 situations, the efficiency of CGAL is a little higher than 389 that of MG and RG. Because NP varies from 3 to 5 and 383 384

Figure 11: The efficiency of 4 algorithms for various distributions of 2 M points

7

Point distributions Uniform Line Cross Circle Spiral Mixed

Point distributions Uniform Line Cross Circle Spiral Mixed

Table 2: The efficiency of 4 algorithms for various distributions of 2 M points RG MG Time (ms) Time (ms) Total time (ms) NP NR NP NR NP Grid Insert Total Grid Insert Total 2808 5.34 5.89 219 2574 2840 5.34 5.89 265 2589 2870 3.91 2901 20.23 4.23 109 3525 3650 4.85 4.98 374 2309 2699 4.43 2855 14.87 4.41 94 3073 3198 4.44 5.21 405 2340 2761 3.93 2855 8.68 4.78 93 2480 2605 4.54 5.34 281 2277 2574 3.77 2870 7.38 6.35 94 2683 2793 4.46 5.83 265 2340 2605 3.73 2855 6.22 5.65 203 2574 2840 5.24 5.72 280 2527 2823 3.98

HG Time (ms) NR Grid Insert Total 4.15 219 2012 2309 4.54 374 2168 2574 4.54 390 2090 2527 4.36 296 1919 2262 4.73 265 1982 2278 4.15 234 1981 2293

Table 3: The efficiency of 4 algorithms for various distributions of 4 M points RG MG Time (ms) Time (ms) Total time (ms) NP NR NP NR NP Grid Insert Total Grid Insert Total 5647 5.85 5.69 327 5179 5569 5.85 5.69 406 5211 5632 4.51 5803 27.87 4.18 187 8549 8752 4.49 5.23 780 4711 5507 3.97 5679 20.29 4.34 187 7082 7285 4.89 5.17 718 4773 5507 4.39 5694 11.28 4.61 172 5367 5570 4.87 5.26 546 4602 5163 4.20 5725 9.04 6.52 172 5694 5897 4.45 5.84 484 4726 5210 3.77 5741 7.53 5.48 328 5289 5679 5.63 5.57 452 5288 5756 4.41

HG Time (ms) Grid Insert Total 4.12 359 4134 4571 4.73 795 4353 5195 4.55 718 4337 5086 4.41 546 4040 4633 4.87 499 4056 4586 4.13 405 4072 4555

CGAL

CGAL

Table 4: The efficiency of 4 algorithms for various distributions of 6 M points RG MG Point Time (ms) Time (ms) distributions Total time (ms) NP NR NP NR NP Grid Insert Total Grid Insert Total Uniform 8626 5.47 6.22 936 8003 9157 5.47 6.22 1155 8066 9236 3.79 Line 8798 18.59 4.32 374 10171 10592 4.66 5.19 1264 7114 8393 4.11 Cross 8626 13.79 4.56 375 9344 9766 4.50 5.42 1357 7254 8627 3.90 Circle 8612 8.74 5.03 328 7722 8112 4.65 5.59 1076 7129 8221 3.74 Spiral 8658 7.11 6.48 297 8159 8518 4.74 6.05 811 7316 8127 3.84 Mixed 8626 6.73 5.93 873 8065 9157 5.35 5.99 1217 7785 9017 3.88

NR

CGAL

CGAL Point distributions Total time (ms) Uniform Line Cross Circle Spiral Mixed

11623 11918 11590 11606 11668 11622

Table 5: The efficiency of 4 algorithms for various distributions of 8 M points RG MG Time (ms) Time (ms) NP NR NP NR NP Grid Insert Total Grid Insert Total 5.59 6.14 1170 10686 12059 5.59 6.14 1389 10749 12153 3.91 21.18 4.29 468 14618 15148 4.90 5.14 1607 9563 11186 4.36 15.61 4.52 468 12667 13198 4.45 5.45 1747 9719 11482 3.83 9.65 4.94 406 10608 11076 4.64 5.56 1357 9485 10858 3.74 7.62 6.53 359 11139 11560 4.55 6.00 1155 9672 10842 3.70 7.24 5.87 1092 10858 12152 5.44 5.94 1466 10358 11840 3.95

Table 6: The efficiency of 4 algorithms for various distributions of 10 M points RG MG Point Time (ms) Time (ms) distributions Total time (ms) NP NR NP NR NP Grid Insert Total Grid Insert Total Uniform 14461 5.72 6.07 1373 13323 14914 5.72 6.07 1622 13432 15070 4.05 Line 14929 23.46 4.27 577 19406 20030 4.80 5.22 2059 11965 14040 4.24 Cross 14508 17.21 4.49 577 16567 17207 4.58 5.44 2075 12200 14290 3.96 Circle 14508 10.47 4.87 530 13510 14102 4.76 5.51 1638 12028 13682 3.90 Spiral 14586 7.09 6.57 436 14149 14648 4.51 5.98 1419 12168 13603 3.68 Mixed 14539 7.71 5.81 1279 13821 15319 5.54 5.89 1685 13057 14757 4.06

NR 4.21 4.57 4.53 4.33 4.73 4.22

NR 4.17 4.56 4.57 4.37 4.75 4.19

CGAL

8

NR 4.15 4.62 4.57 4.38 4.78 4.17

HG Time (ms) Grid Insert Total 1014 6131 7363 1279 6630 7987 1373 6349 7784 1076 5804 6942 811 6006 6895 1045 6099 7378

HG Time (ms) Grid Insert Total 1216 8159 9609 1638 8736 10436 1763 8487 10312 1342 7972 9392 1139 8003 9220 1279 8127 9625

HG Time (ms) Grid Insert Total 1435 10218 11887 2106 10952 13120 2106 10639 12823 1607 9844 11529 1420 10280 11778 1528 10156 11918

A spiral distribution is locally similar to the circle distribution. CGAL and RG of spiral distribution have sim423 ilar efficiency. The efficiency of MG is approximately 424 5% higher than that of CGAL and RG. Compared with 425 MG, the NP and NR of HG are both roughly reduced by 426 20%, so the efficiency is about 15% higher than MG. 427 The test result of mixed point distribution is similar to 428 that of the uniform point distribution because every grid 429 has points for both distributions. Thus, triangles will 430 fill the whole plane after adding control points, which 431 results in less elongated triangles. The efficiency of 432 CGAL, RG and MG is similar to each other. Because 433 NP and NR of HG are reduced apparently, its efficiency 434 is 20% higher compared with RG and MG. 435 From the perspective comparison of different algo436 rithms, the total time of CGAL scarcely changes with 437 various point distributions, and the speed of construct438 ing Delaunay triangulation is approximately stable at 439 690 points/ms. The efficiency of RG alters dramati440 cally with different point distributions. The speed can 441 rise up to 710 points/ms for circle point distribution, and 442 drop down to 500 points/ms for line point distribution. 443 Similar to CGAL, the efficiency of MG changes mildly 444 with different point distributions. It is a little slower 445 for uniform and mixed point distributions (about 670 446 points/ms), and a little quicker for non-uniform point 447 distributions (about 720 points/ms). For HG, the effi448 ciency improvement for line and cross point distribu449 tions (about 770 points/ms) is slightly lower than that 450 for other point distributions (about 850 points/ms). 451 Figure 16 shows the total time of HG applied in dif452 ferent point distributions from 2 to 10 million points. 453 It can be seen that the time increase is nearly linear 454 with the point number, which is consistent with the time 455 complexity analysis in Section 3.5. It can be observed 456 that there are very small changes for the efficiency of 457 HG with various point distributions, illustrating that HG 458 is widely applicable for point cloud dataset no matter 459 uniformly or non-uniformly distributed. 460 Taking mixed distribution as an example, the effi461 ciency of four different algorithms for various point 462 numbers is shown in Figure 17. It can be seen that 463 the running time of HG has the slowest increasing rate, 464 which makes the superiority of HG more and more ob465 vious as the point number increase. 421 422

Figure 16: The efficiency for various distributions by HG

Figure 17: The efficiency of 4 algorithms for mixed distribution

NR varies from 4.1 to 4.2, which are obviously the lowest among these four algorithms, HG has the highest 392 efficiency compared with the other three algorithms for 393 various point number. The total time of HG on con394 structing Delaunay triangulation can decrease by 20% 395 compared with MG. 396 There is a great difference between RG and other 397 three algorithms for points distributing on a line. The 398 main cause for such a large difference is due to the time 399 consuming on finding the base triangle. NP of RG is 400 more than 18 and has the lowest performance in time 401 consumption compared with the other three algorithms. 402 The efficiency of MG is slightly higher than that of 403 CGAL. HG is the most efficient algorithm and the total 404 time is approximately 5% to 10% less compared with 405 MG for different point number. 406 The results for the cross distribution are very similar 407 to those of the line distribution. RG has the lowest effi408 ciency with the highest NP, which results in about 30% 409 more time compared with the most efficient algorithm 410 HG. With approximately 0.5 higher in NP and 1 higher 411 in NR, the total time of MG is about 10% longer than 412 HG. In addition, the efficiency of CGAL is between MG 413 and HG. 414 For circle point distribution, CGAL spends the 415 longest time for various point number. The total time 416 varies with different point number for RG but scarcely 417 changes for MG, because MG has a higher NP but a 418 lower NR. The NP and NR of HG vary from 3.7 to 4.4 419 and 4.3 to 4.5 respectively, leading to the highest effi420 ciency compared with the other three algorithms. 390 391

466

5. Conclusion

In terms of the problems that a large quantity of conflicting elongated triangles are generated and points in 469 adjacent grids are far from each other during construct470 ing Delaunay triangulation, an adaptive Hilbert curve 467 468

9

algorithm is proposed to improve the efficiency of Delaunay triangulation in this paper. By dividing multiple 473 grids, adopting adjusted Hilbert curve to traverse grids 474 and adding control points, the efficiency of Delaunay 475 triangulation is further improved. The experimental re476 sults indicate that the NP and NR values of the algo477 rithm proposed in this paper are lower than regular grid 478 insertion algorithm and multi-grid insertion algorithm 479 regardless of point distributions, which results in less 480 time consumption for Delaunay triangulation. Addi481 tionally, the efficiency of the algorithm proposed in this 482 paper is 10%-20% higher than CGAL, which is gener483 ally used as a fast method in Delaunay triangulation. 484 As a result, the grid traversing sequence can influ485 ence the efficiency of incremental Delaunay algorithm 486 in grid insertion scheme and the algorithm proposed in 487 this paper optimizes the sequence to have a better per488 formance for large-scale and randomly distributed point 489 cloud data. 490 As discussed by previous researches [36], there are 491 some degenerate cases during Delaunay triangulation 492 for a large-scale data. Actually, we have tested our al493 gorithm on the point sets which have many four-points494 co-circular cases or low-precision points, and the results 495 indicate that no degenerate case appears. Certainly, we 496 can’t guarantee that there are no other degenerate cases 497 in our algorithm while handling all kinds of point set. 498 Thorough analysis of degenerate cases of our algorithm 499 should be carried out in the future. 471 472

521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558

500

6. Acknowledgements

This research is supported by the project (Grant No. 502 201205001) of ”the Public Science and Technology Re503 search Funds Projects of Ocean” and sub-project (Grant 504 No. 2011ZX05056-001) of ”The Important Project of 505 Science and Technology in Developing Great Oil & Gas 506 Field and Coal Bed Gas”. 501

559 560 561 562 563 564 565 566 567 568 569 570

507

508 509 510 511 512 513 514 515 516 517 518 519 520

References [1] Y. J. Liu, Z. Q. Chen, K. Tang, Construction of Iso-contours, Bisectors and Voronoi Diagrams on Triangulated Surfaces. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(8) (2011) 1502-1517. [2] B. A. Lewis, J. S. Robinson, Triangulation of planar regions with applications. Computer Journal, 21(4)(1978) 324-332. [3] C. L. Lawson, Software for C1 surface interpolation, in: Proceedings of the Symposium on Mathematical Software, ACM, New York, NY, USA, 1977, pp. 161-194. [4] P. J. Green, R. Sibson, Computing Dirichlet tessellations in the plane. Computer Journal, 21(2)(1978) 168-173. [5] S. Fortune, A sweepline algorithm for voronoi diagrams. Algorithmica, 2(1987) 153–174.

571 572 573 574 575 576 577 578 579 580 581 582 583 584 585

10

[6] R. A. Dwyer, Higher-dimensional voronoi diagrams in linear expected time. Discrete & Computational Geometry, 6(1991) 343– 367. [7] C. Brad Barber. Computational geometry with imprecise data and arithmetic. PhD thesis, Princeton, 1993 [8] P. Su, R. L. S. Drysdale, A comparison of sequential Delaunay triangulation algorithms, Comput. Ge-om. 7(1997) 361385. [9] Y. K. Rui, J. C. Wang, A new study of compound algorithm based on sweepline and divide-and-conquer algorithms for constructing Delaunay triangulation. Acta Geodaetica et Cartographica Sinica, 36(3)(2007) 358-362. [10] S. H. Lo, Parallel Delaunay triangulation—application to two dimensions. Finite Elements in Analysis and Design, 62(2012) 37–48. [11] S. B. Bi, D. Q. Chen, J. Yan, Y. Guo, Planar Delaunay triangulation algorithm based on 2D convex hull. Computer Science, 41(10)(2014) 317–320. [12] J. Li, D. R. Li, Z. F. Shao, A streaming data Delaunay triangulation algorithm based on parallel computing. Geomatics and Information Science of Wuhan University, 38(7)(2013) 794-798. [13] Y. L. Cai, S. M. Xiao, L. Qi, Locking paralleled GPU-based method research for unstructured mesh generation. Computer Engineering and Applications, 50(6)(2014) 56-60. [14] D. T. Lee, B. J. Schachter, Two algorithms for constructing a Delaunay triangulation. International Journal of Computer and Information Sciences, 9(3)(1980) 219-242. [15] A. Bowyer. Computing Dirichlet tessellations. Computer Journal, 24(2)(1981) 162-166. [16] D. F. Watson. Computing the n-dimension Delaunay tessellation with application to Voronoi polytopes. Computer Journal, 24(2)(1981) 167-172. [17] S. W. Sloan, A fast algorithm for constructing Delaunay triangulations in the plane. Advances in Engineering Software, 9(1)(1987) 34-55. [18] G. Macedonia, M. T. Pareschi, An algorithm for the triangulation of arbitrarily distributed points: applications to volume estimate and terrain fitting. Computers & Geosciences, 17(7)(1991) 859-874. [19] L. De Floriani, E. Puppo, An on-line algorithm for constrained Delaunay triangulation. CVGIP: Graphical Models and Image Processing, 54(4)(1992) 290-300. [20] V. J. D. Tsai, Delaunay triangulations in TIN creation: an overview and a linear-time algorithm. International Journal of Geographical Information Systems, 7(6)(1993) 501-524. [21] S. W. Cheng, T. K. Dey, and J. R. Shewchuk. Delaunay Mesh Generation. CRC Press, 2012. [22] X. Wang, Study on an algorithm for fast constructing Delaunay triangulation and 3D Visualization in OpenGL Environment. Science Technology and Engineering, 11(9)(2011) 2070-2074. [23] N. Amenta, S. Choi, G. Rote, Incremental constructions con BRIO, in: Proceedings of the Nineteenth Annual Symposium on Computational Geometry, ACM, New York, NY, USA, 2003, pp. 211-219. [24] Y. Liu, J. Snoeyink, A comparison of five implementations of 3D Delaunay tessellation. Combinatorial and Computational Geometry, 52(2005) 439–458. [25] S. Zhou, C. B. Jones, HCPO: an efficient insertion order for incremental Delaunay triangulation. Information Processing Letters, 93(1)(2005) 37–42. [26] J. D. Boissonnat, O. Devillers, S. Hornus, Incremental construction of the Delaunay triangulation and the Delaunay graph in medium dimension, in: Proceedings of the 25th Annual Symposium on Computational Geometry, ACM, New York, NY, USA, 2009, pp. 208–216. [27] K. Buchin, Organizing point sets: space-filling curves, Delau-

586 587 588

[28]

589 590 591

[29]

592 593

[30]

594 595 596

[31]

597 598 599

[32]

600 601

[33]

602 603 604

[34]

605 606

[35]

607 608 609 610 611 612 613

[36]

nay tessellations of random point sets, and flow complexes, Ph.D. Thesis, Free University, Berlin, 2007. K. Buchin, Constructing Delaunay triangulations along spacefilling curves, in: Proc. 2nd Internat. Sympos. Voronoi Diagrams in Science and Engineering, Seoul, Korea, 2005, pp. 184–195. L. Devroye, J. Jabbour, C. Zamora-Cura, Squarish k-d trees. SIAM Journal on Computing, 30(5)(2000) 1678–1700. L. Devroye, C. Lemaire, J. M. Moreau, Expected time analysis for Delaunay point location. Computational Geometry, 29(2)(2004) 61–89. S. H. Lo, Delaunay triangulation of non-uniform point distributions by means of multi-grid insertion. Finite Elements in Analysis and Design, 63(2013) 8–22. D. Hilbert, Ueber die stetige Abbildung einer Line auf ein Fl¨achenst¨uck. Mathematische Annalen, 38(3)(1891) 459-460. S. F. Frisken, R. N. Perry, Simple and efficient traversal methods for quadtrees and octrees. Journal of Graphics Tools, 7(3)(2002) 1-11. Peitgen, H. O. and Saupe, D. (Eds.). The Science of Fractal Images. New York: Springer-Verlag, pp. 278 and 284, 1988. J. A. Quinn, F. C. Langbein, R. R. Martin, G. Elber, DensityControlled Sampling of Parametric Surfaces Using Adaptive Space-Filling Curves, In: Proceedings of Geometric Modeling and Processing, No. 4077 in LNCS, Springer, 2006, pp. 658678. K. Sugihara, M. IRI, Construction of the Voronoi diagram for “one million” generators in single-precision arithmetic, Proceedings of the IEEE, 80(9)(1992) 1471-1484.

11