GCC Code Coverage Report

Directory: ./
File: lib/geogram/points/colocate.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 92 103 89.3%
Functions: 11 11 100.0%
Branches: 68 122 55.7%
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/points/colocate.h>
41 #include <geogram/points/nn_search.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/algorithm.h>
46
47 namespace {
48
49 using namespace GEO;
50
51 /**
52 * \brief Implements the colocate() algorithm when a tolerance is used.
53 * \details Uses multiple threads to speedup computations.
54 */
55 class Colocate {
56 public:
57 /**
58 * \brief Creates a new Colocate object.
59 * \param[in] NN the nearest neighbors search data structure
60 * \param[out] old2new where to store the relation between old indices
61 * and colocated vertices
62 * \param[in] tolerance maximum distance for colocated vertices
63 */
64 41 Colocate(
65 NearestNeighborSearch* NN,
66 vector<index_t>& old2new,
67 double tolerance
68 41 ) :
69 41 NN_(NN),
70 41 old2new_(old2new),
71 41 sq_tolerance_(geo_sqr(tolerance)) {
72 41 }
73
74 /**
75 * \brief Returns the number of points.
76 */
77 32622 index_t nb_points() const {
78 32622 return NN_->nb_points();
79 }
80
81 /**
82 * \brief Finds the nearest neighbors of a given point.
83 * \param[in] i index of the query point
84 * \param[in] nb maximum number of neighbors
85 * \return true when all the neighbors nearer than
86 * tolerance have been found, false otherwise.
87 */
88 32622 bool find_nearest_neighbors(index_t i, index_t nb) {
89 // allocated on the stack, more multithread-friendly
90 // and no need to deallocate (and VC++ does not support
91 // int neighbors[nb] where nb is a variable)
92 32622 index_t* neighbors = (index_t*) alloca(sizeof(index_t) * nb);
93 32622 double* dist = (double*) alloca(sizeof(double) * nb);
94
95
1/2
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
 32622 NN_->get_nearest_neighbors(
96
1/2
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
 32622 nb, NN_->point_ptr(i), neighbors, dist
97 );
98
99 32622 index_t smallest = i;
100
1/2
✓ Branch 0 taken 68498 times.
✗ Branch 1 not taken.
 68498 for(index_t jj = 0; jj < nb; jj++) {
101
2/2
✓ Branch 0 taken 32622 times.
✓ Branch 1 taken 35876 times.
 68498 if(dist[jj] > sq_tolerance_) {
102
1/2
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
 32622 old2new_[i] = smallest;
103 32622 return true;
104 }
105 35876 smallest = std::min(smallest, neighbors[jj]);
106 }
107 ✗ old2new_[i] = smallest;
108 ✗ return false;
109 }
110
111 /**
112 * \brief Finds all the neighbors nearer than tolerance from
113 * a given point.
114 * \details Called in parallel using parallel_for().
115 * \param[in] i index of the query point
116 */
117 32622 void do_it(index_t i) {
118 32622 index_t nb = std::min(index_t(6),nb_points());
119
3/8
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 32622 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32622 times.
 32622 while(!find_nearest_neighbors(i, nb) && nb < nb_points()) {
120 ✗ if(nb == nb_points()) {
121 ✗ break;
122 }
123 ✗ nb += nb / 2;
124 ✗ nb = std::min(nb, nb_points());
125 }
126 32622 }
127
128 private:
129 NearestNeighborSearch* NN_;
130 vector<index_t>& old2new_;
131 double sq_tolerance_;
132 };
133
134 /************************************************************************/
135
136 /**
137 * \brief A comparator for sorting points in lexicographic order.
138 */
139 class ComparePoints {
140 public:
141 /**
142 * \brief Constructs a new ComparePoints
143 * \param[in] points pointer to the array of points
144 * \param[in] dim number of coordinates
145 * \param[in] stride number of doubles between two consecutive points
146 * (= dim if the array is packed).
147 */
148 178 ComparePoints(
149 const double* points, coord_index_t dim, index_t stride
150 178 ) :
151 178 points_(points),
152 178 dim_(dim),
153 178 stride_(stride) {
154 178 }
155
156 /**
157 * \brief Compares two points given their indices.
158 * \param[in] i1 index of the first point
159 * \param[in] i2 index of the second point
160 */
161 6553818 bool is_before(index_t i1, index_t i2) const {
162 6553818 const double* p1 = points_ + i1 * stride_;
163 6553818 const double* p2 = points_ + i2 * stride_;
164
2/2
✓ Branch 0 taken 7559658 times.
✓ Branch 1 taken 83134 times.
 7642792 for(coord_index_t c = 0; c < dim_; c++) {
165
2/2
✓ Branch 0 taken 4179815 times.
✓ Branch 1 taken 3379843 times.
 7559658 if(p1[c] < p2[c]) {
166 4179815 return true;
167 }
168
2/2
✓ Branch 0 taken 2290869 times.
✓ Branch 1 taken 1088974 times.
 3379843 if(p1[c] > p2[c]) {
169 2290869 return false;
170 }
171 }
172 83134 return false;
173 }
174
175 /**
176 * \brief Tests whether two points are identical.
177 * \param[in] i1 index of the first point
178 * \param[in] i2 index of the second point
179 */
180 386211 bool is_same(index_t i1, index_t i2) const {
181 386211 const double* p1 = points_ + i1 * stride_;
182 386211 const double* p2 = points_ + i2 * stride_;
183
2/2
✓ Branch 0 taken 625842 times.
✓ Branch 1 taken 47946 times.
 673788 for(coord_index_t c = 0; c < dim_; c++) {
184
2/2
✓ Branch 0 taken 338265 times.
✓ Branch 1 taken 287577 times.
 625842 if(p1[c] != p2[c]) {
185 338265 return false;
186 }
187 }
188 47946 return true;
189 }
190
191 /**
192 * \brief Compares two points given their indices.
193 * \param[in] i1 index of the first point
194 * \param[in] i2 index of the second point
195 * \return true if point \p i1 is before point \p i2, false
196 * otherwise.
197 */
198 6553818 bool operator() (index_t i1, index_t i2) const {
199 6553818 return is_before(i1, i2);
200 }
201
202 private:
203 const double* points_;
204 coord_index_t dim_;
205 index_t stride_;
206 };
207 }
208
209 /****************************************************************************/
210
211 namespace GEO {
212
213 namespace Geom {
214
215 41 index_t colocate(
216 const double* points,
217 coord_index_t dim,
218 index_t nb_points,
219 vector<index_t>& old2new,
220 double tolerance,
221 index_t stride,
222 const std::string& nn_algo
223 ) {
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
 41 if(nb_points == 0) {
225 ✗ return 0;
226 }
227
228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
 41 if(stride == 0) {
229 ✗ stride = dim;
230 }
231 NearestNeighborSearch_var NN = NearestNeighborSearch::create(
232 dim, nn_algo
233
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 );
234
2/4
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 41 times.
✗ Branch 5 not taken.
 41 NN->set_points(nb_points, points, stride);
235
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 old2new.resize(nb_points, NO_INDEX);
236
1/2
✓ Branch 2 taken 41 times.
✗ Branch 3 not taken.
 41 Colocate colocate_obj(NN, old2new, tolerance);
237
238
3/6
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 41 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 41 times.
✗ Branch 8 not taken.
 82 if(CmdLine::get_arg_bool("sys:multithread")) {
239
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 parallel_for(
240 0, nb_points,
241 32704 [&colocate_obj](index_t i){ colocate_obj.do_it(i); },
242 1, true
243 );
244 } else {
245 ✗ for(index_t i = 0; i < nb_points; i++) {
246 ✗ colocate_obj.do_it(i);
247 }
248 }
249 41 index_t result = 0;
250
2/2
✓ Branch 1 taken 32622 times.
✓ Branch 2 taken 41 times.
 32663 for(index_t i = 0; i < old2new.size(); i++) {
251
5/14
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32622 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32622 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 32622 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
 32622 geo_assert(
252 signed_index_t(old2new[i]) >= 0 &&
253 old2new[i] < nb_points
254 );
255 32622 index_t j = i;
256 // colocate clusters of identical vertices onto smallest index
257
3/4
✓ Branch 1 taken 33977 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1355 times.
✓ Branch 4 taken 32622 times.
 33977 while(old2new[j] != j) {
258
1/2
✓ Branch 1 taken 1355 times.
✗ Branch 2 not taken.
 1355 j = old2new[j];
259 }
260
1/2
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
 32622 old2new[i] = j;
261
3/4
✓ Branch 1 taken 32622 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 31267 times.
✓ Branch 4 taken 1355 times.
 32622 if(old2new[i] == i) {
262 31267 result++;
263 }
264 }
265 41 return result;
266 41 }
267
268 178 index_t colocate_by_lexico_sort(
269 const double* points,
270 coord_index_t dim,
271 index_t nb_points,
272 vector<index_t>& old2new,
273 index_t stride
274 ) {
275
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 178 times.
 178 if(nb_points == 0) {
276 ✗ return 0;
277 }
278
279 178 ComparePoints compare_points(points, dim, stride);
280
1/2
✓ Branch 1 taken 178 times.
✗ Branch 2 not taken.
 178 vector<index_t> sorted_indices(nb_points);
281
2/2
✓ Branch 0 taken 386389 times.
✓ Branch 1 taken 178 times.
 386567 for(index_t i = 0; i < nb_points; i++) {
282
1/2
✓ Branch 1 taken 386389 times.
✗ Branch 2 not taken.
 386389 sorted_indices[i] = i;
283 }
284
1/2
✓ Branch 1 taken 178 times.
✗ Branch 2 not taken.
 178 GEO::sort(
285 178 sorted_indices.begin(), sorted_indices.end(), compare_points
286 );
287
1/2
✓ Branch 1 taken 178 times.
✗ Branch 2 not taken.
 178 old2new.assign(nb_points, NO_INDEX);
288
289 178 index_t nb_distinct = 0;
290
291 178 index_t iv1 = 0;
292
2/2
✓ Branch 0 taken 338443 times.
✓ Branch 1 taken 178 times.
 338621 while(iv1 < nb_points) {
293 338443 nb_distinct++;
294
3/6
✓ Branch 1 taken 338443 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 338443 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 338443 times.
✗ Branch 8 not taken.
 338443 old2new[sorted_indices[iv1]] = sorted_indices[iv1];
295 338443 index_t iv2 = iv1 + 1;
296 338443 while(
297
6/6
✓ Branch 0 taken 386211 times.
✓ Branch 1 taken 178 times.
✓ Branch 2 taken 47946 times.
✓ Branch 3 taken 338265 times.
✓ Branch 4 taken 47946 times.
✓ Branch 5 taken 338443 times.
 772600 iv2 < nb_points &&
298 772422 compare_points.is_same(
299
2/4
✓ Branch 1 taken 386211 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 386211 times.
✗ Branch 5 not taken.
 386211 sorted_indices[iv1], sorted_indices[iv2]
300 )
301 ) {
302
3/6
✓ Branch 1 taken 47946 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 47946 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 47946 times.
✗ Branch 8 not taken.
 47946 old2new[sorted_indices[iv2]] = sorted_indices[iv1];
303 47946 iv2++;
304 }
305 338443 iv1 = iv2;
306 }
307 178 return nb_distinct;
308 178 }
309 }
310 }
311