From 28bf1d4037496397e3bc5d69ce51d8ac9b8a2738 Mon Sep 17 00:00:00 2001 From: Howard Trickey Date: Sun, 2 May 2021 16:37:05 -0400 Subject: Fix T87554 Exact Boolean performance bug. There was a quadratic algorithm extracting triangles from a coplanar cluster. This is now linear. Also found and fixed a bug in the same area related to the triangulator added recently: it didn't get the right correspondence between new edges and original edges. --- source/blender/blenlib/BLI_mesh_intersect.hh | 3 + source/blender/blenlib/intern/mesh_boolean.cc | 217 +---------- source/blender/blenlib/intern/mesh_intersect.cc | 483 +++++++++++++++++++----- 3 files changed, 395 insertions(+), 308 deletions(-) diff --git a/source/blender/blenlib/BLI_mesh_intersect.hh b/source/blender/blenlib/BLI_mesh_intersect.hh index a7996939bb1..7000349c5af 100644 --- a/source/blender/blenlib/BLI_mesh_intersect.hh +++ b/source/blender/blenlib/BLI_mesh_intersect.hh @@ -357,6 +357,9 @@ IMesh trimesh_nary_intersect(const IMesh &tm_in, bool use_self, IMeshArena *arena); +/** Return an IMesh that is a triangulation of a mesh with general polygonal faces. */ +IMesh triangulate_polymesh(IMesh &imesh, IMeshArena *arena); + /** This has the side effect of populating verts in the #IMesh. */ void write_obj_mesh(IMesh &m, const std::string &objname); diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc index 16f2220af4c..25291b8c3b0 100644 --- a/source/blender/blenlib/intern/mesh_boolean.cc +++ b/source/blender/blenlib/intern/mesh_boolean.cc @@ -38,7 +38,6 @@ # include "BLI_math_mpq.hh" # include "BLI_mesh_intersect.hh" # include "BLI_mpq3.hh" -# include "BLI_polyfill_2d.h" # include "BLI_set.hh" # include "BLI_span.hh" # include "BLI_stack.hh" @@ -2680,224 +2679,10 @@ static IMesh raycast_patches_boolean(const IMesh &tm, } } BLI_bvhtree_free(tree); - ans.set_faces(out_faces); - return ans; -} -/** - * Tessellate face f into triangles and return an array of `const Face *` - * giving that triangulation. Intended to be used when f has > 4 vertices. - * Care is taken so that the original edge index associated with - * each edge in the output triangles either matches the original edge - * for the (identical) edge of f, or else is -1. So diagonals added - * for triangulation can later be identified by having #NO_INDEX for original. - * - * This code uses Blenlib's BLI_polyfill_calc to do triangulation, and is therefore quite fast. - * Unfortunately, it can product degenerate triangles that mesh_intersect will remove, leaving - * the mesh non-PWN. - */ -static Array polyfill_triangulate_poly(Face *f, IMeshArena *arena) -{ - /* Similar to loop body in BM_mesh_calc_tesselation. */ - int flen = f->size(); - BLI_assert(flen > 4); - if (!f->plane_populated()) { - f->populate_plane(false); - } - /* Project along negative face normal so (x,y) can be used in 2d. */ - const double3 &poly_normal = f->plane->norm; - float no[3] = {float(poly_normal[0]), float(poly_normal[1]), float(poly_normal[2])}; - normalize_v3(no); - float axis_mat[3][3]; - float(*projverts)[2]; - unsigned int(*tris)[3]; - const int totfilltri = flen - 2; - /* Prepare projected vertices and array to receive triangles in tesselation. */ - tris = static_cast(MEM_malloc_arrayN(totfilltri, sizeof(*tris), __func__)); - projverts = static_cast(MEM_malloc_arrayN(flen, sizeof(*projverts), __func__)); - axis_dominant_v3_to_m3_negate(axis_mat, no); - for (int j = 0; j < flen; ++j) { - const double3 &dco = (*f)[j]->co; - float co[3] = {float(dco[0]), float(dco[1]), float(dco[2])}; - mul_v2_m3v3(projverts[j], axis_mat, co); - } - BLI_polyfill_calc(projverts, flen, 1, tris); - /* Put tesselation triangles into Face form. Record original edges where they exist. */ - Array ans(totfilltri); - for (int t = 0; t < totfilltri; ++t) { - unsigned int *tri = tris[t]; - int eo[3]; - const Vert *v[3]; - for (int k = 0; k < 3; k++) { - BLI_assert(tri[k] < flen); - v[k] = (*f)[tri[k]]; - /* If tri edge goes between two successive indices in - * the original face, then it is an original edge. */ - if ((tri[k] + 1) % flen == tri[(k + 1) % 3]) { - eo[k] = f->edge_orig[tri[k]]; - } - else { - eo[k] = NO_INDEX; - } - ans[t] = arena->add_face( - {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false}); - } - } - - MEM_freeN(tris); - MEM_freeN(projverts); - - return ans; -} - -/** - * Which CDT output edge index is for an edge between output verts - * v1 and v2 (in either order)? - * \return -1 if none. - */ -static int find_cdt_edge(const CDT_result &cdt_out, int v1, int v2) -{ - for (int e : cdt_out.edge.index_range()) { - const std::pair &edge = cdt_out.edge[e]; - if ((edge.first == v1 && edge.second == v2) || (edge.first == v2 && edge.second == v1)) { - return e; - } - } - return -1; -} - -/** - * Tessellate face f into triangles and return an array of `const Face *` - * giving that triangulation. - * Care is taken so that the original edge index associated with - * each edge in the output triangles either matches the original edge - * for the (identical) edge of f, or else is -1. So diagonals added - * for triangulation can later be identified by having #NO_INDEX for original. - * - * The method used is to use the CDT triangulation. Usually that triangulation - * will only use the existing vertices. However, if the face self-intersects - * then the CDT triangulation will include the intersection points. - * If this happens, we use the polyfill triangulator instead. We don't - * use the polyfill triangulator by default because it can create degenerate - * triangles (which we can handle but they'll create non-manifold meshes). - */ -static Array triangulate_poly(Face *f, IMeshArena *arena) -{ - int flen = f->size(); - CDT_input cdt_in; - cdt_in.vert = Array(flen); - cdt_in.face = Array>(1); - cdt_in.face[0].reserve(flen); - for (int i : f->index_range()) { - cdt_in.face[0].append(i); - } - /* Project poly along dominant axis of normal to get 2d coords. */ - if (!f->plane_populated()) { - f->populate_plane(false); - } - const double3 &poly_normal = f->plane->norm; - int axis = double3::dominant_axis(poly_normal); - /* If project down y axis as opposed to x or z, the orientation - * of the polygon will be reversed. - * Yet another reversal happens if the poly normal in the dominant - * direction is opposite that of the positive dominant axis. */ - bool rev1 = (axis == 1); - bool rev2 = poly_normal[axis] < 0; - bool rev = rev1 ^ rev2; - for (int i = 0; i < flen; ++i) { - int ii = rev ? flen - i - 1 : i; - mpq2 &p2d = cdt_in.vert[ii]; - int k = 0; - for (int j = 0; j < 3; ++j) { - if (j != axis) { - p2d[k++] = (*f)[ii]->co_exact[j]; - } - } - } - CDT_result cdt_out = delaunay_2d_calc(cdt_in, CDT_INSIDE); - int n_tris = cdt_out.face.size(); - Array ans(n_tris); - for (int t = 0; t < n_tris; ++t) { - int i_v_out[3]; - const Vert *v[3]; - int eo[3]; - bool needs_steiner = false; - for (int i = 0; i < 3; ++i) { - i_v_out[i] = cdt_out.face[t][i]; - if (cdt_out.vert_orig[i_v_out[i]].size() == 0) { - needs_steiner = true; - break; - } - v[i] = (*f)[cdt_out.vert_orig[i_v_out[i]][0]]; - } - if (needs_steiner) { - /* Fall back on the polyfill triangulator. */ - return polyfill_triangulate_poly(f, arena); - } - for (int i = 0; i < 3; ++i) { - int e_out = find_cdt_edge(cdt_out, i_v_out[i], i_v_out[(i + 1) % 3]); - BLI_assert(e_out != -1); - eo[i] = NO_INDEX; - for (int orig : cdt_out.edge_orig[e_out]) { - if (orig != NO_INDEX) { - eo[i] = orig; - break; - } - } - } - if (rev) { - ans[t] = arena->add_face( - {v[0], v[2], v[1]}, f->orig, {eo[2], eo[1], eo[0]}, {false, false, false}); - } - else { - ans[t] = arena->add_face( - {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false}); - } - } + ans.set_faces(out_faces); return ans; } - -/** - * Return an #IMesh that is a triangulation of a mesh with general - * polygonal faces, #IMesh. - * Added diagonals will be distinguishable by having edge original - * indices of #NO_INDEX. - */ -static IMesh triangulate_polymesh(IMesh &imesh, IMeshArena *arena) -{ - Vector face_tris; - constexpr int estimated_tris_per_face = 3; - face_tris.reserve(estimated_tris_per_face * imesh.face_size()); - for (Face *f : imesh.faces()) { - /* Tessellate face f, following plan similar to #BM_face_calc_tesselation. */ - int flen = f->size(); - if (flen == 3) { - face_tris.append(f); - } - else if (flen == 4) { - const Vert *v0 = (*f)[0]; - const Vert *v1 = (*f)[1]; - const Vert *v2 = (*f)[2]; - const Vert *v3 = (*f)[3]; - int eo_01 = f->edge_orig[0]; - int eo_12 = f->edge_orig[1]; - int eo_23 = f->edge_orig[2]; - int eo_30 = f->edge_orig[3]; - Face *f0 = arena->add_face({v0, v1, v2}, f->orig, {eo_01, eo_12, -1}, {false, false, false}); - Face *f1 = arena->add_face({v0, v2, v3}, f->orig, {-1, eo_23, eo_30}, {false, false, false}); - face_tris.append(f0); - face_tris.append(f1); - } - else { - Array tris = triangulate_poly(f, arena); - for (Face *tri : tris) { - face_tris.append(tri); - } - } - } - return IMesh(face_tris); -} - /** * If \a tri1 and \a tri2 have a common edge (in opposite orientation), * return the indices into \a tri1 and \a tri2 where that common edge starts. Else return (-1,-1). diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc index ce3a5b55f98..3558dfad81d 100644 --- a/source/blender/blenlib/intern/mesh_intersect.cc +++ b/source/blender/blenlib/intern/mesh_intersect.cc @@ -39,6 +39,7 @@ # include "BLI_math_mpq.hh" # include "BLI_mpq2.hh" # include "BLI_mpq3.hh" +# include "BLI_polyfill_2d.h" # include "BLI_set.hh" # include "BLI_span.hh" # include "BLI_task.h" @@ -1576,6 +1577,9 @@ struct CDT_data { Vector is_reversed; /** Result of running CDT on input with (vert, edge, face). */ CDT_result cdt_out; + /** To speed up get_cdt_edge_orig, sometimes populate this map from vertex pair to output edge. + */ + Map, int> verts_to_edge; int proj_axis; }; @@ -1734,6 +1738,32 @@ static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm, return ans; } +/* Return a copy of the argument with the integers ordered in ascending order. */ +static inline std::pair sorted_int_pair(std::pair pair) +{ + if (pair.first <= pair.second) { + return pair; + } + else { + return std::pair(pair.second, pair.first); + } +} + +/** + * Build cd.verts_to_edge to map from a pair of cdt output indices to an index in cd.cdt_out.edge. + * Order the vertex indices so that the smaller one is first in the pair. + */ +static void populate_cdt_edge_map(Map, int> &verts_to_edge, + const CDT_result &cdt_out) +{ + verts_to_edge.reserve(cdt_out.edge.size()); + for (int e : cdt_out.edge.index_range()) { + std::pair vpair = sorted_int_pair(cdt_out.edge[e]); + /* There should be only one edge for each vertex pair. */ + verts_to_edge.add(vpair, e); + } +} + /** * Fills in cd.cdt_out with result of doing the cdt calculation on (vert, edge, face). */ @@ -1766,6 +1796,10 @@ static void do_cdt(CDT_data &cd) } cdt_in.epsilon = 0; /* TODO: needs attention for non-exact T. */ cd.cdt_out = blender::meshintersect::delaunay_2d_calc(cdt_in, CDT_INSIDE); + constexpr int make_edge_map_threshold = 15; + if (cd.cdt_out.edge.size() >= make_edge_map_threshold) { + populate_cdt_edge_map(cd.verts_to_edge, cd.cdt_out); + } if (dbg_level > 0) { std::cout << "\nCDT result\nVerts:\n"; for (int i : cd.cdt_out.vert.index_range()) { @@ -1802,45 +1836,308 @@ static int get_cdt_edge_orig( { int foff = cd.cdt_out.face_edge_offset; *r_is_intersect = false; - for (int e : cd.cdt_out.edge.index_range()) { - std::pair edge = cd.cdt_out.edge[e]; - if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) { - /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */ - /* TODO: if edge has origs from more than on part of the nary input, - * then want to set *r_is_intersect to true. */ - for (int orig_index : cd.cdt_out.edge_orig[e]) { - /* orig_index encodes the triangle and pos within the triangle of the input edge. */ - if (orig_index >= foff) { - int in_face_index = (orig_index / foff) - 1; - int pos = orig_index % foff; - /* We need to retrieve the edge orig field from the Face used to populate the - * in_face_index'th face of the CDT, at the pos'th position of the face. */ - int in_tm_face_index = cd.input_face[in_face_index]; - BLI_assert(in_tm_face_index < in_tm.face_size()); - const Face *facep = in_tm.face(in_tm_face_index); - BLI_assert(pos < facep->size()); - bool is_rev = cd.is_reversed[in_face_index]; - int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos]; - if (eorig != NO_INDEX) { - return eorig; - } - } - else { - /* This edge came from an edge input to the CDT problem, - * so it is an intersect edge. */ - *r_is_intersect = true; - /* TODO: maybe there is an orig index: - * This happens if an input edge was formed by an input face having - * an edge that is co-planar with the cluster, while the face as a whole is not. */ - return NO_INDEX; - } + int e = NO_INDEX; + if (cd.verts_to_edge.size() > 0) { + /* Use the populated map to find the edge, if any, between vertices i0 and i1. */ + std::pair vpair = sorted_int_pair(std::pair(i0, i1)); + e = cd.verts_to_edge.lookup_default(vpair, NO_INDEX); + } + else { + for (int ee : cd.cdt_out.edge.index_range()) { + std::pair edge = cd.cdt_out.edge[ee]; + if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) { + e = ee; + break; } + } + } + if (e == NO_INDEX) { + return NO_INDEX; + } + + /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */ + /* TODO: if edge has origs from more than on part of the nary input, + * then want to set *r_is_intersect to true. */ + for (int orig_index : cd.cdt_out.edge_orig[e]) { + /* orig_index encodes the triangle and pos within the triangle of the input edge. */ + if (orig_index >= foff) { + int in_face_index = (orig_index / foff) - 1; + int pos = orig_index % foff; + /* We need to retrieve the edge orig field from the Face used to populate the + * in_face_index'th face of the CDT, at the pos'th position of the face. */ + int in_tm_face_index = cd.input_face[in_face_index]; + BLI_assert(in_tm_face_index < in_tm.face_size()); + const Face *facep = in_tm.face(in_tm_face_index); + BLI_assert(pos < facep->size()); + bool is_rev = cd.is_reversed[in_face_index]; + int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos]; + if (eorig != NO_INDEX) { + return eorig; + } + } + else { + /* This edge came from an edge input to the CDT problem, + * so it is an intersect edge. */ + *r_is_intersect = true; + /* TODO: maybe there is an orig index: + * This happens if an input edge was formed by an input face having + * an edge that is co-planar with the cluster, while the face as a whole is not. */ return NO_INDEX; } } return NO_INDEX; } +/** + * Make a Face from the CDT output triangle cdt_out_t, which has corresponding input triangle + * cdt_in_t. The triangle indices that are in cd.input_face are from IMesh tm. + * Return a pointer to the newly constructed Face. + */ +static Face *cdt_tri_as_imesh_face( + int cdt_out_t, int cdt_in_t, const CDT_data &cd, const IMesh &tm, IMeshArena *arena) +{ + const CDT_result &cdt_out = cd.cdt_out; + int t_orig = tm.face(cd.input_face[cdt_in_t])->orig; + BLI_assert(cdt_out.face[cdt_out_t].size() == 3); + int i0 = cdt_out.face[cdt_out_t][0]; + int i1 = cdt_out.face[cdt_out_t][1]; + int i2 = cdt_out.face[cdt_out_t][2]; + mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]); + mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]); + mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]); + /* No need to provide an original index: if coord matches + * an original one, then it will already be in the arena + * with the correct orig field. */ + const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX); + const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX); + const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX); + Face *facep; + bool is_isect0; + bool is_isect1; + bool is_isect2; + if (cd.is_reversed[cdt_in_t]) { + int oe0 = get_cdt_edge_orig(i0, i2, cd, tm, &is_isect0); + int oe1 = get_cdt_edge_orig(i2, i1, cd, tm, &is_isect1); + int oe2 = get_cdt_edge_orig(i1, i0, cd, tm, &is_isect2); + facep = arena->add_face( + {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); + } + else { + int oe0 = get_cdt_edge_orig(i0, i1, cd, tm, &is_isect0); + int oe1 = get_cdt_edge_orig(i1, i2, cd, tm, &is_isect1); + int oe2 = get_cdt_edge_orig(i2, i0, cd, tm, &is_isect2); + facep = arena->add_face( + {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); + } + facep->populate_plane(false); + return facep; +} + +/** + * Tessellate face f into triangles and return an array of `const Face *` + * giving that triangulation. Intended to be used when f has > 4 vertices. + * Care is taken so that the original edge index associated with + * each edge in the output triangles either matches the original edge + * for the (identical) edge of f, or else is -1. So diagonals added + * for triangulation can later be identified by having #NO_INDEX for original. + * + * This code uses Blenlib's BLI_polyfill_calc to do triangulation, and is therefore quite fast. + * Unfortunately, it can product degenerate triangles that mesh_intersect will remove, leaving + * the mesh non-PWN. + */ +static Array polyfill_triangulate_poly(Face *f, IMeshArena *arena) +{ + /* Similar to loop body in BM_mesh_calc_tesselation. */ + int flen = f->size(); + BLI_assert(flen > 4); + if (!f->plane_populated()) { + f->populate_plane(false); + } + /* Project along negative face normal so (x,y) can be used in 2d. */ + const double3 &poly_normal = f->plane->norm; + float no[3] = {float(poly_normal[0]), float(poly_normal[1]), float(poly_normal[2])}; + normalize_v3(no); + float axis_mat[3][3]; + float(*projverts)[2]; + unsigned int(*tris)[3]; + const int totfilltri = flen - 2; + /* Prepare projected vertices and array to receive triangles in tesselation. */ + tris = static_cast(MEM_malloc_arrayN(totfilltri, sizeof(*tris), __func__)); + projverts = static_cast(MEM_malloc_arrayN(flen, sizeof(*projverts), __func__)); + axis_dominant_v3_to_m3_negate(axis_mat, no); + for (int j = 0; j < flen; ++j) { + const double3 &dco = (*f)[j]->co; + float co[3] = {float(dco[0]), float(dco[1]), float(dco[2])}; + mul_v2_m3v3(projverts[j], axis_mat, co); + } + BLI_polyfill_calc(projverts, flen, 1, tris); + /* Put tesselation triangles into Face form. Record original edges where they exist. */ + Array ans(totfilltri); + for (int t = 0; t < totfilltri; ++t) { + unsigned int *tri = tris[t]; + int eo[3]; + const Vert *v[3]; + for (int k = 0; k < 3; k++) { + BLI_assert(tri[k] < flen); + v[k] = (*f)[tri[k]]; + /* If tri edge goes between two successive indices in + * the original face, then it is an original edge. */ + if ((tri[k] + 1) % flen == tri[(k + 1) % 3]) { + eo[k] = f->edge_orig[tri[k]]; + } + else { + eo[k] = NO_INDEX; + } + ans[t] = arena->add_face( + {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false}); + } + } + + MEM_freeN(tris); + MEM_freeN(projverts); + + return ans; +} + +/** + * Tessellate face f into triangles and return an array of `const Face *` + * giving that triangulation. + * Care is taken so that the original edge index associated with + * each edge in the output triangles either matches the original edge + * for the (identical) edge of f, or else is -1. So diagonals added + * for triangulation can later be identified by having #NO_INDEX for original. + * + * The method used is to use the CDT triangulation. Usually that triangulation + * will only use the existing vertices. However, if the face self-intersects + * then the CDT triangulation will include the intersection points. + * If this happens, we use the polyfill triangulator instead. We don't + * use the polyfill triangulator by default because it can create degenerate + * triangles (which we can handle but they'll create non-manifold meshes). + */ +static Array triangulate_poly(Face *f, IMeshArena *arena) +{ + int flen = f->size(); + CDT_input cdt_in; + cdt_in.vert = Array(flen); + cdt_in.face = Array>(1); + cdt_in.face[0].reserve(flen); + for (int i : f->index_range()) { + cdt_in.face[0].append(i); + } + /* Project poly along dominant axis of normal to get 2d coords. */ + if (!f->plane_populated()) { + f->populate_plane(false); + } + const double3 &poly_normal = f->plane->norm; + int axis = double3::dominant_axis(poly_normal); + /* If project down y axis as opposed to x or z, the orientation + * of the polygon will be reversed. + * Yet another reversal happens if the poly normal in the dominant + * direction is opposite that of the positive dominant axis. */ + bool rev1 = (axis == 1); + bool rev2 = poly_normal[axis] < 0; + bool rev = rev1 ^ rev2; + for (int i = 0; i < flen; ++i) { + int ii = rev ? flen - i - 1 : i; + mpq2 &p2d = cdt_in.vert[ii]; + int k = 0; + for (int j = 0; j < 3; ++j) { + if (j != axis) { + p2d[k++] = (*f)[ii]->co_exact[j]; + } + } + } + CDT_result cdt_out = delaunay_2d_calc(cdt_in, CDT_INSIDE); + int n_tris = cdt_out.face.size(); + Array ans(n_tris); + for (int t = 0; t < n_tris; ++t) { + int i_v_out[3]; + const Vert *v[3]; + int eo[3]; + bool needs_steiner = false; + for (int i = 0; i < 3; ++i) { + i_v_out[i] = cdt_out.face[t][i]; + if (cdt_out.vert_orig[i_v_out[i]].size() == 0) { + needs_steiner = true; + break; + } + v[i] = (*f)[cdt_out.vert_orig[i_v_out[i]][0]]; + } + if (needs_steiner) { + /* Fall back on the polyfill triangulator. */ + return polyfill_triangulate_poly(f, arena); + } + Map, int> verts_to_edge; + populate_cdt_edge_map(verts_to_edge, cdt_out); + int foff = cdt_out.face_edge_offset; + for (int i = 0; i < 3; ++i) { + std::pair vpair(i_v_out[i], i_v_out[(i + 1) % 3]); + std::pair vpair_canon = sorted_int_pair(vpair); + int e_out = verts_to_edge.lookup_default(vpair_canon, NO_INDEX); + BLI_assert(e_out != NO_INDEX); + eo[i] = NO_INDEX; + for (int orig : cdt_out.edge_orig[e_out]) { + if (orig >= foff) { + int pos = orig % foff; + BLI_assert(pos < f->size()); + eo[i] = f->edge_orig[pos]; + break; + } + } + } + if (rev) { + ans[t] = arena->add_face( + {v[0], v[2], v[1]}, f->orig, {eo[2], eo[1], eo[0]}, {false, false, false}); + } + else { + ans[t] = arena->add_face( + {v[0], v[1], v[2]}, f->orig, {eo[0], eo[1], eo[2]}, {false, false, false}); + } + } + return ans; +} + +/** + * Return an #IMesh that is a triangulation of a mesh with general + * polygonal faces, #IMesh. + * Added diagonals will be distinguishable by having edge original + * indices of #NO_INDEX. + */ +IMesh triangulate_polymesh(IMesh &imesh, IMeshArena *arena) +{ + Vector face_tris; + constexpr int estimated_tris_per_face = 3; + face_tris.reserve(estimated_tris_per_face * imesh.face_size()); + for (Face *f : imesh.faces()) { + /* Tessellate face f, following plan similar to #BM_face_calc_tesselation. */ + int flen = f->size(); + if (flen == 3) { + face_tris.append(f); + } + else if (flen == 4) { + const Vert *v0 = (*f)[0]; + const Vert *v1 = (*f)[1]; + const Vert *v2 = (*f)[2]; + const Vert *v3 = (*f)[3]; + int eo_01 = f->edge_orig[0]; + int eo_12 = f->edge_orig[1]; + int eo_23 = f->edge_orig[2]; + int eo_30 = f->edge_orig[3]; + Face *f0 = arena->add_face({v0, v1, v2}, f->orig, {eo_01, eo_12, -1}, {false, false, false}); + Face *f1 = arena->add_face({v0, v2, v3}, f->orig, {-1, eo_23, eo_30}, {false, false, false}); + face_tris.append(f0); + face_tris.append(f1); + } + else { + Array tris = triangulate_poly(f, arena); + for (Face *tri : tris) { + face_tris.append(tri); + } + } + } + return IMesh(face_tris); +} + /** * Using the result of CDT in cd.cdt_out, extract an #IMesh representing the subdivision * of input triangle t, which should be an element of cd.input_face. @@ -1862,55 +2159,17 @@ static IMesh extract_subdivided_tri(const CDT_data &cd, BLI_assert(false); return IMesh(); } - int t_orig = in_tm.face(t)->orig; constexpr int inline_buf_size = 20; Vector faces; for (int f : cdt_out.face.index_range()) { if (cdt_out.face_orig[f].contains(t_in_cdt)) { - BLI_assert(cdt_out.face[f].size() == 3); - int i0 = cdt_out.face[f][0]; - int i1 = cdt_out.face[f][1]; - int i2 = cdt_out.face[f][2]; - mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]); - mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]); - mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]); - /* No need to provide an original index: if coord matches - * an original one, then it will already be in the arena - * with the correct orig field. */ - const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX); - const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX); - const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX); - Face *facep; - bool is_isect0; - bool is_isect1; - bool is_isect2; - if (cd.is_reversed[t_in_cdt]) { - int oe0 = get_cdt_edge_orig(i0, i2, cd, in_tm, &is_isect0); - int oe1 = get_cdt_edge_orig(i2, i1, cd, in_tm, &is_isect1); - int oe2 = get_cdt_edge_orig(i1, i0, cd, in_tm, &is_isect2); - facep = arena->add_face( - {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); - } - else { - int oe0 = get_cdt_edge_orig(i0, i1, cd, in_tm, &is_isect0); - int oe1 = get_cdt_edge_orig(i1, i2, cd, in_tm, &is_isect1); - int oe2 = get_cdt_edge_orig(i2, i0, cd, in_tm, &is_isect2); - facep = arena->add_face( - {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); - } - facep->populate_plane(false); + Face *facep = cdt_tri_as_imesh_face(f, t_in_cdt, cd, in_tm, arena); faces.append(facep); } } return IMesh(faces); } -static IMesh extract_single_tri(const IMesh &tm, int t) -{ - Face *f = tm.face(t); - return IMesh({f}); -} - static bool bvhtreeverlap_cmp(const BVHTreeOverlap &a, const BVHTreeOverlap &b) { if (a.indexA < b.indexA) { @@ -2126,7 +2385,7 @@ static void calc_overlap_itts(Map, ITT_value> &itt_map, } /** - * Data needed for parallelization of calc_subdivided_tris. + * Data needed for parallelization of calc_subdivided_non_cluster_tris. */ struct OverlapTriRange { int tri_index; @@ -2203,14 +2462,13 @@ static void calc_subdivided_tri_range_func(void *__restrict userdata, * r_tri_subdivided with the result of intersecting it with * all the other triangles in the mesh, if it intersects any others. * But don't do this for triangles that are part of a cluster. - * Also, do nothing here if the answer is just the triangle itself. */ -static void calc_subdivided_tris(Array &r_tri_subdivided, - const IMesh &tm, - const Map, ITT_value> &itt_map, - const CoplanarClusterInfo &clinfo, - const TriOverlaps &ov, - IMeshArena *arena) +static void calc_subdivided_non_cluster_tris(Array &r_tri_subdivided, + const IMesh &tm, + const Map, ITT_value> &itt_map, + const CoplanarClusterInfo &clinfo, + const TriOverlaps &ov, + IMeshArena *arena) { const int dbg_level = 0; if (dbg_level > 0) { @@ -2250,6 +2508,54 @@ static void calc_subdivided_tris(Array &r_tri_subdivided, settings.use_threading = intersect_use_threading; BLI_task_parallel_range( 0, overlap_tri_range_tot, &data, calc_subdivided_tri_range_func, &settings); + /* Now have to put in the triangles that are the same as the input ones, and not in clusters. + */ + for (int t : tm.face_index_range()) { + if (r_tri_subdivided[t].face_size() == 0 && clinfo.tri_cluster(t) == NO_INDEX) { + r_tri_subdivided[t] = IMesh({tm.face(t)}); + } + } +} + +/** + * For each cluster in clinfo, extract the triangles from the cluster + * that correspond to each original triangle t that is part of the cluster, + * and put the resulting triangles into an IMesh in tri_subdivided[t]. + * We have already done the CDT for the triangles in the cluster, whose + * result is in cluster_subdivided[c] for each cluster c. + */ +static void calc_cluster_tris(Array &tri_subdivided, + const IMesh &tm, + const CoplanarClusterInfo &clinfo, + const Array &cluster_subdivided, + IMeshArena *arena) +{ + for (int c : clinfo.index_range()) { + const CoplanarCluster &cl = clinfo.cluster(c); + const CDT_data &cd = cluster_subdivided[c]; + /* Each triangle in cluster c should be an input triangle in cd.input_faces. + * (See prepare_cdt_input_for_cluster.) + * So accumulate a Vector of Face* for each input face by going through the + * output faces and making a Face for each input face that it is part of. + * (The Boolean algorithm wants duplicates if a given output triangle is part + * of more than one input triangle.) + */ + int n_cluster_tris = cl.tot_tri(); + const CDT_result &cdt_out = cd.cdt_out; + BLI_assert(cd.input_face.size() == n_cluster_tris); + Array> face_vec(n_cluster_tris); + for (int cdt_out_t : cdt_out.face.index_range()) { + for (int cdt_in_t : cdt_out.face_orig[cdt_out_t]) { + Face *f = cdt_tri_as_imesh_face(cdt_out_t, cdt_in_t, cd, tm, arena); + face_vec[cdt_in_t].append(f); + } + } + for (int cdt_in_t : cd.input_face.index_range()) { + int tm_t = cd.input_face[cdt_in_t]; + BLI_assert(tri_subdivided[tm_t].face_size() == 0); + tri_subdivided[tm_t] = IMesh(face_vec[cdt_in_t]); + } + } } static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo, @@ -2622,10 +2928,11 @@ IMesh trimesh_nary_intersect(const IMesh &tm_in, doperfmax(2, tri_ov.overlap().size()); # endif Array tri_subdivided(tm_clean->face_size()); - calc_subdivided_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena); + calc_subdivided_non_cluster_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena); # ifdef PERFDEBUG double subdivided_tris_time = PIL_check_seconds_timer(); - std::cout << "subdivided tris found, time = " << subdivided_tris_time - itt_time << "\n"; + std::cout << "subdivided non-cluster tris found, time = " << subdivided_tris_time - itt_time + << "\n"; # endif Array cluster_subdivided(clinfo.tot_cluster()); for (int c : clinfo.index_range()) { @@ -2636,19 +2943,11 @@ IMesh trimesh_nary_intersect(const IMesh &tm_in, std::cout << "subdivided clusters found, time = " << cluster_subdivide_time - subdivided_tris_time << "\n"; # endif - for (int t : tm_clean->face_index_range()) { - int c = clinfo.tri_cluster(t); - if (c != NO_INDEX) { - BLI_assert(tri_subdivided[t].face_size() == 0); - tri_subdivided[t] = extract_subdivided_tri(cluster_subdivided[c], *tm_clean, t, arena); - } - else if (tri_subdivided[t].face_size() == 0) { - tri_subdivided[t] = extract_single_tri(*tm_clean, t); - } - } + calc_cluster_tris(tri_subdivided, *tm_clean, clinfo, cluster_subdivided, arena); # ifdef PERFDEBUG double extract_time = PIL_check_seconds_timer(); - std::cout << "triangles extracted, time = " << extract_time - cluster_subdivide_time << "\n"; + std::cout << "subdivided cluster tris found, time = " << extract_time - cluster_subdivide_time + << "\n"; # endif IMesh combined = union_tri_subdivides(tri_subdivided); if (dbg_level > 1) { -- cgit v1.2.3