| 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_surface_intersection.h> | ||
| 41 | #include <geogram/mesh/mesh_surface_intersection_internal.h> | ||
| 42 | #include <geogram/mesh/mesh_AABB.h> | ||
| 43 | #include <geogram/mesh/mesh_repair.h> | ||
| 44 | #include <geogram/mesh/mesh_fill_holes.h> | ||
| 45 | #include <geogram/mesh/mesh_geometry.h> | ||
| 46 | #include <geogram/mesh/mesh_topology.h> | ||
| 47 | #include <geogram/mesh/mesh_io.h> | ||
| 48 | #include <geogram/mesh/boxes_intersections.h> | ||
| 49 | #include <geogram/mesh/index.h> | ||
| 50 | #include <geogram/delaunay/CDT_2d.h> | ||
| 51 | #include <geogram/numerics/predicates.h> | ||
| 52 | #include <geogram/numerics/expansion_nt.h> | ||
| 53 | #include <geogram/basic/stopwatch.h> | ||
| 54 | #include <geogram/basic/permutation.h> | ||
| 55 | #include <geogram/basic/boolean_expression.h> | ||
| 56 | #include <geogram/basic/debug_stream.h> | ||
| 57 | #include <geogram/basic/algorithm.h> | ||
| 58 | |||
| 59 | #include <sstream> | ||
| 60 | #include <stack> | ||
| 61 | |||
| 62 | #ifdef GEO_COMPILER_CLANG | ||
| 63 | // I'm using long long | ||
| 64 | #pragma GCC diagnostic ignored "-Wc++98-compat-pedantic" | ||
| 65 | #endif | ||
| 66 | |||
| 67 | // If defined, displays status messages and saves files whenever some | ||
| 68 | // error conditions are met. | ||
| 69 | // #define MESH_SURFACE_INTERSECTION_DEBUG | ||
| 70 | |||
| 71 | namespace { | ||
| 72 | using namespace GEO; | ||
| 73 | /** | ||
| 74 | * \brief Removes all facets that have their three vertices aligned | ||
| 75 | */ | ||
| 76 | 60 | void remove_degenerate_triangles(Mesh& M) { | |
| 77 | vector<index_t> remove_f(M.facets.nb()); | ||
| 78 |
1/4✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
60 | parallel_for(0, M.facets.nb(), [&M, &remove_f](index_t f) { |
| 79 |
1/2✓ Branch 0 taken 114536 times.
✗ Branch 1 not taken.
|
114536 | const vec3& p1 = M.facets.point(f,0); |
| 80 | const vec3& p2 = M.facets.point(f,1); | ||
| 81 | const vec3& p3 = M.facets.point(f,2); | ||
| 82 | 114536 | remove_f[f] = PCK::aligned_3d(p1,p2,p3); | |
| 83 | 114536 | }); | |
| 84 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | M.facets.delete_elements(remove_f); |
| 85 | 60 | } | |
| 86 | |||
| 87 | /** | ||
| 88 | * \brief Enumerates the connected components in a facet attribute | ||
| 89 | * \param[in] M a reference to the mesh | ||
| 90 | * \param[in] attribute the name of the facet attribute | ||
| 91 | * \return the number of found connected components | ||
| 92 | */ | ||
| 93 | 128 | index_t get_surface_connected_components( | |
| 94 | Mesh& M, const std::string& attribute = "chart" | ||
| 95 | ) { | ||
| 96 | Attribute<index_t> chart(M.facets.attributes(), attribute); | ||
| 97 |
1/2✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
|
256 | return GEO::get_connected_components(M, chart); |
| 98 | } | ||
| 99 | |||
| 100 | /** | ||
| 101 | * \brief Computes the intersection between two mesh triangular facets | ||
| 102 | * \details This function is just a wrapper around triangles_intersections() | ||
| 103 | * for Mesh facets. | ||
| 104 | * \param[in] M the mesh | ||
| 105 | * \param[in] f1 , f2 the two facets | ||
| 106 | * \param[out] I a vector of triangle intersections | ||
| 107 | * \retval true if there was an intersection | ||
| 108 | * \retval false otherwise | ||
| 109 | */ | ||
| 110 |
1/2✓ Branch 0 taken 981351 times.
✗ Branch 1 not taken.
|
981351 | bool mesh_facets_intersect( |
| 111 | Mesh& M, index_t f1, index_t f2, TriangleIsects& I | ||
| 112 | ) { | ||
| 113 | geo_debug_assert(M.facets.nb_vertices(f1) == 3); | ||
| 114 | geo_debug_assert(M.facets.nb_vertices(f2) == 3); | ||
| 115 | index_t v1 = M.facets.vertex(f1,0); | ||
| 116 | index_t v2 = M.facets.vertex(f1,1); | ||
| 117 | index_t v3 = M.facets.vertex(f1,2); | ||
| 118 | index_t w1 = M.facets.vertex(f2,0); | ||
| 119 | index_t w2 = M.facets.vertex(f2,1); | ||
| 120 | index_t w3 = M.facets.vertex(f2,2); | ||
| 121 | 981351 | return triangles_intersections( | |
| 122 | M.vertices.point(v1), | ||
| 123 | M.vertices.point(v2), | ||
| 124 | M.vertices.point(v3), | ||
| 125 | M.vertices.point(w1), | ||
| 126 | M.vertices.point(w2), | ||
| 127 | M.vertices.point(w3), | ||
| 128 | v1,v2,v3, | ||
| 129 | w1,w2,w3, | ||
| 130 | I | ||
| 131 | 981351 | ); | |
| 132 | } | ||
| 133 | |||
| 134 | /** | ||
| 135 | * \brief Tests whether a segment intersects a triangle | ||
| 136 | * \param[in] q1 , q2 the two extremities of the segment | ||
| 137 | * \param[in] p1 , p2 , p3 the three verties of the triangle | ||
| 138 | * \retval true if the segment has an intersection with the | ||
| 139 | * interior of the triangle | ||
| 140 | * \retval false otherwise | ||
| 141 | * \details Degenerate configurations (segment passing through vertex, | ||
| 142 | * edge, or co-planar with triangle) are symbolically perturbed. | ||
| 143 | */ | ||
| 144 | 2104 | template <class POINT> bool segment_triangle_intersection_SOS( | |
| 145 | const POINT& q1, const POINT& q2, | ||
| 146 | const POINT& p1, const POINT& p2, const POINT& p3 | ||
| 147 | ) { | ||
| 148 | ✗ | Sign o1 = PCK::orient_3d_SOS(q1,p1,p2,p3); | |
| 149 | ✗ | Sign o2 = PCK::orient_3d_SOS(q2,p1,p2,p3); | |
| 150 | |||
| 151 | // There is no intersection if q1 and q2 are on the | ||
| 152 | // same side of the supporting plane of (p1,p2,p3) | ||
| 153 |
2/2✓ Branch 0 taken 752 times.
✓ Branch 1 taken 300 times.
|
2104 | if(o1 == o2) { |
| 154 | return false; | ||
| 155 | } | ||
| 156 | |||
| 157 | // There is an intersection if the three tetrahedra | ||
| 158 | // formed by [q1,q2] and the three edges of the triangle | ||
| 159 | // have the same orientation | ||
| 160 | ✗ | Sign s1 = PCK::orient_3d_SOS(q1,q2,p1,p2); | |
| 161 | ✗ | Sign s2 = PCK::orient_3d_SOS(q1,q2,p2,p3); | |
| 162 |
2/2✓ Branch 0 taken 207 times.
✓ Branch 1 taken 93 times.
|
600 | if(s1*s2 < 0) { |
| 163 | return false; | ||
| 164 | } | ||
| 165 | ✗ | Sign s3 = PCK::orient_3d_SOS(q1,q2,p3,p1); | |
| 166 |
3/4✓ Branch 0 taken 91 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
186 | return(s2*s3 > 0 && s3*s1 > 0); |
| 167 | } | ||
| 168 | } | ||
| 169 | |||
| 170 | |||
| 171 | namespace GEO { | ||
| 172 | |||
| 173 | 60 | MeshSurfaceIntersection::MeshSurfaceIntersection(Mesh& M) : | |
| 174 | lock_(GEOGRAM_SPINLOCK_INIT), | ||
| 175 | 60 | mesh_(M), | |
| 176 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
120 | vertex_to_exact_point_(M.vertices.attributes(), "exact_point"), |
| 177 | 60 | dry_run_(false), | |
| 178 | halfedges_(*this), | ||
| 179 | radial_bundles_(*this), | ||
| 180 | 60 | radial_polylines_(*this) | |
| 181 | { | ||
| 182 |
2/2✓ Branch 0 taken 57904 times.
✓ Branch 1 taken 60 times.
|
57964 | for(index_t v: mesh_.vertices) { |
| 183 | 57904 | vertex_to_exact_point_[v] = nullptr; | |
| 184 | } | ||
| 185 | 60 | verbose_ = false; | |
| 186 | 60 | fine_verbose_ = false; | |
| 187 | 60 | delaunay_ = true; | |
| 188 | 60 | detect_intersecting_neighbors_ = true; | |
| 189 | 60 | use_radial_sort_ = true; | |
| 190 | 60 | monster_threshold_ = NO_INDEX; | |
| 191 | 60 | has_operand_bits_ = false; | |
| 192 | 60 | skeleton_ = nullptr; | |
| 193 | 60 | skeleton_trim_fins_ = false; | |
| 194 | 60 | interpolate_attributes_ = false; | |
| 195 | 60 | } | |
| 196 | |||
| 197 | 60 | MeshSurfaceIntersection::~MeshSurfaceIntersection() { | |
| 198 | 60 | vertex_to_exact_point_.destroy(); | |
| 199 | 120 | } | |
| 200 | |||
| 201 | ✗ | void MeshSurfaceIntersection::remove_external_shell() { | |
| 202 | vector<index_t> remove_f; | ||
| 203 | ✗ | mark_external_shell(remove_f); | |
| 204 | ✗ | mesh_.facets.delete_elements(remove_f); | |
| 205 | ✗ | mesh_.facets.connect(); | |
| 206 | ✗ | } | |
| 207 | |||
| 208 | 27 | void MeshSurfaceIntersection::remove_internal_shells() { | |
| 209 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
|
27 | if(mesh_.facets.nb() == 0) { |
| 210 | ✗ | return; | |
| 211 | } | ||
| 212 | vector<index_t> remove_f; | ||
| 213 |
1/2✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
|
27 | mark_external_shell(remove_f); |
| 214 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 27 times.
|
217093 | for(index_t& i: remove_f) { |
| 215 | 217066 | i=1-i; | |
| 216 | } | ||
| 217 |
1/2✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
|
27 | mesh_.facets.delete_elements(remove_f); |
| 218 | } | ||
| 219 | |||
| 220 | ✗ | void MeshSurfaceIntersection::remove_fins() { | |
| 221 | // TODO | ||
| 222 | ✗ | } | |
| 223 | |||
| 224 | 60 | void MeshSurfaceIntersection::intersect_prologue() { | |
| 225 |
1/2✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
|
60 | Stopwatch W("Prologue", verbose_); |
| 226 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 55 times.
|
60 | if(!mesh_.facets.are_simplices()) { |
| 227 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | tessellate_facets(mesh_,3); |
| 228 | } | ||
| 229 | |||
| 230 | Attribute<index_t> operand_bit; | ||
| 231 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | operand_bit.bind_if_is_defined( |
| 232 |
4/6✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 46 times.
✓ Branch 7 taken 14 times.
|
120 | mesh_.facets.attributes(), "operand_bit" |
| 233 | ); | ||
| 234 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 14 times.
|
60 | has_operand_bits_ = operand_bit.is_bound(); |
| 235 | if(!operand_bit.is_bound()) { | ||
| 236 | // TODO: not good, there might be more than 32 components | ||
| 237 | // (and why do we need to do that BTW?) | ||
| 238 |
2/4✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
|
14 | get_surface_connected_components(mesh_,"operand_bit"); |
| 239 |
2/6✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
14 | operand_bit.bind(mesh_.facets.attributes(), "operand_bit"); |
| 240 |
2/2✓ Branch 0 taken 50108 times.
✓ Branch 1 taken 14 times.
|
50122 | for(index_t f: mesh_.facets) { |
| 241 | 50108 | operand_bit[f] = (index_t(1) << operand_bit[f]) ; | |
| 242 | } | ||
| 243 | } | ||
| 244 | |||
| 245 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | remove_degenerate_triangles(mesh_); |
| 246 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | mesh_colocate_vertices_no_check(mesh_); |
| 247 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | mesh_remove_bad_facets_no_check(mesh_); |
| 248 | |||
| 249 | // Set symbolic perturbation mode to lexicographic order | ||
| 250 | // on point coordinates instead of point indices only, | ||
| 251 | // Needed to get compatible triangulations on coplanar faces | ||
| 252 | // (example, cubes that touch on a facet). | ||
| 253 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | SOS_bkp_ = PCK::get_SOS_mode(); |
| 254 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | PCK::set_SOS_mode(PCK::SOS_LEXICO); |
| 255 | 60 | } | |
| 256 | |||
| 257 | 60 | void MeshSurfaceIntersection::intersect_get_intersections( | |
| 258 | vector<IsectInfo>& intersections | ||
| 259 | ) { | ||
| 260 |
2/4✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 60 times.
✗ Branch 6 not taken.
|
120 | Stopwatch Wtot("Find isects", verbose_); |
| 261 | { | ||
| 262 | vector<std::pair<index_t, index_t> > FF; | ||
| 263 | 981351 | auto report_BB = [this,&FF](index_t f1, index_t f2) { | |
| 264 | // Optionally skip facet pairs that share a vertex or an edge | ||
| 265 | if( | ||
| 266 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 981351 times.
|
981351 | !detect_intersecting_neighbors_ && ( |
| 267 | ✗ | (mesh_.facets.find_adjacent(f1,f2)!=NO_INDEX) || | |
| 268 | ✗ | (mesh_.facets.find_common_vertex(f1,f2)!=NO_INDEX) | |
| 269 | ) | ||
| 270 | ) { | ||
| 271 | return; | ||
| 272 | } | ||
| 273 | 981351 | FF.push_back(std::make_pair(f1,f2)); | |
| 274 | }; | ||
| 275 | |||
| 276 | // For now, keeping parallel AABB construction and evaluation, | ||
| 277 | // ZE seems to be slower, see stats in: | ||
| 278 | // https://github.com/BrunoLevy/geogram/issues/362 | ||
| 279 | static constexpr bool use_ZE = false; // Zoromodian-Edelsbrunner | ||
| 280 | static constexpr bool ZE_groups = false; | ||
| 281 | static constexpr bool ZE_groups_auto_isect = true; | ||
| 282 | |||
| 283 | // Implementation using AABB | ||
| 284 | if(!use_ZE) { | ||
| 285 |
3/6✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
|
60 | Stopwatch* W = new Stopwatch("AABB build", verbose_); |
| 286 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | MeshFacetsAABB AABB(mesh_, AABB_INDIRECT); |
| 287 | 60 | delete W; | |
| 288 | // Get candidate pairs of intersecting facets | ||
| 289 |
4/8✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 60 times.
✗ Branch 11 not taken.
|
120 | W = new Stopwatch("AABB box-box", verbose_); |
| 290 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | AABB.compute_facet_bbox_intersections(report_BB); |
| 291 | 60 | delete W; | |
| 292 | } | ||
| 293 | |||
| 294 | // Alternative implementation using Zoromodian-Edelsbrunner | ||
| 295 | if(use_ZE) { | ||
| 296 | Stopwatch W("AABB box-box", verbose_); | ||
| 297 | vector<Box3d> boxes(mesh_.facets.nb()); | ||
| 298 | parallel_for( | ||
| 299 | 0, mesh_.facets.nb(), [&](index_t f) { | ||
| 300 | vec3 p0 = mesh_.facets.point(f,0); | ||
| 301 | vec3 p1 = mesh_.facets.point(f,1); | ||
| 302 | vec3 p2 = mesh_.facets.point(f,2); | ||
| 303 | boxes[f].xyz_min[0] = std::min(std::min(p0.x,p1.x),p2.x); | ||
| 304 | boxes[f].xyz_min[1] = std::min(std::min(p0.y,p1.y),p2.y); | ||
| 305 | boxes[f].xyz_min[2] = std::min(std::min(p0.z,p1.z),p2.z); | ||
| 306 | boxes[f].xyz_max[0] = std::max(std::max(p0.x,p1.x),p2.x); | ||
| 307 | boxes[f].xyz_max[1] = std::max(std::max(p0.y,p1.y),p2.y); | ||
| 308 | boxes[f].xyz_max[2] = std::max(std::max(p0.z,p1.z),p2.z); | ||
| 309 | } | ||
| 310 | ); | ||
| 311 | |||
| 312 | |||
| 313 | index_t nb_groups = 0; | ||
| 314 | vector<index_t> indices; | ||
| 315 | vector<index_t> group_ptr; | ||
| 316 | bool has_facets_in_several_groups = false; | ||
| 317 | |||
| 318 | auto has_one_bit = [](index_t x)->bool { | ||
| 319 | return ((x != 0) && ((x & (x - 1))) == 0); | ||
| 320 | }; | ||
| 321 | |||
| 322 | // Compute groups | ||
| 323 | if(ZE_groups && has_operand_bits_) { | ||
| 324 | Attribute<index_t> operand_bit; | ||
| 325 | operand_bit.bind_if_is_defined( | ||
| 326 | mesh_.facets.attributes(), "operand_bit" | ||
| 327 | ); | ||
| 328 | if(operand_bit.is_bound()) { | ||
| 329 | vector<index_t> group_size(32,0); | ||
| 330 | for(index_t f: mesh_.facets) { | ||
| 331 | for(index_t g=0; g<32; ++g) { | ||
| 332 | if(operand_bit[f] == (1u << g)) { | ||
| 333 | ++group_size[g]; | ||
| 334 | } else if (!has_one_bit(operand_bit[f])) { | ||
| 335 | has_facets_in_several_groups = true; | ||
| 336 | break; | ||
| 337 | } | ||
| 338 | } | ||
| 339 | } | ||
| 340 | // If some facets belong to several groups (this | ||
| 341 | // can happen when some input facets are exactly | ||
| 342 | // co-planar), then we do not use group-based | ||
| 343 | // computation (it would be possible, but it | ||
| 344 | // would require to split the set of facets in | ||
| 345 | // several groups according to the groups they | ||
| 346 | // belong to). It is also possible to consider | ||
| 347 | // all the facets in multiple groups as a single | ||
| 348 | // group, but doing so is slow for high-arity | ||
| 349 | // operators (example: dragon_bas.scad). | ||
| 350 | if(!has_facets_in_several_groups) { | ||
| 351 | for(index_t g=0; g<32; ++g) { | ||
| 352 | if(group_size[g] != 0) { | ||
| 353 | nb_groups = std::max(nb_groups,g); | ||
| 354 | } | ||
| 355 | } | ||
| 356 | ++nb_groups; | ||
| 357 | group_ptr.resize(nb_groups+1); | ||
| 358 | group_ptr[0] = 0; | ||
| 359 | for(index_t g=0; g<nb_groups; ++g) { | ||
| 360 | group_ptr[g+1] = group_ptr[g] + group_size[g]; | ||
| 361 | } | ||
| 362 | indices.resize(group_ptr[nb_groups]); | ||
| 363 | vector<index_t> group_insert = group_ptr; | ||
| 364 | for(index_t f: mesh_.facets) { | ||
| 365 | for(index_t g=0; g<32; ++g) { | ||
| 366 | if(operand_bit[f] == (1u << g)) { | ||
| 367 | indices[group_insert[g]] = f; | ||
| 368 | ++group_insert[g]; | ||
| 369 | } | ||
| 370 | } | ||
| 371 | } | ||
| 372 | } | ||
| 373 | } | ||
| 374 | } | ||
| 375 | |||
| 376 | if(nb_groups == 0) { | ||
| 377 | boxes_intersections(boxes, report_BB); | ||
| 378 | } else { | ||
| 379 | for(index_t g1 = 0; g1 < nb_groups; ++g1) { | ||
| 380 | for(index_t g2 = 0; g2 < nb_groups; ++g2) { | ||
| 381 | if(ZE_groups_auto_isect || g1 != g2) { | ||
| 382 | index_t b1 = group_ptr[g1]; | ||
| 383 | index_t e1 = group_ptr[g1+1]; | ||
| 384 | index_t b2 = group_ptr[g2]; | ||
| 385 | index_t e2 = group_ptr[g2+1]; | ||
| 386 | boxes_intersections_hybrid_impl( | ||
| 387 | boxes.data(), | ||
| 388 | indices.data()+b1, indices.data()+e1, | ||
| 389 | boxes.data(), | ||
| 390 | indices.data()+b2, indices.data()+e2, | ||
| 391 | report_BB | ||
| 392 | ); | ||
| 393 | } | ||
| 394 | } | ||
| 395 | } | ||
| 396 | } | ||
| 397 | } | ||
| 398 | |||
| 399 | // Compute facet-facet intersections in parallel | ||
| 400 |
3/6✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
|
60 | Stopwatch* W = new Stopwatch("AABB tri-tri", verbose_); |
| 401 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | Process::spinlock lock = GEOGRAM_SPINLOCK_INIT; |
| 402 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | parallel_for( |
| 403 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
981351 | 0, FF.size(), [&](index_t i) { |
| 404 | TriangleIsects I; | ||
| 405 | 981351 | index_t f1 = FF[i].first; | |
| 406 | 981351 | index_t f2 = FF[i].second; | |
| 407 |
2/2✓ Branch 1 taken 944550 times.
✓ Branch 2 taken 36801 times.
|
981351 | if(!mesh_facets_intersect(mesh_,f1, f2, I)) { |
| 408 | 944550 | return; | |
| 409 | } | ||
| 410 | |||
| 411 | 36801 | Process::acquire_spinlock(lock); | |
| 412 | |||
| 413 |
2/2✓ Branch 0 taken 10907 times.
✓ Branch 1 taken 25894 times.
|
36801 | if(I.size() > 2) { |
| 414 | // Coplanar intersection: to generate the edges, | ||
| 415 | // test validity of all possible pairs of vertices. | ||
| 416 |
2/2✓ Branch 0 taken 41988 times.
✓ Branch 1 taken 10907 times.
|
52895 | for(index_t i1=0; i1< I.size(); ++i1) { |
| 417 |
2/2✓ Branch 0 taken 61612 times.
✓ Branch 1 taken 41988 times.
|
103600 | for(index_t i2=0; i2<i1; ++i2) { |
| 418 | IsectInfo II = { | ||
| 419 | f1, f2, | ||
| 420 | 61612 | I[i1].first, I[i1].second, | |
| 421 | 61612 | I[i2].first, I[i2].second | |
| 422 | 61612 | }; | |
| 423 | |||
| 424 | // Valid edges are the ones where both | ||
| 425 | // extremities are on the same edge | ||
| 426 | // of f1 or on the same edge of f2 | ||
| 427 | // (note: it is a *combinatorial* | ||
| 428 | // convex hull). | ||
| 429 | 61612 | TriangleRegion AB1=regions_convex_hull( | |
| 430 | II.A_rgn_f1,II.B_rgn_f1 | ||
| 431 | ); | ||
| 432 | |||
| 433 | 61612 | TriangleRegion AB2=regions_convex_hull( | |
| 434 | II.A_rgn_f2, II.B_rgn_f2 | ||
| 435 | ); | ||
| 436 | |||
| 437 |
4/4✓ Branch 1 taken 40335 times.
✓ Branch 2 taken 21277 times.
✓ Branch 4 taken 20711 times.
✓ Branch 5 taken 19624 times.
|
61612 | if(region_dim(AB1)==1 || region_dim(AB2)==1) { |
| 438 |
2/2✓ Branch 0 taken 41885 times.
✓ Branch 1 taken 103 times.
|
41988 | intersections.push_back(II); |
| 439 | 41988 | II.flip(); | |
| 440 |
2/2✓ Branch 0 taken 41984 times.
✓ Branch 1 taken 4 times.
|
41988 | intersections.push_back(II); |
| 441 | } | ||
| 442 | } | ||
| 443 | } | ||
| 444 | } else { | ||
| 445 | // Intersection is either a segment | ||
| 446 | // or a vertex of f2. | ||
| 447 | 25894 | TriangleRegion A_rgn_f1 = I[0].first; | |
| 448 | 25894 | TriangleRegion A_rgn_f2 = I[0].second; | |
| 449 | |||
| 450 | TriangleRegion B_rgn_f1 = A_rgn_f1; | ||
| 451 | TriangleRegion B_rgn_f2 = A_rgn_f2; | ||
| 452 | |||
| 453 |
2/2✓ Branch 0 taken 21591 times.
✓ Branch 1 taken 4303 times.
|
25894 | if(I.size() == 2) { |
| 454 | 21591 | B_rgn_f1 = I[1].first; | |
| 455 | 21591 | B_rgn_f2 = I[1].second; | |
| 456 | } | ||
| 457 | |||
| 458 | IsectInfo II = { | ||
| 459 | f1, f2, | ||
| 460 | A_rgn_f1, A_rgn_f2, | ||
| 461 | B_rgn_f1, B_rgn_f2 | ||
| 462 | 25894 | }; | |
| 463 |
2/2✓ Branch 0 taken 25476 times.
✓ Branch 1 taken 418 times.
|
25894 | intersections.push_back(II); |
| 464 | 25894 | II.flip(); | |
| 465 |
2/2✓ Branch 0 taken 25842 times.
✓ Branch 1 taken 52 times.
|
25894 | intersections.push_back(II); |
| 466 | } | ||
| 467 | 36801 | Process::release_spinlock(lock); | |
| 468 | } | ||
| 469 | ); | ||
| 470 | 60 | delete W; | |
| 471 | } | ||
| 472 | 60 | } | |
| 473 | |||
| 474 | 60 | void MeshSurfaceIntersection::intersect_remesh_intersections( | |
| 475 | vector<IsectInfo>& intersections | ||
| 476 | ) { | ||
| 477 |
1/2✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
|
60 | Stopwatch W("CDT",verbose_); |
| 478 | |||
| 479 | // Keep track of original facet ids: they will be copied to the subfacets | ||
| 480 | // whenever a facet is split. | ||
| 481 | // NOTE: of course, facet ids are no longer valids once coplanar facets | ||
| 482 | // are merged. | ||
| 483 |
2/4✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
|
60 | original_facet_id_.bind(mesh_.facets.attributes(), "original_facet_id"); |
| 484 |
2/2✓ Branch 0 taken 114500 times.
✓ Branch 1 taken 60 times.
|
114560 | for(index_t f: mesh_.facets) { |
| 485 | 114500 | original_facet_id_[f] = f; | |
| 486 | } | ||
| 487 | |||
| 488 | // We need to copy the initial mesh, because MeshInTriangle needs | ||
| 489 | // to access it in parallel threads, and without a copy, the internal | ||
| 490 | // arrays of the mesh can be modified whenever there is a | ||
| 491 | // reallocation. Without copying, we would need to insert many | ||
| 492 | // locks (each time the mesh is accessed). | ||
| 493 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | mesh_copy_.copy(mesh_); |
| 494 | |||
| 495 | { | ||
| 496 | // Sort intersections by f1, so that all intersections between f1 | ||
| 497 | // and another facet appear as a contiguous sequence. | ||
| 498 | 60 | GEO::sort( | |
| 499 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | intersections.begin(), intersections.end(), |
| 500 | 60 | [](const IsectInfo& a, const IsectInfo& b) -> bool { | |
| 501 |
14/32✗ 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 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 99939 times.
✓ Branch 17 taken 135667 times.
✓ Branch 18 taken 387914 times.
✓ Branch 19 taken 394731 times.
✓ Branch 20 taken 372655 times.
✓ Branch 21 taken 394731 times.
✓ Branch 22 taken 3219 times.
✓ Branch 23 taken 9383 times.
✓ Branch 24 taken 929 times.
✓ Branch 25 taken 2290 times.
✓ Branch 26 taken 1961 times.
✓ Branch 27 taken 7422 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 41 times.
✓ Branch 31 taken 799 times.
|
1799079 | return (a.f1 < b.f1); |
| 502 | } | ||
| 503 | ); | ||
| 504 | |||
| 505 | // Now iterate on all intersections, and identify | ||
| 506 | // the [b,e[ intervals that correspond to the same f1 facet. | ||
| 507 | // Get starting indices of intersections in same facet. | ||
| 508 | vector<index_t> start; | ||
| 509 | { | ||
| 510 | 60 | index_t b=0; | |
| 511 |
2/2✓ Branch 0 taken 15505 times.
✓ Branch 1 taken 60 times.
|
15625 | while(b < intersections.size()) { |
| 512 | start.push_back(b); | ||
| 513 | index_t e = b; | ||
| 514 | while( | ||
| 515 |
2/2✓ Branch 0 taken 151213 times.
✓ Branch 1 taken 56 times.
|
151269 | e < intersections.size() && |
| 516 |
2/2✓ Branch 0 taken 135764 times.
✓ Branch 1 taken 15449 times.
|
151213 | intersections[e].f1 == intersections[b].f1 |
| 517 | ) { | ||
| 518 | 135764 | ++e; | |
| 519 | } | ||
| 520 | 15505 | b = e; | |
| 521 | } | ||
| 522 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | start.push_back(intersections.size()); |
| 523 | } | ||
| 524 | |||
| 525 | // Display intersection stats | ||
| 526 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 47 times.
|
60 | if(verbose_) { |
| 527 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
|
13 | index_t nb_intersections = intersections.size()/2; |
| 528 | 13 | index_t nb_intersected_triangles = (start.size()-1)/2; | |
| 529 | 13 | index_t max_intersections_in_triangle = 0; | |
| 530 |
2/2✓ Branch 0 taken 4527 times.
✓ Branch 1 taken 13 times.
|
4540 | for(index_t i=0; i+1<start.size(); ++i) { |
| 531 | 4527 | max_intersections_in_triangle = std::max( | |
| 532 | max_intersections_in_triangle, | ||
| 533 | 4527 | index_t(start[i+1]-start[i]) | |
| 534 | ); | ||
| 535 | } | ||
| 536 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
13 | Logger::out("CDT") << "Intersections: " |
| 537 | << nb_intersections << std::endl; | ||
| 538 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
13 | Logger::out("CDT") << "Intersected triangles: " |
| 539 | << nb_intersected_triangles | ||
| 540 | << std::endl; | ||
| 541 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
26 | Logger::out("CDT") << "Max intersections in triangle: " |
| 542 |
1/2✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
|
13 | << max_intersections_in_triangle |
| 543 | << std::endl; | ||
| 544 | } | ||
| 545 | |||
| 546 | |||
| 547 | #define TRIANGULATE_IN_PARALLEL | ||
| 548 | |||
| 549 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
|
60 | index_t f_done = 0; |
| 550 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | index_t f_tot = (start.size()-1); |
| 551 | |||
| 552 | #ifdef TRIANGULATE_IN_PARALLEL | ||
| 553 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | parallel_for_slice( |
| 554 |
1/4✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
60 | 0,start.size()-1, [&](index_t k1, index_t k2) { |
| 555 | #else | ||
| 556 | index_t k1 = 0; | ||
| 557 | index_t k2 = start.size()-1; | ||
| 558 | #endif | ||
| 559 | |||
| 560 | 228 | MeshInTriangle MIT(*this); | |
| 561 |
1/2✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
|
228 | MIT.set_delaunay(delaunay_); |
| 562 |
1/2✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
|
228 | MIT.set_dry_run(dry_run_); |
| 563 | |||
| 564 |
1/2✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
|
228 | index_t tid = Thread::current_id(); |
| 565 | |||
| 566 |
2/2✓ Branch 0 taken 15505 times.
✓ Branch 1 taken 228 times.
|
15733 | for(index_t k=k1; k<k2; ++k) { |
| 567 |
2/2✓ Branch 0 taken 4527 times.
✓ Branch 1 taken 10978 times.
|
15505 | index_t b = start[k]; |
| 568 | 15505 | index_t e = start[k+1]; | |
| 569 | |||
| 570 |
2/2✓ Branch 0 taken 4527 times.
✓ Branch 1 taken 10978 times.
|
15505 | if(fine_verbose_) { |
| 571 | 4527 | ++f_done; | |
| 572 |
2/4✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4527 times.
✗ Branch 5 not taken.
|
9054 | Logger::out("Isect") |
| 573 |
1/4✓ Branch 0 taken 4527 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
4527 | << String::format( |
| 574 | "[%2d] %5d/%5d %6d:%3d", | ||
| 575 |
1/2✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
|
4527 | int(tid), int(f_done), int(f_tot), |
| 576 |
1/2✓ Branch 1 taken 4527 times.
✗ Branch 2 not taken.
|
4527 | int(intersections[b].f1), int(e-b) |
| 577 | ) | ||
| 578 | << std::endl; | ||
| 579 | } | ||
| 580 | |||
| 581 |
1/2✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
|
15505 | MIT.begin_facet(intersections[b].f1); |
| 582 |
2/2✓ Branch 0 taken 135764 times.
✓ Branch 1 taken 15505 times.
|
151269 | for(index_t i=b; i<e; ++i) { |
| 583 |
2/2✓ Branch 0 taken 54380 times.
✓ Branch 1 taken 81384 times.
|
135764 | const IsectInfo& II = intersections[i]; |
| 584 | |||
| 585 | // Each IsectInfo is either an individual vertex | ||
| 586 | // or a segment with two vertices. | ||
| 587 | // Each vertex is represented combinatorially. | ||
| 588 | // The MeshInTriangle knows how to compute the | ||
| 589 | // geometry from the combinatorial information. | ||
| 590 | |||
| 591 | if(II.is_point()) { | ||
| 592 | 8606 | MIT.add_vertex( | |
| 593 |
1/2✓ Branch 1 taken 8606 times.
✗ Branch 2 not taken.
|
8606 | II.f2, |
| 594 | II.A_rgn_f1, II.A_rgn_f2 | ||
| 595 | ); | ||
| 596 | } else { | ||
| 597 | 127158 | MIT.add_edge( | |
| 598 | 127158 | II.f2, | |
| 599 | 127158 | II.A_rgn_f1, II.A_rgn_f2, | |
| 600 |
1/2✓ Branch 1 taken 127158 times.
✗ Branch 2 not taken.
|
127158 | II.B_rgn_f1, II.B_rgn_f2 |
| 601 | ); | ||
| 602 | } | ||
| 603 | } | ||
| 604 | |||
| 605 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15505 times.
|
15505 | if(e-b >= monster_threshold_) { |
| 606 | ✗ | index_t f = intersections[b].f1; | |
| 607 | ✗ | MIT.save_constraints( | |
| 608 | ✗ | "constraints_"+String::to_string(f)+".geogram" | |
| 609 | ); | ||
| 610 | } | ||
| 611 | |||
| 612 | // Inserts constraints | ||
| 613 | // and creates new vertices in shared mesh | ||
| 614 |
1/2✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
|
15505 | MIT.commit(); |
| 615 | |||
| 616 | // For debugging, optionally save "monsters" | ||
| 617 | // (that is, triangles that have a huge number | ||
| 618 | // of intersections). | ||
| 619 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15505 times.
|
15505 | if(e-b >= monster_threshold_) { |
| 620 | ✗ | index_t f = intersections[b].f1; | |
| 621 | ✗ | MIT.save( | |
| 622 | ✗ | "triangulation_"+String::to_string(f)+".geogram" | |
| 623 | ); | ||
| 624 | //MIT.save_constraints( | ||
| 625 | // "constraints_"+String::to_string(f)+".geogram" | ||
| 626 | //); | ||
| 627 | std::ofstream out( | ||
| 628 | ✗ | "facet_"+String::to_string(f)+".obj" | |
| 629 | ✗ | ); | |
| 630 | ✗ | for(index_t v: mesh_copy_.facets.vertices(f)) { | |
| 631 | ✗ | vec3 p = mesh_copy_.vertices.point(v); | |
| 632 | ✗ | out << "v " << p << std::endl; | |
| 633 | } | ||
| 634 | ✗ | out << "f "; | |
| 635 | for( | ||
| 636 | index_t lv=0; | ||
| 637 | ✗ | lv<mesh_copy_.facets.nb_vertices(f); ++lv | |
| 638 | ) { | ||
| 639 | ✗ | out << lv+1 << " "; | |
| 640 | } | ||
| 641 | out << std::endl; | ||
| 642 | ✗ | } | |
| 643 | |||
| 644 | // Clear it so that it is clean for next triangle. | ||
| 645 |
1/2✓ Branch 1 taken 15505 times.
✗ Branch 2 not taken.
|
15505 | MIT.clear(); |
| 646 | } | ||
| 647 |
2/2✓ Branch 0 taken 52 times.
✓ Branch 1 taken 176 times.
|
228 | if(fine_verbose_) { |
| 648 |
2/4✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
|
104 | Logger::out("Isect") |
| 649 |
2/6✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 52 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
104 | << String::format("[%2d] done",int(tid)) |
| 650 | << std::endl; | ||
| 651 | } | ||
| 652 | #ifdef TRIANGULATE_IN_PARALLEL | ||
| 653 | 228 | }); | |
| 654 | #endif | ||
| 655 | } | ||
| 656 | 60 | } | |
| 657 | |||
| 658 | 60 | void MeshSurfaceIntersection::intersect_epilogue( | |
| 659 | const vector<IsectInfo>& intersections | ||
| 660 | ) { | ||
| 661 |
1/2✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
|
60 | Stopwatch W_epilogue("Epilogue", verbose_); |
| 662 | // Vertices coming from intersections may land exactly | ||
| 663 | // on an existing vertex (see #111) | ||
| 664 | { | ||
| 665 |
2/4✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
|
60 | Stopwatch W("I on v", verbose_); |
| 666 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | vector<index_t> v2v(mesh_.vertices.nb()); |
| 667 |
2/2✓ Branch 0 taken 82496 times.
✓ Branch 1 taken 60 times.
|
82556 | for(index_t v : mesh_.vertices) { |
| 668 | 82496 | v2v[v] = v; | |
| 669 | } | ||
| 670 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | parallel_for( |
| 671 | 0, mesh_.vertices.nb(), | ||
| 672 | ✗ | [&](index_t v) { | |
| 673 | // If the vertex is an original vertex, not | ||
| 674 | // coming from an intersection, check whether | ||
| 675 | // it also exists as an intersection | ||
| 676 |
2/2✓ Branch 0 taken 57599 times.
✓ Branch 1 taken 24897 times.
|
82496 | if(vertex_to_exact_point_[v] == nullptr) { |
| 677 | 57599 | const vec3& xyz = mesh_.vertices.point(v); | |
| 678 | 57599 | ExactPoint K(xyz[0], xyz[1], xyz[2], 1.0); | |
| 679 | auto it = exact_point_to_vertex_.find(K); | ||
| 680 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 57593 times.
|
57599 | if(it != exact_point_to_vertex_.end()) { |
| 681 | 6 | v2v[v] = it->second; | |
| 682 | } | ||
| 683 | 57599 | } | |
| 684 | 82496 | } | |
| 685 | ); | ||
| 686 |
2/2✓ Branch 0 taken 737661 times.
✓ Branch 1 taken 60 times.
|
737721 | for(index_t c : mesh_.facet_corners) { |
| 687 | 737661 | index_t v = v2v[mesh_.facet_corners.vertex(c)]; | |
| 688 | mesh_.facet_corners.set_vertex(c, v); | ||
| 689 | } | ||
| 690 | 60 | } | |
| 691 | |||
| 692 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
|
60 | if(interpolate_attributes_) { |
| 693 | ✗ | Stopwatch W("Attrib", verbose_); | |
| 694 | Attribute<index_t> original_facet_id( | ||
| 695 | ✗ | mesh_.facets.attributes(), "original_facet_id" | |
| 696 | ); | ||
| 697 | |||
| 698 | ✗ | for(index_t c: mesh_.facet_corners) { | |
| 699 | ✗ | index_t f = c/3; | |
| 700 | ✗ | index_t f0 = original_facet_id[f]; | |
| 701 | // If this is an original facet, there is nothing to do ! | ||
| 702 | ✗ | if(f == f0) { | |
| 703 | ✗ | continue; | |
| 704 | } | ||
| 705 | |||
| 706 | ✗ | index_t v = mesh_.facet_corners.vertex(c); | |
| 707 | index_t v1 = mesh_.facets.vertex(f0,0); | ||
| 708 | index_t v2 = mesh_.facets.vertex(f0,1); | ||
| 709 | index_t v3 = mesh_.facets.vertex(f0,2); | ||
| 710 | |||
| 711 | index_t c1 = mesh_.facets.corner(f0,0); | ||
| 712 | index_t c2 = mesh_.facets.corner(f0,1); | ||
| 713 | index_t c3 = mesh_.facets.corner(f0,2); | ||
| 714 | |||
| 715 | ✗ | if(v == v1) { | |
| 716 | ✗ | mesh_.facet_corners.attributes().copy_item(c,c1); | |
| 717 | ✗ | continue; | |
| 718 | } | ||
| 719 | |||
| 720 | ✗ | if(v == v2) { | |
| 721 | ✗ | mesh_.facet_corners.attributes().copy_item(c,c2); | |
| 722 | ✗ | continue; | |
| 723 | } | ||
| 724 | |||
| 725 | ✗ | if(v == v3) { | |
| 726 | ✗ | mesh_.facet_corners.attributes().copy_item(c,c3); | |
| 727 | ✗ | continue; | |
| 728 | } | ||
| 729 | |||
| 730 | ✗ | vec3 p1=mesh_.vertices.point(v1); | |
| 731 | ✗ | vec3 p2=mesh_.vertices.point(v2); | |
| 732 | ✗ | vec3 p3=mesh_.vertices.point(v3); | |
| 733 | ✗ | vec3 p =mesh_.vertices.point(v); | |
| 734 | |||
| 735 | double a = Geom::triangle_area(p1,p2,p3); | ||
| 736 | double a1 = Geom::triangle_area(p ,p2,p3); | ||
| 737 | double a2 = Geom::triangle_area(p1,p ,p3); | ||
| 738 | double a3 = Geom::triangle_area(p1,p2,p ); | ||
| 739 | |||
| 740 | ✗ | mesh_.facet_corners.attributes().zero_item(c); | |
| 741 | ✗ | mesh_.facet_corners.attributes().madd_item(c, a1/a, c1); | |
| 742 | ✗ | mesh_.facet_corners.attributes().madd_item(c, a2/a, c2); | |
| 743 | ✗ | mesh_.facet_corners.attributes().madd_item(c, a3/a, c3); | |
| 744 | } | ||
| 745 | ✗ | } | |
| 746 | |||
| 747 | |||
| 748 | // Remove original facets that have intersections. | ||
| 749 | { | ||
| 750 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | vector<index_t> has_intersections(mesh_.facets.nb(), 0); |
| 751 |
2/2✓ Branch 0 taken 135764 times.
✓ Branch 1 taken 60 times.
|
135824 | for(const IsectInfo& II: intersections) { |
| 752 | 135764 | has_intersections[II.f1] = 1; | |
| 753 | 135764 | has_intersections[II.f2] = 1; | |
| 754 | } | ||
| 755 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | mesh_.facets.delete_elements(has_intersections); |
| 756 | } | ||
| 757 | |||
| 758 | // There can be duplicated facets coming from | ||
| 759 | // tesselated co-planar facets. | ||
| 760 | // Note: this updates operand_bit attribute | ||
| 761 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | mesh_remove_bad_facets_no_check(mesh_); |
| 762 | |||
| 763 | #ifdef MESH_SURFACE_INTERSECTION_DEBUG | ||
| 764 | std::cerr << "Sanity check: verify that there is no degenerate triangle" | ||
| 765 | << std::endl; | ||
| 766 | // Sanity check: do we have facets with their three vertices | ||
| 767 | // aligned ? Normally cannot happen since we have eliminated | ||
| 768 | // them during intersection, but who knows ? | ||
| 769 | // Actually this happens sometimes... | ||
| 770 | { | ||
| 771 | Attribute<bool> selected(mesh_.facets.attributes(), "selection"); | ||
| 772 | for(index_t t: mesh_.facets) { | ||
| 773 | if(PCK::aligned_3d( | ||
| 774 | exact_vertex(mesh_.facets.vertex(t,0)), | ||
| 775 | exact_vertex(mesh_.facets.vertex(t,1)), | ||
| 776 | exact_vertex(mesh_.facets.vertex(t,2)) | ||
| 777 | )) { | ||
| 778 | selected[t] = true; | ||
| 779 | std::cerr << "FACET HAS 3 ALIGNED VERTICES" << std::endl; | ||
| 780 | } else { | ||
| 781 | selected[t] = false; | ||
| 782 | } | ||
| 783 | |||
| 784 | } | ||
| 785 | } | ||
| 786 | std::cerr << "There is no degenerate triangle" << std::endl; | ||
| 787 | #endif | ||
| 788 | |||
| 789 |
2/2✓ Branch 0 taken 57 times.
✓ Branch 1 taken 3 times.
|
60 | if(use_radial_sort_) { |
| 790 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | build_Weiler_model(); |
| 791 | } | ||
| 792 | |||
| 793 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | PCK::set_SOS_mode(SOS_bkp_); |
| 794 | 60 | } | |
| 795 | |||
| 796 | 60 | void MeshSurfaceIntersection::intersect() { | |
| 797 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
|
60 | if(mesh_.facets.nb() == 0) { |
| 798 | ✗ | return; | |
| 799 | } | ||
| 800 |
1/2✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
|
60 | Stopwatch W("Intersect", verbose_); |
| 801 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | intersect_prologue(); |
| 802 | vector<IsectInfo> intersections; | ||
| 803 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | intersect_get_intersections(intersections); |
| 804 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | intersect_remesh_intersections(intersections); |
| 805 |
1/2✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
|
60 | intersect_epilogue(intersections); |
| 806 | 60 | } | |
| 807 | |||
| 808 |
2/2✓ Branch 0 taken 150739 times.
✓ Branch 1 taken 657051 times.
|
807790 | MeshSurfaceIntersection::ExactPoint MeshSurfaceIntersection::exact_vertex( |
| 809 | index_t v | ||
| 810 | ) const { | ||
| 811 | geo_debug_assert(v < mesh_.vertices.nb()); | ||
| 812 | 807790 | const ExactPoint* p = vertex_to_exact_point_[v]; | |
| 813 |
2/2✓ Branch 0 taken 150739 times.
✓ Branch 1 taken 657051 times.
|
807790 | if(p != nullptr) { |
| 814 | 150739 | return *p; | |
| 815 | } | ||
| 816 | 657051 | const vec3& xyz = mesh_.vertices.point(v); | |
| 817 | 657051 | return ExactPoint(xyz[0], xyz[1], xyz[2], 1.0); | |
| 818 | } | ||
| 819 | |||
| 820 | 91237 | index_t MeshSurfaceIntersection::find_or_create_exact_vertex( | |
| 821 | const ExactPoint& p | ||
| 822 | ) { | ||
| 823 | std::map<ExactPoint,index_t,ExactPointCompare>::iterator it; | ||
| 824 | bool inserted; | ||
| 825 |
1/2✓ Branch 2 taken 91237 times.
✗ Branch 3 not taken.
|
91237 | std::tie(it, inserted) = exact_point_to_vertex_.insert( |
| 826 | 91237 | std::make_pair(p,NO_INDEX) | |
| 827 | ); | ||
| 828 |
2/2✓ Branch 0 taken 66340 times.
✓ Branch 1 taken 24897 times.
|
91237 | if(!inserted) { |
| 829 | 66340 | return it->second; | |
| 830 | } | ||
| 831 | 24897 | vec3 p_inexact = PCK::approximate(p); | |
| 832 | 24897 | index_t v = mesh_.vertices.create_vertex(p_inexact.data()); | |
| 833 | 24897 | it->second = v; | |
| 834 | 24897 | vertex_to_exact_point_[v] = &(it->first); | |
| 835 | 24897 | return v; | |
| 836 | } | ||
| 837 | |||
| 838 | /************************ Radial sort ***********************************/ | ||
| 839 | |||
| 840 | 3820 | void MeshSurfaceIntersection::RadialSort::init(index_t h_ref) { | |
| 841 | 3820 | degenerate_ = false; | |
| 842 | 3820 | h_ref_ = NO_INDEX; // so that normal() computes N_ref_ | |
| 843 | 3820 | N_ref_ = normal(h_ref); | |
| 844 | 3820 | h_ref_ = h_ref; | |
| 845 | 3820 | } | |
| 846 | |||
| 847 | 24513 | bool MeshSurfaceIntersection::RadialSort::operator()( | |
| 848 | index_t h1, index_t h2 | ||
| 849 | ) const { | ||
| 850 | |||
| 851 | // To order facets around radial edges, we use two different things: | ||
| 852 | // | ||
| 853 | // 1) a 2D coordinate system. For each halfedge h, this gives | ||
| 854 | // two signs. The 2D coordinate system is defined by: | ||
| 855 | // - a reference halfedge h_ref passed to init(h_ref) | ||
| 856 | // - u_sign (su): location of opposite vertex relative to facet | ||
| 857 | // incident to h_ref: h_orient(h_ref_, h) | ||
| 858 | // - v_sign (sv): dot product between the normals to the facets | ||
| 859 | // incident to h and incident to h_ref: h_refNorient(h) | ||
| 860 | // u_sign and v_sign determine 4 quadrants (in fact 8 regions if one | ||
| 861 | // counts when they vanish). If two halfedges land in two different | ||
| 862 | // regions then their order is determined. The order of the regions | ||
| 863 | // is stored in su_sv_to_linear_index[3][3]. | ||
| 864 | // | ||
| 865 | // 2) now if two halfedges h1,h2 land in the same region, their order is | ||
| 866 | // given by orient3d(p1,p2,q1,q2) = horient(h1,h2) where p1,p2 are | ||
| 867 | // the vertices of the radial edge and q1 (resp q2) the opposite | ||
| 868 | // vertices of h1 (resp h2). | ||
| 869 | |||
| 870 | 24513 | Sign su1 = h_orient(h_ref_, h1); | |
| 871 | 24513 | Sign su2 = h_orient(h_ref_, h2); | |
| 872 | |||
| 873 | // Optimization: quick exit, h1 appears first if it is the | ||
| 874 | // only one in positive u half. | ||
| 875 |
2/2✓ Branch 0 taken 2121 times.
✓ Branch 1 taken 22392 times.
|
24513 | if(su1 * su2 < 0) { |
| 876 | 2121 | return (su1 > 0); | |
| 877 | } | ||
| 878 | |||
| 879 | 22392 | Sign sv1 = h_refNorient(h1); | |
| 880 | 22392 | Sign sv2 = h_refNorient(h2); | |
| 881 | |||
| 882 | // Map (su,sv) sign pair into a region linear index ("pseudo-angle") | ||
| 883 | // | ||
| 884 | // su | ||
| 885 | // ^ | ||
| 886 | // | + 2 | ||
| 887 | // | 3 | 1 | ||
| 888 | // | 0 4---+---0 | ||
| 889 | // | 5 | 7 | ||
| 890 | // | - 6 | ||
| 891 | // | | ||
| 892 | // | - 0 + | ||
| 893 | // *--------------> sv | ||
| 894 | static int su_sv_to_linear_index[3][3] = { | ||
| 895 | // -, 0, + <--------- su = h_orient(h_ref, h) | ||
| 896 | {5, 4, 3}, // - <--- sv = h_Norient(h) | ||
| 897 | {6,-1, 2}, // 0 | ||
| 898 | {7, 0, 1} // + | ||
| 899 | }; | ||
| 900 | |||
| 901 | 22392 | int theta1 = su_sv_to_linear_index[sv1+1][su1+1]; | |
| 902 | 22392 | int theta2 = su_sv_to_linear_index[sv2+1][su2+1]; | |
| 903 | |||
| 904 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22392 times.
|
22392 | if(theta1 == -1 || theta2 == -1) { |
| 905 | ✗ | report_problem("Triangle with both zero orient and zero Norient"); | |
| 906 | ✗ | return false; | |
| 907 | } | ||
| 908 | |||
| 909 |
2/2✓ Branch 0 taken 22390 times.
✓ Branch 1 taken 2 times.
|
22392 | if(theta1 != theta2) { // Different thetas: then we know the order |
| 910 | 22390 | return (theta2 > theta1); | |
| 911 | } | ||
| 912 | |||
| 913 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if((theta1 & 1) == 0) { // Same thetas, = to 0,2,4,6 (should not happen) |
| 914 | ✗ | report_problem("Both triangles in same reference half-plane"); | |
| 915 | ✗ | return false; | |
| 916 | } | ||
| 917 | |||
| 918 | // Same thetas, both = to 1,3,5 or 7, measure relative orientation | ||
| 919 | 2 | Sign o_12 = h_orient(h1,h2); | |
| 920 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(o_12 == ZERO) { |
| 921 | ✗ | report_problem("Both triangles in same half-plane"); | |
| 922 | ✗ | return false; | |
| 923 | } | ||
| 924 | 2 | return o_12 > 0; | |
| 925 | } | ||
| 926 | |||
| 927 | 49028 | Sign MeshSurfaceIntersection::RadialSort::h_orient( | |
| 928 | index_t h1, index_t h2 | ||
| 929 | ) const { | ||
| 930 | |||
| 931 | // Normally we compute -orient3d(p1,p2,q1,q2) where | ||
| 932 | // p1,p2 denote the radial edge and q1 (resp q2) the opposite | ||
| 933 | // vertex to h1 (resp h2). | ||
| 934 | // Replacing a facet of a tetrahedron with a co-planar triangle | ||
| 935 | // of same orientation does not change signed volume, so instead | ||
| 936 | // of the facet adjacent to h1, we use the facet of the original | ||
| 937 | // mesh (that contains the facet adjacent to h1), that has simpler | ||
| 938 | // coordinates. | ||
| 939 | |||
| 940 |
2/2✓ Branch 0 taken 32427 times.
✓ Branch 1 taken 16601 times.
|
49028 | if(h1 == h2) { |
| 941 | return ZERO; | ||
| 942 | } | ||
| 943 |
1/2✓ Branch 0 taken 32427 times.
✗ Branch 1 not taken.
|
32427 | index_t f1 = I_.halfedges_.facet(h1); |
| 944 | index_t f2 = I_.halfedges_.facet(h2); | ||
| 945 | index_t w1 = I_.halfedges_.vertex(h1,2); | ||
| 946 | index_t w2 = I_.halfedges_.vertex(h2,2); | ||
| 947 | |||
| 948 | // Optimization (using original points as often as possible) | ||
| 949 | // If w1 is an original vertex, use w1 and original facet of h2 | ||
| 950 | // (there is no minus sign, because args are swapped, we have h2 first) | ||
| 951 |
2/2✓ Branch 0 taken 5353 times.
✓ Branch 1 taken 27074 times.
|
32427 | if(I_.is_original_vertex(w1)) { |
| 952 | 5353 | vec3 p0 = mesh_.vertices.point(w1); | |
| 953 | 5353 | auto [q0, q1, q2] = I_.get_initial_facet_vertices(f2); | |
| 954 | return Sign(PCK::orient_3d(q0,q1,q2,p0)); | ||
| 955 | } | ||
| 956 | |||
| 957 | // Optimization (using original points as often as possible) | ||
| 958 | // If w2 is an original vertex, use w2 and original facet of h1 | ||
| 959 |
2/2✓ Branch 0 taken 3493 times.
✓ Branch 1 taken 23581 times.
|
27074 | if(I_.is_original_vertex(w2)) { |
| 960 | 3493 | vec3 q0 = mesh_.vertices.point(w2); | |
| 961 | 3493 | auto [p0, p1, p2] = I_.get_initial_facet_vertices(f1); | |
| 962 | 3493 | return Sign(-PCK::orient_3d(p0,p1,p2,q0)); | |
| 963 | } | ||
| 964 | |||
| 965 | // General case: use w2 (it's an intersection point) and original | ||
| 966 | // facet of h1 converted to exact points. | ||
| 967 | 23581 | auto [p0, p1, p2] = I_.get_initial_facet_vertices(f1); | |
| 968 | |||
| 969 | 23581 | ExactPoint pp0(p0.x, p0.y, p0.z, 1.0); | |
| 970 | 23581 | ExactPoint pp1(p1.x, p1.y, p1.z, 1.0); | |
| 971 | 23581 | ExactPoint pp2(p2.x, p2.y, p2.z, 1.0); | |
| 972 |
1/2✓ Branch 1 taken 23581 times.
✗ Branch 2 not taken.
|
23581 | const ExactPoint& q2 = I_.exact_vertex(w2); |
| 973 | |||
| 974 |
1/2✓ Branch 1 taken 23581 times.
✗ Branch 2 not taken.
|
23581 | return Sign(-PCK::orient_3d(pp0,pp1,pp2,q2)); |
| 975 | |||
| 976 | 23581 | } | |
| 977 | |||
| 978 | 44784 | Sign MeshSurfaceIntersection::RadialSort::h_refNorient(index_t h2) const { | |
| 979 |
2/2✓ Branch 0 taken 28183 times.
✓ Branch 1 taken 16601 times.
|
44784 | if(h2 == h_ref_) { |
| 980 | return POSITIVE; | ||
| 981 | } | ||
| 982 |
3/4✓ Branch 0 taken 33 times.
✓ Branch 1 taken 28150 times.
✓ Branch 3 taken 33 times.
✗ Branch 4 not taken.
|
28183 | static PCK::PredicateStats stats("h_refNorient"); |
| 983 | stats.log_invoke(); | ||
| 984 | 28183 | exact::vec3 N2 = normal(h2); | |
| 985 |
1/2✓ Branch 1 taken 28183 times.
✗ Branch 2 not taken.
|
28183 | return dot(N_ref_,N2).sign(); |
| 986 | 28183 | } | |
| 987 | |||
| 988 | 32003 | exact::vec3 MeshSurfaceIntersection::RadialSort::normal(index_t h) const { | |
| 989 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32003 times.
|
32003 | if(h == h_ref_) { |
| 990 | ✗ | return N_ref_; | |
| 991 | } | ||
| 992 | 32003 | auto [p1, p2, p3] = I_.get_initial_facet_vertices( | |
| 993 | I_.halfedges_.facet(h) | ||
| 994 | 32003 | ); | |
| 995 | exact::vec3 N = cross( | ||
| 996 | 64006 | make_vec3<exact::vec3>(p1,p2), | |
| 997 |
1/2✓ Branch 3 taken 32003 times.
✗ Branch 4 not taken.
|
64006 | make_vec3<exact::vec3>(p1,p3) |
| 998 |
1/2✓ Branch 1 taken 32003 times.
✗ Branch 2 not taken.
|
32003 | ); |
| 999 | Numeric::optimize_number_representation(N); | ||
| 1000 | return N; | ||
| 1001 | 32003 | } | |
| 1002 | |||
| 1003 | ✗ | void MeshSurfaceIntersection::RadialSort::report_problem( | |
| 1004 | const char* message | ||
| 1005 | ) const { | ||
| 1006 | ✗ | degenerate_ = true; | |
| 1007 | ✗ | Logger::err("RadialSort") << message << std::endl; | |
| 1008 | ✗ | } | |
| 1009 | |||
| 1010 | /*****************************************************************/ | ||
| 1011 | |||
| 1012 | 27 | void MeshSurfaceIntersection::mark_external_shell( | |
| 1013 | vector<index_t>& on_external_shell | ||
| 1014 | ) { | ||
| 1015 | // Chart attribute corresponds to volumetric regions | ||
| 1016 |
1/2✓ Branch 2 taken 27 times.
✗ Branch 3 not taken.
|
27 | Attribute<index_t> chart(mesh_.facets.attributes(), "chart"); |
| 1017 | |||
| 1018 | // Get nb charts | ||
| 1019 | 27 | index_t nb_charts = 0; | |
| 1020 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 27 times.
|
217093 | for(index_t f: mesh_.facets) { |
| 1021 | 217066 | nb_charts = std::max(nb_charts, chart[f]+1); | |
| 1022 | } | ||
| 1023 | |||
| 1024 | // Get connected components by traversing both alpha2 and alpha3 links, | ||
| 1025 | // and orient facets coherently | ||
| 1026 | index_t nb_components = 0; | ||
| 1027 | vector<index_t> component(mesh_.facets.nb(), NO_INDEX); | ||
| 1028 | { | ||
| 1029 |
4/4✓ Branch 0 taken 30 times.
✓ Branch 1 taken 217036 times.
✓ Branch 2 taken 217066 times.
✓ Branch 3 taken 27 times.
|
217093 | for(index_t f:mesh_.facets) { |
| 1030 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 217036 times.
|
217066 | if(component[f] == NO_INDEX) { |
| 1031 | std::stack<index_t> S; | ||
| 1032 |
1/2✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
|
30 | component[f] = nb_components; |
| 1033 | S.push(f); | ||
| 1034 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 30 times.
|
217096 | while(!S.empty()) { |
| 1035 |
2/2✓ Branch 0 taken 213745 times.
✓ Branch 1 taken 3321 times.
|
217066 | index_t f1 = S.top(); |
| 1036 | S.pop(); | ||
| 1037 | |||
| 1038 | { | ||
| 1039 |
2/2✓ Branch 0 taken 69653 times.
✓ Branch 1 taken 147413 times.
|
217066 | index_t f2 = halfedges_.facet_alpha3(f1); |
| 1040 |
2/2✓ Branch 0 taken 69653 times.
✓ Branch 1 taken 147413 times.
|
217066 | if(component[f2] == NO_INDEX) { |
| 1041 |
2/2✓ Branch 0 taken 68589 times.
✓ Branch 1 taken 1064 times.
|
69653 | component[f2]=component[f1]; |
| 1042 | S.push(f2); | ||
| 1043 | } | ||
| 1044 | } | ||
| 1045 | |||
| 1046 |
2/2✓ Branch 0 taken 651198 times.
✓ Branch 1 taken 217066 times.
|
868264 | for(index_t le1=0; le1<3; ++le1) { |
| 1047 |
1/2✓ Branch 0 taken 651198 times.
✗ Branch 1 not taken.
|
651198 | index_t f2 = mesh_.facets.adjacent(f1,le1); |
| 1048 |
3/4✓ Branch 0 taken 651198 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 147383 times.
✓ Branch 3 taken 503815 times.
|
651198 | if(f2!=NO_INDEX && component[f2]==NO_INDEX) { |
| 1049 | #ifdef GEO_DEBUG | ||
| 1050 | index_t le2 = mesh_.facets.find_adjacent(f2,f1); | ||
| 1051 | geo_debug_assert( | ||
| 1052 | mesh_.facets.vertex(f1,le1) != | ||
| 1053 | mesh_.facets.vertex(f2,le2) | ||
| 1054 | ); | ||
| 1055 | #endif | ||
| 1056 |
2/2✓ Branch 0 taken 145126 times.
✓ Branch 1 taken 2257 times.
|
147383 | component[f2] = component[f1]; |
| 1057 | S.push(f2); | ||
| 1058 | } | ||
| 1059 | } | ||
| 1060 | } | ||
| 1061 | } | ||
| 1062 | 217066 | ++nb_components; | |
| 1063 | } | ||
| 1064 | } | ||
| 1065 | |||
| 1066 | // Compute the volume enclosed by each chart | ||
| 1067 | |||
| 1068 |
1/4✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
27 | vector<double> chart_volume(nb_charts,0.0); |
| 1069 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 27 times.
|
217093 | for(index_t f: mesh_.facets) { |
| 1070 | 217066 | vec3 p1 = mesh_.facets.point(f,0); | |
| 1071 | 217066 | vec3 p2 = mesh_.facets.point(f,1); | |
| 1072 | 217066 | vec3 p3 = mesh_.facets.point(f,2); | |
| 1073 | 217066 | chart_volume[chart[f]] += dot(p1,cross(p2,p3)) / 6.0; | |
| 1074 | } | ||
| 1075 | |||
| 1076 |
2/2✓ Branch 0 taken 1122 times.
✓ Branch 1 taken 27 times.
|
1149 | for(index_t c=0; c<chart_volume.size(); ++c) { |
| 1077 | 1122 | chart_volume[c] = ::fabs(chart_volume[c]); | |
| 1078 | } | ||
| 1079 | |||
| 1080 | // For each component, find the chart that encloses the largest | ||
| 1081 | // volume (it is the external boundary of the component) | ||
| 1082 | |||
| 1083 |
2/6✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
27 | vector<double> max_chart_volume_in_component(nb_components, 0.0); |
| 1084 | vector<index_t> chart_with_max_volume_in_component( | ||
| 1085 | nb_components, NO_INDEX | ||
| 1086 | ); | ||
| 1087 | |||
| 1088 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 27 times.
|
217093 | for(index_t f: mesh_.facets) { |
| 1089 |
2/2✓ Branch 0 taken 104878 times.
✓ Branch 1 taken 112188 times.
|
217066 | double V = chart_volume[chart[f]]; |
| 1090 |
2/2✓ Branch 0 taken 104878 times.
✓ Branch 1 taken 112188 times.
|
217066 | if( V >= max_chart_volume_in_component[component[f]]) { |
| 1091 | 104878 | max_chart_volume_in_component[component[f]] = V; | |
| 1092 | 104878 | chart_with_max_volume_in_component[component[f]] = chart[f]; | |
| 1093 | } | ||
| 1094 | } | ||
| 1095 | |||
| 1096 |
1/2✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
|
27 | on_external_shell.resize(mesh_.facets.nb()); |
| 1097 |
2/2✓ Branch 0 taken 217066 times.
✓ Branch 1 taken 27 times.
|
217093 | for(index_t f: mesh_.facets) { |
| 1098 | 217066 | on_external_shell[f] = ( | |
| 1099 | 217066 | chart[f] == chart_with_max_volume_in_component[component[f]] | |
| 1100 | ); | ||
| 1101 | } | ||
| 1102 | 27 | } | |
| 1103 | |||
| 1104 | /*************************************************************************/ | ||
| 1105 | |||
| 1106 | 57 | void MeshSurfaceIntersection::build_Weiler_model() { | |
| 1107 | |||
| 1108 |
1/2✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
|
57 | Stopwatch Wweiler("Weiler",verbose_); |
| 1109 | |||
| 1110 | Attribute<bool> corner_is_on_border( | ||
| 1111 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | mesh_.facet_corners.attributes(), "is_on_border" |
| 1112 | ); | ||
| 1113 |
2/2✓ Branch 0 taken 589587 times.
✓ Branch 1 taken 57 times.
|
589644 | for(index_t c: mesh_.facet_corners) { |
| 1114 | corner_is_on_border[c] = false; | ||
| 1115 | } | ||
| 1116 | |||
| 1117 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | halfedges_.initialize(); |
| 1118 | |||
| 1119 | // Step 1: Duplicate all surfaces and create alpha3 links | ||
| 1120 | { | ||
| 1121 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | index_t nf = mesh_.facets.nb(); |
| 1122 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | mesh_.facets.create_triangles(nf); |
| 1123 |
2/2✓ Branch 0 taken 196529 times.
✓ Branch 1 taken 57 times.
|
196586 | for(index_t f1=0; f1<nf; ++f1) { |
| 1124 | 196529 | index_t f2 = f1+nf; | |
| 1125 |
1/2✓ Branch 0 taken 196529 times.
✗ Branch 1 not taken.
|
196529 | mesh_.facets.set_vertex(f2,0,mesh_.facets.vertex(f1,2)); |
| 1126 | mesh_.facets.set_vertex(f2,1,mesh_.facets.vertex(f1,1)); | ||
| 1127 | mesh_.facets.set_vertex(f2,2,mesh_.facets.vertex(f1,0)); | ||
| 1128 | |||
| 1129 | // Copy attributes | ||
| 1130 |
1/2✓ Branch 1 taken 196529 times.
✗ Branch 2 not taken.
|
196529 | mesh_.facets.attributes().copy_item(f2,f1); |
| 1131 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 196529 times.
|
196529 | if(interpolate_attributes_) { |
| 1132 | ✗ | mesh_.facet_corners.attributes().copy_item(3*f2, 3*f1+2); | |
| 1133 | ✗ | mesh_.facet_corners.attributes().copy_item(3*f2+1,3*f1+1); | |
| 1134 | ✗ | mesh_.facet_corners.attributes().copy_item(3*f2+2,3*f1 ); | |
| 1135 | } | ||
| 1136 | |||
| 1137 | // Sew halfedges (need to be done *after* copy attributes, | ||
| 1138 | // since alpha3 is stored as an .. **attribute** !!! | ||
| 1139 | 196529 | halfedges_.sew3(3*f1, 3*f2+1); | |
| 1140 | 196529 | halfedges_.sew3(3*f1+1,3*f2 ); | |
| 1141 | 196529 | halfedges_.sew3(3*f1+2,3*f2+2); | |
| 1142 | } | ||
| 1143 | } | ||
| 1144 | |||
| 1145 | // Memorize flipped status for facets before classification | ||
| 1146 | // (that changes total number of facets). Used by get_initial_facet() | ||
| 1147 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | f_is_flipped_.bind(mesh_.facets.attributes(), "flipped"); |
| 1148 |
2/2✓ Branch 0 taken 393058 times.
✓ Branch 1 taken 57 times.
|
393115 | for(index_t f: mesh_.facets) { |
| 1149 | 393058 | f_is_flipped_[f] = (f >= mesh_.facets.nb() / 2); | |
| 1150 | } | ||
| 1151 | |||
| 1152 | // Step 2: Clear all facet-facet links | ||
| 1153 |
2/2✓ Branch 0 taken 1179174 times.
✓ Branch 1 taken 57 times.
|
1179231 | for(index_t c: mesh_.facet_corners) { |
| 1154 | mesh_.facet_corners.set_adjacent_facet(c,NO_INDEX); | ||
| 1155 | } | ||
| 1156 | |||
| 1157 | // Step 3: Compute halfedges bundles and radial polylines | ||
| 1158 | { | ||
| 1159 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | Stopwatch W("Bundles",verbose_); |
| 1160 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | radial_bundles_.initialize(); |
| 1161 | 57 | } | |
| 1162 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | radial_polylines_.initialize(); |
| 1163 | |||
| 1164 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | if(skeleton_ != nullptr) { |
| 1165 | ✗ | geo_assert(!dry_run_); | |
| 1166 | ✗ | radial_polylines_.get_skeleton(*skeleton_, skeleton_trim_fins_); | |
| 1167 | } | ||
| 1168 | |||
| 1169 | // Step 4: Connect manifold edges | ||
| 1170 |
2/2✓ Branch 0 taken 551370 times.
✓ Branch 1 taken 57 times.
|
551427 | for(index_t bndl: radial_bundles_) { |
| 1171 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 551370 times.
|
551370 | if(radial_bundles_.nb_halfedges(bndl) == 1) { |
| 1172 | // If there is only one halfedge in the radial bundle, | ||
| 1173 | // this is the original surface border, then we | ||
| 1174 | // create the "hem" by connecting the halfedge on the | ||
| 1175 | // border to its mate in the duplicated facet | ||
| 1176 | // | ||
| 1177 | // Note: if an object has a tiny hole, it has the effect | ||
| 1178 | // of creating a single shell instead of two shells with | ||
| 1179 | // opposite orientation, it may be a problem sometimes... | ||
| 1180 | index_t h = radial_bundles_.halfedge(bndl,0); | ||
| 1181 | halfedges_.sew2(h,halfedges_.alpha3(h)); | ||
| 1182 | corner_is_on_border[h] = true; | ||
| 1183 | corner_is_on_border[halfedges_.alpha3(h)] = true; | ||
| 1184 |
2/2✓ Branch 0 taken 507698 times.
✓ Branch 1 taken 43672 times.
|
551370 | } else if(radial_bundles_.nb_halfedges(bndl) == 2) { |
| 1185 | // If there are two halfedges in the radial bundle, | ||
| 1186 | // then it is a standard manifold edge | ||
| 1187 | index_t h1 = radial_bundles_.halfedge(bndl,0); | ||
| 1188 | index_t h2 = radial_bundles_.halfedge(bndl,1); | ||
| 1189 | halfedges_.sew2(h1,halfedges_.alpha3(h2)); | ||
| 1190 | } | ||
| 1191 | } | ||
| 1192 | |||
| 1193 | // Step 5: get charts | ||
| 1194 | // After this step, charts are interconnected triangles with coherent | ||
| 1195 | // orientation bordered by non-manifold radial edges | ||
| 1196 | index_t nb_charts = 0; | ||
| 1197 | { | ||
| 1198 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | Stopwatch W("Charts",verbose_); |
| 1199 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | nb_charts = get_surface_connected_components(mesh_, "chart"); |
| 1200 | 57 | } | |
| 1201 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 44 times.
|
57 | if(verbose_) { |
| 1202 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
26 | Logger::out("Weiler") |
| 1203 | << "Found " << nb_charts << " charts" << std::endl; | ||
| 1204 | } | ||
| 1205 | |||
| 1206 | // I had before an assertion check that the number of charts is even | ||
| 1207 | // but I removed it: | ||
| 1208 | // the number of charts is not necessarily even, since there | ||
| 1209 | // can be "fins" (see e.g. "saturn" example that creates 5 charts). | ||
| 1210 | |||
| 1211 | // Step 6: Radial sort | ||
| 1212 | { | ||
| 1213 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | Stopwatch("Radial sort", verbose_); |
| 1214 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | radial_polylines_.radial_sort(); |
| 1215 | } | ||
| 1216 | |||
| 1217 | // Step 7: Create alpha2 links | ||
| 1218 |
2/2✓ Branch 0 taken 3820 times.
✓ Branch 1 taken 57 times.
|
3877 | for(index_t P: radial_polylines_) { |
| 1219 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 21836 times.
✓ Branch 2 taken 21836 times.
✓ Branch 3 taken 3820 times.
|
25656 | for(index_t bndl: radial_polylines_.bundles(P)) { |
| 1220 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21836 times.
|
21836 | if(!radial_bundles_.is_sorted(bndl)) { |
| 1221 | ✗ | continue; | |
| 1222 | } | ||
| 1223 | index_t N = radial_bundles_.nb_halfedges(bndl); | ||
| 1224 |
2/2✓ Branch 0 taken 81889 times.
✓ Branch 1 taken 21836 times.
|
103725 | for(index_t i=0; i<N; ++i) { |
| 1225 |
2/2✓ Branch 0 taken 60053 times.
✓ Branch 1 taken 21836 times.
|
81889 | index_t i_next = (i == N-1) ? 0 : i+1; |
| 1226 |
2/2✓ Branch 0 taken 60053 times.
✓ Branch 1 taken 21836 times.
|
81889 | index_t i_prev = (i == 0 ) ? N-1 : i-1; |
| 1227 | index_t h = radial_bundles_.halfedge(bndl,i); | ||
| 1228 | index_t h_next = radial_bundles_.halfedge(bndl,i_next); | ||
| 1229 | index_t h_prev = radial_bundles_.halfedge(bndl,i_prev); | ||
| 1230 | halfedges_.sew2(h,halfedges_.alpha3(h_next)); | ||
| 1231 | halfedges_.sew2(h_prev,halfedges_.alpha3(h)); | ||
| 1232 | } | ||
| 1233 | } | ||
| 1234 | } | ||
| 1235 | |||
| 1236 | // Step 8: Identify regions | ||
| 1237 | { | ||
| 1238 |
2/6✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
57 | Stopwatch W("Regions",verbose_); |
| 1239 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
57 | index_t nb_regions = get_surface_connected_components(mesh_,"chart"); |
| 1240 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 44 times.
|
57 | if(verbose_) { |
| 1241 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
26 | Logger::out("Weiler") |
| 1242 | << "Found " << nb_regions << " regions" << std::endl; | ||
| 1243 | } | ||
| 1244 | 57 | } | |
| 1245 | |||
| 1246 | 57 | } | |
| 1247 | |||
| 1248 | |||
| 1249 | /***********************************************************************/ | ||
| 1250 | |||
| 1251 | 57 | void MeshSurfaceIntersection::RadialBundles::initialize() { | |
| 1252 | // Sorted vector of halfedges | ||
| 1253 | // In this step, we only insert the ones such that v2 < v1 | ||
| 1254 | // For the ones such that v1 > v2, we can deduce the information | ||
| 1255 | // from the first half. | ||
| 1256 | // Note: some of them come from the original halfedges and some of | ||
| 1257 | // them from the ones we just created by duplicating the faces. | ||
| 1258 | 57 | H_.reserve(mesh_.facet_corners.nb()); | |
| 1259 | 57 | H_.resize(0); | |
| 1260 |
2/2✓ Branch 0 taken 1179174 times.
✓ Branch 1 taken 57 times.
|
1179231 | for(index_t h: mesh_.facet_corners) { |
| 1261 |
3/4✓ Branch 0 taken 1179174 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 589587 times.
✓ Branch 3 taken 589587 times.
|
2358348 | if(I_.halfedges_.vertex(h,0) < I_.halfedges_.vertex(h,1)) { |
| 1262 | H_.push_back(h); | ||
| 1263 | } | ||
| 1264 | } | ||
| 1265 | |||
| 1266 | |||
| 1267 | // Step 2: Lexicographic sort of the H array, so that "bundles" will be | ||
| 1268 | // contiguous. By "bundle", I mean the set of halfedges having the same | ||
| 1269 | // extremities in the same order. | ||
| 1270 | 57 | GEO::sort( | |
| 1271 | 57 | H_.begin(), H_.end(), | |
| 1272 | 57 | [&](index_t h1, index_t h2)->bool { | |
| 1273 | geo_debug_assert(h1 < mesh_.facet_corners.nb()); | ||
| 1274 | geo_debug_assert(h2 < mesh_.facet_corners.nb()); | ||
| 1275 |
1/2✓ Branch 0 taken 9669269 times.
✗ Branch 1 not taken.
|
9669269 | index_t v10 = I_.halfedges_.vertex(h1,0); |
| 1276 | index_t v20 = I_.halfedges_.vertex(h2,0); | ||
| 1277 |
2/2✓ Branch 0 taken 4526848 times.
✓ Branch 1 taken 5142421 times.
|
9669269 | if(v10 < v20) { |
| 1278 | return true; | ||
| 1279 | } | ||
| 1280 |
2/2✓ Branch 0 taken 1597479 times.
✓ Branch 1 taken 2929369 times.
|
4526848 | if(v10 > v20) { |
| 1281 | return false; | ||
| 1282 | } | ||
| 1283 | 1597479 | return (I_.halfedges_.vertex(h1,1) < I_.halfedges_.vertex(h2,1)); | |
| 1284 | } | ||
| 1285 | ); | ||
| 1286 | |||
| 1287 | // Step 3: find "bundles" of halfedges in the sorted array. Two halfedges | ||
| 1288 | // are in the same bundle if they connect the same pair of vertices (in | ||
| 1289 | // the same order). | ||
| 1290 | // Halfedges of bundle B are H[b],H[b+1]...H[e-1] | ||
| 1291 | // where b = bndl_start[B] and e = bndl_start[B+1] | ||
| 1292 | // (there are bndl_start.size()-1 bundles) | ||
| 1293 | 57 | bndl_start_.resize(0); | |
| 1294 |
2/2✓ Branch 0 taken 275685 times.
✓ Branch 1 taken 57 times.
|
551484 | for(index_t b=0, e=0; b<H_.size(); b=e) { |
| 1295 | bndl_start_.push_back(b); | ||
| 1296 |
1/2✓ Branch 0 taken 275685 times.
✗ Branch 1 not taken.
|
275685 | index_t v0 = I_.halfedges_.vertex(H_[b],0); |
| 1297 | index_t v1 = I_.halfedges_.vertex(H_[b],1); | ||
| 1298 | 275685 | e = b+1; | |
| 1299 | 275685 | while( | |
| 1300 |
2/2✓ Branch 0 taken 514087 times.
✓ Branch 1 taken 75443 times.
|
589530 | e < H_.size() && |
| 1301 |
5/6✓ Branch 0 taken 589530 times.
✓ Branch 1 taken 57 times.
✓ Branch 2 taken 589530 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 313902 times.
✓ Branch 5 taken 200185 times.
|
1103674 | I_.halfedges_.vertex(H_[e],0) == v0 && |
| 1302 | I_.halfedges_.vertex(H_[e],1) == v1 | ||
| 1303 | ) { | ||
| 1304 | 313902 | ++e; | |
| 1305 | } | ||
| 1306 | } | ||
| 1307 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 57 times.
|
57 | bndl_start_.push_back(H_.size()); |
| 1308 | |||
| 1309 | // Step 4: construct second half of the sorted halfedges array from | ||
| 1310 | // the first half. Remember, at this point, H contains a mixture of | ||
| 1311 | // original halfedges and new ones, but they are all such that v1 < v2. | ||
| 1312 | // We now store all the bundles such that v1 > v2 (by "mirroring" the | ||
| 1313 | // initial ones, that is, traversing alpha3). | ||
| 1314 | { | ||
| 1315 | index_t bndl_start_size = bndl_start_.size(); | ||
| 1316 |
2/2✓ Branch 0 taken 275685 times.
✓ Branch 1 taken 57 times.
|
275742 | for(index_t bndl=0; bndl+1<bndl_start_size; ++bndl) { |
| 1317 | 275685 | index_t b = bndl_start_[bndl]; | |
| 1318 | 275685 | index_t e = bndl_start_[bndl+1]; | |
| 1319 |
2/2✓ Branch 0 taken 589587 times.
✓ Branch 1 taken 275685 times.
|
865272 | for(index_t i=b; i<e; ++i) { |
| 1320 | 589587 | H_.push_back(I_.halfedges_.alpha3(H_[i])); | |
| 1321 | } | ||
| 1322 | 275685 | bndl_start_.push_back(H_.size()); | |
| 1323 | } | ||
| 1324 | } | ||
| 1325 | |||
| 1326 | // Step 5: chain bundles around vertices | ||
| 1327 | 57 | v_first_bndl_.assign(mesh_.vertices.nb(), NO_INDEX); | |
| 1328 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | bndl_next_around_v_.assign(bndl_start_.size()-1, NO_INDEX); |
| 1329 |
2/2✓ Branch 0 taken 551370 times.
✓ Branch 1 taken 57 times.
|
551484 | for(index_t bndl = 0; bndl < nb(); ++bndl) { |
| 1330 | // Skip regular bundles (that are inside charts), | ||
| 1331 | // we only need chaining along non-manifold radial polylines | ||
| 1332 |
2/2✓ Branch 0 taken 507698 times.
✓ Branch 1 taken 43672 times.
|
551370 | if(nb_halfedges(bndl) == 2) { |
| 1333 | 507698 | continue; | |
| 1334 | } | ||
| 1335 | 43672 | index_t v1 = vertex(bndl,0); | |
| 1336 | 43672 | bndl_next_around_v_[bndl] = v_first_bndl_[v1]; | |
| 1337 | 43672 | v_first_bndl_[v1] = bndl; | |
| 1338 | } | ||
| 1339 | |||
| 1340 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 57 times.
|
57 | bndl_is_sorted_.assign(nb(), false); |
| 1341 | |||
| 1342 | if(!facet_chart_.is_bound()) { | ||
| 1343 |
1/2✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
|
114 | facet_chart_.bind(mesh_.facets.attributes(), "chart"); |
| 1344 | } | ||
| 1345 | |||
| 1346 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 44 times.
|
57 | if(I_.verbose_) { |
| 1347 |
1/2✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
|
26 | Logger::out("Weiler") << "Found " << nb() << " bundles" |
| 1348 | << std::endl; | ||
| 1349 | } | ||
| 1350 | 57 | } | |
| 1351 | |||
| 1352 | 21836 | void MeshSurfaceIntersection::RadialBundles::get_sorted_incident_charts( | |
| 1353 | index_t bndl, vector<ChartPos>& chart_pos | ||
| 1354 | ) { | ||
| 1355 | 21836 | chart_pos.resize(0); | |
| 1356 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 81889 times.
✓ Branch 2 taken 81889 times.
✓ Branch 3 taken 21836 times.
|
103725 | for(index_t h: halfedges(bndl)) { |
| 1357 | chart_pos.push_back( | ||
| 1358 | 81889 | std::make_pair( | |
| 1359 | facet_chart_[I_.halfedges_.facet(h)], chart_pos.size() | ||
| 1360 | ) | ||
| 1361 | ); | ||
| 1362 | } | ||
| 1363 | 21836 | std::sort( | |
| 1364 | chart_pos.begin(), chart_pos.end(), | ||
| 1365 | [](const ChartPos& a, const ChartPos& b)->bool { | ||
| 1366 |
4/22✗ 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 22001 times.
✓ Branch 13 taken 38052 times.
✓ Branch 14 taken 19012 times.
✓ Branch 15 taken 38052 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
117117 | return a.first < b.first; |
| 1367 | } | ||
| 1368 | ); | ||
| 1369 | 21836 | } | |
| 1370 | |||
| 1371 | 57 | void MeshSurfaceIntersection::RadialPolylines::initialize() { | |
| 1372 | 57 | B_.resize(0); | |
| 1373 | 57 | polyline_start_.resize(0); | |
| 1374 | 57 | polyline_start_.push_back(0); | |
| 1375 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | vector<bool> bndl_visited(I_.radial_bundles_.nb(),false); |
| 1376 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
✓ Branch 2 taken 551370 times.
✓ Branch 3 taken 57 times.
|
551484 | for(index_t bndl: I_.radial_bundles_) { |
| 1377 | |||
| 1378 | // Bundle already visited, or regular internal edgge | ||
| 1379 |
4/4✓ Branch 0 taken 511518 times.
✓ Branch 1 taken 39852 times.
✓ Branch 2 taken 507698 times.
✓ Branch 3 taken 3820 times.
|
551370 | if(bndl_visited[bndl] || I_.radial_bundles_.nb_halfedges(bndl)==2) { |
| 1380 | 547550 | continue; | |
| 1381 | } | ||
| 1382 | |||
| 1383 | // Find first bundle of the poly line: | ||
| 1384 | // Traverse bundles on border backwards until a | ||
| 1385 | // non-manifold vertex is reached or until we have | ||
| 1386 | // looped to our starting point. | ||
| 1387 | index_t bndl_first = bndl; | ||
| 1388 | for(;;) { | ||
| 1389 | index_t bndl_p = | ||
| 1390 |
1/2✓ Branch 1 taken 14571 times.
✗ Branch 2 not taken.
|
14571 | I_.radial_bundles_.prev_along_polyline(bndl_first); |
| 1391 |
2/2✓ Branch 0 taken 10751 times.
✓ Branch 1 taken 3820 times.
|
14571 | if(bndl_p == NO_INDEX || bndl_p == bndl) { |
| 1392 | break; | ||
| 1393 | } | ||
| 1394 | bndl_first = bndl_p; | ||
| 1395 | } | ||
| 1396 | |||
| 1397 | // Traverse polyline by traversing bundles on border forward | ||
| 1398 | // until a non-manifold vertex is reached or until we have looped | ||
| 1399 | // to our starting point. | ||
| 1400 | 3820 | index_t bndl_cur = bndl_first; | |
| 1401 | for(;;) { | ||
| 1402 | B_.push_back(bndl_cur); | ||
| 1403 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21836 times.
|
21836 | bndl_visited[bndl_cur] = true; |
| 1404 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21836 times.
|
21836 | bndl_visited[I_.radial_bundles_.opposite(bndl_cur)] = true; |
| 1405 |
1/2✓ Branch 1 taken 21836 times.
✗ Branch 2 not taken.
|
21836 | index_t bndl_n=I_.radial_bundles_.next_along_polyline(bndl_cur); |
| 1406 |
2/2✓ Branch 0 taken 18016 times.
✓ Branch 1 taken 3820 times.
|
21836 | if(bndl_n == NO_INDEX || bndl_n == bndl_first) { |
| 1407 | break; | ||
| 1408 | } | ||
| 1409 | 18016 | bndl_cur = bndl_n; | |
| 1410 | 18016 | } | |
| 1411 |
1/2✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
|
3820 | polyline_start_.push_back(B_.size()); |
| 1412 | } | ||
| 1413 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 44 times.
|
57 | if(I_.verbose_) { |
| 1414 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
26 | Logger::out("Weiler") << "Found " <<nb()<< " polylines" |
| 1415 | << std::endl; | ||
| 1416 | } | ||
| 1417 | 57 | } | |
| 1418 | |||
| 1419 | ✗ | void MeshSurfaceIntersection::RadialPolylines::get_skeleton( | |
| 1420 | Mesh& skeleton, bool trim_fins | ||
| 1421 | ) { | ||
| 1422 | ✗ | skeleton.clear(); | |
| 1423 | ✗ | skeleton.vertices.set_dimension(3); | |
| 1424 | Attribute<bool> new_v_selection( | ||
| 1425 | ✗ | skeleton.vertices.attributes(), "selection" | |
| 1426 | ); | ||
| 1427 | Attribute<bool> v_selection( | ||
| 1428 | ✗ | mesh_.vertices.attributes(), "selection" | |
| 1429 | ); | ||
| 1430 | ✗ | vector<index_t> v_id(mesh_.vertices.nb(), NO_INDEX); | |
| 1431 | ✗ | for(index_t bndl: B_) { | |
| 1432 | ✗ | for(index_t lv=0; lv<2; ++lv) { | |
| 1433 | ✗ | index_t v = I_.radial_bundles_.vertex(bndl,lv); | |
| 1434 | ✗ | if(v_id[v] == NO_INDEX) { | |
| 1435 | ✗ | v_id[v] = skeleton.vertices.create_vertex( | |
| 1436 | ✗ | mesh_.vertices.point(v) | |
| 1437 | ); | ||
| 1438 | ✗ | if(I_.radial_bundles_.nb_bundles_around_vertex(v) != 2) { | |
| 1439 | v_selection[v] = true; | ||
| 1440 | ✗ | new_v_selection[v_id[v]] = true; | |
| 1441 | } | ||
| 1442 | } | ||
| 1443 | } | ||
| 1444 | } | ||
| 1445 | ✗ | for(index_t polyline: *this) { | |
| 1446 | ✗ | for(index_t bndl: bundles(polyline)) { | |
| 1447 | ✗ | if(trim_fins && I_.radial_bundles_.nb_halfedges(bndl) < 3) { | |
| 1448 | ✗ | continue; | |
| 1449 | } | ||
| 1450 | ✗ | index_t v1 = I_.radial_bundles_.vertex(bndl,0); | |
| 1451 | ✗ | index_t v2 = I_.radial_bundles_.vertex(bndl,1); | |
| 1452 | geo_debug_assert(v_id[v1] != NO_INDEX); | ||
| 1453 | geo_debug_assert(v_id[v2] != NO_INDEX); | ||
| 1454 | ✗ | skeleton.edges.create_edge(v_id[v1], v_id[v2]); | |
| 1455 | } | ||
| 1456 | } | ||
| 1457 | ✗ | skeleton.vertices.remove_isolated(); | |
| 1458 | ✗ | } | |
| 1459 | |||
| 1460 | 57 | void MeshSurfaceIntersection::RadialPolylines::radial_sort() { | |
| 1461 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 44 times.
|
57 | if(I_.verbose_) { |
| 1462 |
1/2✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
|
26 | Logger::out("Radial sort") |
| 1463 | << "Nb radial polylines:" << nb() << std::endl; | ||
| 1464 | } | ||
| 1465 |
1/2✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
|
57 | Stopwatch W("Radial sort",I_.verbose_); |
| 1466 | |||
| 1467 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | index_t nb_sorted = 0; |
| 1468 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | index_t nb_to_sort = nb(); |
| 1469 | |||
| 1470 | // For each polyline in parallel | ||
| 1471 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | parallel_for_slice( |
| 1472 | 0, nb(), | ||
| 1473 | 246 | [&](index_t b, index_t e) { | |
| 1474 | 189 | index_t tid = Thread::current_id(); | |
| 1475 | 189 | RadialSort RS(I_); | |
| 1476 | |||
| 1477 | vector<RadialBundles::ChartPos> ref_chart_to_radial_id; | ||
| 1478 | vector<RadialBundles::ChartPos> cur_chart_to_radial_id; | ||
| 1479 | vector<index_t> bndl_h; | ||
| 1480 | |||
| 1481 | // P: current polyline | ||
| 1482 |
2/2✓ Branch 0 taken 3820 times.
✓ Branch 1 taken 189 times.
|
4009 | for(index_t P = b; P < e; ++P) { |
| 1483 | |||
| 1484 | index_t bndl_ref = NO_INDEX; // reference bundle | ||
| 1485 | index_t N = NO_INDEX; // nb halfedges in ref bundle | ||
| 1486 | |||
| 1487 |
1/2✓ Branch 0 taken 3820 times.
✗ Branch 1 not taken.
|
3820 | for(index_t bndl: bundles(P)) { |
| 1488 |
1/2✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
|
3820 | bool OK = I_.radial_bundles_.radial_sort(bndl, RS); |
| 1489 |
1/2✓ Branch 0 taken 3820 times.
✗ Branch 1 not taken.
|
3820 | if(OK) { |
| 1490 |
1/2✓ Branch 1 taken 3820 times.
✗ Branch 2 not taken.
|
3820 | I_.radial_bundles_.get_sorted_incident_charts( |
| 1491 | bndl, ref_chart_to_radial_id | ||
| 1492 | ); | ||
| 1493 | for( | ||
| 1494 | index_t i=0; | ||
| 1495 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 13909 times.
✓ Branch 2 taken 10089 times.
✓ Branch 3 taken 3820 times.
|
13909 | i+1<ref_chart_to_radial_id.size(); ++i |
| 1496 | ) { | ||
| 1497 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10089 times.
|
10089 | OK = OK && ( |
| 1498 | 10089 | ref_chart_to_radial_id[i].first != | |
| 1499 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10089 times.
|
10089 | ref_chart_to_radial_id[i+1].first |
| 1500 | ); | ||
| 1501 | } | ||
| 1502 | } | ||
| 1503 |
1/2✓ Branch 0 taken 3820 times.
✗ Branch 1 not taken.
|
3820 | if(OK) { |
| 1504 | bndl_ref = bndl; | ||
| 1505 | 3820 | N = I_.radial_bundles_.nb_halfedges(bndl); | |
| 1506 | 3820 | break; | |
| 1507 | } | ||
| 1508 | } | ||
| 1509 | |||
| 1510 | // Copy order to all other bundles in polyline | ||
| 1511 |
2/2✓ Branch 0 taken 21836 times.
✓ Branch 1 taken 3820 times.
|
25656 | for(index_t bndl: bundles(P)) { |
| 1512 |
2/2✓ Branch 0 taken 3820 times.
✓ Branch 1 taken 18016 times.
|
21836 | if(bndl == bndl_ref) { |
| 1513 | 3820 | continue; | |
| 1514 | } | ||
| 1515 | |||
| 1516 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18016 times.
|
18016 | if(I_.radial_bundles_.nb_halfedges(bndl)==1) { |
| 1517 | ✗ | continue; | |
| 1518 | } | ||
| 1519 | |||
| 1520 | // Necessary condition: reference bndl should be OK | ||
| 1521 | // and current bndl should have same nbr of halfedges | ||
| 1522 |
2/4✓ Branch 0 taken 18016 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18016 times.
✗ Branch 3 not taken.
|
18016 | bool OK = (bndl_ref != NO_INDEX) && ( |
| 1523 | I_.radial_bundles_.nb_halfedges(bndl)==N | ||
| 1524 | ); | ||
| 1525 | |||
| 1526 | // Now check that charts around current bndl are the | ||
| 1527 | // same as charts around reference bndl. It can happen | ||
| 1528 | // that they differ even if N matches (example21.csg) | ||
| 1529 | if(OK) { | ||
| 1530 |
1/2✓ Branch 1 taken 18016 times.
✗ Branch 2 not taken.
|
18016 | I_.radial_bundles_.get_sorted_incident_charts( |
| 1531 | bndl, cur_chart_to_radial_id | ||
| 1532 | ); | ||
| 1533 |
2/2✓ Branch 0 taken 67980 times.
✓ Branch 1 taken 18016 times.
|
85996 | for(index_t i=0; i<N; ++i) { |
| 1534 |
1/2✓ Branch 0 taken 67980 times.
✗ Branch 1 not taken.
|
67980 | OK = OK && ( |
| 1535 | 67980 | cur_chart_to_radial_id[i].first == | |
| 1536 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 67980 times.
|
67980 | ref_chart_to_radial_id[i].first |
| 1537 | ); | ||
| 1538 | } | ||
| 1539 | } | ||
| 1540 | |||
| 1541 |
1/2✓ Branch 0 taken 18016 times.
✗ Branch 1 not taken.
|
18016 | if(OK) { |
| 1542 | // Here ref_bnfl has a valid order, and bndl has | ||
| 1543 | // the same surrounding charts as ref_bndl, | ||
| 1544 | // so we can reuse the radial order computed | ||
| 1545 | // in ref_bndl | ||
| 1546 | bndl_h.assign(N, NO_INDEX); | ||
| 1547 |
2/2✓ Branch 0 taken 67980 times.
✓ Branch 1 taken 18016 times.
|
85996 | for(index_t i=0; i<N; ++i) { |
| 1548 | // This halfedge ... | ||
| 1549 | 67980 | index_t h = I_.radial_bundles_.halfedge( | |
| 1550 | 67980 | bndl, cur_chart_to_radial_id[i].second | |
| 1551 | ); | ||
| 1552 | // ... goes here. | ||
| 1553 | 67980 | bndl_h[ref_chart_to_radial_id[i].second] = h; | |
| 1554 | } | ||
| 1555 | 18016 | I_.radial_bundles_.set_sorted_halfedges(bndl,bndl_h); | |
| 1556 | } else { | ||
| 1557 | // Else compute the radial sort geometrically | ||
| 1558 | ✗ | OK = I_.radial_bundles_.radial_sort(bndl, RS); | |
| 1559 | // May return !OK (if not using geogram_plus) when it | ||
| 1560 | // cannot sort (example_022.csg and example_024.csg) | ||
| 1561 | ✗ | if(!OK) { | |
| 1562 | std::cerr << std::endl | ||
| 1563 | << "FATAL ERROR: " | ||
| 1564 | << "Did not manage to sort a bundle" | ||
| 1565 | << " in Polyline of length " | ||
| 1566 | ✗ | << nb_bundles(P) | |
| 1567 | << std::endl; | ||
| 1568 | std::cerr | ||
| 1569 | << "(if you reached this point, you may" | ||
| 1570 | << " need geogramplus, contact TESSAEL)" | ||
| 1571 | << std::endl; | ||
| 1572 | ✗ | geo_assert_not_reached; | |
| 1573 | } | ||
| 1574 | } | ||
| 1575 | } | ||
| 1576 |
2/2✓ Branch 0 taken 1500 times.
✓ Branch 1 taken 2320 times.
|
3820 | if(I_.fine_verbose_) { |
| 1577 | 1500 | ++nb_sorted; | |
| 1578 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1488 times.
|
1500 | if(!(nb_sorted%100)) { |
| 1579 |
2/4✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
|
24 | Logger::out("Radial sort") |
| 1580 |
1/4✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12 | << String::format( |
| 1581 | "[%2d] %6d/%6d", | ||
| 1582 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | int(tid),int(nb_sorted),int(nb_to_sort) |
| 1583 | ) | ||
| 1584 | << std::endl; | ||
| 1585 | } | ||
| 1586 | } | ||
| 1587 | } | ||
| 1588 |
2/2✓ Branch 0 taken 49 times.
✓ Branch 1 taken 140 times.
|
189 | if(I_.fine_verbose_) { |
| 1589 |
2/6✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
98 | Logger::out("Radial sort") |
| 1590 |
2/6✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 49 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
98 | << String::format("[%2d] done",int(tid)) |
| 1591 | << std::endl; | ||
| 1592 | } | ||
| 1593 | 189 | } | |
| 1594 | ); | ||
| 1595 | 57 | } | |
| 1596 | |||
| 1597 | /***********************************************************************/ | ||
| 1598 | } | ||
| 1599 | |||
| 1600 | |||
| 1601 | namespace { | ||
| 1602 | using namespace GEO; | ||
| 1603 | |||
| 1604 | /** | ||
| 1605 | * \brief Gets the position of the leftmost | ||
| 1606 | * bit set in a 32 bits integer | ||
| 1607 | * \param[in] x the integer | ||
| 1608 | * \return the position of the leftmost bit | ||
| 1609 | * set, or NO_INDEX if the specified integer | ||
| 1610 | * is zero. | ||
| 1611 | */ | ||
| 1612 | inline index_t leftmost_bit_set(index_t x) { | ||
| 1613 | index_t result = NO_INDEX; | ||
| 1614 |
2/2✓ Branch 0 taken 960 times.
✓ Branch 1 taken 30 times.
|
990 | for(index_t i=0; i<32; ++i) { |
| 1615 |
2/2✓ Branch 0 taken 67 times.
✓ Branch 1 taken 893 times.
|
960 | if((x&1) != 0) { |
| 1616 | result = i; | ||
| 1617 | } | ||
| 1618 | 960 | x = x >> 1; | |
| 1619 | } | ||
| 1620 | return result; | ||
| 1621 | } | ||
| 1622 | } | ||
| 1623 | |||
| 1624 | /***************************************************/ | ||
| 1625 | |||
| 1626 | namespace GEO { | ||
| 1627 | |||
| 1628 | 30 | void MeshSurfaceIntersection::classify(const std::string& expr) { | |
| 1629 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
|
30 | if(mesh_.facets.nb() == 0) { |
| 1630 | ✗ | return; | |
| 1631 | } | ||
| 1632 |
1/2✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
|
30 | Stopwatch W("Classify", verbose_); |
| 1633 | |||
| 1634 | // Takes as input a Weiler model, with duplicated interfaces, | ||
| 1635 | // operand bit (that indices for each triangle the set of operands | ||
| 1636 | // it corresponds to), volumetric alpha3 links and correct facet | ||
| 1637 | // adjacency links. It computes the operand_inclusion_bits attribute, | ||
| 1638 | // that indicates for each triangle the set of operands that contains | ||
| 1639 | // it, then evalues the boolean expression \p expr on all facets, | ||
| 1640 | // and keeps the facets on the boundary of the region where \p expr | ||
| 1641 | // evaluates to true. | ||
| 1642 | |||
| 1643 | // Chart attribute corresponds to volumetric regions | ||
| 1644 |
2/4✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
|
30 | Attribute<index_t> chart(mesh_.facets.attributes(), "chart"); |
| 1645 | |||
| 1646 | // For each facet, bit n set if facet belongs to the boundary | ||
| 1647 | // of operand n. There can be several bit sets if two operands | ||
| 1648 | // are tangent (and share facets). | ||
| 1649 | // For each facet pair (f, g=halfedges_.facet_alpha3(f)), | ||
| 1650 | // we have operand_bit[f] = operand_bit[g] | ||
| 1651 |
2/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
30 | Attribute<index_t> operand_bit(mesh_.facets.attributes(), "operand_bit"); |
| 1652 | |||
| 1653 | // For each facet, bit n set if facet is inside operand n. | ||
| 1654 | // For each facet pair (f, g=halfedges_.facet_alpha3(f)), | ||
| 1655 | // we have nth-bit(operand_bit[f]) != nth-bit(operand_bit[g]), | ||
| 1656 | // that is, among a facet pair (f,g) on the boundary of operand n, | ||
| 1657 | // one of f,g is considered to be "outside" the operand, and the | ||
| 1658 | // other one is considered to be "inside". | ||
| 1659 | // The one that is considered to be "outside" is the one that belongs | ||
| 1660 | // to the chart that encloses the largest volume. | ||
| 1661 | Attribute<index_t> operand_inclusion_bits( | ||
| 1662 |
2/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
30 | mesh_.facets.attributes(), "operand_inclusion_bits" |
| 1663 | ); | ||
| 1664 | |||
| 1665 | // Get nb charts and nb operands | ||
| 1666 | 30 | index_t nb_charts = 0; | |
| 1667 | index_t nb_operands = 0; | ||
| 1668 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1669 | 175992 | nb_charts = std::max(nb_charts, chart[f]+1); | |
| 1670 | 175992 | nb_operands = nb_operands | operand_bit[f]; | |
| 1671 | } | ||
| 1672 | |||
| 1673 |
1/2✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
|
30 | if(nb_operands != 0) { |
| 1674 | 30 | nb_operands = leftmost_bit_set(nb_operands) + 1; | |
| 1675 | } | ||
| 1676 | |||
| 1677 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
|
30 | if(verbose_) { |
| 1678 | ✗ | Logger::out("Weiler") << "nb operands=" << nb_operands << std::endl; | |
| 1679 | } | ||
| 1680 | |||
| 1681 | // Get connected components, obtained by traversing all facet | ||
| 1682 | // adjacency links and volumetric alpha3 links | ||
| 1683 | |||
| 1684 | index_t nb_components = 0; | ||
| 1685 | |||
| 1686 | Attribute<index_t> facet_component( | ||
| 1687 |
2/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
30 | mesh_.facets.attributes(), "component" |
| 1688 | ); | ||
| 1689 | |||
| 1690 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1691 | 175992 | facet_component[f] = NO_INDEX; | |
| 1692 | } | ||
| 1693 | // TODO: Switch to vector when I will not need to visualize anymore | ||
| 1694 | // vector<index_t> facet_component(mesh_.facets.nb(), NO_INDEX); | ||
| 1695 | |||
| 1696 | // one vertex per component | ||
| 1697 | vector<index_t> component_vertex; | ||
| 1698 | |||
| 1699 | // ith bit is 1 if component is inside ith operand | ||
| 1700 | vector<index_t> component_inclusion_bits; | ||
| 1701 | |||
| 1702 | // get the connected components obtained by traversing both alpha2 | ||
| 1703 | // and alpha3 links | ||
| 1704 | { | ||
| 1705 |
4/4✓ Branch 0 taken 36 times.
✓ Branch 1 taken 175956 times.
✓ Branch 2 taken 175992 times.
✓ Branch 3 taken 30 times.
|
176022 | for(index_t f:mesh_.facets) { |
| 1706 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 175956 times.
|
175992 | if(facet_component[f] == NO_INDEX) { |
| 1707 | |||
| 1708 |
2/4✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
|
72 | component_vertex.push_back(mesh_.facets.vertex(f,0)); |
| 1709 |
2/4✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
|
36 | component_inclusion_bits.push_back(0); |
| 1710 | |||
| 1711 | std::stack<index_t> S; | ||
| 1712 |
1/2✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
|
36 | facet_component[f] = nb_components; |
| 1713 | S.push(f); | ||
| 1714 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 36 times.
|
176028 | while(!S.empty()) { |
| 1715 |
2/2✓ Branch 0 taken 173284 times.
✓ Branch 1 taken 2708 times.
|
175992 | index_t f1 = S.top(); |
| 1716 | S.pop(); | ||
| 1717 | { | ||
| 1718 |
2/2✓ Branch 0 taken 55809 times.
✓ Branch 1 taken 120183 times.
|
175992 | index_t f2 = halfedges_.facet_alpha3(f1); |
| 1719 | geo_debug_assert(f2 != NO_INDEX); | ||
| 1720 |
2/2✓ Branch 0 taken 55809 times.
✓ Branch 1 taken 120183 times.
|
175992 | if(facet_component[f2] == NO_INDEX) { |
| 1721 |
2/2✓ Branch 0 taken 54962 times.
✓ Branch 1 taken 847 times.
|
55809 | facet_component[f2]=facet_component[f1]; |
| 1722 | S.push(f2); | ||
| 1723 | } | ||
| 1724 | } | ||
| 1725 | |||
| 1726 |
2/2✓ Branch 0 taken 527976 times.
✓ Branch 1 taken 175992 times.
|
703968 | for(index_t le1=0; le1<3; ++le1) { |
| 1727 |
1/2✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
|
527976 | index_t f2 = mesh_.facets.adjacent(f1,le1); |
| 1728 | if( | ||
| 1729 |
1/2✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
|
527976 | f2 != NO_INDEX && |
| 1730 |
2/2✓ Branch 0 taken 120147 times.
✓ Branch 1 taken 407829 times.
|
527976 | facet_component[f2] == NO_INDEX |
| 1731 | ) { | ||
| 1732 | #ifdef GEO_DEBUG | ||
| 1733 | index_t le2 = mesh_.facets.find_adjacent(f2,f1); | ||
| 1734 | geo_debug_assert( | ||
| 1735 | mesh_.facets.vertex(f1,le1) != | ||
| 1736 | mesh_.facets.vertex(f2,le2) | ||
| 1737 | ); | ||
| 1738 | #endif | ||
| 1739 |
2/2✓ Branch 0 taken 118286 times.
✓ Branch 1 taken 1861 times.
|
120147 | facet_component[f2] = facet_component[f1]; |
| 1740 | S.push(f2); | ||
| 1741 | } | ||
| 1742 | } | ||
| 1743 | } | ||
| 1744 | 36 | ++nb_components; | |
| 1745 | } | ||
| 1746 | } | ||
| 1747 | |||
| 1748 | // Prefer original vertices for starting raytracing | ||
| 1749 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1750 |
1/2✓ Branch 0 taken 175992 times.
✗ Branch 1 not taken.
|
175992 | index_t component = facet_component[f]; |
| 1751 |
1/2✓ Branch 0 taken 175992 times.
✗ Branch 1 not taken.
|
175992 | index_t v = component_vertex[component]; |
| 1752 | // If component's vertex is already an original vertex | ||
| 1753 | // we are done. | ||
| 1754 |
1/2✓ Branch 0 taken 175992 times.
✗ Branch 1 not taken.
|
175992 | if(is_original_vertex(v)) { |
| 1755 | 175992 | continue; | |
| 1756 | } | ||
| 1757 | // See if we can find an original vertex in facet's vertices, | ||
| 1758 | // if yes, replace component's vertex. | ||
| 1759 | ✗ | for(index_t w: mesh_.facets.vertices(f)) { | |
| 1760 | ✗ | if(is_original_vertex(w)) { | |
| 1761 | ✗ | component_vertex[component] = w; | |
| 1762 | ✗ | break; | |
| 1763 | } | ||
| 1764 | } | ||
| 1765 | } | ||
| 1766 | } | ||
| 1767 | |||
| 1768 | // Compute the volume enclosed by each chart | ||
| 1769 | |||
| 1770 |
1/4✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
30 | vector<double> chart_volume(nb_charts,0.0); |
| 1771 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1772 | 175992 | vec3 p1 = mesh_.facets.point(f,0); | |
| 1773 | 175992 | vec3 p2 = mesh_.facets.point(f,1); | |
| 1774 | 175992 | vec3 p3 = mesh_.facets.point(f,2); | |
| 1775 | 175992 | chart_volume[chart[f]] += dot(p1,cross(p2,p3)) / 6.0; | |
| 1776 | } | ||
| 1777 | |||
| 1778 | // For each component, find the chart that encloses the largest | ||
| 1779 | // volume (it is the external boundary of the component) | ||
| 1780 | |||
| 1781 |
2/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
30 | vector<double> max_chart_volume_in_component(nb_components, 0.0); |
| 1782 | vector<index_t> chart_with_max_volume_in_component( | ||
| 1783 | nb_components, NO_INDEX | ||
| 1784 | ); | ||
| 1785 | |||
| 1786 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1787 |
2/2✓ Branch 0 taken 76264 times.
✓ Branch 1 taken 99728 times.
|
175992 | double V = chart_volume[chart[f]]; |
| 1788 | 175992 | if( ::fabs(V) >= | |
| 1789 |
2/2✓ Branch 0 taken 76264 times.
✓ Branch 1 taken 99728 times.
|
175992 | ::fabs(max_chart_volume_in_component[facet_component[f]]) |
| 1790 | ) { | ||
| 1791 | 76264 | max_chart_volume_in_component[facet_component[f]] = V; | |
| 1792 | 76264 | chart_with_max_volume_in_component[facet_component[f]] = | |
| 1793 | chart[f]; | ||
| 1794 | } | ||
| 1795 | } | ||
| 1796 | |||
| 1797 | |||
| 1798 | // If there is more than one component, one needs to check | ||
| 1799 | // whether some components "float" inside other ones. | ||
| 1800 | // To do that, we determine the component inclusion bits | ||
| 1801 | // by launching a ray from a vertex of the component, and | ||
| 1802 | // checking parity of the number of intersections for each operand. | ||
| 1803 | |||
| 1804 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 25 times.
|
30 | if(nb_components > 1) { |
| 1805 | |||
| 1806 | // Copy facet_component[] from mesh_ to mesh_copy_, | ||
| 1807 | // following original_facet_id_ links. | ||
| 1808 | // Note: when some input facet are overlapping and co-planar, | ||
| 1809 | // some original facets may become orphan (facet_component_copy_ | ||
| 1810 | // remains NO_INDEX for those facets). | ||
| 1811 | // They are ignored in raytracing. This may be a problem in very | ||
| 1812 | // specific configurations (to be investigated). | ||
| 1813 | // TODO: find a way of keeping original facet <-> intersection facet | ||
| 1814 | // relations even when facets are merged. | ||
| 1815 | { | ||
| 1816 | Attribute<index_t> facet_component_copy( | ||
| 1817 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
10 | mesh_copy_.facets.attributes(), "component" |
| 1818 | ); | ||
| 1819 |
2/2✓ Branch 0 taken 830 times.
✓ Branch 1 taken 5 times.
|
835 | for(index_t f: mesh_copy_.facets) { |
| 1820 | 830 | facet_component_copy[f] = NO_INDEX; | |
| 1821 | } | ||
| 1822 |
2/2✓ Branch 0 taken 3836 times.
✓ Branch 1 taken 5 times.
|
3841 | for(index_t f: mesh_.facets) { |
| 1823 | 3836 | index_t original_f = original_facet_id_[f]; | |
| 1824 | 3836 | index_t component = facet_component[f]; | |
| 1825 | 3836 | facet_component_copy[original_f] = component; | |
| 1826 | } | ||
| 1827 | } | ||
| 1828 | |||
| 1829 | |||
| 1830 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(verbose_) { |
| 1831 | ✗ | Logger::out("Weiler") << "Classifying " << nb_components | |
| 1832 | << " components using ray tracing" | ||
| 1833 | << std::endl; | ||
| 1834 | } | ||
| 1835 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | parallel_for( |
| 1836 | 16 | 0, nb_components, [&](index_t component) { | |
| 1837 | 11 | component_inclusion_bits[component] = | |
| 1838 | 11 | compute_component_inclusion_bits( | |
| 1839 | 11 | component, component_vertex[component] | |
| 1840 | ); | ||
| 1841 | 11 | } | |
| 1842 | ); | ||
| 1843 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(verbose_) { |
| 1844 | ✗ | Logger::out("Weiler") << "Done." << std::endl; | |
| 1845 | } | ||
| 1846 | } | ||
| 1847 | |||
| 1848 | // Compute operand inclusion bits for each facet, | ||
| 1849 | // by propagating component's inclusion bits | ||
| 1850 | // from component's external shell | ||
| 1851 | |||
| 1852 | { | ||
| 1853 |
2/4✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
|
30 | vector<index_t> visited(mesh_.facets.nb(), false); |
| 1854 | std::stack<index_t> S; | ||
| 1855 | |||
| 1856 |
4/4✓ Branch 0 taken 49804 times.
✓ Branch 1 taken 126188 times.
✓ Branch 2 taken 175992 times.
✓ Branch 3 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1857 | 175992 | if(chart[f] == | |
| 1858 |
2/2✓ Branch 0 taken 49804 times.
✓ Branch 1 taken 126188 times.
|
175992 | chart_with_max_volume_in_component[facet_component[f]] |
| 1859 | ) { | ||
| 1860 | 49804 | visited[f] = 1; | |
| 1861 | 49804 | operand_inclusion_bits[f] = | |
| 1862 |
2/2✓ Branch 0 taken 49038 times.
✓ Branch 1 taken 766 times.
|
49804 | component_inclusion_bits[facet_component[f]]; |
| 1863 | S.push(f); | ||
| 1864 | } | ||
| 1865 | } | ||
| 1866 | |||
| 1867 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | while(!S.empty()) { |
| 1868 |
2/2✓ Branch 0 taken 173233 times.
✓ Branch 1 taken 2759 times.
|
175992 | index_t f1 = S.top(); |
| 1869 | S.pop(); | ||
| 1870 | { | ||
| 1871 |
2/2✓ Branch 0 taken 23877 times.
✓ Branch 1 taken 152115 times.
|
175992 | index_t f2 = halfedges_.facet_alpha3(f1); |
| 1872 |
2/2✓ Branch 0 taken 23877 times.
✓ Branch 1 taken 152115 times.
|
175992 | if(f2 != NO_INDEX && !visited[f2]) { |
| 1873 |
2/2✓ Branch 0 taken 23501 times.
✓ Branch 1 taken 376 times.
|
23877 | visited[f2] = true; |
| 1874 | S.push(f2); | ||
| 1875 | 23877 | operand_inclusion_bits[f2] = | |
| 1876 | 23877 | operand_inclusion_bits[f1] ^ operand_bit[f1]; | |
| 1877 | } | ||
| 1878 | } | ||
| 1879 |
2/2✓ Branch 0 taken 527976 times.
✓ Branch 1 taken 175992 times.
|
703968 | for(index_t le=0; le<3; ++le) { |
| 1880 |
1/2✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
|
527976 | index_t f2 = mesh_.facets.adjacent(f1,le); |
| 1881 |
3/4✓ Branch 0 taken 527976 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 102311 times.
✓ Branch 3 taken 425665 times.
|
527976 | if(f2 != NO_INDEX && !visited[f2]) { |
| 1882 |
2/2✓ Branch 0 taken 100694 times.
✓ Branch 1 taken 1617 times.
|
102311 | visited[f2] = true; |
| 1883 | S.push(f2); | ||
| 1884 | 102311 | operand_inclusion_bits[f2] = operand_inclusion_bits[f1]; | |
| 1885 | } | ||
| 1886 | } | ||
| 1887 | } | ||
| 1888 | } | ||
| 1889 | |||
| 1890 | // Classify facets based on ther operand inclusion bits and on the | ||
| 1891 | // boolean expression | ||
| 1892 | |||
| 1893 |
1/4✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
30 | vector<index_t> classify_facet(mesh_.facets.nb(), 0); |
| 1894 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 21 times.
|
30 | if(expr == "intersection") { |
| 1895 | // If operation is an intersection, return the neighbors of | ||
| 1896 | // the facets that have all their operand inclusion bit sets. | ||
| 1897 | 9 | index_t all_bits_set = (1u << nb_operands)-1u; | |
| 1898 |
2/2✓ Branch 0 taken 81288 times.
✓ Branch 1 taken 9 times.
|
81297 | for(index_t f: mesh_.facets) { |
| 1899 | bool flipped = | ||
| 1900 |
1/2✓ Branch 0 taken 81288 times.
✗ Branch 1 not taken.
|
81288 | (max_chart_volume_in_component[facet_component[f]] < 0.0); |
| 1901 |
1/2✓ Branch 0 taken 81288 times.
✗ Branch 1 not taken.
|
81288 | if(flipped) { |
| 1902 | 81288 | classify_facet[f] = | |
| 1903 | 81288 | (operand_inclusion_bits[f] == all_bits_set); | |
| 1904 | } else { | ||
| 1905 | ✗ | classify_facet[f] = | |
| 1906 | ( | ||
| 1907 | ✗ | operand_inclusion_bits[halfedges_.facet_alpha3(f)] == | |
| 1908 | all_bits_set | ||
| 1909 | ); | ||
| 1910 | } | ||
| 1911 | } | ||
| 1912 | } else { | ||
| 1913 | |||
| 1914 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
|
42 | BooleanExpression E(expr == "union" ? "*" : expr); |
| 1915 | |||
| 1916 | // For a general operation, return the facets f for which the | ||
| 1917 | // expression evaluates to false on f and to true on the neighbors | ||
| 1918 | // of f. Rember: what we want to compute is the *boundary* of the | ||
| 1919 | // region defined by the boolean expression, that is, the facets | ||
| 1920 | // for which the result of the boolean expression changes when they | ||
| 1921 | // are traversed by alpha3. | ||
| 1922 | |||
| 1923 | 94704 | auto classify = [&](index_t f)->bool { | |
| 1924 | bool flipped = | ||
| 1925 |
1/2✓ Branch 0 taken 94704 times.
✗ Branch 1 not taken.
|
94704 | (max_chart_volume_in_component[facet_component[f]]<0.0); |
| 1926 | 94704 | index_t f_in_sets = operand_inclusion_bits[f]; | |
| 1927 | index_t g_in_sets = operand_inclusion_bits[ | ||
| 1928 |
1/2✓ Branch 0 taken 94704 times.
✗ Branch 1 not taken.
|
94704 | halfedges_.facet_alpha3(f) |
| 1929 | 94704 | ]; | |
| 1930 | return flipped | ||
| 1931 |
4/6✓ Branch 0 taken 94704 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 23506 times.
✓ Branch 4 taken 71198 times.
✓ Branch 6 taken 23506 times.
✗ Branch 7 not taken.
|
94704 | ? (E(f_in_sets) && !E(g_in_sets)) |
| 1932 | ✗ | : (E(g_in_sets) && !E(f_in_sets)); | |
| 1933 | 21 | }; | |
| 1934 | |||
| 1935 | try { | ||
| 1936 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 13 times.
|
21 | if(mesh_.facets.nb() > 1024) { |
| 1937 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | parallel_for( |
| 1938 | 8 | 0, mesh_.facets.nb(), [&](index_t f) { | |
| 1939 | 88700 | classify_facet[f] = classify(f); | |
| 1940 | } | ||
| 1941 | ); | ||
| 1942 | } else { | ||
| 1943 |
2/2✓ Branch 0 taken 6004 times.
✓ Branch 1 taken 13 times.
|
6017 | for(index_t f: mesh_.facets) { |
| 1944 |
1/2✓ Branch 1 taken 6004 times.
✗ Branch 2 not taken.
|
6004 | classify_facet[f] = classify(f); |
| 1945 | } | ||
| 1946 | } | ||
| 1947 | ✗ | } catch(...) { | |
| 1948 | ✗ | } | |
| 1949 | } | ||
| 1950 | |||
| 1951 | // invert classification (we mark the facets that | ||
| 1952 | // we want to delete) | ||
| 1953 |
2/2✓ Branch 0 taken 175992 times.
✓ Branch 1 taken 30 times.
|
176022 | for(index_t f: mesh_.facets) { |
| 1954 | 175992 | classify_facet[f] = 1u - classify_facet[f]; | |
| 1955 | } | ||
| 1956 | |||
| 1957 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
30 | mesh_.facets.delete_elements(classify_facet); |
| 1958 | //remove_fins(); | ||
| 1959 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
30 | mesh_.facets.connect(); |
| 1960 | |||
| 1961 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
|
30 | if(verbose_) { |
| 1962 | ✗ | Logger::out("Weiler") << "Facets classified" << std::endl; | |
| 1963 | } | ||
| 1964 | 30 | } | |
| 1965 | |||
| 1966 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | index_t MeshSurfaceIntersection::compute_component_inclusion_bits( |
| 1967 | index_t component, index_t v | ||
| 1968 | ) { | ||
| 1969 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | if(!is_original_vertex(v)) { |
| 1970 | ✗ | return compute_component_inclusion_bits_exact(component, v); | |
| 1971 | } | ||
| 1972 | |||
| 1973 | Attribute<index_t> operand_bit_copy( | ||
| 1974 | 11 | mesh_copy_.facets.attributes(), "operand_bit" | |
| 1975 | ); | ||
| 1976 | |||
| 1977 | Attribute<index_t> facet_component_copy( | ||
| 1978 |
1/4✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
11 | mesh_copy_.facets.attributes(), "component" |
| 1979 | ); | ||
| 1980 | |||
| 1981 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | if(verbose_) { |
| 1982 | ✗ | Logger::out("Weiler") << " component" << component << std::endl; | |
| 1983 | } | ||
| 1984 | |||
| 1985 | index_t component_inclusion_bits = 0; | ||
| 1986 | |||
| 1987 | 11 | vec3 q1 = mesh_.vertices.point(v); | |
| 1988 | 22 | vec3 q2{q1.x, q1.y, q1.z + 1e6}; | |
| 1989 | |||
| 1990 |
2/2✓ Branch 0 taken 1882 times.
✓ Branch 1 taken 11 times.
|
1893 | for(index_t f: mesh_copy_.facets) { |
| 1991 | |||
| 1992 | // TODO: understand when it can happen (probably when merging | ||
| 1993 | // overlapping facets, one of the merged facets is no longer seen | ||
| 1994 | // by the mesh_ -> mesh_copy_ original_facet relation). | ||
| 1995 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1882 times.
|
1882 | if(facet_component_copy[f] == NO_INDEX) { |
| 1996 | ✗ | continue; | |
| 1997 | } | ||
| 1998 | |||
| 1999 |
2/2✓ Branch 0 taken 830 times.
✓ Branch 1 taken 1052 times.
|
1882 | if(facet_component_copy[f] == component) { |
| 2000 | 830 | continue; | |
| 2001 | } | ||
| 2002 | |||
| 2003 |
1/2✓ Branch 0 taken 1052 times.
✗ Branch 1 not taken.
|
1052 | index_t facet_operand_bits = operand_bit_copy[f]; |
| 2004 | const vec3& p1 = mesh_copy_.facets.point(f,0); | ||
| 2005 | const vec3& p2 = mesh_copy_.facets.point(f,1); | ||
| 2006 | const vec3& p3 = mesh_copy_.facets.point(f,2); | ||
| 2007 |
3/4✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1050 times.
✓ Branch 4 taken 2 times.
|
1052 | if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) { |
| 2008 | 2 | component_inclusion_bits ^= facet_operand_bits; | |
| 2009 | } | ||
| 2010 | } | ||
| 2011 | return component_inclusion_bits; | ||
| 2012 | } | ||
| 2013 | |||
| 2014 | ✗ | index_t MeshSurfaceIntersection::compute_component_inclusion_bits_exact( | |
| 2015 | index_t component, index_t v | ||
| 2016 | ) { | ||
| 2017 | Attribute<index_t> operand_bit_copy( | ||
| 2018 | ✗ | mesh_copy_.facets.attributes(), "operand_bit" | |
| 2019 | ); | ||
| 2020 | |||
| 2021 | Attribute<index_t> facet_component_copy( | ||
| 2022 | ✗ | mesh_copy_.facets.attributes(), "component" | |
| 2023 | ); | ||
| 2024 | |||
| 2025 | ✗ | if(verbose_) { | |
| 2026 | ✗ | Logger::out("Weiler") << " componentE" << component << std::endl; | |
| 2027 | } | ||
| 2028 | |||
| 2029 | index_t component_inclusion_bits = 0; | ||
| 2030 | |||
| 2031 | ✗ | ExactPoint q1 = exact_vertex(v); | |
| 2032 | ✗ | ExactPoint q2 = q1; | |
| 2033 | ✗ | if(q2.w.is_one()) { | |
| 2034 | ✗ | q2.z += exact::scalar(1e6); | |
| 2035 | } else { | ||
| 2036 | ✗ | q2.z += (q2.w * exact::scalar(1e6)); | |
| 2037 | } | ||
| 2038 | |||
| 2039 | ✗ | for(index_t f: mesh_copy_.facets) { | |
| 2040 | ✗ | if(facet_component_copy[f] == NO_INDEX) { | |
| 2041 | ✗ | continue; | |
| 2042 | } | ||
| 2043 | |||
| 2044 | ✗ | if(facet_component_copy[f] == component) { | |
| 2045 | ✗ | continue; | |
| 2046 | } | ||
| 2047 | |||
| 2048 | ✗ | index_t facet_operand_bits = operand_bit_copy[f]; | |
| 2049 | ✗ | ExactPoint p1(mesh_copy_.facets.point(f,0)); | |
| 2050 | ✗ | ExactPoint p2(mesh_copy_.facets.point(f,1)); | |
| 2051 | ✗ | ExactPoint p3(mesh_copy_.facets.point(f,2)); | |
| 2052 | ✗ | if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) { | |
| 2053 | ✗ | component_inclusion_bits ^= facet_operand_bits; | |
| 2054 | } | ||
| 2055 | ✗ | } | |
| 2056 | ✗ | return component_inclusion_bits; | |
| 2057 | ✗ | } | |
| 2058 | |||
| 2059 | index_t | ||
| 2060 | ✗ | MeshSurfaceIntersection::compute_component_inclusion_bits_exact_exact( | |
| 2061 | index_t component, index_t v | ||
| 2062 | ) { | ||
| 2063 | Attribute<index_t> operand_bit( | ||
| 2064 | ✗ | mesh_.facets.attributes(), "operand_bit" | |
| 2065 | ); | ||
| 2066 | |||
| 2067 | Attribute<index_t> facet_component( | ||
| 2068 | ✗ | mesh_.facets.attributes(), "component" | |
| 2069 | ); | ||
| 2070 | |||
| 2071 | ✗ | if(verbose_) { | |
| 2072 | ✗ | Logger::out("Weiler") << " componentEE" << component << std::endl; | |
| 2073 | } | ||
| 2074 | |||
| 2075 | index_t component_inclusion_bits = 0; | ||
| 2076 | |||
| 2077 | ✗ | ExactPoint q1 = exact_vertex(v); | |
| 2078 | ✗ | ExactPoint q2 = q1; | |
| 2079 | ✗ | if(q2.w.is_one()) { | |
| 2080 | ✗ | q2.z += exact::scalar(1e6); | |
| 2081 | } else { | ||
| 2082 | ✗ | q2.z += (q2.w * exact::scalar(1e6)); | |
| 2083 | } | ||
| 2084 | |||
| 2085 | ✗ | for(index_t f: mesh_.facets) { | |
| 2086 | ✗ | if(facet_component[f] == NO_INDEX) { | |
| 2087 | ✗ | continue; | |
| 2088 | } | ||
| 2089 | |||
| 2090 | ✗ | if(facet_component[f] == component) { | |
| 2091 | ✗ | continue; | |
| 2092 | } | ||
| 2093 | |||
| 2094 | ✗ | index_t facet_operand_bits = operand_bit[f]; | |
| 2095 | ✗ | ExactPoint p1 = exact_vertex(mesh_.facets.vertex(f,0)); | |
| 2096 | ✗ | ExactPoint p2 = exact_vertex(mesh_.facets.vertex(f,1)); | |
| 2097 | ✗ | ExactPoint p3 = exact_vertex(mesh_.facets.vertex(f,2)); | |
| 2098 | ✗ | if(segment_triangle_intersection_SOS(q1,q2,p1,p2,p3)) { | |
| 2099 | ✗ | component_inclusion_bits ^= facet_operand_bits; | |
| 2100 | } | ||
| 2101 | ✗ | } | |
| 2102 | ✗ | return component_inclusion_bits; | |
| 2103 | ✗ | } | |
| 2104 | |||
| 2105 | /*************************************************************************/ | ||
| 2106 | |||
| 2107 | 57 | void MeshSurfaceIntersection::simplify_coplanar_facets( | |
| 2108 | double angle_tolerance | ||
| 2109 | ) { | ||
| 2110 |
1/2✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
|
57 | if(mesh_.facets.nb() == 0) { |
| 2111 | ✗ | return; | |
| 2112 | } | ||
| 2113 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | if(interpolate_attributes_) { |
| 2114 | ✗ | Logger::warn("Intersect") | |
| 2115 | << "Cannot simplify coplanar facets with interpolated attributes" | ||
| 2116 | << std::endl; | ||
| 2117 | ✗ | return; | |
| 2118 | } | ||
| 2119 | |||
| 2120 |
1/2✓ Branch 2 taken 57 times.
✗ Branch 3 not taken.
|
57 | Stopwatch W("Coplanar",verbose_); |
| 2121 |
2/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
|
114 | Attribute<index_t> facet_group(mesh_.facets.attributes(), "group"); |
| 2122 | vector<index_t> group_facet; // one facet per group | ||
| 2123 |
2/2✓ Branch 0 taken 123492 times.
✓ Branch 1 taken 57 times.
|
123549 | for(index_t f: mesh_.facets) { |
| 2124 | 123492 | facet_group[f] = NO_INDEX; | |
| 2125 | } | ||
| 2126 |
2/6✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
57 | Attribute<bool> keep_vertex(mesh_.vertices.attributes(), "keep"); |
| 2127 |
2/2✓ Branch 0 taken 61752 times.
✓ Branch 1 taken 57 times.
|
61809 | for(index_t v: mesh_.vertices) { |
| 2128 | keep_vertex[v] = false; | ||
| 2129 | } | ||
| 2130 | 57 | index_t current_group = 0; | |
| 2131 | { | ||
| 2132 | // clear attributes -------------v | ||
| 2133 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | CoplanarFacets coplanar(*this, true, angle_tolerance); |
| 2134 |
4/4✓ Branch 0 taken 73084 times.
✓ Branch 1 taken 50408 times.
✓ Branch 2 taken 123492 times.
✓ Branch 3 taken 57 times.
|
123549 | for(index_t f: mesh_.facets) { |
| 2135 |
2/2✓ Branch 0 taken 73084 times.
✓ Branch 1 taken 50408 times.
|
123492 | if(facet_group[f] == NO_INDEX) { |
| 2136 |
1/2✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
|
73084 | coplanar.get(f, current_group); // This sets facet_group_[f] |
| 2137 |
1/2✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
|
73084 | coplanar.mark_vertices_to_keep(); |
| 2138 | group_facet.push_back(f); | ||
| 2139 | 73084 | ++current_group; | |
| 2140 | } | ||
| 2141 | } | ||
| 2142 | 57 | } | |
| 2143 | |||
| 2144 |
1/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
57 | vector<index_t> remove_f(mesh_.facets.nb(), 0); |
| 2145 | 57 | index_t nb_groups = current_group; | |
| 2146 | geo_debug_assert(nb_groups == group_facet.size()); | ||
| 2147 | |||
| 2148 | // Avoid to have reallocations in parallel with access by | ||
| 2149 | // preallocating facets (we are going to create a maximum | ||
| 2150 | // nb of facets that corresponds to the actual nb of facets) | ||
| 2151 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | mesh_.facets.reserve(mesh_.facets.nb()); |
| 2152 | |||
| 2153 | // Triangulate coplanar facet groups in parallel | ||
| 2154 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | Process::spinlock lock = GEOGRAM_SPINLOCK_INIT; |
| 2155 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | parallel_for_slice( |
| 2156 | 285 | 0, nb_groups, [&](index_t b, index_t e) { | |
| 2157 | // do not clear attributes -----v | ||
| 2158 | 228 | CoplanarFacets coplanar(*this,false,angle_tolerance); | |
| 2159 | |||
| 2160 |
2/2✓ Branch 0 taken 73084 times.
✓ Branch 1 taken 228 times.
|
73312 | for(index_t group=b; group<e; ++group) { |
| 2161 |
1/2✓ Branch 1 taken 73084 times.
✗ Branch 2 not taken.
|
73084 | coplanar.get(group_facet[group],group); |
| 2162 | |||
| 2163 |
2/2✓ Branch 0 taken 65017 times.
✓ Branch 1 taken 8067 times.
|
73084 | if(coplanar.nb_facets() < 2) { |
| 2164 | 65017 | continue; | |
| 2165 | } | ||
| 2166 | |||
| 2167 |
1/2✓ Branch 1 taken 8067 times.
✗ Branch 2 not taken.
|
8067 | coplanar.triangulate(); |
| 2168 | |||
| 2169 | bool OK = true; | ||
| 2170 |
2/2✓ Branch 0 taken 27543 times.
✓ Branch 1 taken 8067 times.
|
35610 | for(index_t t=0; t<coplanar.CDT.nT(); ++t) { |
| 2171 | // If one of v1,v2,v3 is NO_INDEX, | ||
| 2172 | // it means that v1,v2 or v3 was one of the four | ||
| 2173 | // vertices of the external quad. | ||
| 2174 | // It means that there was probably an | ||
| 2175 | // inside/outside classification error. | ||
| 2176 | |||
| 2177 | index_t v1=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,0)); | ||
| 2178 | index_t v2=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,1)); | ||
| 2179 | index_t v3=coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,2)); | ||
| 2180 | 27543 | OK = OK && (v1 != NO_INDEX); | |
| 2181 | 27543 | OK = OK && (v2 != NO_INDEX); | |
| 2182 | 27543 | OK = OK && (v3 != NO_INDEX); | |
| 2183 | } | ||
| 2184 | |||
| 2185 |
1/2✓ Branch 0 taken 8067 times.
✗ Branch 1 not taken.
|
8067 | if(OK) { |
| 2186 | 8067 | coplanar.mark_facets(remove_f); | |
| 2187 | 8067 | Process::acquire_spinlock(lock); | |
| 2188 |
2/2✓ Branch 0 taken 27543 times.
✓ Branch 1 taken 8067 times.
|
63153 | for(index_t t=0; t<coplanar.CDT.nT(); ++t) { |
| 2189 | index_t v1= | ||
| 2190 | coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,0)); | ||
| 2191 | index_t v2= | ||
| 2192 | coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,1)); | ||
| 2193 | index_t v3= | ||
| 2194 | coplanar.CDT.vertex_id(coplanar.CDT.Tv(t,2)); | ||
| 2195 | index_t new_f = | ||
| 2196 |
1/2✓ Branch 1 taken 27543 times.
✗ Branch 2 not taken.
|
27543 | mesh_.facets.create_triangle(v1,v2,v3); |
| 2197 | 27543 | facet_group[new_f] = current_group; | |
| 2198 | } | ||
| 2199 | 8067 | Process::release_spinlock(lock); | |
| 2200 | } | ||
| 2201 | } | ||
| 2202 | 228 | } | |
| 2203 | ); | ||
| 2204 | |||
| 2205 | // Delete temporary attributes | ||
| 2206 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | facet_group.destroy(); |
| 2207 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | keep_vertex.destroy(); |
| 2208 | |||
| 2209 |
1/4✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
57 | remove_f.resize(mesh_.facets.nb(),0); |
| 2210 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | mesh_.facets.delete_elements(remove_f); |
| 2211 |
1/2✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
|
57 | mesh_.facets.connect(); |
| 2212 | 57 | } | |
| 2213 | |||
| 2214 | |||
| 2215 | |||
| 2216 | } | ||
| 2217 | |||
| 2218 | /******************************************************************************/ | ||
| 2219 | |||
| 2220 | namespace { | ||
| 2221 | using namespace GEO; | ||
| 2222 | |||
| 2223 | ✗ | void copy_operand(Mesh& result, const Mesh& operand, index_t operand_id) { | |
| 2224 | Attribute<index_t> operand_bit( | ||
| 2225 | ✗ | result.facets.attributes(), "operand_bit" | |
| 2226 | ); | ||
| 2227 | index_t v_ofs = result.vertices.create_vertices(operand.vertices.nb()); | ||
| 2228 | ✗ | for(index_t v: operand.vertices) { | |
| 2229 | ✗ | result.vertices.point(v + v_ofs) = operand.vertices.point(v); | |
| 2230 | } | ||
| 2231 | ✗ | for(index_t f1: operand.facets) { | |
| 2232 | index_t N = operand.facets.nb_vertices(f1); | ||
| 2233 | ✗ | index_t f2 = result.facets.create_polygon(N); | |
| 2234 | ✗ | for(index_t lv=0; lv<N; ++lv) { | |
| 2235 | ✗ | result.facets.set_vertex( | |
| 2236 | f2,lv,operand.facets.vertex(f1,lv) + v_ofs | ||
| 2237 | ); | ||
| 2238 | } | ||
| 2239 | ✗ | operand_bit[f2] = index_t(1) << operand_id; | |
| 2240 | } | ||
| 2241 | ✗ | } | |
| 2242 | } | ||
| 2243 | |||
| 2244 | namespace GEO { | ||
| 2245 | |||
| 2246 | ✗ | void mesh_boolean_operation( | |
| 2247 | Mesh& result, const Mesh& A, const Mesh& B, | ||
| 2248 | const std::string& operation, | ||
| 2249 | MeshBooleanOperationFlags flags | ||
| 2250 | ) { | ||
| 2251 | bool verbose = ((flags & MESH_BOOL_OPS_VERBOSE) != 0); | ||
| 2252 | ✗ | if(&result == &A) { | |
| 2253 | Attribute<index_t> operand_bit( | ||
| 2254 | ✗ | result.facets.attributes(), "operand_bit" | |
| 2255 | ); | ||
| 2256 | ✗ | for(index_t f: A.facets) { | |
| 2257 | ✗ | operand_bit[f] = index_t(1); | |
| 2258 | } | ||
| 2259 | ✗ | copy_operand(result,B,1); | |
| 2260 | ✗ | } else if(&result == &B) { | |
| 2261 | ✗ | mesh_boolean_operation( | |
| 2262 | ✗ | result, B, A, (operation=="A-B") ? "B-A" : operation, verbose | |
| 2263 | ); | ||
| 2264 | ✗ | return; | |
| 2265 | } else { | ||
| 2266 | ✗ | result.clear(); | |
| 2267 | ✗ | result.vertices.set_dimension(3); | |
| 2268 | ✗ | copy_operand(result,A,0); | |
| 2269 | ✗ | copy_operand(result,B,1); | |
| 2270 | } | ||
| 2271 | ✗ | MeshSurfaceIntersection I(result); | |
| 2272 | I.set_radial_sort(true); | ||
| 2273 | I.set_verbose(verbose); | ||
| 2274 | ✗ | if((flags & MESH_BOOL_OPS_ATTRIBS) != 0) { | |
| 2275 | I.set_interpolate_attributes(true); | ||
| 2276 | } | ||
| 2277 | ✗ | if((flags & MESH_BOOL_OPS_NO_CHECK_NEIGHBORS) != 0) { | |
| 2278 | I.set_detect_intersecting_neighbors(false); | ||
| 2279 | } | ||
| 2280 | ✗ | I.intersect(); | |
| 2281 | ✗ | I.classify(operation); | |
| 2282 | ✗ | if( | |
| 2283 | ✗ | (flags & MESH_BOOL_OPS_ATTRIBS) == 0 && | |
| 2284 | (flags & MESH_BOOL_OPS_NO_SIMPLIFY) == 0 | ||
| 2285 | ) { | ||
| 2286 | ✗ | I.simplify_coplanar_facets(); | |
| 2287 | } | ||
| 2288 | ✗ | } | |
| 2289 | |||
| 2290 | 1 | void mesh_remove_intersections(Mesh& M, index_t max_iter, bool verbose) { | |
| 2291 | // TODO: same as tet_meshing() (compute union) ? | ||
| 2292 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
|
4 | for(index_t k=0; k<max_iter; ++k) { |
| 2293 | 3 | MeshSurfaceIntersection I(M); | |
| 2294 | I.set_radial_sort(false); | ||
| 2295 | I.set_verbose(verbose); | ||
| 2296 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | I.intersect(); |
| 2297 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | mesh_repair(M); |
| 2298 | 3 | } | |
| 2299 | 1 | } | |
| 2300 | |||
| 2301 | ✗ | bool mesh_facets_have_intersection(Mesh& M, index_t f1, index_t f2) { | |
| 2302 | ✗ | for(auto [ p1, p2, p3] : M.facets.triangle_points(f1)) { | |
| 2303 | ✗ | for(auto [ q1, q2, q3] : M.facets.triangle_points(f2)) { | |
| 2304 | ✗ | if(triangles_intersections(p1,p2,p3,q1,q2,q3)) { | |
| 2305 | ✗ | return true; | |
| 2306 | } | ||
| 2307 | } | ||
| 2308 | } | ||
| 2309 | ✗ | return false; | |
| 2310 | } | ||
| 2311 | } | ||
| 2312 |