added 3d convex hull

examples/CMakeLists.txt

@@ -12,7 +12,7 @@ set(EXAMPLES "primes" "BFS" "word_counts" "tokens" "filter" "linefit"
   "knuth_morris_pratt" "huffman_tree" "decision_tree_c45" "karatsuba" "suffix_tree"
   "2d_linear_program" "box_kdtree" "radix_tree" "ray_trace" "hash_map" "oct_tree"
   "3d_range" "rectangle_intersection" "star_connectivity" "ldd_connectivity" "boruvka"
-  "counting_sort" "integer_sort" "lubys")
+  "counting_sort" "integer_sort" "lubys" "convex_hull_3d")
 
 function(add_example NAME)
   add_executable(${NAME} ${NAME}.cpp)

examples/README.md

@@ -77,7 +77,8 @@ The examples are as follows, broken down into categories.
 ## Geometry and Graphics
 
 - 2d_linear_program : solve linear program max c^Tx, Ax <= b, in 2 dimensions.
-- delaunay : Delaunay triangulation in 2 dimensions
+- convex_hull_3d : convex hull of points in 3 dimensions using parallel incremental algorithm
+- delaunay : Delaunay triangulation in 2 dimensions using parallel incremental algorithm
 - oct_tree : the z-tree variant of oct trees in any constant dimension
 - knn : k-nearest-neighbors of a set of points in constant dimension
 - quickhull : convex hull of points in 2 dimensions using quickhull algorithm

examples/convex_hull_3d.cpp

@@ -0,0 +1,39 @@
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <random>
+
+#include <parlay/primitives.h>
+#include <parlay/sequence.h>
+#include <parlay/random.h>
+
+#include "convex_hull_3d.h"
+
+// **************************************************************
+// Driver
+// **************************************************************
+int main(int argc, char* argv[]) {
+  auto usage = "Usage: <n>";
+  if (argc != 2) std::cout << usage << std::endl;
+  else {
+    point_id n;
+    try {n = std::stoi(argv[1]);}
+    catch (...) {std::cout << usage << std::endl; return 1;}
+    parlay::random_generator gen(0);
+    std::uniform_real_distribution<real> dis(0.0,1.0);
+
+    // generate n random points in a unit square
+    auto points = parlay::tabulate(n, [&] (point_id i) -> point {
+      auto r = gen[i];
+      return point{i, dis(r), dis(r), dis(r)};});
+
+    parlay::internal::timer t("Time");
+    parlay::sequence<tri> result;
+    for (int i=0; i < 5; i++) {
+      result = convex_hull_3d(points);
+      t.next("convex_hull_3d");
+    }
+
+    std::cout << "number of triangles in the mesh = " << result.size() << std::endl;
+  }
+}

examples/convex_hull_3d.h

@@ -0,0 +1,177 @@
+#include <utility>
+#include <memory>
+#include <optional>
+
+#include <parlay/primitives.h>
+#include <parlay/sequence.h>
+#include <parlay/utilities.h>
+
+#include "hash_map.h"
+
+// **************************************************************
+// Parallel Convel Hull in 3D
+// From the paper:
+// Randomized Incremental Convex Hull is Highly Parallel
+// Blelloch, Gu, Shun and Sun
+//
+// Author: Jekyeom Jeon
+// **************************************************************
+
+using real = float;
+using point_id = int;
+
+using tri = std::array<point_id,3>;
+using edge = std::array<point_id,2>;
+
+struct point {
+  point_id id; real x; real y; real z;
+  const bool operator<(const point p) const {return id < p.id;}
+};
+
+struct vect {
+  double x; double y; double z;
+  vect(double x, double y, double z) : x(x), y(y), z(z) {}
+};
+
+struct triangle {
+  tri t;
+  point_id pid; // point in convex hull not in triangle
+  parlay::sequence<point> conflicts;
+  triangle();
+  triangle(tri t, point_id pid, parlay::sequence<point> c) : t(t), pid(pid), conflicts(c) {}
+};
+
+// A "staged" version of a test if two points (p and q) are on the
+// opposite side of a plane defined by three other points (a, b, c).
+// It first takes the 4 points (a, b, c, p) and returns a function
+// which takes the point q, answering the query.  Since we often check
+// multiple points against a single circle this reduces the work.
+// Works by taking the cross product of c->a and c->b and checking that
+// the dot of this with c->p and c->q have opposite signs.
+inline auto is_opposite (point a, point b, point c, point p) {
+  auto v_from_c = [=] (point p) { return vect(p.x - c.x, p.y - c.y, p.z - c.z);};
+  auto dot = [] (vect v1, vect v2) { return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;};
+  auto cross = [] (vect v1, vect v2) {
+    return vect{v1.y*v2.z - v1.z*v2.y,
+                v1.z*v2.x - v1.x*v2.z,
+                v1.x*v2.y - v1.y*v2.x};};
+  vect cp = cross(v_from_c(a), v_from_c(b));
+  bool p_side = dot(cp, v_from_c(p)) > 0.0;
+  return [=] (point q) -> bool {
+    return ((dot(cp, v_from_c(q)) > 0.0) != p_side);};
+}
+
+// **************************************************************
+// The main body
+// **************************************************************
+
+using Points = parlay::sequence<point>;
+
+struct Convex_Hull_3d {
+  using triangle_ptr = std::shared_ptr<triangle>;
+  hash_map<edge, triangle_ptr> map_facets;
+  Points points;
+  hash_map<tri, bool> convex_hull;
+  point_id n;
+
+  point_id min_conflicts(triangle_ptr& t) {
+    return (t->conflicts.size() == 0) ? n : t->conflicts[0].id;}
+
+  // recursive routine
+  // convex hull in 3d space
+  void process_ridge(triangle_ptr& t1, edge r, triangle_ptr& t2) {
+    if (t1->conflicts.size() == 0 && t2->conflicts.size() == 0) {
+      return;
+    } else if (min_conflicts(t2) == min_conflicts(t1)) {
+      // H \ {t1, t2}
+      convex_hull.remove(t1->t);
+      convex_hull.remove(t2->t);
+    } else if (min_conflicts(t2) < min_conflicts(t1)) {
+      process_ridge(t2, r, t1);
+    } else {
+      point_id pid = min_conflicts(t1);
+      tri t{r[0], r[1], pid};
+
+      // C(t) <- v' in C(t1) U C(t2) | visible(v', t)
+      auto t1_U_t2 = parlay::merge(t1->conflicts, t2->conflicts);
+
+      // removes first point, removes duplicates, and checks for visibility
+      auto test = is_opposite(points[t[0]], points[t[1]], points[t[2]], points[t1->pid]);
+      auto keep = parlay::tabulate(t1_U_t2.size(), [&] (long i) {
+          return (((i != 0) && (t1_U_t2[i].id != t1_U_t2[i-1].id)) &&
+                  test(t1_U_t2[i]));});
+      auto p = parlay::pack(t1_U_t2, keep);
+
+      auto t_new = std::make_shared<triangle>(t, t1->pid, p);
+
+      // H <- (H \ {t1}) U {t}
+      convex_hull.remove(t1->t);
+      convex_hull.insert(t, true);
+
+      auto check_edge = [&] (edge e, triangle_ptr& tp) {
+        auto key = (e[0] < e[1]) ? e : edge{e[1], e[0]};
+        if (map_facets.insert(key, tp)) return;
+        auto tt = map_facets.remove(key).value();
+        process_ridge(tp, e, tt);};
+
+      parlay::par_do3([&] {process_ridge(t_new, r, t2);},
+                      [&] {check_edge(edge{r[0], pid}, t_new);},
+                      [&] {check_edge(edge{r[1], pid}, t_new);});
+    }
+  }
+
+  // toplevel routine
+  // The result is set of facets: result_facets
+  // assum that P has more than 4 points
+  Convex_Hull_3d(const Points& P) :
+      map_facets(hash_map<edge, triangle_ptr>(6*P.size())),
+      convex_hull(hash_map<tri, bool>(6*P.size())),
+      n(P.size()) {
+
+    points = P;
+    // first 4 points are used to define an initial tetrahedron with 4 faces
+    tri init_tri[4] = {{0, 1, 2},
+                       {1, 2, 3},
+                       {0, 2, 3},
+                       {0, 1, 3}};
+    // points on the inside of each face (i.e. the missing point)
+    point_id remain[4] = {3, 0, 1, 2}; 
+
+    // insert the initial convex hull
+    for (auto t : init_tri)
+      convex_hull.insert(t, true);
+
+    // remove first 4 points
+    Points target_points = P.subseq(4, P.size());
+
+    // initialize each of the 4 faces of the tetrahedron
+    triangle_ptr t[4];
+    parlay::parallel_for(0, 4, [&] (int i) {
+        point p0 = points[init_tri[i][0]];
+        point p1 = points[init_tri[i][1]];
+        point p2 = points[init_tri[i][2]];
+        point pc = points[remain[i]];
+        auto test = is_opposite(p0, p1, p2, pc);
+        t[i] = std::make_shared<triangle>(init_tri[i], remain[i],
+                                          parlay::filter(target_points, test));
+      });
+
+    // The 6 ridges of the initial tetrahedron along with their neighboring triangles
+    std::tuple<int, int, edge> ridges[6] =
+      {{0, 1, edge{1, 2}},
+       {0, 2, edge{0, 2}},
+       {0, 3, edge{0, 1}},
+       {1, 2, edge{2, 3}},
+       {1, 3, edge{1, 3}},
+       {2, 3, edge{0, 3}}};
+
+    parlay::parallel_for(0, 6, [&](auto i) {
+        auto [t1_idx, t2_idx, e] = ridges[i];
+        process_ridge(t[t1_idx], e, t[t2_idx]); });
+  }
+};
+
+parlay::sequence<tri> convex_hull_3d(const Points& P) {
+  Convex_Hull_3d ch3d(P);
+  return ch3d.convex_hull.keys();
+}

examples/delaunay.h

@@ -77,7 +77,9 @@ struct Delaunay {
     auto a = parlay::merge(t1->conflicts, t2->conflicts);
     auto is_in_circle = in_circle(points[t[0]],points[t[1]],points[t[2]]);
     auto keep = parlay::tabulate(a.size(), [&] (long i) {
+        // throw out the first point and one of the copies if appears in both
       return ((i != 0) && (a[i].id != a[i-1].id) &&
+              // allways keep a point if it appears in both
               ((i+1 < a.size() && a[i].id == a[i+1].id) ||
                is_in_circle(a[i])));},500);
     return parlay::pack(a, keep);
@@ -88,6 +90,7 @@ struct Delaunay {
   // 2d Delaunay instead of arbitrary convex hull
   void process_edge(triangle_ptr& t1, edge e, triangle_ptr& t2) {
     if (t1->conflicts.size() == 0 && t2->conflicts.size() == 0) {
+      // only add triangles to mesh when they are definitely in there
       mesh.insert(t1->t,true); mesh.insert(t2->t,true);
       t1 = t2 = nullptr;
     } else if (earliest(t2) == earliest(t1)) {
@@ -119,14 +122,14 @@ struct Delaunay {
       mesh(hash_map<tri,bool>(2*P.size())),
       edges(hash_map<edge,triangle_ptr>(6*P.size())),
       n(P.size()) {
-    points = P;
     // enclosing triangle
     point p0{n,0.0,100.0};
     point p1{n+1,100.0,-100.0};
     point p2{n+2,-100.0,-100.0};
     points = parlay::append(P, Points({p0, p1, p2}));
     auto t = std::make_shared<triangle>(tri{n, n+1, n+2}, P);
     n += 3;
+    // a dummy triangle with no points, for the boundary
     auto te = std::make_shared<triangle>(tri{-1, -1, -1}, Points());
     std::shared_ptr<triangle> te2, te3, t2, t3;
     t2 = t3 = t;