Earclip (elalish#564)

* earclip compiles

* some tests passing

* 46 tests passing

* refactor

* cleanup

* new cost

* cleanup

* fixed TouchingHole

* more fixes

* only 5 tests failing

* fixed short sides

* fixed Sponge2

* better cost

* improved keyholing

* Moved Rect to public

* added top/bottom check

* fixed bBox

* Polygon tests pass

* fixed Sponge3a

* another fix, another test

* keyhole refactor

* Added ClipDegenerates

* Refactor FindStarts

* fixed bugs

* fixed Sponge4

* refactor clipped

* another bug, another test

* another bug, another test

* Sponge4 passing

* all tests passing

* removed old triangulator

* fix isfinite

* debug checks

* remove 2-manifold check

* added comments

* fixed missing point bug

* another test, another fix

* two tests added and fixed

* fix build

* fixed CoplanarProp

* fixed seg fault

* reenable 2-manifold check

* refactor for polygon area

* factored out loop

* tests passing

* another test, another fix

* another test, another fix

* fixed dedupe segfault

* ignore duplicated verts

* re-enabled turn180

* re-enabled duplicate tests

* another test, another fix

* refactor keyholing

* more precise area computation

* disabled a test

* added test

* another test, another test

* fixed bridge bound

* fewer short edges

* replaced intersection sorting with InsideEdge

* another test, another fix

* allow a few degenerates in Sponge4

.vscode/settings.json

@@ -144,4 +144,5 @@
     "editor.defaultFormatter": "ms-python.black-formatter"
   },
   "python.formatting.provider": "none",
+  "clang-format.executable": "clang-format",
 }

bindings/c/include/manifoldc.h

@@ -301,8 +301,6 @@ ManifoldRect *manifold_rect_mul(void *mem, ManifoldRect *r, float x, float y);
 int manifold_rect_does_overlap_rect(ManifoldRect *a, ManifoldRect *r);
 int manifold_rect_is_empty(ManifoldRect *r);
 int manifold_rect_is_finite(ManifoldRect *r);
-ManifoldCrossSection *manifold_rect_as_cross_section(void *mem,
-                                                     ManifoldRect *r);
 
 // Bounding Box
 

bindings/c/rect.cpp

@@ -98,9 +98,3 @@ int manifold_rect_does_overlap_rect(ManifoldRect *a, ManifoldRect *r) {
 int manifold_rect_is_empty(ManifoldRect *r) { return from_c(r)->IsEmpty(); }
 
 int manifold_rect_is_finite(ManifoldRect *r) { return from_c(r)->IsFinite(); }
-
-ManifoldCrossSection *manifold_rect_as_cross_section(void *mem,
-                                                     ManifoldRect *r) {
-  auto cs = from_c(r)->AsCrossSection();
-  return to_c(new (mem) CrossSection(cs));
-}

src/cross_section/include/cross_section.h

@@ -24,8 +24,6 @@
 
 namespace manifold {
 
-class Rect;
-
 /** @addtogroup Core
  *  @{
  */
@@ -172,64 +170,4 @@ class CrossSection {
   std::shared_ptr<const PathImpl> GetPaths() const;
 };
 /** @} */
-
-/** @addtogroup Connections
- *  @{
- */
-
-/**
- * Axis-aligned rectangular bounds.
- */
-class Rect {
- public:
-  glm::vec2 min = glm::vec2(0);
-  glm::vec2 max = glm::vec2(0);
-
-  /** @name Creation
-   *  Constructors
-   */
-  ///@{
-  Rect();
-  ~Rect();
-  Rect(const Rect& other);
-  Rect& operator=(const Rect& other);
-  Rect(Rect&&) noexcept;
-  Rect& operator=(Rect&&) noexcept;
-  Rect(const glm::vec2 a, const glm::vec2 b);
-  ///@}
-
-  /** @name Information
-   *  Details of the rectangle
-   */
-  ///@{
-  glm::vec2 Size() const;
-  float Area() const;
-  float Scale() const;
-  glm::vec2 Center() const;
-  bool Contains(const glm::vec2& pt) const;
-  bool Contains(const Rect& other) const;
-  bool DoesOverlap(const Rect& other) const;
-  bool IsEmpty() const;
-  bool IsFinite() const;
-  ///@}
-
-  /** @name Modification
-   */
-  ///@{
-  void Union(const glm::vec2 p);
-  Rect Union(const Rect& other) const;
-  Rect operator+(const glm::vec2 shift) const;
-  Rect& operator+=(const glm::vec2 shift);
-  Rect operator*(const glm::vec2 scale) const;
-  Rect& operator*=(const glm::vec2 scale);
-  Rect Transform(const glm::mat3x2& m) const;
-  ///@}
-
-  /** @name Conversion
-   */
-  ///@{
-  CrossSection AsCrossSection() const;
-  ///@}
-};
-/** @} */
 }  // namespace manifold

src/cross_section/src/cross_section.cpp

@@ -742,167 +742,4 @@ Polygons CrossSection::ToPolygons() const {
   }
   return polys;
 }
-
-// Rect
-
-/**
- * Default constructor is an empty rectangle..
- */
-Rect::Rect() {}
-
-Rect::~Rect() = default;
-Rect::Rect(Rect&&) noexcept = default;
-Rect& Rect::operator=(Rect&&) noexcept = default;
-Rect::Rect(const Rect& other) {
-  min = glm::vec2(other.min);
-  max = glm::vec2(other.max);
-}
-
-/**
- * Create a rectangle that contains the two given points.
- */
-Rect::Rect(const glm::vec2 a, const glm::vec2 b) {
-  min = glm::min(a, b);
-  max = glm::max(a, b);
-}
-
-/**
- * Return the dimensions of the rectangle.
- */
-glm::vec2 Rect::Size() const { return max - min; }
-
-/**
- * Return the area of the rectangle.
- */
-float Rect::Area() const {
-  auto sz = Size();
-  return sz.x * sz.y;
-}
-
-/**
- * Returns the absolute-largest coordinate value of any contained
- * point.
- */
-float Rect::Scale() const {
-  glm::vec2 absMax = glm::max(glm::abs(min), glm::abs(max));
-  return glm::max(absMax.x, absMax.y);
-}
-
-/**
- * Returns the center point of the rectangle.
- */
-glm::vec2 Rect::Center() const { return 0.5f * (max + min); }
-
-/**
- * Does this rectangle contain (includes on border) the given point?
- */
-bool Rect::Contains(const glm::vec2& p) const {
-  return glm::all(glm::greaterThanEqual(p, min)) &&
-         glm::all(glm::greaterThanEqual(max, p));
-}
-
-/**
- * Does this rectangle contain (includes equal) the given rectangle?
- */
-bool Rect::Contains(const Rect& rect) const {
-  return glm::all(glm::greaterThanEqual(rect.min, min)) &&
-         glm::all(glm::greaterThanEqual(max, rect.max));
-}
-
-/**
- * Does this rectangle overlap the one given (including equality)?
- */
-bool Rect::DoesOverlap(const Rect& rect) const {
-  return min.x <= rect.max.x && min.y <= rect.max.y && max.x >= rect.min.x &&
-         max.y >= rect.min.y;
-}
-
-/**
- * Is the rectangle empty (containing no space)?
- */
-bool Rect::IsEmpty() const { return max.y <= min.y || max.x <= min.x; };
-
-/**
- * Does this recangle have finite bounds?
- */
-bool Rect::IsFinite() const {
-  return glm::all(glm::isfinite(min)) && glm::all(glm::isfinite(max));
-}
-
-/**
- * Expand this rectangle (in place) to include the given point.
- */
-void Rect::Union(const glm::vec2 p) {
-  min = glm::min(min, p);
-  max = glm::max(max, p);
-}
-
-/**
- * Expand this rectangle to include the given Rect.
- */
-Rect Rect::Union(const Rect& rect) const {
-  Rect out;
-  out.min = glm::min(min, rect.min);
-  out.max = glm::max(max, rect.max);
-  return out;
-}
-
-/**
- * Shift this rectangle by the given vector.
- */
-Rect Rect::operator+(const glm::vec2 shift) const {
-  Rect out;
-  out.min = min + shift;
-  out.max = max + shift;
-  return out;
-}
-
-/**
- * Shift this rectangle in-place by the given vector.
- */
-Rect& Rect::operator+=(const glm::vec2 shift) {
-  min += shift;
-  max += shift;
-  return *this;
-}
-
-/**
- * Scale this rectangle by the given vector.
- */
-Rect Rect::operator*(const glm::vec2 scale) const {
-  Rect out;
-  out.min = min * scale;
-  out.max = max * scale;
-  return out;
-}
-
-/**
- * Scale this rectangle in-place by the given vector.
- */
-Rect& Rect::operator*=(const glm::vec2 scale) {
-  min *= scale;
-  max *= scale;
-  return *this;
-}
-
-/**
- * Transform the rectangle by the given axis-aligned affine transform.
- *
- * Ensure the transform passed in is axis-aligned (rotations are all
- * multiples of 90 degrees), or else the resulting rectangle will no longer
- * bound properly.
- */
-Rect Rect::Transform(const glm::mat3x2& m) const {
-  Rect rect;
-  rect.min = m * glm::vec3(min, 1);
-  rect.max = m * glm::vec3(max, 1);
-  return rect;
-}
-
-/**
- * Return a CrossSection with an outline defined by this axis-aligned
- * rectangle.
- */
-CrossSection Rect::AsCrossSection() const { return CrossSection(*this); }
-
 }  // namespace manifold

src/manifold/src/edge_op.cpp

@@ -145,24 +145,29 @@ namespace manifold {
 void Manifold::Impl::SimplifyTopology() {
   if (!halfedge_.size()) return;
 
-  int nbEdges = halfedge_.size();
+  const int nbEdges = halfedge_.size();
+  auto policy = autoPolicy(nbEdges);
   int numFlagged = 0;
   Vec<uint8_t> bflags(nbEdges);
 
-  Vec<SortEntry> entries(nbEdges);
-  auto policy = autoPolicy(halfedge_.size());
-  for_each_n(policy, countAt(0), nbEdges, [&](int i) {
-    entries[i].start = halfedge_[i].startVert;
-    entries[i].end = halfedge_[i].endVert;
-    entries[i].index = i;
-  });
+  // In the case of a very bad triangulation, it is possible to create pinched
+  // verts. They must be removed before edge collapse.
+  SplitPinchedVerts();
+
+  Vec<SortEntry> entries;
+  entries.reserve(nbEdges / 2);
+  for (int i = 0; i < nbEdges; ++i) {
+    if (halfedge_[i].IsForward()) {
+      entries.push_back({halfedge_[i].startVert, halfedge_[i].endVert, i});
+    }
+  }
 
   stable_sort(policy, entries.begin(), entries.end());
-  for (int i = 0; i < nbEdges - 1; ++i) {
+  for (int i = 0; i < entries.size() - 1; ++i) {
     if (entries[i].start == entries[i + 1].start &&
         entries[i].end == entries[i + 1].end) {
-      DedupeEdge(std::min(entries[i].index, entries[i + 1].index));
-      entries[i + 1].index = std::max(entries[i].index, entries[i + 1].index);
+      DedupeEdge(entries[i].index);
+      numFlagged++;
     }
   }
 
@@ -243,6 +248,8 @@ void Manifold::Impl::SimplifyTopology() {
 #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;
@@ -251,6 +258,7 @@ void Manifold::Impl::DedupeEdge(const int edge) {
   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]);
@@ -297,6 +305,45 @@ void Manifold::Impl::DedupeEdge(const int edge) {
 
     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]);
+
+    do {
+      halfedge_[current].endVert = newVert;
+      current = NextHalfedge(current);
+      halfedge_[current].startVert = newVert;
+      current = halfedge_[current].pairedHalfedge;
+    } while (current != edge);
+  }
+
+  // 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]);
+
+    do {
+      halfedge_[current].endVert = newVert;
+      current = NextHalfedge(current);
+      halfedge_[current].startVert = newVert;
+      current = halfedge_[current].pairedHalfedge;
+    } while (current != pair);
+  }
 }
 
 void Manifold::Impl::PairUp(int edge0, int edge1) {
@@ -380,6 +427,10 @@ void Manifold::Impl::RemoveIfFolded(int edge) {
   }
 }
 
+// Collapses the given edge by removing startVert. 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 will be marked NaN,
+// but other edges may still be pointing to it.
 void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
   Vec<TriRef>& triRef = meshRelation_.triRef;
   Vec<glm::ivec3>& triProp = meshRelation_.triProperties;