ochafik/__ Coplanar remeshing in OpenSCAD.md Secret

## __ Coplanar remeshing in OpenSCAD.md

      
    Raw
  

              __ Coplanar remeshing in OpenSCAD.md
            
          
    Note: the code below is slightly ahead of openscad/openscad#4117 as of this writing, needs cleanup!
Manual remesh to cleanup the output of CGAL corefinement functions. Hooking in as corefinement visitors, as suggested by @sloriot here.
If this works out well we won't need to wait for CGAL 5.5 to fix the hanging cases of fast-csg (see openscad/openscad#4107)
Note: Somehow visitors get copied so I added a delegate to share the state between all visitors.


## _cgalutils-corefinement-visitor.h
// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
#pragma once

#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include "cgalutils-coplanar-faces-remesher.h"
#include "cgalutils-mesh-edits.h"

namespace CGALUtils {

namespace internal {

template <typename Pt>
void forceExact(Pt& point)
{
  CGAL::exact(point.x());
  CGAL::exact(point.y());
  CGAL::exact(point.z());
}

/*! Implements corefinement callbacks to support remeshing coplanar faces and forcing new numbers to
 * be exact. */
template <typename TriangleMesh>
struct CorefinementVisitorDelegate_ { // : public PMP::Corefinement::Default_visitor<TriangleMesh> {
private:
  typedef boost::graph_traits<TriangleMesh> GT;
  typedef typename GT::face_descriptor face_descriptor;
  typedef typename GT::halfedge_descriptor halfedge_descriptor;
  typedef typename GT::edge_descriptor edge_descriptor;
  typedef typename GT::vertex_descriptor vertex_descriptor;

  typedef size_t PatchId; // Starting from 1

  bool forceNewLazyNumbersToExact_;
  TriangleMesh *mesh1_, *mesh2_, *meshOut_;

  // The following fields are scoped by mesh ([mesh1_, mesh2_, meshOut_]).
  // If meshOut_ == mesh1_ (or mesh2_), the third element of each of these
  // fields will go unused, as all accesses are using the index given by indexOf

  std::unordered_map<face_descriptor, PatchId> faceToPatchId_[3];
  std::unordered_set<PatchId> splitPatchIds_;

  face_descriptor faceBeingSplit_[3];
  std::vector<face_descriptor> facesBeingCreated_[3];

  PatchId getPatchId(const face_descriptor& f, const TriangleMesh& tm) {
    auto id = faceToPatchId_[indexOf(tm)][f];
    CGAL_assertion(id);
    return id;
  }

  void setPatchId(const face_descriptor& f, const TriangleMesh& tm, PatchId id) {
    CGAL_assertion(id);
    faceToPatchId_[indexOf(tm)][f] = id;
  }

  size_t indexOf(const TriangleMesh& tm) {
    if (&tm == mesh1_) return 0;
    if (&tm == mesh2_) return 1;
    if (&tm == meshOut_) return 2;
    CGAL_assertion(false && "Unknown mesh");
  }

public:
  CorefinementVisitorDelegate_(TriangleMesh& m1, TriangleMesh& m2, TriangleMesh& mout,
                               bool forceNewLazyNumbersToExact)
    : mesh1_(&m1), mesh2_(&m2), meshOut_(&mout), forceNewLazyNumbersToExact_(forceNewLazyNumbersToExact)
  {
    TriangleMesh *meshes[3] = {&m1, &m2, &mout};

    // Give each existing face its own global patch id.
    PatchId nextPatchId = 1;
    for (auto mi = 0; mi < 3; mi++) {
      auto& m = *meshes[mi];

      if (mi != indexOf(m)) {
        // &mout is often == &m1, so don't reassign its patch ids!
        continue;
      }

      for (auto f : m.faces()) setPatchId(f, m, nextPatchId++);
    }
  }

  void remeshSplitFaces(TriangleMesh& tm)
  {
    auto mi = indexOf(tm);
    CoplanarFacesRemesher<TriangleMesh> remesher(tm);
    remesher.remesh(faceToPatchId_[mi], splitPatchIds_);
  }

  void after_face_copy(face_descriptor f_old, const TriangleMesh& m_old, face_descriptor f_new,
                       TriangleMesh& m_new)
  {
    setPatchId(f_new, m_new, getPatchId(f_old, m_old));
  }

  void before_subface_creations(face_descriptor f_split, TriangleMesh& tm)
  {
    auto mi = indexOf(tm);
    faceBeingSplit_[mi] = f_split;
    facesBeingCreated_[mi].clear();
  }

  void after_subface_created(face_descriptor fi, TriangleMesh& tm)
  {
    auto mi = indexOf(tm);
    auto id = getPatchId(faceBeingSplit_[mi], tm);
    setPatchId(fi, tm, getPatchId(faceBeingSplit_[mi], tm));
    // splitPatchIds_.insert(id); // Doing it once in after_subface_creations instead.

#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
    if (forceNewLazyNumbersToExact_) {
      facesBeingCreated_[mi].push_back(fi);
    }
#endif
  }

  void after_subface_creations(TriangleMesh& tm)
  {
    auto mi = indexOf(tm);
    auto id = getPatchId(faceBeingSplit_[mi], tm);
    splitPatchIds_.insert(id);

    // From CGAL 5.4 on we rely on new_vertex_added instead.
#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
    if (forceNewLazyNumbersToExact_) {
      for (auto& fi : facesBeingCreated_[mi]) {
        CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend;
        for (boost::tie(vbegin, vend) = vertices_around_face(tm.halfedge(fi), tm); vbegin != vend;
             ++vbegin) {
          forceExact(tm.point(*vbegin));
        }
      }
    }
#endif // if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
  }

  // Note: only called in CGAL 5.4+
  void new_vertex_added(std::size_t node_id, vertex_descriptor vh, const TriangleMesh& tm)
  {
    if (forceNewLazyNumbersToExact_) {
      forceExact(tm.point(vh));
    }
  }

};

} // namespace internal

template <typename TriangleMesh>
struct CorefinementVisitor : public PMP::Corefinement::Default_visitor<TriangleMesh> {
  typedef boost::graph_traits<TriangleMesh> GT;
  typedef typename GT::face_descriptor face_descriptor;
  typedef typename GT::vertex_descriptor vertex_descriptor;

  typedef internal::CorefinementVisitorDelegate_<TriangleMesh> Delegate;

  std::shared_ptr<Delegate> delegate_;

  CorefinementVisitor(TriangleMesh& m1, TriangleMesh& m2, TriangleMesh& mout,
                      bool forceNewLazyNumbersToExact)
    : delegate_(std::make_shared<Delegate>(m1, m2, mout, forceNewLazyNumbersToExact))
  {
  }

  void remeshSplitFaces(TriangleMesh& tm) {
    delegate_->remeshSplitFaces(tm);
  }

  void before_subface_creations(face_descriptor f_old, TriangleMesh& tm) {
    delegate_->before_subface_creations(f_old, tm);
  }
  void after_subface_creations(TriangleMesh& tm) {
    delegate_->after_subface_creations(tm);
  }
  void after_subface_created(face_descriptor f_new, TriangleMesh& tm) {
    delegate_->after_subface_created(f_new, tm);
  }
  void after_face_copy(face_descriptor f_old, const TriangleMesh& old_tm, face_descriptor f_new, TriangleMesh& new_tm) {
    delegate_->after_face_copy(f_old, old_tm, f_new, new_tm);
  }
  void new_vertex_added(std::size_t node_id, vertex_descriptor vh, const TriangleMesh& tm) {
    delegate_->new_vertex_added(node_id, vh, tm);
  }
};

} // namespace CGALUtils

## cgal-faces-remesher.h
// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
#pragma once

#include <iterator>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Kernel/global_functions.h>
#include <CGAL/Kernel/function_objects.h>
// #include <CGAL/Polygon_mesh_processing/stitch_borders.h>

#include "cgalutils-mesh-edits.h"

namespace CGALUtils {

/*! Helper that remeshes coplanar faces.
 *
 * This is useful for the results of corefinement operations, which result in
 * suboptimal meshes which growth compounds over several operations and makes
 * some models exponentially slow.
 *
 * Two faces with the same patch id *MUST* be coplanar, but need not be continuous
 * (technically the same patch id could result in different patches, e.g. if
 * during corefinement the same original face got split in disjoint sets of
 * descendants).
 */
template <typename TriangleMesh>
struct CoplanarFacesRemesher {
  typedef boost::graph_traits<TriangleMesh> GT;
  typedef typename GT::face_descriptor face_descriptor;
  typedef typename GT::halfedge_descriptor halfedge_descriptor;
  typedef typename GT::edge_descriptor edge_descriptor;
  typedef typename GT::vertex_descriptor vertex_descriptor;

  typedef std::function<bool (const face_descriptor&)> face_predicate;

  TriangleMesh& tm;

public:

  CoplanarFacesRemesher(TriangleMesh& mesh) : tm(mesh) {}

  template <typename PatchId>
  void remesh(
    const std::unordered_map<face_descriptor, PatchId>& faceToCoplanarPatchId,
    const std::unordered_set<PatchId>& patchesToRemesh)
  {
    if (patchesToRemesh.empty()) return;

    auto facesBefore = tm.number_of_faces();
    auto verbose = Feature::ExperimentalFastCsgDebug.is_enabled();

    // PMP::stitch_borders(tm);

    auto printVertex = [&](auto v) {
        auto p = tm.point(v);
        //std::cerr << p.x() <<  " -> ";
        printf("[%g, %g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()));
        // printf("[%g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()));
      };

    auto printVertices = [&](auto& path) {
        for (auto v : path) {
          printVertex(v);
          printf(", ");
        }
        printf("\n");
      };

    auto printHEdges = [&](auto& path) {
        for (auto he : path) {
          std::cout << he << " (";
          printVertex(tm.source(he));
          std::cout << " -> ";
          printVertex(tm.target(he));
          std::cout << "), ";
        }
        printf("\n");
      };

    auto patchToPolyhedronStr = [&](auto &patchFaces) {
      std::vector<vertex_descriptor> vertices;
      std::map<vertex_descriptor, size_t> vertexIndex;

      std::ostringstream verticesOut;
      std::ostringstream facesOut;

      facesOut << "[\n";
      for (auto& face : patchFaces) {
        CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
        facesOut << "  [";
        auto first = true;
        for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
          auto he = *heIt;
          auto v = tm.source(he);
          auto it = vertexIndex.find(v);
          size_t idx;
          if (it == vertexIndex.end()) {
            idx = vertices.size();
            vertices.push_back(v);
            vertexIndex[v] = idx;
          } else {
            idx = it->second;
          }
          if (first) first = false;
          else facesOut << ", ";
          facesOut << idx;
        }
        facesOut << "],\n";
      }
      facesOut << "]";

      verticesOut << "[\n";
      for (auto v : vertices) {
        auto p = tm.point(v);
        verticesOut << "[" << CGAL::to_double(p.x()) << ", " << CGAL::to_double(p.y()) << ", " << CGAL::to_double(p.z()) << "],\n";
      }
      verticesOut << "]";

      std::ostringstream out;
      out << "polyhedron(" << verticesOut.str().c_str() << ", " << facesOut.str().c_str() << ");";
      return out.str();
    };

    auto patchBorderToPolyhedronStr = [&](auto &borderPathVertices) {
      std::vector<vertex_descriptor> vertices;
      std::map<vertex_descriptor, size_t> vertexIndex;

      std::ostringstream verticesOut;
      std::ostringstream facesOut;

      facesOut << "[[";
      auto first = true;
      for (auto& v : borderPathVertices) {
        auto it = vertexIndex.find(v);
        size_t idx;
        if (it == vertexIndex.end()) {
          idx = vertices.size();
          vertices.push_back(v);
          vertexIndex[v] = idx;
        } else {
          idx = it->second;
        }
        if (first) first = false;
        else facesOut << ", ";
        facesOut << idx;
      }
      facesOut << "]]\n";

      verticesOut << "[\n";
      for (auto v : vertices) {
        auto p = tm.point(v);
        verticesOut << "[" << CGAL::to_double(p.x()) << ", " << CGAL::to_double(p.y()) << ", " << CGAL::to_double(p.z()) << "],\n";
      }
      verticesOut << "]";

      std::ostringstream out;
      out << "polyhedron(" << verticesOut.str().c_str() << ", " << facesOut.str().c_str() << ");";
      return out.str();
    };

    try {
      TriangleMeshEdits<TriangleMesh> edits;
      std::unordered_set<PatchId> patchesProcessed;
      patchesProcessed.reserve(patchesToRemesh.size());

      std::unordered_set<face_descriptor> facesProcessed;
      // patchesProcessed.reserve(patchesToRemesh.size());

      std::set<face_descriptor> loopLocalPatchFaces;
      std::map<PatchId, bool> loopLocalIsPatchIdMap;
      std::vector<halfedge_descriptor> loopLocalBorderEdges;
      std::vector<halfedge_descriptor> loopLocalBorderPath;
      std::vector<vertex_descriptor> loopLocalBorderPathVertices;
      std::unordered_set<halfedge_descriptor> loopLocalBorderPathEdges;

      std::unordered_map<PatchId, std::string> allPatchPolyStrings;
      std::unordered_map<PatchId, std::string> allPatchReplacementsPolyStrings;

      for (auto face : tm.faces()) {
        if (tm.is_removed(face)) {
          continue;
        }
        auto idIt = faceToCoplanarPatchId.find(face);
        if (idIt == faceToCoplanarPatchId.end()) {
          continue;
        }
        const auto id = idIt->second;

        if (!getenv("COPL") && (patchesToRemesh.find(id) == patchesToRemesh.end())) {
          continue;
        }
        if (!patchesProcessed.insert(id).second) {
          continue;
        }
        if (facesProcessed.find(face) != facesProcessed.end()) {
          continue;
        }

        loopLocalIsPatchIdMap.clear();
        auto& isPatchIdMap = loopLocalIsPatchIdMap;
        isPatchIdMap[id] = true;

        loopLocalPatchFaces.clear();
        auto& patchFaces = loopLocalPatchFaces;
        patchFaces.insert(face);

        auto isHalfedgeOnBorder = [&](auto& he) {
            auto neighbourFace = tm.face(tm.opposite(he));
            if (!neighbourFace.is_valid()) {
              return true;
            }

            auto neighbourIdIt = faceToCoplanarPatchId.find(neighbourFace);
            if (neighbourIdIt == faceToCoplanarPatchId.end()) {
              std::cerr << "Unknown face, weird!";
              return true;
            }
            auto neighbourId = neighbourIdIt->second;

            auto isPatchIdIt = isPatchIdMap.find(neighbourId);
            if (isPatchIdIt != isPatchIdMap.end()) {
              auto isPatchId = isPatchIdIt->second;
              return !isPatchId;
            }

            if (getenv("COPL")) {
              auto isCoplanar = isEdgeBetweenCoplanarFaces(he);
              // std::cerr << "isCoplanar(" << face << ", " << neighbourFace << "): " << (isCoplanar ? "true" : "false") << "\n";
              isPatchIdMap[neighbourId] = isCoplanar;
              return !isCoplanar;
            }
            // We're on the border
            return true;
          };

        floodFillPatch(patchFaces, isHalfedgeOnBorder);

        allPatchPolyStrings[id] = patchToPolyhedronStr(patchFaces);

        auto conflicting = false;
        for (auto& face : patchFaces) {
          if (!facesProcessed.insert(face).second) {
            std::cerr << "Patch " << id << " trying to process already processed face " << face << "\n";
            conflicting = true;
            break;
          }
        }
        if (conflicting) {
          continue;
        }

        if (patchFaces.size() < 2) {
          continue;
        }

        auto isFaceOnPatch = [&](auto& f) {
            return patchFaces.find(f) != patchFaces.end();
          };

        loopLocalBorderEdges.clear();
        auto &borderEdges = loopLocalBorderEdges;

        for (auto& face : patchFaces) {
          CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
          for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
            auto he = *heIt;
            if (isHalfedgeOnBorder(he)) {
              borderEdges.push_back(he);
            }
          }
        }

        if (borderEdges.empty()) {
          std::cerr << "Failed to find a border edge for patch " << id << "!\n";
          continue;
        }
        halfedge_descriptor borderEdge = *borderEdges.begin();

        loopLocalBorderPath.clear();
        auto& borderPath = loopLocalBorderPath;

        if (!walkAroundPatch(borderEdge, isFaceOnPatch, borderPath)) {
          LOG(message_group::Error, Location::NONE, "",
              "[fast-csg-remesh] Failed to collect path around patch faces, invalid mesh!");
          return;
        }

        if (borderPath.size() <= 3) {
          continue;
        }

        loopLocalBorderPathVertices.clear();
        auto& borderPathVertices = loopLocalBorderPathVertices;

        loopLocalBorderPathEdges.clear();
        auto &borderPathEdges = loopLocalBorderPathEdges;

        for (auto he : borderPath) {
          borderPathEdges.insert(he);
          borderPathVertices.push_back(tm.target(he));
        }

        auto hasHoles = false;
        for (auto he : borderEdges) {
          if (borderPathEdges.find(he) == borderPathEdges.end()) {
            // Found a border halfedge that's not in the border path we've walked around the patch.
            // This means the patch is holed / has more than one border: abort!!
            hasHoles = true;
            break;
          }
        }
        if (hasHoles) {
          if (verbose) {
            LOG(message_group::None, Location::NONE, "",
              "[fast-csg-remesh] Skipping remeshing of patch with %1$lu faces as it seems to have holes.", borderPathEdges.size());
          }
          continue;
        }

        // auto veryVerbose = verbose;
        // // Only on z=1 plane to reduce clutter
        // for (auto he : borderPath) if (CGAL::to_double(tm.point(tm.target(he)).z()) != 1) veryVerbose = false;
        // if (veryVerbose) {

        //   std::cerr << borderPath.size() << " edges around patch " << id << ":\n";
        //   printHEdges(borderPath);
        //     std::cout << "\n";

        //   std::cerr << patchFaces.size() << " faces in patch " << id << ":\n";
        //   for (auto f : patchFaces) {
        //     std::cout << "\tFace " << f << ": ";
        //     std::vector<vertex_descriptor> vertices;
        //     std::vector<halfedge_descriptor> hedges;
        //     CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
        //     for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(f), tm); heIt != heEnd; ++heIt) {
        //       auto he = *heIt;
        //       vertices.push_back(tm.target(he));
        //       hedges.push_back(he);
        //     }
        //     printVertices(vertices);
        //     printHEdges(hedges);
        //     std::cout << "\n";
        //   }

        //   std::cerr << "Patch " << id << ": removing " << patchFaces.size() << " faces\n";

        //   std::cout << "Border of length " << borderPathVertices.size() << " for patch " << id << ": ";
        //   printVertices(borderPathVertices);
        // }

        // Only remesh patches that have consecutive collinear edges.
        auto lengthBefore = borderPathVertices.size();
        auto collapsed = TriangleMeshEdits<TriangleMesh>::collapsePathWithConsecutiveCollinearEdges(borderPathVertices, tm);

          auto lengthAfter = borderPathVertices.size();
          if (collapsed) {
            std::cerr << "Collapsed path around patch " << id << " (" << patchFaces.size() << " faces) from " << lengthBefore << " to " << lengthAfter << " vertices\n";
            if (lengthBefore == lengthAfter) {
              std::cerr << "# WTF??\n";
              // continue;
            }
          }

          allPatchReplacementsPolyStrings[id] = patchBorderToPolyhedronStr(borderPathVertices);

          auto onlyPatchIdStr = getenv("PATCHID");
          if (onlyPatchIdStr && (atoi(onlyPatchIdStr) != id)) {
            std::cerr << "Skipping patch " << id << " (PATCHID = " << atoi(onlyPatchIdStr)  << ")\n";
            continue;
          }

          // Cover the patch with a polygon. It will be triangulated when the
          // edits are applied.

          // edits.replaceFaces(patchFaces, borderPathVertices);

          for (auto& face : patchFaces) edits.removeFace(face);
          edits.addFace(borderPathVertices);

      }

      {
        std::ofstream fout("patches.scad");
        fout << "before=true;\npatchIndex=-1;\n";
        size_t patchIndex = 0;
        for (auto &p : allPatchPolyStrings) {
          auto id = p.first;
          fout << "// Patch id " << id << " (index " << patchIndex << "):\n";
          if (allPatchReplacementsPolyStrings[id].empty()) fout << "%";
          fout << "if (patchIndex < 0 || patchIndex == " << patchIndex << ") {\n";
          fout << "  if (before) { " << p.second.c_str() << " }\n";
          fout << "  else { " << allPatchReplacementsPolyStrings[id] << " }\n";
          fout << "}\n";
          patchIndex++;
        }
        std::ofstream fout("patches.scad");
      }

      if (!edits.isEmpty()) {
        edits.apply(tm);
      }

      auto facesAfter = tm.number_of_faces();
      if (verbose) {
        LOG(message_group::None, Location::NONE, "",
            "[fast-csg-remesh] Remeshed from %1$lu to %2$lu faces (%3$lu % improvement)", facesBefore, facesAfter,
            (facesBefore - facesAfter) * 100 / facesBefore);
      }
    } catch (const CGAL::Assertion_exception& e) {
      LOG(message_group::Error, Location::NONE, "", "CGAL error in remeshSplitFaces: %1$s", e.what());
    }
  }

private:

  bool isEdgeBetweenCoplanarFaces(const halfedge_descriptor& h) {
    auto& p = tm.point(tm.source(h));
    auto& q = tm.point(tm.target(h));
    auto& r = tm.point(tm.target(tm.next(h)));
    auto& s = tm.point(tm.target(tm.next(tm.opposite(h))));

    return CGAL::coplanar(p, q, r, s);
    // auto coplanarCosThreshold = getenv("COPL") ? atof(getenv("COPL")) : -1;
    // return coplanarCosThreshold < 0
    //   ? CGAL::coplanar(p, q, r, s)
    //   : CGAL::compare_dihedral_angle(p, q, r, s, coplanarCosThreshold) == CGAL::LARGER;
  }

  /*! Expand the patch of known faces to surrounding neighbours that pass the (fast) predicate. */
  bool floodFillPatch(
    std::set<face_descriptor>& facePatch,
    const std::function<bool(const halfedge_descriptor&)>& isHalfedgeOnBorder)
  {
    // Avoid recursion as its depth would be model-dependent.
    std::vector<face_descriptor> unprocessedFaces(facePatch.begin(), facePatch.end());


    while (!unprocessedFaces.empty()) {
      auto face = unprocessedFaces.back();
      unprocessedFaces.pop_back();

      CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
      for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
        auto he = *heIt;
        if (isHalfedgeOnBorder(he)) {
          continue;
        }

        auto neighbourFace = tm.face(tm.opposite(he));
        if (!facePatch.insert(neighbourFace).second) {
          continue;
        }

        unprocessedFaces.push_back(neighbourFace);
      }
    }

    return true;
  }

  /*! Returns true if the border is closed. */
  bool walkAroundPatch(halfedge_descriptor startingEdge,
                       const face_predicate& isFaceOnPatch,
                       std::vector<halfedge_descriptor>& borderOut)
  {
    auto currentEdge = startingEdge;
    borderOut.push_back(startingEdge);
    // std::set<vertex_descriptor> visitedVertices;

    std::unordered_set<vertex_descriptor> visitedVertices;

    while (tm.target(currentEdge) != tm.source(startingEdge)) {
      // visitedVertices.insert(tm.target(currentEdge));

      auto foundNext = false;
      CGAL::Halfedge_around_source_iterator<TriangleMesh> heIt, heEnd;
      for (boost::tie(heIt, heEnd) = halfedges_around_source(tm.target(currentEdge), tm); heIt != heEnd; ++heIt) {
        auto he = *heIt;
        if (he == tm.opposite(currentEdge)) {
          // Don't go back to where we just came from so quickly!
          continue;
        }
        // if (visitedVertices.find(tm.target(he)) != visitedVertices.end()) {
        //   continue;
        // }

        if (!isFaceOnPatch(tm.face(he))) {
          // This halfedge's face isn't on the patch, ignore it.
          continue;
        }

        if (isFaceOnPatch(tm.face(tm.opposite(he)))) {
          // Not a border edge as both faces are on the patch
          continue;
        }

        auto v = tm.target(he);
        if (!visitedVertices.insert(v).second) {
          continue;
        }

        // std::cerr << "Edge " << he << " between " << tm.face(he) << " and " << tm.face(tm.opposite(he)) << " is on the border!\n";

        // Edge is on the border of the patch.
        currentEdge = he;
        borderOut.push_back(currentEdge);
        foundNext = true;
        break;
      }
      if (!foundNext) {
        std::cerr << "Error: Failed to find the next edge to walk to along the patch border!\n";
        return false;
      }
    }

    return true;
  }
};

} // namespace CGALUtils

## cgalutils-corefine.cc
// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
#include "cgalutils.h"

#include "cgalutils-corefinement-visitor.h"

namespace CGALUtils {

typedef CGAL::Surface_mesh<CGAL::Point_3<CGAL_HybridKernel3>> CGAL_HybridSurfaceMesh;

#if FAST_CSG_KERNEL_IS_LAZY

/*! Visitor that forces exact numbers for the vertices of all the faces created during corefinement.
 */
template <typename TriangleMesh>
struct ExactLazyNumbersVisitor
  : public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<TriangleMesh> {
  typedef boost::graph_traits<TriangleMesh> GT;
  typedef typename GT::face_descriptor face_descriptor;
  typedef typename GT::halfedge_descriptor halfedge_descriptor;
  typedef typename GT::vertex_descriptor vertex_descriptor;

#if CGAL_VERSION_NR >= CGAL_VERSION_NUMBER(5, 4, 0)
  void new_vertex_added(std::size_t i_id, vertex_descriptor v, const TriangleMesh& tm) {
    auto& pt = tm.point(v);
    CGAL::exact(pt.x());
    CGAL::exact(pt.y());
    CGAL::exact(pt.z());
  }
#else
  face_descriptor split_face;
  std::vector<face_descriptor> created_faces;

  void before_subface_creations(face_descriptor f_split, TriangleMesh& tm)
  {
    split_face = f_split;
    created_faces.clear();
  }

  void after_subface_creations(TriangleMesh& mesh)
  {
    for (auto& fi : created_faces) {
      CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend;
      for (boost::tie(vbegin, vend) = vertices_around_face(mesh.halfedge(fi), mesh); vbegin != vend;
           ++vbegin) {
        auto& v = mesh.point(*vbegin);
        CGAL::exact(v.x());
        CGAL::exact(v.y());
        CGAL::exact(v.z());
      }
    }
    created_faces.clear();
  }
  void after_subface_created(face_descriptor fi, TriangleMesh& tm) { created_faces.push_back(fi); }
#endif // if CGAL_VERSION_NR >= CGAL_VERSION_NUMBER(5, 4, 0)
};

#endif// FAST_CSG_KERNEL_IS_LAZY

template <typename K>
bool corefineAndComputeUnion(
  CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
{
  auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
  auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
  if (exactCallback && !remesh) {
    auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
    return PMP::corefine_and_compute_union(lhs, rhs, out, param, param);
  } else if (remesh) {
    CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
    auto param = PMP::parameters::visitor(visitor);
    auto result = PMP::corefine_and_compute_union(lhs, rhs, out, param, param, param);
    visitor.remeshSplitFaces(out);
    return result;
  } else {
    return PMP::corefine_and_compute_union(lhs, rhs, out);
  }
}

template <typename K>
bool corefineAndComputeIntersection(
  CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
{
  auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
  auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
  if (exactCallback && !remesh) {
    auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
    return PMP::corefine_and_compute_intersection(lhs, rhs, out, param, param);
  } else if (remesh) {
    CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
    auto param = PMP::parameters::visitor(visitor);
    auto result = PMP::corefine_and_compute_intersection(lhs, rhs, out, param, param, param);
    visitor.remeshSplitFaces(out);
    return result;
  } else {
    return PMP::corefine_and_compute_intersection(lhs, rhs, out);
  }
}

template <typename K>
bool corefineAndComputeDifference(
  CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
  CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
{
  auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
  auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
  if (exactCallback && !remesh) {
    auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
    return PMP::corefine_and_compute_difference(lhs, rhs, out, param, param);
  } else if (remesh) {
    CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
    auto param = PMP::parameters::visitor(visitor);
    auto result = PMP::corefine_and_compute_difference(lhs, rhs, out, param, param, param);
    visitor.remeshSplitFaces(out);
    return result;
  } else {
    return PMP::corefine_and_compute_difference(lhs, rhs, out);
  }
}

template bool corefineAndComputeUnion(
  CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);
template bool corefineAndComputeIntersection(
  CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);
template bool corefineAndComputeDifference(
  CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);

} // namespace CGALUtils


## cgalutils-mesh-edits.h
// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
#pragma once

#include <iterator>
#include <list>
#include <unordered_set>
#include <CGAL/Kernel/global_functions.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include "grid.h"

template <typename T, typename IndicesContainer>
bool removeItemsAtSortedIndices(std::vector<T>& items, const IndicesContainer& removalIndices) {
  if (removalIndices.empty()) {
    return false;
  }
  size_t outputIndex = 0;
  auto removalIndexIt = removalIndices.begin();

  for (size_t i = 0, n = items.size(); i < n; i++) {
    if (removalIndexIt != removalIndices.end() && i == *removalIndexIt) {
      ++removalIndexIt;
      continue;
    }
    if (i != outputIndex) {
      items[outputIndex] = items[i];
    }
    outputIndex++;
  }
  items.resize(outputIndex);
  return true;
}

namespace CGALUtils {

namespace PMP = CGAL::Polygon_mesh_processing;

/*! Buffer of changes to be applied to a triangle mesh. */
template <typename TriangleMesh>
class TriangleMeshEdits
{

private:
  typedef boost::graph_traits<TriangleMesh> GT;
  typedef typename GT::face_descriptor face_descriptor;
  typedef typename GT::halfedge_descriptor halfedge_descriptor;
  typedef typename GT::edge_descriptor edge_descriptor;
  typedef typename GT::vertex_descriptor vertex_descriptor;

  std::unordered_set<face_descriptor> facesToRemove;
  std::unordered_set<vertex_descriptor> verticesToRemove;
  std::vector<std::vector<vertex_descriptor>> facesToAdd;
  std::vector<std::pair<std::vector<face_descriptor>, std::vector<vertex_descriptor>>> faceReplacements;
  std::unordered_map<vertex_descriptor, vertex_descriptor> vertexReplacements;

public:

  bool isEmpty() {
    return facesToRemove.empty() &&
        verticesToRemove.empty() &&
        facesToAdd.empty() &&
        faceReplacements.empty() &&
        vertexReplacements.empty();
  }

  void removeFace(const face_descriptor& f) {
    facesToRemove.insert(f);
  }

  void removeVertex(const vertex_descriptor& v) {
    verticesToRemove.insert(v);
  }

  void addFace(const std::vector<vertex_descriptor>& vertices) {
    facesToAdd.push_back(vertices);
  }

  template <class FaceContainer, class VertexContainer>
  void replaceFaces(const FaceContainer& originalFaces, const VertexContainer& replacementFaceVertices) {
    faceReplacements.emplace_back(
      make_pair(
        std::vector<face_descriptor>(originalFaces.begin(), originalFaces.end()),
        std::vector<vertex_descriptor>(replacementFaceVertices.begin(), replacementFaceVertices.end())));
  }

  void replaceVertex(const vertex_descriptor& original, const vertex_descriptor& replacement) {
    vertexReplacements[original] = replacement;
  }

  static bool collapsePathWithConsecutiveCollinearEdges(std::vector<vertex_descriptor>& path, const TriangleMesh& tm) {
    if (path.size() <= 3) {
      return false;
    }

    std::vector<size_t> indicesToRemove;

    const auto *p1 = &tm.point(path[0]);
    const auto *p2 = &tm.point(path[1]);
    const auto *p3 = &tm.point(path[2]);

    for (size_t i = 0, n = path.size(); i < n; i++) {
      if (CGAL::are_ordered_along_line(*p1, *p2, *p3)) {
        // p2 (at index i + 1) can be removed.
        if (indicesToRemove.empty()) {
          // Late reservation: when we start to have stuff to remove, think big and reserve the worst case.
          indicesToRemove.reserve(std::min(n - i, n - 3));
        }
        if (i == n - 1) {
          indicesToRemove.insert(indicesToRemove.begin(), 0);
        } else {
          indicesToRemove.push_back(i + 1);
        }
      }
      p1 = p2;
      p2 = p3;
      p3 = &tm.point(path[(i + 3) % n]);
    }

    if (path.size() - indicesToRemove.size() < 3) {
      return false;
    }

    return removeItemsAtSortedIndices(path, indicesToRemove);
  }

  /*! Mutating in place is tricky, to say the least, so this creates a new mesh
   * and overwrites the original to it at the end for now. */
  void apply(TriangleMesh& src) const
  {
//     {
//       TriangleMeshEdits<TriangleMesh> edits;
//       edits.facesToRemove = facesToRemove;
//       edits.verticesToRemove = verticesToRemove;
//       edits.vertexReplacements = vertexReplacements;
//       doApply(src);
//     }
//     {
//       TriangleMeshEdits<TriangleMesh> edits;
//       edits.facesToAdd = facesToAdd;
//       edits.vertexReplacements = vertexReplacements;
//       doApply(src);
//     }
//   }

// private:

//   void doApply(TriangleMesh& src)
//   {
    auto printVertex = [&](auto v) {
        auto p = src.point(v);
        //std::cerr << p.x() <<  " -> ";
        printf("[%g, %g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()));
        // printf("[%g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()));
      };

    auto printVertices = [&](auto& path) {
        for (auto v : path) {
          printVertex(v);
          printf(", ");
        }
        printf("\n");
      };

    TriangleMesh copy;

    auto edgesAdded = 0;
    for (auto& vs : facesToAdd) edgesAdded += vs.size();

    auto projectedVertexCount = src.number_of_vertices() - verticesToRemove.size();
    auto projectedHalfedgeCount = src.number_of_halfedges() + edgesAdded * 2; // This is crude
    auto projectedFaceCount = src.number_of_faces() - facesToRemove.size() + facesToAdd.size();
    copy.reserve(copy.number_of_vertices() + projectedVertexCount,
                 copy.number_of_halfedges() + projectedHalfedgeCount,
                 copy.number_of_faces() + projectedFaceCount);

    std::unordered_map<vertex_descriptor, vertex_descriptor> vertexMap;
    vertexMap.reserve(projectedVertexCount);

    auto getDestinationVertex = [&](auto srcVertex) {
        auto repIt = vertexReplacements.find(srcVertex);
        if (repIt != vertexReplacements.end()) {
          srcVertex = repIt->second;
        }
        auto it = vertexMap.find(srcVertex);
        if (it == vertexMap.end()) {
          auto v = copy.add_vertex(src.point(srcVertex));
          vertexMap[srcVertex] = v;
          return v;
        }
        return it->second;
      };

    std::vector<vertex_descriptor> polygon;


    std::unordered_set<face_descriptor> facesBeingReplaced;
    for (auto &rep : faceReplacements) {
      std::copy(rep.first.begin(), rep.first.end(), inserter(facesBeingReplaced));
    }

    auto copyFace = [&](auto &f) {
      polygon.clear();

      CGAL::Vertex_around_face_iterator<TriangleMesh> vit, vend;
      for (boost::tie(vit, vend) = vertices_around_face(src.halfedge(f), src); vit != vend; ++vit) {
        auto v = *vit;
        polygon.push_back(getDestinationVertex(v));
      }

      auto face = copy.add_face(polygon);
      if (face.is_valid()) {
        if (polygon.size() > 3) {
          PMP::triangulate_face(face, copy);
        }
        // std::cerr << "Succesfully added face with " << polygon.size() << " vertices:\n";
        return true;

      } else {
        std::cerr << "Failed to add face with " << polygon.size() << " vertices:\n";
        printVertices(polygon);
        return false;
      }
    };

    for (auto f : src.faces()) {
      if (src.is_removed(f)) {
        continue;
      }
      if (facesToRemove.find(f) != facesToRemove.end()) {
        continue;
      }
      if (facesBeingReplaced.find(f) != facesBeingReplaced.end()) {
        continue;
      }
      copyFace(f);
    }

    for (auto& originalPolygon : facesToAdd) {
      polygon.clear();

      for (auto v : originalPolygon) {
        polygon.push_back(getDestinationVertex(v));
      }

      auto face = copy.add_face(polygon);
      if (face.is_valid()) {
        if (polygon.size() > 3) {
          PMP::triangulate_face(face, copy);
        }
      } else {
        std::cerr << "Failed to add face with " << polygon.size() << " vertices:\n";
        printVertices(polygon);
      }
    }

    for (auto &rep : faceReplacements) {
      polygon.clear();

      auto &originalFaces = rep.first;
      auto &rawReplacementPolygon = rep.second;

      for (auto v : rawReplacementPolygon) {
        polygon.push_back(getDestinationVertex(v));
      }

      face_descriptor face;
      try {
        face = copy.add_face(polygon);
      } catch (const CGAL::Failure_exception& e) {
        LOG(message_group::Warning, Location::NONE, "", "[fast-csg] Remesh add_face error: %1$s\n", e.what());
      }

      if (face.is_valid()) {
        if (polygon.size() > 3) {
          PMP::triangulate_face(face, copy);
        }
      } else {

#ifdef DEBUG
        auto isValid = CGAL::is_valid_polygon_mesh(src);
        try {
          face = copy.add_face(polygon);
        } catch (const CGAL::Failure_exception& e) {
          LOG(message_group::Warning, Location::NONE, "", "[fast-csg] Remesh add_face error: %1$s\n", e.what());
        }
#endif // DEBUG

        std::cerr << "Failed to replace " << originalFaces.size() << " faces with a " << polygon.size() << " vertices poly. Keeping the originals\n";
        printVertices(polygon);

        for (auto &f : originalFaces) {
          copyFace(f);
        }
      }
    }
    // std::vector<std::pair<std::vector<face_descriptor>, std::vector<vertex_descriptor>>> faceReplacements;

    src = copy;
  }
};

} // namespace CGALUtils
	// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
	#pragma once

	#include <CGAL/Polygon_mesh_processing/corefinement.h>
	#include "cgalutils-coplanar-faces-remesher.h"
	#include "cgalutils-mesh-edits.h"

	namespace CGALUtils {

	namespace internal {

	template <typename Pt>
	void forceExact(Pt& point)
	{
	CGAL::exact(point.x());
	CGAL::exact(point.y());
	CGAL::exact(point.z());
	}

	/*! Implements corefinement callbacks to support remeshing coplanar faces and forcing new numbers to
	* be exact. */
	template <typename TriangleMesh>
	struct CorefinementVisitorDelegate_ { // : public PMP::Corefinement::Default_visitor<TriangleMesh> {
	private:
	typedef boost::graph_traits<TriangleMesh> GT;
	typedef typename GT::face_descriptor face_descriptor;
	typedef typename GT::halfedge_descriptor halfedge_descriptor;
	typedef typename GT::edge_descriptor edge_descriptor;
	typedef typename GT::vertex_descriptor vertex_descriptor;

	typedef size_t PatchId; // Starting from 1

	bool forceNewLazyNumbersToExact_;
	TriangleMesh mesh1_, mesh2_, *meshOut_;

	// The following fields are scoped by mesh ([mesh1_, mesh2_, meshOut_]).
	// If meshOut_ == mesh1_ (or mesh2_), the third element of each of these
	// fields will go unused, as all accesses are using the index given by indexOf

	std::unordered_map<face_descriptor, PatchId> faceToPatchId_[3];
	std::unordered_set<PatchId> splitPatchIds_;

	face_descriptor faceBeingSplit_[3];
	std::vector<face_descriptor> facesBeingCreated_[3];

	PatchId getPatchId(const face_descriptor& f, const TriangleMesh& tm) {
	auto id = faceToPatchId_[indexOf(tm)][f];
	CGAL_assertion(id);
	return id;
	}

	void setPatchId(const face_descriptor& f, const TriangleMesh& tm, PatchId id) {
	CGAL_assertion(id);
	faceToPatchId_[indexOf(tm)][f] = id;
	}

	size_t indexOf(const TriangleMesh& tm) {
	if (&tm == mesh1_) return 0;
	if (&tm == mesh2_) return 1;
	if (&tm == meshOut_) return 2;
	CGAL_assertion(false && "Unknown mesh");
	}

	public:
	CorefinementVisitorDelegate_(TriangleMesh& m1, TriangleMesh& m2, TriangleMesh& mout,
	bool forceNewLazyNumbersToExact)
	: mesh1_(&m1), mesh2_(&m2), meshOut_(&mout), forceNewLazyNumbersToExact_(forceNewLazyNumbersToExact)
	{
	TriangleMesh *meshes[3] = {&m1, &m2, &mout};

	// Give each existing face its own global patch id.
	PatchId nextPatchId = 1;
	for (auto mi = 0; mi < 3; mi++) {
	auto& m = *meshes[mi];

	if (mi != indexOf(m)) {
	// &mout is often == &m1, so don't reassign its patch ids!
	continue;
	}

	for (auto f : m.faces()) setPatchId(f, m, nextPatchId++);
	}
	}

	void remeshSplitFaces(TriangleMesh& tm)
	{
	auto mi = indexOf(tm);
	CoplanarFacesRemesher<TriangleMesh> remesher(tm);
	remesher.remesh(faceToPatchId_[mi], splitPatchIds_);
	}

	void after_face_copy(face_descriptor f_old, const TriangleMesh& m_old, face_descriptor f_new,
	TriangleMesh& m_new)
	{
	setPatchId(f_new, m_new, getPatchId(f_old, m_old));
	}

	void before_subface_creations(face_descriptor f_split, TriangleMesh& tm)
	{
	auto mi = indexOf(tm);
	faceBeingSplit_[mi] = f_split;
	facesBeingCreated_[mi].clear();
	}

	void after_subface_created(face_descriptor fi, TriangleMesh& tm)
	{
	auto mi = indexOf(tm);
	auto id = getPatchId(faceBeingSplit_[mi], tm);
	setPatchId(fi, tm, getPatchId(faceBeingSplit_[mi], tm));
	// splitPatchIds_.insert(id); // Doing it once in after_subface_creations instead.

	#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
	if (forceNewLazyNumbersToExact_) {
	facesBeingCreated_[mi].push_back(fi);
	}
	#endif
	}

	void after_subface_creations(TriangleMesh& tm)
	{
	auto mi = indexOf(tm);
	auto id = getPatchId(faceBeingSplit_[mi], tm);
	splitPatchIds_.insert(id);

	// From CGAL 5.4 on we rely on new_vertex_added instead.
	#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
	if (forceNewLazyNumbersToExact_) {
	for (auto& fi : facesBeingCreated_[mi]) {
	CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend;
	for (boost::tie(vbegin, vend) = vertices_around_face(tm.halfedge(fi), tm); vbegin != vend;
	++vbegin) {
	forceExact(tm.point(*vbegin));
	}
	}
	}
	#endif // if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(5, 4, 0)
	}

	// Note: only called in CGAL 5.4+
	void new_vertex_added(std::size_t node_id, vertex_descriptor vh, const TriangleMesh& tm)
	{
	if (forceNewLazyNumbersToExact_) {
	forceExact(tm.point(vh));
	}
	}

	};

	} // namespace internal

	template <typename TriangleMesh>
	struct CorefinementVisitor : public PMP::Corefinement::Default_visitor<TriangleMesh> {
	typedef boost::graph_traits<TriangleMesh> GT;
	typedef typename GT::face_descriptor face_descriptor;
	typedef typename GT::vertex_descriptor vertex_descriptor;

	typedef internal::CorefinementVisitorDelegate_<TriangleMesh> Delegate;

	std::shared_ptr<Delegate> delegate_;

	CorefinementVisitor(TriangleMesh& m1, TriangleMesh& m2, TriangleMesh& mout,
	bool forceNewLazyNumbersToExact)
	: delegate_(std::make_shared<Delegate>(m1, m2, mout, forceNewLazyNumbersToExact))
	{
	}

	void remeshSplitFaces(TriangleMesh& tm) {
	delegate_->remeshSplitFaces(tm);
	}

	void before_subface_creations(face_descriptor f_old, TriangleMesh& tm) {
	delegate_->before_subface_creations(f_old, tm);
	}
	void after_subface_creations(TriangleMesh& tm) {
	delegate_->after_subface_creations(tm);
	}
	void after_subface_created(face_descriptor f_new, TriangleMesh& tm) {
	delegate_->after_subface_created(f_new, tm);
	}
	void after_face_copy(face_descriptor f_old, const TriangleMesh& old_tm, face_descriptor f_new, TriangleMesh& new_tm) {
	delegate_->after_face_copy(f_old, old_tm, f_new, new_tm);
	}
	void new_vertex_added(std::size_t node_id, vertex_descriptor vh, const TriangleMesh& tm) {
	delegate_->new_vertex_added(node_id, vh, tm);
	}
	};

	} // namespace CGALUtils
	// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
	#pragma once

	#include <iterator>
	#include <CGAL/Surface_mesh.h>
	#include <CGAL/Kernel/global_functions.h>
	#include <CGAL/Kernel/function_objects.h>
	// #include <CGAL/Polygon_mesh_processing/stitch_borders.h>

	#include "cgalutils-mesh-edits.h"

	namespace CGALUtils {

	/*! Helper that remeshes coplanar faces.
	*
	* This is useful for the results of corefinement operations, which result in
	* suboptimal meshes which growth compounds over several operations and makes
	* some models exponentially slow.
	*
	* Two faces with the same patch id MUST be coplanar, but need not be continuous
	* (technically the same patch id could result in different patches, e.g. if
	* during corefinement the same original face got split in disjoint sets of
	* descendants).
	*/
	template <typename TriangleMesh>
	struct CoplanarFacesRemesher {
	typedef boost::graph_traits<TriangleMesh> GT;
	typedef typename GT::face_descriptor face_descriptor;
	typedef typename GT::halfedge_descriptor halfedge_descriptor;
	typedef typename GT::edge_descriptor edge_descriptor;
	typedef typename GT::vertex_descriptor vertex_descriptor;

	typedef std::function<bool (const face_descriptor&)> face_predicate;

	TriangleMesh& tm;

	public:

	CoplanarFacesRemesher(TriangleMesh& mesh) : tm(mesh) {}

	template <typename PatchId>
	void remesh(
	const std::unordered_map<face_descriptor, PatchId>& faceToCoplanarPatchId,
	const std::unordered_set<PatchId>& patchesToRemesh)
	{
	if (patchesToRemesh.empty()) return;

	auto facesBefore = tm.number_of_faces();
	auto verbose = Feature::ExperimentalFastCsgDebug.is_enabled();

	// PMP::stitch_borders(tm);

	auto printVertex = [&](auto v) {
	auto p = tm.point(v);
	//std::cerr << p.x() << " -> ";
	printf("[%g, %g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()));
	// printf("[%g, %g]", CGAL::to_double(p.x()), CGAL::to_double(p.y()));
	};

	auto printVertices = [&](auto& path) {
	for (auto v : path) {
	printVertex(v);
	printf(", ");
	}
	printf("\n");
	};

	auto printHEdges = [&](auto& path) {
	for (auto he : path) {
	std::cout << he << " (";
	printVertex(tm.source(he));
	std::cout << " -> ";
	printVertex(tm.target(he));
	std::cout << "), ";
	}
	printf("\n");
	};

	auto patchToPolyhedronStr = [&](auto &patchFaces) {
	std::vector<vertex_descriptor> vertices;
	std::map<vertex_descriptor, size_t> vertexIndex;

	std::ostringstream verticesOut;
	std::ostringstream facesOut;

	facesOut << "[\n";
	for (auto& face : patchFaces) {
	CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
	facesOut << " [";
	auto first = true;
	for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
	auto he = *heIt;
	auto v = tm.source(he);
	auto it = vertexIndex.find(v);
	size_t idx;
	if (it == vertexIndex.end()) {
	idx = vertices.size();
	vertices.push_back(v);
	vertexIndex[v] = idx;
	} else {
	idx = it->second;
	}
	if (first) first = false;
	else facesOut << ", ";
	facesOut << idx;
	}
	facesOut << "],\n";
	}
	facesOut << "]";

	verticesOut << "[\n";
	for (auto v : vertices) {
	auto p = tm.point(v);
	verticesOut << "[" << CGAL::to_double(p.x()) << ", " << CGAL::to_double(p.y()) << ", " << CGAL::to_double(p.z()) << "],\n";
	}
	verticesOut << "]";

	std::ostringstream out;
	out << "polyhedron(" << verticesOut.str().c_str() << ", " << facesOut.str().c_str() << ");";
	return out.str();
	};

	auto patchBorderToPolyhedronStr = [&](auto &borderPathVertices) {
	std::vector<vertex_descriptor> vertices;
	std::map<vertex_descriptor, size_t> vertexIndex;

	std::ostringstream verticesOut;
	std::ostringstream facesOut;

	facesOut << "[[";
	auto first = true;
	for (auto& v : borderPathVertices) {
	auto it = vertexIndex.find(v);
	size_t idx;
	if (it == vertexIndex.end()) {
	idx = vertices.size();
	vertices.push_back(v);
	vertexIndex[v] = idx;
	} else {
	idx = it->second;
	}
	if (first) first = false;
	else facesOut << ", ";
	facesOut << idx;
	}
	facesOut << "]]\n";

	verticesOut << "[\n";
	for (auto v : vertices) {
	auto p = tm.point(v);
	verticesOut << "[" << CGAL::to_double(p.x()) << ", " << CGAL::to_double(p.y()) << ", " << CGAL::to_double(p.z()) << "],\n";
	}
	verticesOut << "]";

	std::ostringstream out;
	out << "polyhedron(" << verticesOut.str().c_str() << ", " << facesOut.str().c_str() << ");";
	return out.str();
	};

	try {
	TriangleMeshEdits<TriangleMesh> edits;
	std::unordered_set<PatchId> patchesProcessed;
	patchesProcessed.reserve(patchesToRemesh.size());

	std::unordered_set<face_descriptor> facesProcessed;
	// patchesProcessed.reserve(patchesToRemesh.size());

	std::set<face_descriptor> loopLocalPatchFaces;
	std::map<PatchId, bool> loopLocalIsPatchIdMap;
	std::vector<halfedge_descriptor> loopLocalBorderEdges;
	std::vector<halfedge_descriptor> loopLocalBorderPath;
	std::vector<vertex_descriptor> loopLocalBorderPathVertices;
	std::unordered_set<halfedge_descriptor> loopLocalBorderPathEdges;

	std::unordered_map<PatchId, std::string> allPatchPolyStrings;
	std::unordered_map<PatchId, std::string> allPatchReplacementsPolyStrings;

	for (auto face : tm.faces()) {
	if (tm.is_removed(face)) {
	continue;
	}
	auto idIt = faceToCoplanarPatchId.find(face);
	if (idIt == faceToCoplanarPatchId.end()) {
	continue;
	}
	const auto id = idIt->second;

	if (!getenv("COPL") && (patchesToRemesh.find(id) == patchesToRemesh.end())) {
	continue;
	}
	if (!patchesProcessed.insert(id).second) {
	continue;
	}
	if (facesProcessed.find(face) != facesProcessed.end()) {
	continue;
	}

	loopLocalIsPatchIdMap.clear();
	auto& isPatchIdMap = loopLocalIsPatchIdMap;
	isPatchIdMap[id] = true;

	loopLocalPatchFaces.clear();
	auto& patchFaces = loopLocalPatchFaces;
	patchFaces.insert(face);

	auto isHalfedgeOnBorder = [&](auto& he) {
	auto neighbourFace = tm.face(tm.opposite(he));
	if (!neighbourFace.is_valid()) {
	return true;
	}

	auto neighbourIdIt = faceToCoplanarPatchId.find(neighbourFace);
	if (neighbourIdIt == faceToCoplanarPatchId.end()) {
	std::cerr << "Unknown face, weird!";
	return true;
	}
	auto neighbourId = neighbourIdIt->second;

	auto isPatchIdIt = isPatchIdMap.find(neighbourId);
	if (isPatchIdIt != isPatchIdMap.end()) {
	auto isPatchId = isPatchIdIt->second;
	return !isPatchId;
	}

	if (getenv("COPL")) {
	auto isCoplanar = isEdgeBetweenCoplanarFaces(he);
	// std::cerr << "isCoplanar(" << face << ", " << neighbourFace << "): " << (isCoplanar ? "true" : "false") << "\n";
	isPatchIdMap[neighbourId] = isCoplanar;
	return !isCoplanar;
	}
	// We're on the border
	return true;
	};

	floodFillPatch(patchFaces, isHalfedgeOnBorder);

	allPatchPolyStrings[id] = patchToPolyhedronStr(patchFaces);

	auto conflicting = false;
	for (auto& face : patchFaces) {
	if (!facesProcessed.insert(face).second) {
	std::cerr << "Patch " << id << " trying to process already processed face " << face << "\n";
	conflicting = true;
	break;
	}
	}
	if (conflicting) {
	continue;
	}

	if (patchFaces.size() < 2) {
	continue;
	}

	auto isFaceOnPatch = [&](auto& f) {
	return patchFaces.find(f) != patchFaces.end();
	};

	loopLocalBorderEdges.clear();
	auto &borderEdges = loopLocalBorderEdges;

	for (auto& face : patchFaces) {
	CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
	for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
	auto he = *heIt;
	if (isHalfedgeOnBorder(he)) {
	borderEdges.push_back(he);
	}
	}
	}

	if (borderEdges.empty()) {
	std::cerr << "Failed to find a border edge for patch " << id << "!\n";
	continue;
	}
	halfedge_descriptor borderEdge = *borderEdges.begin();

	loopLocalBorderPath.clear();
	auto& borderPath = loopLocalBorderPath;

	if (!walkAroundPatch(borderEdge, isFaceOnPatch, borderPath)) {
	LOG(message_group::Error, Location::NONE, "",
	"[fast-csg-remesh] Failed to collect path around patch faces, invalid mesh!");
	return;
	}

	if (borderPath.size() <= 3) {
	continue;
	}

	loopLocalBorderPathVertices.clear();
	auto& borderPathVertices = loopLocalBorderPathVertices;

	loopLocalBorderPathEdges.clear();
	auto &borderPathEdges = loopLocalBorderPathEdges;

	for (auto he : borderPath) {
	borderPathEdges.insert(he);
	borderPathVertices.push_back(tm.target(he));
	}

	auto hasHoles = false;
	for (auto he : borderEdges) {
	if (borderPathEdges.find(he) == borderPathEdges.end()) {
	// Found a border halfedge that's not in the border path we've walked around the patch.
	// This means the patch is holed / has more than one border: abort!!
	hasHoles = true;
	break;
	}
	}
	if (hasHoles) {
	if (verbose) {
	LOG(message_group::None, Location::NONE, "",
	"[fast-csg-remesh] Skipping remeshing of patch with %1$lu faces as it seems to have holes.", borderPathEdges.size());
	}
	continue;
	}

	// auto veryVerbose = verbose;
	// // Only on z=1 plane to reduce clutter
	// for (auto he : borderPath) if (CGAL::to_double(tm.point(tm.target(he)).z()) != 1) veryVerbose = false;
	// if (veryVerbose) {

	// std::cerr << borderPath.size() << " edges around patch " << id << ":\n";
	// printHEdges(borderPath);
	// std::cout << "\n";

	// std::cerr << patchFaces.size() << " faces in patch " << id << ":\n";
	// for (auto f : patchFaces) {
	// std::cout << "\tFace " << f << ": ";
	// std::vector<vertex_descriptor> vertices;
	// std::vector<halfedge_descriptor> hedges;
	// CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
	// for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(f), tm); heIt != heEnd; ++heIt) {
	// auto he = *heIt;
	// vertices.push_back(tm.target(he));
	// hedges.push_back(he);
	// }
	// printVertices(vertices);
	// printHEdges(hedges);
	// std::cout << "\n";
	// }

	// std::cerr << "Patch " << id << ": removing " << patchFaces.size() << " faces\n";

	// std::cout << "Border of length " << borderPathVertices.size() << " for patch " << id << ": ";
	// printVertices(borderPathVertices);
	// }

	// Only remesh patches that have consecutive collinear edges.
	auto lengthBefore = borderPathVertices.size();
	auto collapsed = TriangleMeshEdits<TriangleMesh>::collapsePathWithConsecutiveCollinearEdges(borderPathVertices, tm);

	auto lengthAfter = borderPathVertices.size();
	if (collapsed) {
	std::cerr << "Collapsed path around patch " << id << " (" << patchFaces.size() << " faces) from " << lengthBefore << " to " << lengthAfter << " vertices\n";
	if (lengthBefore == lengthAfter) {
	std::cerr << "# WTF??\n";
	// continue;
	}
	}

	allPatchReplacementsPolyStrings[id] = patchBorderToPolyhedronStr(borderPathVertices);

	auto onlyPatchIdStr = getenv("PATCHID");
	if (onlyPatchIdStr && (atoi(onlyPatchIdStr) != id)) {
	std::cerr << "Skipping patch " << id << " (PATCHID = " << atoi(onlyPatchIdStr) << ")\n";
	continue;
	}

	// Cover the patch with a polygon. It will be triangulated when the
	// edits are applied.

	// edits.replaceFaces(patchFaces, borderPathVertices);

	for (auto& face : patchFaces) edits.removeFace(face);
	edits.addFace(borderPathVertices);

	}

	{
	std::ofstream fout("patches.scad");
	fout << "before=true;\npatchIndex=-1;\n";
	size_t patchIndex = 0;
	for (auto &p : allPatchPolyStrings) {
	auto id = p.first;
	fout << "// Patch id " << id << " (index " << patchIndex << "):\n";
	if (allPatchReplacementsPolyStrings[id].empty()) fout << "%";
	fout << "if (patchIndex < 0 \|\| patchIndex == " << patchIndex << ") {\n";
	fout << " if (before) { " << p.second.c_str() << " }\n";
	fout << " else { " << allPatchReplacementsPolyStrings[id] << " }\n";
	fout << "}\n";
	patchIndex++;
	}
	std::ofstream fout("patches.scad");
	}

	if (!edits.isEmpty()) {
	edits.apply(tm);
	}

	auto facesAfter = tm.number_of_faces();
	if (verbose) {
	LOG(message_group::None, Location::NONE, "",
	"[fast-csg-remesh] Remeshed from %1$lu to %2$lu faces (%3$lu % improvement)", facesBefore, facesAfter,
	(facesBefore - facesAfter) * 100 / facesBefore);
	}
	} catch (const CGAL::Assertion_exception& e) {
	LOG(message_group::Error, Location::NONE, "", "CGAL error in remeshSplitFaces: %1$s", e.what());
	}
	}

	private:

	bool isEdgeBetweenCoplanarFaces(const halfedge_descriptor& h) {
	auto& p = tm.point(tm.source(h));
	auto& q = tm.point(tm.target(h));
	auto& r = tm.point(tm.target(tm.next(h)));
	auto& s = tm.point(tm.target(tm.next(tm.opposite(h))));

	return CGAL::coplanar(p, q, r, s);
	// auto coplanarCosThreshold = getenv("COPL") ? atof(getenv("COPL")) : -1;
	// return coplanarCosThreshold < 0
	// ? CGAL::coplanar(p, q, r, s)
	// : CGAL::compare_dihedral_angle(p, q, r, s, coplanarCosThreshold) == CGAL::LARGER;
	}

	/! Expand the patch of known faces to surrounding neighbours that pass the (fast) predicate. /
	bool floodFillPatch(
	std::set<face_descriptor>& facePatch,
	const std::function<bool(const halfedge_descriptor&)>& isHalfedgeOnBorder)
	{
	// Avoid recursion as its depth would be model-dependent.
	std::vector<face_descriptor> unprocessedFaces(facePatch.begin(), facePatch.end());


	while (!unprocessedFaces.empty()) {
	auto face = unprocessedFaces.back();
	unprocessedFaces.pop_back();

	CGAL::Halfedge_around_face_iterator<TriangleMesh> heIt, heEnd;
	for (boost::tie(heIt, heEnd) = halfedges_around_face(tm.halfedge(face), tm); heIt != heEnd; ++heIt) {
	auto he = *heIt;
	if (isHalfedgeOnBorder(he)) {
	continue;
	}

	auto neighbourFace = tm.face(tm.opposite(he));
	if (!facePatch.insert(neighbourFace).second) {
	continue;
	}

	unprocessedFaces.push_back(neighbourFace);
	}
	}

	return true;
	}

	/! Returns true if the border is closed. /
	bool walkAroundPatch(halfedge_descriptor startingEdge,
	const face_predicate& isFaceOnPatch,
	std::vector<halfedge_descriptor>& borderOut)
	{
	auto currentEdge = startingEdge;
	borderOut.push_back(startingEdge);
	// std::set<vertex_descriptor> visitedVertices;

	std::unordered_set<vertex_descriptor> visitedVertices;

	while (tm.target(currentEdge) != tm.source(startingEdge)) {
	// visitedVertices.insert(tm.target(currentEdge));

	auto foundNext = false;
	CGAL::Halfedge_around_source_iterator<TriangleMesh> heIt, heEnd;
	for (boost::tie(heIt, heEnd) = halfedges_around_source(tm.target(currentEdge), tm); heIt != heEnd; ++heIt) {
	auto he = *heIt;
	if (he == tm.opposite(currentEdge)) {
	// Don't go back to where we just came from so quickly!
	continue;
	}
	// if (visitedVertices.find(tm.target(he)) != visitedVertices.end()) {
	// continue;
	// }

	if (!isFaceOnPatch(tm.face(he))) {
	// This halfedge's face isn't on the patch, ignore it.
	continue;
	}

	if (isFaceOnPatch(tm.face(tm.opposite(he)))) {
	// Not a border edge as both faces are on the patch
	continue;
	}

	auto v = tm.target(he);
	if (!visitedVertices.insert(v).second) {
	continue;
	}

	// std::cerr << "Edge " << he << " between " << tm.face(he) << " and " << tm.face(tm.opposite(he)) << " is on the border!\n";

	// Edge is on the border of the patch.
	currentEdge = he;
	borderOut.push_back(currentEdge);
	foundNext = true;
	break;
	}
	if (!foundNext) {
	std::cerr << "Error: Failed to find the next edge to walk to along the patch border!\n";
	return false;
	}
	}

	return true;
	}
	};

	} // namespace CGALUtils
	// Portions of this file are Copyright 2021 Google LLC, and licensed under GPL2+. See COPYING.
	#include "cgalutils.h"

	#include "cgalutils-corefinement-visitor.h"

	namespace CGALUtils {

	typedef CGAL::Surface_mesh<CGAL::Point_3<CGAL_HybridKernel3>> CGAL_HybridSurfaceMesh;

	#if FAST_CSG_KERNEL_IS_LAZY

	/*! Visitor that forces exact numbers for the vertices of all the faces created during corefinement.
	*/
	template <typename TriangleMesh>
	struct ExactLazyNumbersVisitor
	: public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<TriangleMesh> {
	typedef boost::graph_traits<TriangleMesh> GT;
	typedef typename GT::face_descriptor face_descriptor;
	typedef typename GT::halfedge_descriptor halfedge_descriptor;
	typedef typename GT::vertex_descriptor vertex_descriptor;

	#if CGAL_VERSION_NR >= CGAL_VERSION_NUMBER(5, 4, 0)
	void new_vertex_added(std::size_t i_id, vertex_descriptor v, const TriangleMesh& tm) {
	auto& pt = tm.point(v);
	CGAL::exact(pt.x());
	CGAL::exact(pt.y());
	CGAL::exact(pt.z());
	}
	#else
	face_descriptor split_face;
	std::vector<face_descriptor> created_faces;

	void before_subface_creations(face_descriptor f_split, TriangleMesh& tm)
	{
	split_face = f_split;
	created_faces.clear();
	}

	void after_subface_creations(TriangleMesh& mesh)
	{
	for (auto& fi : created_faces) {
	CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend;
	for (boost::tie(vbegin, vend) = vertices_around_face(mesh.halfedge(fi), mesh); vbegin != vend;
	++vbegin) {
	auto& v = mesh.point(*vbegin);
	CGAL::exact(v.x());
	CGAL::exact(v.y());
	CGAL::exact(v.z());
	}
	}
	created_faces.clear();
	}
	void after_subface_created(face_descriptor fi, TriangleMesh& tm) { created_faces.push_back(fi); }
	#endif // if CGAL_VERSION_NR >= CGAL_VERSION_NUMBER(5, 4, 0)
	};

	#endif// FAST_CSG_KERNEL_IS_LAZY

	template <typename K>
	bool corefineAndComputeUnion(
	CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
	{
	auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
	auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
	if (exactCallback && !remesh) {
	auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
	return PMP::corefine_and_compute_union(lhs, rhs, out, param, param);
	} else if (remesh) {
	CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
	auto param = PMP::parameters::visitor(visitor);
	auto result = PMP::corefine_and_compute_union(lhs, rhs, out, param, param, param);
	visitor.remeshSplitFaces(out);
	return result;
	} else {
	return PMP::corefine_and_compute_union(lhs, rhs, out);
	}
	}

	template <typename K>
	bool corefineAndComputeIntersection(
	CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
	{
	auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
	auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
	if (exactCallback && !remesh) {
	auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
	return PMP::corefine_and_compute_intersection(lhs, rhs, out, param, param);
	} else if (remesh) {
	CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
	auto param = PMP::parameters::visitor(visitor);
	auto result = PMP::corefine_and_compute_intersection(lhs, rhs, out, param, param, param);
	visitor.remeshSplitFaces(out);
	return result;
	} else {
	return PMP::corefine_and_compute_intersection(lhs, rhs, out);
	}
	}

	template <typename K>
	bool corefineAndComputeDifference(
	CGAL::Surface_mesh<CGAL::Point_3<K>>& lhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& rhs,
	CGAL::Surface_mesh<CGAL::Point_3<K>>& out)
	{
	auto remesh = Feature::ExperimentalFastCsgRemesh.is_enabled();
	auto exactCallback = Feature::ExperimentalFastCsgExactCorefinementCallback.is_enabled();
	if (exactCallback && !remesh) {
	auto param = PMP::parameters::visitor(ExactLazyNumbersVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>>());
	return PMP::corefine_and_compute_difference(lhs, rhs, out, param, param);
	} else if (remesh) {
	CorefinementVisitor<CGAL::Surface_mesh<CGAL::Point_3<K>>> visitor(lhs, rhs, out, exactCallback);
	auto param = PMP::parameters::visitor(visitor);
	auto result = PMP::corefine_and_compute_difference(lhs, rhs, out, param, param, param);
	visitor.remeshSplitFaces(out);
	return result;
	} else {
	return PMP::corefine_and_compute_difference(lhs, rhs, out);
	}
	}

	template bool corefineAndComputeUnion(
	CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);
	template bool corefineAndComputeIntersection(
	CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);
	template bool corefineAndComputeDifference(
	CGAL_HybridSurfaceMesh& lhs, CGAL_HybridSurfaceMesh& rhs, CGAL_HybridSurfaceMesh& out);

	} // namespace CGALUtils