text: /*
   * 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);
    }
  
    Span<int> patches() const
    {
      return Span<int>(patches_);
    }
  
    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());
      if (shape >= 0) {
        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 << ": " << tm.face(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] >= 1) {
            return false;
          }
        }
        return true;
      }
      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.
   * If \a high_confidence is true then we want a higher degree
   * of "insideness" than if it is false.
   */
  static bool point_is_inside_shape(const IMesh &tm,
                                    std::function<int(int)> shape_fn,
                                    const double3 &testp,
                                    int shape,
                                    bool high_confidence)
  {
    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. */
    if (high_confidence) {
      return (gwn > 0.9);
    }
  
    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());
    Array<int> winding(nshapes, 0);
    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";
      }
      for (int other_shape = 0; other_shape < nshapes; ++other_shape) {
        if (other_shape == shape) {
          continue;
        }
        /* The point_is_inside_shape function has to approximate if the other
         * shape is not PWN. For most operations, even a hint of being inside
         * gives good results, but when shape is a cutter in a Difference
         * operation, we want to be pretty sure that the point is inside other_shape.
         * E.g., T75827.
         */
        bool need_high_confidence = (op == BoolOpType::Difference) && (shape != 0);
        bool inside = point_is_inside_shape(
            tm, shape_fn, test_point_db, other_shape, need_high_confidence);
        if (dbg_level > 0) {
          std::cout << "test point is " << (inside ? "inside" : "outside") << " other_shape "
                    << other_shape << "\n";
        }
        winding[other_shape] = inside;
      }
      /* Find out the "in the output volume" flag for each of the cases of winding[shape] == 0
       * and winding[shape] == 1. If the flags are different, this patch should be in the output.
       * Also, if this is a Difference and the shape isn't the first one, need to flip the normals.
       */
      winding[shape] = 0;
      bool in_output_volume_0 = apply_bool_op(op, winding);
      winding[shape] = 1;
      bool in_output_volume_1 = apply_bool_op(op, winding);
      bool do_remove = in_output_volume_0 == in_output_volume_1;
      bool do_flip = !do_remove && op == BoolOpType::Difference && shape != 0;
      if (dbg_level > 0) {
        std::cout << "winding = ";
        for (int i = 0; i < nshapes; ++i) {
          std::cout << winding[i] << " ";
        }
        std::cout << "\niv0=" << in_output_volume_0 << ", iv1=" << in_output_volume_1 << "\n";
        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);
          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;
    bool any_faces_erased = false;
    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) {
        any_faces_erased |= imesh->erase_face_positions(f, face_pos_erase, arena);
      }
    }
    imesh->set_dirty_verts();
    if (any_faces_erased) {
      imesh->remove_null_faces();
    }
  }
  
  /**
   * 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 = 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;
      if (dbg_level > 2) {
        ans.populate_vert();
        dump_test_spec(ans);
      }
    }
    return ans;
  }
  
  }  // namespace blender::meshintersect
  
  #endif  // WITH_GMP
  