GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_distance.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 47 63 74.6%
Functions: 3 4 75.0%
Branches: 34 84 40.5%
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/mesh/mesh_distance.h>
41 #include <geogram/mesh/mesh_AABB.h>
42 #include <geogram/mesh/mesh_repair.h>
43 #include <geogram/mesh/mesh_sampling.h>
44 #include <geogram/basic/process.h>
45 #include <geogram/basic/stopwatch.h>
46
47 #include <algorithm>
48
49 namespace {
50
51 using namespace GEO;
52
53 /**
54 * \brief Computes largest distance between
55 * an interval of an array of points and
56 * a surface stored in a MeshFacetsAABB.
57 * \details Used by compute_max_distance() that
58 * runs multiple instances of this class
59 * (one per core).
60 */
61 class DistanceThread : public Thread {
62 public:
63 /**
64 * \brief Constructs a new DistanceThread
65 * \details This DistanceThread will compute the distances
66 * between a batch of points and an axis-aligned bounding box
67 * tree.
68 * \param[in] AABB The axis-aligned bounding box tree used to compute
69 * the distance
70 * \param[in] from first point index
71 * \param[in] to one position past the last point index
72 * \param[in] points_ptr pointer to the points
73 * \param[in] points_stride number of doubles between two
74 * consecutive points
75 */
76 DistanceThread(
77 const MeshFacetsAABB& AABB,
78 index_t from, index_t to,
79 const double* points_ptr, index_t points_stride = 3
80 128 ) :
81 128 AABB_(AABB),
82 128 from_(from),
83 128 to_(to),
84 128 max_sq_dist_(0.0),
85 128 points_ptr_(points_ptr),
86
1/2
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
 128 points_stride_(points_stride) {
87 }
88
89 /**
90 * \brief Runs the thread.
91 * \details Computes the distances for the batch of points associated
92 * with this thread.
93 */
94 128 void run() override {
95
2/2
✓ Branch 0 taken 800714 times.
✓ Branch 1 taken 128 times.
 800842 for(index_t v = from_; v < to_; v++) {
96 // TODO: optimization
97 // if we know that the points are spatially
98 // sorted, then we can use AABB_.squared_distance_with_hint()
99 800714 double sq_dist = AABB_.squared_distance(
100 *reinterpret_cast<const vec3*>(
101 800714 points_ptr_ + v * points_stride_
102 )
103 );
104 800714 max_sq_dist_ = std::max(
105 max_sq_dist_, sq_dist
106 );
107 }
108 128 }
109
110 /**
111 * \brief Gets the computed max squared distance.
112 * \return the maximum squared distance computed so far
113 */
114 double max_squared_distance() const {
115 128 return max_sq_dist_;
116 }
117
118 private:
119 const MeshFacetsAABB& AABB_;
120 index_t from_;
121 index_t to_;
122 double max_sq_dist_;
123 const double* points_ptr_;
124 index_t points_stride_;
125 };
126
127 /**
128 * \brief Computes largest distance between
129 * an array of points and a mesh stored in
130 * a MeshFacetsAABB.
131 * \details Uses a multithread implementation.
132 * \param[out] result the maximum squared distance
133 * \param[in] AABB the mesh stored in an axis-aligned bounding box
134 * \param[in] nb_points number of query points
135 * \param[in] points_ptr pointer to the points coordinates
136 * \param[in] points_stride number of doubles between two consecutive
137 * points
138 */
139 32 void compute_max_distance(
140 double& result,
141 const MeshFacetsAABB& AABB,
142 index_t nb_points,
143 const double* points_ptr,
144 index_t points_stride = 3
145 ) {
146 32 Stopwatch W;
147 TypedThreadGroup<DistanceThread> threads;
148
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 index_t nb_threads = Process::maximum_concurrent_threads();
149 32 index_t batch_size = nb_points / nb_threads;
150 index_t cur = 0;
151 index_t remaining = nb_points;
152
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 32 times.
 160 for(index_t i = 0; i < nb_threads; i++) {
153 index_t this_batch_size = batch_size;
154
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 96 times.
 128 if(i == nb_threads - 1) {
155 this_batch_size = remaining;
156 }
157 threads.push_back(
158 128 new DistanceThread(
159 AABB,
160 cur, cur + this_batch_size,
161 points_ptr, points_stride
162
2/4
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 128 times.
✗ Branch 5 not taken.
 128 )
163 );
164 cur += this_batch_size;
165 128 remaining -= this_batch_size;
166 }
167
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 32 geo_assert(remaining == 0);
168
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 Process::run_threads(threads);
169
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 32 times.
 160 for(index_t t = 0; t < threads.size(); t++) {
170
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 83 times.
 173 result = std::max(result, threads[t]->max_squared_distance());
171 }
172
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 double elapsed = W.elapsed_time();
173
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 27 times.
 32 if(elapsed == 0.0) {
174
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 Logger::out("AABB")
175 << "???? Mqueries / s (too fast to be measured)"
176 << std::endl;
177 } else {
178
2/4
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
 54 Logger::out("AABB")
179
1/2
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
 27 << double(nb_points) / (1.0e6 * elapsed)
180 << " Mqueries / s"
181 << std::endl;
182 }
183 32 }
184 }
185
186 /****************************************************************************/
187
188 namespace GEO {
189
190 16 double mesh_one_sided_Hausdorff_distance(
191 Mesh& m1, Mesh& m2, double sampling_step
192 ) {
193 16 double result = 0.0;
194 16 MeshFacetsAABB AABB(m2);
195
196 index_t nb_points = 0;
197
198
2/4
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 if(m1.cells.nb() == 0 && m1.edges.nb() == 0) {
199 nb_points = m1.vertices.nb();
200
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 compute_max_distance(
201 result, AABB, m1.vertices.nb(),
202 m1.vertices.point_ptr(0), m1.vertices.dimension()
203 );
204 } else {
205 // If the mesh has cells, then we need to remove the vertices
206 // that are not on the surface, else their distance to the
207 // other surface will be included in the computation !
208 ✗ vector<bool> on_surface(m1.vertices.nb(),false);
209 ✗ for(index_t f: m1.facets) {
210 ✗ for(index_t v: m1.facets.vertices(f)) {
211 on_surface[v] = true;
212 }
213 }
214 ✗ for(index_t v: m1.vertices) {
215 ✗ if(on_surface[v]) {
216 ✗ ++nb_points;
217 }
218 }
219 vector<double> points;
220 ✗ points.reserve(nb_points*m1.vertices.dimension());
221 ✗ for(index_t v: m1.vertices) {
222 ✗ if(on_surface[v]) {
223 ✗ for(index_t c=0; c<m1.vertices.dimension(); ++c) {
224 ✗ points.push_back(m1.vertices.point_ptr(v)[c]);
225 }
226 }
227 }
228
229 ✗ compute_max_distance(
230 result, AABB, nb_points, points.data(),
231 m1.vertices.dimension()
232 );
233 }
234
235 index_t nb_samples = index_t(
236
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 Geom::mesh_area(m1, 3) / geo_sqr(sampling_step)
237 16 );
238
239
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if(nb_samples > nb_points) {
240
241 16 nb_samples -= nb_points;
242
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 Logger::out("Distance") << "Using " << nb_samples
243 << " additional samples"
244 << std::endl;
245
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 vector<double> samples(nb_samples * 3);
246 Attribute<double> weight; // left unbound
247
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 mesh_generate_random_samples_on_surface<3>(
248 m1, samples.data(), nb_samples, weight
249 );
250
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 compute_max_distance(
251 result, AABB, nb_samples, samples.data(), 3
252 );
253 }
254
255 16 return ::sqrt(result);
256 }
257
258 ✗ double mesh_symmetric_Hausdorff_distance(
259 Mesh& m1, Mesh& m2, double sampling_step
260 ) {
261 return std::max(
262 ✗ mesh_one_sided_Hausdorff_distance(m1, m2, sampling_step),
263 ✗ mesh_one_sided_Hausdorff_distance(m2, m1, sampling_step)
264 ✗ );
265 }
266 }
267