GCC Code Coverage Report

Directory: ./
File: delaunay/LFS.cpp
Date: 2026-09-27 03:12:47
Exec Total Coverage
Lines: 96 97 99.0%
Functions: 3 3 100.0%
Branches: 107 162 66.0%
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/LFS.h>
41 #include <geogram/basic/geometry.h>
42 #include <geogram/basic/logger.h>
43 #include <geogram/basic/command_line.h>
44
45 namespace {
46
47 using namespace GEO;
48
49 const index_t pid_size = 6;
50 int pid[pid_size][4] = {
51 {0, 1, 2, 3},
52 {0, 2, 1, 3},
53 {0, 3, 1, 2},
54 {1, 2, 0, 3},
55 {1, 3, 2, 0},
56 {2, 3, 0, 1}
57 };
58
59 /**
60 * \brief Computes the circumcenter and squared radius of a tetrahedron.
61 * \param[in] p first vertex of the tetrahedron
62 * \param[in] q second vertex of the tetrahedron
63 * \param[in] r third vertex of the tetrahedron
64 * \param[in] s fourth vertex of the tetrahedron
65 * \param[out] t_circumcenter computed circumcenter of the tetrahedron
66 * \param[out] squared_radius computed squared radius of the tetrahedron
67 * \param[out] dihedral_angles if non-nullptr,
68 * computed dihedral angles of the tetrahedron
69 * \return true if successful, false otherwise
70 * (for instance, if the tetrahedron is flat).
71 */
72
1/2
✓ Branch 0 taken 55452 times.
✗ Branch 1 not taken.
 55452 bool tetra_circumcenter_squaredradius(
73 const vec3& p, const vec3& q,
74 const vec3& r, const vec3& s,
75 vec3& t_circumcenter, double& squared_radius,
76 double* dihedral_angles = nullptr
77 ) {
78
79 vec3 qp = q - p;
80 vec3 rp = r - p;
81 vec3 sp = s - p;
82
83 55452 double tet_vol_x_12 = 2.0 * dot(qp, cross(rp, sp));
84
85 // it is not a safe check, one should avoid
86 // to input a degenerated tetra.
87
1/2
✓ Branch 0 taken 55452 times.
✗ Branch 1 not taken.
 55452 if(tet_vol_x_12 == 0.0) {
88 return false;
89 }
90
91 vec3 d_t_circumcenter =
92 length2(qp) * cross(rp, sp) +
93 length2(rp) * cross(sp, qp) +
94 length2(sp) * cross(qp, rp);
95
96 55452 squared_radius = length2(
97 55452 d_t_circumcenter) / (tet_vol_x_12 * tet_vol_x_12
98 );
99 55452 d_t_circumcenter = (1.0 / tet_vol_x_12) * d_t_circumcenter;
100 55452 t_circumcenter = p + d_t_circumcenter;
101
102
1/2
✓ Branch 0 taken 55452 times.
✗ Branch 1 not taken.
 55452 if(dihedral_angles) {
103 // compute Dihedral angle directly
104 // the following computation is not optimized.
105
2/2
✓ Branch 0 taken 221808 times.
✓ Branch 1 taken 55452 times.
 277260 vec3 point[4];
106 55452 point[0] = p;
107 55452 point[1] = q;
108 55452 point[2] = r;
109 55452 point[3] = s;
110
2/2
✓ Branch 0 taken 332712 times.
✓ Branch 1 taken 55452 times.
 388164 for(index_t i = 0; i < pid_size; i++) {
111 332712 vec3 b1 = point[pid[i][0]] - point[pid[i][2]];
112 332712 vec3 b2 = point[pid[i][1]] - point[pid[i][0]];
113 332712 vec3 b3 = point[pid[i][3]] - point[pid[i][1]];
114 vec3 b2b3 = cross(b2, b3);
115 332712 dihedral_angles[i] = ::fabs(
116 332712 ::atan2(
117 length(b2) * dot(b1, b2b3),
118 dot(cross(b1, b2), b2b3)
119 )
120 );
121 }
122 }
123 return true;
124 }
125
126 /**
127 * \brief Gets a Delaunay vertex by global index.
128 * \param[in] delaunay the Delaunay triangulation
129 * \param[in] v the index of the vertex
130 * \return a const reference to the vertex, as a vec3
131 */
132 inline const vec3& delaunay_vertex(Delaunay* delaunay, index_t v) {
133 return *(const vec3*) delaunay->vertex_ptr(v);
134 }
135
136 /**
137 * \brief Gets a Delaunay vertex by tetrahedron index
138 * and local vertex index.
139 * \param[in] delaunay the Delaunay triangulation
140 * \param[in] c the index of the tetrahedron
141 * \param[in] lv the local index of the vertex (0,1,2 or 3)
142 * in tetrahedron \p c.
143 * \return a const reference to the vertex, as a vec3.
144 */
145 inline const vec3& delaunay_tet_vertex(
146 Delaunay* delaunay, index_t c, index_t lv
147 ) {
148 return delaunay_vertex(delaunay, index_t(delaunay->cell_vertex(c, lv)));
149 }
150
151 index_t tet_facet_vertex[4][3] = {
152 {1, 2, 3},
153 {0, 3, 2},
154 {3, 0, 1},
155 {2, 1, 0}
156 };
157
158 /**
159 * \brief Computes the normal to a tetrahedron facet.
160 * \param[in] delaunay the Delaunay triangulation
161 * \param[in] t the index of the tetrahedron
162 * \param[in] f the local index (0,1,2 or 3) of
163 * the facet in the tetrahedron \p t
164 * \return the normal to the facet \p f or tetrahedron \p t
165 */
166 6249 vec3 delaunay_facet_normal(
167 Delaunay* delaunay, index_t t, index_t f
168 ) {
169 geo_debug_assert(f < 4);
170 6249 index_t v1 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][0]));
171 6249 index_t v2 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][1]));
172 6249 index_t v3 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][2]));
173 const vec3& p1 = delaunay_vertex(delaunay, v1);
174 const vec3& p2 = delaunay_vertex(delaunay, v2);
175 const vec3& p3 = delaunay_vertex(delaunay, v3);
176 6249 return cross(p2 - p1, p3 - p1);
177 }
178 }
179
180 /****************************************************************************/
181
182 namespace GEO {
183
184 1 void LocalFeatureSize::init(index_t nb_pts, const double* pts) {
185 // Note: I need PDEL here instead of ANN/BNN since I need
186 // to access the cells (and ANN/BNN only give me the
187 // neighbors)
188
189
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Logger::out("LFS") << "Delaunay" << std::endl;
190
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Delaunay_var delaunay = Delaunay::create(3, "PDEL");
191
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 delaunay->set_vertices(nb_pts, pts);
192
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1 Logger::out("LFS") << "Done Delaunay" << std::endl;
193
194
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector<vec3> circumcenter(delaunay->nb_cells());
195 vector<bool> voronoi_cell_is_infinite(
196 delaunay->nb_vertices(), false
197
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 );
198
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector<bool> is_sliver(delaunay->nb_cells());
199
200 static constexpr index_t NO_POLE = NO_INDEX;
201 // Index of incident tet whose circumcenter is the pole
202
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector<index_t> positive_pole(delaunay->nb_vertices(),NO_POLE);
203
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector<index_t> negative_pole(delaunay->nb_vertices(),NO_POLE);
204
205 vector<vec3> avg_infinite_dir(
206 ✗ delaunay->nb_vertices(), vec3(0, 0, 0)
207
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 );
208
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1 vector<double> dist(delaunay->nb_vertices(), 0.0);
209
210
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Logger::out("LFS") << "(1) Circumcenters and slivers" << std::endl;
211 // Step 1: compute circumcenters and check for slivers
212 1 const double sliver_quality = sliver_angle_threshold_ / 180.0 * M_PI;
213
3/4
✓ Branch 1 taken 55453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55452 times.
✓ Branch 4 taken 1 times.
 55453 for(index_t t = 0; t < delaunay->nb_cells(); t++) {
214 const vec3& p = delaunay_tet_vertex(delaunay, t, 0);
215 const vec3& q = delaunay_tet_vertex(delaunay, t, 1);
216 const vec3& r = delaunay_tet_vertex(delaunay, t, 2);
217 const vec3& s = delaunay_tet_vertex(delaunay, t, 3);
218
219 double dihedral_angle[6];
220 double tet_radius;
221 55452 bool ok = tetra_circumcenter_squaredradius(
222 p, q, r, s,
223 circumcenter[t], tet_radius, dihedral_angle
224 );
225
226
2/2
✓ Branch 0 taken 332454 times.
✓ Branch 1 taken 55376 times.
 387830 for(index_t a = 0; a < 6; a++) {
227
2/2
✓ Branch 0 taken 332378 times.
✓ Branch 1 taken 76 times.
 332454 if(!ok) {
228 break;
229 }
230 332378 ok = ok &&
231
2/2
✓ Branch 0 taken 332298 times.
✓ Branch 1 taken 80 times.
 332378 dihedral_angle[a] >= sliver_quality &&
232
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 332291 times.
 332298 dihedral_angle[a] <= M_PI - sliver_quality;
233 }
234 is_sliver[t] = !ok;
235 }
236
237
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Logger::out("LFS") << "(2) Positive poles" << std::endl;
238 // Step 2: compute positive poles
239
3/4
✓ Branch 1 taken 55453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55452 times.
✓ Branch 4 taken 1 times.
 55453 for(index_t t = 0; t < delaunay->nb_cells(); t++) {
240
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 55365 times.
 55452 if(is_sliver[t]) {
241 87 continue;
242 }
243 index_t f_inf = NO_INDEX;
244
2/2
✓ Branch 0 taken 209300 times.
✓ Branch 1 taken 49116 times.
 258416 for(index_t f = 0; f < 4; f++) {
245
3/4
✓ Branch 1 taken 209300 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 203051 times.
✓ Branch 4 taken 6249 times.
 209300 if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
246 f_inf = index_t(f);
247 break;
248 }
249 }
250
2/2
✓ Branch 0 taken 49116 times.
✓ Branch 1 taken 6249 times.
 55365 if(f_inf == NO_INDEX) {
251 // tet t does not have facet on border
252 const vec3& c = circumcenter[t];
253
2/2
✓ Branch 0 taken 196464 times.
✓ Branch 1 taken 49116 times.
 245580 for(index_t lv = 0; lv < 4; lv++) {
254
1/2
✓ Branch 1 taken 196464 times.
✗ Branch 2 not taken.
 196464 index_t iv = index_t(delaunay->cell_vertex(t, lv));
255 const vec3& v = delaunay_vertex(delaunay, iv);
256 double d = length2(c - v);
257
2/2
✓ Branch 0 taken 34757 times.
✓ Branch 1 taken 161707 times.
 196464 if(d > dist[iv]) {
258 34757 dist[iv] = d;
259 34757 positive_pole[iv] = t;
260 }
261 }
262 } else {
263 // tet t has facet f_inf on border
264 // positive pole = average direction of infinite Voronoi edges
265
1/2
✓ Branch 1 taken 6249 times.
✗ Branch 2 not taken.
 6249 index_t v = index_t(delaunay->cell_vertex(t, index_t(f_inf)));
266 6249 voronoi_cell_is_infinite[v] = true;
267 avg_infinite_dir[v] +=
268 6249 delaunay_facet_normal(delaunay, t, index_t(f_inf));
269 }
270 }
271
272
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Logger::out("LFS") << "(3) Negative poles" << std::endl;
273 // Step 3: compute negative poles
274 std::fill(dist.begin(), dist.end(), 0.0);
275
3/4
✓ Branch 1 taken 55453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55452 times.
✓ Branch 4 taken 1 times.
 55453 for(index_t t = 0; t < delaunay->nb_cells(); t++) {
276
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 55365 times.
 55452 if(is_sliver[t]) {
277 87 continue;
278 }
279 bool t_is_infinite = false;
280
2/2
✓ Branch 0 taken 209300 times.
✓ Branch 1 taken 49116 times.
 258416 for(index_t f = 0; f < 4; f++) {
281
3/4
✓ Branch 1 taken 209300 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 203051 times.
✓ Branch 4 taken 6249 times.
 209300 if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
282 t_is_infinite = true;
283 break;
284 }
285 }
286
2/2
✓ Branch 0 taken 6249 times.
✓ Branch 1 taken 49116 times.
 55365 if(t_is_infinite) {
287 6249 continue;
288 }
289 const vec3& c = circumcenter[t];
290
2/2
✓ Branch 0 taken 196464 times.
✓ Branch 1 taken 49116 times.
 245580 for(index_t lv = 0; lv < 4; lv++) {
291
1/2
✓ Branch 1 taken 196464 times.
✗ Branch 2 not taken.
 196464 index_t iv = index_t(delaunay->cell_vertex(t, lv));
292 const vec3& v = delaunay_vertex(delaunay, iv);
293 vec3 N;
294
2/2
✓ Branch 0 taken 40326 times.
✓ Branch 1 taken 156138 times.
 196464 if(voronoi_cell_is_infinite[iv]) {
295 N = -avg_infinite_dir[iv];
296 } else {
297 156138 N = v - circumcenter[positive_pole[iv]];
298 }
299 double d = dot(c - v, N);
300
2/2
✓ Branch 0 taken 25664 times.
✓ Branch 1 taken 170800 times.
 196464 if(d > dist[iv]) {
301 25664 dist[iv] = d;
302 25664 negative_pole[iv] = t;
303 }
304 }
305 }
306
307
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1 Logger::out("LFS") << "(4) Kd-tree" << std::endl;
308 // Step 4: create search structure
309
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector<bool> is_pole(delaunay->nb_cells(), false);
310
3/4
✓ Branch 1 taken 10001 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10000 times.
✓ Branch 4 taken 1 times.
 10001 for(index_t iv = 0; iv < delaunay->nb_vertices(); iv++) {
311
2/2
✓ Branch 0 taken 8475 times.
✓ Branch 1 taken 1525 times.
 10000 if(negative_pole[iv] != NO_POLE) {
312 8475 is_pole[index_t(negative_pole[iv])] = true;
313 }
314
2/2
✓ Branch 0 taken 7107 times.
✓ Branch 1 taken 2893 times.
 10000 if(!voronoi_cell_is_infinite[iv]) {
315
2/2
✓ Branch 0 taken 6899 times.
✓ Branch 1 taken 208 times.
 7107 if(positive_pole[iv] != NO_POLE) {
316 6899 is_pole[index_t(positive_pole[iv])] = true;
317 }
318 }
319 }
320
321 index_t nb_poles = 0;
322
3/4
✓ Branch 1 taken 55453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55452 times.
✓ Branch 4 taken 1 times.
 55453 for(index_t t = 0; t < delaunay->nb_cells(); t++) {
323
2/2
✓ Branch 0 taken 15155 times.
✓ Branch 1 taken 40297 times.
 55452 if(is_pole[t]) {
324 15155 nb_poles++;
325 }
326 }
327
328
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 poles_.reserve(nb_poles * 3);
329
3/4
✓ Branch 1 taken 55453 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55452 times.
✓ Branch 4 taken 1 times.
 55453 for(index_t t = 0; t < delaunay->nb_cells(); t++) {
330
2/2
✓ Branch 0 taken 15155 times.
✓ Branch 1 taken 40297 times.
 55452 if(is_pole[t]) {
331
1/2
✓ Branch 0 taken 15155 times.
✗ Branch 1 not taken.
 15155 poles_.push_back(circumcenter[t].x);
332
1/2
✓ Branch 0 taken 15155 times.
✗ Branch 1 not taken.
 15155 poles_.push_back(circumcenter[t].y);
333
1/2
✓ Branch 0 taken 15155 times.
✗ Branch 1 not taken.
 15155 poles_.push_back(circumcenter[t].z);
334 }
335 }
336
337
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 spatial_search_ = Delaunay::create(3, "NN");
338
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 spatial_search_->set_vertices(poles_.size() / 3, poles_.data());
339
340
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Logger::out("LFS") << "Done init." << std::endl;
341 1 }
342 }
343