From 1fcdde0b50fc384a8b16643bdaae3140c2c28219 Mon Sep 17 00:00:00 2001 From: Guy Blelloch Date: Tue, 7 May 2024 18:03:40 -0400 Subject: [PATCH] added 3d convex hull --- examples/CMakeLists.txt | 2 +- examples/README.md | 3 +- examples/convex_hull_3d.cpp | 39 ++++++++ examples/convex_hull_3d.h | 177 ++++++++++++++++++++++++++++++++++++ examples/delaunay.h | 5 +- 5 files changed, 223 insertions(+), 3 deletions(-) create mode 100644 examples/convex_hull_3d.cpp create mode 100644 examples/convex_hull_3d.h diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a980a096..f224044d 100644 --- a/examples/CMakeLists.txt +++ b/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) diff --git a/examples/README.md b/examples/README.md index 86b897dd..e55181f8 100644 --- a/examples/README.md +++ b/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 diff --git a/examples/convex_hull_3d.cpp b/examples/convex_hull_3d.cpp new file mode 100644 index 00000000..6c3c60a5 --- /dev/null +++ b/examples/convex_hull_3d.cpp @@ -0,0 +1,39 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include "convex_hull_3d.h" + +// ************************************************************** +// Driver +// ************************************************************** +int main(int argc, char* argv[]) { + auto usage = "Usage: "; + 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 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 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; + } +} diff --git a/examples/convex_hull_3d.h b/examples/convex_hull_3d.h new file mode 100644 index 00000000..7c69d28d --- /dev/null +++ b/examples/convex_hull_3d.h @@ -0,0 +1,177 @@ +#include +#include +#include + +#include +#include +#include + +#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; +using edge = std::array; + +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 conflicts; + triangle(); + triangle(tri t, point_id pid, parlay::sequence 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; + +struct Convex_Hull_3d { + using triangle_ptr = std::shared_ptr; + hash_map map_facets; + Points points; + hash_map 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(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(6*P.size())), + convex_hull(hash_map(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(init_tri[i], remain[i], + parlay::filter(target_points, test)); + }); + + // The 6 ridges of the initial tetrahedron along with their neighboring triangles + std::tuple 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 convex_hull_3d(const Points& P) { + Convex_Hull_3d ch3d(P); + return ch3d.convex_hull.keys(); +} diff --git a/examples/delaunay.h b/examples/delaunay.h index d37d9f81..cd051905 100644 --- a/examples/delaunay.h +++ b/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,7 +122,6 @@ struct Delaunay { mesh(hash_map(2*P.size())), edges(hash_map(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}; @@ -127,6 +129,7 @@ struct Delaunay { points = parlay::append(P, Points({p0, p1, p2})); auto t = std::make_shared(tri{n, n+1, n+2}, P); n += 3; + // a dummy triangle with no points, for the boundary auto te = std::make_shared(tri{-1, -1, -1}, Points()); std::shared_ptr te2, te3, t2, t3; t2 = t3 = t;