| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2025 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/en/bruno-levy-1 | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #include <geogram/mesh/mesh_minkowski.h> | ||
| 41 | #include <geogram/mesh/mesh.h> | ||
| 42 | #include <geogram/mesh/mesh_geometry.h> | ||
| 43 | #include <geogram/mesh/mesh_convex_hull.h> | ||
| 44 | #include <geogram/mesh/mesh_surface_intersection.h> | ||
| 45 | #include <geogram/mesh/mesh_topology.h> | ||
| 46 | #include <geogram/mesh/mesh_fill_holes.h> | ||
| 47 | #include <geogram/delaunay/delaunay.h> | ||
| 48 | #include <geogram/numerics/predicates.h> | ||
| 49 | |||
| 50 | // Implement an exact version of Hichem Barki's algorithm: | ||
| 51 | // 1) 2D convex hull: implement QuickHull, make a generic version that has | ||
| 52 | // its own orient2d_compare() predicate has a "plug-in" | ||
| 53 | |||
| 54 | namespace { | ||
| 55 | using namespace GEO; | ||
| 56 | |||
| 57 | /** | ||
| 58 | * \brief Convexity of an edge in a 3D mesh | ||
| 59 | * \param[in] M a const reference to the mesh | ||
| 60 | * \param[in] f1 an edge facet | ||
| 61 | * \param[in] le1 a local edge index in \p f1 | ||
| 62 | * \param[in] angle_tolerance maximum normal angle deviation for | ||
| 63 | * considering that the two adjacent facets are co-planar | ||
| 64 | * \retval POSITIVE if the edge is convex | ||
| 65 | * \retval ZERO if the two facets adjacent to the edge are coplanar, | ||
| 66 | * up to \p angle_tolerance | ||
| 67 | * \retval NEGATIVE if the edge is concave | ||
| 68 | */ | ||
| 69 | ✗ | Sign edge_convexity( | |
| 70 | const Mesh& M, index_t f1, index_t le1, | ||
| 71 | double angle_tolerance = 0.0 | ||
| 72 | ) { | ||
| 73 | // Tiny tolerance (0.01 degrees) for nearly co-planar facets (like | ||
| 74 | // square facets in CSGBuilder's tesselated spheres). | ||
| 75 | ✗ | index_t n1 = M.facets.nb_vertices(f1); | |
| 76 | ✗ | index_t f2 = M.facets.adjacent(f1,le1); | |
| 77 | ✗ | geo_debug_assert(f2 != NO_INDEX); | |
| 78 | ✗ | index_t n2 = M.facets.nb_vertices(f2); | |
| 79 | ✗ | index_t le2 = M.facets.find_adjacent(f2,f1); | |
| 80 | ✗ | vec3 p1 = M.vertices.point(M.facets.vertex(f1,le1)); | |
| 81 | ✗ | vec3 p2 = M.vertices.point(M.facets.vertex(f1,(le1+1)%n1)); | |
| 82 | ✗ | vec3 p3 = M.vertices.point(M.facets.vertex(f1,(le1+2)%n1)); | |
| 83 | ✗ | vec3 p4 = M.vertices.point(M.facets.vertex(f2,(le2+2)%n2)); | |
| 84 | ✗ | if(angle_tolerance != 0.0) { | |
| 85 | ✗ | vec3 N1 = cross(p2-p1, p3-p1); | |
| 86 | ✗ | vec3 N2 = cross(p4-p1, p2-p1); | |
| 87 | ✗ | double a1 = Geom::angle(N1,N2); | |
| 88 | ✗ | if(::fabs(a1) < angle_tolerance) { | |
| 89 | ✗ | return ZERO; | |
| 90 | } | ||
| 91 | } | ||
| 92 | ✗ | return PCK::orient_3d(p1,p2,p3,p4); | |
| 93 | } | ||
| 94 | |||
| 95 | ✗ | bool mesh_is_convex_3d(const Mesh& M, double angle_tolerance = 0.0) { | |
| 96 | ✗ | index_t nb_positive = 0; | |
| 97 | ✗ | index_t nb_negative = 0; | |
| 98 | ✗ | for(index_t f: M.facets) { | |
| 99 | ✗ | for(index_t le=0; le<M.facets.nb_vertices(f); ++le) { | |
| 100 | ✗ | index_t g = M.facets.adjacent(f,le); | |
| 101 | ✗ | geo_debug_assert(g != NO_INDEX); | |
| 102 | ✗ | if(f > g) { | |
| 103 | ✗ | Sign s = edge_convexity(M,f,le,angle_tolerance); | |
| 104 | ✗ | nb_positive += (s > 0); | |
| 105 | ✗ | nb_negative += (s < 0); | |
| 106 | } | ||
| 107 | } | ||
| 108 | } | ||
| 109 | ✗ | return (nb_positive == 0) || (nb_negative == 0); | |
| 110 | } | ||
| 111 | |||
| 112 | ✗ | bool mesh_is_convex_2d(const Mesh& M) { | |
| 113 | ✗ | vector<index_t> nxt(M.vertices.nb(), NO_INDEX); | |
| 114 | ✗ | for(index_t e: M.edges) { | |
| 115 | ✗ | index_t v1 = M.edges.vertex(e,0); | |
| 116 | ✗ | index_t v2 = M.edges.vertex(e,1); | |
| 117 | ✗ | geo_debug_assert(nxt[v1] == NO_INDEX); | |
| 118 | ✗ | nxt[v1] = v2; | |
| 119 | } | ||
| 120 | |||
| 121 | ✗ | index_t nb_positive = 0; | |
| 122 | ✗ | index_t nb_negative = 0; | |
| 123 | ✗ | for(index_t v1: M.vertices) { | |
| 124 | ✗ | index_t v2 = nxt[v1]; | |
| 125 | ✗ | index_t v3 = nxt[v2]; | |
| 126 | ✗ | Sign s = PCK::orient_2d( | |
| 127 | M.vertices.point_ptr(v1), | ||
| 128 | M.vertices.point_ptr(v2), | ||
| 129 | M.vertices.point_ptr(v3) | ||
| 130 | ); | ||
| 131 | ✗ | nb_positive += (s > 0); | |
| 132 | ✗ | nb_negative += (s < 0); | |
| 133 | } | ||
| 134 | ✗ | return (nb_positive == 0) || (nb_negative == 0); | |
| 135 | ✗ | } | |
| 136 | |||
| 137 | ✗ | void compute_minkowski_sum_convex_convex_3d( | |
| 138 | Mesh& result, const Mesh& op1, const Mesh& op2 | ||
| 139 | ) { | ||
| 140 | ✗ | result.clear(); | |
| 141 | ✗ | result.vertices.set_dimension(3); | |
| 142 | |||
| 143 | ✗ | result.vertices.create_vertices(op1.vertices.nb()*op2.vertices.nb()); | |
| 144 | ✗ | for(index_t v1: op1.vertices) { | |
| 145 | ✗ | for(index_t v2: op2.vertices) { | |
| 146 | ✗ | result.vertices.point(v1*op2.vertices.nb()+v2) = | |
| 147 | ✗ | op1.vertices.point(v1) + op2.vertices.point(v2) ; | |
| 148 | } | ||
| 149 | } | ||
| 150 | ✗ | compute_convex_hull_3d(result); | |
| 151 | ✗ | MeshSurfaceIntersection I(result); | |
| 152 | // TODO: make it work without binding original_facet_id | ||
| 153 | Attribute<index_t> original_facet_id( | ||
| 154 | ✗ | result.facets.attributes(), "original_facet_id" | |
| 155 | ✗ | ); | |
| 156 | ✗ | for(index_t f: result.facets) { | |
| 157 | ✗ | original_facet_id[f] = f; | |
| 158 | } | ||
| 159 | ✗ | I.simplify_coplanar_facets(); | |
| 160 | ✗ | } | |
| 161 | |||
| 162 | /** | ||
| 163 | * \brief Reference: Hichem Barki Ph.D. thesis, Algo 2 P. 90 | ||
| 164 | * A is possibly non-convex | ||
| 165 | * B is convex | ||
| 166 | */ | ||
| 167 | class Minkovski { | ||
| 168 | public: | ||
| 169 | ✗ | Minkovski( | |
| 170 | const Mesh& A, const Mesh& B, Mesh& result | ||
| 171 | ✗ | ) : A_(A), B_(B), result_(result) { | |
| 172 | |||
| 173 | ✗ | Av_f_.assign(A_.facets.nb(), NO_INDEX); | |
| 174 | ✗ | Af_to_Bvcontrib_.resize(A_.facets.nb()); | |
| 175 | |||
| 176 | ✗ | A_has_borders_ = false; | |
| 177 | // Store one facet incident to each vertex for vertex_is_elevated() | ||
| 178 | ✗ | for(index_t f: A_.facets) { | |
| 179 | ✗ | for(index_t v: A_.facets.vertices(f)) { | |
| 180 | ✗ | Av_f_[v] = f; | |
| 181 | } | ||
| 182 | ✗ | for(index_t g: A_.facets.adjacent(f)) { | |
| 183 | ✗ | A_has_borders_ = (A_has_borders_ || g == NO_INDEX); | |
| 184 | } | ||
| 185 | } | ||
| 186 | ✗ | } | |
| 187 | |||
| 188 | ✗ | void compute() { | |
| 189 | ✗ | result_.clear(); | |
| 190 | ✗ | result_.vertices.set_dimension(3); | |
| 191 | |||
| 192 | // Translated facets | ||
| 193 | ✗ | vector<index_t> afv; | |
| 194 | ✗ | for(index_t af: A_.facets) { | |
| 195 | ✗ | afv.resize(0); | |
| 196 | ✗ | for(index_t v: A_.facets.vertices(af)) { | |
| 197 | ✗ | afv.push_back(v); | |
| 198 | } | ||
| 199 | |||
| 200 | ✗ | vec3 ap1 = A_.facets.point(af, 0); | |
| 201 | ✗ | vec3 ap2 = A_.facets.point(af, 1); | |
| 202 | ✗ | vec3 ap3 = A_.facets.point(af, 2); | |
| 203 | |||
| 204 | ✗ | vec3 N = cross(ap2-ap1,ap3-ap1); | |
| 205 | ✗ | if(length(N) < 1e-6) { | |
| 206 | ✗ | continue; | |
| 207 | } | ||
| 208 | |||
| 209 | ✗ | vector<index_t>& af_bvcontrib = Af_to_Bvcontrib_[af]; | |
| 210 | ✗ | for(index_t bv: B_.vertices) { | |
| 211 | ✗ | if(af_bvcontrib.size() == 0) { | |
| 212 | ✗ | af_bvcontrib.push_back(bv); | |
| 213 | } else { | ||
| 214 | ✗ | vec3 bp = B_.vertices.point(bv); | |
| 215 | ✗ | vec3 best_bp_so_far = B_.vertices.point(af_bvcontrib[0]); | |
| 216 | ✗ | Sign s = N_dot_compare( | |
| 217 | ap1, ap2, ap3, bp, best_bp_so_far | ||
| 218 | ); | ||
| 219 | ✗ | if(s >= 0) { | |
| 220 | ✗ | if(s > 0) { | |
| 221 | ✗ | af_bvcontrib.resize(0); | |
| 222 | } | ||
| 223 | ✗ | af_bvcontrib.push_back(bv); | |
| 224 | } | ||
| 225 | } | ||
| 226 | } | ||
| 227 | ✗ | minkowski_2d(afv, af_bvcontrib); | |
| 228 | } | ||
| 229 | |||
| 230 | // Corner facets | ||
| 231 | ✗ | for(index_t bf: B_.facets) { | |
| 232 | ✗ | vec3 bp1 = B_.facets.point(bf, 0); | |
| 233 | ✗ | vec3 bp2 = B_.facets.point(bf, 1); | |
| 234 | ✗ | vec3 bp3 = B_.facets.point(bf, 2); | |
| 235 | ✗ | vec3 bfN = cross(bp2-bp1,bp3-bp1); | |
| 236 | ✗ | if(length(bfN) < 1e-6) { | |
| 237 | ✗ | continue; | |
| 238 | } | ||
| 239 | ✗ | for(index_t av: A_.vertices) { | |
| 240 | ✗ | if(A_vertex_is_elevated(av, bp1, bp2, bp3)) { | |
| 241 | ✗ | vec3 T = A_.vertices.point(av); | |
| 242 | ✗ | index_t deg = B_.facets.nb_vertices(bf); | |
| 243 | ✗ | index_t first_v = result_.vertices.create_vertices(deg); | |
| 244 | ✗ | index_t new_f = result_.facets.create_polygon(deg); | |
| 245 | ✗ | for(index_t lv=0; lv<deg; ++lv) { | |
| 246 | ✗ | result_.facets.set_vertex(new_f, lv, first_v+lv); | |
| 247 | ✗ | result_.vertices.point(first_v+lv) = | |
| 248 | ✗ | B_.facets.point(bf,lv) + T; | |
| 249 | } | ||
| 250 | } | ||
| 251 | } | ||
| 252 | // TODO: edges of op1 elevated w.r.t. V | ||
| 253 | } | ||
| 254 | |||
| 255 | // Edge facets | ||
| 256 | ✗ | for_each_edge(A_, [&](const Edge& aE) { | |
| 257 | ✗ | if(!A_facets_have_distinct_contributing_vertices(aE.f, aE.g)) { | |
| 258 | ✗ | return; | |
| 259 | } | ||
| 260 | ✗ | if(edge_convexity(A_, aE.f, aE.le) == POSITIVE) { | |
| 261 | ✗ | return; | |
| 262 | } | ||
| 263 | ✗ | for_each_edge(B_, [&](const Edge& bE) { | |
| 264 | // B is almost always a tessellated sphere with | ||
| 265 | // quads made of two triangles. Ignore the diagonal | ||
| 266 | // edges of those quads. | ||
| 267 | //if(Geom::angle(bE.Nf, bE.Ng) < 0.01 * M_PI / 180.0) { | ||
| 268 | // return; | ||
| 269 | //} | ||
| 270 | |||
| 271 | ✗ | Sign s1 = geo_sgn(dot(aE.dir,bE.Nf)); | |
| 272 | ✗ | Sign s2 = geo_sgn(dot(aE.dir,bE.Ng)); | |
| 273 | |||
| 274 | ✗ | if(s1 == s2) { | |
| 275 | ✗ | return; | |
| 276 | } | ||
| 277 | |||
| 278 | ✗ | vec3 abN = cross(aE.dir,bE.dir); | |
| 279 | ✗ | if( | |
| 280 | ✗ | double(s1)*dot(cross(aE.Nf,abN),aE.dir) < 0 && | |
| 281 | ✗ | double(s1)*dot(cross(abN,aE.Ng),aE.dir) < 0 | |
| 282 | ) { | ||
| 283 | ✗ | create_quad( | |
| 284 | ✗ | aE.p1+bE.p1, aE.p1+bE.p2, aE.p2+bE.p2, aE.p2+bE.p1 | |
| 285 | ); | ||
| 286 | } | ||
| 287 | }); | ||
| 288 | }); | ||
| 289 | |||
| 290 | ✗ | tessellate_facets(result_, 3); | |
| 291 | ✗ | } | |
| 292 | |||
| 293 | |||
| 294 | protected: | ||
| 295 | |||
| 296 | struct Edge { | ||
| 297 | ✗ | Edge( | |
| 298 | const Mesh& mesh_in, index_t f_in, index_t le_in | ||
| 299 | ✗ | ) : mesh(mesh_in), f(f_in), le(le_in) { | |
| 300 | ✗ | g = mesh.facets.adjacent(f,le); | |
| 301 | ✗ | if(!valid()) { | |
| 302 | ✗ | return; | |
| 303 | } | ||
| 304 | ✗ | v1 = mesh.facets.vertex(f,le); | |
| 305 | ✗ | v2 = mesh.facets.vertex(f,(le+1)%mesh.facets.nb_vertices(f)); | |
| 306 | ✗ | p1 = mesh.vertices.point(v1); | |
| 307 | ✗ | p2 = mesh.vertices.point(v2); | |
| 308 | ✗ | dir = p2-p1; | |
| 309 | ✗ | Nf = Geom::mesh_facet_normal(mesh,f); | |
| 310 | ✗ | Ng = Geom::mesh_facet_normal(mesh,g); | |
| 311 | } | ||
| 312 | ✗ | bool valid() { | |
| 313 | ✗ | return g != NO_INDEX && f > g; | |
| 314 | } | ||
| 315 | const Mesh& mesh; | ||
| 316 | index_t f; | ||
| 317 | index_t le; | ||
| 318 | index_t g; | ||
| 319 | index_t v1; | ||
| 320 | index_t v2; | ||
| 321 | vec3 p1; | ||
| 322 | vec3 p2; | ||
| 323 | vec3 dir; | ||
| 324 | vec3 Nf; | ||
| 325 | vec3 Ng; | ||
| 326 | }; | ||
| 327 | |||
| 328 | ✗ | void for_each_edge( | |
| 329 | const Mesh& M, std::function<void(const Edge& E)> doit | ||
| 330 | ) { | ||
| 331 | ✗ | for(index_t f: M.facets) { | |
| 332 | ✗ | for(index_t le=0; le<M.facets.nb_vertices(f); ++le) { | |
| 333 | ✗ | Edge E(M,f,le); | |
| 334 | ✗ | if(E.valid()) { | |
| 335 | ✗ | doit(E); | |
| 336 | } | ||
| 337 | } | ||
| 338 | } | ||
| 339 | ✗ | } | |
| 340 | |||
| 341 | ✗ | void create_quad(vec3 p1, vec3 p2, vec3 p3, vec3 p4) { | |
| 342 | ✗ | index_t first_v =result_.vertices.create_vertices(4); | |
| 343 | ✗ | index_t new_f = result_.facets.create_polygon(4); | |
| 344 | ✗ | for(index_t lv=0; lv<4; ++lv) { | |
| 345 | ✗ | result_.facets.set_vertex(new_f, lv, first_v+lv); | |
| 346 | } | ||
| 347 | ✗ | result_.vertices.point(first_v ) = p1; | |
| 348 | ✗ | result_.vertices.point(first_v+1) = p2; | |
| 349 | ✗ | result_.vertices.point(first_v+2) = p3; | |
| 350 | ✗ | result_.vertices.point(first_v+3) = p4; | |
| 351 | ✗ | } | |
| 352 | |||
| 353 | ✗ | bool A_facets_have_distinct_contributing_vertices( | |
| 354 | index_t Af1, index_t Af2 | ||
| 355 | ) const { | ||
| 356 | ✗ | const vector<index_t>& c1 = Af_to_Bvcontrib_[Af1]; | |
| 357 | ✗ | const vector<index_t>& c2 = Af_to_Bvcontrib_[Af2]; | |
| 358 | ✗ | for(index_t v1: c1) { | |
| 359 | ✗ | if(std::find(c2.begin(), c2.end(), v1) == c2.end()) { | |
| 360 | ✗ | return true; | |
| 361 | } | ||
| 362 | } | ||
| 363 | ✗ | return false; | |
| 364 | } | ||
| 365 | |||
| 366 | ✗ | bool A_vertex_is_elevated( | |
| 367 | index_t av, vec3 bp1, vec3 bp2, vec3 bp3 | ||
| 368 | ) const { | ||
| 369 | ✗ | vec3 ap1 = A_.vertices.point(av); | |
| 370 | ✗ | index_t first_f = Av_f_[av]; | |
| 371 | ✗ | index_t f = first_f; | |
| 372 | ✗ | index_t lv = A_.facets.find_vertex(f,av); | |
| 373 | do { | ||
| 374 | ✗ | index_t N = A_.facets.nb_vertices(f); | |
| 375 | ✗ | index_t v2 = A_.facets.vertex(f, (lv + 1) % N); | |
| 376 | ✗ | vec3 ap2 = A_.vertices.point(v2); | |
| 377 | ✗ | Sign s = N_dot_compare(bp1,bp2,bp3,ap1,ap2); | |
| 378 | ✗ | if(s < 0) { | |
| 379 | ✗ | return false; | |
| 380 | } | ||
| 381 | ✗ | f = A_.facets.adjacent(f,lv); | |
| 382 | ✗ | lv = A_.facets.find_vertex(f,av); | |
| 383 | ✗ | } while(f != first_f); | |
| 384 | ✗ | return true; | |
| 385 | } | ||
| 386 | |||
| 387 | ✗ | void minkowski_2d( | |
| 388 | const vector<index_t>& av, const vector<index_t>& bv | ||
| 389 | ) { | ||
| 390 | ✗ | geo_debug_assert(av.size() >= 3); | |
| 391 | ✗ | geo_debug_assert(bv.size() > 0); | |
| 392 | ✗ | if(bv.size() == 1) { | |
| 393 | // simple translation | ||
| 394 | ✗ | vec3 T = B_.vertices.point(bv[0]); | |
| 395 | ✗ | index_t N = av.size(); | |
| 396 | ✗ | index_t first_v = result_.vertices.create_vertices(N); | |
| 397 | ✗ | index_t new_f = result_.facets.create_polygon(N); | |
| 398 | ✗ | for(index_t lv=0; lv<N; ++lv) { | |
| 399 | ✗ | result_.facets.set_vertex(new_f, lv, first_v+lv); | |
| 400 | ✗ | result_.vertices.point(first_v+lv) = | |
| 401 | ✗ | A_.vertices.point(av[lv]) + T; | |
| 402 | } | ||
| 403 | } else { | ||
| 404 | // Super ugly: use Delaunay in 2D orthogonal plane to | ||
| 405 | // compute convex hull. | ||
| 406 | ✗ | vec3 p1 = A_.vertices.point(av[0]); | |
| 407 | ✗ | vec3 p2 = A_.vertices.point(av[1]); | |
| 408 | ✗ | vec3 p3 = A_.vertices.point(av[2]); | |
| 409 | ✗ | vec3 NN = normalize(cross(p2-p1,p3-p1)); | |
| 410 | ✗ | vec3 U = normalize(Geom::perpendicular(NN)); | |
| 411 | ✗ | vec3 V = cross(NN,U); | |
| 412 | ✗ | vector<vec2> uv; | |
| 413 | ✗ | vector<vec3> xyz; | |
| 414 | ✗ | uv.reserve(av.size()*bv.size()); | |
| 415 | ✗ | xyz.reserve(av.size()*bv.size()); | |
| 416 | ✗ | for(index_t curav: av) { | |
| 417 | ✗ | vec3 ap = A_.vertices.point(curav); | |
| 418 | ✗ | for(index_t curbv: bv) { | |
| 419 | ✗ | vec3 bp = B_.vertices.point(curbv); | |
| 420 | ✗ | vec3 abp = ap+bp-p1; | |
| 421 | ✗ | uv.emplace_back(dot(abp,U), dot(abp,V)); | |
| 422 | ✗ | xyz.emplace_back(ap+bp); | |
| 423 | } | ||
| 424 | } | ||
| 425 | Delaunay_var delaunay = Delaunay::create( | ||
| 426 | ✗ | coord_index_t(2), "BDEL2d" | |
| 427 | ✗ | ); | |
| 428 | ✗ | delaunay->set_keeps_infinite(true); | |
| 429 | ✗ | delaunay->set_vertices(uv.size(), uv[0].data()); | |
| 430 | ✗ | vector<index_t> nxt(uv.size(), NO_INDEX); | |
| 431 | ✗ | index_t first = NO_INDEX; | |
| 432 | ✗ | for(index_t t=delaunay->nb_finite_cells(); | |
| 433 | ✗ | t<delaunay->nb_cells(); ++t) { | |
| 434 | ✗ | index_t v1= NO_INDEX, v2=NO_INDEX; | |
| 435 | ✗ | for(index_t lv=0; lv<3; ++lv) { | |
| 436 | ✗ | if(delaunay->cell_vertex(t,lv) == NO_INDEX) { | |
| 437 | ✗ | v1 = delaunay->cell_vertex(t,(lv+1)%3); | |
| 438 | ✗ | v2 = delaunay->cell_vertex(t,(lv+2)%3); | |
| 439 | } | ||
| 440 | } | ||
| 441 | ✗ | first = v2; | |
| 442 | ✗ | nxt[v2] = v1; | |
| 443 | } | ||
| 444 | ✗ | vector<vec3> f_xyz; | |
| 445 | ✗ | index_t v = first; | |
| 446 | do { | ||
| 447 | ✗ | f_xyz.push_back(xyz[v]); | |
| 448 | ✗ | v = nxt[v]; | |
| 449 | ✗ | } while(v != first); | |
| 450 | |||
| 451 | ✗ | if(f_xyz.size() < 3) { | |
| 452 | ✗ | return; | |
| 453 | } | ||
| 454 | |||
| 455 | ✗ | index_t N = f_xyz.size(); | |
| 456 | ✗ | index_t first_v = result_.vertices.create_vertices(N); | |
| 457 | ✗ | index_t new_f = result_.facets.create_polygon(N); | |
| 458 | ✗ | for(index_t lv=0; lv<N; ++lv) { | |
| 459 | ✗ | result_.facets.set_vertex(new_f, lv, first_v+lv); | |
| 460 | ✗ | result_.vertices.point(first_v+lv) = f_xyz[lv]; | |
| 461 | } | ||
| 462 | ✗ | } | |
| 463 | } | ||
| 464 | |||
| 465 | /** | ||
| 466 | * \brief Compares the dot product between the normal to a triangle | ||
| 467 | * and two vectors | ||
| 468 | * \param[in] p1 , p2 , p3 the three vertices of the triangle, that | ||
| 469 | * define the normal vector N = cross(p2-p1,p3-p1) | ||
| 470 | * \param[in] q1 , q2 the two points to be compared relative to N | ||
| 471 | * \retval POSITIVE if dot(N,q2) > dot(N,q1) | ||
| 472 | * \retval ZERO if dot(N,q2) = dot(N,q1) | ||
| 473 | * \retval NEGATIVE if dot(N,q2) < dot(N,q1) | ||
| 474 | */ | ||
| 475 | ✗ | Sign N_dot_compare( | |
| 476 | const vec3& p1, const vec3& p2, const vec3& p3, | ||
| 477 | const vec3& q1, const vec3& q2 | ||
| 478 | ) const { | ||
| 479 | // TODO: new specialized predicate | ||
| 480 | ✗ | return PCK::det_3d(p3-p1, p2-p1, q2-q1); | |
| 481 | } | ||
| 482 | |||
| 483 | private: | ||
| 484 | const Mesh& A_; | ||
| 485 | const Mesh& B_; | ||
| 486 | bool A_has_borders_; | ||
| 487 | vector<index_t> Av_f_; | ||
| 488 | vector<vector<index_t>> Af_to_Bvcontrib_; | ||
| 489 | Mesh& result_; | ||
| 490 | }; | ||
| 491 | |||
| 492 | ✗ | void compute_minkowski_sum_non_convex_convex_3d( | |
| 493 | Mesh& result, const Mesh& A, const Mesh& B | ||
| 494 | ) { | ||
| 495 | // TODO: why needed ? | ||
| 496 | ✗ | reorient_connected_components(const_cast<Mesh&>(A)); | |
| 497 | ✗ | reorient_connected_components(const_cast<Mesh&>(B)); | |
| 498 | ✗ | Minkovski mink(A,B,result); | |
| 499 | ✗ | mink.compute(); | |
| 500 | ✗ | } | |
| 501 | |||
| 502 | ✗ | void compute_minkowski_sum_convex_convex_2d( | |
| 503 | Mesh& result, const Mesh& op1, const Mesh& op2 | ||
| 504 | ) { | ||
| 505 | ✗ | result.clear(); | |
| 506 | ✗ | result.vertices.set_dimension(2); | |
| 507 | ✗ | result.vertices.create_vertices(op1.vertices.nb()*op2.vertices.nb()); | |
| 508 | ✗ | for(index_t v1: op1.vertices) { | |
| 509 | ✗ | for(index_t v2: op2.vertices) { | |
| 510 | ✗ | result.vertices.point<2>(v1*op2.vertices.nb()+v2) = | |
| 511 | ✗ | op1.vertices.point<2>(v1) + op2.vertices.point<2>(v2) ; | |
| 512 | } | ||
| 513 | } | ||
| 514 | ✗ | compute_convex_hull_2d(result); | |
| 515 | ✗ | } | |
| 516 | |||
| 517 | |||
| 518 | } | ||
| 519 | |||
| 520 | namespace GEO { | ||
| 521 | |||
| 522 | ✗ | void compute_minkowski_sum_3d( | |
| 523 | Mesh& result, const Mesh& op1, const Mesh& op2 | ||
| 524 | ) { | ||
| 525 | ✗ | const double angle_tol = 0.01 * M_PI / 180.0; | |
| 526 | ✗ | bool op1_is_convex = mesh_is_convex_3d(op1, angle_tol); | |
| 527 | ✗ | bool op2_is_convex = mesh_is_convex_3d(op2, angle_tol); | |
| 528 | |||
| 529 | ✗ | if(op1_is_convex && op2_is_convex) { | |
| 530 | ✗ | compute_minkowski_sum_convex_convex_3d(result, op1, op2); | |
| 531 | ✗ | return; | |
| 532 | } | ||
| 533 | |||
| 534 | ✗ | if(!op1_is_convex && op2_is_convex) { | |
| 535 | // op2 may be not exactly convex -> convexify it | ||
| 536 | ✗ | Mesh op2_convex; | |
| 537 | ✗ | op2_convex.copy(op2); | |
| 538 | ✗ | compute_convex_hull_3d(op2_convex); | |
| 539 | ✗ | compute_minkowski_sum_non_convex_convex_3d(result, op1, op2_convex); | |
| 540 | ✗ | return; | |
| 541 | ✗ | } | |
| 542 | |||
| 543 | ✗ | if(op1_is_convex && !op2_is_convex) { | |
| 544 | // op1 may be not exactly convex -> convexify it | ||
| 545 | ✗ | Mesh op1_convex; | |
| 546 | ✗ | op1_convex.copy(op2); | |
| 547 | ✗ | compute_convex_hull_3d(op1_convex); | |
| 548 | ✗ | compute_minkowski_sum_non_convex_convex_3d(result, op2, op1_convex); | |
| 549 | ✗ | return; | |
| 550 | ✗ | } | |
| 551 | |||
| 552 | throw( | ||
| 553 | std::logic_error( | ||
| 554 | "compute_Minkowski_sum_3d: " | ||
| 555 | "non-convex case not implemented yet" | ||
| 556 | ) | ||
| 557 | ✗ | ); | |
| 558 | } | ||
| 559 | |||
| 560 | ✗ | void compute_minkowski_sum_2d( | |
| 561 | Mesh& result, const Mesh& op1, const Mesh& op2 | ||
| 562 | ) { | ||
| 563 | ✗ | if(mesh_is_convex_2d(op1) && mesh_is_convex_2d(op2)) { | |
| 564 | ✗ | compute_minkowski_sum_convex_convex_2d(result, op1, op2); | |
| 565 | ✗ | return; | |
| 566 | } | ||
| 567 | |||
| 568 | throw( | ||
| 569 | std::logic_error( | ||
| 570 | "compute_Minkowski_sum_2d: " | ||
| 571 | "non-convex case not implemented yet" | ||
| 572 | ) | ||
| 573 | ✗ | ); | |
| 574 | } | ||
| 575 | |||
| 576 | } | ||
| 577 |