From 96739c3ea038703e582b3170ceca4ad5dead7173 Mon Sep 17 00:00:00 2001 From: Revar Desmera Date: Tue, 16 Apr 2019 15:34:54 -0700 Subject: [PATCH] Added convex_hull() --- convex_hull.scad | 200 +++++++++++++++++++++++++++ math.scad | 264 ++++++++++++++++++++++++++++++++++++ paths.scad | 29 +--- tests/test_convex_hull.scad | 73 ++++++++++ tests/test_math.scad | 15 ++ 5 files changed, 556 insertions(+), 25 deletions(-) create mode 100644 convex_hull.scad create mode 100644 tests/test_convex_hull.scad diff --git a/convex_hull.scad b/convex_hull.scad new file mode 100644 index 0000000..4b5270e --- /dev/null +++ b/convex_hull.scad @@ -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 +// ``` +// Derived from Linde's Hull: +// - https://github.com/openscad/scad-utils +////////////////////////////////////////////////////////////////////// + +include + + + +// 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 diff --git a/math.scad b/math.scad index 33fd301..54b664f 100644 --- a/math.scad +++ b/math.scad @@ -44,6 +44,8 @@ include 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)=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 diff --git a/paths.scad b/paths.scad index e64cb28..8dbf2fb 100644 --- a/paths.scad +++ b/paths.scad @@ -55,18 +55,7 @@ use // 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]) diff --git a/tests/test_convex_hull.scad b/tests/test_convex_hull.scad new file mode 100644 index 0000000..d43ffac --- /dev/null +++ b/tests/test_convex_hull.scad @@ -0,0 +1,73 @@ +include +include + + +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 diff --git a/tests/test_math.scad b/tests/test_math.scad index 99d39d2..dcbad7d 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -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