GCC Code Coverage Report

Directory: ./
File: lib/geogram/points/colocate.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 64 72 88.9%
Functions: 5 5 100.0%
Branches: 67 116 57.8%
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 Colocate(
65 NearestNeighborSearch* NN,
66 vector<index_t>& old2new,
67 double tolerance
68 43 ) :
69 43 NN_(NN),
70
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 43 old2new_(old2new),
71 43 sq_tolerance_(geo_sqr(tolerance)) {
72 }
73
74 /**
75 * \brief Returns the number of points.
76 */
77 index_t nb_points() const {
78 ✗ 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 36614 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 36614 index_t* neighbors = (index_t*) alloca(sizeof(index_t) * nb);
93 36614 double* dist = (double*) alloca(sizeof(double) * nb);
94
95 36614 NN_->get_nearest_neighbors(
96 nb, NN_->point_ptr(i), neighbors, dist
97 );
98
99 36614 index_t smallest = i;
100
1/2
✓ Branch 0 taken 76482 times.
✗ Branch 1 not taken.
 76482 for(index_t jj = 0; jj < nb; jj++) {
101
2/2
✓ Branch 0 taken 36614 times.
✓ Branch 1 taken 39868 times.
 76482 if(dist[jj] > sq_tolerance_) {
102 36614 old2new_[i] = smallest;
103 36614 return true;
104 }
105
2/2
✓ Branch 0 taken 1497 times.
✓ Branch 1 taken 38371 times.
 41365 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 36614 void do_it(index_t i) {
118 36614 index_t nb = std::min(index_t(6),nb_points());
119
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 36614 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 36614 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 36614 }
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 ComparePoints(
149 const double* points, coord_index_t dim, index_t stride
150 179 ) :
151 179 points_(points),
152 179 dim_(dim),
153 179 stride_(stride) {
154 }
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 6575662 bool is_before(index_t i1, index_t i2) const {
162 6575662 const double* p1 = points_ + i1 * stride_;
163 6575662 const double* p2 = points_ + i2 * stride_;
164
2/2
✓ Branch 0 taken 7583238 times.
✓ Branch 1 taken 83134 times.
 7666372 for(coord_index_t c = 0; c < dim_; c++) {
165
2/2
✓ Branch 0 taken 3392277 times.
✓ Branch 1 taken 4190961 times.
 7583238 if(p1[c] < p2[c]) {
166 return true;
167 }
168
2/2
✓ Branch 0 taken 1090710 times.
✓ Branch 1 taken 2301567 times.
 3392277 if(p1[c] > p2[c]) {
169 return false;
170 }
171 }
172 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 bool is_same(index_t i1, index_t i2) const {
181 388206 const double* p1 = points_ + i1 * stride_;
182 388206 const double* p2 = points_ + i2 * stride_;
183
2/2
✓ Branch 0 taken 629092 times.
✓ Branch 1 taken 47946 times.
 677038 for(coord_index_t c = 0; c < dim_; c++) {
184
2/2
✓ Branch 0 taken 288832 times.
✓ Branch 1 taken 340260 times.
 629092 if(p1[c] != p2[c]) {
185 return false;
186 }
187 }
188 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 bool operator() (index_t i1, index_t i2) const {
199
22/36
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 8 times.
✓ Branch 13 taken 37 times.
✓ Branch 14 taken 54 times.
✓ Branch 15 taken 42 times.
✓ Branch 16 taken 988388 times.
✓ Branch 17 taken 387780 times.
✓ Branch 18 taken 1769627 times.
✓ Branch 19 taken 971238 times.
✓ Branch 20 taken 1383676 times.
✓ Branch 21 taken 971238 times.
✓ Branch 22 taken 19563 times.
✓ Branch 23 taken 20237 times.
✓ Branch 24 taken 9098 times.
✓ Branch 25 taken 10465 times.
✓ Branch 26 taken 5478 times.
✓ Branch 27 taken 4987 times.
✓ Branch 28 taken 9228 times.
✓ Branch 29 taken 11009 times.
✓ Branch 30 taken 5415 times.
✓ Branch 31 taken 5594 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 426 times.
✓ Branch 35 taken 2074 times.
 6514388 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 43 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 taken 43 times.
✗ Branch 1 not taken.
 43 if(nb_points == 0) {
225 return 0;
226 }
227
228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43 times.
 43 if(stride == 0) {
229 ✗ stride = dim;
230 }
231 NearestNeighborSearch_var NN = NearestNeighborSearch::create(
232 dim, nn_algo
233 43 );
234
2/4
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
 43 NN->set_points(nb_points, points, stride);
235
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 43 old2new.resize(nb_points, NO_INDEX);
236 Colocate colocate_obj(NN, old2new, tolerance);
237
238
3/6
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 86 if(CmdLine::get_arg_bool("sys:multithread")) {
239
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 43 parallel_for(
240 0, nb_points,
241 36657 [&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 index_t result = 0;
250
2/2
✓ Branch 0 taken 36614 times.
✓ Branch 1 taken 43 times.
 73271 for(index_t i = 0; i < old2new.size(); i++) {
251
2/12
✓ Branch 0 taken 36614 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 36614 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 36614 geo_assert(
252 signed_index_t(old2new[i]) >= 0 &&
253 old2new[i] < nb_points
254 );
255 index_t j = i;
256 // colocate clusters of identical vertices onto smallest index
257
2/2
✓ Branch 0 taken 1355 times.
✓ Branch 1 taken 36614 times.
 37969 while(old2new[j] != j) {
258 j = old2new[j];
259 }
260 36614 old2new[i] = j;
261
2/2
✓ Branch 0 taken 35259 times.
✓ Branch 1 taken 1355 times.
 36614 if(old2new[i] == i) {
262 35259 result++;
263 }
264 }
265 return result;
266 }
267
268 179 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 taken 179 times.
✗ Branch 1 not taken.
 179 if(nb_points == 0) {
276 return 0;
277 }
278
279 ComparePoints compare_points(points, dim, stride);
280 vector<index_t> sorted_indices(nb_points);
281
2/2
✓ Branch 0 taken 388385 times.
✓ Branch 1 taken 179 times.
 388564 for(index_t i = 0; i < nb_points; i++) {
282 388385 sorted_indices[i] = i;
283 }
284 179 GEO::sort(
285
1/4
✓ Branch 1 taken 179 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 179 sorted_indices.begin(), sorted_indices.end(), compare_points
286 );
287
1/2
✓ Branch 1 taken 179 times.
✗ Branch 2 not taken.
 179 old2new.assign(nb_points, NO_INDEX);
288
289 index_t nb_distinct = 0;
290
291 index_t iv1 = 0;
292
2/2
✓ Branch 0 taken 340439 times.
✓ Branch 1 taken 179 times.
 340618 while(iv1 < nb_points) {
293 340439 nb_distinct++;
294 340439 old2new[sorted_indices[iv1]] = sorted_indices[iv1];
295 340439 index_t iv2 = iv1 + 1;
296 340439 while(
297
4/4
✓ Branch 0 taken 388206 times.
✓ Branch 1 taken 179 times.
✓ Branch 2 taken 47946 times.
✓ Branch 3 taken 340260 times.
 776591 iv2 < nb_points &&
298 388206 compare_points.is_same(
299 sorted_indices[iv1], sorted_indices[iv2]
300 )
301 ) {
302 47946 old2new[sorted_indices[iv2]] = sorted_indices[iv1];
303 47946 iv2++;
304 }
305 iv1 = iv2;
306 }
307 return nb_distinct;
308 }
309 }
310 }
311