Add 2d hull operation (#486)

This adds a 2d convex hull operation to CrossSection (ported from OCADml, which I ported based on BOSL2).

src/cross_section/include/cross_section.h

@@ -153,6 +153,15 @@ class CrossSection {
   std::vector<CrossSection> Decompose() const;
   ///@}
 
+  /** @name Convex Hulling
+   */
+  ///@{
+  CrossSection Hull() const;
+  static CrossSection Hull(const std::vector<CrossSection>& crossSections);
+  static CrossSection Hull(const SimplePolygon poly);
+  static CrossSection Hull(const Polygons polys);
+  ///@}
+  ///
   /** @name Conversion
    */
   ///@{

src/cross_section/src/cross_section.cpp

@@ -132,6 +132,58 @@ void decompose_hole(const C2::PolyTreeD* outline,
   }
 }
 
+bool V2Lesser(glm::vec2 a, glm::vec2 b) {
+  if (a.x == b.x) return a.y < b.y;
+  return a.x < b.x;
+}
+
+double Cross(glm::vec2 a, glm::vec2 b) { return a.x * b.y - a.y * b.x; }
+
+double SquaredDistance(glm::vec2 a, glm::vec2 b) {
+  auto d = a - b;
+  return glm::dot(d, d);
+}
+
+bool IsCw(glm::vec2 a, glm::vec2 b, glm::vec2 c) {
+  double lhs = Cross(a - c, b - c);
+  double rhs = 1e-18 * SquaredDistance(a, c) * SquaredDistance(b, c);
+  return lhs * std::abs(lhs) <= rhs;
+}
+
+void HullBacktrack(const SimplePolygon& pts, const int idx,
+                   std::vector<int>& keep, const int hold) {
+  const int stop = keep.size() - hold;
+  int i = 0;
+  while (i < stop && !IsCw(pts[idx], pts[keep[keep.size() - 1]],
+                           pts[keep[keep.size() - 2]])) {
+    keep.pop_back();
+    i++;
+  }
+}
+
+// Based on method described here:
+// https://www.hackerearth.com/practice/math/geometry/line-sweep-technique/tutorial/
+C2::PathD HullImpl(SimplePolygon& pts) {
+  int len = pts.size();
+  if (len < 3) return C2::PathD();  // not enough points to create a polygon
+  std::sort(pts.begin(), pts.end(), V2Lesser);
+  auto keep = std::vector<int>{0, 1};
+  for (int i = 2; i < len; i++) {
+    HullBacktrack(pts, i, keep, 1);
+    keep.push_back(i);
+  }
+  int nLower = keep.size();
+  for (int i = 0; i < len - 1; i++) {
+    int idx = len - 2 - i;
+    HullBacktrack(pts, idx, keep, nLower);
+    if (idx > 0) keep.push_back(idx);
+  }
+  auto path = C2::PathD(keep.size());
+  for (int i = 0; i < keep.size(); i++) {
+    path[i] = v2_to_pd(pts[keep[i]]);
+  }
+  return path;
+}
 }  // namespace
 
 namespace manifold {
@@ -560,6 +612,41 @@ CrossSection CrossSection::Offset(double delta, JoinType jointype,
   return CrossSection(ps);
 }
 
+CrossSection CrossSection::Hull(
+    const std::vector<CrossSection>& crossSections) {
+  int n = 0;
+  for (auto cs : crossSections) n += cs.NumVert();
+  SimplePolygon pts;
+  pts.reserve(n);
+  for (auto cs : crossSections) {
+    auto paths = cs.GetPaths();
+    for (auto path : paths) {
+      for (auto p : path) {
+        pts.push_back(v2_of_pd(p));
+      }
+    }
+  }
+  return CrossSection(C2::PathsD{HullImpl(pts)});
+}
+
+CrossSection CrossSection::Hull() const {
+  return Hull(std::vector<CrossSection>{*this});
+}
+
+CrossSection CrossSection::Hull(SimplePolygon pts) {
+  return CrossSection(C2::PathsD{HullImpl(pts)});
+}
+
+CrossSection CrossSection::Hull(const Polygons polys) {
+  SimplePolygon pts;
+  for (auto poly : polys) {
+    for (auto p : poly) {
+      pts.push_back(p);
+    }
+  }
+  return Hull(pts);
+}
+
 /**
  * Return the total area covered by complex polygons making up the
  * CrossSection.

test/cross_section_test.cpp

@@ -16,6 +16,8 @@
 
 #include <gtest/gtest.h>
 
+#include <vector>
+
 #include "glm/geometric.hpp"
 #include "manifold.h"
 #include "polygon.h"
@@ -167,3 +169,25 @@ TEST(CrossSection, FillRule) {
   CrossSection nonZero(polygon, CrossSection::FillRule::NonZero);
   EXPECT_NEAR(nonZero.Area(), 0.875, 0.001);
 }
+
+TEST(CrossSection, Hull) {
+  auto circ = CrossSection::Circle(10, 360);
+  auto circs = std::vector<CrossSection>{circ, circ.Translate({0, 30}),
+                                         circ.Translate({30, 0})};
+  auto circ_tri = CrossSection::Hull(circs);
+  auto centres = SimplePolygon{{0, 0}, {0, 30}, {30, 0}, {15, 5}};
+  auto tri = CrossSection::Hull(centres);
+
+#ifdef MANIFOLD_EXPORT
+  if (options.exportModels) {
+    auto circ_tri_ex = Manifold::Extrude(circ_tri, 10);
+    ExportMesh("cross_section_hull_circ_tri.glb", circ_tri_ex.GetMesh(), {});
+  }
+#endif
+
+  auto circ_area = circ.Area();
+  EXPECT_FLOAT_EQ(circ_area, (circ - circ.Scale({0.8, 0.8})).Hull().Area());
+  EXPECT_FLOAT_EQ(
+      circ_area * 2.5,
+      (CrossSection::BatchBoolean(circs, OpType::Add) - tri).Area());
+}