| 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_topology.h> | ||
| 41 | #include <geogram/mesh/mesh.h> | ||
| 42 | #include <geogram/basic/memory.h> | ||
| 43 | #include <geogram/basic/logger.h> | ||
| 44 | #include <stack> | ||
| 45 | |||
| 46 | namespace { | ||
| 47 | |||
| 48 | using namespace GEO; | ||
| 49 | |||
| 50 | /** | ||
| 51 | * \brief Computes the number of surface vertices that are not isolated. | ||
| 52 | * \param[in] M the mesh | ||
| 53 | * \return the number of vertices of a mesh with at least | ||
| 54 | * one incident surface facet | ||
| 55 | */ | ||
| 56 | 151 | index_t nb_non_isolated_surface_vertices(const Mesh& M) { | |
| 57 | index_t result = 0; | ||
| 58 | 151 | std::vector<bool> visited(M.vertices.nb(), false); | |
| 59 |
2/2✓ Branch 0 taken 1410285 times.
✓ Branch 1 taken 151 times.
|
1410436 | for(index_t c: M.facet_corners) { |
| 60 | visited[M.facet_corners.vertex(c)] = true; | ||
| 61 | } | ||
| 62 |
2/2✓ Branch 0 taken 315933 times.
✓ Branch 1 taken 151 times.
|
316084 | for(index_t v: M.vertices) { |
| 63 |
1/2✓ Branch 0 taken 315933 times.
✗ Branch 1 not taken.
|
315933 | if(visited[v]) { |
| 64 | 315933 | ++result; | |
| 65 | } | ||
| 66 | } | ||
| 67 | 151 | return result; | |
| 68 | } | ||
| 69 | } | ||
| 70 | |||
| 71 | namespace GEO { | ||
| 72 | |||
| 73 | 390 | index_t get_connected_components( | |
| 74 | const Mesh& M, vector<index_t>& component | ||
| 75 | ) { | ||
| 76 | index_t nb_components = 0; | ||
| 77 | 390 | component.assign(M.facets.nb(), NO_INDEX); | |
| 78 |
4/4✓ Branch 0 taken 1340 times.
✓ Branch 1 taken 1592789 times.
✓ Branch 2 taken 1594129 times.
✓ Branch 3 taken 390 times.
|
1594519 | for(index_t f: M.facets) { |
| 79 |
2/2✓ Branch 0 taken 1340 times.
✓ Branch 1 taken 1592789 times.
|
1594129 | if(component[f] == NO_INDEX) { |
| 80 | std::stack<index_t> S; | ||
| 81 | S.push(f); | ||
| 82 | 1340 | component[f] = nb_components; | |
| 83 | do { | ||
| 84 |
2/2✓ Branch 0 taken 1569869 times.
✓ Branch 1 taken 24260 times.
|
1594129 | index_t cur_f = S.top(); |
| 85 | S.pop(); | ||
| 86 |
4/4✓ Branch 0 taken 1376943 times.
✓ Branch 1 taken 217186 times.
✓ Branch 2 taken 5012070 times.
✓ Branch 3 taken 1594129 times.
|
8200328 | for(index_t adj_f: M.facets.adjacent(cur_f)) { |
| 87 |
4/4✓ Branch 0 taken 4986732 times.
✓ Branch 1 taken 25338 times.
✓ Branch 2 taken 1592789 times.
✓ Branch 3 taken 3393943 times.
|
5012070 | if(adj_f != NO_FACET && component[adj_f] == NO_INDEX) { |
| 88 |
1/2✓ Branch 1 taken 1592789 times.
✗ Branch 2 not taken.
|
3185578 | S.push(index_t(adj_f)); |
| 89 | 1592789 | component[adj_f] = nb_components; | |
| 90 | } | ||
| 91 | } | ||
| 92 |
2/2✓ Branch 0 taken 1592789 times.
✓ Branch 1 taken 1340 times.
|
1594129 | } while(!S.empty()); |
| 93 | 1340 | nb_components++; | |
| 94 | } | ||
| 95 | } | ||
| 96 | 390 | return nb_components; | |
| 97 | } | ||
| 98 | |||
| 99 | 128 | index_t GEOGRAM_API get_connected_components( | |
| 100 | const Mesh& M, Attribute<index_t>& component | ||
| 101 | ) { | ||
| 102 | index_t nb_components = 0; | ||
| 103 |
2/2✓ Branch 0 taken 836224 times.
✓ Branch 1 taken 128 times.
|
836352 | for(index_t f: M.facets) { |
| 104 | 836224 | component[f] = NO_INDEX; | |
| 105 | } | ||
| 106 |
4/4✓ Branch 0 taken 9061 times.
✓ Branch 1 taken 827163 times.
✓ Branch 2 taken 836224 times.
✓ Branch 3 taken 128 times.
|
836352 | for(index_t f: M.facets) { |
| 107 |
2/2✓ Branch 0 taken 9061 times.
✓ Branch 1 taken 827163 times.
|
836224 | if(component[f] == NO_INDEX) { |
| 108 | std::stack<index_t> S; | ||
| 109 | S.push(f); | ||
| 110 | 9061 | component[f] = nb_components; | |
| 111 | do { | ||
| 112 |
2/2✓ Branch 0 taken 826628 times.
✓ Branch 1 taken 9596 times.
|
836224 | index_t cur_f = S.top(); |
| 113 | S.pop(); | ||
| 114 |
3/4✓ Branch 0 taken 836224 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2508672 times.
✓ Branch 3 taken 836224 times.
|
4181120 | for(index_t adj_f: M.facets.adjacent(cur_f)) { |
| 115 |
4/4✓ Branch 0 taken 2344714 times.
✓ Branch 1 taken 163958 times.
✓ Branch 2 taken 827163 times.
✓ Branch 3 taken 1517551 times.
|
2508672 | if(adj_f != NO_FACET && component[adj_f] == NO_INDEX) { |
| 116 |
1/2✓ Branch 1 taken 827163 times.
✗ Branch 2 not taken.
|
1654326 | S.push(index_t(adj_f)); |
| 117 | 827163 | component[adj_f] = nb_components; | |
| 118 | } | ||
| 119 | } | ||
| 120 |
2/2✓ Branch 0 taken 827163 times.
✓ Branch 1 taken 9061 times.
|
836224 | } while(!S.empty()); |
| 121 | 9061 | nb_components++; | |
| 122 | } | ||
| 123 | } | ||
| 124 | 128 | return nb_components; | |
| 125 | } | ||
| 126 | |||
| 127 |
1/2✓ Branch 1 taken 153 times.
✗ Branch 2 not taken.
|
153 | index_t mesh_nb_connected_components(const Mesh& M) { |
| 128 | vector<index_t> component; | ||
| 129 |
1/2✓ Branch 1 taken 153 times.
✗ Branch 2 not taken.
|
306 | return get_connected_components(M, component); |
| 130 | } | ||
| 131 | |||
| 132 | 151 | signed_index_t mesh_Xi(const Mesh& M) { | |
| 133 | 151 | index_t nb_v = nb_non_isolated_surface_vertices(M); | |
| 134 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 151 times.
|
151 | if(nb_v != M.vertices.nb()) { |
| 135 | ✗ | if(M.cells.nb() == 0) { | |
| 136 | ✗ | Logger::warn("Topology") | |
| 137 | ✗ | << "Surface mesh has " << M.vertices.nb() - nb_v | |
| 138 | << " isolated vertices" | ||
| 139 | << std::endl; | ||
| 140 | } else { | ||
| 141 | ✗ | Logger::out("Topology") | |
| 142 | ✗ | << "Surface mesh has " << M.vertices.nb() - nb_v | |
| 143 | << " isolated vertices " | ||
| 144 | << " (but they may be attached to tetrahedra)" | ||
| 145 | << std::endl; | ||
| 146 | } | ||
| 147 | } | ||
| 148 | 151 | signed_index_t result = signed_index_t(nb_v + M.facets.nb()); | |
| 149 |
2/2✓ Branch 0 taken 393534 times.
✓ Branch 1 taken 151 times.
|
393685 | for(index_t f: M.facets) { |
| 150 |
2/2✓ Branch 0 taken 1410285 times.
✓ Branch 1 taken 393534 times.
|
1803819 | for(index_t c: M.facets.corners(f)) { |
| 151 | index_t f2 = M.facet_corners.adjacent_facet(c); | ||
| 152 |
2/2✓ Branch 0 taken 709418 times.
✓ Branch 1 taken 700867 times.
|
1410285 | if(f2 == NO_FACET || f > f2) { |
| 153 | 709418 | --result; | |
| 154 | } // We count each edge once, | ||
| 155 | } | ||
| 156 | } | ||
| 157 | 151 | return result; | |
| 158 | } | ||
| 159 | |||
| 160 | 153 | signed_index_t mesh_nb_borders(const Mesh& M) { | |
| 161 | // Step 1: chain vertices around borders | ||
| 162 | 153 | std::vector<index_t> next_around_border(M.vertices.nb(),NO_VERTEX); | |
| 163 |
2/2✓ Branch 0 taken 517045 times.
✓ Branch 1 taken 152 times.
|
517197 | for(index_t f: M.facets) { |
| 164 |
2/2✓ Branch 0 taken 1781162 times.
✓ Branch 1 taken 517044 times.
|
2298206 | for(index_t c1: M.facets.corners(f)) { |
| 165 |
2/2✓ Branch 0 taken 7723 times.
✓ Branch 1 taken 1773439 times.
|
1781162 | if(M.facet_corners.adjacent_facet(c1) == NO_FACET) { |
| 166 | index_t c2 = M.facets.next_corner_around_facet(f, c1); | ||
| 167 | index_t v1 = M.facet_corners.vertex(c1); | ||
| 168 | index_t v2 = M.facet_corners.vertex(c2); | ||
| 169 |
2/2✓ Branch 0 taken 7722 times.
✓ Branch 1 taken 1 times.
|
7723 | if(next_around_border[v1] != NO_VERTEX) { |
| 170 | // If this happens, then the same vertex | ||
| 171 | // is incident to more than two edges on | ||
| 172 | // the border (non-manifold configuration, | ||
| 173 | // return "error value" -1) | ||
| 174 | return -1; | ||
| 175 | } | ||
| 176 | // May happen with non-manifold configuration, | ||
| 177 | // where several connected component of the | ||
| 178 | // border can touch the same vertex several | ||
| 179 | // times. | ||
| 180 | 7722 | next_around_border[v1] = v2; | |
| 181 | } | ||
| 182 | } | ||
| 183 | } | ||
| 184 | // Step 2: count connected components of the borders | ||
| 185 | index_t result = 0; | ||
| 186 |
2/2✓ Branch 0 taken 377456 times.
✓ Branch 1 taken 152 times.
|
377608 | for(index_t v: M.vertices) { |
| 187 |
2/2✓ Branch 0 taken 159 times.
✓ Branch 1 taken 377297 times.
|
377456 | if(next_around_border[v] != NO_VERTEX) { |
| 188 | 159 | result++; | |
| 189 | index_t cur = v; | ||
| 190 |
2/2✓ Branch 0 taken 7716 times.
✓ Branch 1 taken 159 times.
|
7875 | while(next_around_border[cur] != NO_VERTEX) { |
| 191 | index_t next = next_around_border[cur]; | ||
| 192 | 7716 | next_around_border[cur] = NO_VERTEX; | |
| 193 | cur = next; | ||
| 194 | } | ||
| 195 | } | ||
| 196 | } | ||
| 197 | 152 | return signed_index_t(result); | |
| 198 | } | ||
| 199 | |||
| 200 | 69 | bool meshes_have_same_topology( | |
| 201 | const Mesh& M1, const Mesh& M2, bool verbose | ||
| 202 | ) { | ||
| 203 | 69 | signed_index_t Xi1 = mesh_Xi(M1); | |
| 204 | 69 | signed_index_t Xi2 = mesh_Xi(M2); | |
| 205 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
|
69 | if(!verbose && Xi1 != Xi2) { |
| 206 | return false; | ||
| 207 | } | ||
| 208 | |||
| 209 | 69 | signed_index_t nbB1 = mesh_nb_borders(M1); | |
| 210 | 69 | signed_index_t nbB2 = mesh_nb_borders(M2); | |
| 211 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
|
69 | if(!verbose && nbB1 != nbB2) { |
| 212 | return false; | ||
| 213 | } | ||
| 214 | |||
| 215 | 69 | index_t nb_conn1 = mesh_nb_connected_components(M1); | |
| 216 | 69 | index_t nb_conn2 = mesh_nb_connected_components(M2); | |
| 217 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
|
69 | if(!verbose && nb_conn1 != nb_conn2) { |
| 218 | return false; | ||
| 219 | } | ||
| 220 |
1/2✓ Branch 0 taken 69 times.
✗ Branch 1 not taken.
|
69 | if(!verbose) { |
| 221 | return true; | ||
| 222 | } | ||
| 223 | |||
| 224 |
3/4✓ Branch 0 taken 67 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 67 times.
|
69 | bool result = (Xi1 == Xi2 && nbB1 == nbB2 && nb_conn1 == nb_conn2); |
| 225 |
1/2✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
|
69 | Logger::out("Topology") |
| 226 | << "M1: Xi=" << Xi1 << " nbB=" << nbB1 | ||
| 227 | << " nbConn=" << nb_conn1 << std::endl; | ||
| 228 | |||
| 229 |
1/2✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
|
69 | Logger::out("Topology") |
| 230 | << "M2: Xi=" << Xi2 << " nbB=" << nbB2 | ||
| 231 | << " nbConn=" << nb_conn2 << std::endl; | ||
| 232 | |||
| 233 |
4/6✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 67 times.
✓ Branch 7 taken 69 times.
✗ Branch 8 not taken.
|
71 | Logger::out("Topology") << (result ? "match." : "mismatch.") |
| 234 | << std::endl; | ||
| 235 | 69 | return result; | |
| 236 | } | ||
| 237 | |||
| 238 |
1/2✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
|
19 | void reorient_connected_components(Mesh& surf) { |
| 239 | vector<index_t> component; | ||
| 240 |
1/2✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
|
19 | index_t nb_components = get_connected_components(surf, component); |
| 241 | |||
| 242 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
19 | vector<vec3> comp_G(nb_components,{0.0,0.0,0.0}); |
| 243 |
1/4✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
19 | vector<index_t> comp_N(nb_components,0); |
| 244 |
1/4✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
19 | vector<double> comp_signed_V(nb_components,0.0); |
| 245 | |||
| 246 |
2/2✓ Branch 0 taken 48562 times.
✓ Branch 1 taken 19 times.
|
48581 | for(index_t f: surf.facets) { |
| 247 | 48562 | index_t comp = component[f]; | |
| 248 |
2/2✓ Branch 0 taken 145686 times.
✓ Branch 1 taken 48562 times.
|
242810 | for(index_t lv=0; lv<surf.facets.nb_vertices(f); ++lv) { |
| 249 | index_t v = surf.facets.vertex(f,lv); | ||
| 250 | comp_G[comp] += surf.vertices.point(v); | ||
| 251 | 145686 | ++comp_N[comp]; | |
| 252 | } | ||
| 253 | } | ||
| 254 | |||
| 255 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 19 times.
|
55 | for(index_t comp=0; comp<nb_components; ++comp) { |
| 256 | 36 | comp_G[comp] /= double(comp_N[comp]); | |
| 257 | } | ||
| 258 | |||
| 259 |
2/2✓ Branch 0 taken 48562 times.
✓ Branch 1 taken 19 times.
|
48581 | for(index_t f: surf.facets) { |
| 260 | 48562 | index_t comp = component[f]; | |
| 261 |
2/2✓ Branch 0 taken 48562 times.
✓ Branch 1 taken 48562 times.
|
97124 | for(auto [ p1, p2, p3] : surf.facets.triangle_points(f)) { |
| 262 | 48562 | comp_signed_V[comp] += Geom::tetra_signed_volume( | |
| 263 | comp_G[comp],p1,p2,p3 | ||
| 264 | ); | ||
| 265 | } | ||
| 266 | } | ||
| 267 | |||
| 268 |
2/2✓ Branch 0 taken 48562 times.
✓ Branch 1 taken 19 times.
|
48581 | for(index_t f: surf.facets) { |
| 269 | 48562 | index_t comp = component[f]; | |
| 270 |
2/2✓ Branch 0 taken 9776 times.
✓ Branch 1 taken 38786 times.
|
48562 | if(comp_signed_V[comp] < 0.0) { |
| 271 |
1/2✓ Branch 1 taken 9776 times.
✗ Branch 2 not taken.
|
9776 | surf.facets.flip(f); |
| 272 | } | ||
| 273 | } | ||
| 274 | 19 | } | |
| 275 | |||
| 276 | } | ||
| 277 |