From 319ef14e6c319508133c6961f6c5b5d0406c9d24 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Mon, 21 Jun 2021 18:25:35 +0100 Subject: [PATCH 01/12] convex collision and distance --- geometry.scad | 417 ++++++++++++++++++++++++++++++--------- tests/test_geometry.scad | 68 +++++-- vnf.scad | 16 +- 3 files changed, 385 insertions(+), 116 deletions(-) diff --git a/geometry.scad b/geometry.scad index 24fbac6..36e30dc 100644 --- a/geometry.scad +++ b/geometry.scad @@ -19,15 +19,8 @@ // edge = Array of two points forming the line segment to test against. // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) function point_on_segment2d(point, edge, eps=EPSILON) = - assert( is_vector(point,2), "Invalid point." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - assert( _valid_line(edge,2,eps=eps), "Invalid segment." ) - let( dp = point-edge[0], - de = edge[1]-edge[0], - ne = norm(de) ) - ( dp*de >= -eps*ne ) - && ( (dp-de)*de <= eps*ne ) // point projects on the segment - && _dist2line(point-edge[0],unit(de))eps*max(norm(line[1]),norm(line[0])); //Internal function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); @@ -85,22 +78,57 @@ function collinear(a, b, c, eps=EPSILON) = : noncollinear_triple(points,error=false,eps=eps)==[]; -// Function: distance_from_line() +// Function: point_line_distance() // Usage: -// distance_from_line(line, pt); +// point_line_distance(line, pt); // Description: // Finds the perpendicular distance of a point `pt` from the line `line`. // Arguments: // line = A list of two points, defining a line that both are on. // pt = A point to find the distance of from the line. // Example: -// distance_from_line([[-10,0], [10,0]], [3,8]); // Returns: 8 -function distance_from_line(line, pt) = +// dist = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +function point_line_distance(pt, line) = assert( _valid_line(line) && is_vector(pt,len(line[0])), "Invalid line, invalid point or incompatible dimensions." ) _dist2line(pt-line[0],unit(line[1]-line[0])); +// Function: point_segment_distance() +// Usage: +// dist = point_segment_distance(pt, seg); +// Description: +// Returns the closest distance of the given point to the given line segment. +// Arguments: +// pt = The point to check the distance of. +// seg = The two points representing the line segment to check the distance of. +// Example: +// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5 +function point_segment_distance(pt, seg) = + assert( is_matrix(concat([pt],seg),3), + "Input should be a point and a valid segment with the dimension equal to the point." ) + norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) : + norm(pt-segment_closest_point(seg,pt)); + + +// Function: segment_distance() +// Usage: +// dist = segment_distance(seg1, seg2); +// Description: +// Returns the closest distance of the two given line segments. +// Arguments: +// seg1 = The list of two points representing the first line segment to check the distance of. +// seg2 = The list of two points representing the second line segment to check the distance of. +// Example: +// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5 +// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0 +function segment_distance(seg1, seg2) = + assert( is_matrix(concat(seg1,seg2),4), + "Inputs should be two valid segments." ) + convex_distance(seg1,seg2); + + // Function: line_normal() // Usage: // line_normal([P1,P2]) @@ -436,17 +464,9 @@ function ray_closest_point(ray,pt) = // color("blue") translate(pt) sphere(r=1,$fn=12); // color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = - assert(_valid_line(seg), "Invalid segment." ) - assert(len(pt)==len(seg[0]), "Incompatible dimensions." ) - approx(seg[0],seg[1])? seg[0] : - let( - seglen = norm(seg[1]-seg[0]), - segvec = (seg[1]-seg[0])/seglen, - projection = (pt-seg[0]) * segvec - ) - projection<=0 ? seg[0] : - projection>=seglen ? seg[1] : - seg[0] + projection*segvec; + assert( is_matrix(concat([pt],seg),3) , + "Invalid point or segment or incompatible dimensions." ) + pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0]; // Function: line_from_points() @@ -454,7 +474,7 @@ function segment_closest_point(seg,pt) = // line_from_points(points, [fast], [eps]); // Description: // Given a list of 2 or more collinear points, returns a line containing them. -// If `fast` is false and the points are coincident, then `undef` is returned. +// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: // points = The list of points to find the line through. @@ -464,7 +484,7 @@ function line_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=undef), "Improper point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) let( pb = furthest_point(points[0],points) ) - approx(norm(points[pb]-points[0]),0) ? undef : + norm(points[pb]-points[0])=0, "The tolerance should be a positive number." ) - assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.") + assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.") let( bounded = is_list(bounded)? bounded : [bounded, bounded], @@ -1182,7 +1202,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) assert(is_path(poly,dim=3), "Invalid polygon." ) assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).") - assert(_valid_line(line,dim=3,eps=eps), "Invalid line." ) + assert(_valid_line(line,dim=3,eps=eps), "Invalid 3d line." ) let( bounded = is_list(bounded)? bounded : [bounded, bounded], poly = deduplicate(poly), @@ -1310,7 +1330,7 @@ function points_on_plane(points, plane, eps=EPSILON) = // plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`. // point = The 3D point to test. function in_front_of_plane(plane, point) = - distance_from_plane(plane, point) > EPSILON; + point_plane_distance(plane, point) > EPSILON; @@ -1425,6 +1445,7 @@ module circle_2tangents(pt1, pt2, pt3, r, d, h, center=false) { } } + // Function&Module: circle_3points() // Usage: As Function // circ = circle_3points(pt1, pt2, pt3); @@ -1636,7 +1657,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // eps = epsilon used for identifying the case with one solution. Default: 1e-9 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) = let(r=get_radius(r=r,d=d,dflt=undef)) - assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") + assert(_valid_line(line,2), "Invalid 2d line.") assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_num(r) && r>0, "Radius must be positive") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition") @@ -1680,14 +1701,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = pb = points[b], nrm = norm(pa-pb) ) - approx(nrm, 0) + nrm <= eps*max(norm(pa),norm(pb)) ? assert(!error, "Cannot find three noncollinear points in pointlist.") [] : let( n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist)0 && len(pts[0])>0 , "Invalid pointlist." ) - let(ptsT = transpose(pts)) - [ - [for(row=ptsT) min(row)], - [for(row=ptsT) max(row)] - ]; - + assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." ) + let( + select = ident(len(pts[0])), + spread = [for(i=[0:len(pts[0])-1]) + let( spreadi = pts*select[i] ) + [min(spreadi), max(spreadi)] ] ) + transpose(spread); // Function: closest_point() // Usage: @@ -1747,7 +1768,7 @@ function furthest_point(pt, points) = // area = polygon_area(poly); // Description: // Given a 2D or 3D planar polygon, returns the area of that polygon. -// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is []. +// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon. // Arguments: // poly = Polygon to compute the area of. @@ -1759,53 +1780,16 @@ function polygon_area(poly, signed=false) = ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ) signed ? total : abs(total) : let( plane = plane_from_polygon(poly) ) - plane==[]? [] : + plane==[]? undef : let( n = plane_normal(plane), - total = sum([ - for(i=[1:1:len(poly)-2]) - let( - v1 = poly[i] - poly[0], - v2 = poly[i+1] - poly[0] - ) - cross(v1,v2) - ])* n/2 + total = sum([ for(i=[1:1:len(poly)-2]) + cross(poly[i]-poly[0], poly[i+1]-poly[0]) + ]) * n/2 ) signed ? total : abs(total); -// Function: is_convex_polygon() -// Usage: -// is_convex_polygon(poly); -// Description: -// Returns true if the given 2D or 3D polygon is convex. -// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. -// If the points are collinear an error is generated. -// Arguments: -// poly = Polygon to check. -// eps = Tolerance for the collinearity test. Default: EPSILON. -// Example: -// is_convex_polygon(circle(d=50)); // Returns: true -// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true -// Example: -// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; -// is_convex_polygon(spiral); // Returns: false -function is_convex_polygon(poly,eps=EPSILON) = - assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly), - p0 = poly[0] ) - assert( lp>=3 , "A polygon must have at least 3 points" ) - let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(p0)==2 - ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) - min(crosses) >=0 || max(crosses)<=0 - : let( prod = crosses*sum(crosses), - minc = min(prod), - maxc = max(prod) ) - assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) - minc>=0 || maxc<=0; - - // Function: polygon_shift() // Usage: // polygon_shift(poly, i); @@ -1960,7 +1944,6 @@ function centroid(poly, eps=EPSILON) = val[1]/val[0]/3; - // Function: point_in_polygon() // Usage: // point_in_polygon(point, poly, ) @@ -1972,9 +1955,9 @@ function centroid(poly, eps=EPSILON) = // Returns -1 if the point is outside the polygon. // Returns 0 if the point is on the boundary. // Returns 1 if the point lies in the interior. -// The polygon does not need to be simple: it can have self-intersections. +// The polygon does not need to be simple: it may have self-intersections. // But the polygon cannot have holes (it must be simply connected). -// Rounding error may give mixed results for points on or near the boundary. +// Rounding errors may give mixed results for points on or near the boundary. // Arguments: // point = The 2D point to check position of. // poly = The list of 2D path points forming the perimeter of the polygon. @@ -2067,7 +2050,7 @@ function ccw_polygon(poly) = // poly = The list of the path points for the perimeter of the polygon. function reverse_polygon(poly) = assert(is_path(poly), "Input should be a polygon") - let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]]; + [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; // Function: polygon_normal() @@ -2075,7 +2058,7 @@ function reverse_polygon(poly) = // n = polygon_normal(poly); // Description: // Given a 3D planar polygon, returns a unit-length normal vector for the -// clockwise orientation of the polygon. If the polygon points are collinear, returns []. +// clockwise orientation of the polygon. If the polygon points are collinear, returns undef. // It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. @@ -2083,7 +2066,7 @@ function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? [] : + triple==[] ? undef : point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; @@ -2236,5 +2219,253 @@ function split_polygons_at_each_z(polys, zs, _i=0) = ], zs, _i=_i+1 ); +// Section: Convex Sets + +// Function: is_convex_polygon() +// Usage: +// is_convex_polygon(poly); +// Description: +// Returns true if the given 2D or 3D polygon is convex. +// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. +// If the points are collinear an error is generated. +// Arguments: +// poly = Polygon to check. +// eps = Tolerance for the collinearity test. Default: EPSILON. +// Example: +// is_convex_polygon(circle(d=50)); // Returns: true +// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true +// Example: +// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; +// is_convex_polygon(spiral); // Returns: false +function is_convex_polygon(poly,eps=EPSILON) = + assert(is_path(poly), "The input should be a 2D or 3D polygon." ) + let( lp = len(poly) ) + assert( lp>=3 , "A polygon must have at least 3 points" ) + let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) + len(poly[0])==2 + ? assert( max(max(crosses),-min(crosses))>eps, "The points are collinear" ) + min(crosses) >=0 || max(crosses)<=0 + : let( prod = crosses*sum(crosses), + minc = min(prod), + maxc = max(prod) ) + assert( max(maxc,-minc)>eps, "The points are collinear" ) + minc>=0 || maxc<=0; + + +// Function: convex_distance() +// Usage: +// convex_distance(pts1, pts2,); +// See also: +// convex_collision +// Descrition: +// Returns the smallest distance between a point in convex hull of `points1` +// and a point in the convex hull of `points2`. All the points in the lists +// should have the same dimension, either 2D or 3D. +// A zero result means the hulls intercept whithin a tolerance `eps`. +// Arguments: +// points1 - first list of 2d or 3d points. +// points2 - second list of 2d or 3d points. +// eps - tolerance in distance evaluations. Default: EPSILON. +// Example(2D): +// pts1 = move([-3,0], p=square(3,center=true)); +// pts2 = rot(a=45, p=square(2,center=true)); +// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; +// polygon(pts1); +// polygon(pts2); +// polygon(pts3); +// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864 +// echo(convex_distance(pts2,pts3)); // Returns: 0 +// Example(3D): +// sphr1 = sphere(2,$fn=10); +// sphr2 = move([4,0,0], p=sphr1); +// sphr3 = move([4.5,0,0], p=sphr1); +// vnf_polyhedron(sphr1); +// vnf_polyhedron(sphr2); +// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0 +// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5 +function convex_distance(points1, points2, eps=EPSILON) = + assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), + "The input list should be a consistent non empty list of points of same dimension.") + assert(len(points1[0])==2 || len(points1[0])==3 , + "The input points should be 2d or 3d points.") + let( d = points1[0]-points2[0] ) + norm(d)); +// See also: +// convex_distance +// Descrition: +// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2` +// otherwise, `false`. +// All the points in the lists should have the same dimension, either 2D or 3D. +// This function is tipically faster than `convex_distance` to find a non-collision. +// Arguments: +// points1 - first list of 2d or 3d points. +// points2 - second list of 2d or 3d points. +// eps - tolerance for the intersection tests. Default: EPSILON. +// Example(2D): +// pts1 = move([-3,0], p=square(3,center=true)); +// pts2 = rot(a=45, p=square(2,center=true)); +// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; +// polygon(pts1); +// polygon(pts2); +// polygon(pts3); +// echo(convex_collision(pts1,pts2)); // Returns: false +// echo(convex_collision(pts2,pts3)); // Returns: true +// Example(3D): +// sphr1 = sphere(2,$fn=10); +// sphr2 = move([4,0,0], p=sphr1); +// sphr3 = move([4.5,0,0], p=sphr1); +// vnf_polyhedron(sphr1); +// vnf_polyhedron(sphr2); +// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true +// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false +// +function convex_collision(points1, points2, eps=EPSILON) = + assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), + "The input list should be a consistent non empty list of points of same dimension.") + assert(len(points1[0])==2 || len(points1[0])==3 , + "The input points should be 2d or 3d points.") + let( d = points1[0]-points2[0] ) + norm(d) eps ? false : // no collision + let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) + _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps); + + +// given a simplex s, returns a pair: +// - the point of the s closest to the origin +// - the smallest sub-simplex of s that contains that point +function _closest_simplex(s,eps=EPSILON) = + assert(len(s)>=2 && len(s)<=4, "Internal error.") + len(s)==2 ? _closest_s1(s,eps) : + len(s)==3 ? _closest_s2(s,eps) + : _closest_s3(s,eps); + + +// find the closest to a 1-simplex +// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf +function _closest_s1(s,eps=EPSILON) = + norm(s[1]-s[0])1 ? [ s[1], [s[1]] ] : + [ s[0]+t*c, s ]; + + +// find the closest to a 2-simplex +// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf +function _closest_s2(s,eps=EPSILON) = + let( + dim = len(s[0]), + a = dim==3 ? s[0]: [ each s[0], 0] , + b = dim==3 ? s[1]: [ each s[1], 0] , + c = dim==3 ? s[2]: [ each s[2], 0] , + ab = norm(a-b), + bc = norm(b-c), + ca = norm(c-a), + nr = cross(b-a,c-a) + ) + norm(nr) <= eps*max(ab,bc,ca) // degenerate case + ? let( i = max_index([ab, bc, ca]) ) + _closest_s1([s[i],s[(i+1)%3]],eps) +// considering that s[2] was the last inserted vertex in s, +// the only possible outcomes are : +// s, [s[0],s[2]] and [s[1],s[2]] + : let( + class = (cross(nr,a-b)*a<0 ? 1 : 0 ) + + (cross(nr,c-a)*a<0 ? 2 : 0 ) + + (cross(nr,b-c)*b<0 ? 4 : 0 ) + ) + assert( class!=1, "Internal error" ) + class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri +// class==1 ? _closest_s1([s[0],s[1]]) : + class==2 ? _closest_s1([s[0],s[2]],eps) : + class==4 ? _closest_s1([s[1],s[2]],eps) : +// class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) : + class==3 ? _closest_s1([s[0],s[2]],eps) : +// class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) : + class==5 ? _closest_s1([s[1],s[2]],eps) : + c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps); + + +// find the closest to a 3-simplex +// it seems that degenerate 3-simplices are correctly manage without extra code +function _closest_s3(s,eps=EPSILON) = + assert( len(s[0])==3 && len(s)==4, "Internal error." ) + let( nr = cross(s[1]-s[0],s[2]-s[0]), + sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] ) + norm(nr)0)==(nrm*s[i]<0) ) i ] + ) + len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex + len(facing)==1 ? _closest_s2(tris[facing[0]], eps) : + let( // look for the origin-facing tri closest to the origin + closest = [for(i=facing) _closest_s2(tris[i], eps) ], + dist = [for(cl=closest) norm(cl[0]) ], + nearest = min_index(dist) + ) + closest[nearest]; + + +function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]); + + +function _support_diff(p1,p2,d) = + let( p1d = p1*d, p2d = p2*d ) + p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]]; + + + // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 04f5126..31fc8ac 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -9,7 +9,9 @@ include <../std.scad> test_point_on_segment2d(); test_point_left_of_line2d(); test_collinear(); -test_distance_from_line(); +test_point_line_distance(); +test_point_segment_distance(); +test_segment_distance(); test_line_normal(); test_line_intersection(); //test_line_ray_intersection(); @@ -44,7 +46,7 @@ test_plane_normal(); test_plane_offset(); test_projection_on_plane(); test_plane_point_nearest_origin(); -test_distance_from_plane(); +test_point_plane_distance(); test__general_plane_line_intersection(); test_plane_line_angle(); @@ -88,7 +90,7 @@ test_cleanup_path(); test_simplify_path(); test_simplify_path_indexed(); test_is_region(); - +test_convex_distance(); // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -230,7 +232,7 @@ module test__general_plane_line_intersection() { interspoint = line1[0]+inters1[1]*(line1[1]-line1[0]); assert_approx(inters1[0],interspoint, info1); assert_approx(point3d(plane1)*inters1[0], plane1[3], info1); // interspoint on the plane - assert_approx(distance_from_plane(plane1, inters1[0]), 0, info1); // inters1[0] on the plane + assert_approx(point_plane_distance(plane1, inters1[0]), 0, info1); // inters1[0] on the plane } // line parallel to the plane @@ -351,13 +353,35 @@ module test_collinear() { *test_collinear(); -module test_distance_from_line() { - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,1,1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [-1,-1,-1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON); +module test_point_line_distance() { + assert_approx(point_line_distance([1,1,1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([-1,-1,-1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([1,-1,0], [[-10,-10,-10], [10,10,10]]), sqrt(2)); + assert_approx(point_line_distance([8,-8,0], [[-10,-10,-10], [10,10,10]]), 8*sqrt(2)); } -*test_distance_from_line(); +*test_point_line_distance(); + + +module test_point_segment_distance() { + assert_approx(point_segment_distance([3,8], [[-10,0], [10,0]]), 8); + assert_approx(point_segment_distance([14,3], [[-10,0], [10,0]]), 5); +} +*test_point_segment_distance(); + + +module test_segment_distance() { + assert_approx(segment_distance([[-14,3], [-14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-15,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,-3], [14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-14,3], [-14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-6,5], [4,-5]], [[-10,0], [10,0]]), 0); + assert_approx(segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]), 0); +} +*test_segment_distance(); module test_line_normal() { @@ -713,12 +737,12 @@ module test_plane_normal() { *test_plane_normal(); -module test_distance_from_plane() { +module test_point_plane_distance() { plane1 = plane3pt([-10,0,0], [0,10,0], [10,0,0]); - assert(distance_from_plane(plane1, [0,0,5]) == 5); - assert(distance_from_plane(plane1, [5,5,8]) == 8); + assert(point_plane_distance(plane1, [0,0,5]) == 5); + assert(point_plane_distance(plane1, [5,5,8]) == 8); } -*test_distance_from_plane(); +*test_point_plane_distance(); module test_polygon_line_intersection() { @@ -1051,6 +1075,20 @@ module test_is_region() { } *test_is_region(); - +module test_convex_distance() { + c1 = circle(10,$fn=24); + c2 = move([15,0], p=c1); + assert(convex_distance(c1, c2)==0); + c3 = move([22,0],c1); + assert(abs(convex_distance(c1, c3)-2)2) + _vnf_validate_err("MULTCONN", [for (i=uniq_edges[i]) varr[i]]) + ]), + issues = concat(issues, multconn_edges) + ) multconn_edges? issues : let( repeated_faces = [ for (i=idx(dfaces), j=idx(dfaces)) @@ -864,14 +872,6 @@ function vnf_validate(vnf, show_warns=true, check_isects=false) = ], issues = concat(issues, repeated_faces) ) repeated_faces? issues : - let( - multconn_edges = unique([ - for (i = idx(uniq_edges)) - if (edgecnts[1][i]>2) - _vnf_validate_err("MULTCONN", [for (i=uniq_edges[i]) varr[i]]) - ]), - issues = concat(issues, multconn_edges) - ) multconn_edges? issues : let( reversals = unique([ for(i = idx(dfaces), j = idx(dfaces)) if(i != j) From 3a857d89ecd25809deb1097b72a97a52a70a9955 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Mon, 21 Jun 2021 18:43:51 +0100 Subject: [PATCH 02/12] Revert "convex collision and distance" This reverts commit 319ef14e6c319508133c6961f6c5b5d0406c9d24. --- geometry.scad | 417 +++++++++------------------------------ tests/test_geometry.scad | 68 ++----- vnf.scad | 16 +- 3 files changed, 116 insertions(+), 385 deletions(-) diff --git a/geometry.scad b/geometry.scad index 36e30dc..24fbac6 100644 --- a/geometry.scad +++ b/geometry.scad @@ -19,8 +19,15 @@ // edge = Array of two points forming the line segment to test against. // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) function point_on_segment2d(point, edge, eps=EPSILON) = + assert( is_vector(point,2), "Invalid point." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - point_segment_distance(point, edge)= -eps*ne ) + && ( (dp-de)*de <= eps*ne ) // point projects on the segment + && _dist2line(point-edge[0],unit(de))eps*max(norm(line[1]),norm(line[0])); + && ! approx(norm(line[1]-line[0]), 0, eps); //Internal function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); @@ -78,57 +85,22 @@ function collinear(a, b, c, eps=EPSILON) = : noncollinear_triple(points,error=false,eps=eps)==[]; -// Function: point_line_distance() +// Function: distance_from_line() // Usage: -// point_line_distance(line, pt); +// distance_from_line(line, pt); // Description: // Finds the perpendicular distance of a point `pt` from the line `line`. // Arguments: // line = A list of two points, defining a line that both are on. // pt = A point to find the distance of from the line. // Example: -// dist = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 -function point_line_distance(pt, line) = +// distance_from_line([[-10,0], [10,0]], [3,8]); // Returns: 8 +function distance_from_line(line, pt) = assert( _valid_line(line) && is_vector(pt,len(line[0])), "Invalid line, invalid point or incompatible dimensions." ) _dist2line(pt-line[0],unit(line[1]-line[0])); -// Function: point_segment_distance() -// Usage: -// dist = point_segment_distance(pt, seg); -// Description: -// Returns the closest distance of the given point to the given line segment. -// Arguments: -// pt = The point to check the distance of. -// seg = The two points representing the line segment to check the distance of. -// Example: -// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 -// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5 -function point_segment_distance(pt, seg) = - assert( is_matrix(concat([pt],seg),3), - "Input should be a point and a valid segment with the dimension equal to the point." ) - norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) : - norm(pt-segment_closest_point(seg,pt)); - - -// Function: segment_distance() -// Usage: -// dist = segment_distance(seg1, seg2); -// Description: -// Returns the closest distance of the two given line segments. -// Arguments: -// seg1 = The list of two points representing the first line segment to check the distance of. -// seg2 = The list of two points representing the second line segment to check the distance of. -// Example: -// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5 -// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0 -function segment_distance(seg1, seg2) = - assert( is_matrix(concat(seg1,seg2),4), - "Inputs should be two valid segments." ) - convex_distance(seg1,seg2); - - // Function: line_normal() // Usage: // line_normal([P1,P2]) @@ -464,9 +436,17 @@ function ray_closest_point(ray,pt) = // color("blue") translate(pt) sphere(r=1,$fn=12); // color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = - assert( is_matrix(concat([pt],seg),3) , - "Invalid point or segment or incompatible dimensions." ) - pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0]; + assert(_valid_line(seg), "Invalid segment." ) + assert(len(pt)==len(seg[0]), "Incompatible dimensions." ) + approx(seg[0],seg[1])? seg[0] : + let( + seglen = norm(seg[1]-seg[0]), + segvec = (seg[1]-seg[0])/seglen, + projection = (pt-seg[0]) * segvec + ) + projection<=0 ? seg[0] : + projection>=seglen ? seg[1] : + seg[0] + projection*segvec; // Function: line_from_points() @@ -474,7 +454,7 @@ function segment_closest_point(seg,pt) = // line_from_points(points, [fast], [eps]); // Description: // Given a list of 2 or more collinear points, returns a line containing them. -// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned. +// If `fast` is false and the points are coincident, then `undef` is returned. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: // points = The list of points to find the line through. @@ -484,7 +464,7 @@ function line_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=undef), "Improper point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) let( pb = furthest_point(points[0],points) ) - norm(points[pb]-points[0])=0, "The tolerance should be a positive number." ) - assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.") + assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.") let( bounded = is_list(bounded)? bounded : [bounded, bounded], @@ -1202,7 +1182,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) assert(is_path(poly,dim=3), "Invalid polygon." ) assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).") - assert(_valid_line(line,dim=3,eps=eps), "Invalid 3d line." ) + assert(_valid_line(line,dim=3,eps=eps), "Invalid line." ) let( bounded = is_list(bounded)? bounded : [bounded, bounded], poly = deduplicate(poly), @@ -1330,7 +1310,7 @@ function points_on_plane(points, plane, eps=EPSILON) = // plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`. // point = The 3D point to test. function in_front_of_plane(plane, point) = - point_plane_distance(plane, point) > EPSILON; + distance_from_plane(plane, point) > EPSILON; @@ -1445,7 +1425,6 @@ module circle_2tangents(pt1, pt2, pt3, r, d, h, center=false) { } } - // Function&Module: circle_3points() // Usage: As Function // circ = circle_3points(pt1, pt2, pt3); @@ -1657,7 +1636,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // eps = epsilon used for identifying the case with one solution. Default: 1e-9 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) = let(r=get_radius(r=r,d=d,dflt=undef)) - assert(_valid_line(line,2), "Invalid 2d line.") + assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_num(r) && r>0, "Radius must be positive") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition") @@ -1701,14 +1680,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = pb = points[b], nrm = norm(pa-pb) ) - nrm <= eps*max(norm(pa),norm(pb)) + approx(nrm, 0) ? assert(!error, "Cannot find three noncollinear points in pointlist.") [] : let( n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist) < eps*nrm + max(distlist)0 && len(pts[0])>0 , "Invalid pointlist." ) + let(ptsT = transpose(pts)) + [ + [for(row=ptsT) min(row)], + [for(row=ptsT) max(row)] + ]; + // Function: closest_point() // Usage: @@ -1768,7 +1747,7 @@ function furthest_point(pt, points) = // area = polygon_area(poly); // Description: // Given a 2D or 3D planar polygon, returns the area of that polygon. -// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`. +// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is []. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon. // Arguments: // poly = Polygon to compute the area of. @@ -1780,16 +1759,53 @@ function polygon_area(poly, signed=false) = ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ) signed ? total : abs(total) : let( plane = plane_from_polygon(poly) ) - plane==[]? undef : + plane==[]? [] : let( n = plane_normal(plane), - total = sum([ for(i=[1:1:len(poly)-2]) - cross(poly[i]-poly[0], poly[i+1]-poly[0]) - ]) * n/2 + total = sum([ + for(i=[1:1:len(poly)-2]) + let( + v1 = poly[i] - poly[0], + v2 = poly[i+1] - poly[0] + ) + cross(v1,v2) + ])* n/2 ) signed ? total : abs(total); +// Function: is_convex_polygon() +// Usage: +// is_convex_polygon(poly); +// Description: +// Returns true if the given 2D or 3D polygon is convex. +// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. +// If the points are collinear an error is generated. +// Arguments: +// poly = Polygon to check. +// eps = Tolerance for the collinearity test. Default: EPSILON. +// Example: +// is_convex_polygon(circle(d=50)); // Returns: true +// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true +// Example: +// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; +// is_convex_polygon(spiral); // Returns: false +function is_convex_polygon(poly,eps=EPSILON) = + assert(is_path(poly), "The input should be a 2D or 3D polygon." ) + let( lp = len(poly), + p0 = poly[0] ) + assert( lp>=3 , "A polygon must have at least 3 points" ) + let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) + len(p0)==2 + ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) + min(crosses) >=0 || max(crosses)<=0 + : let( prod = crosses*sum(crosses), + minc = min(prod), + maxc = max(prod) ) + assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) + minc>=0 || maxc<=0; + + // Function: polygon_shift() // Usage: // polygon_shift(poly, i); @@ -1944,6 +1960,7 @@ function centroid(poly, eps=EPSILON) = val[1]/val[0]/3; + // Function: point_in_polygon() // Usage: // point_in_polygon(point, poly, ) @@ -1955,9 +1972,9 @@ function centroid(poly, eps=EPSILON) = // Returns -1 if the point is outside the polygon. // Returns 0 if the point is on the boundary. // Returns 1 if the point lies in the interior. -// The polygon does not need to be simple: it may have self-intersections. +// The polygon does not need to be simple: it can have self-intersections. // But the polygon cannot have holes (it must be simply connected). -// Rounding errors may give mixed results for points on or near the boundary. +// Rounding error may give mixed results for points on or near the boundary. // Arguments: // point = The 2D point to check position of. // poly = The list of 2D path points forming the perimeter of the polygon. @@ -2050,7 +2067,7 @@ function ccw_polygon(poly) = // poly = The list of the path points for the perimeter of the polygon. function reverse_polygon(poly) = assert(is_path(poly), "Input should be a polygon") - [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; + let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]]; // Function: polygon_normal() @@ -2058,7 +2075,7 @@ function reverse_polygon(poly) = // n = polygon_normal(poly); // Description: // Given a 3D planar polygon, returns a unit-length normal vector for the -// clockwise orientation of the polygon. If the polygon points are collinear, returns undef. +// clockwise orientation of the polygon. If the polygon points are collinear, returns []. // It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. @@ -2066,7 +2083,7 @@ function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? undef : + triple==[] ? [] : point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; @@ -2219,253 +2236,5 @@ function split_polygons_at_each_z(polys, zs, _i=0) = ], zs, _i=_i+1 ); -// Section: Convex Sets - -// Function: is_convex_polygon() -// Usage: -// is_convex_polygon(poly); -// Description: -// Returns true if the given 2D or 3D polygon is convex. -// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. -// If the points are collinear an error is generated. -// Arguments: -// poly = Polygon to check. -// eps = Tolerance for the collinearity test. Default: EPSILON. -// Example: -// is_convex_polygon(circle(d=50)); // Returns: true -// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true -// Example: -// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; -// is_convex_polygon(spiral); // Returns: false -function is_convex_polygon(poly,eps=EPSILON) = - assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly) ) - assert( lp>=3 , "A polygon must have at least 3 points" ) - let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(poly[0])==2 - ? assert( max(max(crosses),-min(crosses))>eps, "The points are collinear" ) - min(crosses) >=0 || max(crosses)<=0 - : let( prod = crosses*sum(crosses), - minc = min(prod), - maxc = max(prod) ) - assert( max(maxc,-minc)>eps, "The points are collinear" ) - minc>=0 || maxc<=0; - - -// Function: convex_distance() -// Usage: -// convex_distance(pts1, pts2,); -// See also: -// convex_collision -// Descrition: -// Returns the smallest distance between a point in convex hull of `points1` -// and a point in the convex hull of `points2`. All the points in the lists -// should have the same dimension, either 2D or 3D. -// A zero result means the hulls intercept whithin a tolerance `eps`. -// Arguments: -// points1 - first list of 2d or 3d points. -// points2 - second list of 2d or 3d points. -// eps - tolerance in distance evaluations. Default: EPSILON. -// Example(2D): -// pts1 = move([-3,0], p=square(3,center=true)); -// pts2 = rot(a=45, p=square(2,center=true)); -// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; -// polygon(pts1); -// polygon(pts2); -// polygon(pts3); -// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864 -// echo(convex_distance(pts2,pts3)); // Returns: 0 -// Example(3D): -// sphr1 = sphere(2,$fn=10); -// sphr2 = move([4,0,0], p=sphr1); -// sphr3 = move([4.5,0,0], p=sphr1); -// vnf_polyhedron(sphr1); -// vnf_polyhedron(sphr2); -// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0 -// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5 -function convex_distance(points1, points2, eps=EPSILON) = - assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), - "The input list should be a consistent non empty list of points of same dimension.") - assert(len(points1[0])==2 || len(points1[0])==3 , - "The input points should be 2d or 3d points.") - let( d = points1[0]-points2[0] ) - norm(d)); -// See also: -// convex_distance -// Descrition: -// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2` -// otherwise, `false`. -// All the points in the lists should have the same dimension, either 2D or 3D. -// This function is tipically faster than `convex_distance` to find a non-collision. -// Arguments: -// points1 - first list of 2d or 3d points. -// points2 - second list of 2d or 3d points. -// eps - tolerance for the intersection tests. Default: EPSILON. -// Example(2D): -// pts1 = move([-3,0], p=square(3,center=true)); -// pts2 = rot(a=45, p=square(2,center=true)); -// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; -// polygon(pts1); -// polygon(pts2); -// polygon(pts3); -// echo(convex_collision(pts1,pts2)); // Returns: false -// echo(convex_collision(pts2,pts3)); // Returns: true -// Example(3D): -// sphr1 = sphere(2,$fn=10); -// sphr2 = move([4,0,0], p=sphr1); -// sphr3 = move([4.5,0,0], p=sphr1); -// vnf_polyhedron(sphr1); -// vnf_polyhedron(sphr2); -// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true -// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false -// -function convex_collision(points1, points2, eps=EPSILON) = - assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), - "The input list should be a consistent non empty list of points of same dimension.") - assert(len(points1[0])==2 || len(points1[0])==3 , - "The input points should be 2d or 3d points.") - let( d = points1[0]-points2[0] ) - norm(d) eps ? false : // no collision - let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) - _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps); - - -// given a simplex s, returns a pair: -// - the point of the s closest to the origin -// - the smallest sub-simplex of s that contains that point -function _closest_simplex(s,eps=EPSILON) = - assert(len(s)>=2 && len(s)<=4, "Internal error.") - len(s)==2 ? _closest_s1(s,eps) : - len(s)==3 ? _closest_s2(s,eps) - : _closest_s3(s,eps); - - -// find the closest to a 1-simplex -// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf -function _closest_s1(s,eps=EPSILON) = - norm(s[1]-s[0])1 ? [ s[1], [s[1]] ] : - [ s[0]+t*c, s ]; - - -// find the closest to a 2-simplex -// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf -function _closest_s2(s,eps=EPSILON) = - let( - dim = len(s[0]), - a = dim==3 ? s[0]: [ each s[0], 0] , - b = dim==3 ? s[1]: [ each s[1], 0] , - c = dim==3 ? s[2]: [ each s[2], 0] , - ab = norm(a-b), - bc = norm(b-c), - ca = norm(c-a), - nr = cross(b-a,c-a) - ) - norm(nr) <= eps*max(ab,bc,ca) // degenerate case - ? let( i = max_index([ab, bc, ca]) ) - _closest_s1([s[i],s[(i+1)%3]],eps) -// considering that s[2] was the last inserted vertex in s, -// the only possible outcomes are : -// s, [s[0],s[2]] and [s[1],s[2]] - : let( - class = (cross(nr,a-b)*a<0 ? 1 : 0 ) - + (cross(nr,c-a)*a<0 ? 2 : 0 ) - + (cross(nr,b-c)*b<0 ? 4 : 0 ) - ) - assert( class!=1, "Internal error" ) - class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri -// class==1 ? _closest_s1([s[0],s[1]]) : - class==2 ? _closest_s1([s[0],s[2]],eps) : - class==4 ? _closest_s1([s[1],s[2]],eps) : -// class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) : - class==3 ? _closest_s1([s[0],s[2]],eps) : -// class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) : - class==5 ? _closest_s1([s[1],s[2]],eps) : - c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps); - - -// find the closest to a 3-simplex -// it seems that degenerate 3-simplices are correctly manage without extra code -function _closest_s3(s,eps=EPSILON) = - assert( len(s[0])==3 && len(s)==4, "Internal error." ) - let( nr = cross(s[1]-s[0],s[2]-s[0]), - sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] ) - norm(nr)0)==(nrm*s[i]<0) ) i ] - ) - len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex - len(facing)==1 ? _closest_s2(tris[facing[0]], eps) : - let( // look for the origin-facing tri closest to the origin - closest = [for(i=facing) _closest_s2(tris[i], eps) ], - dist = [for(cl=closest) norm(cl[0]) ], - nearest = min_index(dist) - ) - closest[nearest]; - - -function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]); - - -function _support_diff(p1,p2,d) = - let( p1d = p1*d, p2d = p2*d ) - p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]]; - - - // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 31fc8ac..04f5126 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -9,9 +9,7 @@ include <../std.scad> test_point_on_segment2d(); test_point_left_of_line2d(); test_collinear(); -test_point_line_distance(); -test_point_segment_distance(); -test_segment_distance(); +test_distance_from_line(); test_line_normal(); test_line_intersection(); //test_line_ray_intersection(); @@ -46,7 +44,7 @@ test_plane_normal(); test_plane_offset(); test_projection_on_plane(); test_plane_point_nearest_origin(); -test_point_plane_distance(); +test_distance_from_plane(); test__general_plane_line_intersection(); test_plane_line_angle(); @@ -90,7 +88,7 @@ test_cleanup_path(); test_simplify_path(); test_simplify_path_indexed(); test_is_region(); -test_convex_distance(); + // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -232,7 +230,7 @@ module test__general_plane_line_intersection() { interspoint = line1[0]+inters1[1]*(line1[1]-line1[0]); assert_approx(inters1[0],interspoint, info1); assert_approx(point3d(plane1)*inters1[0], plane1[3], info1); // interspoint on the plane - assert_approx(point_plane_distance(plane1, inters1[0]), 0, info1); // inters1[0] on the plane + assert_approx(distance_from_plane(plane1, inters1[0]), 0, info1); // inters1[0] on the plane } // line parallel to the plane @@ -353,35 +351,13 @@ module test_collinear() { *test_collinear(); -module test_point_line_distance() { - assert_approx(point_line_distance([1,1,1], [[-10,-10,-10], [10,10,10]]), 0); - assert_approx(point_line_distance([-1,-1,-1], [[-10,-10,-10], [10,10,10]]), 0); - assert_approx(point_line_distance([1,-1,0], [[-10,-10,-10], [10,10,10]]), sqrt(2)); - assert_approx(point_line_distance([8,-8,0], [[-10,-10,-10], [10,10,10]]), 8*sqrt(2)); +module test_distance_from_line() { + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,1,1])) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [-1,-1,-1])) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON); } -*test_point_line_distance(); - - -module test_point_segment_distance() { - assert_approx(point_segment_distance([3,8], [[-10,0], [10,0]]), 8); - assert_approx(point_segment_distance([14,3], [[-10,0], [10,0]]), 5); -} -*test_point_segment_distance(); - - -module test_segment_distance() { - assert_approx(segment_distance([[-14,3], [-14,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,3], [14,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,-3], [-14,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,-3], [-15,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,-3], [14,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,3], [14,-3]], [[-10,0], [10,0]]), 4); - assert_approx(segment_distance([[-14,3], [-14,-3]], [[-10,0], [10,0]]), 4); - assert_approx(segment_distance([[-6,5], [4,-5]], [[-10,0], [10,0]]), 0); - assert_approx(segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]), 0); -} -*test_segment_distance(); +*test_distance_from_line(); module test_line_normal() { @@ -737,12 +713,12 @@ module test_plane_normal() { *test_plane_normal(); -module test_point_plane_distance() { +module test_distance_from_plane() { plane1 = plane3pt([-10,0,0], [0,10,0], [10,0,0]); - assert(point_plane_distance(plane1, [0,0,5]) == 5); - assert(point_plane_distance(plane1, [5,5,8]) == 8); + assert(distance_from_plane(plane1, [0,0,5]) == 5); + assert(distance_from_plane(plane1, [5,5,8]) == 8); } -*test_point_plane_distance(); +*test_distance_from_plane(); module test_polygon_line_intersection() { @@ -1075,20 +1051,6 @@ module test_is_region() { } *test_is_region(); -module test_convex_distance() { - c1 = circle(10,$fn=24); - c2 = move([15,0], p=c1); - assert(convex_distance(c1, c2)==0); - c3 = move([22,0],c1); - assert(abs(convex_distance(c1, c3)-2)2) - _vnf_validate_err("MULTCONN", [for (i=uniq_edges[i]) varr[i]]) - ]), - issues = concat(issues, multconn_edges) - ) multconn_edges? issues : let( repeated_faces = [ for (i=idx(dfaces), j=idx(dfaces)) @@ -872,6 +864,14 @@ function vnf_validate(vnf, show_warns=true, check_isects=false) = ], issues = concat(issues, repeated_faces) ) repeated_faces? issues : + let( + multconn_edges = unique([ + for (i = idx(uniq_edges)) + if (edgecnts[1][i]>2) + _vnf_validate_err("MULTCONN", [for (i=uniq_edges[i]) varr[i]]) + ]), + issues = concat(issues, multconn_edges) + ) multconn_edges? issues : let( reversals = unique([ for(i = idx(dfaces), j = idx(dfaces)) if(i != j) From eccb006f850b9becb816b0133aecc5a3c378bdf0 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Mon, 21 Jun 2021 19:08:13 +0100 Subject: [PATCH 03/12] Introduces convex collision and distance --- geometry.scad | 418 ++++++++++++++++++++++++++++++--------- tests/test_geometry.scad | 68 +++++-- 2 files changed, 380 insertions(+), 106 deletions(-) diff --git a/geometry.scad b/geometry.scad index 24fbac6..0221a2b 100644 --- a/geometry.scad +++ b/geometry.scad @@ -19,15 +19,8 @@ // edge = Array of two points forming the line segment to test against. // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) function point_on_segment2d(point, edge, eps=EPSILON) = - assert( is_vector(point,2), "Invalid point." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - assert( _valid_line(edge,2,eps=eps), "Invalid segment." ) - let( dp = point-edge[0], - de = edge[1]-edge[0], - ne = norm(de) ) - ( dp*de >= -eps*ne ) - && ( (dp-de)*de <= eps*ne ) // point projects on the segment - && _dist2line(point-edge[0],unit(de))eps*max(norm(line[1]),norm(line[0])); //Internal function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); @@ -85,22 +78,57 @@ function collinear(a, b, c, eps=EPSILON) = : noncollinear_triple(points,error=false,eps=eps)==[]; -// Function: distance_from_line() +// Function: point_line_distance() // Usage: -// distance_from_line(line, pt); +// point_line_distance(line, pt); // Description: // Finds the perpendicular distance of a point `pt` from the line `line`. // Arguments: // line = A list of two points, defining a line that both are on. // pt = A point to find the distance of from the line. // Example: -// distance_from_line([[-10,0], [10,0]], [3,8]); // Returns: 8 -function distance_from_line(line, pt) = +// dist = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +function point_line_distance(pt, line) = assert( _valid_line(line) && is_vector(pt,len(line[0])), "Invalid line, invalid point or incompatible dimensions." ) _dist2line(pt-line[0],unit(line[1]-line[0])); +// Function: point_segment_distance() +// Usage: +// dist = point_segment_distance(pt, seg); +// Description: +// Returns the closest distance of the given point to the given line segment. +// Arguments: +// pt = The point to check the distance of. +// seg = The two points representing the line segment to check the distance of. +// Example: +// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5 +function point_segment_distance(pt, seg) = + assert( is_matrix(concat([pt],seg),3), + "Input should be a point and a valid segment with the dimension equal to the point." ) + norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) : + norm(pt-segment_closest_point(seg,pt)); + + +// Function: segment_distance() +// Usage: +// dist = segment_distance(seg1, seg2); +// Description: +// Returns the closest distance of the two given line segments. +// Arguments: +// seg1 = The list of two points representing the first line segment to check the distance of. +// seg2 = The list of two points representing the second line segment to check the distance of. +// Example: +// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5 +// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0 +function segment_distance(seg1, seg2) = + assert( is_matrix(concat(seg1,seg2),4), + "Inputs should be two valid segments." ) + convex_distance(seg1,seg2); + + // Function: line_normal() // Usage: // line_normal([P1,P2]) @@ -436,17 +464,9 @@ function ray_closest_point(ray,pt) = // color("blue") translate(pt) sphere(r=1,$fn=12); // color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = - assert(_valid_line(seg), "Invalid segment." ) - assert(len(pt)==len(seg[0]), "Incompatible dimensions." ) - approx(seg[0],seg[1])? seg[0] : - let( - seglen = norm(seg[1]-seg[0]), - segvec = (seg[1]-seg[0])/seglen, - projection = (pt-seg[0]) * segvec - ) - projection<=0 ? seg[0] : - projection>=seglen ? seg[1] : - seg[0] + projection*segvec; + assert( is_matrix(concat([pt],seg),3) , + "Invalid point or segment or incompatible dimensions." ) + pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0]; // Function: line_from_points() @@ -454,7 +474,7 @@ function segment_closest_point(seg,pt) = // line_from_points(points, [fast], [eps]); // Description: // Given a list of 2 or more collinear points, returns a line containing them. -// If `fast` is false and the points are coincident, then `undef` is returned. +// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: // points = The list of points to find the line through. @@ -464,7 +484,7 @@ function line_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=undef), "Improper point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) let( pb = furthest_point(points[0],points) ) - approx(norm(points[pb]-points[0]),0) ? undef : + norm(points[pb]-points[0])=0, "The tolerance should be a positive number." ) - assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.") + assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.") let( bounded = is_list(bounded)? bounded : [bounded, bounded], @@ -1182,7 +1202,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) assert(is_path(poly,dim=3), "Invalid polygon." ) assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).") - assert(_valid_line(line,dim=3,eps=eps), "Invalid line." ) + assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." ) let( bounded = is_list(bounded)? bounded : [bounded, bounded], poly = deduplicate(poly), @@ -1310,7 +1330,7 @@ function points_on_plane(points, plane, eps=EPSILON) = // plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`. // point = The 3D point to test. function in_front_of_plane(plane, point) = - distance_from_plane(plane, point) > EPSILON; + point_plane_distance(plane, point) > EPSILON; @@ -1636,7 +1656,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // eps = epsilon used for identifying the case with one solution. Default: 1e-9 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) = let(r=get_radius(r=r,d=d,dflt=undef)) - assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") + assert(_valid_line(line,2), "Invalid 2d line.") assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_num(r) && r>0, "Radius must be positive") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition") @@ -1680,14 +1700,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = pb = points[b], nrm = norm(pa-pb) ) - approx(nrm, 0) + nrm <= eps*max(norm(pa),norm(pb)) ? assert(!error, "Cannot find three noncollinear points in pointlist.") [] : let( n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist)0 && len(pts[0])>0 , "Invalid pointlist." ) - let(ptsT = transpose(pts)) - [ - [for(row=ptsT) min(row)], - [for(row=ptsT) max(row)] - ]; + assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." ) + let( + select = ident(len(pts[0])), + spread = [for(i=[0:len(pts[0])-1]) + let( spreadi = pts*select[i] ) + [min(spreadi), max(spreadi)] ] ) + transpose(spread); // Function: closest_point() @@ -1747,7 +1768,7 @@ function furthest_point(pt, points) = // area = polygon_area(poly); // Description: // Given a 2D or 3D planar polygon, returns the area of that polygon. -// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is []. +// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon. // Arguments: // poly = Polygon to compute the area of. @@ -1759,53 +1780,17 @@ function polygon_area(poly, signed=false) = ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ) signed ? total : abs(total) : let( plane = plane_from_polygon(poly) ) - plane==[]? [] : + plane==[]? undef : let( n = plane_normal(plane), - total = sum([ - for(i=[1:1:len(poly)-2]) - let( - v1 = poly[i] - poly[0], - v2 = poly[i+1] - poly[0] - ) - cross(v1,v2) - ])* n/2 + total = + sum([ for(i=[1:1:len(poly)-2]) + cross(poly[i]-poly[0], poly[i+1]-poly[0]) + ]) * n/2 ) signed ? total : abs(total); -// Function: is_convex_polygon() -// Usage: -// is_convex_polygon(poly); -// Description: -// Returns true if the given 2D or 3D polygon is convex. -// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. -// If the points are collinear an error is generated. -// Arguments: -// poly = Polygon to check. -// eps = Tolerance for the collinearity test. Default: EPSILON. -// Example: -// is_convex_polygon(circle(d=50)); // Returns: true -// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true -// Example: -// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; -// is_convex_polygon(spiral); // Returns: false -function is_convex_polygon(poly,eps=EPSILON) = - assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly), - p0 = poly[0] ) - assert( lp>=3 , "A polygon must have at least 3 points" ) - let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(p0)==2 - ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) - min(crosses) >=0 || max(crosses)<=0 - : let( prod = crosses*sum(crosses), - minc = min(prod), - maxc = max(prod) ) - assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) - minc>=0 || maxc<=0; - - // Function: polygon_shift() // Usage: // polygon_shift(poly, i); @@ -1972,9 +1957,9 @@ function centroid(poly, eps=EPSILON) = // Returns -1 if the point is outside the polygon. // Returns 0 if the point is on the boundary. // Returns 1 if the point lies in the interior. -// The polygon does not need to be simple: it can have self-intersections. +// The polygon does not need to be simple: it may have self-intersections. // But the polygon cannot have holes (it must be simply connected). -// Rounding error may give mixed results for points on or near the boundary. +// Rounding errors may give mixed results for points on or near the boundary. // Arguments: // point = The 2D point to check position of. // poly = The list of 2D path points forming the perimeter of the polygon. @@ -2067,7 +2052,7 @@ function ccw_polygon(poly) = // poly = The list of the path points for the perimeter of the polygon. function reverse_polygon(poly) = assert(is_path(poly), "Input should be a polygon") - let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]]; + [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; // Function: polygon_normal() @@ -2075,7 +2060,7 @@ function reverse_polygon(poly) = // n = polygon_normal(poly); // Description: // Given a 3D planar polygon, returns a unit-length normal vector for the -// clockwise orientation of the polygon. If the polygon points are collinear, returns []. +// clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`. // It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. @@ -2083,7 +2068,7 @@ function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? [] : + triple==[] ? undef : point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; @@ -2237,4 +2222,255 @@ function split_polygons_at_each_z(polys, zs, _i=0) = ); + +// Section: Convex Sets + + +// Function: is_convex_polygon() +// Usage: +// is_convex_polygon(poly); +// Description: +// Returns true if the given 2D or 3D polygon is convex. +// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. +// If the points are collinear an error is generated. +// Arguments: +// poly = Polygon to check. +// eps = Tolerance for the collinearity test. Default: EPSILON. +// Example: +// is_convex_polygon(circle(d=50)); // Returns: true +// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true +// Example: +// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; +// is_convex_polygon(spiral); // Returns: false +function is_convex_polygon(poly,eps=EPSILON) = + assert(is_path(poly), "The input should be a 2D or 3D polygon." ) + let( lp = len(poly), + p0 = poly[0] ) + assert( lp>=3 , "A polygon must have at least 3 points" ) + let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) + len(p0)==2 + ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) + min(crosses) >=0 || max(crosses)<=0 + : let( prod = crosses*sum(crosses), + minc = min(prod), + maxc = max(prod) ) + assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) + minc>=0 || maxc<=0; + + +// Function: convex_distance() +// Usage: +// convex_distance(pts1, pts2,); +// See also: +// convex_collision +// Description: +// Returns the smallest distance between a point in convex hull of `points1` +// and a point in the convex hull of `points2`. All the points in the lists +// should have the same dimension, either 2D or 3D. +// A zero result means the hulls intercept whithin a tolerance `eps`. +// Arguments: +// points1 - first list of 2d or 3d points. +// points2 - second list of 2d or 3d points. +// eps - tolerance in distance evaluations. Default: EPSILON. +// Example(2D): +// pts1 = move([-3,0], p=square(3,center=true)); +// pts2 = rot(a=45, p=square(2,center=true)); +// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; +// polygon(pts1); +// polygon(pts2); +// polygon(pts3); +// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864 +// echo(convex_distance(pts2,pts3)); // Returns: 0 +// Example(3D): +// sphr1 = sphere(2,$fn=10); +// sphr2 = move([4,0,0], p=sphr1); +// sphr3 = move([4.5,0,0], p=sphr1); +// vnf_polyhedron(sphr1); +// vnf_polyhedron(sphr2); +// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0 +// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5 +function convex_distance(points1, points2, eps=EPSILON) = + assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), + "The input list should be a consistent non empty list of points of same dimension.") + assert(len(points1[0])==2 || len(points1[0])==3 , + "The input points should be 2d or 3d points.") + let( d = points1[0]-points2[0] ) + norm(d)); +// See also: +// convex_distance +// Description: +// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2` +// otherwise, `false`. +// All the points in the lists should have the same dimension, either 2D or 3D. +// This function is tipically faster than `convex_distance` to find a non-collision. +// Arguments: +// points1 - first list of 2d or 3d points. +// points2 - second list of 2d or 3d points. +// eps - tolerance for the intersection tests. Default: EPSILON. +// Example(2D): +// pts1 = move([-3,0], p=square(3,center=true)); +// pts2 = rot(a=45, p=square(2,center=true)); +// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; +// polygon(pts1); +// polygon(pts2); +// polygon(pts3); +// echo(convex_collision(pts1,pts2)); // Returns: false +// echo(convex_collision(pts2,pts3)); // Returns: true +// Example(3D): +// sphr1 = sphere(2,$fn=10); +// sphr2 = move([4,0,0], p=sphr1); +// sphr3 = move([4.5,0,0], p=sphr1); +// vnf_polyhedron(sphr1); +// vnf_polyhedron(sphr2); +// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true +// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false +// +function convex_collision(points1, points2, eps=EPSILON) = + assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), + "The input list should be a consistent non empty list of points of same dimension.") + assert(len(points1[0])==2 || len(points1[0])==3 , + "The input points should be 2d or 3d points.") + let( d = points1[0]-points2[0] ) + norm(d) eps ? false : // no collision + let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) + _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps); + + +// given a simplex s, returns a pair: +// - the point of the s closest to the origin +// - the smallest sub-simplex of s that contains that point +function _closest_simplex(s,eps=EPSILON) = + assert(len(s)>=2 && len(s)<=4, "Internal error.") + len(s)==2 ? _closest_s1(s,eps) : + len(s)==3 ? _closest_s2(s,eps) + : _closest_s3(s,eps); + + +// find the closest to a 1-simplex +// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf +function _closest_s1(s,eps=EPSILON) = + norm(s[1]-s[0])1 ? [ s[1], [s[1]] ] : + [ s[0]+t*c, s ]; + + +// find the closest to a 2-simplex +// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf +function _closest_s2(s,eps=EPSILON) = + let( + dim = len(s[0]), + a = dim==3 ? s[0]: [ each s[0], 0] , + b = dim==3 ? s[1]: [ each s[1], 0] , + c = dim==3 ? s[2]: [ each s[2], 0] , + ab = norm(a-b), + bc = norm(b-c), + ca = norm(c-a), + nr = cross(b-a,c-a) + ) + norm(nr) <= eps*max(ab,bc,ca) // degenerate case + ? let( i = max_index([ab, bc, ca]) ) + _closest_s1([s[i],s[(i+1)%3]],eps) +// considering that s[2] was the last inserted vertex in s, +// the only possible outcomes are : +// s, [s[0],s[2]] and [s[1],s[2]] + : let( + class = (cross(nr,a-b)*a<0 ? 1 : 0 ) + + (cross(nr,c-a)*a<0 ? 2 : 0 ) + + (cross(nr,b-c)*b<0 ? 4 : 0 ) + ) + assert( class!=1, "Internal error" ) + class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri +// class==1 ? _closest_s1([s[0],s[1]]) : + class==2 ? _closest_s1([s[0],s[2]],eps) : + class==4 ? _closest_s1([s[1],s[2]],eps) : +// class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) : + class==3 ? _closest_s1([s[0],s[2]],eps) : +// class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) : + class==5 ? _closest_s1([s[1],s[2]],eps) : + c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps); + + +// find the closest to a 3-simplex +// it seems that degenerate 3-simplices are correctly manage without extra code +function _closest_s3(s,eps=EPSILON) = + assert( len(s[0])==3 && len(s)==4, "Internal error." ) + let( nr = cross(s[1]-s[0],s[2]-s[0]), + sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] ) + norm(nr)0)==(nrm*s[i]<0) ) i ] + ) + len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex + len(facing)==1 ? _closest_s2(tris[facing[0]], eps) : + let( // look for the origin-facing tri closest to the origin + closest = [for(i=facing) _closest_s2(tris[i], eps) ], + dist = [for(cl=closest) norm(cl[0]) ], + nearest = min_index(dist) + ) + closest[nearest]; + + +function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]); + + +function _support_diff(p1,p2,d) = + let( p1d = p1*d, p2d = p2*d ) + p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]]; + + + // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 04f5126..31fc8ac 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -9,7 +9,9 @@ include <../std.scad> test_point_on_segment2d(); test_point_left_of_line2d(); test_collinear(); -test_distance_from_line(); +test_point_line_distance(); +test_point_segment_distance(); +test_segment_distance(); test_line_normal(); test_line_intersection(); //test_line_ray_intersection(); @@ -44,7 +46,7 @@ test_plane_normal(); test_plane_offset(); test_projection_on_plane(); test_plane_point_nearest_origin(); -test_distance_from_plane(); +test_point_plane_distance(); test__general_plane_line_intersection(); test_plane_line_angle(); @@ -88,7 +90,7 @@ test_cleanup_path(); test_simplify_path(); test_simplify_path_indexed(); test_is_region(); - +test_convex_distance(); // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -230,7 +232,7 @@ module test__general_plane_line_intersection() { interspoint = line1[0]+inters1[1]*(line1[1]-line1[0]); assert_approx(inters1[0],interspoint, info1); assert_approx(point3d(plane1)*inters1[0], plane1[3], info1); // interspoint on the plane - assert_approx(distance_from_plane(plane1, inters1[0]), 0, info1); // inters1[0] on the plane + assert_approx(point_plane_distance(plane1, inters1[0]), 0, info1); // inters1[0] on the plane } // line parallel to the plane @@ -351,13 +353,35 @@ module test_collinear() { *test_collinear(); -module test_distance_from_line() { - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,1,1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [-1,-1,-1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON); +module test_point_line_distance() { + assert_approx(point_line_distance([1,1,1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([-1,-1,-1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([1,-1,0], [[-10,-10,-10], [10,10,10]]), sqrt(2)); + assert_approx(point_line_distance([8,-8,0], [[-10,-10,-10], [10,10,10]]), 8*sqrt(2)); } -*test_distance_from_line(); +*test_point_line_distance(); + + +module test_point_segment_distance() { + assert_approx(point_segment_distance([3,8], [[-10,0], [10,0]]), 8); + assert_approx(point_segment_distance([14,3], [[-10,0], [10,0]]), 5); +} +*test_point_segment_distance(); + + +module test_segment_distance() { + assert_approx(segment_distance([[-14,3], [-14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-15,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,-3], [14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-14,3], [-14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-6,5], [4,-5]], [[-10,0], [10,0]]), 0); + assert_approx(segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]), 0); +} +*test_segment_distance(); module test_line_normal() { @@ -713,12 +737,12 @@ module test_plane_normal() { *test_plane_normal(); -module test_distance_from_plane() { +module test_point_plane_distance() { plane1 = plane3pt([-10,0,0], [0,10,0], [10,0,0]); - assert(distance_from_plane(plane1, [0,0,5]) == 5); - assert(distance_from_plane(plane1, [5,5,8]) == 8); + assert(point_plane_distance(plane1, [0,0,5]) == 5); + assert(point_plane_distance(plane1, [5,5,8]) == 8); } -*test_distance_from_plane(); +*test_point_plane_distance(); module test_polygon_line_intersection() { @@ -1051,6 +1075,20 @@ module test_is_region() { } *test_is_region(); - +module test_convex_distance() { + c1 = circle(10,$fn=24); + c2 = move([15,0], p=c1); + assert(convex_distance(c1, c2)==0); + c3 = move([22,0],c1); + assert(abs(convex_distance(c1, c3)-2) Date: Mon, 21 Jun 2021 19:51:26 +0100 Subject: [PATCH 04/12] Revert "Introduces convex collision and distance" This reverts commit eccb006f850b9becb816b0133aecc5a3c378bdf0. --- geometry.scad | 418 +++++++++------------------------------ tests/test_geometry.scad | 68 ++----- 2 files changed, 106 insertions(+), 380 deletions(-) diff --git a/geometry.scad b/geometry.scad index 0221a2b..24fbac6 100644 --- a/geometry.scad +++ b/geometry.scad @@ -19,8 +19,15 @@ // edge = Array of two points forming the line segment to test against. // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) function point_on_segment2d(point, edge, eps=EPSILON) = + assert( is_vector(point,2), "Invalid point." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - point_segment_distance(point, edge)= -eps*ne ) + && ( (dp-de)*de <= eps*ne ) // point projects on the segment + && _dist2line(point-edge[0],unit(de))eps*max(norm(line[1]),norm(line[0])); + && ! approx(norm(line[1]-line[0]), 0, eps); //Internal function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); @@ -78,57 +85,22 @@ function collinear(a, b, c, eps=EPSILON) = : noncollinear_triple(points,error=false,eps=eps)==[]; -// Function: point_line_distance() +// Function: distance_from_line() // Usage: -// point_line_distance(line, pt); +// distance_from_line(line, pt); // Description: // Finds the perpendicular distance of a point `pt` from the line `line`. // Arguments: // line = A list of two points, defining a line that both are on. // pt = A point to find the distance of from the line. // Example: -// dist = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 -function point_line_distance(pt, line) = +// distance_from_line([[-10,0], [10,0]], [3,8]); // Returns: 8 +function distance_from_line(line, pt) = assert( _valid_line(line) && is_vector(pt,len(line[0])), "Invalid line, invalid point or incompatible dimensions." ) _dist2line(pt-line[0],unit(line[1]-line[0])); -// Function: point_segment_distance() -// Usage: -// dist = point_segment_distance(pt, seg); -// Description: -// Returns the closest distance of the given point to the given line segment. -// Arguments: -// pt = The point to check the distance of. -// seg = The two points representing the line segment to check the distance of. -// Example: -// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 -// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5 -function point_segment_distance(pt, seg) = - assert( is_matrix(concat([pt],seg),3), - "Input should be a point and a valid segment with the dimension equal to the point." ) - norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) : - norm(pt-segment_closest_point(seg,pt)); - - -// Function: segment_distance() -// Usage: -// dist = segment_distance(seg1, seg2); -// Description: -// Returns the closest distance of the two given line segments. -// Arguments: -// seg1 = The list of two points representing the first line segment to check the distance of. -// seg2 = The list of two points representing the second line segment to check the distance of. -// Example: -// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5 -// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0 -function segment_distance(seg1, seg2) = - assert( is_matrix(concat(seg1,seg2),4), - "Inputs should be two valid segments." ) - convex_distance(seg1,seg2); - - // Function: line_normal() // Usage: // line_normal([P1,P2]) @@ -464,9 +436,17 @@ function ray_closest_point(ray,pt) = // color("blue") translate(pt) sphere(r=1,$fn=12); // color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = - assert( is_matrix(concat([pt],seg),3) , - "Invalid point or segment or incompatible dimensions." ) - pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0]; + assert(_valid_line(seg), "Invalid segment." ) + assert(len(pt)==len(seg[0]), "Incompatible dimensions." ) + approx(seg[0],seg[1])? seg[0] : + let( + seglen = norm(seg[1]-seg[0]), + segvec = (seg[1]-seg[0])/seglen, + projection = (pt-seg[0]) * segvec + ) + projection<=0 ? seg[0] : + projection>=seglen ? seg[1] : + seg[0] + projection*segvec; // Function: line_from_points() @@ -474,7 +454,7 @@ function segment_closest_point(seg,pt) = // line_from_points(points, [fast], [eps]); // Description: // Given a list of 2 or more collinear points, returns a line containing them. -// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned. +// If `fast` is false and the points are coincident, then `undef` is returned. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: // points = The list of points to find the line through. @@ -484,7 +464,7 @@ function line_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=undef), "Improper point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) let( pb = furthest_point(points[0],points) ) - norm(points[pb]-points[0])=0, "The tolerance should be a positive number." ) - assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.") + assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.") let( bounded = is_list(bounded)? bounded : [bounded, bounded], @@ -1202,7 +1182,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) assert(is_path(poly,dim=3), "Invalid polygon." ) assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).") - assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." ) + assert(_valid_line(line,dim=3,eps=eps), "Invalid line." ) let( bounded = is_list(bounded)? bounded : [bounded, bounded], poly = deduplicate(poly), @@ -1330,7 +1310,7 @@ function points_on_plane(points, plane, eps=EPSILON) = // plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`. // point = The 3D point to test. function in_front_of_plane(plane, point) = - point_plane_distance(plane, point) > EPSILON; + distance_from_plane(plane, point) > EPSILON; @@ -1656,7 +1636,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // eps = epsilon used for identifying the case with one solution. Default: 1e-9 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) = let(r=get_radius(r=r,d=d,dflt=undef)) - assert(_valid_line(line,2), "Invalid 2d line.") + assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_num(r) && r>0, "Radius must be positive") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition") @@ -1700,14 +1680,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = pb = points[b], nrm = norm(pa-pb) ) - nrm <= eps*max(norm(pa),norm(pb)) + approx(nrm, 0) ? assert(!error, "Cannot find three noncollinear points in pointlist.") [] : let( n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist) < eps*nrm + max(distlist)0 && len(pts[0])>0 , "Invalid pointlist." ) + let(ptsT = transpose(pts)) + [ + [for(row=ptsT) min(row)], + [for(row=ptsT) max(row)] + ]; // Function: closest_point() @@ -1768,7 +1747,7 @@ function furthest_point(pt, points) = // area = polygon_area(poly); // Description: // Given a 2D or 3D planar polygon, returns the area of that polygon. -// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`. +// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is []. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon. // Arguments: // poly = Polygon to compute the area of. @@ -1780,17 +1759,53 @@ function polygon_area(poly, signed=false) = ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ) signed ? total : abs(total) : let( plane = plane_from_polygon(poly) ) - plane==[]? undef : + plane==[]? [] : let( n = plane_normal(plane), - total = - sum([ for(i=[1:1:len(poly)-2]) - cross(poly[i]-poly[0], poly[i+1]-poly[0]) - ]) * n/2 + total = sum([ + for(i=[1:1:len(poly)-2]) + let( + v1 = poly[i] - poly[0], + v2 = poly[i+1] - poly[0] + ) + cross(v1,v2) + ])* n/2 ) signed ? total : abs(total); +// Function: is_convex_polygon() +// Usage: +// is_convex_polygon(poly); +// Description: +// Returns true if the given 2D or 3D polygon is convex. +// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. +// If the points are collinear an error is generated. +// Arguments: +// poly = Polygon to check. +// eps = Tolerance for the collinearity test. Default: EPSILON. +// Example: +// is_convex_polygon(circle(d=50)); // Returns: true +// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true +// Example: +// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; +// is_convex_polygon(spiral); // Returns: false +function is_convex_polygon(poly,eps=EPSILON) = + assert(is_path(poly), "The input should be a 2D or 3D polygon." ) + let( lp = len(poly), + p0 = poly[0] ) + assert( lp>=3 , "A polygon must have at least 3 points" ) + let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) + len(p0)==2 + ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) + min(crosses) >=0 || max(crosses)<=0 + : let( prod = crosses*sum(crosses), + minc = min(prod), + maxc = max(prod) ) + assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) + minc>=0 || maxc<=0; + + // Function: polygon_shift() // Usage: // polygon_shift(poly, i); @@ -1957,9 +1972,9 @@ function centroid(poly, eps=EPSILON) = // Returns -1 if the point is outside the polygon. // Returns 0 if the point is on the boundary. // Returns 1 if the point lies in the interior. -// The polygon does not need to be simple: it may have self-intersections. +// The polygon does not need to be simple: it can have self-intersections. // But the polygon cannot have holes (it must be simply connected). -// Rounding errors may give mixed results for points on or near the boundary. +// Rounding error may give mixed results for points on or near the boundary. // Arguments: // point = The 2D point to check position of. // poly = The list of 2D path points forming the perimeter of the polygon. @@ -2052,7 +2067,7 @@ function ccw_polygon(poly) = // poly = The list of the path points for the perimeter of the polygon. function reverse_polygon(poly) = assert(is_path(poly), "Input should be a polygon") - [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; + let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]]; // Function: polygon_normal() @@ -2060,7 +2075,7 @@ function reverse_polygon(poly) = // n = polygon_normal(poly); // Description: // Given a 3D planar polygon, returns a unit-length normal vector for the -// clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`. +// clockwise orientation of the polygon. If the polygon points are collinear, returns []. // It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. @@ -2068,7 +2083,7 @@ function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? undef : + triple==[] ? [] : point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; @@ -2222,255 +2237,4 @@ function split_polygons_at_each_z(polys, zs, _i=0) = ); - -// Section: Convex Sets - - -// Function: is_convex_polygon() -// Usage: -// is_convex_polygon(poly); -// Description: -// Returns true if the given 2D or 3D polygon is convex. -// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. -// If the points are collinear an error is generated. -// Arguments: -// poly = Polygon to check. -// eps = Tolerance for the collinearity test. Default: EPSILON. -// Example: -// is_convex_polygon(circle(d=50)); // Returns: true -// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true -// Example: -// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; -// is_convex_polygon(spiral); // Returns: false -function is_convex_polygon(poly,eps=EPSILON) = - assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly), - p0 = poly[0] ) - assert( lp>=3 , "A polygon must have at least 3 points" ) - let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(p0)==2 - ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) - min(crosses) >=0 || max(crosses)<=0 - : let( prod = crosses*sum(crosses), - minc = min(prod), - maxc = max(prod) ) - assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) - minc>=0 || maxc<=0; - - -// Function: convex_distance() -// Usage: -// convex_distance(pts1, pts2,); -// See also: -// convex_collision -// Description: -// Returns the smallest distance between a point in convex hull of `points1` -// and a point in the convex hull of `points2`. All the points in the lists -// should have the same dimension, either 2D or 3D. -// A zero result means the hulls intercept whithin a tolerance `eps`. -// Arguments: -// points1 - first list of 2d or 3d points. -// points2 - second list of 2d or 3d points. -// eps - tolerance in distance evaluations. Default: EPSILON. -// Example(2D): -// pts1 = move([-3,0], p=square(3,center=true)); -// pts2 = rot(a=45, p=square(2,center=true)); -// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; -// polygon(pts1); -// polygon(pts2); -// polygon(pts3); -// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864 -// echo(convex_distance(pts2,pts3)); // Returns: 0 -// Example(3D): -// sphr1 = sphere(2,$fn=10); -// sphr2 = move([4,0,0], p=sphr1); -// sphr3 = move([4.5,0,0], p=sphr1); -// vnf_polyhedron(sphr1); -// vnf_polyhedron(sphr2); -// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0 -// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5 -function convex_distance(points1, points2, eps=EPSILON) = - assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), - "The input list should be a consistent non empty list of points of same dimension.") - assert(len(points1[0])==2 || len(points1[0])==3 , - "The input points should be 2d or 3d points.") - let( d = points1[0]-points2[0] ) - norm(d)); -// See also: -// convex_distance -// Description: -// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2` -// otherwise, `false`. -// All the points in the lists should have the same dimension, either 2D or 3D. -// This function is tipically faster than `convex_distance` to find a non-collision. -// Arguments: -// points1 - first list of 2d or 3d points. -// points2 - second list of 2d or 3d points. -// eps - tolerance for the intersection tests. Default: EPSILON. -// Example(2D): -// pts1 = move([-3,0], p=square(3,center=true)); -// pts2 = rot(a=45, p=square(2,center=true)); -// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ]; -// polygon(pts1); -// polygon(pts2); -// polygon(pts3); -// echo(convex_collision(pts1,pts2)); // Returns: false -// echo(convex_collision(pts2,pts3)); // Returns: true -// Example(3D): -// sphr1 = sphere(2,$fn=10); -// sphr2 = move([4,0,0], p=sphr1); -// sphr3 = move([4.5,0,0], p=sphr1); -// vnf_polyhedron(sphr1); -// vnf_polyhedron(sphr2); -// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true -// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false -// -function convex_collision(points1, points2, eps=EPSILON) = - assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), - "The input list should be a consistent non empty list of points of same dimension.") - assert(len(points1[0])==2 || len(points1[0])==3 , - "The input points should be 2d or 3d points.") - let( d = points1[0]-points2[0] ) - norm(d) eps ? false : // no collision - let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) - _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps); - - -// given a simplex s, returns a pair: -// - the point of the s closest to the origin -// - the smallest sub-simplex of s that contains that point -function _closest_simplex(s,eps=EPSILON) = - assert(len(s)>=2 && len(s)<=4, "Internal error.") - len(s)==2 ? _closest_s1(s,eps) : - len(s)==3 ? _closest_s2(s,eps) - : _closest_s3(s,eps); - - -// find the closest to a 1-simplex -// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf -function _closest_s1(s,eps=EPSILON) = - norm(s[1]-s[0])1 ? [ s[1], [s[1]] ] : - [ s[0]+t*c, s ]; - - -// find the closest to a 2-simplex -// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf -function _closest_s2(s,eps=EPSILON) = - let( - dim = len(s[0]), - a = dim==3 ? s[0]: [ each s[0], 0] , - b = dim==3 ? s[1]: [ each s[1], 0] , - c = dim==3 ? s[2]: [ each s[2], 0] , - ab = norm(a-b), - bc = norm(b-c), - ca = norm(c-a), - nr = cross(b-a,c-a) - ) - norm(nr) <= eps*max(ab,bc,ca) // degenerate case - ? let( i = max_index([ab, bc, ca]) ) - _closest_s1([s[i],s[(i+1)%3]],eps) -// considering that s[2] was the last inserted vertex in s, -// the only possible outcomes are : -// s, [s[0],s[2]] and [s[1],s[2]] - : let( - class = (cross(nr,a-b)*a<0 ? 1 : 0 ) - + (cross(nr,c-a)*a<0 ? 2 : 0 ) - + (cross(nr,b-c)*b<0 ? 4 : 0 ) - ) - assert( class!=1, "Internal error" ) - class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri -// class==1 ? _closest_s1([s[0],s[1]]) : - class==2 ? _closest_s1([s[0],s[2]],eps) : - class==4 ? _closest_s1([s[1],s[2]],eps) : -// class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) : - class==3 ? _closest_s1([s[0],s[2]],eps) : -// class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) : - class==5 ? _closest_s1([s[1],s[2]],eps) : - c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps); - - -// find the closest to a 3-simplex -// it seems that degenerate 3-simplices are correctly manage without extra code -function _closest_s3(s,eps=EPSILON) = - assert( len(s[0])==3 && len(s)==4, "Internal error." ) - let( nr = cross(s[1]-s[0],s[2]-s[0]), - sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] ) - norm(nr)0)==(nrm*s[i]<0) ) i ] - ) - len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex - len(facing)==1 ? _closest_s2(tris[facing[0]], eps) : - let( // look for the origin-facing tri closest to the origin - closest = [for(i=facing) _closest_s2(tris[i], eps) ], - dist = [for(cl=closest) norm(cl[0]) ], - nearest = min_index(dist) - ) - closest[nearest]; - - -function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]); - - -function _support_diff(p1,p2,d) = - let( p1d = p1*d, p2d = p2*d ) - p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]]; - - - // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 31fc8ac..04f5126 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -9,9 +9,7 @@ include <../std.scad> test_point_on_segment2d(); test_point_left_of_line2d(); test_collinear(); -test_point_line_distance(); -test_point_segment_distance(); -test_segment_distance(); +test_distance_from_line(); test_line_normal(); test_line_intersection(); //test_line_ray_intersection(); @@ -46,7 +44,7 @@ test_plane_normal(); test_plane_offset(); test_projection_on_plane(); test_plane_point_nearest_origin(); -test_point_plane_distance(); +test_distance_from_plane(); test__general_plane_line_intersection(); test_plane_line_angle(); @@ -90,7 +88,7 @@ test_cleanup_path(); test_simplify_path(); test_simplify_path_indexed(); test_is_region(); -test_convex_distance(); + // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -232,7 +230,7 @@ module test__general_plane_line_intersection() { interspoint = line1[0]+inters1[1]*(line1[1]-line1[0]); assert_approx(inters1[0],interspoint, info1); assert_approx(point3d(plane1)*inters1[0], plane1[3], info1); // interspoint on the plane - assert_approx(point_plane_distance(plane1, inters1[0]), 0, info1); // inters1[0] on the plane + assert_approx(distance_from_plane(plane1, inters1[0]), 0, info1); // inters1[0] on the plane } // line parallel to the plane @@ -353,35 +351,13 @@ module test_collinear() { *test_collinear(); -module test_point_line_distance() { - assert_approx(point_line_distance([1,1,1], [[-10,-10,-10], [10,10,10]]), 0); - assert_approx(point_line_distance([-1,-1,-1], [[-10,-10,-10], [10,10,10]]), 0); - assert_approx(point_line_distance([1,-1,0], [[-10,-10,-10], [10,10,10]]), sqrt(2)); - assert_approx(point_line_distance([8,-8,0], [[-10,-10,-10], [10,10,10]]), 8*sqrt(2)); +module test_distance_from_line() { + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,1,1])) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [-1,-1,-1])) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON); + assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON); } -*test_point_line_distance(); - - -module test_point_segment_distance() { - assert_approx(point_segment_distance([3,8], [[-10,0], [10,0]]), 8); - assert_approx(point_segment_distance([14,3], [[-10,0], [10,0]]), 5); -} -*test_point_segment_distance(); - - -module test_segment_distance() { - assert_approx(segment_distance([[-14,3], [-14,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,3], [14,9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,-3], [-14,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[-14,-3], [-15,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,-3], [14,-9]], [[-10,0], [10,0]]), 5); - assert_approx(segment_distance([[14,3], [14,-3]], [[-10,0], [10,0]]), 4); - assert_approx(segment_distance([[-14,3], [-14,-3]], [[-10,0], [10,0]]), 4); - assert_approx(segment_distance([[-6,5], [4,-5]], [[-10,0], [10,0]]), 0); - assert_approx(segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]), 0); -} -*test_segment_distance(); +*test_distance_from_line(); module test_line_normal() { @@ -737,12 +713,12 @@ module test_plane_normal() { *test_plane_normal(); -module test_point_plane_distance() { +module test_distance_from_plane() { plane1 = plane3pt([-10,0,0], [0,10,0], [10,0,0]); - assert(point_plane_distance(plane1, [0,0,5]) == 5); - assert(point_plane_distance(plane1, [5,5,8]) == 8); + assert(distance_from_plane(plane1, [0,0,5]) == 5); + assert(distance_from_plane(plane1, [5,5,8]) == 8); } -*test_point_plane_distance(); +*test_distance_from_plane(); module test_polygon_line_intersection() { @@ -1075,20 +1051,6 @@ module test_is_region() { } *test_is_region(); -module test_convex_distance() { - c1 = circle(10,$fn=24); - c2 = move([15,0], p=c1); - assert(convex_distance(c1, c2)==0); - c3 = move([22,0],c1); - assert(abs(convex_distance(c1, c3)-2) Date: Mon, 21 Jun 2021 19:57:28 +0100 Subject: [PATCH 05/12] minor doc correction --- geometry.scad | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geometry.scad b/geometry.scad index e2569b9..46cff38 100644 --- a/geometry.scad +++ b/geometry.scad @@ -2257,7 +2257,7 @@ function is_convex_polygon(poly,eps=EPSILON) = // convex_distance(pts1, pts2,); // See also: // convex_collision -// Descrition: +// Description: // Returns the smallest distance between a point in convex hull of `points1` // and a point in the convex hull of `points2`. All the points in the lists // should have the same dimension, either 2D or 3D. @@ -2317,7 +2317,7 @@ function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) = // convex_collision(pts1, pts2,); // See also: // convex_distance -// Descrition: +// Description: // Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2` // otherwise, `false`. // All the points in the lists should have the same dimension, either 2D or 3D. From 9c2e4c23ac4b4197e0b048a2b8cb5a7da1efb8f9 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Mon, 21 Jun 2021 20:19:07 +0100 Subject: [PATCH 06/12] introduce convex collision and distance --- geometry.scad | 196 ++++++++++++++++++--------------------- tests/test_geometry.scad | 70 ++++++++++---- 2 files changed, 145 insertions(+), 121 deletions(-) diff --git a/geometry.scad b/geometry.scad index cf429c8..0221a2b 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1,4 +1,4 @@ -////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////// // LibFile: geometry.scad // Geometry helpers. // Includes: @@ -19,15 +19,8 @@ // edge = Array of two points forming the line segment to test against. // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) function point_on_segment2d(point, edge, eps=EPSILON) = - assert( is_vector(point,2), "Invalid point." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - assert( _valid_line(edge,2,eps=eps), "Invalid segment." ) - let( dp = point-edge[0], - de = edge[1]-edge[0], - ne = norm(de) ) - ( dp*de >= -eps*ne ) - && ( (dp-de)*de <= eps*ne ) // point projects on the segment - && _dist2line(point-edge[0],unit(de))eps*max(norm(line[1]),norm(line[0])); //Internal function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); @@ -85,22 +78,57 @@ function collinear(a, b, c, eps=EPSILON) = : noncollinear_triple(points,error=false,eps=eps)==[]; -// Function: distance_from_line() +// Function: point_line_distance() // Usage: -// distance_from_line(line, pt); +// point_line_distance(line, pt); // Description: // Finds the perpendicular distance of a point `pt` from the line `line`. // Arguments: // line = A list of two points, defining a line that both are on. // pt = A point to find the distance of from the line. // Example: -// distance_from_line([[-10,0], [10,0]], [3,8]); // Returns: 8 -function distance_from_line(line, pt) = +// dist = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +function point_line_distance(pt, line) = assert( _valid_line(line) && is_vector(pt,len(line[0])), "Invalid line, invalid point or incompatible dimensions." ) _dist2line(pt-line[0],unit(line[1]-line[0])); +// Function: point_segment_distance() +// Usage: +// dist = point_segment_distance(pt, seg); +// Description: +// Returns the closest distance of the given point to the given line segment. +// Arguments: +// pt = The point to check the distance of. +// seg = The two points representing the line segment to check the distance of. +// Example: +// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8 +// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5 +function point_segment_distance(pt, seg) = + assert( is_matrix(concat([pt],seg),3), + "Input should be a point and a valid segment with the dimension equal to the point." ) + norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) : + norm(pt-segment_closest_point(seg,pt)); + + +// Function: segment_distance() +// Usage: +// dist = segment_distance(seg1, seg2); +// Description: +// Returns the closest distance of the two given line segments. +// Arguments: +// seg1 = The list of two points representing the first line segment to check the distance of. +// seg2 = The list of two points representing the second line segment to check the distance of. +// Example: +// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5 +// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0 +function segment_distance(seg1, seg2) = + assert( is_matrix(concat(seg1,seg2),4), + "Inputs should be two valid segments." ) + convex_distance(seg1,seg2); + + // Function: line_normal() // Usage: // line_normal([P1,P2]) @@ -436,17 +464,9 @@ function ray_closest_point(ray,pt) = // color("blue") translate(pt) sphere(r=1,$fn=12); // color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = - assert(_valid_line(seg), "Invalid segment." ) - assert(len(pt)==len(seg[0]), "Incompatible dimensions." ) - approx(seg[0],seg[1])? seg[0] : - let( - seglen = norm(seg[1]-seg[0]), - segvec = (seg[1]-seg[0])/seglen, - projection = (pt-seg[0]) * segvec - ) - projection<=0 ? seg[0] : - projection>=seglen ? seg[1] : - seg[0] + projection*segvec; + assert( is_matrix(concat([pt],seg),3) , + "Invalid point or segment or incompatible dimensions." ) + pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0]; // Function: line_from_points() @@ -454,7 +474,7 @@ function segment_closest_point(seg,pt) = // line_from_points(points, [fast], [eps]); // Description: // Given a list of 2 or more collinear points, returns a line containing them. -// If `fast` is false and the points are coincident, then `undef` is returned. +// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: // points = The list of points to find the line through. @@ -464,7 +484,7 @@ function line_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=undef), "Improper point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) let( pb = furthest_point(points[0],points) ) - approx(norm(points[pb]-points[0]),0) ? undef : + norm(points[pb]-points[0])=0, "The tolerance should be a positive number." ) - assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.") + assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.") let( bounded = is_list(bounded)? bounded : [bounded, bounded], @@ -1182,7 +1202,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) assert(is_path(poly,dim=3), "Invalid polygon." ) assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).") - assert(_valid_line(line,dim=3,eps=eps), "Invalid line." ) + assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." ) let( bounded = is_list(bounded)? bounded : [bounded, bounded], poly = deduplicate(poly), @@ -1310,7 +1330,7 @@ function points_on_plane(points, plane, eps=EPSILON) = // plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`. // point = The 3D point to test. function in_front_of_plane(plane, point) = - distance_from_plane(plane, point) > EPSILON; + point_plane_distance(plane, point) > EPSILON; @@ -1636,7 +1656,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // eps = epsilon used for identifying the case with one solution. Default: 1e-9 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) = let(r=get_radius(r=r,d=d,dflt=undef)) - assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") + assert(_valid_line(line,2), "Invalid 2d line.") assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_num(r) && r>0, "Radius must be positive") assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition") @@ -1680,14 +1700,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = pb = points[b], nrm = norm(pa-pb) ) - approx(nrm, 0) + nrm <= eps*max(norm(pa),norm(pb)) ? assert(!error, "Cannot find three noncollinear points in pointlist.") [] : let( n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist)0 && len(pts[0])>0 , "Invalid pointlist." ) - let(ptsT = transpose(pts)) - [ - [for(row=ptsT) min(row)], - [for(row=ptsT) max(row)] - ]; + assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." ) + let( + select = ident(len(pts[0])), + spread = [for(i=[0:len(pts[0])-1]) + let( spreadi = pts*select[i] ) + [min(spreadi), max(spreadi)] ] ) + transpose(spread); // Function: closest_point() @@ -1747,7 +1768,7 @@ function furthest_point(pt, points) = // area = polygon_area(poly); // Description: // Given a 2D or 3D planar polygon, returns the area of that polygon. -// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is []. +// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon. // Arguments: // poly = Polygon to compute the area of. @@ -1759,53 +1780,17 @@ function polygon_area(poly, signed=false) = ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ) signed ? total : abs(total) : let( plane = plane_from_polygon(poly) ) - plane==[]? [] : + plane==[]? undef : let( n = plane_normal(plane), - total = sum([ - for(i=[1:1:len(poly)-2]) - let( - v1 = poly[i] - poly[0], - v2 = poly[i+1] - poly[0] - ) - cross(v1,v2) - ])* n/2 + total = + sum([ for(i=[1:1:len(poly)-2]) + cross(poly[i]-poly[0], poly[i+1]-poly[0]) + ]) * n/2 ) signed ? total : abs(total); -// Function: is_convex_polygon() -// Usage: -// is_convex_polygon(poly); -// Description: -// Returns true if the given 2D or 3D polygon is convex. -// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar. -// If the points are collinear an error is generated. -// Arguments: -// poly = Polygon to check. -// eps = Tolerance for the collinearity test. Default: EPSILON. -// Example: -// is_convex_polygon(circle(d=50)); // Returns: true -// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true -// Example: -// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; -// is_convex_polygon(spiral); // Returns: false -function is_convex_polygon(poly,eps=EPSILON) = - assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly), - p0 = poly[0] ) - assert( lp>=3 , "A polygon must have at least 3 points" ) - let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(p0)==2 - ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) - min(crosses) >=0 || max(crosses)<=0 - : let( prod = crosses*sum(crosses), - minc = min(prod), - maxc = max(prod) ) - assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) - minc>=0 || maxc<=0; - - // Function: polygon_shift() // Usage: // polygon_shift(poly, i); @@ -1972,9 +1957,9 @@ function centroid(poly, eps=EPSILON) = // Returns -1 if the point is outside the polygon. // Returns 0 if the point is on the boundary. // Returns 1 if the point lies in the interior. -// The polygon does not need to be simple: it can have self-intersections. +// The polygon does not need to be simple: it may have self-intersections. // But the polygon cannot have holes (it must be simply connected). -// Rounding error may give mixed results for points on or near the boundary. +// Rounding errors may give mixed results for points on or near the boundary. // Arguments: // point = The 2D point to check position of. // poly = The list of 2D path points forming the perimeter of the polygon. @@ -2067,7 +2052,7 @@ function ccw_polygon(poly) = // poly = The list of the path points for the perimeter of the polygon. function reverse_polygon(poly) = assert(is_path(poly), "Input should be a polygon") - let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]]; + [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; // Function: polygon_normal() @@ -2075,7 +2060,7 @@ function reverse_polygon(poly) = // n = polygon_normal(poly); // Description: // Given a 3D planar polygon, returns a unit-length normal vector for the -// clockwise orientation of the polygon. If the polygon points are collinear, returns []. +// clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`. // It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. @@ -2083,7 +2068,7 @@ function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? [] : + triple==[] ? undef : point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; @@ -2237,8 +2222,10 @@ function split_polygons_at_each_z(polys, zs, _i=0) = ); + // Section: Convex Sets + // Function: is_convex_polygon() // Usage: // is_convex_polygon(poly); @@ -2257,16 +2244,17 @@ function split_polygons_at_each_z(polys, zs, _i=0) = // is_convex_polygon(spiral); // Returns: false function is_convex_polygon(poly,eps=EPSILON) = assert(is_path(poly), "The input should be a 2D or 3D polygon." ) - let( lp = len(poly) ) + let( lp = len(poly), + p0 = poly[0] ) assert( lp>=3 , "A polygon must have at least 3 points" ) let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] ) - len(poly[0])==2 - ? assert( max(max(crosses),-min(crosses))>eps, "The points are collinear" ) + len(p0)==2 + ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" ) min(crosses) >=0 || max(crosses)<=0 : let( prod = crosses*sum(crosses), minc = min(prod), maxc = max(prod) ) - assert( max(maxc,-minc)>eps, "The points are collinear" ) + assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) minc>=0 || maxc<=0; @@ -2328,7 +2316,7 @@ function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) = close || [for(nv=norm(v), s=simplex) if(norm(s-v)<=eps*nv) 1]!=[] ? d : let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) _GJK_distance(points1, points2, eps, lbd, newsplx[0], newsplx[1]); - + // Function: convex_collision() // Usage: @@ -2371,7 +2359,7 @@ function convex_collision(points1, points2, eps=EPSILON) = norm(d)1 ? [ s[1], [s[1]] ] : [ s[0]+t*c, s ]; - - + + // find the closest to a 2-simplex // Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf function _closest_s2(s,eps=EPSILON) = @@ -2474,10 +2462,10 @@ function _closest_s3(s,eps=EPSILON) = nearest = min_index(dist) ) closest[nearest]; - + function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]); - + function _support_diff(p1,p2,d) = let( p1d = p1*d, p2d = p2*d ) @@ -2485,6 +2473,4 @@ function _support_diff(p1,p2,d) = - - // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 04f5126..a496b0f 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -9,7 +9,9 @@ include <../std.scad> test_point_on_segment2d(); test_point_left_of_line2d(); test_collinear(); -test_distance_from_line(); +test_point_line_distance(); +test_point_segment_distance(); +test_segment_distance(); test_line_normal(); test_line_intersection(); //test_line_ray_intersection(); @@ -44,7 +46,7 @@ test_plane_normal(); test_plane_offset(); test_projection_on_plane(); test_plane_point_nearest_origin(); -test_distance_from_plane(); +test_point_plane_distance(); test__general_plane_line_intersection(); test_plane_line_angle(); @@ -88,7 +90,7 @@ test_cleanup_path(); test_simplify_path(); test_simplify_path_indexed(); test_is_region(); - +test_convex_distance(); // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -230,7 +232,7 @@ module test__general_plane_line_intersection() { interspoint = line1[0]+inters1[1]*(line1[1]-line1[0]); assert_approx(inters1[0],interspoint, info1); assert_approx(point3d(plane1)*inters1[0], plane1[3], info1); // interspoint on the plane - assert_approx(distance_from_plane(plane1, inters1[0]), 0, info1); // inters1[0] on the plane + assert_approx(point_plane_distance(plane1, inters1[0]), 0, info1); // inters1[0] on the plane } // line parallel to the plane @@ -299,7 +301,7 @@ module test_line_from_points() { } *test_line_from_points(); -module test_point_on_segment2d() { +module test_point_on_segment2d() { assert(point_on_segment2d([-15,0], [[-10,0], [10,0]]) == false); assert(point_on_segment2d([-10,0], [[-10,0], [10,0]]) == true); assert(point_on_segment2d([-5,0], [[-10,0], [10,0]]) == true); @@ -351,13 +353,35 @@ module test_collinear() { *test_collinear(); -module test_distance_from_line() { - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,1,1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [-1,-1,-1])) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON); - assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON); +module test_point_line_distance() { + assert_approx(point_line_distance([1,1,1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([-1,-1,-1], [[-10,-10,-10], [10,10,10]]), 0); + assert_approx(point_line_distance([1,-1,0], [[-10,-10,-10], [10,10,10]]), sqrt(2)); + assert_approx(point_line_distance([8,-8,0], [[-10,-10,-10], [10,10,10]]), 8*sqrt(2)); } -*test_distance_from_line(); +*test_point_line_distance(); + + +module test_point_segment_distance() { + assert_approx(point_segment_distance([3,8], [[-10,0], [10,0]]), 8); + assert_approx(point_segment_distance([14,3], [[-10,0], [10,0]]), 5); +} +*test_point_segment_distance(); + + +module test_segment_distance() { + assert_approx(segment_distance([[-14,3], [-14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[-14,-3], [-15,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,-3], [14,-9]], [[-10,0], [10,0]]), 5); + assert_approx(segment_distance([[14,3], [14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-14,3], [-14,-3]], [[-10,0], [10,0]]), 4); + assert_approx(segment_distance([[-6,5], [4,-5]], [[-10,0], [10,0]]), 0); + assert_approx(segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]), 0); +} +*test_segment_distance(); module test_line_normal() { @@ -713,12 +737,12 @@ module test_plane_normal() { *test_plane_normal(); -module test_distance_from_plane() { +module test_point_plane_distance() { plane1 = plane3pt([-10,0,0], [0,10,0], [10,0,0]); - assert(distance_from_plane(plane1, [0,0,5]) == 5); - assert(distance_from_plane(plane1, [5,5,8]) == 8); + assert(point_plane_distance(plane1, [0,0,5]) == 5); + assert(point_plane_distance(plane1, [5,5,8]) == 8); } -*test_distance_from_plane(); +*test_point_plane_distance(); module test_polygon_line_intersection() { @@ -1051,6 +1075,20 @@ module test_is_region() { } *test_is_region(); - +module test_convex_distance() { + c1 = circle(10,$fn=24); + c2 = move([15,0], p=c1); + assert(convex_distance(c1, c2)==0); + c3 = move([22,0],c1); + assert(abs(convex_distance(c1, c3)-2) Date: Mon, 21 Jun 2021 20:43:12 +0100 Subject: [PATCH 07/12] test convex collision and distance --- tests/test_geometry.scad | 45 +++++++++++++++++++++++++++++++--------- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index a496b0f..992a681 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -91,6 +91,7 @@ test_simplify_path(); test_simplify_path_indexed(); test_is_region(); test_convex_distance(); +test_convex_collision(); // to be used when there are two alternative symmetrical outcomes // from a function like a plane output; v must be a vector @@ -1076,18 +1077,42 @@ module test_is_region() { *test_is_region(); module test_convex_distance() { +// 2D c1 = circle(10,$fn=24); - c2 = move([15,0], p=c1); - assert(convex_distance(c1, c2)==0); - c3 = move([22,0],c1); - assert(abs(convex_distance(c1, c3)-2) Date: Mon, 21 Jun 2021 23:24:54 +0100 Subject: [PATCH 08/12] introduce group_sort and enhance unique --- arrays.scad | 102 ++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 88 insertions(+), 14 deletions(-) diff --git a/arrays.scad b/arrays.scad index 40b232f..3972b0e 100644 --- a/arrays.scad +++ b/arrays.scad @@ -59,13 +59,16 @@ function _same_type(a,b, depth) = // Returns a portion of a list, wrapping around past the beginning, if end pivot) li ] + ) + concat( _group_sort_by_index(lesser,idx), + [equal], + _group_sort_by_index(greater,idx) ) ; + +function _group_sort(l) = + len(l) == 0 ? [] : + len(l) == 1 ? [l] : + let( pivot = l[floor(len(l)/2)], + equal = [ for(li=l) if( li==pivot) li ] , + lesser = [ for(li=l) if( li< pivot) li ] , + greater = [ for(li=l) if( li> pivot) li ] + ) + concat( _group_sort(lesser), + [equal], + _group_sort(greater) ) ; + // Sort a vector of scalar values with the native comparison operator // all elements should have the same type. @@ -1007,7 +1035,7 @@ function _sort_general(arr, idx=undef, indexed=false) = // lexical sort using compare_vals() function _lexical_sort(arr) = - arr==[] ? [] : len(arr)==1? arr : + len(arr)<=1? arr : let( pivot = arr[floor(len(arr)/2)] ) let( lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], @@ -1123,26 +1151,69 @@ function sortidx(list, idx=undef) = : _sort_general(list,idx,indexed=true); +// Function: group_sort() +// Usage: +// ulist = group_sort(list); +// Topics: List Handling +// See Also: shuffle(), sort(), sortidx(), unique(), unique_count() +// Description: +// Given a list of values, returns the sorted list with all repeated items grouped in a list. +// When the list entries are themselves lists, the sorting may be done based on the `idx` entry +// of those entries, that should be numbers. +// The result is always a list of lists. +// Arguments: +// list = The list to sort. +// idx = If given, do the comparison based just on the specified index. Default: zero. +// Example: +// sorted = group_sort([5,2,8,3,1,3,8,7,5]); // Returns: [[1],[2],[3,3],[5,5],[7],[8,8]] +// sorted2 = group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0); // Returns: [[[2,"b"],[2,"e"]], [[5,"a"],[5,"c"]], [[3,"d"]] ] +function group_sort(list, idx=undef) = + assert(is_list(list), "Input should be a list." ) + assert(is_undef(idx) || (is_finite(idx) && idx>=0) , "Invalid index." ) + len(list)<=1 ? [list] : + is_vector(list)? _group_sort(list) : + let( idx = is_undef(idx) ? 0 : idx ) + assert( [for(i=idx(list)) if(!is_list(list[i]) || len(list[i])>=idx || !is_num(list[idx])) 1]==[], + "Some entry of the list is a list shorter than `idx` or the indexed entry of it is not a number." ) + _group_sort_by_index(list,idx); + + // Function: unique() // Usage: // ulist = unique(list); // Topics: List Handling // See Also: shuffle(), sort(), sortidx(), unique_count() // Description: -// Returns a sorted list with all repeated items removed. +// Given a string or a list returns the sorted string or the sorted list with all repeated items removed. +// The sorting order of non homogeneous lists is the function `sort` order. // Arguments: // list = The list to uniquify. // Example: // sorted = unique([5,2,8,3,1,3,8,7,5]); // Returns: [1,2,3,5,7,8] +// sorted = unique("axdbxxc"); // Returns: "abcdx" +// sorted = unique([true,2,"xba",[1,0],true,[0,0],3,"a",[0,0],2]); // Returns: [true,2,3,"a","xba",[0,0],[1,0]] function unique(list) = assert(is_list(list)||is_string(list), "Invalid input." ) is_string(list)? str_join(unique([for (x = list) x])) : len(list)<=1? list : - let( sorted = sort(list)) - [ for (i=[0:1:len(sorted)-1]) - if (i==0 || (sorted[i] != sorted[i-1])) - sorted[i] - ]; + is_homogeneous(list,1) && ! is_list(list[0]) + ? _unique_sort(list) + : let( sorted = sort(list)) + [ for (i=[0:1:len(sorted)-1]) + if (i==0 || (sorted[i] != sorted[i-1])) + sorted[i] + ]; + +function _unique_sort(l) = + len(l) <= 1 ? l : + let( pivot = l[floor(len(l)/2)], + equal = [ for(li=l) if( li==pivot) li ] , + lesser = [ for(li=l) if( lipivot) li ] + ) + concat( _unique_sort(lesser), + equal[0], + _unique_sort(greater) ) ; // Function: unique_count() @@ -1160,11 +1231,14 @@ function unique(list) = function unique_count(list) = assert(is_list(list) || is_string(list), "Invalid input." ) list == [] ? [[],[]] : - let( list=sort(list) ) - let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] ) - [ select(list,ind), deltas( concat(ind,[len(list)]) ) ]; - + is_homogeneous(list,1) && ! is_list(list[0]) + ? let( sorted = _group_sort_scalars(list) ) + [ [for(s=sorted) s[0] ], [for(s=sorted) len(s) ] ] + : let( list=sort(list) ) + let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] ) + [ select(list,ind), deltas( concat(ind,[len(list)]) ) ]; + // Section: List Iteration Helpers // Function: idx() From b924eaa303000dc5465978cb288d2cc772c7fa88 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Tue, 22 Jun 2021 01:13:40 +0100 Subject: [PATCH 09/12] test group_sort --- arrays.scad | 18 +++++++++--------- tests/test_arrays.scad | 16 +++++++++++++--- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/arrays.scad b/arrays.scad index 3972b0e..0df54af 100644 --- a/arrays.scad +++ b/arrays.scad @@ -68,8 +68,8 @@ function _same_type(a,b, depth) = // list = select(list,start,end); // Arguments: // list = The list to get the portion of. -// start = The index of the first item or an index range or a list of indices. -// end = The index of the last item. +// start = Either the index of the first item or an index range or a list of indices. +// end = The index of the last item when `start` is a number. When `start` is a list or a range, `end` should not be given. // See Also: slice(), subindex(), last() // Example: // l = [3,4,5,6,7,8,9]; @@ -92,7 +92,7 @@ function select(list, start, end) = ? list[ (start%l+l)%l ] : assert( is_list(start) || is_range(start), "Invalid start parameter") [for (i=start) list[ (i%l+l)%l ] ] - : assert(is_finite(start), "Invalid start parameter.") + : assert(is_finite(start), "When `end` is given, `start` parameter should be a number.") assert(is_finite(end), "Invalid end parameter.") let( s = (start%l+l)%l, e = (end%l+l)%l ) (s <= e) @@ -1062,7 +1062,7 @@ function _indexed_sort(arrind) = // Usage: // slist = sort(list, ); // Topics: List Handling -// See Also: shuffle(), sortidx(), unique(), unique_count() +// See Also: shuffle(), sortidx(), unique(), unique_count(), group_sort() // Description: // Sorts the given list in lexicographic order. If the input is a homogeneous simple list or a homogeneous // list of vectors (see function is_homogeneous), the sorting method uses the native comparison operator and is faster. @@ -1103,7 +1103,7 @@ function sort(list, idx=undef) = // Usage: // idxlist = sortidx(list, ); // Topics: List Handling -// See Also: shuffle(), sort(), unique(), unique_count() +// See Also: shuffle(), sort(), group_sort(), unique(), unique_count() // Description: // Given a list, sort it as function `sort()`, and returns // a list of indexes into the original list in that sorted order. @@ -1167,13 +1167,13 @@ function sortidx(list, idx=undef) = // Example: // sorted = group_sort([5,2,8,3,1,3,8,7,5]); // Returns: [[1],[2],[3,3],[5,5],[7],[8,8]] // sorted2 = group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0); // Returns: [[[2,"b"],[2,"e"]], [[5,"a"],[5,"c"]], [[3,"d"]] ] -function group_sort(list, idx=undef) = +function group_sort(list, idx) = assert(is_list(list), "Input should be a list." ) assert(is_undef(idx) || (is_finite(idx) && idx>=0) , "Invalid index." ) len(list)<=1 ? [list] : is_vector(list)? _group_sort(list) : let( idx = is_undef(idx) ? 0 : idx ) - assert( [for(i=idx(list)) if(!is_list(list[i]) || len(list[i])>=idx || !is_num(list[idx])) 1]==[], + assert( [for(entry=list) if(!is_list(entry) || len(entry) Date: Tue, 22 Jun 2021 01:20:05 +0100 Subject: [PATCH 10/12] Extend convolve and deltas; do minor doc corrections --- coords.scad | 2 +- distributors.scad | 2 +- math.scad | 38 +++++++++++++++++++++++++------------- paths.scad | 3 ++- 4 files changed, 29 insertions(+), 16 deletions(-) diff --git a/coords.scad b/coords.scad index b32ad73..f0b4272 100644 --- a/coords.scad +++ b/coords.scad @@ -237,7 +237,7 @@ function project_plane(plane,p) = let(plane = plane_from_points(plane)) assert(is_def(plane), "Point list is not coplanar") project_plane(plane) - : assert(is_def(p), str("Invalid plane specification",plane)) + : assert(is_def(p), str("Invalid plane specification: ",plane)) is_vnf(p) ? [project_plane(plane,p[0]), p[1]] : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ? // bezier patch or region [for(plist=p) project_plane(plane,plist)] diff --git a/distributors.scad b/distributors.scad index 6010b12..ee60d6b 100644 --- a/distributors.scad +++ b/distributors.scad @@ -34,7 +34,7 @@ module move_copies(a=[[0,0,0]]) assert(is_list(a)); for ($idx = idx(a)) { $pos = a[$idx]; - assert(is_vector($pos)); + assert(is_vector($pos),"move_copies offsets should be a 2d or 3d vector."); translate($pos) children(); } } diff --git a/math.scad b/math.scad index 96bd475..30b585e 100644 --- a/math.scad +++ b/math.scad @@ -640,17 +640,19 @@ function sum_of_sines(a, sines) = // Usage: // delts = deltas(v); // Description: -// Returns a list with the deltas of adjacent entries in the given list. +// Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front. // The list should be a consistent list of numeric components (numbers, vectors, matrix, etc). // Given [a,b,c,d], returns [b-a,c-b,d-c]. +// // Arguments: // v = The list to get the deltas of. +// wrap = If true, wrap back to the start from the end. ie: return the difference between the last and first items as the last delta. Default: false // Example: // deltas([2,5,9,17]); // returns [3,4,8]. // deltas([[1,2,3], [3,6,8], [4,8,11]]); // returns [[2,4,5], [1,2,3]] -function deltas(v) = +function deltas(v, wrap=false) = assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.") - [for (p=pair(v)) p[1]-p[0]] ; + [for (p=pair(v,wrap)) p[1]-p[0]] ; // Function: product() @@ -771,21 +773,31 @@ function _med3(a,b,c) = // Usage: // x = convolve(p,q); // Description: -// Given two vectors, finds the convolution of them. -// The length of the returned vector is len(p)+len(q)-1 . +// Given two vectors, or one vector and a path or +// two paths of the same dimension, finds the convolution of them. +// If both parameter are vectors, returns the vector convolution. +// If one parameter is a vector and the other a path, +// convolves using products by scalars and returns a path. +// If both parameters are paths, convolve using scalar products +// and returns a vector. +// The returned vector or path has length len(p)+len(q)-1. // Arguments: -// p = The first vector. -// q = The second vector. +// p = The first vector or path. +// q = The second vector or path. // Example: // a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1] // b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3] +// c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]] +// d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns: [3,9,11,7] function convolve(p,q) = p==[] || q==[] ? [] : - assert( is_vector(p) && is_vector(q), "The inputs should be vectors.") + assert( (is_vector(p) || is_matrix(p)) + && ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) , + "The inputs should be vectors or paths all of the same dimension.") let( n = len(p), m = len(q)) [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) ) - [for(j=[k1:k2]) p[i-j] ] * [for(j=[k1:k2]) q[j] ] + sum([for(j=[k1:k2]) p[i-j]*q[j] ]) ]; @@ -1694,7 +1706,7 @@ function polynomial(p,z,k,total) = // x = polymult(p,q) // x = polymult([p1,p2,p3,...]) // Description: -// Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, +// Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first, // computes the coefficient list of the product polynomial. function poly_mult(p,q) = is_undef(q) ? @@ -1714,8 +1726,8 @@ function poly_mult(p,q) = // Description: // Computes division of the numerator polynomial by the denominator polynomial and returns // a list of two polynomials, [quotient, remainder]. If the division has no remainder then -// the zero polynomial [] is returned for the remainder. Similarly if the quotient is zero -// the returned quotient will be []. +// the zero polynomial [0] is returned for the remainder. Similarly if the quotient is zero +// the returned quotient will be [0]. function poly_div(n,d) = assert( is_vector(n) && is_vector(d) , "Invalid polynomials." ) let( d = _poly_trim(d), @@ -1740,7 +1752,7 @@ function _poly_div(n,d,q) = /// _poly_trim(p,[eps]) /// Description: /// Removes leading zero terms of a polynomial. By default zeros must be exact, -/// or give epsilon for approximate zeros. +/// or give epsilon for approximate zeros. Returns [0] for a zero polynomial. function _poly_trim(p,eps=0) = let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i]) len(nz)==0 ? [0] : list_tail(p,nz[0]); diff --git a/paths.scad b/paths.scad index 88d9ec7..3c79b97 100644 --- a/paths.scad +++ b/paths.scad @@ -997,7 +997,7 @@ module jittered_poly(path, dist=1/512) { // Module: extrude_from_to() // Description: -// Extrudes a 2D shape between the points pt1 and pt2. Takes as children a set of 2D shapes to extrude. +// Extrudes a 2D shape between the 3d points pt1 and pt2. Takes as children a set of 2D shapes to extrude. // Arguments: // pt1 = starting point of extrusion. // pt2 = ending point of extrusion. @@ -1010,6 +1010,7 @@ module jittered_poly(path, dist=1/512) { // xcopies(3) circle(3, $fn=32); // } module extrude_from_to(pt1, pt2, convexity, twist, scale, slices) { + assert( is_path([pt1,pt2],3), "The points should be 3d points"); rtp = xyz_to_spherical(pt2-pt1); translate(pt1) { rotate([0, rtp[2], rtp[1]]) { From 0b38bc463fb445bd65a4a11ace5b2ddf120f3c7a Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Tue, 22 Jun 2021 01:34:39 +0100 Subject: [PATCH 11/12] update test_math --- tests/test_math.scad | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_math.scad b/tests/test_math.scad index adcfb6d..b5c25b8 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -546,7 +546,9 @@ test_sum_of_sines(); module test_deltas() { assert_equal(deltas([2,5,9,17]), [3,4,8]); + assert_equal(deltas([2,5,9,17],wrap=true), [3,4,8,-15]); assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]]), [[2,4,5], [1,2,3]]); + assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]],wrap=true), [[2,4,5], [1,2,3], [-3,-6,-8]]); } test_deltas(); @@ -582,6 +584,10 @@ module test_convolve() { assert_equal(convolve([1,1],[]), []); assert_equal(convolve([1,1],[1,2,1]), [1,3,3,1]); assert_equal(convolve([1,2,3],[1,2,1]), [1,4,8,8,3]); + assert_equal(convolve([1,2,3],[[1],[2],[1]]), [[1], [4], [8], [8], [3]]); + assert_equal(convolve([[1],[2],[3]],[[1],[2],[1]]), [1,4,8,8,3]); + assert_equal(convolve([[1,0],[2,1],[3,2]],[[1,0],[2,1],[1,2]]), [1,4,9,12,7]); + assert_equal(convolve([1,2,3],[[1,0],[2,1],[1,2]]), [[1,0],[4,1],[8,4],[8,7],[3,6]]); } test_convolve(); From cb5cd93f2530696642751893f5d1ac0700fbd917 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Tue, 22 Jun 2021 04:12:14 +0100 Subject: [PATCH 12/12] cleaning --- hull.scad | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/hull.scad b/hull.scad index 1dc147c..4adb7c3 100644 --- a/hull.scad +++ b/hull.scad @@ -202,17 +202,15 @@ function hull3d_faces(points) = // Adds the remaining points one by one to the convex hull -function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = +function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = //let( EPSILON=1e-12 ) _i >= len(remaining) ? triangles : let ( // pick a point i = remaining[_i], // evaluate the triangle plane equations at point i -// xx=[for(i=[0:len(planes)-1],p=[planes[i]]) if(len(p)!=4) echo(i=i,len(p))0 ],//echo([each points[i], -1]), -// planeq_val = [for(p=planes) p*[each points[i], -1]], planeq_val = planes*[each points[i], -1], // find the triangles that are in conflict with the point (point not inside) - conflicts = [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]>EPSILON) i ], + conflicts = [for (i = [0:1:len(planeq_val)-1]) if (planeq_val[i]>EPSILON) i ], // collect the halfedges of all triangles that are in conflict halfedges = [ for(c = conflicts, i = [0:2]) @@ -222,18 +220,15 @@ function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = horizon = _remove_internal_edges(halfedges), // generate new triangles connecting point i to each horizon halfedge vertices tri2add = [ for (h = horizon) concat(h,i) ], -// w=[for(t=tri2add) if(collinear(points[t[0]],points[t[1]],points[t[2]])) echo(t)], // add tria2add and remove conflict triangles new_triangles = concat( tri2add, [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) triangles[i] ] ), -// y=[for(t=tri2add) if([]==plane3pt_indexed(points,t[0],t[1],t[2])) echo(tri2add=t,pts=[for(ti=t) points[ti]])], // add the plane equations of new added triangles and remove the plane equations of the conflict ones new_planes = - concat( [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ], - [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ] - ) + [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) , + for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ] ) _hull3d_iterative( points, new_triangles, @@ -249,7 +244,6 @@ function _remove_internal_edges(halfedges) = [ h ]; - function _find_first_noncoplanar(plane, points, i=0) = (i >= len(points) || !points_on_plane([points[i]],plane))? i : _find_first_noncoplanar(plane, points, i+1);