grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2024-12-28 15:59:45 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Added convex_hull()

Browse source
This commit is contained in:
Revar Desmera 2019-04-16 15:34:54 -07:00
parent fde93d9991
commit 96739c3ea0
5 changed files with 556 additions and 25 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

200
convex_hull.scad Normal file
View file

@ -0,0 +1,200 @@
//////////////////////////////////////////////////////////////////////
// LibFile: convex_hull.scad
// Functions to create 2D and 3D convex hulls.
// To use, add the following line to the beginning of your file:
// ```
// include <BOSL/convex_hull.scad>
// ```
// Derived from Linde's Hull:
// - https://github.com/openscad/scad-utils
//////////////////////////////////////////////////////////////////////
include <BOSL/math.scad>
// Section: Generalized Hull
// Function: convex_hull()
// Usage:
// convex_hull(points)
// Description:
// When given a list of 3D points, returns a list of faces for
// the minimal convex hull polyhedron of those points. Each face
// is a list of indexes into `points`.
// When given a list of 2D points, or 3D points that are all
// coplanar, returns a list of indices into `points` for the path
// that forms the minimal convex hull polygon of those points.
// Arguments:
// points = The list of points to find the minimal convex hull of.
function convex_hull(points) =
!(len(points) > 0) ? [] :
len(points[0]) == 2 ? convex_hull2d(points) :
len(points[0]) == 3 ? convex_hull3d(points) : [];
// Section: 2D Hull
// Function: convex_hull2d()
// Usage:
// convex_hull2d(points)
// Description:
// Takes a list of arbitrary 2D points, and finds the minimal convex
// hull polygon to enclose them. Returns a path as a list of indices
// into `points`.
function convex_hull2d(points) =
(len(points) < 3)? [] : let(
a=0, b=1,
c = _find_first_noncollinear([a,b], points, 2)
) (c == len(points))? _convex_hull_collinear(points) : let(
remaining = [ for (i = [2:len(points)-1]) if (i != c) i ],
ccw = triangle_area2d(points[a], points[b], points[c]) > 0,
polygon = ccw? [a,b,c] : [a,c,b]
) _convex_hull_iterative_2d(points, polygon, remaining);
// Adds the remaining points one by one to the convex hull
function _convex_hull_iterative_2d(points, polygon, remaining, _i=0) =
(_i >= len(remaining))? polygon : let (
// pick a point
i = remaining[_i],
// find the segments that are in conflict with the point (point not inside)
conflicts = _find_conflicting_segments(points, polygon, points[i])
// no conflicts, skip point and move on
) (len(conflicts) == 0)? _convex_hull_iterative_2d(points, polygon, remaining, _i+1) : let(
// find the first conflicting segment and the first not conflicting
// conflict will be sorted, if not wrapping around, do it the easy way
polygon = _remove_conflicts_and_insert_point(polygon, conflicts, i)
) _convex_hull_iterative_2d(points, polygon, remaining, _i+1);
function _find_first_noncollinear(line, points, i) =
(i>=len(points) || !collinear_indexed(points, line[0], line[1], i))? i :
_find_first_noncollinear(line, points, i+1);
function _find_conflicting_segments(points, polygon, point) = [
for (i = [0:len(polygon)-1]) let(
j = (i+1) % len(polygon),
p1 = points[polygon[i]],
p2 = points[polygon[j]],
area = triangle_area2d(p1, p2, point)
) if (area < 0) i
];
// remove the conflicting segments from the polygon
function _remove_conflicts_and_insert_point(polygon, conflicts, point) =
(conflicts[0] == 0)? let(
nonconflicting = [ for(i = [0:len(polygon)-1]) if (!in_list(i, conflicts)) i ],
new_indices = concat(nonconflicting, (nonconflicting[len(nonconflicting)-1]+1) % len(polygon)),
polygon = concat([ for (i = new_indices) polygon[i] ], point)
) polygon : let(
before_conflicts = [ for(i = [0:min(conflicts)]) polygon[i] ],
after_conflicts = (max(conflicts) >= (len(polygon)-1))? [] : [ for(i = [max(conflicts)+1:len(polygon)-1]) polygon[i] ],
polygon = concat(before_conflicts, point, after_conflicts)
) polygon;
// Section: 3D Hull
// Function: convex_hull3d()
// Usage:
// convex_hull3d(points)
// Description:
// Takes a list of arbitrary 3D points, and finds the minimal convex
// hull polyhedron to enclose them. Returns a list of faces, where
// each face is a list of indexes into the given `points` list.
// If all points passed to it are coplanar, then the return is the
// list of indices of points forming the minimal convex hull polygon.
function convex_hull3d(points) =
(len(points) < 3)? list_range(len(points)) : let (
// start with a single triangle
a=0, b=1, c=2,
plane = plane3pt_indexed(points, a, b, c),
d = _find_first_noncoplanar(plane, points, 3)
) (d == len(points))? /* all coplanar*/ let (
pts2d = [ for (p = points) xyz_to_planar(p, points[a], points[b], points[c]) ],
hull2d = convex_hull2d(pts2d)
) hull2d : let(
remaining = [for (i = [3:len(points)-1]) if (i != d) i],
// Build an initial tetrahedron.
// Swap b, c if d is in front of triangle t.
ifop = in_front_of_plane(plane, points[d]),
bc = ifop? [c,b] : [b,c],
b = bc[0],
c = bc[1],
triangles = [
[a,b,c],
[d,b,a],
[c,d,a],
[b,d,c]
],
// calculate the plane equations
planes = [ for (t = triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
) _convex_hull_iterative(points, triangles, planes, remaining);
// Adds the remaining points one by one to the convex hull
function _convex_hull_iterative(points, triangles, planes, remaining, _i=0) =
_i >= len(remaining) ? triangles :
let (
// pick a point
i = remaining[_i],
// find the triangles that are in conflict with the point (point not inside)
conflicts = _find_conflicts(points[i], planes),
// for all triangles that are in conflict, collect their halfedges
halfedges = [
for(c = conflicts, i = [0:2]) let(
j = (i+1)%3
) [triangles[c][i], triangles[c][j]]
],
// find the outer perimeter of the set of conflicting triangles
horizon = _remove_internal_edges(halfedges),
// generate a new triangle for each horizon halfedge together with the picked point i
new_triangles = [ for (h = horizon) concat(h,i) ],
// calculate the corresponding plane equations
new_planes = [ for (t = new_triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
) _convex_hull_iterative(
points,
// remove the conflicting triangles and add the new ones
concat(list_remove(triangles, conflicts), new_triangles),
concat(list_remove(planes, conflicts), new_planes),
remaining,
_i+1
);
function _convex_hull_collinear(points) =
let(
a = points[0],
n = points[1] - a,
points1d = [ for(p = points) (p-a)*n ],
min_i = min_index(points1d),
max_i = max_index(points1d)
) [min_i, max_i];
function _remove_internal_edges(halfedges) = [
for (h = halfedges)
if (!in_list(reverse(h), halfedges))
h
];
function _find_conflicts(point, planes) = [
for (i = [0:len(planes)-1])
if (in_front_of_plane(planes[i], point))
i
];
function _find_first_noncoplanar(plane, points, i) =
(i >= len(points) || !coplanar(plane, points[i]))? i :
_find_first_noncoplanar(plane, points, i+1);
// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

264
math.scad
View file

@ -44,6 +44,8 @@ include <compat.scad>
PHI = (1+sqrt(5))/2; // The golden ratio phi.
EPSILON = 1e-9; // A really small value useful in comparing FP numbers. ie: abs(a-b)<EPSILON
// Function: Cpi()
@ -94,6 +96,36 @@ function quantup(x,y) = ceil(x/y)*y;
function constrain(v, minval, maxval) = min(maxval, max(minval, v));
// Function: min_index()
// Usage:
// min_index(vals);
// Description:
// Returns the index of the minimal value in the given list.
function min_index(vals, _minval, _minidx, _i=0) =
_i>=len(vals)? _minidx :
min_index(
vals,
((_minval == undef || vals[_i] < _minval)? vals[_i] : _minval),
((_minval == undef || vals[_i] < _minval)? _i : _minidx),
_i+1
);
// Function: max_index()
// Usage:
// max_index(vals);
// Description:
// Returns the index of the maximum value in the given list.
function max_index(vals, _maxval, _maxidx, _i=0) =
_i>=len(vals)? _maxidx :
max_index(
vals,
((_maxval == undef || vals[_i] > _maxval)? vals[_i] : _maxval),
((_maxval == undef || vals[_i] > _maxval)? _i : _maxidx),
_i+1
);
// Function: posmod()
// Usage:
// posmod(x,m)
@ -788,6 +820,20 @@ function unique(arr) =
// Function: list_remove()
// Usage:
// list_remove(list, elements)
// Description:
// Remove all items from `list` whose indexes are in `elements`.
// Arguments:
// list = The list to remove items from.
// elements = The list of indexes of items to remove.
function list_remove(list, elements) = [
for (i = [0:len(list)-1]) if (!search(i, elements)) list[i]
];
// Internal. Not exposed.
function _array_dim_recurse(v) =
!is_list(v[0])? (
@ -1114,6 +1160,42 @@ function xy_to_polar(x,y=undef) = let(
) [norm([xx,yy]), atan2(yy,xx)];
// Function: xyz_to_planar()
// Usage:
// xyz_to_planar(point, a, b, c);
// Description:
// Given three points defining a plane, returns the projected planar
// [X,Y] coordinates of the closest point to a 3D `point`. The origin
// of the planar coordinate system [0,0] will be at point `a`, and the
// Y+ axis direction will be towards point `b`. This coordinate system
// can be useful in taking a set of nearly coplanar points, and converting
// them to a pure XY set of coordinates for manipulation, before convering
// them back to the original 3D plane.
function xyz_to_planar(point, a, b, c) = let(
u = normalize(b-a),
v = normalize(c-a),
n = normalize(cross(u,v)),
w = normalize(cross(n,u)),
relpoint = point-a
) [relpoint * w, relpoint * u];
// Function: planar_to_xyz()
// Usage:
// planar_to_xyz(point, a, b, c);
// Description:
// Given three points defining a plane, converts a planar [X,Y]
// coordinate to the actual corresponding 3D point on the plane.
// The origin of the planar coordinate system [0,0] will be at point
// `a`, and the Y+ axis direction will be towards point `b`.
function planar_to_xyz(point, a, b, c) = let(
u = normalize(b-a),
v = normalize(c-a),
n = normalize(cross(u,v)),
w = normalize(cross(n,u))
) a + point.x * w + point.y * u;
// Function: cylindrical_to_xyz()
// Usage:
// cylindrical_to_xyz(r, theta, z)
@ -1618,4 +1700,186 @@ function pointlist_bounds(pts) = [
];
// Function: triangle_area2d()
// Usage:
// triangle_area2d(a,b,c);
// Description:
// Returns the area of a triangle formed between three vertices.
// Result will be negative if the points are in clockwise order.
// Examples:
// triangle_area2d([0,0], [5,10], [10,0]); // Returns -50
// triangle_area2d([10,0], [5,10], [0,0]); // Returns 50
function triangle_area2d(a,b,c) =
(
a.x * (b.y - c.y) +
b.x * (c.y - a.y) +
c.x * (a.y - b.y)
) / 2;
// Function: right_of_line2d()
// Usage:
// right_of_line2d(line, pt)
// Description:
// Returns true if the given point is to the left of the given line.
// Arguments:
// line = A list of two points.
// pt = The point to test.
function right_of_line2d(line, pt) =
triangle_area2d(line[0], line[1], pt) < 0;
// Function: collinear()
// Usage:
// collinear(a, b, c, [eps]);
// Description:
// Returns true if three points are co-linear.
// Arguments:
// a = First point.
// b = Second point.
// c = Third point.
// eps = Acceptable max angle variance. Default: EPSILON (1e-9) degrees.
function collinear(a, b, c, eps=EPSILON) =
abs(vector_angle(b-a,c-a)) < eps;
// Function: collinear_indexed()
// Usage:
// collinear_indexed(points, a, b, c, [eps]);
// Description:
// Returns true if three points are co-linear.
// Arguments:
// points = A list of points.
// a = Index in `points` of first point.
// b = Index in `points` of second point.
// c = Index in `points` of third point.
// eps = Acceptable max angle variance. Default: EPSILON (1e-9) degrees.
function collinear_indexed(points, a, b, c, eps=EPSILON) =
let(
p1=points[a],
p2=points[b],
p3=points[c]
) abs(vector_angle(p2-p1,p3-p1)) < eps;
// Function: plane3pt()
// Usage:
// plane3pt(p1, p2, p3);
// Description:
// Generates the cartesian equation of a plane from three non-colinear points on the plane.
// Returns [A,B,C,D] where Ax+By+Cz+D=0 is the equation of a plane.
// Arguments:
// p1 = The first point on the plane.
// p2 = The second point on the plane.
// p3 = The third point on the plane.
function plane3pt(p1, p2, p3) =
let(normal = normalize(cross(p3-p1, p2-p1))) concat(normal, [normal*p1]);
// Function: plane3pt_indexed()
// Usage:
// plane3pt_indexed(points, i1, i2, i3);
// Description:
// Given a list of points, and the indexes of three of those points,
// generates the cartesian equation of a plane that those points all
// lie on. Requires that the three indexed points be non-collinear.
// Returns [A,B,C,D] where Ax+By+Cz+D=0 is the equation of a plane.
// Arguments:
// points = A list of points.
// i1 = The index into `points` of the first point on the plane.
// i2 = The index into `points` of the second point on the plane.
// i3 = The index into `points` of the third point on the plane.
function plane3pt_indexed(points, i1, i2, i3) =
let(
p1 = points[i1],
p2 = points[i2],
p3 = points[i3],
normal = normalize(cross(p3-p1, p2-p1))
) concat(normal, [normal*p1]);
// Function: distance_from_plane()
// Usage:
// distance_from_plane(plane, point)
// Description:
// Given a plane as [A,B,C,D] where the cartesian equation for that plane
// is Ax+By+Cz+D=0, determines how far from that plane the given point is.
// The returned distance will be positive if the point is in front of the
// plane; on the same side of the plane as the normal of that plane points
// towards. If the point is behind the plane, then the distance returned
// will be negative. The normal of the plane is the same as [A,B,C].
// Arguments:
// plane = The [A,B,C,D] values for the equation of the plane.
// point = The point to test.
function distance_from_plane(plane, point) =
[plane.x, plane.y, plane.z] * point - plane[3];
// Function: coplanar()
// Usage:
// coplanar(plane, point);
// Description:
// Given a plane as [A,B,C,D] where the cartesian equation for that plane
// is Ax+By+Cz+D=0, determines if the given point is on that plane.
// Returns true if the point is on that plane.
// Arguments:
// plane = The [A,B,C,D] values for the equation of the plane.
// point = The point to test.
function coplanar(plane, point) =
abs(distance_from_plane(plane, point)) <= EPSILON;
// Function: in_front_of_plane()
// Usage:
// in_front_of_plane(plane, point);
// Description:
// Given a plane as [A,B,C,D] where the cartesian equation for that plane
// is Ax+By+Cz+D=0, determines if the given point is on the side of that
// plane that the normal points towards. The normal of the plane is the
// same as [A,B,C].
// Arguments:
// plane = The [A,B,C,D] values for the equation of the plane.
// point = The point to test.
function in_front_of_plane(plane, point) =
distance_from_plane(plane, point) > EPSILON;
// Function: simplify_path()
// Description:
// Takes a path and removes unnecessary collinear points.
// Usage:
// simplify_path(path, [eps])
// Arguments:
// path = A list of 2D path points.
// eps = Largest angle variance allowed. Default: EPSILON (1-e9) degrees.
function simplify_path(path, eps=EPSILON, _a=0, _b=2, _acc=[]) =
(_b >= len(path))? concat([path[0]], _acc, [path[len(path)-1]]) :
simplify_path(
path, eps,
(collinear_indexed(path, _a, _b-1, _b, eps=eps)? _a : _b-1),
_b+1,
(collinear_indexed(path, _a, _b-1, _b, eps=eps)? _acc : concat(_acc, [path[_b-1]]))
);
// Function: simplify_path_indexed()
// Description:
// Takes a list of points, and a path as a list of indexes into `points`,
// and removes all path points that are unecessarily collinear.
// Usage:
// simplify_path_indexed(path, eps)
// Arguments:
// points = A list of points.
// path = A list of indexes into `points` that forms a path.
// eps = Largest angle variance allowed. Default: EPSILON (1-e9) degrees.
function simplify_path_indexed(points, path, eps=EPSILON, _a=0, _b=2, _acc=[]) =
(_b >= len(path))? concat([path[0]], _acc, [path[len(path)-1]]) :
simplify_path_indexed(
points, path, eps,
(collinear_indexed(points, path[_a], path[_b-1], path[_b], eps=eps)? _a : _b-1),
_b+1,
(collinear_indexed(points, path[_a], path[_b-1], path[_b], eps=eps)? _acc : concat(_acc, [path[_b-1]]))
);
// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

29
paths.scad
View file

@ -55,18 +55,7 @@ use <triangulation.scad>
// Arguments:
// path = A list of 2D path points.
// eps = Largest angle delta between segments to count as colinear. Default: 1e-6
function simplify2d_path(path, eps=1e-6) = concat(
[path[0]],
[
for (
i = [1:len(path)-2]
) let (
v1 = path[i] - path[i-1],
v2 = path[i+1] - path[i-1]
) if (abs(cross(v1,v2)) > eps) path[i]
],
[path[len(path)-1]]
);
function simplify2d_path(path, eps=1e-6) = simplify_path(path, eps=eps);
// Function: simplify3d_path()
@ -77,18 +66,7 @@ function simplify2d_path(path, eps=1e-6) = concat(
// Arguments:
// path = A list of 3D path points.
// eps = Largest angle delta between segments to count as colinear. Default: 1e-6
function simplify3d_path(path, eps=1e-6) = concat(
[path[0]],
[
for (
i = [1:len(path)-2]
) let (
v1 = path[i] - path[i-1],
v2 = path[i+1] - path[i-1]
) if (vector_angle(v1,v2) > eps) path[i]
],
[path[len(path)-1]]
);
function simplify3d_path(path, eps=1e-6) = simplify_path(path, eps=eps);
// Function: path_length()
@ -119,7 +97,8 @@ function path_length(path) =
// scale = [X,Y] scaling factors for each axis. Default: `[1,1]`
// Example(2D):
// trace_polyline(path2d_regular_ngon(n=12, r=50), N=1, showpts=true);
function path2d_regular_ngon(n=6, r=undef, d=undef, cp=[0,0], scale=[1,1]) = let(
function path2d_regular_ngon(n=6, r=undef, d=undef, cp=[0,0], scale=[1,1]) =
let(
rr=get_radius(r=r, d=d, dflt=100)
) [
for (i=[0:n-1])

73
tests/test_convex_hull.scad Normal file
View file

@ -0,0 +1,73 @@
include <BOSL/math.scad>
include <BOSL/convex_hull.scad>
testpoints_on_sphere = [ for(p =
[
[1,PHI,0], [-1,PHI,0], [1,-PHI,0], [-1,-PHI,0],
[0,1,PHI], [0,-1,PHI], [0,1,-PHI], [0,-1,-PHI],
[PHI,0,1], [-PHI,0,1], [PHI,0,-1], [-PHI,0,-1]
])
normalize(p)
];
testpoints_circular = [ for(a = [0:15:360-EPSILON]) [cos(a),sin(a)] ];
testpoints_coplanar = let(u = normalize([1,3,7]), v = normalize([-2,1,-2])) [ for(i = [1:10]) rands(-1,1,1)[0] * u + rands(-1,1,1)[0] * v ];
testpoints_collinear_2d = let(u = normalize([5,3])) [ for(i = [1:20]) rands(-1,1,1)[0] * u ];
testpoints_collinear_3d = let(u = normalize([5,3,-5])) [ for(i = [1:20]) rands(-1,1,1)[0] * u ];
testpoints2d = 20 * [for (i = [1:10]) concat(rands(-1,1,2))];
testpoints3d = 20 * [for (i = [1:50]) concat(rands(-1,1,3))];
// All points are on the sphere, no point should be red
translate([-50,0]) visualize_hull(20*testpoints_on_sphere);
// 2D points
translate([50,0]) visualize_hull(testpoints2d);
// All points on a circle, no point should be red
translate([0,50]) visualize_hull(20*testpoints_circular);
// All points 3d but collinear
translate([0,-50]) visualize_hull(20*testpoints_coplanar);
// Collinear
translate([50,50]) visualize_hull(20*testpoints_collinear_2d);
// Collinear
translate([-50,50]) visualize_hull(20*testpoints_collinear_3d);
// 3D points
visualize_hull(testpoints3d);
module visualize_hull(points) {
hull = convex_hull(points);
%if (len(hull) > 0 && is_list(hull[0]) && len(hull[0]) > 0)
polyhedron(points=points, faces = hull);
else
polyhedron(points=points, faces = [hull]);
for (i = [0:len(points)-1]) {
p = points[i];
$fn = 16;
translate(p) {
if (hull_contains_index(hull,i)) {
color("blue") sphere(1);
} else {
color("red") sphere(1);
}
}
}
function hull_contains_index(hull, index) =
search(index,hull,1,0) ||
search(index,hull,1,1) ||
search(index,hull,1,2);
}
// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

15
tests/test_math.scad
View file

@ -630,6 +630,21 @@ module test_rotate_points3d() {
test_rotate_points3d();
module test_simplify_path()
{
path = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
assert(simplify_path(path) == [[-20,10],[-10,0],[10,0], [10,10]]);
}
test_simplify_path();
module test_simplify_path_indexed()
{
points = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
path = list_range(len(points));
assert(simplify_path_indexed(points, path) == [0,1,5,6]);
}
test_simplify_path_indexed();
// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap