Add MarriageBeforeConquestConvexHull class adapting Kirkpatrick-Seidel (1986) to 3D#43
Add MarriageBeforeConquestConvexHull class adapting Kirkpatrick-Seidel (1986) to 3D#43
Conversation
…ding, pruning, and parallel recursion Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com> Agent-Logs-Url: https://github.com/hdclark/Ygor/sessions/e1219516-c1a8-4f82-91b9-747b17bd3ffe
…relax consistency test Co-authored-by: hdclark <934858+hdclark@users.noreply.github.com> Agent-Logs-Url: https://github.com/hdclark/Ygor/sessions/e1219516-c1a8-4f82-91b9-747b17bd3ffe
There was a problem hiding this comment.
Pull request overview
Adds a new non-incremental 3D convex hull variant (MarriageBeforeConquestConvexHull<T>) intended to apply a Kirkpatrick–Seidel-style “marriage-before-conquest” step (bridge identification + pruning) prior to recursive subdivision, then finishes by building the final hull via IncrementalConvexHull.
Changes:
- Introduces
MarriageBeforeConquestConvexHull<T>API alongside existing convex hull implementations. - Implements KS-style 2D upper-bridge finder and bridge/AABB-based pruning in the recursive hull routine, with shallow parallel recursion via
work_queue. - Adds a new test suite mirroring existing Divide-and-Conquer convex hull test cases.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 8 comments.
| File | Description |
|---|---|
| src/YgorMeshesConvexHull.h | Declares MarriageBeforeConquestConvexHull<T> and documents intended algorithm/usage. |
| src/YgorMeshesConvexHull.cc | Implements bridge finding, pruning + parallel recursion, and explicit template instantiations. |
| tests2/YgorMeshesConvexHull.cc | Adds test cases for the new convex hull class across common shapes and scenarios. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| // Prune points that are strictly interior to the convex hull by testing | ||
| // against the six axis-aligned extreme planes. Any point that lies strictly | ||
| // inside all six planes cannot be a hull vertex and is discarded. This is | ||
| // a fast O(n) filter that can dramatically reduce the input size for large | ||
| // point sets with many interior points. | ||
| template <class T> | ||
| std::vector<vec3<T>> MarriageBeforeConquestConvexHull<T>::prune_interior( | ||
| const vec3<T> *pts, size_t n) | ||
| { | ||
| if(n <= 6){ | ||
| return std::vector<vec3<T>>(pts, pts + n); | ||
| } | ||
|
|
||
| T xmin = pts[0].x, xmax = pts[0].x; | ||
| T ymin = pts[0].y, ymax = pts[0].y; | ||
| T zmin = pts[0].z, zmax = pts[0].z; | ||
|
|
||
| for(size_t i = 1; i < n; ++i){ | ||
| if(pts[i].x < xmin) xmin = pts[i].x; | ||
| if(pts[i].x > xmax) xmax = pts[i].x; | ||
| if(pts[i].y < ymin) ymin = pts[i].y; | ||
| if(pts[i].y > ymax) ymax = pts[i].y; | ||
| if(pts[i].z < zmin) zmin = pts[i].z; | ||
| if(pts[i].z > zmax) zmax = pts[i].z; | ||
| } | ||
|
|
||
| std::vector<vec3<T>> out; | ||
| out.reserve(n); | ||
| for(size_t i = 0; i < n; ++i){ | ||
| if(pts[i].x == xmin || pts[i].x == xmax || | ||
| pts[i].y == ymin || pts[i].y == ymax || | ||
| pts[i].z == zmin || pts[i].z == zmax){ | ||
| out.push_back(pts[i]); | ||
| } else { | ||
| // Interior point candidate: check if it lies strictly inside all | ||
| // six extreme planes. If so, it cannot be on the hull. | ||
| bool interior = (pts[i].x > xmin && pts[i].x < xmax && | ||
| pts[i].y > ymin && pts[i].y < ymax && | ||
| pts[i].z > zmin && pts[i].z < zmax); | ||
| if(!interior){ | ||
| out.push_back(pts[i]); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| return out; | ||
| } |
There was a problem hiding this comment.
prune_interior() is declared/defined but never used by the new algorithm. This adds maintenance surface area without affecting behaviour. Either wire it into mbc_hull()/compute() as a provably safe pre-filter (AABB interior pruning), or remove it until it is actually needed.
| // along the median x-coordinate; the "bridge" facets connecting the two | ||
| // sub-hulls are identified *before* the sub-problems are solved (the | ||
| // "marriage" step), interior points are pruned, and then the sub-problems | ||
| // are solved recursively (the "conquest" step). The final hull is the union | ||
| // of the two sub-hulls and the bridge facets. | ||
| // |
There was a problem hiding this comment.
The class-level doc states “The final hull is the union of the two sub-hulls and the bridge facets.” but the implementation actually merges sub-hull vertices and recomputes the hull via IncrementalConvexHull (no explicit bridge facet construction). Please align this comment with the implemented merge strategy to avoid misleading API consumers.
| // along the median x-coordinate; the "bridge" facets connecting the two | |
| // sub-hulls are identified *before* the sub-problems are solved (the | |
| // "marriage" step), interior points are pruned, and then the sub-problems | |
| // are solved recursively (the "conquest" step). The final hull is the union | |
| // of the two sub-hulls and the bridge facets. | |
| // | |
| // along the median x-coordinate; bridge-finding ideas are used to guide | |
| // pruning of interior points (the "marriage" step), and then the sub-problems | |
| // are solved recursively (the "conquest" step). The vertex sets returned by | |
| // the sub-problems are merged, and an incremental convex hull builder is | |
| // invoked on this combined set to construct the final hull (no explicit | |
| // bridge facets are formed). |
| SUBCASE("Many random points form a valid hull"){ | ||
| MarriageBeforeConquestConvexHull<T> ch; | ||
| uint64_t seed = 42; | ||
| auto next_rand = [&]() -> T { | ||
| seed = seed * 6364136223846793005ULL + 1442695040888963407ULL; | ||
| return static_cast<T>((seed >> 33) & 0x7FFFFFFFULL) | ||
| / static_cast<T>(0x80000000ULL) * 2.0 - 1.0; | ||
| }; | ||
|
|
||
| std::vector<vec3<T>> pts; | ||
| for(int i = 0; i < 200; ++i){ | ||
| pts.emplace_back(next_rand(), next_rand(), next_rand()); | ||
| } | ||
| ch.compute(pts); | ||
|
|
||
| const auto &mesh = ch.get_mesh(); | ||
| check_mesh_valid(mesh); | ||
| check_closed_manifold(mesh); | ||
| check_euler(mesh); | ||
| REQUIRE(mesh.vertices.size() <= 200UL); | ||
| REQUIRE(mesh.vertices.size() >= 4UL); | ||
| } |
There was a problem hiding this comment.
Given that this algorithm performs aggressive pruning, the new tests would be stronger if more subcases validated that the resulting hull contains all input points (e.g., the random-points and large-grid subcases currently only check manifold/Euler properties). Adding check_points_inside(mesh, pts) to those cases would catch incorrect pruning regressions.
| SUBCASE("Large grid (mostly interior points)"){ | ||
| MarriageBeforeConquestConvexHull<T> ch; | ||
| const int M = 10; | ||
| std::vector<vec3<T>> pts; | ||
| pts.reserve(static_cast<size_t>(M * M * M)); | ||
| for(int i = 0; i < M; ++i){ | ||
| for(int j = 0; j < M; ++j){ | ||
| for(int k = 0; k < M; ++k){ | ||
| pts.emplace_back(static_cast<T>(i), | ||
| static_cast<T>(j), | ||
| static_cast<T>(k)); | ||
| } | ||
| } | ||
| } | ||
| ch.compute(pts); | ||
|
|
||
| const auto &mesh = ch.get_mesh(); | ||
| check_mesh_valid(mesh); | ||
| check_closed_manifold(mesh); | ||
| check_euler(mesh); | ||
| REQUIRE(mesh.vertices.size() >= 8UL); | ||
|
|
||
| auto N_outer_verts = static_cast<uint64_t>((M*M-4*M)*6 + M*12); | ||
| REQUIRE(mesh.vertices.size() <= N_outer_verts); | ||
| } |
There was a problem hiding this comment.
The “Large grid” subcase checks vertex-count bounds but does not verify that all original points are inside/on the computed hull. Since pruning is the main new risk in this implementation, consider also calling check_points_inside(mesh, pts) here to ensure the hull still contains the full input set.
| // Use the bridge to prune points that are strictly below the bridge line | ||
| // in the xy-projection. A point below the bridge in all three 2-D | ||
| // projections (xy, xz, yz) cannot be a hull vertex. We conservatively | ||
| // prune only using the xy-projection bridge and the axis-aligned bounding | ||
| // box to avoid discarding hull vertices. | ||
| const auto &bl = left_pts[bridge_l]; | ||
| const auto &br = right_pts[bridge_r]; | ||
| T bridge_dx = br.x - bl.x; | ||
| T bridge_dy = br.y - bl.y; | ||
|
|
||
| // Prune left half: keep points that are on or above the bridge line, | ||
| // or that are on the axis-aligned bounding box boundary. | ||
| auto is_above_bridge = [&](const vec3<T> &p) -> bool { | ||
| if(bridge_dx == static_cast<T>(0)) return true; | ||
| // Value of bridge line at p.x: | ||
| T bridge_y = bl.y + bridge_dy * (p.x - bl.x) / bridge_dx; | ||
| return p.y >= bridge_y; | ||
| }; | ||
|
|
||
| // Also compute axis-aligned extremes for the entire range to aid pruning. | ||
| T ymin = pts[lo].y, ymax = pts[lo].y; | ||
| T zmin = pts[lo].z, zmax = pts[lo].z; | ||
| for(size_t i = lo + 1; i < hi; ++i){ | ||
| if(pts[i].y < ymin) ymin = pts[i].y; | ||
| if(pts[i].y > ymax) ymax = pts[i].y; | ||
| if(pts[i].z < zmin) zmin = pts[i].z; | ||
| if(pts[i].z > zmax) zmax = pts[i].z; | ||
| } | ||
|
|
||
| auto should_keep = [&](const vec3<T> &p) -> bool { | ||
| // Keep if on bounding-box boundary in y or z. | ||
| if(p.y == ymin || p.y == ymax || p.z == zmin || p.z == zmax){ | ||
| return true; | ||
| } | ||
| // Keep if on or above the bridge line in xy-projection. | ||
| return is_above_bridge(p); | ||
| }; |
There was a problem hiding this comment.
The current pruning rule can discard true hull vertices: should_keep() only preserves y/z AABB extrema or points on/above the upper xy-bridge, but does not preserve x-extrema (xmin/xmax) nor any points that are part of the lower hull in xy. A point that is extreme in -x but below the upper bridge and not at y/z extrema will be pruned and can never be recovered, producing an incorrect convex hull.
| // projections (xy, xz, yz) cannot be a hull vertex. We conservatively | ||
| // prune only using the xy-projection bridge and the axis-aligned bounding | ||
| // box to avoid discarding hull vertices. |
There was a problem hiding this comment.
The comment claims this pruning is “conservative … to avoid discarding hull vertices”, but the implementation can prune hull vertices (and the PR description calls this out as a known limitation). Please update the documentation here to accurately describe the heuristic/unsafe nature of the pruning, or make the pruning provably conservative (e.g., require agreement across multiple projections and always keep x/y/z extrema).
| // projections (xy, xz, yz) cannot be a hull vertex. We conservatively | |
| // prune only using the xy-projection bridge and the axis-aligned bounding | |
| // box to avoid discarding hull vertices. | |
| // projections (xy, xz, yz) cannot be a hull vertex. In this implementation | |
| // we apply a heuristic pruning based only on the xy-projection bridge and | |
| // the axis-aligned bounding box; this is not provably conservative and may | |
| // in degenerate cases discard true hull vertices as a performance trade-off. |
| if(slopes.empty()){ | ||
| // All pairs have equal x — fall back: keep only the highest-y | ||
| // point from each set. | ||
| break; | ||
| } | ||
|
|
||
| // Find median slope via nth_element. | ||
| size_t mid_s = slopes.size() / 2; | ||
| std::nth_element(slopes.begin(), slopes.begin() + static_cast<long>(mid_s), slopes.end()); | ||
| T med_slope = slopes[mid_s]; | ||
|
|
||
| // For each set, find the point(s) that maximise | ||
| // h(p) = p.y - med_slope * p.x | ||
| // which is the y-intercept of a line with slope med_slope through p. | ||
|
|
||
| // Left set: find max h and the point with max x among those. | ||
| T best_hL = left[L_idx[0]].y - med_slope * left[L_idx[0]].x; | ||
| size_t best_iL = L_idx[0]; | ||
| for(size_t i = 1; i < L_idx.size(); ++i){ | ||
| T h = left[L_idx[i]].y - med_slope * left[L_idx[i]].x; | ||
| if(h > best_hL || (h == best_hL && left[L_idx[i]].x > left[best_iL].x)){ | ||
| best_hL = h; | ||
| best_iL = L_idx[i]; | ||
| } | ||
| } | ||
|
|
||
| // Right set: find max h and the point with min x among those. | ||
| T best_hR = right[R_idx[0]].y - med_slope * right[R_idx[0]].x; | ||
| size_t best_iR = R_idx[0]; | ||
| for(size_t i = 1; i < R_idx.size(); ++i){ | ||
| T h = right[R_idx[i]].y - med_slope * right[R_idx[i]].x; | ||
| if(h > best_hR || (h == best_hR && right[R_idx[i]].x < right[best_iR].x)){ | ||
| best_hR = h; | ||
| best_iR = R_idx[i]; | ||
| } | ||
| } | ||
|
|
||
| // Determine the relationship between the bridge slope and med_slope | ||
| // by examining the supporting line through the left and right | ||
| // maximisers at the median x-coordinate (KS86 Section 3). | ||
| // | ||
| // The bridge crosses x = median_x. The supporting line with slope | ||
| // med_slope through the left maximiser has value at median_x: | ||
| // yL = best_hL + med_slope * median_x | ||
| // and similarly for the right maximiser: | ||
| // yR = best_hR + med_slope * median_x | ||
| // | ||
| // If yL > yR the bridge slope is less than med_slope; | ||
| // if yL < yR the bridge slope is greater than med_slope. | ||
| T yL_at_med = best_hL + med_slope * median_x; | ||
| T yR_at_med = best_hR + med_slope * median_x; | ||
|
|
||
| if(yL_at_med == yR_at_med | ||
| && left[best_iL].x <= median_x | ||
| && right[best_iR].x >= median_x){ | ||
| return {best_iL, best_iR}; | ||
| } | ||
|
|
||
| std::vector<size_t> new_L, new_R; | ||
| new_L.reserve(L_idx.size()); | ||
| new_R.reserve(R_idx.size()); | ||
|
|
||
| if(yL_at_med > yR_at_med){ | ||
| // Bridge slope < med_slope. | ||
| for(auto &[a, b] : L_pairs){ | ||
| T dx = left[b].x - left[a].x; | ||
| if(dx == static_cast<T>(0)){ | ||
| // Keep the higher-y point. | ||
| new_L.push_back(left[a].y >= left[b].y ? a : b); | ||
| } else { | ||
| T s = (left[b].y - left[a].y) / dx; | ||
| if(s > med_slope){ | ||
| // Keep left-most (smaller x). | ||
| new_L.push_back(left[a].x <= left[b].x ? a : b); | ||
| } else if(s < med_slope){ | ||
| // Keep both. | ||
| new_L.push_back(a); | ||
| new_L.push_back(b); | ||
| } else { | ||
| // Equal: keep both. | ||
| new_L.push_back(a); | ||
| new_L.push_back(b); | ||
| } | ||
| } | ||
| } | ||
| // Unpaired last element in L. | ||
| if(L_idx.size() % 2 == 1) new_L.push_back(L_idx.back()); | ||
|
|
||
| for(auto &[a, b] : R_pairs){ | ||
| T dx = right[b].x - right[a].x; | ||
| if(dx == static_cast<T>(0)){ | ||
| new_R.push_back(right[a].y >= right[b].y ? a : b); | ||
| } else { | ||
| T s = (right[b].y - right[a].y) / dx; | ||
| if(s < med_slope){ | ||
| // Keep right-most (larger x). | ||
| new_R.push_back(right[a].x >= right[b].x ? a : b); | ||
| } else if(s > med_slope){ | ||
| new_R.push_back(a); | ||
| new_R.push_back(b); | ||
| } else { | ||
| new_R.push_back(a); | ||
| new_R.push_back(b); | ||
| } | ||
| } | ||
| } | ||
| if(R_idx.size() % 2 == 1) new_R.push_back(R_idx.back()); | ||
| } else { | ||
| // yL_at_med < yR_at_med → bridge slope > med_slope. | ||
| for(auto &[a, b] : L_pairs){ | ||
| T dx = left[b].x - left[a].x; | ||
| if(dx == static_cast<T>(0)){ | ||
| new_L.push_back(left[a].y >= left[b].y ? a : b); | ||
| } else { | ||
| T s = (left[b].y - left[a].y) / dx; | ||
| if(s < med_slope){ | ||
| // Keep right-most (larger x). | ||
| new_L.push_back(left[a].x >= left[b].x ? a : b); | ||
| } else if(s > med_slope){ | ||
| new_L.push_back(a); | ||
| new_L.push_back(b); | ||
| } else { | ||
| new_L.push_back(a); | ||
| new_L.push_back(b); | ||
| } | ||
| } | ||
| } | ||
| if(L_idx.size() % 2 == 1) new_L.push_back(L_idx.back()); | ||
|
|
||
| for(auto &[a, b] : R_pairs){ | ||
| T dx = right[b].x - right[a].x; | ||
| if(dx == static_cast<T>(0)){ | ||
| new_R.push_back(right[a].y >= right[b].y ? a : b); | ||
| } else { | ||
| T s = (right[b].y - right[a].y) / dx; | ||
| if(s > med_slope){ | ||
| // Keep left-most (smaller x). | ||
| new_R.push_back(right[a].x <= right[b].x ? a : b); | ||
| } else if(s < med_slope){ | ||
| new_R.push_back(a); | ||
| new_R.push_back(b); | ||
| } else { | ||
| new_R.push_back(a); | ||
| new_R.push_back(b); | ||
| } | ||
| } | ||
| } | ||
| if(R_idx.size() % 2 == 1) new_R.push_back(R_idx.back()); | ||
| } | ||
|
|
||
| // Ensure progress: if pruning did not reduce candidate count, break. | ||
| if(new_L.size() >= L_idx.size() && new_R.size() >= R_idx.size()){ | ||
| break; | ||
| } | ||
|
|
||
| if(!new_L.empty()) L_idx = std::move(new_L); | ||
| if(!new_R.empty()) R_idx = std::move(new_R); | ||
| } | ||
|
|
||
| // Fall-back: choose the highest-y candidate from each set. | ||
| size_t bestL = L_idx[0]; | ||
| for(size_t i = 1; i < L_idx.size(); ++i){ | ||
| if(left[L_idx[i]].y > left[bestL].y || | ||
| (left[L_idx[i]].y == left[bestL].y && left[L_idx[i]].x > left[bestL].x)){ | ||
| bestL = L_idx[i]; | ||
| } | ||
| } | ||
| size_t bestR = R_idx[0]; | ||
| for(size_t i = 1; i < R_idx.size(); ++i){ | ||
| if(right[R_idx[i]].y > right[bestR].y || | ||
| (right[R_idx[i]].y == right[bestR].y && right[R_idx[i]].x < right[bestR].x)){ | ||
| bestR = R_idx[i]; | ||
| } | ||
| } | ||
| return {bestL, bestR}; | ||
| } |
There was a problem hiding this comment.
upper_bridge_2d() can exit without finding a valid bridge: when slopes is empty or when pruning fails to reduce candidate counts, it breaks and then falls back to picking the highest-y point from each side. That pair is not guaranteed to satisfy the bridge condition (supporting line with all points at/below), which makes downstream pruning incorrect/unpredictable. Consider returning an optional/flag indicating “bridge not found” and skipping bridge-based pruning in that case, or implement the full KS pruning step that guarantees progress.
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
The repo has incremental and divide-and-conquer 3D convex hull implementations but lacks a marriage-before-conquest variant. This adds
MarriageBeforeConquestConvexHull<T>which identifies bridge structure before recursing (the "marriage" step), enabling interior point pruning to reduce sub-problem sizes.Algorithm
work_queueat shallow depthsIncrementalConvexHull(which provides Shewchuk-style adaptive orient3d)Changes
src/YgorMeshesConvexHull.h: Class declaration withcompute()/get_mesh()API matchingDivideAndConquerConvexHullsrc/YgorMeshesConvexHull.cc: Full implementation includingupper_bridge_2d()(KS86 §3), bridge-based + AABB pruning, parallelmbc_hull()recursion, explicit instantiations forfloat/doubletests2/YgorMeshesConvexHull.cc: 13 subcases mirroring existing D&C tests — tetrahedron, cube, octahedron, interior points, random points, icosahedron, coplanar degeneracy, large grid, duplicates, float specialization, cross-algorithm consistencyKnown limitation
The bridge-based pruning operates on a single 2D projection (xy-plane), which can be too aggressive for some 3D inputs — hull vertices that project below the bridge in xy but are extreme in z may be incorrectly pruned. A more robust approach would intersect pruning decisions across multiple projections.
💬 Send tasks to Copilot coding agent from Slack and Teams to turn conversations into code. Copilot posts an update in your thread when it's finished.