edge_op.cpp

// Copyright 2021 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <unordered_map>

#include "impl.h"
#include "parallel.h"
#include "shared.h"

namespace {
using namespace manifold;

ivec3 TriOf(int edge) {
  ivec3 triEdge;
  triEdge[0] = edge;
  triEdge[1] = NextHalfedge(triEdge[0]);
  triEdge[2] = NextHalfedge(triEdge[1]);
  return triEdge;
}

bool Is01Longest(vec2 v0, vec2 v1, vec2 v2) {
  const vec2 e[3] = {v1 - v0, v2 - v1, v0 - v2};
  double l[3];
  for (int i : {0, 1, 2}) l[i] = la::dot(e[i], e[i]);
  return l[0] > l[1] && l[0] > l[2];
}

struct DuplicateEdge {
  const Halfedge* sortedHalfedge;

  bool operator()(int edge) {
    const Halfedge& halfedge = sortedHalfedge[edge];
    const Halfedge& nextHalfedge = sortedHalfedge[edge + 1];
    return halfedge.startVert == nextHalfedge.startVert &&
           halfedge.endVert == nextHalfedge.endVert;
  }
};

struct ShortEdge {
  VecView<const Halfedge> halfedge;
  VecView<const vec3> vertPos;
  const double epsilon;
  const int firstNewVert;

  bool operator()(int edge) const {
    const Halfedge& half = halfedge[edge];
    if (half.pairedHalfedge < 0 ||
        (half.startVert < firstNewVert && half.endVert < firstNewVert))
      return false;
    // Flag short edges
    const vec3 delta = vertPos[half.endVert] - vertPos[half.startVert];
    return la::dot(delta, delta) < epsilon * epsilon;
  }
};

struct FlagEdge {
  VecView<const Halfedge> halfedge;
  VecView<const TriRef> triRef;
  const int firstNewVert;

  bool operator()(int edge) const {
    const Halfedge& half = halfedge[edge];
    if (half.pairedHalfedge < 0 || half.startVert < firstNewVert) return false;
    // Flag redundant edges - those where the startVert is surrounded by only
    // two original triangles.
    const TriRef ref0 = triRef[edge / 3];
    int current = NextHalfedge(half.pairedHalfedge);
    TriRef ref1 = triRef[current / 3];
    bool ref1Updated = !ref0.SameFace(ref1);
    while (current != edge) {
      current = NextHalfedge(halfedge[current].pairedHalfedge);
      int tri = current / 3;
      const TriRef ref = triRef[tri];
      if (!ref.SameFace(ref0) && !ref.SameFace(ref1)) {
        if (!ref1Updated) {
          ref1 = ref;
          ref1Updated = true;
        } else {
          return false;
        }
      }
    }
    return true;
  }
};

struct SwappableEdge {
  VecView<const Halfedge> halfedge;
  VecView<const vec3> vertPos;
  VecView<const vec3> triNormal;
  const double tolerance;
  const int firstNewVert;

  bool operator()(int edge) const {
    const Halfedge& half = halfedge[edge];
    if (half.pairedHalfedge < 0) return false;
    if (half.startVert < firstNewVert && half.endVert < firstNewVert &&
        halfedge[NextHalfedge(edge)].endVert < firstNewVert &&
        halfedge[NextHalfedge(half.pairedHalfedge)].endVert < firstNewVert)
      return false;

    int tri = edge / 3;
    ivec3 triEdge = TriOf(edge);
    mat2x3 projection = GetAxisAlignedProjection(triNormal[tri]);
    vec2 v[3];
    for (int i : {0, 1, 2})
      v[i] = projection * vertPos[halfedge[triEdge[i]].startVert];
    if (CCW(v[0], v[1], v[2], tolerance) > 0 || !Is01Longest(v[0], v[1], v[2]))
      return false;

    // Switch to neighbor's projection.
    edge = half.pairedHalfedge;
    tri = edge / 3;
    triEdge = TriOf(edge);
    projection = GetAxisAlignedProjection(triNormal[tri]);
    for (int i : {0, 1, 2})
      v[i] = projection * vertPos[halfedge[triEdge[i]].startVert];
    return CCW(v[0], v[1], v[2], tolerance) > 0 ||
           Is01Longest(v[0], v[1], v[2]);
  }
};

struct FlagStore {
#if MANIFOLD_PAR == 1
  tbb::combinable<std::vector<size_t>> store;
#endif
  std::vector<size_t> s;

  template <typename Pred, typename F>
  void run_seq(size_t n, Pred pred, F f) {
    for (size_t i = 0; i < n; ++i)
      if (pred(i)) s.push_back(i);
    for (size_t i : s) f(i);
    s.clear();
  }

#if MANIFOLD_PAR == 1
  template <typename Pred, typename F>
  void run_par(size_t n, Pred pred, F f) {
    // Test pred in parallel, store i into thread-local vectors when pred(i) is
    // true. After testing pred, iterate and call f over the indices in
    // ascending order by using a heap in a single thread
    auto& store = this->store;
    tbb::parallel_for(tbb::blocked_range<size_t>(0, n),
                      [&store, &pred](const auto& r) {
                        auto& local = store.local();
                        for (auto i = r.begin(); i < r.end(); ++i) {
                          if (pred(i)) local.push_back(i);
                        }
                      });

    std::vector<std::vector<size_t>> stores;
    std::vector<size_t> result;
    store.combine_each(
        [&](auto& data) { stores.emplace_back(std::move(data)); });
    std::vector<size_t> sizes;
    size_t total_size = 0;
    for (const auto& tmp : stores) {
      sizes.push_back(total_size);
      total_size += tmp.size();
    }
    result.resize(total_size);
    for_each_n(ExecutionPolicy::Seq, countAt(0), stores.size(), [&](size_t i) {
      std::copy(stores[i].begin(), stores[i].end(), result.begin() + sizes[i]);
    });
    stable_sort(autoPolicy(result.size()), result.begin(), result.end());
    for (size_t x : result) f(x);
  }
#endif

  template <typename Pred, typename F>
  void run(size_t n, Pred pred, F f) {
#if MANIFOLD_PAR == 1
    if (n > 1e5) {
      run_par(n, pred, f);
    } else
#endif
    {
      run_seq(n, pred, f);
    }
  }
};

}  // namespace

namespace manifold {

/**
 * Duplicates just enough verts to covert an even-manifold to a proper
 * 2-manifold, splitting non-manifold verts and edges with too many triangles.
 */
void Manifold::Impl::CleanupTopology() {
  if (!halfedge_.size()) return;
  DEBUG_ASSERT(IsManifold(), logicErr, "polygon mesh is not manifold!");

  // In the case of a very bad triangulation, it is possible to create pinched
  // verts. They must be removed before edge collapse.
  SplitPinchedVerts();
  DedupeEdges();
}

/**
 * Collapses degenerate triangles by removing edges shorter than tolerance_ and
 * any edge that is preceeded by an edge that joins the same two face relations.
 * It also performs edge swaps on the long edges of degenerate triangles, though
 * there are some configurations of degenerates that cannot be removed this way.
 *
 * Before collapsing edges, the mesh is checked for duplicate edges (more than
 * one pair of triangles sharing the same edge), which are removed by
 * duplicating one vert and adding two triangles. These degenerate triangles are
 * likely to be collapsed again in the subsequent simplification.
 *
 * Note when an edge collapse would result in something non-manifold, the
 * vertices are duplicated in such a way as to remove handles or separate
 * meshes, thus decreasing the Genus(). It only increases when meshes that have
 * collapsed to just a pair of triangles are removed entirely.
 *
 * Verts with index less than firstNewVert will be left uncollapsed. This is
 * zero by default so that everything can be collapsed.
 *
 * Rather than actually removing the edges, this step merely marks them for
 * removal, by setting vertPos to NaN and halfedge to {-1, -1, -1, -1}.
 */
void Manifold::Impl::SimplifyTopology(int firstNewVert) {
  if (!halfedge_.size()) return;

  CleanupTopology();
  CollapseShortEdges(firstNewVert);
  CollapseColinearEdges(firstNewVert);
  SwapDegenerates(firstNewVert);
}

void Manifold::Impl::RemoveDegenerates(int firstNewVert) {
  if (!halfedge_.size()) return;

  CleanupTopology();
  CollapseShortEdges(firstNewVert);
  SwapDegenerates(firstNewVert);
}

void Manifold::Impl::CollapseShortEdges(int firstNewVert) {
  ZoneScopedN("CollapseShortEdge");
  FlagStore s;
  size_t numFlagged = 0;
  const size_t nbEdges = halfedge_.size();

  std::vector<int> scratchBuffer;
  scratchBuffer.reserve(10);
  // Short edges get to skip several checks and hence remove more classes of
  // degenerate triangles than flagged edges do, but this could in theory lead
  // to error stacking where a vertex moves too far. For this reason this is
  // restricted to epsilon, rather than tolerance.
  ShortEdge se{halfedge_, vertPos_, epsilon_, firstNewVert};
  s.run(nbEdges, se, [&](size_t i) {
    const bool didCollapse = CollapseEdge(i, scratchBuffer);
    if (didCollapse) numFlagged++;
    scratchBuffer.resize(0);
  });

#ifdef MANIFOLD_DEBUG
  if (ManifoldParams().verbose > 0 && numFlagged > 0) {
    std::cout << "collapsed " << numFlagged << " short edges" << std::endl;
  }
#endif
}

void Manifold::Impl::CollapseColinearEdges(int firstNewVert) {
  FlagStore s;
  size_t numFlagged = 0;
  const size_t nbEdges = halfedge_.size();
  std::vector<int> scratchBuffer;
  scratchBuffer.reserve(10);
  while (1) {
    ZoneScopedN("CollapseFlaggedEdge");
    numFlagged = 0;
    // Collapse colinear edges, but only remove new verts, i.e. verts with
    // index
    // >= firstNewVert. This is used to keep the Boolean from changing the
    // non-intersecting parts of the input meshes. Colinear is defined not by a
    // local check, but by the global MarkCoplanar function, which keeps this
    // from being vulnerable to error stacking.
    FlagEdge se{halfedge_, meshRelation_.triRef, firstNewVert};
    s.run(nbEdges, se, [&](size_t i) {
      const bool didCollapse = CollapseEdge(i, scratchBuffer);
      if (didCollapse) numFlagged++;
      scratchBuffer.resize(0);
    });
    if (numFlagged == 0) break;

#ifdef MANIFOLD_DEBUG
    if (ManifoldParams().verbose > 0 && numFlagged > 0) {
      std::cout << "collapsed " << numFlagged << " colinear edges" << std::endl;
    }
#endif
  }
}

void Manifold::Impl::SwapDegenerates(int firstNewVert) {
  ZoneScopedN("RecursiveEdgeSwap");
  FlagStore s;
  size_t numFlagged = 0;
  const size_t nbEdges = halfedge_.size();
  std::vector<int> scratchBuffer;
  scratchBuffer.reserve(10);

  SwappableEdge se{halfedge_, vertPos_, faceNormal_, tolerance_, firstNewVert};
  std::vector<int> edgeSwapStack;
  std::vector<int> visited(halfedge_.size(), -1);
  int tag = 0;
  s.run(nbEdges, se, [&](size_t i) {
    numFlagged++;
    tag++;
    RecursiveEdgeSwap(i, tag, visited, edgeSwapStack, scratchBuffer);
    while (!edgeSwapStack.empty()) {
      int last = edgeSwapStack.back();
      edgeSwapStack.pop_back();
      RecursiveEdgeSwap(last, tag, visited, edgeSwapStack, scratchBuffer);
    }
  });

#ifdef MANIFOLD_DEBUG
  if (ManifoldParams().verbose > 0 && numFlagged > 0) {
    std::cout << "swapped " << numFlagged << " edges" << std::endl;
  }
#endif
}

// Deduplicate the given 4-manifold edge by duplicating endVert, thus making the
// edges distinct. Also duplicates startVert if it becomes pinched.
void Manifold::Impl::DedupeEdge(const int edge) {
  // Orbit endVert
  const int startVert = halfedge_[edge].startVert;
  const int endVert = halfedge_[edge].endVert;
  const int endProp = halfedge_[NextHalfedge(edge)].propVert;
  int current = halfedge_[NextHalfedge(edge)].pairedHalfedge;
  while (current != edge) {
    const int vert = halfedge_[current].startVert;
    if (vert == startVert) {
      // Single topological unit needs 2 faces added to be split
      const int newVert = vertPos_.size();
      vertPos_.push_back(vertPos_[endVert]);
      if (vertNormal_.size() > 0) vertNormal_.push_back(vertNormal_[endVert]);
      current = halfedge_[NextHalfedge(current)].pairedHalfedge;
      const int opposite = halfedge_[NextHalfedge(edge)].pairedHalfedge;

      UpdateVert(newVert, current, opposite);

      int newHalfedge = halfedge_.size();
      int oldFace = current / 3;
      int outsideVert = halfedge_[current].startVert;
      halfedge_.push_back({endVert, newVert, -1, endProp});
      halfedge_.push_back({newVert, outsideVert, -1, endProp});
      halfedge_.push_back(
          {outsideVert, endVert, -1, halfedge_[current].propVert});
      PairUp(newHalfedge + 2, halfedge_[current].pairedHalfedge);
      PairUp(newHalfedge + 1, current);
      if (meshRelation_.triRef.size() > 0)
        meshRelation_.triRef.push_back(meshRelation_.triRef[oldFace]);
      if (faceNormal_.size() > 0) faceNormal_.push_back(faceNormal_[oldFace]);

      newHalfedge += 3;
      oldFace = opposite / 3;
      outsideVert = halfedge_[opposite].startVert;
      halfedge_.push_back({newVert, endVert, -1, endProp});  // fix prop
      halfedge_.push_back({endVert, outsideVert, -1, endProp});
      halfedge_.push_back(
          {outsideVert, newVert, -1, halfedge_[opposite].propVert});
      PairUp(newHalfedge + 2, halfedge_[opposite].pairedHalfedge);
      PairUp(newHalfedge + 1, opposite);
      PairUp(newHalfedge, newHalfedge - 3);
      if (meshRelation_.triRef.size() > 0)
        meshRelation_.triRef.push_back(meshRelation_.triRef[oldFace]);
      if (faceNormal_.size() > 0) faceNormal_.push_back(faceNormal_[oldFace]);

      break;
    }

    current = halfedge_[NextHalfedge(current)].pairedHalfedge;
  }

  if (current == edge) {
    // Separate topological unit needs no new faces to be split
    const int newVert = vertPos_.size();
    vertPos_.push_back(vertPos_[endVert]);
    if (vertNormal_.size() > 0) vertNormal_.push_back(vertNormal_[endVert]);

    ForVert(NextHalfedge(current), [this, newVert](int e) {
      halfedge_[e].startVert = newVert;
      halfedge_[halfedge_[e].pairedHalfedge].endVert = newVert;
    });
  }

  // Orbit startVert
  const int pair = halfedge_[edge].pairedHalfedge;
  current = halfedge_[NextHalfedge(pair)].pairedHalfedge;
  while (current != pair) {
    const int vert = halfedge_[current].startVert;
    if (vert == endVert) {
      break;  // Connected: not a pinched vert
    }
    current = halfedge_[NextHalfedge(current)].pairedHalfedge;
  }

  if (current == pair) {
    // Split the pinched vert the previous split created.
    const int newVert = vertPos_.size();
    vertPos_.push_back(vertPos_[endVert]);
    if (vertNormal_.size() > 0) vertNormal_.push_back(vertNormal_[endVert]);

    ForVert(NextHalfedge(current), [this, newVert](int e) {
      halfedge_[e].startVert = newVert;
      halfedge_[halfedge_[e].pairedHalfedge].endVert = newVert;
    });
  }
}

void Manifold::Impl::PairUp(int edge0, int edge1) {
  halfedge_[edge0].pairedHalfedge = edge1;
  halfedge_[edge1].pairedHalfedge = edge0;
}

// Traverses CW around startEdge.endVert from startEdge to endEdge
// (edgeEdge.endVert must == startEdge.endVert), updating each edge to point
// to vert instead.
void Manifold::Impl::UpdateVert(int vert, int startEdge, int endEdge) {
  int current = startEdge;
  while (current != endEdge) {
    halfedge_[current].endVert = vert;
    current = NextHalfedge(current);
    halfedge_[current].startVert = vert;
    current = halfedge_[current].pairedHalfedge;
    DEBUG_ASSERT(current != startEdge, logicErr, "infinite loop in decimator!");
  }
}

// In the event that the edge collapse would create a non-manifold edge,
// instead we duplicate the two verts and attach the manifolds the other way
// across this edge.
void Manifold::Impl::FormLoop(int current, int end) {
  int startVert = vertPos_.size();
  vertPos_.push_back(vertPos_[halfedge_[current].startVert]);
  int endVert = vertPos_.size();
  vertPos_.push_back(vertPos_[halfedge_[current].endVert]);

  int oldMatch = halfedge_[current].pairedHalfedge;
  int newMatch = halfedge_[end].pairedHalfedge;

  UpdateVert(startVert, oldMatch, newMatch);
  UpdateVert(endVert, end, current);

  halfedge_[current].pairedHalfedge = newMatch;
  halfedge_[newMatch].pairedHalfedge = current;
  halfedge_[end].pairedHalfedge = oldMatch;
  halfedge_[oldMatch].pairedHalfedge = end;

  RemoveIfFolded(end);
}

void Manifold::Impl::CollapseTri(const ivec3& triEdge) {
  if (halfedge_[triEdge[1]].pairedHalfedge == -1) return;
  int pair1 = halfedge_[triEdge[1]].pairedHalfedge;
  int pair2 = halfedge_[triEdge[2]].pairedHalfedge;
  halfedge_[pair1].pairedHalfedge = pair2;
  halfedge_[pair2].pairedHalfedge = pair1;
  for (int i : {0, 1, 2}) {
    halfedge_[triEdge[i]] = {-1, -1, -1, halfedge_[triEdge[i]].propVert};
  }
}

void Manifold::Impl::RemoveIfFolded(int edge) {
  const ivec3 tri0edge = TriOf(edge);
  const ivec3 tri1edge = TriOf(halfedge_[edge].pairedHalfedge);
  if (halfedge_[tri0edge[1]].pairedHalfedge == -1) return;
  if (halfedge_[tri0edge[1]].endVert == halfedge_[tri1edge[1]].endVert) {
    if (halfedge_[tri0edge[1]].pairedHalfedge == tri1edge[2]) {
      if (halfedge_[tri0edge[2]].pairedHalfedge == tri1edge[1]) {
        for (int i : {0, 1, 2})
          vertPos_[halfedge_[tri0edge[i]].startVert] = vec3(NAN);
      } else {
        vertPos_[halfedge_[tri0edge[1]].startVert] = vec3(NAN);
      }
    } else {
      if (halfedge_[tri0edge[2]].pairedHalfedge == tri1edge[1]) {
        vertPos_[halfedge_[tri1edge[1]].startVert] = vec3(NAN);
      }
    }
    PairUp(halfedge_[tri0edge[1]].pairedHalfedge,
           halfedge_[tri1edge[2]].pairedHalfedge);
    PairUp(halfedge_[tri0edge[2]].pairedHalfedge,
           halfedge_[tri1edge[1]].pairedHalfedge);
    for (int i : {0, 1, 2}) {
      halfedge_[tri0edge[i]] = {-1, -1, -1};
      halfedge_[tri1edge[i]] = {-1, -1, -1};
    }
  }
}

// Collapses the given edge by removing startVert - returns false if the edge
// cannot be collapsed. May split the mesh topologically if the collapse would
// have resulted in a 4-manifold edge. Do not collapse an edge if startVert is
// pinched - the vert would be marked NaN, but other edges could still be
// pointing to it.
bool Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
  Vec<TriRef>& triRef = meshRelation_.triRef;

  const Halfedge toRemove = halfedge_[edge];
  if (toRemove.pairedHalfedge < 0) return false;

  const int endVert = toRemove.endVert;
  const ivec3 tri0edge = TriOf(edge);
  const ivec3 tri1edge = TriOf(toRemove.pairedHalfedge);

  const vec3 pNew = vertPos_[endVert];
  const vec3 pOld = vertPos_[toRemove.startVert];
  const vec3 delta = pNew - pOld;
  const bool shortEdge = la::dot(delta, delta) < epsilon_ * epsilon_;

  // Orbit startVert
  int start = halfedge_[tri1edge[1]].pairedHalfedge;
  int current = tri1edge[2];
  if (!shortEdge) {
    current = start;
    TriRef refCheck = triRef[toRemove.pairedHalfedge / 3];
    vec3 pLast = vertPos_[halfedge_[tri1edge[1]].endVert];
    while (current != tri1edge[0]) {
      current = NextHalfedge(current);
      vec3 pNext = vertPos_[halfedge_[current].endVert];
      const int tri = current / 3;
      const TriRef ref = triRef[tri];
      const mat2x3 projection = GetAxisAlignedProjection(faceNormal_[tri]);
      // Don't collapse if the edge is not redundant (this may have changed due
      // to the collapse of neighbors).
      if (!ref.SameFace(refCheck)) {
        const TriRef oldRef = refCheck;
        refCheck = triRef[edge / 3];
        if (!ref.SameFace(refCheck)) {
          return false;
        }
        if (ref.meshID != oldRef.meshID || ref.faceID != oldRef.faceID ||
            la::dot(faceNormal_[toRemove.pairedHalfedge / 3],
                    faceNormal_[tri]) < -0.5) {
          // Restrict collapse to colinear edges when the edge separates faces
          // or the edge is sharp. This ensures large shifts are not introduced
          // parallel to the tangent plane.
          if (CCW(projection * pLast, projection * pOld, projection * pNew,
                  epsilon_) != 0)
            return false;
        }
      }

      // Don't collapse edge if it would cause a triangle to invert.
      if (CCW(projection * pNext, projection * pLast, projection * pNew,
              epsilon_) < 0)
        return false;

      pLast = pNext;
      current = halfedge_[current].pairedHalfedge;
    }
  }

  // Orbit endVert
  {
    int current = halfedge_[tri0edge[1]].pairedHalfedge;
    while (current != tri1edge[2]) {
      current = NextHalfedge(current);
      edges.push_back(current);
      current = halfedge_[current].pairedHalfedge;
    }
  }

  // Remove toRemove.startVert and replace with endVert.
  vertPos_[toRemove.startVert] = vec3(NAN);
  CollapseTri(tri1edge);

  // Orbit startVert
  const int tri0 = edge / 3;
  const int tri1 = toRemove.pairedHalfedge / 3;
  current = start;
  while (current != tri0edge[2]) {
    current = NextHalfedge(current);

    if (NumProp() > 0) {
      // Update the shifted triangles to the vertBary of endVert
      const int tri = current / 3;
      if (triRef[tri].SameFace(triRef[tri0])) {
        halfedge_[current].propVert = halfedge_[NextHalfedge(edge)].propVert;
      } else if (triRef[tri].SameFace(triRef[tri1])) {
        halfedge_[current].propVert =
            halfedge_[toRemove.pairedHalfedge].propVert;
      }
    }

    const int vert = halfedge_[current].endVert;
    const int next = halfedge_[current].pairedHalfedge;
    for (size_t i = 0; i < edges.size(); ++i) {
      if (vert == halfedge_[edges[i]].endVert) {
        FormLoop(edges[i], current);
        start = next;
        edges.resize(i);
        break;
      }
    }
    current = next;
  }

  UpdateVert(endVert, start, tri0edge[2]);
  CollapseTri(tri0edge);
  RemoveIfFolded(start);
  return true;
}

void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag,
                                       std::vector<int>& visited,
                                       std::vector<int>& edgeSwapStack,
                                       std::vector<int>& edges) {
  Vec<TriRef>& triRef = meshRelation_.triRef;

  if (edge < 0) return;
  const int pair = halfedge_[edge].pairedHalfedge;
  if (pair < 0) return;

  // avoid infinite recursion
  if (visited[edge] == tag && visited[pair] == tag) return;

  const ivec3 tri0edge = TriOf(edge);
  const ivec3 tri1edge = TriOf(pair);

  mat2x3 projection = GetAxisAlignedProjection(faceNormal_[edge / 3]);
  vec2 v[4];
  for (int i : {0, 1, 2})
    v[i] = projection * vertPos_[halfedge_[tri0edge[i]].startVert];
  // Only operate on the long edge of a degenerate triangle.
  if (CCW(v[0], v[1], v[2], tolerance_) > 0 || !Is01Longest(v[0], v[1], v[2]))
    return;

  // Switch to neighbor's projection.
  projection = GetAxisAlignedProjection(faceNormal_[pair / 3]);
  for (int i : {0, 1, 2})
    v[i] = projection * vertPos_[halfedge_[tri0edge[i]].startVert];
  v[3] = projection * vertPos_[halfedge_[tri1edge[2]].startVert];

  auto SwapEdge = [&]() {
    // The 0-verts are swapped to the opposite 2-verts.
    const int v0 = halfedge_[tri0edge[2]].startVert;
    const int v1 = halfedge_[tri1edge[2]].startVert;
    halfedge_[tri0edge[0]].startVert = v1;
    halfedge_[tri0edge[2]].endVert = v1;
    halfedge_[tri1edge[0]].startVert = v0;
    halfedge_[tri1edge[2]].endVert = v0;
    PairUp(tri0edge[0], halfedge_[tri1edge[2]].pairedHalfedge);
    PairUp(tri1edge[0], halfedge_[tri0edge[2]].pairedHalfedge);
    PairUp(tri0edge[2], tri1edge[2]);
    // Both triangles are now subsets of the neighboring triangle.
    const int tri0 = tri0edge[0] / 3;
    const int tri1 = tri1edge[0] / 3;
    faceNormal_[tri0] = faceNormal_[tri1];
    triRef[tri0] = triRef[tri1];
    const double l01 = la::length(v[1] - v[0]);
    const double l02 = la::length(v[2] - v[0]);
    const double a = std::max(0.0, std::min(1.0, l02 / l01));
    // Update properties if applicable
    if (properties_.size() > 0) {
      Vec<double>& prop = properties_;
      halfedge_[tri0edge[1]].propVert = halfedge_[tri1edge[0]].propVert;
      halfedge_[tri0edge[0]].propVert = halfedge_[tri1edge[2]].propVert;
      halfedge_[tri0edge[2]].propVert = halfedge_[tri1edge[2]].propVert;
      const int numProp = NumProp();
      const int newProp = prop.size() / numProp;
      const int propIdx0 = halfedge_[tri1edge[0]].propVert;
      const int propIdx1 = halfedge_[tri1edge[1]].propVert;
      for (int p = 0; p < numProp; ++p) {
        prop.push_back(a * prop[numProp * propIdx0 + p] +
                       (1 - a) * prop[numProp * propIdx1 + p]);
      }
      halfedge_[tri1edge[0]].propVert = newProp;
      halfedge_[tri0edge[2]].propVert = newProp;
    }

    // if the new edge already exists, duplicate the verts and split the mesh.
    int current = halfedge_[tri1edge[0]].pairedHalfedge;
    const int endVert = halfedge_[tri1edge[1]].endVert;
    while (current != tri0edge[1]) {
      current = NextHalfedge(current);
      if (halfedge_[current].endVert == endVert) {
        FormLoop(tri0edge[2], current);
        RemoveIfFolded(tri0edge[2]);
        return;
      }
      current = halfedge_[current].pairedHalfedge;
    }
  };

  // Only operate if the other triangles are not degenerate.
  if (CCW(v[1], v[0], v[3], tolerance_) <= 0) {
    if (!Is01Longest(v[1], v[0], v[3])) return;
    // Two facing, long-edge degenerates can swap.
    SwapEdge();
    const vec2 e23 = v[3] - v[2];
    if (la::dot(e23, e23) < tolerance_ * tolerance_) {
      tag++;
      CollapseEdge(tri0edge[2], edges);
      edges.resize(0);
    } else {
      visited[edge] = tag;
      visited[pair] = tag;
      edgeSwapStack.insert(edgeSwapStack.end(), {tri1edge[1], tri1edge[0],
                                                 tri0edge[1], tri0edge[0]});
    }
    return;
  } else if (CCW(v[0], v[3], v[2], tolerance_) <= 0 ||
             CCW(v[1], v[2], v[3], tolerance_) <= 0) {
    return;
  }
  // Normal path
  SwapEdge();
  visited[edge] = tag;
  visited[pair] = tag;
  edgeSwapStack.insert(edgeSwapStack.end(),
                       {halfedge_[tri1edge[0]].pairedHalfedge,
                        halfedge_[tri0edge[1]].pairedHalfedge});
}

void Manifold::Impl::SplitPinchedVerts() {
  ZoneScoped;

  auto nbEdges = halfedge_.size();
#if MANIFOLD_PAR == 1
  if (nbEdges > 1e4) {
    std::mutex mutex;
    std::vector<size_t> pinched;
    // This parallelized version is non-trivial so we can't reuse the code
    //
    // The idea here is to identify cycles of halfedges that can be iterated
    // through using ForVert. Pinched verts are vertices where there are
    // multiple cycles associated with the vertex. Each cycle is identified with
    // the largest halfedge index within the cycle, and when there are multiple
    // cycles associated with the same starting vertex but with different ids,
    // it means we have a pinched vertex. This check is done by using a single
    // atomic cas operation, the expected case is either invalid id (the vertex
    // was not processed) or with the same id.
    //
    // The local store is to store the processed halfedges, so to avoid
    // repetitive processing. Note that it only approximates the processed
    // halfedges because it is thread local. This is why we need a vector to
    // deduplicate the probematic halfedges we found.
    std::vector<std::atomic<size_t>> largestEdge(NumVert());
    for_each(ExecutionPolicy::Par, countAt(0), countAt(NumVert()),
             [&largestEdge](size_t i) {
               largestEdge[i].store(std::numeric_limits<size_t>::max());
             });
    tbb::combinable<std::vector<bool>> store(
        [nbEdges]() { return std::vector<bool>(nbEdges, false); });
    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, nbEdges),
        [&store, &mutex, &pinched, &largestEdge, this](const auto& r) {
          auto& local = store.local();
          std::vector<size_t> pinchedLocal;
          for (auto i = r.begin(); i < r.end(); ++i) {
            if (local[i]) continue;
            local[i] = true;
            const int vert = halfedge_[i].startVert;
            if (vert == -1) continue;
            size_t largest = i;
            ForVert(i, [&local, &largest](int current) {
              local[current] = true;
              largest = std::max(largest, static_cast<size_t>(current));
            });
            size_t expected = std::numeric_limits<size_t>::max();
            if (!largestEdge[vert].compare_exchange_strong(expected, largest) &&
                expected != largest) {
              // we know that there is another loop...
              pinchedLocal.push_back(largest);
            }
          }
          if (!pinchedLocal.empty()) {
            std::lock_guard<std::mutex> lock(mutex);
            pinched.insert(pinched.end(), pinchedLocal.begin(),
                           pinchedLocal.end());
          }
        });

    manifold::stable_sort(pinched.begin(), pinched.end());
    std::vector<bool> halfedgeProcessed(nbEdges, false);
    for (size_t i : pinched) {
      if (halfedgeProcessed[i]) continue;
      vertPos_.push_back(vertPos_[halfedge_[i].startVert]);
      const int vert = NumVert() - 1;
      ForVert(i, [this, vert, &halfedgeProcessed](int current) {
        halfedgeProcessed[current] = true;
        halfedge_[current].startVert = vert;
        halfedge_[halfedge_[current].pairedHalfedge].endVert = vert;
      });
    }
  } else
#endif
  {
    std::vector<bool> vertProcessed(NumVert(), false);
    std::vector<bool> halfedgeProcessed(nbEdges, false);
    for (size_t i = 0; i < nbEdges; ++i) {
      if (halfedgeProcessed[i]) continue;
      int vert = halfedge_[i].startVert;
      if (vert == -1) continue;
      if (vertProcessed[vert]) {
        vertPos_.push_back(vertPos_[vert]);
        vert = NumVert() - 1;
      } else {
        vertProcessed[vert] = true;
      }
      ForVert(i, [this, &halfedgeProcessed, vert](int current) {
        halfedgeProcessed[current] = true;
        halfedge_[current].startVert = vert;
        halfedge_[halfedge_[current].pairedHalfedge].endVert = vert;
      });
    }
  }
}

void Manifold::Impl::DedupeEdges() {
  while (1) {
    ZoneScopedN("DedupeEdge");

    const size_t nbEdges = halfedge_.size();
    std::vector<size_t> duplicates;
    auto localLoop = [&](size_t start, size_t end, std::vector<bool>& local,
                         std::vector<size_t>& results) {
      // Iterate over all halfedges that start with the same vertex, and check
      // for halfedges with the same ending vertex.
      // Note: we use Vec and linear search when the number of neighbor is
      // small because unordered_set requires allocations and is expensive.
      // We switch to unordered_set when the number of neighbor is
      // larger to avoid making things quadratic.
      // We do it in two pass, the first pass to find the minimal halfedges with
      // the target start and end verts, the second pass flag all the duplicated
      // halfedges that are not having the minimal index as duplicates.
      // This ensures deterministic result.
      //
      // The local store is to store the processed halfedges, so to avoid
      // repetitive processing. Note that it only approximates the processed
      // halfedges because it is thread local.
      Vec<std::pair<int, int>> endVerts;
      std::unordered_map<int, int> endVertSet;
      for (auto i = start; i < end; ++i) {
        if (local[i] || halfedge_[i].startVert == -1 ||
            halfedge_[i].endVert == -1)
          continue;
        // we want to keep the allocation
        endVerts.clear(false);
        endVertSet.clear();

        // first iteration, populate entries
        // this makes sure we always report the same set of entries
        ForVert(i, [&local, &endVerts, &endVertSet, this](int current) {
          local[current] = true;
          if (halfedge_[current].startVert == -1 ||
              halfedge_[current].endVert == -1) {
            return;
          }
          int endV = halfedge_[current].endVert;
          if (endVertSet.empty()) {
            auto iter = std::find_if(endVerts.begin(), endVerts.end(),
                                     [endV](const std::pair<int, int>& pair) {
                                       return pair.first == endV;
                                     });
            if (iter != endVerts.end()) {
              iter->second = std::min(iter->second, current);
            } else {
              endVerts.push_back({endV, current});
              if (endVerts.size() > 32) {
                endVertSet.insert(endVerts.begin(), endVerts.end());
                endVerts.clear(false);
              }
            }
          } else {
            auto pair = endVertSet.insert({endV, current});
            if (!pair.second)
              pair.first->second = std::min(pair.first->second, current);
          }
        });
        // second iteration, actually check for duplicates
        // we always report the same set of duplicates, excluding the smallest
        // halfedge in the set of duplicates
        ForVert(i, [&endVerts, &endVertSet, &results, this](int current) {
          if (halfedge_[current].startVert == -1 ||
              halfedge_[current].endVert == -1) {
            return;
          }
          int endV = halfedge_[current].endVert;
          if (endVertSet.empty()) {
            auto iter = std::find_if(endVerts.begin(), endVerts.end(),
                                     [endV](const std::pair<int, int>& pair) {
                                       return pair.first == endV;
                                     });
            if (iter->second != current) results.push_back(current);
          } else {
            auto iter = endVertSet.find(endV);
            if (iter->second != current) results.push_back(current);
          }
        });
      }
    };
#if MANIFOLD_PAR == 1
    if (nbEdges > 1e4) {
      std::mutex mutex;
      tbb::combinable<std::vector<bool>> store(
          [nbEdges]() { return std::vector<bool>(nbEdges, false); });
      tbb::parallel_for(
          tbb::blocked_range<size_t>(0, nbEdges),
          [&store, &mutex, &duplicates, &localLoop](const auto& r) {
            auto& local = store.local();
            std::vector<size_t> duplicatesLocal;
            localLoop(r.begin(), r.end(), local, duplicatesLocal);
            if (!duplicatesLocal.empty()) {
              std::lock_guard<std::mutex> lock(mutex);
              duplicates.insert(duplicates.end(), duplicatesLocal.begin(),
                                duplicatesLocal.end());
            }
          });
      manifold::stable_sort(duplicates.begin(), duplicates.end());
      duplicates.resize(
          std::distance(duplicates.begin(),
                        std::unique(duplicates.begin(), duplicates.end())));
    } else
#endif
    {
      std::vector<bool> local(nbEdges, false);
      localLoop(0, nbEdges, local, duplicates);
    }

    size_t numFlagged = 0;
    for (size_t i : duplicates) {
      DedupeEdge(i);
      numFlagged++;
    }

    if (numFlagged == 0) break;

#ifdef MANIFOLD_DEBUG
    if (ManifoldParams().verbose > 0) {
      std::cout << "found " << numFlagged << " duplicate edges to split"
                << std::endl;
    }
#endif
  }
}
}  // namespace manifold