GCC Code Coverage Report

Directory: ./
File: lib/geogram/delaunay/delaunay_2d.cpp
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 304 507 60.0%
Functions: 11 19 57.9%
Branches: 270 774 34.9%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/delaunay/delaunay_2d.h>
41 #include <geogram/basic/logger.h>
42 #include <geogram/basic/geometry_nd.h>
43 #include <geogram/basic/process.h>
44 #include <geogram/basic/command_line.h>
45 #include <geogram/basic/stopwatch.h>
46 #include <geogram/basic/matrix.h>
47 #include <geogram/basic/permutation.h>
48 #include <geogram/mesh/mesh_reorder.h>
49 #include <geogram/bibliography/bibliography.h>
50
51 #include <stack>
52 #include <algorithm>
53
54 // TODO: optimizations:
55 // - convex hull traversal for nearest_vertex()
56
57 namespace {
58 using namespace GEO;
59
60 /**
61 * \brief Computes the (approximate) orientation predicate in 2d.
62 * \details Computes the sign of the (approximate) signed volume of
63 * the triangle p0, p1, p2
64 * \param[in] p0 first vertex of the triangle
65 * \param[in] p1 second vertex of the triangle
66 * \param[in] p2 third vertex of the triangle
67 * \retval POSITIVE if the triangle is oriented positively
68 * \retval ZERO if the triangle is flat
69 * \retval NEGATIVE if the triangle is oriented negatively
70 * \todo check whether orientation is inverted as compared to
71 * Shewchuk's version.
72 */
73 1799 inline Sign orient_2d_inexact(
74 const double* p0, const double* p1, const double* p2
75 ) {
76 1799 double a11 = p1[0] - p0[0] ;
77 1799 double a12 = p1[1] - p0[1] ;
78
79 1799 double a21 = p2[0] - p0[0] ;
80 1799 double a22 = p2[1] - p0[1] ;
81
82 1799 double Delta = det2x2(
83 a11,a12,
84 a21,a22
85 1799 );
86
87
1/2
✓ Branch 1 taken 1799 times.
✗ Branch 2 not taken.
 3598 return geo_sgn(Delta);
88 }
89 }
90
91 namespace GEO {
92
93 // triangle edge vertex is such that the triangle
94 // formed with:
95 // vertex lv
96 // triangle_edge_vertex[lv][0]
97 // triangle_edge_vertex[lv][1]
98 // has the same orientation as the original triangle for
99 // any vertex lv.
100
101 char Delaunay2d::triangle_edge_vertex_[3][2] = {
102 {1,2},
103 {2,0},
104 {0,1}
105 };
106
107 2 Delaunay2d::Delaunay2d(coord_index_t dimension) :
108
1/2
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 Delaunay(dimension)
109 {
110
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
111 "DBLP:journals/cj/Bowyer81",
112 "One of the two initial references to the algorithm, "
113 "discovered independently and simultaneously by Bowyer and Watson."
114 );
115
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
116 "journals/cj/Watson81",
117 "One of the two initial references to the algorithm, "
118 "discovered independently and simultaneously by Bowyer and Watson."
119 );
120
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
121 "DBLP:conf/compgeom/AmentaCR03",
122 "Using spatial sorting has a dramatic impact on the performances."
123 );
124
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
125 "DBLP:journals/comgeo/FunkeMN05",
126 "Initializing \\verb|locate()| with a non-exact version "
127 " (structural filtering) gains (a bit of) performance."
128 );
129
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
130 "DBLP:journals/comgeo/BoissonnatDPTY02",
131 "The idea of traversing the cavity from inside "
132 " used in GEOGRAM is inspired by the implementation of "
133 " \\verb|Delaunay_triangulation_3| in CGAL."
134 );
135
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
136 "DBLP:conf/imr/Si06",
137 "The triangulation data structure used in GEOGRAM is inspired "
138 "by Tetgen."
139 );
140
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
141 "DBLP:journals/ijfcs/DevillersPT02",
142 "Analysis of the different versions of the line walk algorithm "
143 " used by \\verb|locate()|."
144 );
145
146
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if(dimension != 2 && dimension != 3) {
147 ✗ throw InvalidDimension(dimension, "Delaunay2d", "2 or 3");
148 }
149 2 first_free_ = END_OF_LIST;
150 2 weighted_ = (dimension == 3);
151 // In weighted mode, vertices are 3d but combinatorics is 2d.
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(weighted_) {
153 ✗ cell_size_ = 3;
154 ✗ cell_v_stride_ = 3;
155 ✗ cell_neigh_stride_ = 3;
156 }
157 2 cur_stamp_ = 0;
158
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay");
159
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose");
160
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 2 debug_mode_ = (debug_mode_ || verbose_debug_mode_);
161
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark");
162 2 has_empty_cells_ = false;
163 2 abort_if_empty_cell_ = false;
164 2 }
165
166 8 Delaunay2d::~Delaunay2d() {
167 8 }
168
169 2 void Delaunay2d::set_vertices(index_t nb_vertices, const double* vertices) {
170 2 has_empty_cells_ = false;
171
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 Stopwatch W("DelInternal", benchmark_mode_);
172
173 2 cur_stamp_ = 0;
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(weighted_) {
175 ✗ heights_.resize(nb_vertices);
176 ✗ for(index_t i = 0; i < nb_vertices; ++i) {
177 // Client code uses 3d embedding with ti = sqrt(W - wi)
178 // where W = max(wi)
179 // We recompute the standard "shifted" lifting on
180 // the paraboloid from it.
181 // (we use wi - W, everything is shifted by W, but
182 // we do not care since the power diagram is invariant
183 // by a translation of all weights).
184 ✗ double w = -geo_sqr(vertices[3 * i + 2]);
185 ✗ heights_[i] = -w +
186 ✗ geo_sqr(vertices[3 * i]) +
187 ✗ geo_sqr(vertices[3 * i + 1]);
188 }
189 }
190
191
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Delaunay::set_vertices(nb_vertices, vertices);
192
193 2 index_t expected_triangles = nb_vertices * 2;
194
195
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_v_store_.reserve(expected_triangles * 3);
196
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_cell_store_.reserve(expected_triangles * 3);
197
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_next_.reserve(expected_triangles);
198
199
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_v_store_.resize(0);
200
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_cell_store_.resize(0);
201
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_next_.resize(0);
202 2 first_free_ = END_OF_LIST;
203
204 // Sort the vertices spatially. This makes localisation
205 // faster.
206
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(do_reorder_) {
207
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 compute_BRIO_order(
208
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_vertices, vertex_ptr(0), reorder_, 2, dimension_
209 );
210 } else {
211 ✗ reorder_.resize(nb_vertices);
212 ✗ for(index_t i = 0; i < nb_vertices; ++i) {
213 ✗ reorder_[i] = i;
214 }
215 }
216
217 2 double sorting_time = 0;
218
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(benchmark_mode_) {
219
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 sorting_time = W.elapsed_time();
220
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 6 Logger::out("DelInternal1") << "BRIO sorting:"
221
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << sorting_time
222
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
223 }
224
225 // The indices of the vertices of the first triangle.
226 index_t v0, v1, v2;
227
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 if(!create_first_triangle(v0, v1, v2)) {
228 ✗ Logger::warn("Delaunay2d") << "All the points are colinear"
229 ✗ << std::endl;
230 ✗ return;
231 }
232
233 2 index_t hint = NO_TRIANGLE;
234 // Insert all the vertices incrementally.
235
2/2
✓ Branch 0 taken 229 times.
✓ Branch 1 taken 2 times.
 231 for(index_t i = 0; i < nb_vertices; ++i) {
236
1/2
✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
 229 index_t v = reorder_[i];
237 // Do not re-insert the first four vertices.
238
6/6
✓ Branch 0 taken 227 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 225 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 223 times.
✓ Branch 5 taken 2 times.
 229 if(v != v0 && v != v1 && v != v2) {
239
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 index_t new_hint = insert(v, hint);
240
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 127 times.
 223 if(new_hint == NO_TRIANGLE) {
241 96 has_empty_cells_ = true;
242
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 96 times.
 96 if(abort_if_empty_cell_) {
243 ✗ return;
244 }
245 } else {
246 127 hint = new_hint;
247 }
248 }
249 }
250
251
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(benchmark_mode_) {
252
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 6 Logger::out("DelInternal2") << "Core insertion algo:"
253
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 << W.elapsed_time() - sorting_time
254
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << std::endl;
255 }
256
257
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(debug_mode_) {
258 ✗ check_combinatorics(verbose_debug_mode_);
259 ✗ check_geometry(verbose_debug_mode_);
260 }
261
262 // Compress cell_to_v_store_ and cell_to_cell_store_
263 // (remove free and virtual tetrahedra).
264 // Since cell_next_ is not used at this point,
265 // we reuse it for storing the conversion array that
266 // maps old trgl indices to new trgl indices
267 // Note: trgl_is_real() uses the previous value of
268 // cell_next(), but we are processing indices
269 // in increasing order and since old2new[t] is always
270 // smaller or equal to t, we never overwrite a value
271 // before needing it.
272
273 2 vector<index_t>& old2new = cell_next_;
274 2 index_t nb_triangles = 0;
275 2 index_t nb_triangles_to_delete = 0;
276
277 {
278
3/4
✓ Branch 1 taken 280 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 278 times.
✓ Branch 4 taken 2 times.
 280 for(index_t t = 0; t < max_t(); ++t) {
279 278 if(
280
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 278 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 153 times.
✓ Branch 8 taken 125 times.
 556 (keep_infinite_ && !triangle_is_free(t)) ||
281
3/4
✓ Branch 1 taken 278 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 153 times.
✓ Branch 4 taken 125 times.
 278 triangle_is_real(t)
282 ) {
283
2/2
✓ Branch 0 taken 152 times.
✓ Branch 1 taken 1 times.
 153 if(t != nb_triangles) {
284 152 Memory::copy(
285 152 &cell_to_v_store_[nb_triangles * 3],
286
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
 152 &cell_to_v_store_[t * 3],
287 3 * sizeof(index_t)
288 );
289 152 Memory::copy(
290 152 &cell_to_cell_store_[nb_triangles * 3],
291
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
 152 &cell_to_cell_store_[t * 3],
292 3 * sizeof(index_t)
293 );
294 }
295
1/2
✓ Branch 1 taken 153 times.
✗ Branch 2 not taken.
 153 old2new[t] = nb_triangles;
296 153 ++nb_triangles;
297 } else {
298
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 125 old2new[t] = NO_INDEX;
299 125 ++nb_triangles_to_delete;
300 }
301 }
302
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_v_store_.resize(3 * nb_triangles);
303
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_cell_store_.resize(3 * nb_triangles);
304
2/2
✓ Branch 0 taken 459 times.
✓ Branch 1 taken 2 times.
 461 for(index_t i = 0; i < 3 * nb_triangles; ++i) {
305
1/2
✓ Branch 1 taken 459 times.
✗ Branch 2 not taken.
 459 index_t t = cell_to_cell_store_[i];
306
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 459 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 459 geo_debug_assert(t != NO_INDEX);
307
1/2
✓ Branch 1 taken 459 times.
✗ Branch 2 not taken.
 459 t = old2new[t];
308 // Note: t can be equal to -1 when a real trgl is
309 // adjacent to a virtual one (and this is how the
310 // rest of Vorpaline expects to see trgls on the
311 // border).
312
1/2
✓ Branch 1 taken 459 times.
✗ Branch 2 not taken.
 459 cell_to_cell_store_[i] = t;
313 }
314 }
315
316 // In "keep_infinite" mode, we reorder the cells in such
317 // a way that finite cells have indices [0..nb_finite_cells_-1]
318 // and infinite cells have indices [nb_finite_cells_ .. nb_cells_-1]
319
320
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(keep_infinite_) {
321 ✗ nb_finite_cells_ = 0;
322 ✗ index_t finite_ptr = 0;
323 ✗ index_t infinite_ptr = nb_triangles - 1;
324 for(;;) {
325 ✗ while(triangle_is_finite(finite_ptr)) {
326 ✗ old2new[finite_ptr] = finite_ptr;
327 ✗ ++finite_ptr;
328 ✗ ++nb_finite_cells_;
329 }
330 ✗ while(!triangle_is_finite(infinite_ptr)) {
331 ✗ old2new[infinite_ptr] = infinite_ptr;
332 ✗ --infinite_ptr;
333 }
334 ✗ if(finite_ptr > infinite_ptr) {
335 ✗ break;
336 }
337 ✗ old2new[finite_ptr] = infinite_ptr;
338 ✗ old2new[infinite_ptr] = finite_ptr;
339 ✗ ++nb_finite_cells_;
340 ✗ for(index_t lf=0; lf<3; ++lf) {
341 ✗ std::swap(
342 ✗ cell_to_cell_store_[3*finite_ptr + lf],
343 ✗ cell_to_cell_store_[3*infinite_ptr + lf]
344 );
345 }
346 ✗ for(index_t lv=0; lv<3; ++lv) {
347 ✗ std::swap(
348 ✗ cell_to_v_store_[3*finite_ptr + lv],
349 ✗ cell_to_v_store_[3*infinite_ptr + lv]
350 );
351 }
352 ✗ ++finite_ptr;
353 ✗ --infinite_ptr;
354 ✗ }
355 ✗ for(index_t i = 0; i < 3 * nb_triangles; ++i) {
356 ✗ index_t t = cell_to_cell_store_[i];
357 ✗ geo_debug_assert(t != NO_INDEX);
358 ✗ t = old2new[t];
359 ✗ geo_debug_assert(t != NO_INDEX);
360 ✗ cell_to_cell_store_[i] = t;
361 }
362 }
363
364
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(benchmark_mode_) {
365
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(keep_infinite_) {
366 ✗ Logger::out("DelCompress")
367 ✗ << "Removed " << nb_triangles_to_delete
368 ✗ << " triangles (free list)" << std::endl;
369 } else {
370
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 6 Logger::out("DelCompress")
371
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 << "Removed " << nb_triangles_to_delete
372
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 << " triangles (free list and infinite)" << std::endl;
373 }
374 }
375
376
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 set_arrays(
377 nb_triangles,
378 2 cell_to_v_store_.data(), cell_to_cell_store_.data()
379 );
380
381 // Not mandatory, but doing so makes it possible to
382 // use locate() in derived classes outside of
383 // set_vertices().
384
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 cell_next_.assign(cell_next_.size(),NO_INDEX);
385
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 }
386
387 ✗ index_t Delaunay2d::nearest_vertex(const double* p) const {
388
389 // TODO: For the moment, we fallback to the (unefficient)
390 // baseclass implementation when in weighted mode.
391 ✗ if(weighted_) {
392 ✗ return Delaunay::nearest_vertex(p);
393 }
394
395 // Find a triangle (real or virtual) that contains p
396 ✗ index_t t = locate(p, NO_TRIANGLE, thread_safe());
397
398 // If p is outside the convex hull of the inserted points,
399 // a special traversal is required (not implemented yet).
400 // TODO: implement convex hull boundary traversal
401 // (for now we fallback to linear search implemented
402 // in baseclass)
403 ✗ if(t == NO_TRIANGLE || triangle_is_virtual(t)) {
404 ✗ return Delaunay::nearest_vertex(p);
405 }
406
407 ✗ double sq_dist = 1e30;
408 ✗ index_t result = NO_TRIANGLE;
409
410 // Find the nearest vertex among t's vertices
411 ✗ for(index_t lv = 0; lv < 3; ++lv) {
412 ✗ index_t v = triangle_vertex(t, lv);
413 // If the triangle is virtual, then the first vertex
414 // is the vertex at infinity and is skipped.
415 ✗ if(v == NO_INDEX) {
416 ✗ continue;
417 }
418 ✗ double cur_sq_dist = Geom::distance2(p, vertex_ptr(v), 2);
419 ✗ if(cur_sq_dist < sq_dist) {
420 ✗ sq_dist = cur_sq_dist;
421 ✗ result = v;
422 }
423 }
424 ✗ return result;
425 }
426
427 223 index_t Delaunay2d::locate_inexact(
428 const double* p, index_t hint, index_t max_iter
429 ) const {
430
431 // If no hint specified, find a triangle randomly
432
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 223 times.
 225 while(hint == NO_TRIANGLE) {
433 2 hint = index_t(Numeric::random_int32()) % max_t();
434
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(triangle_is_free(hint)) {
435 ✗ hint = NO_TRIANGLE;
436 }
437 }
438
439
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 223 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 223 geo_debug_assert(!triangle_is_free(hint));
440
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 223 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 223 geo_debug_assert(!triangle_is_in_list(hint));
441
442 // Always start from a real trgl. If the trgl is virtual,
443 // find its real neighbor (always opposite to the
444 // infinite vertex)
445
2/2
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 153 times.
 223 if(triangle_is_virtual(hint)) {
446
1/2
✓ Branch 0 taken 144 times.
✗ Branch 1 not taken.
 144 for(index_t lf = 0; lf < 3; ++lf) {
447
2/2
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 74 times.
 144 if(triangle_vertex(hint, lf) == VERTEX_AT_INFINITY) {
448 70 hint = triangle_adjacent(hint, lf);
449
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 70 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 70 geo_debug_assert(hint != NO_TRIANGLE);
450 70 break;
451 }
452 }
453 }
454
455 223 index_t t = hint;
456 223 index_t t_pred = NO_TRIANGLE;
457
458 997 still_walking:
459 {
460 const double* pv[3];
461
2/4
✓ Branch 1 taken 997 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 997 times.
✗ Branch 5 not taken.
 997 pv[0] = vertex_ptr(finite_triangle_vertex(t,0));
462
2/4
✓ Branch 1 taken 997 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 997 times.
✗ Branch 5 not taken.
 997 pv[1] = vertex_ptr(finite_triangle_vertex(t,1));
463
2/4
✓ Branch 1 taken 997 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 997 times.
✗ Branch 5 not taken.
 997 pv[2] = vertex_ptr(finite_triangle_vertex(t,2));
464
465
2/2
✓ Branch 0 taken 2218 times.
✓ Branch 1 taken 188 times.
 2406 for(index_t le = 0; le < 3; ++le) {
466
467
1/2
✓ Branch 1 taken 2218 times.
✗ Branch 2 not taken.
 2218 index_t t_next = triangle_adjacent(t,le);
468
469 // If the opposite trgl is -1, then it means that
470 // we are trying to locate() (e.g. called from
471 // nearest_vertex) within a triangulation
472 // from which the infinite trgls were removed.
473
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2218 times.
 2218 if(t_next == NO_INDEX) {
474 35 return NO_TRIANGLE;
475 }
476
477 // If the candidate next triangle is the
478 // one we came from, then we know already that
479 // the orientation is positive, thus we examine
480 // the next candidate (or exit the loop if they
481 // are exhausted).
482
2/2
✓ Branch 0 taken 419 times.
✓ Branch 1 taken 1799 times.
 2218 if(t_next == t_pred) {
483 419 continue ;
484 }
485
486 // To test the orientation of p w.r.t. the facet f of
487 // t, we replace vertex number f with p in t (same
488 // convention as in CGAL).
489 1799 const double* pv_bkp = pv[le];
490 1799 pv[le] = p;
491
1/2
✓ Branch 1 taken 1799 times.
✗ Branch 2 not taken.
 1799 Sign ori = orient_2d_inexact(pv[0], pv[1], pv[2]);
492
493 // If the orientation is not negative, then we cannot
494 // walk towards t_next, and examine the next candidate
495 // (or exit the loop if they are exhausted).
496
2/2
✓ Branch 0 taken 990 times.
✓ Branch 1 taken 809 times.
 1799 if(ori != NEGATIVE) {
497 990 pv[le] = pv_bkp;
498 990 continue;
499 }
500
501 // If the opposite trgl is a virtual trgl, then
502 // the point has a positive orientation relative
503 // to the facet on the border of the convex hull,
504 // thus t_next is a trgl in conflict and we are
505 // done.
506
3/4
✓ Branch 1 taken 809 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35 times.
✓ Branch 4 taken 774 times.
 809 if(triangle_is_virtual(t_next)) {
507 35 return t_next;
508 }
509
510 // If we reach this point, then t_next is a valid
511 // successor, thus we are still walking.
512 774 t_pred = t;
513 774 t = t_next;
514
1/2
✓ Branch 0 taken 774 times.
✗ Branch 1 not taken.
 774 if(--max_iter != 0) {
515 774 goto still_walking;
516 }
517 }
518 }
519
520 // If we reach this point, we did not find a valid successor
521 // for walking (a face for which p has negative orientation),
522 // thus we reached the trgl for which p has all positive
523 // face orientations (i.e. the trgl that contains p).
524
525 188 return t;
526 }
527
528
529 223 index_t Delaunay2d::locate(
530 const double* p, index_t hint, bool thread_safe, Sign* orient
531 ) const {
532
533 // Try improving the hint by using the
534 // inexact locate function. This gains
535 // (a little bit) performance (a few
536 // percent in total Delaunay computation
537 // time), but it is better than nothing...
538 // Note: there is a maximum number of trgls
539 // traversed by locate_inexact() (2500)
540 // since there exists configurations in which
541 // locate_inexact() loops forever !
542
543
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 hint = locate_inexact(p, hint, 2500);
544
545 static Process::spinlock locate_lock = GEOGRAM_SPINLOCK_INIT;
546
547 // We need to have this spinlock because
548 // of random() that is not thread-safe
549 // (TODO: implement a random() function with
550 // thread local storage)
551
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 223 times.
 223 if(thread_safe) {
552 ✗ Process::acquire_spinlock(locate_lock);
553 }
554
555 // If no hint specified, find a triangle randomly
556
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 223 times.
 223 while(hint == NO_TRIANGLE) {
557 ✗ hint = index_t(Numeric::random_int32()) % max_t();
558 ✗ if(triangle_is_free(hint)) {
559 ✗ hint = NO_TRIANGLE;
560 }
561 }
562
563
2/8
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 223 geo_debug_assert(!triangle_is_free(hint));
564
2/8
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 223 geo_debug_assert(!triangle_is_in_list(hint));
565
566 // Always start from a real trgl. If the trgl is virtual,
567 // find its real neighbor (always opposite to the
568 // infinite vertex)
569
3/4
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35 times.
✓ Branch 4 taken 188 times.
 223 if(triangle_is_virtual(hint)) {
570
1/2
✓ Branch 0 taken 83 times.
✗ Branch 1 not taken.
 83 for(index_t le = 0; le < 3; ++le) {
571
3/4
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35 times.
✓ Branch 4 taken 48 times.
 83 if(triangle_vertex(hint, le) == VERTEX_AT_INFINITY) {
572
1/2
✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
 35 hint = triangle_adjacent(hint, le);
573
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 35 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 35 geo_debug_assert(hint != NO_TRIANGLE);
574 35 break;
575 }
576 }
577 }
578
579
2/8
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 223 geo_debug_assert(!triangle_is_free(hint));
580
2/8
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 223 geo_debug_assert(!triangle_is_in_list(hint));
581
2/8
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 223 geo_debug_assert(!triangle_is_virtual(hint));
582
583 223 index_t t = hint;
584 223 index_t t_pred = NO_TRIANGLE;
585 Sign orient_local[3];
586
1/2
✓ Branch 0 taken 223 times.
✗ Branch 1 not taken.
 223 if(orient == nullptr) {
587 ✗ orient = orient_local;
588 }
589
590
591 223 still_walking:
592 {
593 const double* pv[3];
594
2/4
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 5 not taken.
 223 pv[0] = vertex_ptr(finite_triangle_vertex(t,0));
595
2/4
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 5 not taken.
 223 pv[1] = vertex_ptr(finite_triangle_vertex(t,1));
596
2/4
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 223 times.
✗ Branch 5 not taken.
 223 pv[2] = vertex_ptr(finite_triangle_vertex(t,2));
597
598 // Start from a random facet
599
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 index_t e0 = index_t(Numeric::random_int32()) % 3;
600
2/2
✓ Branch 0 taken 640 times.
✓ Branch 1 taken 188 times.
 828 for(index_t de = 0; de < 3; ++de) {
601 640 index_t le = (e0 + de) % 3;
602
603
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 640 index_t t_next = triangle_adjacent(t,le);
604
605 // If the opposite triangle is -1, then it means that
606 // we are trying to locate() (e.g. called from
607 // nearest_vertex) within a triangulation
608 // from which the infinite trgls were removed.
609
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
 640 if(t_next == NO_INDEX) {
610 ✗ if(thread_safe) {
611 ✗ Process::release_spinlock(locate_lock);
612 }
613 35 return NO_TRIANGLE;
614 }
615
616
2/8
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 640 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 640 geo_debug_assert(!triangle_is_free(t_next));
617
2/8
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 640 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 640 geo_debug_assert(!triangle_is_in_list(t_next));
618
619 // If the candidate next triangle is the
620 // one we came from, then we know already that
621 // the orientation is positive, thus we examine
622 // the next candidate (or exit the loop if they
623 // are exhausted).
624
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
 640 if(t_next == t_pred) {
625 ✗ orient[le] = POSITIVE ;
626 ✗ continue ;
627 }
628
629 // To test the orientation of p w.r.t. the facet f of
630 // t, we replace vertex number f with p in t (same
631 // convention as in CGAL).
632 // This is equivalent to trgl_facet_point_orient3d(t,f,p)
633 // (but less costly, saves a couple of lookups)
634 640 const double* pv_bkp = pv[le];
635 640 pv[le] = p;
636
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 640 orient[le] = PCK::orient_2d(pv[0], pv[1], pv[2]);
637
638 // If the orientation is not negative, then we cannot
639 // walk towards t_next, and examine the next candidate
640 // (or exit the loop if they are exhausted).
641
2/2
✓ Branch 0 taken 605 times.
✓ Branch 1 taken 35 times.
 640 if(orient[le] != NEGATIVE) {
642 605 pv[le] = pv_bkp;
643 605 continue;
644 }
645
646 // If the opposite trgl is a virtual trgl, then
647 // the point has a positive orientation relative
648 // to the facet on the border of the convex hull,
649 // thus t_next is a trgl in conflict and we are
650 // done.
651
2/4
✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35 times.
✗ Branch 4 not taken.
 35 if(triangle_is_virtual(t_next)) {
652
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 35 times.
 35 if(thread_safe) {
653 ✗ Process::release_spinlock(locate_lock);
654 }
655
2/2
✓ Branch 0 taken 105 times.
✓ Branch 1 taken 35 times.
 140 for(index_t tle = 0; tle < 3; ++tle) {
656 105 orient[tle] = POSITIVE;
657 }
658 35 return t_next;
659 }
660
661 // If we reach this point, then t_next is a valid
662 // successor, thus we are still walking.
663 ✗ t_pred = t;
664 ✗ t = t_next;
665 ✗ goto still_walking;
666 }
667 }
668
669 // If we reach this point, we did not find a valid successor
670 // for walking (a face for which p has negative orientation),
671 // thus we reached the trgl for which p has all positive
672 // face orientations (i.e. the trgl that contains p).
673
674
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 188 times.
 188 if(thread_safe) {
675 ✗ Process::release_spinlock(locate_lock);
676 }
677 188 return t;
678 }
679
680 223 void Delaunay2d::find_conflict_zone(
681 index_t v,
682 index_t t, const Sign* orient,
683 index_t& t_bndry, index_t& e_bndry,
684 index_t& first, index_t& last
685 ) {
686 223 first = last = END_OF_LIST;
687
688 // Generate a unique stamp from current vertex index,
689 // used for marking triangles
690 223 set_triangle_mark_stamp(v);
691
692 // Pointer to the coordinates of the point to be inserted
693 223 const double* p = vertex_ptr(v);
694
695
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 223 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 223 geo_debug_assert(t != NO_TRIANGLE);
696
697 // Test whether the point already exists in
698 // the triangulation. The point already exists
699 // if it's located on three faces of the
700 // triangle returned by locate().
701 223 int nb_zero =
702 223 (orient[0] == ZERO) +
703 223 (orient[1] == ZERO) +
704 223 (orient[2] == ZERO) ;
705
706
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 127 times.
 223 if(nb_zero >= 2) {
707 96 return;
708 }
709
710 // Weighted triangulations can have dangling
711 // vertices. Such vertices p are characterized by
712 // the fact that p is not in conflict with the
713 // triangle returned by locate().
714
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 127 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 127 times.
 127 if(weighted_ && !triangle_is_conflict(t, p)) {
715 ✗ return;
716 }
717
718 // Note: points on edges and on facets are
719 // handled by the way triangle_is_in_conflict()
720 // is implemented, that naturally inserts
721 // the correct triangles in the conflict list.
722
723
724 // Mark t as conflict
725 127 add_triangle_to_list(t, first, last);
726
727 // A small optimization: if the point to be inserted
728 // is on some faces of the located triangle, insert
729 // the neighbors accros those edges in the conflict list.
730 // It saves a couple of calls to the predicates in this
731 // specific case (combinatorics are in general less
732 // expensive than the predicates).
733
3/4
✓ Branch 0 taken 127 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 86 times.
✓ Branch 3 taken 41 times.
 127 if(!weighted_ && nb_zero != 0) {
734
2/2
✓ Branch 0 taken 258 times.
✓ Branch 1 taken 86 times.
 344 for(index_t le = 0; le < 3; ++le) {
735
2/2
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 172 times.
 258 if(orient[le] == ZERO) {
736 86 index_t t2 = triangle_adjacent(t, le);
737 86 add_triangle_to_list(t2, first, last);
738 }
739 }
740
2/2
✓ Branch 0 taken 258 times.
✓ Branch 1 taken 86 times.
 344 for(index_t le = 0; le < 3; ++le) {
741
2/2
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 172 times.
 258 if(orient[le] == ZERO) {
742 86 index_t t2 = triangle_adjacent(t, le);
743 86 find_conflict_zone_iterative(
744 p,t2,t_bndry,e_bndry,first,last
745 );
746 }
747 }
748 }
749
750 // Determine the conflict list by greedy propagation from t.
751 127 find_conflict_zone_iterative(p,t,t_bndry,e_bndry,first,last);
752 }
753
754 213 void Delaunay2d::find_conflict_zone_iterative(
755 const double* p, index_t t_in,
756 index_t& t_bndry, index_t& e_bndry,
757 index_t& first, index_t& last
758 ) {
759
760 213 S_.push(t_in);
761
762
2/2
✓ Branch 1 taken 468 times.
✓ Branch 2 taken 213 times.
 681 while(!S_.empty()) {
763
764 468 index_t t = S_.top();
765 468 S_.pop();
766
767
2/2
✓ Branch 0 taken 1404 times.
✓ Branch 1 taken 468 times.
 1872 for(index_t le = 0; le < 3; ++le) {
768
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 index_t t2 = triangle_adjacent(t, le);
769
770 1846 if(
771
5/6
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 977 times.
✓ Branch 4 taken 427 times.
✓ Branch 5 taken 442 times.
✓ Branch 6 taken 962 times.
 2381 triangle_is_in_list(t2) || // known as conflict
772
3/4
✓ Branch 1 taken 977 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 15 times.
✓ Branch 4 taken 962 times.
 977 triangle_is_marked(t2) // known as non-conflict
773 ) {
774 697 continue;
775 }
776
777
3/4
✓ Branch 1 taken 962 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 255 times.
✓ Branch 4 taken 707 times.
 962 if(triangle_is_conflict(t2, p)) {
778 // Chain t2 in conflict list
779
1/2
✓ Branch 1 taken 255 times.
✗ Branch 2 not taken.
 255 add_triangle_to_list(t2, first, last);
780
1/2
✓ Branch 1 taken 255 times.
✗ Branch 2 not taken.
 255 S_.push(t2);
781 255 continue;
782 }
783
784 // At this point, t is in conflict
785 // and t2 is not in conflict.
786 // We keep a reference to a trgl on the boundary
787 707 t_bndry = t;
788 707 e_bndry = le;
789 // Mark t2 as visited (but not conflict)
790
1/2
✓ Branch 1 taken 707 times.
✗ Branch 2 not taken.
 707 mark_triangle(t2);
791 }
792 }
793 213 }
794
795 127 index_t Delaunay2d::stellate_conflict_zone(
796 index_t v_in, index_t t1, index_t t1ebord
797 ) {
798
799 127 index_t t = t1;
800 127 index_t e = t1ebord;
801 127 index_t t_adj = triangle_adjacent(t,e);
802
803
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 127 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 127 geo_debug_assert(t_adj != NO_INDEX);
804
805
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 127 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 127 geo_debug_assert(triangle_is_in_list(t));
806
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 127 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 127 geo_debug_assert(!triangle_is_in_list(t_adj));
807
808
809 127 index_t new_t_first = NO_INDEX;
810 127 index_t new_t_prev = NO_INDEX;
811
812 do {
813
814 722 index_t v1 = triangle_vertex(t, (e+1)%3);
815 722 index_t v2 = triangle_vertex(t, (e+2)%3);
816
817 // Create new triangle
818 722 index_t new_t = new_triangle(v_in, v1, v2);
819
820 // Connect new triangle to triangle on the other
821 // side of the conflict zone.
822 722 set_triangle_adjacent(new_t, 0, t_adj);
823 722 index_t adj_e = find_triangle_adjacent(t_adj, t);
824 722 set_triangle_adjacent(t_adj, adj_e, new_t);
825
826
827 // Move to next triangle
828 722 e = (e + 1)%3;
829 722 t_adj = triangle_adjacent(t,e);
830
2/2
✓ Branch 1 taken 682 times.
✓ Branch 2 taken 722 times.
 1404 while(triangle_is_in_list(t_adj)) {
831 682 t = t_adj;
832 682 e = (find_triangle_vertex(t,v2) + 2)%3;
833 682 t_adj = triangle_adjacent(t,e);
834
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 682 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 682 geo_debug_assert(t_adj != NO_INDEX);
835 }
836
837
2/2
✓ Branch 0 taken 127 times.
✓ Branch 1 taken 595 times.
 722 if(new_t_prev == NO_INDEX) {
838 127 new_t_first = new_t;
839 } else {
840 595 set_triangle_adjacent(new_t_prev, 1, new_t);
841 595 set_triangle_adjacent(new_t, 2, new_t_prev);
842 }
843
844 722 new_t_prev = new_t;
845
846
4/4
✓ Branch 0 taken 458 times.
✓ Branch 1 taken 264 times.
✓ Branch 2 taken 137 times.
✓ Branch 3 taken 127 times.
 722 } while((t != t1) || (e != t1ebord));
847
848 // Connect last triangle to first triangle
849 127 set_triangle_adjacent(new_t_prev, 1, new_t_first);
850 127 set_triangle_adjacent(new_t_first, 2, new_t_prev);
851
852 127 return new_t_prev;
853 }
854
855 223 index_t Delaunay2d::insert(index_t v, index_t hint) {
856 223 index_t t_bndry = NO_TRIANGLE;
857 223 index_t e_bndry = NO_INDEX;
858 223 index_t first_conflict = NO_TRIANGLE;
859 223 index_t last_conflict = NO_TRIANGLE;
860
861
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 const double* p = vertex_ptr(v);
862
863 Sign orient[3];
864
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 index_t t = locate(p, hint, false, orient);
865
1/2
✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
 223 find_conflict_zone(
866 v,t,orient,t_bndry,e_bndry,first_conflict,last_conflict
867 );
868
869 // The conflict list can be empty if:
870 // - Vertex v already exists in the triangulation
871 // - The triangulation is weighted and v is not visible
872
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 127 times.
 223 if(first_conflict == END_OF_LIST) {
873 96 return NO_TRIANGLE;
874 }
875
876
1/2
✓ Branch 1 taken 127 times.
✗ Branch 2 not taken.
 127 index_t new_triangle = stellate_conflict_zone(v,t_bndry,e_bndry);
877
878 // Recycle the triangles of the conflict zone.
879
1/2
✓ Branch 1 taken 127 times.
✗ Branch 2 not taken.
 127 cell_next_[last_conflict] = first_free_;
880 127 first_free_ = first_conflict;
881
882 // Return one of the newly created triangles
883 127 return new_triangle;
884 }
885
886 2 bool Delaunay2d::create_first_triangle(
887 index_t& iv0, index_t& iv1, index_t& iv2
888 ) {
889
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(nb_vertices() < 3) {
890 ✗ return false;
891 }
892
893 2 iv0 = 0;
894
895 2 iv1 = 1;
896 2 while(
897
3/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 2 times.
 6 iv1 < nb_vertices() &&
898
5/8
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 2 times.
 3 PCK::points_are_identical_2d(
899 vertex_ptr(iv0), vertex_ptr(iv1)
900 )
901 ) {
902 1 ++iv1;
903 }
904
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(iv1 == nb_vertices()) {
905 ✗ return false;
906 }
907
908 2 iv2 = iv1 + 1;
909 2 Sign s = ZERO;
910 2 while(
911
3/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 2 times.
 6 iv2 < nb_vertices() &&
912
6/10
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 2 times.
 3 (s = PCK::orient_2d(
913 vertex_ptr(iv0), vertex_ptr(iv1), vertex_ptr(iv2)
914 )) == ZERO
915 ) {
916 1 ++iv2;
917 }
918
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(iv2 == nb_vertices()) {
919 ✗ return false;
920 }
921
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(s == NEGATIVE) {
922 2 std::swap(iv1,iv2);
923 }
924
925 // Create the first triangle
926
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 index_t t0 = new_triangle(iv0, iv1, iv2);
927
928 // Create the first three virtual triangles surrounding it
929 index_t t[3];
930
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 for(index_t e = 0; e < 3; ++e) {
931 // In reverse order since it is an adjacent triangle
932
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 index_t v1 = triangle_vertex(t0, triangle_edge_vertex(e,1));
933
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 index_t v2 = triangle_vertex(t0, triangle_edge_vertex(e,0));
934
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 t[e] = new_triangle(VERTEX_AT_INFINITY, v1, v2);
935 }
936
937 // Connect the virtual triangles to the real one
938
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 for(index_t e=0; e<3; ++e) {
939
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 set_triangle_adjacent(t[e], 0, t0);
940
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 set_triangle_adjacent(t0, e, t[e]);
941 }
942
943 // Interconnect the three virtual triangles along their common
944 // edges
945
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 for(index_t e = 0; e < 3; ++e) {
946 // In reverse order since it is an adjacent triangle
947
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 index_t lv1 = triangle_edge_vertex(e,1);
948
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 index_t lv2 = triangle_edge_vertex(e,0);
949
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 set_triangle_adjacent(t[e], 1, t[lv1]);
950
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 set_triangle_adjacent(t[e], 2, t[lv2]);
951 }
952
953 2 return true;
954 }
955
956 /************************************************************************/
957
958 ✗ void Delaunay2d::show_triangle(index_t t) const {
959 std::cerr << "tri"
960 ✗ << (triangle_is_in_list(t) ? '*' : ' ')
961 ✗ << t
962 ✗ << ", v=["
963 ✗ << triangle_vertex(t, 0)
964 ✗ << ' '
965 ✗ << triangle_vertex(t, 1)
966 ✗ << ' '
967 ✗ << triangle_vertex(t, 2)
968 ✗ << "] adj=[";
969 ✗ show_triangle_adjacent(t, 0);
970 ✗ show_triangle_adjacent(t, 1);
971 ✗ show_triangle_adjacent(t, 2);
972 ✗ std::cerr << "] ";
973
974 ✗ for(index_t e = 0; e < 3; ++e) {
975 ✗ std::cerr << 'e' << e << ':';
976 ✗ for(index_t v = 0; v < 2; ++v) {
977 ✗ std::cerr << triangle_vertex(t, triangle_edge_vertex(e,v))
978 ✗ << ',';
979 }
980 ✗ std::cerr << ' ';
981 }
982 ✗ std::cerr << std::endl;
983 ✗ }
984
985 ✗ void Delaunay2d::show_triangle_adjacent(index_t t, index_t le) const {
986 ✗ index_t adj = triangle_adjacent(t, le);
987 ✗ if(adj != NO_INDEX) {
988 ✗ std::cerr << (triangle_is_in_list(adj) ? '*' : ' ');
989 }
990 ✗ std::cerr << adj;
991 ✗ std::cerr << ' ';
992 ✗ }
993
994 ✗ void Delaunay2d::show_list(
995 index_t first, const std::string& list_name
996 ) const {
997 ✗ index_t t = first;
998 ✗ std::cerr << "tri list: " << list_name << std::endl;
999 ✗ while(t != END_OF_LIST) {
1000 ✗ show_triangle(t);
1001 ✗ t = triangle_next(t);
1002 }
1003 ✗ std::cerr << "-------------" << std::endl;
1004 ✗ }
1005
1006 ✗ void Delaunay2d::check_combinatorics(bool verbose) const {
1007 ✗ if(verbose) {
1008 ✗ std::cerr << std::endl;
1009 }
1010 ✗ bool ok = true;
1011 ✗ std::vector<bool> v_has_triangle(nb_vertices(), false);
1012 ✗ for(index_t t = 0; t < max_t(); ++t) {
1013 ✗ if(triangle_is_free(t)) {
1014 /*
1015 if(verbose) {
1016 std::cerr << "-Deleted tri: ";
1017 show_tri(t);
1018 }
1019 */
1020 } else {
1021 /*
1022 if(verbose) {
1023 std::cerr << "Checking tri: ";
1024 show_tri(t);
1025 }
1026 */
1027 ✗ for(index_t le = 0; le < 3; ++le) {
1028 ✗ if(triangle_adjacent(t, le) == NO_INDEX) {
1029 ✗ std::cerr << le << ":Missing adjacent tri"
1030 ✗ << std::endl;
1031 ✗ ok = false;
1032 ✗ } else if(triangle_adjacent(t, le) == t) {
1033 ✗ std::cerr << le << ":Tri is adjacent to itself"
1034 ✗ << std::endl;
1035 ✗ ok = false;
1036 } else {
1037 ✗ index_t t2 = triangle_adjacent(t, le);
1038 ✗ bool found = false;
1039 ✗ for(index_t le2 = 0; le2 < 3; ++le2) {
1040 ✗ if(triangle_adjacent(t2, le2) == t) {
1041 ✗ found = true;
1042 }
1043 }
1044 ✗ if(!found) {
1045 std::cerr
1046 ✗ << le << ":Adjacent link is not bidirectional"
1047 ✗ << std::endl;
1048 ✗ ok = false;
1049 }
1050 }
1051 }
1052 ✗ index_t nb_infinite = 0;
1053 ✗ for(index_t lv = 0; lv < 3; ++lv) {
1054 ✗ if(triangle_vertex(t, lv) == NO_INDEX) {
1055 ✗ ++nb_infinite;
1056 }
1057 }
1058 ✗ if(nb_infinite > 1) {
1059 ✗ ok = false;
1060 ✗ std::cerr << "More than one infinite vertex"
1061 ✗ << std::endl;
1062 }
1063 }
1064 ✗ for(index_t lv = 0; lv < 3; ++lv) {
1065 ✗ index_t v = triangle_vertex(t, lv);
1066 ✗ if(v != NO_INDEX) {
1067 ✗ v_has_triangle[v] = true;
1068 }
1069 }
1070 }
1071 ✗ for(index_t v = 0; v < nb_vertices(); ++v) {
1072 ✗ if(!v_has_triangle[v]) {
1073 ✗ if(verbose) {
1074 ✗ std::cerr << "Vertex " << v
1075 ✗ << " is isolated (duplicated ?)" << std::endl;
1076 }
1077 }
1078 }
1079 ✗ geo_assert(ok);
1080 ✗ if(verbose) {
1081 ✗ std::cerr << std::endl;
1082 }
1083 ✗ std::cerr << std::endl << "Delaunay Combi OK" << std::endl;
1084 ✗ }
1085
1086 ✗ void Delaunay2d::check_geometry(bool verbose) const {
1087 ✗ bool ok = true;
1088 ✗ for(index_t t = 0; t < max_t(); ++t) {
1089 ✗ if(!triangle_is_free(t)) {
1090 ✗ index_t v0 = triangle_vertex(t, 0);
1091 ✗ index_t v1 = triangle_vertex(t, 1);
1092 ✗ index_t v2 = triangle_vertex(t, 2);
1093 ✗ for(index_t v = 0; v < nb_vertices(); ++v) {
1094 ✗ if(v == v0 || v == v1 || v == v2) {
1095 ✗ continue;
1096 }
1097 ✗ if(triangle_is_conflict(t, vertex_ptr(v))) {
1098 ✗ ok = false;
1099 ✗ if(verbose) {
1100 ✗ std::cerr << "Tri " << t <<
1101 ✗ " is in conflict with vertex " << v
1102 ✗ << std::endl;
1103
1104 ✗ std::cerr << " offending tri: ";
1105 ✗ show_triangle(t);
1106 }
1107 }
1108 }
1109 }
1110 }
1111 ✗ geo_assert(ok);
1112 ✗ std::cerr << std::endl << "Delaunay Geo OK" << std::endl;
1113 ✗ }
1114
1115 /************************************************************************/
1116
1117 ✗ RegularWeightedDelaunay2d::RegularWeightedDelaunay2d(
1118 coord_index_t dimension
1119 ✗ ) :
1120 ✗ Delaunay2d(3)
1121 {
1122 ✗ if(dimension != 3) {
1123 ✗ throw InvalidDimension(dimension, "RegularWeightedDelaunay2d", "3");
1124 }
1125 ✗ }
1126
1127 ✗ RegularWeightedDelaunay2d::~RegularWeightedDelaunay2d() {
1128 ✗ }
1129 }
1130