diff options
Diffstat (limited to 'source/blender/blenlib/intern/mesh_boolean.cc')
| -rw-r--r-- | source/blender/blenlib/intern/mesh_boolean.cc | 3382 |
1 files changed, 3382 insertions, 0 deletions
diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc new file mode 100644 index 00000000000..387d879c5af --- /dev/null +++ b/source/blender/blenlib/intern/mesh_boolean.cc @@ -0,0 +1,3382 @@ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/** \file + * \ingroup bli + */ + +#ifdef WITH_GMP + +# include <algorithm> +# include <fstream> +# include <iostream> + +# include "BLI_array.hh" +# include "BLI_assert.h" +# include "BLI_delaunay_2d.h" +# include "BLI_hash.hh" +# include "BLI_map.hh" +# include "BLI_math.h" +# include "BLI_math_boolean.hh" +# include "BLI_math_mpq.hh" +# include "BLI_mesh_intersect.hh" +# include "BLI_mpq3.hh" +# include "BLI_set.hh" +# include "BLI_span.hh" +# include "BLI_stack.hh" +# include "BLI_vector.hh" +# include "BLI_vector_set.hh" + +# include "BLI_mesh_boolean.hh" + +namespace blender::meshintersect { + +/** + * Edge as two `const` Vert *'s, in a canonical order (lower vert id first). + * We use the Vert id field for hashing to get algorithms + * that yield predictable results from run-to-run and machine-to-machine. + */ +class Edge { + const Vert *v_[2]{nullptr, nullptr}; + + public: + Edge() = default; + Edge(const Vert *v0, const Vert *v1) + { + if (v0->id <= v1->id) { + v_[0] = v0; + v_[1] = v1; + } + else { + v_[0] = v1; + v_[1] = v0; + } + } + + const Vert *v0() const + { + return v_[0]; + } + + const Vert *v1() const + { + return v_[1]; + } + + const Vert *operator[](int i) const + { + return v_[i]; + } + + bool operator==(Edge other) const + { + return v_[0]->id == other.v_[0]->id && v_[1]->id == other.v_[1]->id; + } + + uint64_t hash() const + { + constexpr uint64_t h1 = 33; + uint64_t v0hash = DefaultHash<int>{}(v_[0]->id); + uint64_t v1hash = DefaultHash<int>{}(v_[1]->id); + return v0hash ^ (v1hash * h1); + } +}; + +static std::ostream &operator<<(std::ostream &os, const Edge &e) +{ + if (e.v0() == nullptr) { + BLI_assert(e.v1() == nullptr); + os << "(null,null)"; + } + else { + os << "(" << e.v0() << "," << e.v1() << ")"; + } + return os; +} + +static std::ostream &operator<<(std::ostream &os, const Span<int> &a) +{ + for (int i : a.index_range()) { + os << a[i]; + if (i != a.size() - 1) { + os << " "; + } + } + return os; +} + +static std::ostream &operator<<(std::ostream &os, const Array<int> &iarr) +{ + os << Span<int>(iarr); + return os; +} + +/** Holds information about topology of an #IMesh that is all triangles. */ +class TriMeshTopology : NonCopyable { + /** Triangles that contain a given Edge (either order). */ + Map<Edge, Vector<int> *> edge_tri_; + /** Edges incident on each vertex. */ + Map<const Vert *, Vector<Edge>> vert_edges_; + + public: + TriMeshTopology(const IMesh &tm); + ~TriMeshTopology(); + + /* If e is manifold, return index of the other triangle (not t) that has it. + * Else return NO_INDEX. */ + int other_tri_if_manifold(Edge e, int t) const + { + if (edge_tri_.contains(e)) { + auto *p = edge_tri_.lookup(e); + if (p->size() == 2) { + return ((*p)[0] == t) ? (*p)[1] : (*p)[0]; + } + } + return NO_INDEX; + } + + /* Which triangles share edge e (in either orientation)? */ + const Vector<int> *edge_tris(Edge e) const + { + return edge_tri_.lookup_default(e, nullptr); + } + + /* Which edges are incident on the given vertex? + * We assume v has some incident edges. */ + const Vector<Edge> &vert_edges(const Vert *v) const + { + return vert_edges_.lookup(v); + } + + Map<Edge, Vector<int> *>::ItemIterator edge_tri_map_items() const + { + return edge_tri_.items(); + } +}; + +TriMeshTopology::TriMeshTopology(const IMesh &tm) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "TRIMESHTOPOLOGY CONSTRUCTION\n"; + } + /* If everything were manifold, `F+V-E=2` and `E=3F/2`. + * So an likely overestimate, allowing for non-manifoldness, is `E=2F` and `V=F`. */ + const int estimate_num_edges = 2 * tm.face_size(); + const int estimate_num_verts = tm.face_size(); + edge_tri_.reserve(estimate_num_edges); + vert_edges_.reserve(estimate_num_verts); + for (int t : tm.face_index_range()) { + const Face &tri = *tm.face(t); + BLI_assert(tri.is_tri()); + for (int i = 0; i < 3; ++i) { + const Vert *v = tri[i]; + const Vert *vnext = tri[(i + 1) % 3]; + Edge e(v, vnext); + Vector<Edge> *edges = vert_edges_.lookup_ptr(v); + if (edges == nullptr) { + vert_edges_.add_new(v, Vector<Edge>()); + edges = vert_edges_.lookup_ptr(v); + BLI_assert(edges != nullptr); + } + edges->append_non_duplicates(e); + auto createf = [t](Vector<int> **pvec) { *pvec = new Vector<int>{t}; }; + auto modifyf = [t](Vector<int> **pvec) { (*pvec)->append_non_duplicates(t); }; + this->edge_tri_.add_or_modify(Edge(v, vnext), createf, modifyf); + } + } + /* Debugging. */ + if (dbg_level > 0) { + std::cout << "After TriMeshTopology construction\n"; + for (auto item : edge_tri_.items()) { + std::cout << "tris for edge " << item.key << ": " << *item.value << "\n"; + constexpr bool print_stats = false; + if (print_stats) { + edge_tri_.print_stats(); + } + } + for (auto item : vert_edges_.items()) { + std::cout << "edges for vert " << item.key << ":\n"; + for (const Edge &e : item.value) { + std::cout << " " << e << "\n"; + } + std::cout << "\n"; + } + } +} + +TriMeshTopology::~TriMeshTopology() +{ + for (const Vector<int> *vec : edge_tri_.values()) { + delete vec; + } +} + +/** A Patch is a maximal set of triangles that share manifold edges only. */ +class Patch { + Vector<int> tri_; /* Indices of triangles in the Patch. */ + + public: + Patch() = default; + + void add_tri(int t) + { + tri_.append(t); + } + + int tot_tri() const + { + return tri_.size(); + } + + int tri(int i) const + { + return tri_[i]; + } + + IndexRange tri_range() const + { + return IndexRange(tri_.size()); + } + + Span<int> tris() const + { + return Span<int>(tri_); + } + + int cell_above{NO_INDEX}; + int cell_below{NO_INDEX}; + int component{NO_INDEX}; +}; + +static std::ostream &operator<<(std::ostream &os, const Patch &patch) +{ + os << "Patch " << patch.tris(); + if (patch.cell_above != NO_INDEX) { + os << " cell_above=" << patch.cell_above; + } + else { + os << " cell_above not set"; + } + if (patch.cell_below != NO_INDEX) { + os << " cell_below=" << patch.cell_below; + } + else { + os << " cell_below not set"; + } + return os; +} + +class PatchesInfo { + /** All of the Patches for a #IMesh. */ + Vector<Patch> patch_; + /** Patch index for corresponding triangle. */ + Array<int> tri_patch_; + /** Shared edge for incident patches; (-1, -1) if none. */ + Map<std::pair<int, int>, Edge> pp_edge_; + + public: + explicit PatchesInfo(int ntri) + { + constexpr int max_expected_patch_patch_incidences = 100; + tri_patch_ = Array<int>(ntri, NO_INDEX); + pp_edge_.reserve(max_expected_patch_patch_incidences); + } + + int tri_patch(int t) const + { + return tri_patch_[t]; + } + + int add_patch() + { + int patch_index = patch_.append_and_get_index(Patch()); + return patch_index; + } + + void grow_patch(int patch_index, int t) + { + tri_patch_[t] = patch_index; + patch_[patch_index].add_tri(t); + } + + bool tri_is_assigned(int t) const + { + return tri_patch_[t] != NO_INDEX; + } + + const Patch &patch(int patch_index) const + { + return patch_[patch_index]; + } + + Patch &patch(int patch_index) + { + return patch_[patch_index]; + } + + int tot_patch() const + { + return patch_.size(); + } + + IndexRange index_range() const + { + return IndexRange(patch_.size()); + } + + const Patch *begin() const + { + return patch_.begin(); + } + + const Patch *end() const + { + return patch_.end(); + } + + Patch *begin() + { + return patch_.begin(); + } + + Patch *end() + { + return patch_.end(); + } + + void add_new_patch_patch_edge(int p1, int p2, Edge e) + { + pp_edge_.add_new(std::pair<int, int>(p1, p2), e); + pp_edge_.add_new(std::pair<int, int>(p2, p1), e); + } + + Edge patch_patch_edge(int p1, int p2) + { + return pp_edge_.lookup_default(std::pair<int, int>(p1, p2), Edge()); + } +}; + +static bool apply_bool_op(BoolOpType bool_optype, const Array<int> &winding); + +/** + * A Cell is a volume of 3-space, surrounded by patches. + * We will partition all 3-space into Cells. + * One cell, the Ambient cell, contains all other cells. + */ +class Cell { + Vector<int> patches_; + Array<int> winding_; + int merged_to_{NO_INDEX}; + bool winding_assigned_{false}; + /* in_output_volume_ will be true when this cell should be in the output volume. */ + bool in_output_volume_{false}; + /* zero_volume_ will be true when this is a zero-volume cell (inside a stack of identical + * triangles). */ + bool zero_volume_{false}; + + public: + Cell() = default; + + void add_patch(int p) + { + patches_.append(p); + } + + void add_patch_non_duplicates(int p) + { + patches_.append_non_duplicates(p); + } + + const Span<int> patches() const + { + return Span<int>(patches_); + } + + const Span<int> winding() const + { + return Span<int>(winding_); + } + + void init_winding(int winding_len) + { + winding_ = Array<int>(winding_len); + } + + void seed_ambient_winding() + { + winding_.fill(0); + winding_assigned_ = true; + } + + void set_winding_and_in_output_volume(const Cell &from_cell, + int shape, + int delta, + BoolOpType bool_optype) + { + std::copy(from_cell.winding().begin(), from_cell.winding().end(), winding_.begin()); + winding_[shape] += delta; + winding_assigned_ = true; + in_output_volume_ = apply_bool_op(bool_optype, winding_); + } + + bool in_output_volume() const + { + return in_output_volume_; + } + + bool winding_assigned() const + { + return winding_assigned_; + } + + bool zero_volume() const + { + return zero_volume_; + } + + int merged_to() const + { + return merged_to_; + } + + void set_merged_to(int c) + { + merged_to_ = c; + } + + /** + * Call this when it is possible that this Cell has zero volume, + * and if it does, set zero_volume_ to true. + */ + void check_for_zero_volume(const PatchesInfo &pinfo, const IMesh &mesh); +}; + +static std::ostream &operator<<(std::ostream &os, const Cell &cell) +{ + os << "Cell patches " << cell.patches(); + if (cell.winding().size() > 0) { + os << " winding=" << cell.winding(); + os << " in_output_volume=" << cell.in_output_volume(); + } + os << " zv=" << cell.zero_volume(); + return os; +} + +static bool tris_have_same_verts(const IMesh &mesh, int t1, int t2) +{ + const Face &tri1 = *mesh.face(t1); + const Face &tri2 = *mesh.face(t2); + BLI_assert(tri1.size() == 3 && tri2.size() == 3); + if (tri1.vert[0] == tri2.vert[0]) { + return ((tri1.vert[1] == tri2.vert[1] && tri1.vert[2] == tri2.vert[2]) || + (tri1.vert[1] == tri2.vert[2] && tri1.vert[2] == tri2.vert[1])); + } + if (tri1.vert[0] == tri2.vert[1]) { + return ((tri1.vert[1] == tri2.vert[0] && tri1.vert[2] == tri2.vert[2]) || + (tri1.vert[1] == tri2.vert[2] && tri1.vert[2] == tri2.vert[0])); + } + if (tri1.vert[0] == tri2.vert[2]) { + return ((tri1.vert[1] == tri2.vert[0] && tri1.vert[2] == tri2.vert[1]) || + (tri1.vert[1] == tri2.vert[1] && tri1.vert[2] == tri2.vert[0])); + } + return false; +} + +/** + * A Cell will have zero volume if it is bounded by exactly two patches and those + * patches are geometrically identical triangles (perhaps flipped versions of each other). + * If this Cell has zero volume, set its zero_volume_ member to true. + */ +void Cell::check_for_zero_volume(const PatchesInfo &pinfo, const IMesh &mesh) +{ + if (patches_.size() == 2) { + const Patch &p1 = pinfo.patch(patches_[0]); + const Patch &p2 = pinfo.patch(patches_[1]); + if (p1.tot_tri() == 1 && p2.tot_tri() == 1) { + if (tris_have_same_verts(mesh, p1.tri(0), p2.tri(0))) { + zero_volume_ = true; + } + } + } +} + +/* Information about all the Cells. */ +class CellsInfo { + Vector<Cell> cell_; + + public: + CellsInfo() = default; + + int add_cell() + { + int index = cell_.append_and_get_index(Cell()); + return index; + } + + Cell &cell(int c) + { + return cell_[c]; + } + + const Cell &cell(int c) const + { + return cell_[c]; + } + + int tot_cell() const + { + return cell_.size(); + } + + IndexRange index_range() const + { + return cell_.index_range(); + } + + const Cell *begin() const + { + return cell_.begin(); + } + + const Cell *end() const + { + return cell_.end(); + } + + Cell *begin() + { + return cell_.begin(); + } + + Cell *end() + { + return cell_.end(); + } + + void init_windings(int winding_len) + { + for (Cell &cell : cell_) { + cell.init_winding(winding_len); + } + } +}; + +/** + * For Debugging: write a .obj file showing the patch/cell structure or just the cells. + */ +static void write_obj_cell_patch(const IMesh &m, + const CellsInfo &cinfo, + const PatchesInfo &pinfo, + bool cells_only, + const std::string &name) +{ + /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library. + * This is just for developer debugging anyway, + * and should never be called in production Blender. */ +# ifdef _WIN_32 + const char *objdir = BLI_getenv("HOME"); +# else + const char *objdir = "/tmp/"; +# endif + + std::string fname = std::string(objdir) + name + std::string("_cellpatch.obj"); + std::ofstream f; + f.open(fname); + if (!f) { + std::cout << "Could not open file " << fname << "\n"; + return; + } + + /* Copy IMesh so can populate verts. */ + IMesh mm = m; + mm.populate_vert(); + f << "o cellpatch\n"; + for (const Vert *v : mm.vertices()) { + const double3 dv = v->co; + f << "v " << dv[0] << " " << dv[1] << " " << dv[2] << "\n"; + } + if (!cells_only) { + for (int p : pinfo.index_range()) { + f << "g patch" << p << "\n"; + const Patch &patch = pinfo.patch(p); + for (int t : patch.tris()) { + const Face &tri = *mm.face(t); + f << "f "; + for (const Vert *v : tri) { + f << mm.lookup_vert(v) + 1 << " "; + } + f << "\n"; + } + } + } + for (int c : cinfo.index_range()) { + f << "g cell" << c << "\n"; + const Cell &cell = cinfo.cell(c); + for (int p : cell.patches()) { + const Patch &patch = pinfo.patch(p); + for (int t : patch.tris()) { + const Face &tri = *mm.face(t); + f << "f "; + for (const Vert *v : tri) { + f << mm.lookup_vert(v) + 1 << " "; + } + f << "\n"; + } + } + } + f.close(); +} + +static void merge_cells(int merge_to, int merge_from, CellsInfo &cinfo, PatchesInfo &pinfo) +{ + if (merge_to == merge_from) { + return; + } + Cell &merge_from_cell = cinfo.cell(merge_from); + Cell &merge_to_cell = cinfo.cell(merge_to); + int final_merge_to = merge_to; + while (merge_to_cell.merged_to() != NO_INDEX) { + final_merge_to = merge_to_cell.merged_to(); + merge_to_cell = cinfo.cell(final_merge_to); + } + for (Patch &patch : pinfo) { + if (patch.cell_above == merge_from) { + patch.cell_above = final_merge_to; + } + if (patch.cell_below == merge_from) { + patch.cell_below = final_merge_to; + } + } + for (int cell_p : merge_from_cell.patches()) { + merge_to_cell.add_patch_non_duplicates(cell_p); + } + merge_from_cell.set_merged_to(final_merge_to); +} + +/** + * Partition the triangles of \a tm into Patches. + */ +static PatchesInfo find_patches(const IMesh &tm, const TriMeshTopology &tmtopo) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nFIND_PATCHES\n"; + } + int ntri = tm.face_size(); + PatchesInfo pinfo(ntri); + /* Algorithm: Grow patches across manifold edges as long as there are unassigned triangles. */ + Stack<int> cur_patch_grow; + for (int t : tm.face_index_range()) { + if (pinfo.tri_patch(t) == -1) { + cur_patch_grow.push(t); + int cur_patch_index = pinfo.add_patch(); + while (!cur_patch_grow.is_empty()) { + int tcand = cur_patch_grow.pop(); + if (dbg_level > 1) { + std::cout << "pop tcand = " << tcand << "; assigned = " << pinfo.tri_is_assigned(tcand) + << "\n"; + } + if (pinfo.tri_is_assigned(tcand)) { + continue; + } + if (dbg_level > 1) { + std::cout << "grow patch from seed tcand=" << tcand << "\n"; + } + pinfo.grow_patch(cur_patch_index, tcand); + const Face &tri = *tm.face(tcand); + for (int i = 0; i < 3; ++i) { + Edge e(tri[i], tri[(i + 1) % 3]); + int t_other = tmtopo.other_tri_if_manifold(e, tcand); + if (dbg_level > 1) { + std::cout << " edge " << e << " generates t_other=" << t_other << "\n"; + } + if (t_other != NO_INDEX) { + if (!pinfo.tri_is_assigned(t_other)) { + if (dbg_level > 1) { + std::cout << " push t_other = " << t_other << "\n"; + } + cur_patch_grow.push(t_other); + } + } + else { + /* e is non-manifold. Set any patch-patch incidences we can. */ + if (dbg_level > 1) { + std::cout << " e non-manifold case\n"; + } + const Vector<int> *etris = tmtopo.edge_tris(e); + if (etris != nullptr) { + for (int i : etris->index_range()) { + int t_other = (*etris)[i]; + if (t_other != tcand && pinfo.tri_is_assigned(t_other)) { + int p_other = pinfo.tri_patch(t_other); + if (p_other == cur_patch_index) { + continue; + } + if (pinfo.patch_patch_edge(cur_patch_index, p_other).v0() == nullptr) { + pinfo.add_new_patch_patch_edge(cur_patch_index, p_other, e); + if (dbg_level > 1) { + std::cout << "added patch_patch_edge (" << cur_patch_index << "," << p_other + << ") = " << e << "\n"; + } + } + } + } + } + } + } + } + } + } + if (dbg_level > 0) { + std::cout << "\nafter FIND_PATCHES: found " << pinfo.tot_patch() << " patches\n"; + for (int p : pinfo.index_range()) { + std::cout << p << ": " << pinfo.patch(p) << "\n"; + } + if (dbg_level > 1) { + std::cout << "\ntriangle map\n"; + for (int t : tm.face_index_range()) { + std::cout << t << ": patch " << pinfo.tri_patch(t) << "\n"; + } + } + std::cout << "\npatch-patch incidences\n"; + for (int p1 : pinfo.index_range()) { + for (int p2 : pinfo.index_range()) { + Edge e = pinfo.patch_patch_edge(p1, p2); + if (e.v0() != nullptr) { + std::cout << "p" << p1 << " and p" << p2 << " share edge " << e << "\n"; + } + } + } + } + return pinfo; +} + +/** + * If e is an edge in tri, return the vertex that isn't part of tri, + * the "flap" vertex, or nullptr if e is not part of tri. + * Also, e may be reversed in tri. + * Set *r_rev to true if it is reversed, else false. + */ +static const Vert *find_flap_vert(const Face &tri, const Edge e, bool *r_rev) +{ + *r_rev = false; + const Vert *flapv; + if (tri[0] == e.v0()) { + if (tri[1] == e.v1()) { + *r_rev = false; + flapv = tri[2]; + } + else { + if (tri[2] != e.v1()) { + return nullptr; + } + *r_rev = true; + flapv = tri[1]; + } + } + else if (tri[1] == e.v0()) { + if (tri[2] == e.v1()) { + *r_rev = false; + flapv = tri[0]; + } + else { + if (tri[0] != e.v1()) { + return nullptr; + } + *r_rev = true; + flapv = tri[2]; + } + } + else { + if (tri[2] != e.v0()) { + return nullptr; + } + if (tri[0] == e.v1()) { + *r_rev = false; + flapv = tri[1]; + } + else { + if (tri[1] != e.v1()) { + return nullptr; + } + *r_rev = true; + flapv = tri[0]; + } + } + return flapv; +} + +/** + * Triangle \a tri and tri0 share edge e. + * Classify \a tri with respect to tri0 as described in + * sort_tris_around_edge, and return 1, 2, 3, or 4 as \a tri is: + * (1) co-planar with tri0 and on same side of e + * (2) co-planar with tri0 and on opposite side of e + * (3) below plane of tri0 + * (4) above plane of tri0 + * For "above" and "below", we use the orientation of non-reversed + * orientation of tri0. + * Because of the way the intersect mesh was made, we can assume + * that if a triangle is in class 1 then it is has the same flap vert + * as tri0. + */ +static int sort_tris_class(const Face &tri, const Face &tri0, const Edge e) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "classify e = " << e << "\n"; + } + mpq3 a0 = tri0[0]->co_exact; + mpq3 a1 = tri0[1]->co_exact; + mpq3 a2 = tri0[2]->co_exact; + bool rev; + bool rev0; + const Vert *flapv0 = find_flap_vert(tri0, e, &rev0); + const Vert *flapv = find_flap_vert(tri, e, &rev); + if (dbg_level > 0) { + std::cout << " t0 = " << tri0[0] << " " << tri0[1] << " " << tri0[2]; + std::cout << " rev0 = " << rev0 << " flapv0 = " << flapv0 << "\n"; + std::cout << " t = " << tri[0] << " " << tri[1] << " " << tri[2]; + std::cout << " rev = " << rev << " flapv = " << flapv << "\n"; + } + BLI_assert(flapv != nullptr && flapv0 != nullptr); + const mpq3 flap = flapv->co_exact; + /* orient will be positive if flap is below oriented plane of a0,a1,a2. */ + int orient = orient3d(a0, a1, a2, flap); + int ans; + if (orient > 0) { + ans = rev0 ? 4 : 3; + } + else if (orient < 0) { + ans = rev0 ? 3 : 4; + } + else { + ans = flapv == flapv0 ? 1 : 2; + } + if (dbg_level > 0) { + std::cout << " orient = " << orient << " ans = " << ans << "\n"; + } + return ans; +} + +constexpr int EXTRA_TRI_INDEX = INT_MAX; + +/** + * To ensure consistent ordering of co-planar triangles if they happen to be sorted around + * more than one edge, sort the triangle indices in g (in place) by their index -- but also apply + * a sign to the index: positive if the triangle has edge e in the same orientation, + * otherwise negative. + */ +static void sort_by_signed_triangle_index(Vector<int> &g, + const Edge e, + const IMesh &tm, + const Face *extra_tri) +{ + Array<int> signed_g(g.size()); + for (int i : g.index_range()) { + const Face &tri = g[i] == EXTRA_TRI_INDEX ? *extra_tri : *tm.face(g[i]); + bool rev; + find_flap_vert(tri, e, &rev); + signed_g[i] = rev ? -g[i] : g[i]; + } + std::sort(signed_g.begin(), signed_g.end()); + + for (int i : g.index_range()) { + g[i] = abs(signed_g[i]); + } +} + +/** + * Sort the triangles \a tris, which all share edge e, as they appear + * geometrically clockwise when looking down edge e. + * Triangle t0 is the first triangle in the top-level call + * to this recursive routine. The merge step below differs + * for the top level call and all the rest, so this distinguishes those cases. + * Care is taken in the case of duplicate triangles to have + * an ordering that is consistent with that which would happen + * if another edge of the triangle were sorted around. + * + * We sometimes need to do this with an extra triangle that is not part of tm. + * To accommodate this: + * If extra_tri is non-null, then an index of EXTRA_TRI_INDEX should use it for the triangle. + */ +static Array<int> sort_tris_around_edge(const IMesh &tm, + const TriMeshTopology &tmtopo, + const Edge e, + const Span<int> tris, + const int t0, + const Face *extra_tri) +{ + /* Divide and conquer, quick-sort-like sort. + * Pick a triangle t0, then partition into groups: + * (1) co-planar with t0 and on same side of e + * (2) co-planar with t0 and on opposite side of e + * (3) below plane of t0 + * (4) above plane of t0 + * Each group is sorted and then the sorts are merged to give the answer. + * We don't expect the input array to be very large - should typically + * be only 3 or 4 - so OK to make copies of arrays instead of swapping + * around in a single array. */ + const int dbg_level = 0; + if (tris.size() == 0) { + return Array<int>(); + } + if (dbg_level > 0) { + if (t0 == tris[0]) { + std::cout << "\n"; + } + std::cout << "sort_tris_around_edge " << e << "\n"; + std::cout << "tris = " << tris << "\n"; + std::cout << "t0 = " << t0 << "\n"; + } + Vector<int> g1{tris[0]}; + Vector<int> g2; + Vector<int> g3; + Vector<int> g4; + std::array<Vector<int> *, 4> groups = {&g1, &g2, &g3, &g4}; + const Face &triref = *tm.face(tris[0]); + for (int i : tris.index_range()) { + if (i == 0) { + continue; + } + int t = tris[i]; + BLI_assert(t < tm.face_size() || (t == EXTRA_TRI_INDEX && extra_tri != nullptr)); + const Face &tri = (t == EXTRA_TRI_INDEX) ? *extra_tri : *tm.face(t); + if (dbg_level > 2) { + std::cout << "classifying tri " << t << " with respect to " << tris[0] << "\n"; + } + int group_num = sort_tris_class(tri, triref, e); + if (dbg_level > 2) { + std::cout << " classify result : " << group_num << "\n"; + } + groups[group_num - 1]->append(t); + } + if (dbg_level > 1) { + std::cout << "g1 = " << g1 << "\n"; + std::cout << "g2 = " << g2 << "\n"; + std::cout << "g3 = " << g3 << "\n"; + std::cout << "g4 = " << g4 << "\n"; + } + if (g1.size() > 1) { + sort_by_signed_triangle_index(g1, e, tm, extra_tri); + if (dbg_level > 1) { + std::cout << "g1 sorted: " << g1 << "\n"; + } + } + if (g2.size() > 1) { + sort_by_signed_triangle_index(g2, e, tm, extra_tri); + if (dbg_level > 1) { + std::cout << "g2 sorted: " << g2 << "\n"; + } + } + if (g3.size() > 1) { + Array<int> g3sorted = sort_tris_around_edge(tm, tmtopo, e, g3, t0, extra_tri); + std::copy(g3sorted.begin(), g3sorted.end(), g3.begin()); + if (dbg_level > 1) { + std::cout << "g3 sorted: " << g3 << "\n"; + } + } + if (g4.size() > 1) { + Array<int> g4sorted = sort_tris_around_edge(tm, tmtopo, e, g4, t0, extra_tri); + std::copy(g4sorted.begin(), g4sorted.end(), g4.begin()); + if (dbg_level > 1) { + std::cout << "g4 sorted: " << g4 << "\n"; + } + } + int group_tot_size = g1.size() + g2.size() + g3.size() + g4.size(); + Array<int> ans(group_tot_size); + int *p = ans.begin(); + if (tris[0] == t0) { + p = std::copy(g1.begin(), g1.end(), p); + p = std::copy(g4.begin(), g4.end(), p); + p = std::copy(g2.begin(), g2.end(), p); + std::copy(g3.begin(), g3.end(), p); + } + else { + p = std::copy(g3.begin(), g3.end(), p); + p = std::copy(g1.begin(), g1.end(), p); + p = std::copy(g4.begin(), g4.end(), p); + std::copy(g2.begin(), g2.end(), p); + } + if (dbg_level > 0) { + std::cout << "sorted tris = " << ans << "\n"; + } + return ans; +} + +/** + * Find the Cells around edge e. + * This possibly makes new cells in \a cinfo, and sets up the + * bipartite graph edges between cells and patches. + * Will modify \a pinfo and \a cinfo and the patches and cells they contain. + */ +static void find_cells_from_edge(const IMesh &tm, + const TriMeshTopology &tmtopo, + PatchesInfo &pinfo, + CellsInfo &cinfo, + const Edge e) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_CELLS_FROM_EDGE " << e << "\n"; + } + const Vector<int> *edge_tris = tmtopo.edge_tris(e); + BLI_assert(edge_tris != nullptr); + Array<int> sorted_tris = sort_tris_around_edge( + tm, tmtopo, e, Span<int>(*edge_tris), (*edge_tris)[0], nullptr); + + int n_edge_tris = edge_tris->size(); + Array<int> edge_patches(n_edge_tris); + for (int i = 0; i < n_edge_tris; ++i) { + edge_patches[i] = pinfo.tri_patch(sorted_tris[i]); + if (dbg_level > 1) { + std::cout << "edge_patches[" << i << "] = " << edge_patches[i] << "\n"; + } + } + for (int i = 0; i < n_edge_tris; ++i) { + int inext = (i + 1) % n_edge_tris; + int r_index = edge_patches[i]; + int rnext_index = edge_patches[inext]; + Patch &r = pinfo.patch(r_index); + Patch &rnext = pinfo.patch(rnext_index); + bool r_flipped; + bool rnext_flipped; + find_flap_vert(*tm.face(sorted_tris[i]), e, &r_flipped); + find_flap_vert(*tm.face(sorted_tris[inext]), e, &rnext_flipped); + int *r_follow_cell = r_flipped ? &r.cell_below : &r.cell_above; + int *rnext_prev_cell = rnext_flipped ? &rnext.cell_above : &rnext.cell_below; + if (dbg_level > 0) { + std::cout << "process patch pair " << r_index << " " << rnext_index << "\n"; + std::cout << " r_flipped = " << r_flipped << " rnext_flipped = " << rnext_flipped << "\n"; + std::cout << " r_follow_cell (" << (r_flipped ? "below" : "above") + << ") = " << *r_follow_cell << "\n"; + std::cout << " rnext_prev_cell (" << (rnext_flipped ? "above" : "below") + << ") = " << *rnext_prev_cell << "\n"; + } + if (*r_follow_cell == NO_INDEX && *rnext_prev_cell == NO_INDEX) { + /* Neither is assigned: make a new cell. */ + int c = cinfo.add_cell(); + *r_follow_cell = c; + *rnext_prev_cell = c; + Cell &cell = cinfo.cell(c); + cell.add_patch(r_index); + cell.add_patch(rnext_index); + cell.check_for_zero_volume(pinfo, tm); + if (dbg_level > 0) { + std::cout << " made new cell " << c << "\n"; + std::cout << " p" << r_index << "." << (r_flipped ? "cell_below" : "cell_above") << " = c" + << c << "\n"; + std::cout << " p" << rnext_index << "." << (rnext_flipped ? "cell_above" : "cell_below") + << " = c" << c << "\n"; + } + } + else if (*r_follow_cell != NO_INDEX && *rnext_prev_cell == NO_INDEX) { + int c = *r_follow_cell; + *rnext_prev_cell = c; + Cell &cell = cinfo.cell(c); + cell.add_patch(rnext_index); + cell.check_for_zero_volume(pinfo, tm); + if (dbg_level > 0) { + std::cout << " reuse r_follow: p" << rnext_index << "." + << (rnext_flipped ? "cell_above" : "cell_below") << " = c" << c << "\n"; + } + } + else if (*r_follow_cell == NO_INDEX && *rnext_prev_cell != NO_INDEX) { + int c = *rnext_prev_cell; + *r_follow_cell = c; + Cell &cell = cinfo.cell(c); + cell.add_patch(r_index); + cell.check_for_zero_volume(pinfo, tm); + if (dbg_level > 0) { + std::cout << " reuse rnext prev: rprev_p" << r_index << "." + << (r_flipped ? "cell_below" : "cell_above") << " = c" << c << "\n"; + } + } + else { + if (*r_follow_cell != *rnext_prev_cell) { + if (dbg_level > 0) { + std::cout << " merge cell " << *rnext_prev_cell << " into cell " << *r_follow_cell + << "\n"; + } + merge_cells(*r_follow_cell, *rnext_prev_cell, cinfo, pinfo); + } + } + } +} + +/** + * Find the partition of 3-space into Cells. + * This assigns the cell_above and cell_below for each Patch. + */ +static CellsInfo find_cells(const IMesh &tm, const TriMeshTopology &tmtopo, PatchesInfo &pinfo) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nFIND_CELLS\n"; + } + CellsInfo cinfo; + /* For each unique edge shared between patch pairs, process it. */ + Set<Edge> processed_edges; + int np = pinfo.tot_patch(); + for (int p = 0; p < np; ++p) { + for (int q = p + 1; q < np; ++q) { + Edge e = pinfo.patch_patch_edge(p, q); + if (e.v0() != nullptr) { + if (!processed_edges.contains(e)) { + processed_edges.add_new(e); + find_cells_from_edge(tm, tmtopo, pinfo, cinfo, e); + } + } + } + } + /* Some patches may have no cells at this point. These are either: + * (a) a closed manifold patch only incident on itself (sphere, torus, klein bottle, etc.). + * (b) an open manifold patch only incident on itself (has non-manifold boundaries). + * Make above and below cells for these patches. This will create a disconnected patch-cell + * bipartite graph, which will have to be fixed later. */ + for (int p : pinfo.index_range()) { + Patch &patch = pinfo.patch(p); + if (patch.cell_above == NO_INDEX) { + int c = cinfo.add_cell(); + patch.cell_above = c; + Cell &cell = cinfo.cell(c); + cell.add_patch(p); + } + if (patch.cell_below == NO_INDEX) { + int c = cinfo.add_cell(); + patch.cell_below = c; + Cell &cell = cinfo.cell(c); + cell.add_patch(p); + } + } + if (dbg_level > 0) { + std::cout << "\nFIND_CELLS found " << cinfo.tot_cell() << " cells\nCells\n"; + for (int i : cinfo.index_range()) { + std::cout << i << ": " << cinfo.cell(i) << "\n"; + } + std::cout << "Patches\n"; + for (int i : pinfo.index_range()) { + std::cout << i << ": " << pinfo.patch(i) << "\n"; + } + if (dbg_level > 1) { + write_obj_cell_patch(tm, cinfo, pinfo, false, "postfindcells"); + } + } + return cinfo; +} + +/** + * Find the connected patch components (connects are via intermediate cells), and put + * component numbers in each patch. + * Return a Vector of components - each a Vector of the patch ids in the component. + */ +static Vector<Vector<int>> find_patch_components(const CellsInfo &cinfo, PatchesInfo &pinfo) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_PATCH_COMPONENTS\n"; + } + if (pinfo.tot_patch() == 0) { + return Vector<Vector<int>>(); + } + int current_component = 0; + Stack<int> stack; /* Patch indices to visit. */ + Vector<Vector<int>> ans; + for (int pstart : pinfo.index_range()) { + Patch &patch_pstart = pinfo.patch(pstart); + if (patch_pstart.component != NO_INDEX) { + continue; + } + ans.append(Vector<int>()); + ans[current_component].append(pstart); + stack.push(pstart); + patch_pstart.component = current_component; + while (!stack.is_empty()) { + int p = stack.pop(); + Patch &patch = pinfo.patch(p); + BLI_assert(patch.component == current_component); + for (int c : {patch.cell_above, patch.cell_below}) { + for (int pn : cinfo.cell(c).patches()) { + Patch &patch_neighbor = pinfo.patch(pn); + if (patch_neighbor.component == NO_INDEX) { + patch_neighbor.component = current_component; + stack.push(pn); + ans[current_component].append(pn); + } + } + } + } + ++current_component; + } + if (dbg_level > 0) { + std::cout << "found " << ans.size() << " components\n"; + for (int comp : ans.index_range()) { + std::cout << comp << ": " << ans[comp] << "\n"; + } + } + return ans; +} + +/** + * Do all patches have cell_above and cell_below set? + * Is the bipartite graph connected? + */ +static bool patch_cell_graph_ok(const CellsInfo &cinfo, const PatchesInfo &pinfo) +{ + for (int c : cinfo.index_range()) { + const Cell &cell = cinfo.cell(c); + if (cell.merged_to() != NO_INDEX) { + continue; + } + if (cell.patches().size() == 0) { + std::cout << "Patch/Cell graph disconnected at Cell " << c << " with no patches\n"; + return false; + } + for (int p : cell.patches()) { + if (p >= pinfo.tot_patch()) { + std::cout << "Patch/Cell graph has bad patch index at Cell " << c << "\n"; + return false; + } + } + } + for (int p : pinfo.index_range()) { + const Patch &patch = pinfo.patch(p); + if (patch.cell_above == NO_INDEX || patch.cell_below == NO_INDEX) { + std::cout << "Patch/Cell graph disconnected at Patch " << p + << " with one or two missing cells\n"; + return false; + } + if (patch.cell_above >= cinfo.tot_cell() || patch.cell_below >= cinfo.tot_cell()) { + std::cout << "Patch/Cell graph has bad cell index at Patch " << p << "\n"; + return false; + } + } + return true; +} + +/** + * Is trimesh tm PWN ("Piece-wise constant Winding Number")? + * See Zhou et al. paper for exact definition, but roughly + * means that the faces connect so as to form closed volumes. + * The actual definition says that if you calculate the + * generalized winding number of every point not exactly on + * the mesh, it will always be an integer. + * Necessary (but not sufficient) conditions that a mesh be PWN: + * No edges with a non-zero sum of incident face directions. + * I think that cases like Klein bottles are likely to satisfy + * this without being PWN. So this routine will be only + * approximately right. + */ +static bool is_pwn(const IMesh &tm, const TriMeshTopology &tmtopo) +{ + constexpr int dbg_level = 0; + for (auto item : tmtopo.edge_tri_map_items()) { + const Edge &edge = item.key; + int tot_orient = 0; + /* For each face t attached to edge, add +1 if the edge + * is positively in t, and -1 if negatively in t. */ + for (int t : *item.value) { + const Face &face = *tm.face(t); + BLI_assert(face.size() == 3); + for (int i : face.index_range()) { + if (face[i] == edge.v0()) { + if (face[(i + 1) % 3] == edge.v1()) { + ++tot_orient; + } + else { + BLI_assert(face[(i + 3 - 1) % 3] == edge.v1()); + --tot_orient; + } + } + } + } + if (tot_orient != 0) { + if (dbg_level > 0) { + std::cout << "edge causing non-pwn: " << edge << "\n"; + } + return false; + } + } + return true; +} + +/** + * Find which of the cells around edge e contains point p. + * Do this by inserting a dummy triangle containing v and sorting the + * triangles around the edge to find out where in the sort order + * the dummy triangle lies, then finding which cell is between + * the two triangles on either side of the dummy. + */ +static int find_cell_for_point_near_edge(mpq3 p, + const Edge &e, + const IMesh &tm, + const TriMeshTopology &tmtopo, + const PatchesInfo &pinfo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_CELL_FOR_POINT_NEAR_EDGE, p=" << p << " e=" << e << "\n"; + } + const Vector<int> *etris = tmtopo.edge_tris(e); + const Vert *dummy_vert = arena->add_or_find_vert(p, NO_INDEX); + const Face *dummy_tri = arena->add_face({e.v0(), e.v1(), dummy_vert}, + NO_INDEX, + {NO_INDEX, NO_INDEX, NO_INDEX}, + {false, false, false}); + BLI_assert(etris != nullptr); + Array<int> edge_tris(etris->size() + 1); + std::copy(etris->begin(), etris->end(), edge_tris.begin()); + edge_tris[edge_tris.size() - 1] = EXTRA_TRI_INDEX; + Array<int> sorted_tris = sort_tris_around_edge( + tm, tmtopo, e, edge_tris, edge_tris[0], dummy_tri); + if (dbg_level > 0) { + std::cout << "sorted tris = " << sorted_tris << "\n"; + } + int *p_sorted_dummy = std::find(sorted_tris.begin(), sorted_tris.end(), EXTRA_TRI_INDEX); + BLI_assert(p_sorted_dummy != sorted_tris.end()); + int dummy_index = p_sorted_dummy - sorted_tris.begin(); + int prev_tri = (dummy_index == 0) ? sorted_tris[sorted_tris.size() - 1] : + sorted_tris[dummy_index - 1]; + int next_tri = (dummy_index == sorted_tris.size() - 1) ? sorted_tris[0] : + sorted_tris[dummy_index + 1]; + if (dbg_level > 0) { + std::cout << "prev tri to dummy = " << prev_tri << "; next tri to dummy = " << next_tri + << "\n"; + } + const Patch &prev_patch = pinfo.patch(pinfo.tri_patch(prev_tri)); + if (dbg_level > 0) { + std::cout << "prev_patch = " << prev_patch << "\n"; + } + bool prev_flipped; + find_flap_vert(*tm.face(prev_tri), e, &prev_flipped); + int c = prev_flipped ? prev_patch.cell_below : prev_patch.cell_above; + if (dbg_level > 0) { + std::cout << "find_cell_for_point_near_edge returns " << c << "\n"; + } + return c; +} + +/** + * Find the ambient cell -- that is, the cell that is outside + * all other cells. + * If component_patches != nullptr, restrict consideration to patches + * in that vector. + * + * The method is to find an edge known to be on the convex hull + * of the mesh, then insert a dummy triangle that has that edge + * and a point known to be outside the whole mesh. Then sorting + * the triangles around the edge will reveal where the dummy triangle + * fits in that sorting order, and hence, the two adjacent patches + * to the dummy triangle - thus revealing the cell that the point + * known to be outside the whole mesh is in. + */ +static int find_ambient_cell(const IMesh &tm, + const Vector<int> *component_patches, + const TriMeshTopology &tmtopo, + const PatchesInfo &pinfo, + IMeshArena *arena) +{ + int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_AMBIENT_CELL\n"; + } + /* First find a vertex with the maximum x value. */ + /* Prefer not to populate the verts in the #IMesh just for this. */ + const Vert *v_extreme; + mpq_class extreme_x; + if (component_patches == nullptr) { + v_extreme = (*tm.face(0))[0]; + extreme_x = v_extreme->co_exact.x; + for (const Face *f : tm.faces()) { + for (const Vert *v : *f) { + const mpq_class &x = v->co_exact.x; + if (x > extreme_x) { + v_extreme = v; + extreme_x = x; + } + } + } + } + else { + if (dbg_level > 0) { + std::cout << "restrict to patches " << *component_patches << "\n"; + } + int p0 = (*component_patches)[0]; + v_extreme = (*tm.face(pinfo.patch(p0).tri(0)))[0]; + extreme_x = v_extreme->co_exact.x; + for (int p : *component_patches) { + for (int t : pinfo.patch(p).tris()) { + const Face *f = tm.face(t); + for (const Vert *v : *f) { + const mpq_class &x = v->co_exact.x; + if (x > extreme_x) { + v_extreme = v; + extreme_x = x; + } + } + } + } + } + if (dbg_level > 0) { + std::cout << "v_extreme = " << v_extreme << "\n"; + } + /* Find edge attached to v_extreme with max absolute slope + * when projected onto the XY plane. That edge is guaranteed to + * be on the convex hull of the mesh. */ + const Vector<Edge> &edges = tmtopo.vert_edges(v_extreme); + const mpq_class extreme_y = v_extreme->co_exact.y; + Edge ehull; + mpq_class max_abs_slope = -1; + for (Edge e : edges) { + const Vert *v_other = (e.v0() == v_extreme) ? e.v1() : e.v0(); + const mpq3 &co_other = v_other->co_exact; + mpq_class delta_x = co_other.x - extreme_x; + if (delta_x == 0) { + /* Vertical slope. */ + ehull = e; + break; + } + mpq_class abs_slope = abs((co_other.y - extreme_y) / delta_x); + if (abs_slope > max_abs_slope) { + ehull = e; + max_abs_slope = abs_slope; + } + } + if (dbg_level > 0) { + std::cout << "ehull = " << ehull << " slope = " << max_abs_slope << "\n"; + } + /* Sort triangles around ehull, including a dummy triangle that include a known point in ambient + * cell. */ + mpq3 p_in_ambient = v_extreme->co_exact; + p_in_ambient.x += 1; + int c_ambient = find_cell_for_point_near_edge(p_in_ambient, ehull, tm, tmtopo, pinfo, arena); + if (dbg_level > 0) { + std::cout << "FIND_AMBIENT_CELL returns " << c_ambient << "\n"; + } + return c_ambient; +} + +/** + * We need an edge on the convex hull of the edges incident on \a closestp + * in order to sort around, including a dummy triangle that has \a testp and + * the sorting edge vertices. So we don't want an edge that is co-linear + * with the line through \a testp and \a closestp. + * The method is to project onto a plane that contains `testp-closestp`, + * and then choose the edge that, when projected, has the maximum absolute + * slope (regarding the line `testp-closestp` as the x-axis for slope computation). + */ +static Edge find_good_sorting_edge(const Vert *testp, + const Vert *closestp, + const TriMeshTopology &tmtopo) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_GOOD_SORTING_EDGE testp = " << testp << ", closestp = " << closestp << "\n"; + } + /* We want to project the edges incident to closestp onto a plane + * whose ordinate direction will be regarded as going from closestp to testp, + * and whose abscissa direction is some perpendicular to that. + * A perpendicular direction can be found by swapping two coordinates + * and negating one, and zeroing out the third, being careful that one + * of the swapped vertices is non-zero. */ + const mpq3 &co_closest = closestp->co_exact; + const mpq3 &co_test = testp->co_exact; + BLI_assert(co_test != co_closest); + mpq3 abscissa = co_test - co_closest; + /* Find a non-zero-component axis of abscissa. */ + int axis; + for (axis = 0; axis < 3; ++axis) { + if (abscissa[axis] != 0) { + break; + } + } + BLI_assert(axis < 3); + int axis_next = (axis + 1) % 3; + int axis_next_next = (axis_next + 1) % 3; + mpq3 ordinate; + ordinate[axis] = abscissa[axis_next]; + ordinate[axis_next] = -abscissa[axis]; + ordinate[axis_next_next] = 0; + /* By construction, dot(abscissa, ordinate) == 0, so they are perpendicular. */ + mpq3 normal = mpq3::cross(abscissa, ordinate); + if (dbg_level > 0) { + std::cout << "abscissa = " << abscissa << "\n"; + std::cout << "ordinate = " << ordinate << "\n"; + std::cout << "normal = " << normal << "\n"; + } + mpq_class nlen2 = normal.length_squared(); + mpq_class max_abs_slope = -1; + Edge esort; + const Vector<Edge> &edges = tmtopo.vert_edges(closestp); + for (Edge e : edges) { + const Vert *v_other = (e.v0() == closestp) ? e.v1() : e.v0(); + const mpq3 &co_other = v_other->co_exact; + mpq3 evec = co_other - co_closest; + /* Get projection of evec onto plane of abscissa and ordinate. */ + mpq3 proj_evec = evec - (mpq3::dot(evec, normal) / nlen2) * normal; + /* The projection calculations along the abscissa and ordinate should + * be scaled by 1/abscissa and 1/ordinate respectively, + * but we can skip: it won't affect which `evec` has the maximum slope. */ + mpq_class evec_a = mpq3::dot(proj_evec, abscissa); + mpq_class evec_o = mpq3::dot(proj_evec, ordinate); + if (dbg_level > 0) { + std::cout << "e = " << e << "\n"; + std::cout << "v_other = " << v_other << "\n"; + std::cout << "evec = " << evec << ", proj_evec = " << proj_evec << "\n"; + std::cout << "evec_a = " << evec_a << ", evec_o = " << evec_o << "\n"; + } + if (evec_a == 0) { + /* evec is perpendicular to abscissa. */ + esort = e; + if (dbg_level > 0) { + std::cout << "perpendicular esort is " << esort << "\n"; + } + break; + } + mpq_class abs_slope = abs(evec_o / evec_a); + if (abs_slope > max_abs_slope) { + esort = e; + max_abs_slope = abs_slope; + if (dbg_level > 0) { + std::cout << "with abs_slope " << abs_slope << " new esort is " << esort << "\n"; + } + } + } + return esort; +} + +/** + * Find the cell that contains v. Consider the cells adjacent to triangle t. + * The close_edge and close_vert values are what were returned by + * #closest_on_tri_to_point when determining that v was close to t. + * They will indicate whether the point of closest approach to t is to + * an edge of t, a vertex of t, or somewhere inside t. + * + * The algorithm is similar to the one for find_ambient_cell, except that + * instead of an arbitrary point known to be outside the whole mesh, we + * have a particular point (v) and we just want to determine the patches + * that that point is between in sorting-around-an-edge order. + */ +static int find_containing_cell(const Vert *v, + int t, + int close_edge, + int close_vert, + const PatchesInfo &pinfo, + const IMesh &tm, + const TriMeshTopology &tmtopo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_CONTAINING_CELL v=" << v << ", t=" << t << "\n"; + } + const Face &tri = *tm.face(t); + Edge etest; + if (close_edge == -1 && close_vert == -1) { + /* Choose any edge if closest point is inside the triangle. */ + close_edge = 0; + } + if (close_edge != -1) { + const Vert *v0 = tri[close_edge]; + const Vert *v1 = tri[(close_edge + 1) % 3]; + const Vector<Edge> &edges = tmtopo.vert_edges(v0); + if (dbg_level > 0) { + std::cout << "look for edge containing " << v0 << " and " << v1 << "\n"; + std::cout << " in edges: "; + for (Edge e : edges) { + std::cout << e << " "; + } + std::cout << "\n"; + } + for (Edge e : edges) { + if ((e.v0() == v0 && e.v1() == v1) || (e.v0() == v1 && e.v1() == v0)) { + etest = e; + break; + } + } + } + else { + int cv = close_vert; + const Vert *vert_cv = tri[cv]; + if (vert_cv == v) { + /* Need to use another one to find sorting edge. */ + vert_cv = tri[(cv + 1) % 3]; + BLI_assert(vert_cv != v); + } + etest = find_good_sorting_edge(v, vert_cv, tmtopo); + } + BLI_assert(etest.v0() != nullptr); + if (dbg_level > 0) { + std::cout << "etest = " << etest << "\n"; + } + int c = find_cell_for_point_near_edge(v->co_exact, etest, tm, tmtopo, pinfo, arena); + if (dbg_level > 0) { + std::cout << "find_containing_cell returns " << c << "\n"; + } + return c; +} + +/** + * Find the closest point in triangle (a, b, c) to point p. + * Return the distance squared to that point. + * Also, if the closest point in the triangle is on a vertex, + * return 0, 1, or 2 for a, b, c in *r_vert; else -1. + * If the closest point is on an edge, return 0, 1, or 2 + * for edges ab, bc, or ca in *r_edge; else -1. + * (Adapted from #closest_on_tri_to_point_v3()). + */ +static mpq_class closest_on_tri_to_point( + const mpq3 &p, const mpq3 &a, const mpq3 &b, const mpq3 &c, int *r_edge, int *r_vert) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "CLOSEST_ON_TRI_TO_POINT p = " << p << "\n"; + std::cout << " a = " << a << ", b = " << b << ", c = " << c << "\n"; + } + /* Check if p in vertex region outside a. */ + mpq3 ab = b - a; + mpq3 ac = c - a; + mpq3 ap = p - a; + mpq_class d1 = mpq3::dot(ab, ap); + mpq_class d2 = mpq3::dot(ac, ap); + if (d1 <= 0 && d2 <= 0) { + /* Barycentric coordinates (1,0,0). */ + *r_edge = -1; + *r_vert = 0; + if (dbg_level > 0) { + std::cout << " answer = a\n"; + } + return mpq3::distance_squared(p, a); + } + /* Check if p in vertex region outside b. */ + mpq3 bp = p - b; + mpq_class d3 = mpq3::dot(ab, bp); + mpq_class d4 = mpq3::dot(ac, bp); + if (d3 >= 0 && d4 <= d3) { + /* Barycentric coordinates (0,1,0). */ + *r_edge = -1; + *r_vert = 1; + if (dbg_level > 0) { + std::cout << " answer = b\n"; + } + return mpq3::distance_squared(p, b); + } + /* Check if p in region of ab. */ + mpq_class vc = d1 * d4 - d3 * d2; + if (vc <= 0 && d1 >= 0 && d3 <= 0) { + mpq_class v = d1 / (d1 - d3); + /* Barycentric coordinates (1-v,v,0). */ + mpq3 r = a + v * ab; + *r_vert = -1; + *r_edge = 0; + if (dbg_level > 0) { + std::cout << " answer = on ab at " << r << "\n"; + } + return mpq3::distance_squared(p, r); + } + /* Check if p in vertex region outside c. */ + mpq3 cp = p - c; + mpq_class d5 = mpq3::dot(ab, cp); + mpq_class d6 = mpq3::dot(ac, cp); + if (d6 >= 0 && d5 <= d6) { + /* Barycentric coordinates (0,0,1). */ + *r_edge = -1; + *r_vert = 2; + if (dbg_level > 0) { + std::cout << " answer = c\n"; + } + return mpq3::distance_squared(p, c); + } + /* Check if p in edge region of ac. */ + mpq_class vb = d5 * d2 - d1 * d6; + if (vb <= 0 && d2 >= 0 && d6 <= 0) { + mpq_class w = d2 / (d2 - d6); + /* Barycentric coordinates (1-w,0,w). */ + mpq3 r = a + w * ac; + *r_vert = -1; + *r_edge = 2; + if (dbg_level > 0) { + std::cout << " answer = on ac at " << r << "\n"; + } + return mpq3::distance_squared(p, r); + } + /* Check if p in edge region of bc. */ + mpq_class va = d3 * d6 - d5 * d4; + if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) { + mpq_class w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); + /* Barycentric coordinates (0,1-w,w). */ + mpq3 r = c - b; + r = w * r; + r = r + b; + *r_vert = -1; + *r_edge = 1; + if (dbg_level > 0) { + std::cout << " answer = on bc at " << r << "\n"; + } + return mpq3::distance_squared(p, r); + } + /* p inside face region. Compute barycentric coordinates (u,v,w). */ + mpq_class denom = 1 / (va + vb + vc); + mpq_class v = vb * denom; + mpq_class w = vc * denom; + ac = w * ac; + mpq3 r = a + v * ab; + r = r + ac; + *r_vert = -1; + *r_edge = -1; + if (dbg_level > 0) { + std::cout << " answer = inside at " << r << "\n"; + } + return mpq3::distance_squared(p, r); +} + +struct ComponentContainer { + int containing_component{NO_INDEX}; + int nearest_cell{NO_INDEX}; + mpq_class dist_to_cell; + + ComponentContainer(int cc, int cell, mpq_class d) + : containing_component(cc), nearest_cell(cell), dist_to_cell(d) + { + } +}; + +/** + * Find out all the components, not equal to comp, that contain a point + * in comp in a non-ambient cell of those components. + * In other words, find the components that comp is nested inside + * (maybe not directly nested, which is why there can be more than one). + */ +static Vector<ComponentContainer> find_component_containers(int comp, + const Vector<Vector<int>> &components, + const Array<int> &ambient_cell, + const IMesh &tm, + const PatchesInfo &pinfo, + const TriMeshTopology &tmtopo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_COMPONENT_CONTAINERS for comp " << comp << "\n"; + } + Vector<ComponentContainer> ans; + int test_p = components[comp][0]; + int test_t = pinfo.patch(test_p).tri(0); + const Vert *test_v = tm.face(test_t)[0].vert[0]; + if (dbg_level > 0) { + std::cout << "test vertex in comp: " << test_v << "\n"; + } + for (int comp_other : components.index_range()) { + if (comp == comp_other) { + continue; + } + if (dbg_level > 0) { + std::cout << "comp_other = " << comp_other << "\n"; + } + int nearest_tri = NO_INDEX; + int nearest_tri_close_vert = -1; + int nearest_tri_close_edge = -1; + mpq_class nearest_tri_dist_squared; + for (int p : components[comp_other]) { + const Patch &patch = pinfo.patch(p); + for (int t : patch.tris()) { + const Face &tri = *tm.face(t); + if (dbg_level > 1) { + std::cout << "tri " << t << " = " << &tri << "\n"; + } + int close_vert; + int close_edge; + mpq_class d2 = closest_on_tri_to_point(test_v->co_exact, + tri[0]->co_exact, + tri[1]->co_exact, + tri[2]->co_exact, + &close_edge, + &close_vert); + if (dbg_level > 1) { + std::cout << " close_edge=" << close_edge << " close_vert=" << close_vert + << " dsquared=" << d2.get_d() << "\n"; + } + if (nearest_tri == NO_INDEX || d2 < nearest_tri_dist_squared) { + nearest_tri = t; + nearest_tri_close_edge = close_edge; + nearest_tri_close_vert = close_vert; + nearest_tri_dist_squared = d2; + } + } + } + if (dbg_level > 0) { + std::cout << "closest tri to comp=" << comp << " in comp_other=" << comp_other << " is t" + << nearest_tri << "\n"; + const Face *tn = tm.face(nearest_tri); + std::cout << "tri = " << tn << "\n"; + std::cout << " (" << tn->vert[0]->co << "," << tn->vert[1]->co << "," << tn->vert[2]->co + << ")\n"; + } + int containing_cell = find_containing_cell(test_v, + nearest_tri, + nearest_tri_close_edge, + nearest_tri_close_vert, + + pinfo, + tm, + tmtopo, + arena); + if (dbg_level > 0) { + std::cout << "containing cell = " << containing_cell << "\n"; + } + if (containing_cell != ambient_cell[comp_other]) { + ans.append(ComponentContainer(comp_other, containing_cell, nearest_tri_dist_squared)); + } + } + return ans; +} + +/** + * The cells and patches are supposed to form a bipartite graph. + * The graph may be disconnected (if parts of meshes are nested or side-by-side + * without intersection with other each other). + * Connect the bipartite graph. This involves discovering the connected components + * of the patches, then the nesting structure of those components. + */ +static void finish_patch_cell_graph(const IMesh &tm, + CellsInfo &cinfo, + PatchesInfo &pinfo, + const TriMeshTopology &tmtopo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FINISH_PATCH_CELL_GRAPH\n"; + } + Vector<Vector<int>> components = find_patch_components(cinfo, pinfo); + if (components.size() <= 1) { + if (dbg_level > 0) { + std::cout << "one component so finish_patch_cell_graph does no work\n"; + } + return; + } + if (dbg_level > 0) { + std::cout << "components:\n"; + for (int comp : components.index_range()) { + std::cout << comp << ": " << components[comp] << "\n"; + } + } + Array<int> ambient_cell(components.size()); + for (int comp : components.index_range()) { + ambient_cell[comp] = find_ambient_cell(tm, &components[comp], tmtopo, pinfo, arena); + } + if (dbg_level > 0) { + std::cout << "ambient cells:\n"; + for (int comp : ambient_cell.index_range()) { + std::cout << comp << ": " << ambient_cell[comp] << "\n"; + } + } + int tot_components = components.size(); + Array<Vector<ComponentContainer>> comp_cont(tot_components); + for (int comp : components.index_range()) { + comp_cont[comp] = find_component_containers( + comp, components, ambient_cell, tm, pinfo, tmtopo, arena); + } + if (dbg_level > 0) { + std::cout << "component containers:\n"; + for (int comp : comp_cont.index_range()) { + std::cout << comp << ": "; + for (const ComponentContainer &cc : comp_cont[comp]) { + std::cout << "[containing_comp=" << cc.containing_component + << ", nearest_cell=" << cc.nearest_cell << ", d2=" << cc.dist_to_cell << "] "; + } + std::cout << "\n"; + } + } + if (dbg_level > 1) { + write_obj_cell_patch(tm, cinfo, pinfo, false, "beforemerge"); + } + /* For nested components, merge their ambient cell with the nearest containing cell. */ + Vector<int> outer_components; + for (int comp : comp_cont.index_range()) { + if (comp_cont[comp].size() == 0) { + outer_components.append(comp); + } + else { + ComponentContainer &closest = comp_cont[comp][0]; + for (int i = 1; i < comp_cont[comp].size(); ++i) { + if (comp_cont[comp][i].dist_to_cell < closest.dist_to_cell) { + closest = comp_cont[comp][i]; + } + } + int comp_ambient = ambient_cell[comp]; + int cont_cell = closest.nearest_cell; + if (dbg_level > 0) { + std::cout << "merge comp " << comp << "'s ambient cell=" << comp_ambient << " to cell " + << cont_cell << "\n"; + } + merge_cells(cont_cell, comp_ambient, cinfo, pinfo); + } + } + /* For outer components (not nested in any other component), merge their ambient cells. */ + if (outer_components.size() > 1) { + int merged_ambient = ambient_cell[outer_components[0]]; + for (int i = 1; i < outer_components.size(); ++i) { + if (dbg_level > 0) { + std::cout << "merge comp " << outer_components[i] + << "'s ambient cell=" << ambient_cell[outer_components[i]] << " to cell " + << merged_ambient << "\n"; + } + merge_cells(merged_ambient, ambient_cell[outer_components[i]], cinfo, pinfo); + } + } + if (dbg_level > 0) { + std::cout << "after FINISH_PATCH_CELL_GRAPH\nCells\n"; + for (int i : cinfo.index_range()) { + if (cinfo.cell(i).merged_to() == NO_INDEX) { + std::cout << i << ": " << cinfo.cell(i) << "\n"; + } + } + std::cout << "Patches\n"; + for (int i : pinfo.index_range()) { + std::cout << i << ": " << pinfo.patch(i) << "\n"; + } + if (dbg_level > 1) { + write_obj_cell_patch(tm, cinfo, pinfo, false, "finished"); + } + } +} + +/** + * Starting with ambient cell c_ambient, with all zeros for winding numbers, + * propagate winding numbers to all the other cells. + * There will be a vector of \a nshapes winding numbers in each cell, one per + * input shape. + * As one crosses a patch into a new cell, the original shape (mesh part) + * that that patch was part of dictates which winding number changes. + * The shape_fn(triangle_number) function should return the shape that the + * triangle is part of. + * Also, as soon as the winding numbers for a cell are set, use bool_optype + * to decide whether that cell is included or excluded from the boolean output. + * If included, the cell's in_output_volume will be set to true. + */ +static void propagate_windings_and_in_output_volume(PatchesInfo &pinfo, + CellsInfo &cinfo, + int c_ambient, + BoolOpType op, + int nshapes, + std::function<int(int)> shape_fn) +{ + int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "PROPAGATE_WINDINGS, ambient cell = " << c_ambient << "\n"; + } + Cell &cell_ambient = cinfo.cell(c_ambient); + cell_ambient.seed_ambient_winding(); + /* Use a vector as a queue. It can't grow bigger than number of cells. */ + Vector<int> queue; + queue.reserve(cinfo.tot_cell()); + int queue_head = 0; + queue.append(c_ambient); + while (queue_head < queue.size()) { + int c = queue[queue_head++]; + if (dbg_level > 1) { + std::cout << "process cell " << c << "\n"; + } + Cell &cell = cinfo.cell(c); + for (int p : cell.patches()) { + Patch &patch = pinfo.patch(p); + bool p_above_c = patch.cell_below == c; + int c_neighbor = p_above_c ? patch.cell_above : patch.cell_below; + if (dbg_level > 1) { + std::cout << " patch " << p << " p_above_c = " << p_above_c << "\n"; + std::cout << " c_neighbor = " << c_neighbor << "\n"; + } + Cell &cell_neighbor = cinfo.cell(c_neighbor); + if (!cell_neighbor.winding_assigned()) { + int winding_delta = p_above_c ? -1 : 1; + int t = patch.tri(0); + int shape = shape_fn(t); + BLI_assert(shape < nshapes); + UNUSED_VARS_NDEBUG(nshapes); + if (dbg_level > 1) { + std::cout << " representative tri " << t << ": in shape " << shape << "\n"; + } + cell_neighbor.set_winding_and_in_output_volume(cell, shape, winding_delta, op); + if (dbg_level > 1) { + std::cout << " now cell_neighbor = " << cell_neighbor << "\n"; + } + queue.append(c_neighbor); + BLI_assert(queue.size() <= cinfo.tot_cell()); + } + } + } + if (dbg_level > 0) { + std::cout << "\nPROPAGATE_WINDINGS result\n"; + for (int i = 0; i < cinfo.tot_cell(); ++i) { + std::cout << i << ": " << cinfo.cell(i) << "\n"; + } + } +} + +/** + * Given an array of winding numbers, where the ith entry is a cell's winding + * number with respect to input shape (mesh part) i, return true if the + * cell should be included in the output of the boolean operation. + * Intersection: all the winding numbers must be nonzero. + * Union: at least one winding number must be nonzero. + * Difference (first shape minus the rest): first winding number must be nonzero + * and the rest must have at least one zero winding number. + */ +static bool apply_bool_op(BoolOpType bool_optype, const Array<int> &winding) +{ + int nw = winding.size(); + BLI_assert(nw > 0); + switch (bool_optype) { + case BoolOpType::Intersect: { + for (int i = 0; i < nw; ++i) { + if (winding[i] == 0) { + return false; + } + } + return true; + } + case BoolOpType::Union: { + for (int i = 0; i < nw; ++i) { + if (winding[i] != 0) { + return true; + } + } + return false; + } + case BoolOpType::Difference: { + /* if nw > 2, make it shape 0 minus the union of the rest. */ + if (winding[0] == 0) { + return false; + } + if (nw == 1) { + return true; + } + for (int i = 1; i < nw; ++i) { + if (winding[i] == 0) { + return true; + } + } + return false; + } + default: + return false; + } +} + +/** + * Special processing for extract_from_in_output_volume_diffs to handle + * triangles that are part of stacks of geometrically identical + * triangles enclosing zero volume cells. + */ +static void extract_zero_volume_cell_tris(Vector<Face *> &r_tris, + const IMesh &tm_subdivided, + const PatchesInfo &pinfo, + const CellsInfo &cinfo, + IMeshArena *arena) +{ + int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "extract_zero_volume_cell_tris\n"; + } + /* Find patches that are adjacent to zero-volume cells. */ + Array<bool> adj_to_zv(pinfo.tot_patch()); + for (int p : pinfo.index_range()) { + const Patch &patch = pinfo.patch(p); + if (cinfo.cell(patch.cell_above).zero_volume() || cinfo.cell(patch.cell_below).zero_volume()) { + adj_to_zv[p] = true; + } + else { + adj_to_zv[p] = false; + } + } + /* Partition the adj_to_zv patches into stacks. */ + Vector<Vector<int>> patch_stacks; + Array<bool> allocated_to_stack(pinfo.tot_patch(), false); + for (int p : pinfo.index_range()) { + if (!adj_to_zv[p] || allocated_to_stack[p]) { + continue; + } + int stack_index = patch_stacks.size(); + patch_stacks.append(Vector<int>{p}); + Vector<int> &stack = patch_stacks[stack_index]; + Vector<bool> flipped{false}; + allocated_to_stack[p] = true; + /* We arbitrarily choose p's above and below directions as above and below for whole stack. + * Triangles in the stack that don't follow that convention are marked with flipped = true. + * The non-zero-volume cell above the whole stack, following this convention, is + * above_stack_cell. The non-zero-volume cell below the whole stack is #below_stack_cell. */ + /* First, walk in the above_cell direction from p. */ + int pwalk = p; + const Patch *pwalk_patch = &pinfo.patch(pwalk); + int c = pwalk_patch->cell_above; + const Cell *cell = &cinfo.cell(c); + while (cell->zero_volume()) { + /* In zero-volume cells, the cell should have exactly two patches. */ + BLI_assert(cell->patches().size() == 2); + int pother = cell->patches()[0] == pwalk ? cell->patches()[1] : cell->patches()[0]; + bool flip = pinfo.patch(pother).cell_above == c; + flipped.append(flip); + stack.append(pother); + allocated_to_stack[pother] = true; + pwalk = pother; + pwalk_patch = &pinfo.patch(pwalk); + c = flip ? pwalk_patch->cell_below : pwalk_patch->cell_above; + cell = &cinfo.cell(c); + } + const Cell *above_stack_cell = cell; + /* Now walk in the below_cell direction from p. */ + pwalk = p; + pwalk_patch = &pinfo.patch(pwalk); + c = pwalk_patch->cell_below; + cell = &cinfo.cell(c); + while (cell->zero_volume()) { + BLI_assert(cell->patches().size() == 2); + int pother = cell->patches()[0] == pwalk ? cell->patches()[1] : cell->patches()[0]; + bool flip = pinfo.patch(pother).cell_below == c; + flipped.append(flip); + stack.append(pother); + allocated_to_stack[pother] = true; + pwalk = pother; + pwalk_patch = &pinfo.patch(pwalk); + c = flip ? pwalk_patch->cell_above : pwalk_patch->cell_below; + cell = &cinfo.cell(c); + } + const Cell *below_stack_cell = cell; + if (dbg_level > 0) { + std::cout << "process zero-volume patch stack " << stack << "\n"; + std::cout << "flipped = "; + for (bool b : flipped) { + std::cout << b << " "; + } + std::cout << "\n"; + } + if (above_stack_cell->in_output_volume() ^ below_stack_cell->in_output_volume()) { + bool need_flipped_tri = above_stack_cell->in_output_volume(); + if (dbg_level > 0) { + std::cout << "need tri " << (need_flipped_tri ? "flipped" : "") << "\n"; + } + int t_to_add = NO_INDEX; + for (int i : stack.index_range()) { + if (flipped[i] == need_flipped_tri) { + t_to_add = pinfo.patch(stack[i]).tri(0); + if (dbg_level > 0) { + std::cout << "using tri " << t_to_add << "\n"; + } + r_tris.append(tm_subdivided.face(t_to_add)); + break; + } + } + if (t_to_add == NO_INDEX) { + const Face *f = tm_subdivided.face(pinfo.patch(p).tri(0)); + const Face &tri = *f; + /* We need flipped version or else we would have found it above. */ + std::array<const Vert *, 3> flipped_vs = {tri[0], tri[2], tri[1]}; + std::array<int, 3> flipped_e_origs = { + tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]}; + std::array<bool, 3> flipped_is_intersect = { + tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]}; + Face *flipped_f = arena->add_face( + flipped_vs, f->orig, flipped_e_origs, flipped_is_intersect); + r_tris.append(flipped_f); + } + } + } +} + +/** + * Extract the output mesh from tm_subdivided and return it as a new mesh. + * The cells in \a cinfo must have cells-to-be-retained with in_output_volume set. + * We keep only triangles between those in the output volume and those not in. + * We flip the normals of any triangle that has an in_output_volume cell above + * and a not-in_output_volume cell below. + * For all stacks of exact duplicate co-planar triangles, we want to + * include either one version of the triangle or none, depending on + * whether the in_output_volume in_output_volumes on either side of the stack are + * different or the same. + */ +static IMesh extract_from_in_output_volume_diffs(const IMesh &tm_subdivided, + const PatchesInfo &pinfo, + const CellsInfo &cinfo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nEXTRACT_FROM_FLAG_DIFFS\n"; + } + Vector<Face *> out_tris; + out_tris.reserve(tm_subdivided.face_size()); + bool any_zero_volume_cell = false; + for (int t : tm_subdivided.face_index_range()) { + int p = pinfo.tri_patch(t); + const Patch &patch = pinfo.patch(p); + const Cell &cell_above = cinfo.cell(patch.cell_above); + const Cell &cell_below = cinfo.cell(patch.cell_below); + if (dbg_level > 0) { + std::cout << "tri " << t << ": cell_above=" << patch.cell_above + << " cell_below=" << patch.cell_below << "\n"; + std::cout << " in_output_volume_above=" << cell_above.in_output_volume() + << " in_output_volume_below=" << cell_below.in_output_volume() << "\n"; + } + bool adjacent_zero_volume_cell = cell_above.zero_volume() || cell_below.zero_volume(); + any_zero_volume_cell |= adjacent_zero_volume_cell; + if (cell_above.in_output_volume() ^ cell_below.in_output_volume() && + !adjacent_zero_volume_cell) { + bool flip = cell_above.in_output_volume(); + if (dbg_level > 0) { + std::cout << "need tri " << t << " flip=" << flip << "\n"; + } + Face *f = tm_subdivided.face(t); + if (flip) { + Face &tri = *f; + std::array<const Vert *, 3> flipped_vs = {tri[0], tri[2], tri[1]}; + std::array<int, 3> flipped_e_origs = { + tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]}; + std::array<bool, 3> flipped_is_intersect = { + tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]}; + Face *flipped_f = arena->add_face( + flipped_vs, f->orig, flipped_e_origs, flipped_is_intersect); + out_tris.append(flipped_f); + } + else { + out_tris.append(f); + } + } + } + if (any_zero_volume_cell) { + extract_zero_volume_cell_tris(out_tris, tm_subdivided, pinfo, cinfo, arena); + } + return IMesh(out_tris); +} + +static const char *bool_optype_name(BoolOpType op) +{ + switch (op) { + case BoolOpType::None: + return "none"; + case BoolOpType::Intersect: + return "intersect"; + case BoolOpType::Union: + return "union"; + case BoolOpType::Difference: + return "difference"; + default: + return "<unknown>"; + } +} + +static mpq3 calc_point_inside_tri(const Face &tri) +{ + const Vert *v0 = tri.vert[0]; + const Vert *v1 = tri.vert[1]; + const Vert *v2 = tri.vert[2]; + mpq3 ans = v0->co_exact / 3 + v1->co_exact / 3 + v2->co_exact / 3; + return ans; +} + +/** + * Return the Generalized Winding Number of point \a testp with respect to the + * volume implied by the faces for which shape_fn returns the value shape. + * See "Robust Inside-Outside Segmentation using Generalized Winding Numbers" + * by Jacobson, Kavan, and Sorkine-Hornung. + * This is like a winding number in that if it is positive, the point + * is inside the volume. But it is tolerant of not-completely-watertight + * volumes, still doing a passable job of classifying inside/outside + * as we intuitively understand that to mean. + * + * TOOD: speed up this calculation using the hierarchical algorithm in that paper. + */ +static double generalized_winding_number(const IMesh &tm, + std::function<int(int)> shape_fn, + const double3 &testp, + int shape) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "GENERALIZED_WINDING_NUMBER testp = " << testp << ", shape = " << shape << "\n"; + } + double gwn = 0; + for (int t : tm.face_index_range()) { + const Face *f = tm.face(t); + const Face &tri = *f; + if (shape_fn(tri.orig) == shape) { + if (dbg_level > 0) { + std::cout << "accumulate for tri t = " << t << " = " << f << "\n"; + } + const Vert *v0 = tri.vert[0]; + const Vert *v1 = tri.vert[1]; + const Vert *v2 = tri.vert[2]; + double3 a = v0->co - testp; + double3 b = v1->co - testp; + double3 c = v2->co - testp; + /* Calculate the solid angle of abc relative to origin. + * See "The Solid Angle of a Plane Triangle" by Oosterom and Strackee + * for the derivation of the formula. */ + double alen = a.length(); + double blen = b.length(); + double clen = c.length(); + double3 bxc = double3::cross_high_precision(b, c); + double num = double3::dot(a, bxc); + double denom = alen * blen * clen + double3::dot(a, b) * clen + double3::dot(a, c) * blen + + double3::dot(b, c) * alen; + if (denom == 0.0) { + if (dbg_level > 0) { + std::cout << "denom == 0, skipping this tri\n"; + } + continue; + } + double x = atan2(num, denom); + double fgwn = 2.0 * x; + if (dbg_level > 0) { + std::cout << "tri contributes " << fgwn << "\n"; + } + gwn += fgwn; + } + } + gwn = gwn / (M_PI * 4.0); + if (dbg_level > 0) { + std::cout << "final gwn = " << gwn << "\n"; + } + return gwn; +} + +/** + * Return true if point \a testp is inside the volume implied by the + * faces for which the shape_fn returns the value shape. + */ +static bool point_is_inside_shape(const IMesh &tm, + std::function<int(int)> shape_fn, + const double3 &testp, + int shape) +{ + double gwn = generalized_winding_number(tm, shape_fn, testp, shape); + /* Due to floating point error, an outside point should get a value + * of zero for gwn, but may have a very slightly positive value instead. + * It is not important to get this epsilon very small, because practical + * cases of interest will have gwn at least 0.2 if it is not zero. */ + return (gwn > 0.01); +} + +/** + * Use the Generalized Winding Number method for deciding if a patch of the + * mesh is supposed to be included or excluded in the boolean result, + * and return the mesh that is the boolean result. + */ +static IMesh gwn_boolean(const IMesh &tm, + BoolOpType op, + int nshapes, + std::function<int(int)> shape_fn, + const PatchesInfo &pinfo, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "GWN_BOOLEAN\n"; + } + IMesh ans; + Vector<Face *> out_faces; + out_faces.reserve(tm.face_size()); + BLI_assert(nshapes == 2); /* TODO: generalize. */ + UNUSED_VARS_NDEBUG(nshapes); + for (int p : pinfo.index_range()) { + const Patch &patch = pinfo.patch(p); + /* For test triangle, choose one in the middle of patch list + * as the ones near the beginning may be very near other patches. */ + int test_t_index = patch.tri(patch.tot_tri() / 2); + Face &tri_test = *tm.face(test_t_index); + /* Assume all triangles in a patch are in the same shape. */ + int shape = shape_fn(tri_test.orig); + if (dbg_level > 0) { + std::cout << "process patch " << p << " = " << patch << "\n"; + std::cout << "test tri = " << test_t_index << " = " << &tri_test << "\n"; + std::cout << "shape = " << shape << "\n"; + } + if (shape == -1) { + continue; + } + mpq3 test_point = calc_point_inside_tri(tri_test); + double3 test_point_db(test_point[0].get_d(), test_point[1].get_d(), test_point[2].get_d()); + if (dbg_level > 0) { + std::cout << "test point = " << test_point_db << "\n"; + } + int other_shape = 1 - shape; + bool inside = point_is_inside_shape(tm, shape_fn, test_point_db, other_shape); + if (dbg_level > 0) { + std::cout << "test point is " << (inside ? "inside\n" : "outside\n"); + } + bool do_remove; + bool do_flip; + switch (op) { + case BoolOpType::Intersect: + do_remove = !inside; + do_flip = false; + break; + case BoolOpType::Union: + do_remove = inside; + do_flip = false; + break; + case BoolOpType::Difference: + do_remove = (shape == 0) ? inside : !inside; + do_flip = (shape == 1); + break; + default: + do_remove = false; + do_flip = false; + BLI_assert(false); + } + if (dbg_level > 0) { + std::cout << "result for patch " << p << ": remove=" << do_remove << ", flip=" << do_flip + << "\n"; + } + if (!do_remove) { + for (int t : patch.tris()) { + Face *f = tm.face(t); + if (!do_flip) { + out_faces.append(f); + } + else { + Face &tri = *f; + /* We need flipped version of f. */ + Array<const Vert *> flipped_vs = {tri[0], tri[2], tri[1]}; + Array<int> flipped_e_origs = {tri.edge_orig[2], tri.edge_orig[1], tri.edge_orig[0]}; + Array<bool> flipped_is_intersect = { + tri.is_intersect[2], tri.is_intersect[1], tri.is_intersect[0]}; + Face *flipped_f = arena->add_face( + flipped_vs, f->orig, flipped_e_origs, flipped_is_intersect); + out_faces.append(flipped_f); + } + } + } + } + ans.set_faces(out_faces); + 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<mpq_class> &cdt_out, int v1, int v2) +{ + for (int e : cdt_out.edge.index_range()) { + const std::pair<int, int> &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. + */ +static Array<Face *> triangulate_poly(Face *f, IMeshArena *arena) +{ + int flen = f->size(); + CDT_input<mpq_class> cdt_in; + cdt_in.vert = Array<mpq2>(flen); + cdt_in.face = Array<Vector<int>>(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<mpq_class> cdt_out = delaunay_2d_calc(cdt_in, CDT_INSIDE); + int n_tris = cdt_out.face.size(); + Array<Face *> ans(n_tris); + for (int t = 0; t < n_tris; ++t) { + int i_v_out[3]; + const Vert *v[3]; + int eo[3]; + for (int i = 0; i < 3; ++i) { + i_v_out[i] = cdt_out.face[t][i]; + v[i] = (*f)[cdt_out.vert_orig[i_v_out[i]][0]]; + } + 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}); + } + } + 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 *> 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<Face *> 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). + */ +static std::pair<int, int> find_tris_common_edge(const Face &tri1, const Face &tri2) +{ + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + if (tri1[(i + 1) % 3] == tri2[j] && tri1[i] == tri2[(j + 1) % 3]) { + return std::pair<int, int>(i, j); + } + } + } + return std::pair<int, int>(-1, -1); +} + +struct MergeEdge { + /** Length (squared) of the edge, used for sorting. */ + double len_squared = 0.0; + /* v1 and v2 are the ends of the edge, ordered so that `v1->id < v2->id` */ + const Vert *v1 = nullptr; + const Vert *v2 = nullptr; + /* left_face and right_face are indices into #FaceMergeState.face. */ + int left_face = -1; + int right_face = -1; + int orig = -1; /* An edge orig index that can be used for this edge. */ + /** Is it allowed to dissolve this edge? */ + bool dissolvable = false; + /** Is this an intersect edge? */ + bool is_intersect = false; + + MergeEdge() = default; + + MergeEdge(const Vert *va, const Vert *vb) + { + if (va->id < vb->id) { + this->v1 = va; + this->v2 = vb; + } + else { + this->v1 = vb; + this->v2 = va; + } + }; +}; + +struct MergeFace { + /** The current sequence of Verts forming this face. */ + Vector<const Vert *> vert; + /** For each position in face, what is index in #FaceMergeState of edge for that position? */ + Vector<int> edge; + /** If not -1, merge_to gives a face index in #FaceMergeState that this is merged to. */ + int merge_to = -1; + /** A face->orig that can be used for the merged face. */ + int orig = -1; +}; +struct FaceMergeState { + /** + * The faces being considered for merging. Some will already have been merge (merge_to != -1). + */ + Vector<MergeFace> face; + /** + * The edges that are part of the faces in face[], together with current topological + * information (their left and right faces) and whether or not they are dissolvable. + */ + Vector<MergeEdge> edge; + /** + * `edge_map` maps a pair of `const Vert *` ids (in canonical order: smaller id first) + * to the index in the above edge vector in which to find the corresponding #MergeEdge. + */ + Map<std::pair<int, int>, int> edge_map; +}; + +static std::ostream &operator<<(std::ostream &os, const FaceMergeState &fms) +{ + os << "faces:\n"; + for (int f : fms.face.index_range()) { + const MergeFace &mf = fms.face[f]; + std::cout << f << ": orig=" << mf.orig << " verts "; + for (const Vert *v : mf.vert) { + std::cout << v << " "; + } + std::cout << "\n"; + std::cout << " edges " << mf.edge << "\n"; + std::cout << " merge_to = " << mf.merge_to << "\n"; + } + os << "\nedges:\n"; + for (int e : fms.edge.index_range()) { + const MergeEdge &me = fms.edge[e]; + std::cout << e << ": (" << me.v1 << "," << me.v2 << ") left=" << me.left_face + << " right=" << me.right_face << " dis=" << me.dissolvable << " orig=" << me.orig + << " is_int=" << me.is_intersect << "\n"; + } + return os; +} + +/** + * \a tris all have the same original face. + * Find the 2d edge/triangle topology for these triangles, but only the ones facing in the + * norm direction, and whether each edge is dissolvable or not. + */ +static void init_face_merge_state(FaceMergeState *fms, + const Vector<int> &tris, + const IMesh &tm, + const double3 &norm) +{ + constexpr int dbg_level = 0; + /* Reserve enough faces and edges so that neither will have to resize. */ + fms->face.reserve(tris.size() + 1); + fms->edge.reserve(3 * tris.size()); + fms->edge_map.reserve(3 * tris.size()); + if (dbg_level > 0) { + std::cout << "\nINIT_FACE_MERGE_STATE\n"; + } + for (int t : tris.index_range()) { + MergeFace mf; + const Face &tri = *tm.face(tris[t]); + if (dbg_level > 0) { + std::cout << "process tri = " << &tri << "\n"; + } + BLI_assert(tri.plane_populated()); + if (double3::dot(norm, tri.plane->norm) <= 0.0) { + if (dbg_level > 0) { + std::cout << "triangle has wrong orientation, skipping\n"; + } + continue; + } + mf.vert.append(tri[0]); + mf.vert.append(tri[1]); + mf.vert.append(tri[2]); + mf.orig = tri.orig; + int f = fms->face.append_and_get_index(mf); + if (dbg_level > 1) { + std::cout << "appended MergeFace for tri at f = " << f << "\n"; + } + for (int i = 0; i < 3; ++i) { + int inext = (i + 1) % 3; + MergeEdge new_me(mf.vert[i], mf.vert[inext]); + std::pair<int, int> canon_vs(new_me.v1->id, new_me.v2->id); + int me_index = fms->edge_map.lookup_default(canon_vs, -1); + if (dbg_level > 1) { + std::cout << "new_me = canon_vs = " << new_me.v1 << ", " << new_me.v2 << "\n"; + std::cout << "me_index lookup = " << me_index << "\n"; + } + if (me_index == -1) { + double3 vec = new_me.v2->co - new_me.v1->co; + new_me.len_squared = vec.length_squared(); + new_me.orig = tri.edge_orig[i]; + new_me.is_intersect = tri.is_intersect[i]; + new_me.dissolvable = (new_me.orig == NO_INDEX && !new_me.is_intersect); + fms->edge.append(new_me); + me_index = fms->edge.size() - 1; + fms->edge_map.add_new(canon_vs, me_index); + if (dbg_level > 1) { + std::cout << "added new me with me_index = " << me_index << "\n"; + std::cout << " len_squared = " << new_me.len_squared << " orig = " << new_me.orig + << ", is_intersect" << new_me.is_intersect + << ", dissolvable = " << new_me.dissolvable << "\n"; + } + } + MergeEdge &me = fms->edge[me_index]; + if (dbg_level > 1) { + std::cout << "retrieved me at index " << me_index << ":\n"; + std::cout << " v1 = " << me.v1 << " v2 = " << me.v2 << "\n"; + std::cout << " dis = " << me.dissolvable << " int = " << me.is_intersect << "\n"; + std::cout << " left_face = " << me.left_face << " right_face = " << me.right_face << "\n"; + } + if (me.dissolvable && tri.edge_orig[i] != NO_INDEX) { + if (dbg_level > 1) { + std::cout << "reassigning orig to " << tri.edge_orig[i] << ", dissolvable = false\n"; + } + me.dissolvable = false; + me.orig = tri.edge_orig[i]; + } + if (me.dissolvable && tri.is_intersect[i]) { + if (dbg_level > 1) { + std::cout << "reassigning dissolvable = false, is_intersect = true\n"; + } + me.dissolvable = false; + me.is_intersect = true; + } + /* This face is left or right depending on orientation of edge. */ + if (me.v1 == mf.vert[i]) { + if (dbg_level > 1) { + std::cout << "me.v1 == mf.vert[i] so set edge[" << me_index << "].left_face = " << f + << "\n"; + } + BLI_assert(me.left_face == -1); + fms->edge[me_index].left_face = f; + } + else { + if (dbg_level > 1) { + std::cout << "me.v1 != mf.vert[i] so set edge[" << me_index << "].right_face = " << f + << "\n"; + } + BLI_assert(me.right_face == -1); + fms->edge[me_index].right_face = f; + } + fms->face[f].edge.append(me_index); + } + } + if (dbg_level > 0) { + std::cout << *fms; + } +} + +/** + * To have a valid #BMesh, there are constraints on what edges can be removed. + * We cannot remove an edge if (a) it would create two disconnected boundary parts + * (which will happen if there's another edge sharing the same two faces); + * or (b) it would create a face with a repeated vertex. + */ +static bool dissolve_leaves_valid_bmesh(FaceMergeState *fms, + const MergeEdge &me, + int me_index, + const MergeFace &mf_left, + const MergeFace &mf_right) +{ + int a_edge_start = mf_left.edge.first_index_of_try(me_index); + BLI_assert(a_edge_start != -1); + int alen = mf_left.vert.size(); + int blen = mf_right.vert.size(); + int b_left_face = me.right_face; + bool ok = true; + /* Is there another edge, not me, in A's face, whose right face is B's left? */ + for (int a_e_index = (a_edge_start + 1) % alen; ok && a_e_index != a_edge_start; + a_e_index = (a_e_index + 1) % alen) { + const MergeEdge &a_me_cur = fms->edge[mf_left.edge[a_e_index]]; + if (a_me_cur.right_face == b_left_face) { + ok = false; + } + } + /* Is there a vert in A, not me.v1 or me.v2, that is also in B? + * One could avoid this O(n^2) algorithm if had a structure + * saying which faces a vertex touches. */ + for (int a_v_index = 0; ok && a_v_index < alen; ++a_v_index) { + const Vert *a_v = mf_left.vert[a_v_index]; + if (a_v != me.v1 && a_v != me.v2) { + for (int b_v_index = 0; b_v_index < blen; ++b_v_index) { + const Vert *b_v = mf_right.vert[b_v_index]; + if (a_v == b_v) { + ok = false; + } + } + } + } + return ok; +} + +/** + * mf_left and mf_right should share a #MergeEdge me, having index me_index. + * We change mf_left to remove edge me and insert the appropriate edges of + * mf_right in between the start and end vertices of that edge. + * We change the left face of the spliced-in edges to be mf_left's index. + * We mark the merge_to property of mf_right, which is now in essence deleted. + */ +static void splice_faces( + FaceMergeState *fms, MergeEdge &me, int me_index, MergeFace &mf_left, MergeFace &mf_right) +{ + int a_edge_start = mf_left.edge.first_index_of_try(me_index); + int b_edge_start = mf_right.edge.first_index_of_try(me_index); + BLI_assert(a_edge_start != -1 && b_edge_start != -1); + int alen = mf_left.vert.size(); + int blen = mf_right.vert.size(); + Vector<const Vert *> splice_vert; + Vector<int> splice_edge; + splice_vert.reserve(alen + blen - 2); + splice_edge.reserve(alen + blen - 2); + int ai = 0; + while (ai < a_edge_start) { + splice_vert.append(mf_left.vert[ai]); + splice_edge.append(mf_left.edge[ai]); + ++ai; + } + int bi = b_edge_start + 1; + while (bi != b_edge_start) { + if (bi >= blen) { + bi = 0; + if (bi == b_edge_start) { + break; + } + } + splice_vert.append(mf_right.vert[bi]); + splice_edge.append(mf_right.edge[bi]); + if (mf_right.vert[bi] == fms->edge[mf_right.edge[bi]].v1) { + fms->edge[mf_right.edge[bi]].left_face = me.left_face; + } + else { + fms->edge[mf_right.edge[bi]].right_face = me.left_face; + } + ++bi; + } + ai = a_edge_start + 1; + while (ai < alen) { + splice_vert.append(mf_left.vert[ai]); + splice_edge.append(mf_left.edge[ai]); + ++ai; + } + mf_right.merge_to = me.left_face; + mf_left.vert = splice_vert; + mf_left.edge = splice_edge; + me.left_face = -1; + me.right_face = -1; +} + +/** + * Given that fms has been properly initialized to contain a set of faces that + * together form a face or part of a face of the original #IMesh, and that + * it has properly recorded with faces are dissolvable, dissolve as many edges as possible. + * We try to dissolve in decreasing order of edge length, so that it is more likely + * that the final output doesn't have awkward looking long edges with extreme angles. + */ +static void do_dissolve(FaceMergeState *fms) +{ + const int dbg_level = 0; + if (dbg_level > 1) { + std::cout << "\nDO_DISSOLVE\n"; + } + Vector<int> dissolve_edges; + for (int e : fms->edge.index_range()) { + if (fms->edge[e].dissolvable) { + dissolve_edges.append(e); + } + } + if (dissolve_edges.size() == 0) { + return; + } + /* Things look nicer if we dissolve the longer edges first. */ + std::sort( + dissolve_edges.begin(), dissolve_edges.end(), [fms](const int &a, const int &b) -> bool { + return (fms->edge[a].len_squared > fms->edge[b].len_squared); + }); + if (dbg_level > 0) { + std::cout << "Sorted dissolvable edges: " << dissolve_edges << "\n"; + } + for (int me_index : dissolve_edges) { + MergeEdge &me = fms->edge[me_index]; + if (me.left_face == -1 || me.right_face == -1) { + continue; + } + MergeFace &mf_left = fms->face[me.left_face]; + MergeFace &mf_right = fms->face[me.right_face]; + if (!dissolve_leaves_valid_bmesh(fms, me, me_index, mf_left, mf_right)) { + continue; + } + if (dbg_level > 0) { + std::cout << "Removing edge " << me_index << "\n"; + } + splice_faces(fms, me, me_index, mf_left, mf_right); + if (dbg_level > 1) { + std::cout << "state after removal:\n"; + std::cout << *fms; + } + } +} + +/** + * Given that \a tris form a triangulation of a face or part of a face that was in \a imesh_in, + * merge as many of the triangles together as possible, by dissolving the edges between them. + * We can only dissolve triangulation edges that don't overlap real input edges, and we + * can only dissolve them if doing so leaves the remaining faces able to create valid #BMesh. + * We can tell edges that don't overlap real input edges because they will have an + * "original edge" that is different from #NO_INDEX. + * \note it is possible that some of the triangles in \a tris have reversed orientation + * to the rest, so we have to handle the two cases separately. + */ +static Vector<Face *> merge_tris_for_face(Vector<int> tris, + const IMesh &tm, + const IMesh &imesh_in, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "merge_tris_for_face\n"; + std::cout << "tris: " << tris << "\n"; + } + Vector<Face *> ans; + if (tris.size() <= 1) { + if (tris.size() == 1) { + ans.append(tm.face(tris[0])); + } + return ans; + } + bool done = false; + double3 first_tri_normal = tm.face(tris[0])->plane->norm; + double3 second_tri_normal = tm.face(tris[1])->plane->norm; + if (tris.size() == 2 && double3::dot(first_tri_normal, second_tri_normal) > 0.0) { + /* Is this a case where quad with one diagonal remained unchanged? + * Worth special handling because this case will be very common. */ + Face &tri1 = *tm.face(tris[0]); + Face &tri2 = *tm.face(tris[1]); + Face *in_face = imesh_in.face(tri1.orig); + if (in_face->size() == 4) { + std::pair<int, int> estarts = find_tris_common_edge(tri1, tri2); + if (estarts.first != -1 && tri1.edge_orig[estarts.first] == NO_INDEX) { + if (dbg_level > 0) { + std::cout << "try recovering orig quad case\n"; + std::cout << "tri1 = " << &tri1 << "\n"; + std::cout << "tri1 = " << &tri2 << "\n"; + } + int i0 = estarts.first; + int i1 = (i0 + 1) % 3; + int i2 = (i0 + 2) % 3; + int j2 = (estarts.second + 2) % 3; + Face tryface({tri1[i1], tri1[i2], tri1[i0], tri2[j2]}, -1, -1, {}, {}); + if (tryface.cyclic_equal(*in_face)) { + if (dbg_level > 0) { + std::cout << "inface = " << in_face << "\n"; + std::cout << "quad recovery worked\n"; + } + ans.append(in_face); + done = true; + } + } + } + } + if (done) { + return ans; + } + + double3 first_tri_normal_rev = -first_tri_normal; + for (const double3 &norm : {first_tri_normal, first_tri_normal_rev}) { + FaceMergeState fms; + init_face_merge_state(&fms, tris, tm, norm); + do_dissolve(&fms); + if (dbg_level > 0) { + std::cout << "faces in merged result:\n"; + } + for (const MergeFace &mf : fms.face) { + if (mf.merge_to == -1) { + Array<int> e_orig(mf.edge.size()); + Array<bool> is_intersect(mf.edge.size()); + for (int i : mf.edge.index_range()) { + e_orig[i] = fms.edge[mf.edge[i]].orig; + is_intersect[i] = fms.edge[mf.edge[i]].is_intersect; + } + Face *facep = arena->add_face(mf.vert, mf.orig, e_orig, is_intersect); + ans.append(facep); + if (dbg_level > 0) { + std::cout << " " << facep << "\n"; + } + } + } + } + return ans; +} + +/** + * Return an array, paralleling imesh_out.vert, saying which vertices can be dissolved. + * A vertex v can be dissolved if (a) it is not an input vertex; (b) it has valence 2; + * and (c) if v's two neighboring vertices are u and w, then (u,v,w) forms a straight line. + * Return the number of dissolvable vertices in r_count_dissolve. + */ +static Array<bool> find_dissolve_verts(IMesh &imesh_out, int *r_count_dissolve) +{ + imesh_out.populate_vert(); + /* dissolve[i] will say whether imesh_out.vert(i) can be dissolved. */ + Array<bool> dissolve(imesh_out.vert_size()); + for (int v_index : imesh_out.vert_index_range()) { + const Vert &vert = *imesh_out.vert(v_index); + dissolve[v_index] = (vert.orig == NO_INDEX); + } + /* neighbors[i] will be a pair giving the up-to-two neighboring vertices + * of the vertex v in position i of imesh_out.vert. + * If we encounter a third, then v will not be dissolvable. */ + Array<std::pair<const Vert *, const Vert *>> neighbors( + imesh_out.vert_size(), std::pair<const Vert *, const Vert *>(nullptr, nullptr)); + for (int f : imesh_out.face_index_range()) { + const Face &face = *imesh_out.face(f); + for (int i : face.index_range()) { + const Vert *v = face[i]; + int v_index = imesh_out.lookup_vert(v); + BLI_assert(v_index != NO_INDEX); + if (dissolve[v_index]) { + const Vert *n1 = face[face.next_pos(i)]; + const Vert *n2 = face[face.prev_pos(i)]; + const Vert *f_n1 = neighbors[v_index].first; + const Vert *f_n2 = neighbors[v_index].second; + if (f_n1 != nullptr) { + /* Already has a neighbor in another face; can't dissolve unless they are the same. */ + if (!((n1 == f_n2 && n2 == f_n1) || (n1 == f_n1 && n2 == f_n2))) { + /* Different neighbors, so can't dissolve. */ + dissolve[v_index] = false; + } + } + else { + /* These are the first-seen neighbors. */ + neighbors[v_index] = std::pair<const Vert *, const Vert *>(n1, n2); + } + } + } + } + int count = 0; + for (int v_out : imesh_out.vert_index_range()) { + if (dissolve[v_out]) { + dissolve[v_out] = false; /* Will set back to true if final condition is satisfied. */ + const std::pair<const Vert *, const Vert *> &nbrs = neighbors[v_out]; + if (nbrs.first != nullptr) { + BLI_assert(nbrs.second != nullptr); + const mpq3 &co1 = nbrs.first->co_exact; + const mpq3 &co2 = nbrs.second->co_exact; + const mpq3 &co = imesh_out.vert(v_out)->co_exact; + mpq3 dir1 = co - co1; + mpq3 dir2 = co2 - co; + mpq3 cross = mpq3::cross(dir1, dir2); + if (cross[0] == 0 && cross[1] == 0 && cross[2] == 0) { + dissolve[v_out] = true; + ++count; + } + } + } + } + if (r_count_dissolve != nullptr) { + *r_count_dissolve = count; + } + return dissolve; +} + +/** + * The dissolve array parallels the `imesh.vert` array. Wherever it is true, + * remove the corresponding vertex from the vertices in the faces of + * `imesh.faces` to account for the close-up of the gaps in `imesh.vert`. + */ +static void dissolve_verts(IMesh *imesh, const Array<bool> dissolve, IMeshArena *arena) +{ + constexpr int inline_face_size = 100; + Vector<bool, inline_face_size> face_pos_erase; + for (int f : imesh->face_index_range()) { + const Face &face = *imesh->face(f); + face_pos_erase.clear(); + int num_erase = 0; + for (const Vert *v : face) { + int v_index = imesh->lookup_vert(v); + BLI_assert(v_index != NO_INDEX); + if (dissolve[v_index]) { + face_pos_erase.append(true); + ++num_erase; + } + else { + face_pos_erase.append(false); + } + } + if (num_erase > 0) { + imesh->erase_face_positions(f, face_pos_erase, arena); + } + } + imesh->set_dirty_verts(); +} + +/** + * The main boolean function operates on a triangle #IMesh and produces a + * Triangle #IMesh as output. + * This function converts back into a general polygonal mesh by removing + * any possible triangulation edges (which can be identified because they + * will have an original edge that is NO_INDEX. + * Not all triangulation edges can be removed: if they ended up non-trivially overlapping a real + * input edge, then we need to keep it. Also, some are necessary to make the output satisfy + * the "valid #BMesh" property: we can't produce output faces that have repeated vertices in them, + * or have several disconnected boundaries (e.g., faces with holes). + */ +static IMesh polymesh_from_trimesh_with_dissolve(const IMesh &tm_out, + const IMesh &imesh_in, + IMeshArena *arena) +{ + const int dbg_level = 0; + if (dbg_level > 1) { + std::cout << "\nPOLYMESH_FROM_TRIMESH_WITH_DISSOLVE\n"; + } + /* For now: need plane normals for all triangles. */ + for (Face *tri : tm_out.faces()) { + tri->populate_plane(false); + } + /* Gather all output triangles that are part of each input face. + * face_output_tris[f] will be indices of triangles in tm_out + * that have f as their original face. */ + int tot_in_face = imesh_in.face_size(); + Array<Vector<int>> face_output_tris(tot_in_face); + for (int t : tm_out.face_index_range()) { + const Face &tri = *tm_out.face(t); + int in_face = tri.orig; + face_output_tris[in_face].append(t); + } + if (dbg_level > 1) { + std::cout << "face_output_tris:\n"; + for (int f : face_output_tris.index_range()) { + std::cout << f << ": " << face_output_tris[f] << "\n"; + } + } + + /* Merge triangles that we can from face_output_tri to make faces for output. + * face_output_face[f] will be new original const Face *'s that + * make up whatever part of the boolean output remains of input face f. */ + Array<Vector<Face *>> face_output_face(tot_in_face); + int tot_out_face = 0; + for (int in_f : imesh_in.face_index_range()) { + if (dbg_level > 1) { + std::cout << "merge tris for face " << in_f << "\n"; + } + int num_out_tris_for_face = face_output_tris.size(); + if (num_out_tris_for_face == 0) { + continue; + } + face_output_face[in_f] = merge_tris_for_face(face_output_tris[in_f], tm_out, imesh_in, arena); + tot_out_face += face_output_face[in_f].size(); + } + Array<Face *> face(tot_out_face); + int out_f_index = 0; + for (int in_f : imesh_in.face_index_range()) { + const Vector<Face *> &f_faces = face_output_face[in_f]; + if (f_faces.size() > 0) { + std::copy(f_faces.begin(), f_faces.end(), &face[out_f_index]); + out_f_index += f_faces.size(); + } + } + IMesh imesh_out(face); + /* Dissolve vertices that were (a) not original; and (b) now have valence 2 and + * are between two other vertices that are exactly in line with them. + * These were created because of triangulation edges that have been dissolved. */ + int count_dissolve; + Array<bool> v_dissolve = find_dissolve_verts(imesh_out, &count_dissolve); + if (count_dissolve > 0) { + dissolve_verts(&imesh_out, v_dissolve, arena); + } + if (dbg_level > 1) { + write_obj_mesh(imesh_out, "boolean_post_dissolve"); + } + + return imesh_out; +} + +/** + * This function does a boolean operation on a TriMesh with nshapes inputs. + * All the shapes are combined in tm_in. + * The shape_fn function should take a triangle index in tm_in and return + * a number in the range 0 to `nshapes-1`, to say which shape that triangle is in. + */ +IMesh boolean_trimesh(IMesh &tm_in, + BoolOpType op, + int nshapes, + std::function<int(int)> shape_fn, + bool use_self, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "BOOLEAN of " << nshapes << " operand" << (nshapes == 1 ? "" : "s") + << " op=" << bool_optype_name(op) << "\n"; + if (dbg_level > 1) { + tm_in.populate_vert(); + std::cout << "boolean_trimesh input:\n" << tm_in; + write_obj_mesh(tm_in, "boolean_in"); + } + } + if (tm_in.face_size() == 0) { + return IMesh(tm_in); + } + IMesh tm_si; + if (use_self) { + tm_si = trimesh_self_intersect(tm_in, arena); + } + else { + tm_si = trimesh_nary_intersect(tm_in, nshapes, shape_fn, use_self, arena); + } + if (dbg_level > 1) { + write_obj_mesh(tm_si, "boolean_tm_si"); + std::cout << "\nboolean_tm_input after intersection:\n" << tm_si; + } + /* It is possible for tm_si to be empty if all the input triangles are bogus/degenerate. */ + if (tm_si.face_size() == 0 || op == BoolOpType::None) { + return tm_si; + } + auto si_shape_fn = [shape_fn, tm_si](int t) { return shape_fn(tm_si.face(t)->orig); }; + TriMeshTopology tm_si_topo(tm_si); + PatchesInfo pinfo = find_patches(tm_si, tm_si_topo); + IMesh tm_out; + if (!is_pwn(tm_si, tm_si_topo)) { + if (dbg_level > 0) { + std::cout << "Input is not PWN, using gwn method\n"; + } + tm_out = gwn_boolean(tm_si, op, nshapes, shape_fn, pinfo, arena); + } + else { + CellsInfo cinfo = find_cells(tm_si, tm_si_topo, pinfo); + if (dbg_level > 0) { + std::cout << "Input is PWN\n"; + } + finish_patch_cell_graph(tm_si, cinfo, pinfo, tm_si_topo, arena); + bool pc_ok = patch_cell_graph_ok(cinfo, pinfo); + if (!pc_ok) { + /* TODO: if bad input can lead to this, diagnose the problem. */ + std::cout << "Something funny about input or a bug in boolean\n"; + return IMesh(tm_in); + } + cinfo.init_windings(nshapes); + int c_ambient = find_ambient_cell(tm_si, nullptr, tm_si_topo, pinfo, arena); + if (c_ambient == NO_INDEX) { + /* TODO: find a way to propagate this error to user properly. */ + std::cout << "Could not find an ambient cell; input not valid?\n"; + return IMesh(tm_si); + } + propagate_windings_and_in_output_volume(pinfo, cinfo, c_ambient, op, nshapes, si_shape_fn); + tm_out = extract_from_in_output_volume_diffs(tm_si, pinfo, cinfo, arena); + if (dbg_level > 0) { + /* Check if output is PWN. */ + TriMeshTopology tm_out_topo(tm_out); + if (!is_pwn(tm_out, tm_out_topo)) { + std::cout << "OUTPUT IS NOT PWN!\n"; + } + } + } + if (dbg_level > 1) { + write_obj_mesh(tm_out, "boolean_tm_output"); + std::cout << "boolean tm output:\n" << tm_out; + } + return tm_out; +} + +static void dump_test_spec(IMesh &imesh) +{ + std::cout << "test spec = " << imesh.vert_size() << " " << imesh.face_size() << "\n"; + for (const Vert *v : imesh.vertices()) { + std::cout << v->co_exact[0] << " " << v->co_exact[1] << " " << v->co_exact[2] << " # " + << v->co[0] << " " << v->co[1] << " " << v->co[2] << "\n"; + } + for (const Face *f : imesh.faces()) { + for (const Vert *fv : *f) { + std::cout << imesh.lookup_vert(fv) << " "; + } + std::cout << "\n"; + } +} + +/** + * Do the boolean operation op on the polygon mesh imesh_in. + * See the header file for a complete description. + */ +IMesh boolean_mesh(IMesh &imesh, + BoolOpType op, + int nshapes, + std::function<int(int)> shape_fn, + bool use_self, + IMesh *imesh_triangulated, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nBOOLEAN_MESH\n" + << nshapes << " operand" << (nshapes == 1 ? "" : "s") + << " op=" << bool_optype_name(op) << "\n"; + if (dbg_level > 1) { + write_obj_mesh(imesh, "boolean_mesh_in"); + std::cout << imesh; + if (dbg_level > 2) { + dump_test_spec(imesh); + } + } + } + IMesh *tm_in = imesh_triangulated; + IMesh our_triangulation; + if (tm_in == nullptr) { + our_triangulation = triangulate_polymesh(imesh, arena); + tm_in = &our_triangulation; + } + if (dbg_level > 1) { + write_obj_mesh(*tm_in, "boolean_tm_in"); + } + IMesh tm_out = boolean_trimesh(*tm_in, op, nshapes, shape_fn, use_self, arena); + if (dbg_level > 1) { + std::cout << "bool_trimesh_output:\n" << tm_out; + write_obj_mesh(tm_out, "bool_trimesh_output"); + } + IMesh ans = polymesh_from_trimesh_with_dissolve(tm_out, imesh, arena); + if (dbg_level > 0) { + std::cout << "boolean_mesh output:\n" << ans; + } + return ans; +} + +} // namespace blender::meshintersect + +#endif // WITH_GMP |