From 49a3a166eb153105c8c667b4b356e3904f88d237 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sat, 10 Apr 2021 20:09:03 +0100 Subject: [PATCH 01/12] Review of geometry.scad for speed --- arrays.scad | 16 +++--- common.scad | 11 ++-- geometry.scad | 112 ++++++++++++++++++++++++++------------- tests/test_geometry.scad | 8 +-- 4 files changed, 95 insertions(+), 52 deletions(-) diff --git a/arrays.scad b/arrays.scad index ac97444..604186b 100644 --- a/arrays.scad +++ b/arrays.scad @@ -41,6 +41,7 @@ function is_homogeneous(l, depth=10) = [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) ) 0 ]; function is_homogenous(l, depth=10) = is_homogeneous(l, depth); + function _same_type(a,b, depth) = (depth==0) || @@ -50,7 +51,7 @@ function _same_type(a,b, depth) = (is_string(a) && is_string(b)) || (is_list(a) && is_list(b) && len(a)==len(b) && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); - + // Function: select() // Topics: List Handling @@ -97,7 +98,6 @@ function select(list, start, end) = // Function: slice() -// Topics: List Handling // Usage: // list = slice(list,s,e); // Description: @@ -476,7 +476,7 @@ function reverse(x) = // l9 = list_rotate([1,2,3,4,5],6); // Returns: [2,3,4,5,1] function list_rotate(list,n=1) = assert(is_list(list)||is_string(list), "Invalid list or string.") - assert(is_finite(n), "Invalid number") + assert(is_int(n), "The rotation number should be integer") let ( ll = len(list), n = ((n % ll) + ll) % ll, @@ -990,7 +990,7 @@ function _sort_vectors(arr, idxlist, _i=0) = _sort_vectors(equal, idxlist, _i+1), _sort_vectors(greater, idxlist, _i ) ); - + // sorting using compare_vals(); returns indexed list when `indexed==true` function _sort_general(arr, idx=undef, indexed=false) = (len(arr)<=1) ? arr : @@ -1332,6 +1332,8 @@ function permutations(l,n=2) = // pairs = zip(a,b); // triples = zip(a,b,c); // quads = zip([LIST1,LIST2,LIST3,LIST4]); +// Topics: List Handling, Iteration +// See Also: zip_long() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1357,6 +1359,8 @@ function zip(a,b,c) = // pairs = zip_long(a,b); // triples = zip_long(a,b,c); // quads = zip_long([LIST1,LIST2,LIST3,LIST4]); +// Topics: List Handling, Iteration +// See Also: zip() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1526,7 +1530,6 @@ function subindex(M, idx) = // [[4,2], 91, false], // [6, [3,4], undef]]; // submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]] - function submatrix(M,idx1,idx2) = [for(i=idx1) [for(j=idx2) M[i][j] ] ]; @@ -1629,7 +1632,6 @@ function block_matrix(M) = assert(badrows==[], "Inconsistent or invalid input") bigM; - // Function: diagonal_matrix() // Usage: // mat = diagonal_matrix(diag, ); @@ -1855,7 +1857,7 @@ function transpose(arr, reverse=false) = // A = matrix to test // eps = epsilon for comparing equality. Default: 1e-12 function is_matrix_symmetric(A,eps=1e-12) = - approx(A,transpose(A)); + approx(A,transpose(A), eps); // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/common.scad b/common.scad index e24c744..6d2d44f 100644 --- a/common.scad +++ b/common.scad @@ -205,7 +205,8 @@ function is_func(x) = version_num()>20210000 && is_function(x); // Description: // Tests whether input is a list of entries which all have the same list structure // and are filled with finite numerical data. You can optionally specify a required -// list structure with the pattern argument. It returns `true` for the empty list. +// list structure with the pattern argument. +// It returns `true` for the empty list regardless the value of the `pattern`. // Arguments: // list = list to check // pattern = optional pattern required to match @@ -293,7 +294,7 @@ function default(v,dflt=undef) = is_undef(v)? dflt : v; // v = The list whose items are being checked. // recursive = If true, sublists are checked recursively for defined values. The first sublist that has a defined item is returned. // Examples: -// val = first_defined([undef,7,undef,true]); // Returns: 1 +// val = first_defined([undef,7,undef,true]); // Returns: 7 function first_defined(v,recursive=false,_i=0) = _i, ); // Description: // Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane, // that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1. -// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned. -// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. +// If the points in the list are collinear or not coplanar, then `undef` is returned. // Arguments: // points = The list of points to find the plane of. -// fast = If true, don't verify that all points in the list are coplanar. Default: false // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) // Example(3D): // xyzpath = rot(45, v=[-0.3,1,0], p=path3d(star(n=6,id=70,d=100), 70)); @@ -911,20 +947,21 @@ function plane_from_normal(normal, pt=[0,0,0]) = // #stroke(xyzpath,closed=true); // cp = centroid(xyzpath); // move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(); -function plane_from_points(points, fast=false, eps=EPSILON) = +function plane_from_points(points,fast=false, eps=EPSILON) = assert( is_path(points,dim=3), "Improper 3d point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - let( - indices = noncollinear_triple(points,error=false) - ) - indices==[] ? undef : - let( - p1 = points[indices[0]], - p2 = points[indices[1]], - p3 = points[indices[2]], - plane = plane3pt(p1,p2,p3) - ) - fast || points_on_plane(points,plane,eps=eps) ? plane : undef; + len(points) == 3 + ? let( plane = plane3pt(points[0],points[1],points[2]) ) + plane==[] ? undef : plane + : let( + cov_evals = _covariance_evals(points), + pm = cov_evals[0], + evals = cov_evals[1], + M = cov_evals[2], + evec = _eigenvec_symm_3(M,evals,i=2) ) +// echo(error_points_plane= abs(max(points*evec)-pm*evec), limit=eps) + !fast && abs(max(points*evec)-pm*evec)>eps*evals[0] ? undef : + [ each evec, pm*evec] ; // Function: plane_from_polygon() @@ -945,17 +982,16 @@ function plane_from_points(points, fast=false, eps=EPSILON) = // #stroke(xyzpath,closed=true); // cp = centroid(xyzpath); // move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(); -function plane_from_polygon(poly, fast=false, eps=EPSILON) = +function plane_from_polygon(poly, fast=false, eps=EPSILON) = assert( is_path(poly,dim=3), "Invalid polygon." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - let( - poly = deduplicate(poly), - n = polygon_normal(poly), - plane = [n.x, n.y, n.z, n*poly[0]] - ) - fast? plane: coplanar(poly,eps=eps)? plane: []; - - + len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) : + let( triple = sort(noncollinear_triple(poly,error=false)) ) + triple==[] ? [] : + let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) + fast? plane: points_on_plane(poly, plane, eps=eps)? plane: []; + + // Function: plane_normal() // Usage: // plane_normal(plane); @@ -1252,9 +1288,11 @@ function coplanar(points, eps=EPSILON) = len(points)<=2 ? false : let( ip = noncollinear_triple(points,error=false,eps=eps) ) ip == [] ? false : - let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), - normal = point3d(plane) ) - max( points*normal ) - plane[3]< eps*norm(normal); + let( + plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), + normal = point3d(plane), + pt_nrm = points*normal ) + abs(max(max(pt_nrm)-plane[3], -min(pt_nrm)+plane[3])) < eps; // Function: points_on_plane() @@ -1665,11 +1703,11 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist) Date: Sat, 10 Apr 2021 20:14:40 +0100 Subject: [PATCH 02/12] Revert "Review of geometry.scad for speed" This reverts commit 49a3a166eb153105c8c667b4b356e3904f88d237. --- arrays.scad | 16 +++--- common.scad | 11 ++-- geometry.scad | 112 +++++++++++++-------------------------- tests/test_geometry.scad | 8 ++- 4 files changed, 52 insertions(+), 95 deletions(-) diff --git a/arrays.scad b/arrays.scad index 604186b..ac97444 100644 --- a/arrays.scad +++ b/arrays.scad @@ -41,7 +41,6 @@ function is_homogeneous(l, depth=10) = [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) ) 0 ]; function is_homogenous(l, depth=10) = is_homogeneous(l, depth); - function _same_type(a,b, depth) = (depth==0) || @@ -51,7 +50,7 @@ function _same_type(a,b, depth) = (is_string(a) && is_string(b)) || (is_list(a) && is_list(b) && len(a)==len(b) && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); - + // Function: select() // Topics: List Handling @@ -98,6 +97,7 @@ function select(list, start, end) = // Function: slice() +// Topics: List Handling // Usage: // list = slice(list,s,e); // Description: @@ -476,7 +476,7 @@ function reverse(x) = // l9 = list_rotate([1,2,3,4,5],6); // Returns: [2,3,4,5,1] function list_rotate(list,n=1) = assert(is_list(list)||is_string(list), "Invalid list or string.") - assert(is_int(n), "The rotation number should be integer") + assert(is_finite(n), "Invalid number") let ( ll = len(list), n = ((n % ll) + ll) % ll, @@ -990,7 +990,7 @@ function _sort_vectors(arr, idxlist, _i=0) = _sort_vectors(equal, idxlist, _i+1), _sort_vectors(greater, idxlist, _i ) ); - + // sorting using compare_vals(); returns indexed list when `indexed==true` function _sort_general(arr, idx=undef, indexed=false) = (len(arr)<=1) ? arr : @@ -1332,8 +1332,6 @@ function permutations(l,n=2) = // pairs = zip(a,b); // triples = zip(a,b,c); // quads = zip([LIST1,LIST2,LIST3,LIST4]); -// Topics: List Handling, Iteration -// See Also: zip_long() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1359,8 +1357,6 @@ function zip(a,b,c) = // pairs = zip_long(a,b); // triples = zip_long(a,b,c); // quads = zip_long([LIST1,LIST2,LIST3,LIST4]); -// Topics: List Handling, Iteration -// See Also: zip() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1530,6 +1526,7 @@ function subindex(M, idx) = // [[4,2], 91, false], // [6, [3,4], undef]]; // submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]] + function submatrix(M,idx1,idx2) = [for(i=idx1) [for(j=idx2) M[i][j] ] ]; @@ -1632,6 +1629,7 @@ function block_matrix(M) = assert(badrows==[], "Inconsistent or invalid input") bigM; + // Function: diagonal_matrix() // Usage: // mat = diagonal_matrix(diag, ); @@ -1857,7 +1855,7 @@ function transpose(arr, reverse=false) = // A = matrix to test // eps = epsilon for comparing equality. Default: 1e-12 function is_matrix_symmetric(A,eps=1e-12) = - approx(A,transpose(A), eps); + approx(A,transpose(A)); // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/common.scad b/common.scad index 6d2d44f..e24c744 100644 --- a/common.scad +++ b/common.scad @@ -205,8 +205,7 @@ function is_func(x) = version_num()>20210000 && is_function(x); // Description: // Tests whether input is a list of entries which all have the same list structure // and are filled with finite numerical data. You can optionally specify a required -// list structure with the pattern argument. -// It returns `true` for the empty list regardless the value of the `pattern`. +// list structure with the pattern argument. It returns `true` for the empty list. // Arguments: // list = list to check // pattern = optional pattern required to match @@ -294,7 +293,7 @@ function default(v,dflt=undef) = is_undef(v)? dflt : v; // v = The list whose items are being checked. // recursive = If true, sublists are checked recursively for defined values. The first sublist that has a defined item is returned. // Examples: -// val = first_defined([undef,7,undef,true]); // Returns: 7 +// val = first_defined([undef,7,undef,true]); // Returns: 1 function first_defined(v,recursive=false,_i=0) = _i, ); // Description: // Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane, // that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1. -// If the points in the list are collinear or not coplanar, then `undef` is returned. +// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned. +// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. // Arguments: // points = The list of points to find the plane of. +// fast = If true, don't verify that all points in the list are coplanar. Default: false // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9) // Example(3D): // xyzpath = rot(45, v=[-0.3,1,0], p=path3d(star(n=6,id=70,d=100), 70)); @@ -947,21 +911,20 @@ function _covariance_evals(points) = // #stroke(xyzpath,closed=true); // cp = centroid(xyzpath); // move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(); -function plane_from_points(points,fast=false, eps=EPSILON) = +function plane_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=3), "Improper 3d point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - len(points) == 3 - ? let( plane = plane3pt(points[0],points[1],points[2]) ) - plane==[] ? undef : plane - : let( - cov_evals = _covariance_evals(points), - pm = cov_evals[0], - evals = cov_evals[1], - M = cov_evals[2], - evec = _eigenvec_symm_3(M,evals,i=2) ) -// echo(error_points_plane= abs(max(points*evec)-pm*evec), limit=eps) - !fast && abs(max(points*evec)-pm*evec)>eps*evals[0] ? undef : - [ each evec, pm*evec] ; + let( + indices = noncollinear_triple(points,error=false) + ) + indices==[] ? undef : + let( + p1 = points[indices[0]], + p2 = points[indices[1]], + p3 = points[indices[2]], + plane = plane3pt(p1,p2,p3) + ) + fast || points_on_plane(points,plane,eps=eps) ? plane : undef; // Function: plane_from_polygon() @@ -982,16 +945,17 @@ function plane_from_points(points,fast=false, eps=EPSILON) = // #stroke(xyzpath,closed=true); // cp = centroid(xyzpath); // move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(); -function plane_from_polygon(poly, fast=false, eps=EPSILON) = +function plane_from_polygon(poly, fast=false, eps=EPSILON) = assert( is_path(poly,dim=3), "Invalid polygon." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) : - let( triple = sort(noncollinear_triple(poly,error=false)) ) - triple==[] ? [] : - let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) - fast? plane: points_on_plane(poly, plane, eps=eps)? plane: []; - - + let( + poly = deduplicate(poly), + n = polygon_normal(poly), + plane = [n.x, n.y, n.z, n*poly[0]] + ) + fast? plane: coplanar(poly,eps=eps)? plane: []; + + // Function: plane_normal() // Usage: // plane_normal(plane); @@ -1288,11 +1252,9 @@ function coplanar(points, eps=EPSILON) = len(points)<=2 ? false : let( ip = noncollinear_triple(points,error=false,eps=eps) ) ip == [] ? false : - let( - plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), - normal = point3d(plane), - pt_nrm = points*normal ) - abs(max(max(pt_nrm)-plane[3], -min(pt_nrm)+plane[3])) < eps; + let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), + normal = point3d(plane) ) + max( points*normal ) - plane[3]< eps*norm(normal); // Function: points_on_plane() @@ -1703,11 +1665,11 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist) Date: Sat, 10 Apr 2021 21:07:23 +0100 Subject: [PATCH 03/12] Review of geometry.scad for speed --- geometry.scad | 107 +++++++++++++++++++++++++++------------ tests/test_geometry.scad | 49 +++++++++--------- 2 files changed, 99 insertions(+), 57 deletions(-) diff --git a/geometry.scad b/geometry.scad index 21b4f29..80f6756 100644 --- a/geometry.scad +++ b/geometry.scad @@ -888,11 +888,49 @@ function plane3pt_indexed(points, i1, i2, i3) = // Example: // plane_from_normal([0,0,1], [2,2,2]); // Returns the xy plane passing through the point (2,2,2) function plane_from_normal(normal, pt=[0,0,0]) = - assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0), - "Inputs `normal` and `pt` should 3d vectors/points and `normal` cannot be zero." ) - concat(normal, normal*pt) / norm(normal); + assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0), + "Inputs `normal` and `pt` should be 3d vectors/points and `normal` cannot be zero." ) + concat(normal, normal*pt) / norm(normal); +// Eigenvalues for a 3x3 symmetrical matrix in decreasing order +// Based on: https://en.wikipedia.org/wiki/Eigenvalue_algorithm +function _eigenvals_symm_3(M) = + let( p1 = pow(M[0][1],2) + pow(M[0][2],2) + pow(M[1][2],2) ) + (p1, ); @@ -900,7 +938,7 @@ function plane_from_normal(normal, pt=[0,0,0]) = // Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane, // that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1. // If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned. -// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. +// If `fast` is true, the polygon coplanarity check is skipped and a best fitted plane is returned. // Arguments: // points = The list of points to find the plane of. // fast = If true, don't verify that all points in the list are coplanar. Default: false @@ -914,17 +952,18 @@ function plane_from_normal(normal, pt=[0,0,0]) = function plane_from_points(points, fast=false, eps=EPSILON) = assert( is_path(points,dim=3), "Improper 3d point list." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - let( - indices = noncollinear_triple(points,error=false) - ) - indices==[] ? undef : - let( - p1 = points[indices[0]], - p2 = points[indices[1]], - p3 = points[indices[2]], - plane = plane3pt(p1,p2,p3) - ) - fast || points_on_plane(points,plane,eps=eps) ? plane : undef; + len(points) == 3 + ? let( plane = plane3pt(points[0],points[1],points[2]) ) + plane==[] ? undef : plane + : let( + cov_evals = _covariance_evals(points), + pm = cov_evals[0], + evals = cov_evals[1], + M = cov_evals[2], + evec = _eigenvec_symm_3(M,evals,i=2) ) +// echo(error_points_plane= abs(max(points*evec)-pm*evec), limit=eps) + !fast && abs(max(points*evec)-pm*evec)>eps*evals[0] ? undef : + [ each evec, pm*evec] ; // Function: plane_from_polygon() @@ -934,7 +973,8 @@ function plane_from_points(points, fast=false, eps=EPSILON) = // Given a 3D planar polygon, returns the normalized cartesian equation of its plane. // Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1. // If not all the points in the polygon are coplanar, then [] is returned. -// If `fast` is true, the polygon coplanarity check is skipped and the plane may not contain all polygon points. +// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned. +// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. // Arguments: // poly = The planar 3D polygon to find the plane of. // fast = If true, doesn't verify that all points in the polygon are coplanar. Default: false @@ -948,14 +988,13 @@ function plane_from_points(points, fast=false, eps=EPSILON) = function plane_from_polygon(poly, fast=false, eps=EPSILON) = assert( is_path(poly,dim=3), "Invalid polygon." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - let( - poly = deduplicate(poly), - n = polygon_normal(poly), - plane = [n.x, n.y, n.z, n*poly[0]] - ) - fast? plane: coplanar(poly,eps=eps)? plane: []; - - + len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) : + let( triple = sort(noncollinear_triple(poly,error=false)) ) + triple==[] ? [] : + let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) + fast? plane: points_on_plane(poly, plane, eps=eps)? plane: []; + + // Function: plane_normal() // Usage: // plane_normal(plane); @@ -1252,9 +1291,11 @@ function coplanar(points, eps=EPSILON) = len(points)<=2 ? false : let( ip = noncollinear_triple(points,error=false,eps=eps) ) ip == [] ? false : - let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), - normal = point3d(plane) ) - max( points*normal ) - plane[3]< eps*norm(normal); + let( + plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), + normal = point3d(plane), + pt_nrm = points*normal ) + abs(max(max(pt_nrm)-plane[3], -min(pt_nrm)+plane[3])) < eps; // Function: points_on_plane() @@ -1665,11 +1706,11 @@ function noncollinear_triple(points,error=true,eps=EPSILON) = n = (pb-pa)/nrm, distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] ) - max(distlist)0, "Invalid pointlist." ) // using is_matrix it accepts len(points)==1 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." ) - let( normal = point3d(plane), - pt_nrm = points*normal ) - abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3]))< eps*norm(normal); + _pointlist_greatest_distance(points,plane) < eps; // Function: in_front_of_plane() From 160e3f3edd390a91c28f1a373d413ca41ab2525e Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 11:35:38 +0100 Subject: [PATCH 05/12] Correction of polygon_area and a better polygon_normal --- geometry.scad | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/geometry.scad b/geometry.scad index 160ecd5..7f06614 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1784,7 +1784,7 @@ function polygon_area(poly, signed=false) = len(poly[0])==2 ? 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_points(poly) ) + : let( plane = plane_from_polygon(poly) ) plane==undef? undef : let( n = plane_normal(plane), @@ -2101,18 +2101,15 @@ function reverse_polygon(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`. +// It doesn't check for coplanarity. // Arguments: // poly = The list of 3D path points for the perimeter of the polygon. function polygon_normal(poly) = assert(is_path(poly,dim=3), "Invalid 3D polygon." ) - let( - poly = cleanup_path(poly), - p0 = poly[0], - n = sum([ - for (i=[1:1:len(poly)-2]) - cross(poly[i+1]-p0, poly[i]-p0) - ]) - ) unit(n,undef); + len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) : + let( triple = sort(noncollinear_triple(poly,error=false)) ) + triple==[] ? undef : + point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; function _split_polygon_at_x(poly, x) = From 80feb93c98fc03bf6bb0b16c03c5628e11637c2c Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 12:32:49 +0100 Subject: [PATCH 06/12] Change undef to [] as return of polygon functions --- geometry.scad | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/geometry.scad b/geometry.scad index 7f06614..d0e98df 100644 --- a/geometry.scad +++ b/geometry.scad @@ -453,7 +453,7 @@ function segment_closest_point(seg,pt) = // Usage: // line_from_points(points, [fast], [eps]); // Description: -// Given a list of 2 or more colinear points, returns a line containing them. +// 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 true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned. // Arguments: @@ -916,12 +916,12 @@ function _eigenvals_symm_3(M) = // using Cayley–Hamilton theorem according to: // https://en.wikipedia.org/wiki/Eigenvalue_algorithm function _eigenvec_symm_3(M,evals,i=0) = - let( - I = ident(3), - A = (M - evals[(i+1)%3]*I) * (M - evals[(i+2)%3]*I) , - k = max_index( [for(i=[0:2]) norm(A[i]) ]) - ) - norm(A[k])=0), "The tolerance should be a non-negative value." ) len(points) == 3 ? let( plane = plane3pt(points[0],points[1],points[2]) ) - plane==[] ? undef : plane + plane==[] ? [] : plane : let( covmix = _covariance_evec_eval(points), pm = covmix[0], @@ -976,7 +976,7 @@ function plane_from_points(points, fast=false, eps=EPSILON) = // Given a 3D planar polygon, returns the normalized cartesian equation of its plane. // Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1. // If not all the points in the polygon are coplanar, then [] is returned. -// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned. +// If `fast` is false and the points in the list are collinear or not coplanar, then [] is returned. // if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. // Arguments: // poly = The planar 3D polygon to find the plane of. @@ -1301,10 +1301,10 @@ function coplanar(points, eps=EPSILON) = // the maximum distance from points to the plane function _pointlist_greatest_distance(points,plane) = let( - normal = point3d(plane), - pt_nrm = points*normal + normal = point3d(plane), + pt_nrm = points*normal ) - abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3])) / norm(normal); + abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3])) / norm(normal); // Function: points_on_plane() @@ -1647,20 +1647,19 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) = // Function: circle_line_intersection() // Usage: -// isect = circle_line_intersection(c,r,line,,); -// isect = circle_line_intersection(c,d,line,,); +// isect = circle_line_intersection(c,,,,); // Description: // Find intersection points between a 2d circle and a line, ray or segment specified by two points. // By default the line is unbounded. // Arguments: // c = center of circle // r = radius of circle +// --- +// d = diameter of circle // line = two points defining the unbounded line // bounded = false for unbounded line, true for a segment, or a vector [false,true] or [true,false] to specify a ray with the first or second end unbounded. Default: false // eps = epsilon used for identifying the case with one solution. Default: 1e-9 -// --- -// d = diameter of circle -function circle_line_intersection(c,r,line,d,bounded=false,eps=EPSILON) = +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(is_vector(c,2), "Circle center must be a 2-vector") @@ -2100,7 +2099,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. @@ -2108,7 +2107,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]])) ; From f9e7baa5157c7970eb5fbd36c46fbe68821875f1 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 12:33:40 +0100 Subject: [PATCH 07/12] Review of is_matrix --- math.scad | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/math.scad b/math.scad index 43b1347..1f6cb92 100644 --- a/math.scad +++ b/math.scad @@ -985,10 +985,8 @@ function determinant(M) = function is_matrix(A,m,n,square=false) = is_list(A) && (( is_undef(m) && len(A) ) || len(A)==m) - && is_list(A[0]) - && (( is_undef(n) && len(A[0]) ) || len(A[0])==n) && (!square || len(A) == len(A[0])) - && is_vector(A[0]) + && is_vector(A[0],n) && is_consistent(A); From 376609d52f1ad754d2aadc1538dcd42bca825461 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 12:34:59 +0100 Subject: [PATCH 08/12] Cosmetics --- common.scad | 11 +-- paths.scad | 216 ++++++++++++++++++++++++++-------------------------- 2 files changed, 114 insertions(+), 113 deletions(-) diff --git a/common.scad b/common.scad index e24c744..6d2d44f 100644 --- a/common.scad +++ b/common.scad @@ -205,7 +205,8 @@ function is_func(x) = version_num()>20210000 && is_function(x); // Description: // Tests whether input is a list of entries which all have the same list structure // and are filled with finite numerical data. You can optionally specify a required -// list structure with the pattern argument. It returns `true` for the empty list. +// list structure with the pattern argument. +// It returns `true` for the empty list regardless the value of the `pattern`. // Arguments: // list = list to check // pattern = optional pattern required to match @@ -293,7 +294,7 @@ function default(v,dflt=undef) = is_undef(v)? dflt : v; // v = The list whose items are being checked. // recursive = If true, sublists are checked recursively for defined values. The first sublist that has a defined item is returned. // Examples: -// val = first_defined([undef,7,undef,true]); // Returns: 1 +// val = first_defined([undef,7,undef,true]); // Returns: 7 function first_defined(v,recursive=false,_i=0) = _i 0? _corner_chamfer_path(p1, p2, p3, side=chamfer[i]) : - rounding[i] > 0? _corner_roundover_path(p1, p2, p3, r=rounding[i]) : - [p2] - ], - out = [ - if (!closed) path[0], - for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let( - p1 = select(path,i-1), - p2 = select(path,i), - crn1 = select(corner_paths,i-1), - crn2 = corner_paths[i], - l1 = norm(last(crn1)-p1), - l2 = norm(crn2[0]-p2), - needed = l1 + l2, - seglen = norm(p2-p1), - check = assert(seglen >= needed, str("Path segment ",i," is too short to fulfill rounding/chamfering for the adjacent corners.")) - ) each crn2, - if (!closed) last(path) - ] - ) deduplicate(out); + let ( + path = deduplicate(path,closed=true), + lp = len(path), + chamfer = is_undef(chamfer)? repeat(0,lp) : + is_vector(chamfer)? list_pad(chamfer,lp,0) : + is_num(chamfer)? repeat(chamfer,lp) : + assert(false, "Bad chamfer value."), + rounding = is_undef(rounding)? repeat(0,lp) : + is_vector(rounding)? list_pad(rounding,lp,0) : + is_num(rounding)? repeat(rounding,lp) : + assert(false, "Bad rounding value."), + corner_paths = [ + for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let( + p1 = select(path,i-1), + p2 = select(path,i), + p3 = select(path,i+1) + ) + chamfer[i] > 0? _corner_chamfer_path(p1, p2, p3, side=chamfer[i]) : + rounding[i] > 0? _corner_roundover_path(p1, p2, p3, r=rounding[i]) : + [p2] + ], + out = [ + if (!closed) path[0], + for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let( + p1 = select(path,i-1), + p2 = select(path,i), + crn1 = select(corner_paths,i-1), + crn2 = corner_paths[i], + l1 = norm(last(crn1)-p1), + l2 = norm(crn2[0]-p2), + needed = l1 + l2, + seglen = norm(p2-p1), + check = assert(seglen >= needed, str("Path segment ",i," is too short to fulfill rounding/chamfering for the adjacent corners.")) + ) each crn2, + if (!closed) last(path) + ] + ) deduplicate(out); function _corner_chamfer_path(p1, p2, p3, dist1, dist2, side, angle) = - let( - v1 = unit(p1 - p2), - v2 = unit(p3 - p2), - n = vector_axis(v1,v2), - ang = vector_angle(v1,v2), - path = (is_num(dist1) && is_undef(dist2) && is_undef(side))? ( - // dist1 & optional angle - assert(dist1 > 0) - let(angle = default(angle,(180-ang)/2)) - assert(is_num(angle)) - assert(angle > 0 && angle < 180) - let( - pta = p2 + dist1*v1, - a3 = 180 - angle - ang - ) assert(a3>0, "Angle too extreme.") - let( - side = sin(angle) * dist1/sin(a3), - ptb = p2 + side*v2 - ) [pta, ptb] - ) : (is_undef(dist1) && is_num(dist2) && is_undef(side))? ( - // dist2 & optional angle - assert(dist2 > 0) - let(angle = default(angle,(180-ang)/2)) - assert(is_num(angle)) - assert(angle > 0 && angle < 180) - let( - ptb = p2 + dist2*v2, - a3 = 180 - angle - ang - ) assert(a3>0, "Angle too extreme.") - let( - side = sin(angle) * dist2/sin(a3), - pta = p2 + side*v1 - ) [pta, ptb] - ) : (is_undef(dist1) && is_undef(dist2) && is_num(side))? ( - // side & optional angle - assert(side > 0) - let(angle = default(angle,(180-ang)/2)) - assert(is_num(angle)) - assert(angle > 0 && angle < 180) - let( - a3 = 180 - angle - ang - ) assert(a3>0, "Angle too extreme.") - let( - dist1 = sin(a3) * side/sin(ang), - dist2 = sin(angle) * side/sin(ang), - pta = p2 + dist1*v1, - ptb = p2 + dist2*v2 - ) [pta, ptb] - ) : (is_num(dist1) && is_num(dist2) && is_undef(side) && is_undef(side))? ( - // dist1 & dist2 - assert(dist1 > 0) - assert(dist2 > 0) - let( - pta = p2 + dist1*v1, - ptb = p2 + dist2*v2 - ) [pta, ptb] - ) : ( - assert(false,"Bad arguments.") - ) - ) path; + let( + v1 = unit(p1 - p2), + v2 = unit(p3 - p2), + n = vector_axis(v1,v2), + ang = vector_angle(v1,v2), + path = (is_num(dist1) && is_undef(dist2) && is_undef(side))? ( + // dist1 & optional angle + assert(dist1 > 0) + let(angle = default(angle,(180-ang)/2)) + assert(is_num(angle)) + assert(angle > 0 && angle < 180) + let( + pta = p2 + dist1*v1, + a3 = 180 - angle - ang + ) assert(a3>0, "Angle too extreme.") + let( + side = sin(angle) * dist1/sin(a3), + ptb = p2 + side*v2 + ) [pta, ptb] + ) : (is_undef(dist1) && is_num(dist2) && is_undef(side))? ( + // dist2 & optional angle + assert(dist2 > 0) + let(angle = default(angle,(180-ang)/2)) + assert(is_num(angle)) + assert(angle > 0 && angle < 180) + let( + ptb = p2 + dist2*v2, + a3 = 180 - angle - ang + ) assert(a3>0, "Angle too extreme.") + let( + side = sin(angle) * dist2/sin(a3), + pta = p2 + side*v1 + ) [pta, ptb] + ) : (is_undef(dist1) && is_undef(dist2) && is_num(side))? ( + // side & optional angle + assert(side > 0) + let(angle = default(angle,(180-ang)/2)) + assert(is_num(angle)) + assert(angle > 0 && angle < 180) + let( + a3 = 180 - angle - ang + ) assert(a3>0, "Angle too extreme.") + let( + dist1 = sin(a3) * side/sin(ang), + dist2 = sin(angle) * side/sin(ang), + pta = p2 + dist1*v1, + ptb = p2 + dist2*v2 + ) [pta, ptb] + ) : (is_num(dist1) && is_num(dist2) && is_undef(side) && is_undef(side))? ( + // dist1 & dist2 + assert(dist1 > 0) + assert(dist2 > 0) + let( + pta = p2 + dist1*v1, + ptb = p2 + dist2*v2 + ) [pta, ptb] + ) : ( + assert(false,"Bad arguments.") + ) + ) path; function _corner_roundover_path(p1, p2, p3, r, d) = - let( - r = get_radius(r=r,d=d,dflt=undef), - res = circle_2tangents(p1, p2, p3, r=r, tangents=true), - cp = res[0], - n = res[1], - tp1 = res[2], - ang = res[4]+res[5], - steps = floor(segs(r)*ang/360+0.5), - step = ang / steps, - path = [for (i=[0:1:steps]) move(cp, p=rot(a=-i*step, v=n, p=tp1-cp))] - ) path; + let( + r = get_radius(r=r,d=d,dflt=undef), + res = circle_2tangents(p1, p2, p3, r=r, tangents=true), + cp = res[0], + n = res[1], + tp1 = res[2], + ang = res[4]+res[5], + steps = floor(segs(r)*ang/360+0.5), + step = ang / steps, + path = [for (i=[0:1:steps]) move(cp, p=rot(a=-i*step, v=n, p=tp1-cp))] + ) path; From c5799d65396a1cb711363f38b6b9f40d2877a2f4 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 12:36:05 +0100 Subject: [PATCH 09/12] Review of vnf_centroid() --- tests/test_vnf.scad | 3 ++- vnf.scad | 16 ++++++---------- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/tests/test_vnf.scad b/tests/test_vnf.scad index bb171d4..b83775e 100644 --- a/tests/test_vnf.scad +++ b/tests/test_vnf.scad @@ -74,7 +74,8 @@ module test_vnf_centroid() { assert_approx(vnf_centroid(cube(100, anchor=TOP)), [0,0,-50]); assert_approx(vnf_centroid(sphere(d=100, anchor=CENTER, $fn=36)), [0,0,0]); assert_approx(vnf_centroid(sphere(d=100, anchor=BOT, $fn=36)), [0,0,50]); -} + ellipse = xscale(2, p=circle($fn=24, r=3)); + assert_approx(vnf_centroid(path_sweep(pentagon(r=1), path3d(ellipse), closed=true)),[0,0,0]);} test_vnf_centroid(); diff --git a/vnf.scad b/vnf.scad index 9c05e0f..7d1c021 100644 --- a/vnf.scad +++ b/vnf.scad @@ -450,18 +450,13 @@ function vnf_volume(vnf) = // Returns the centroid of the given manifold VNF. The VNF must describe a valid polyhedron with consistent face direction and // no holes; otherwise the results are undefined. -// Divide the solid up into tetrahedra with the origin as one vertex. The centroid of a tetrahedron is the average of its vertices. +// Divide the solid up into tetrahedra with the origin as one vertex. +// The centroid of a tetrahedron is the average of its vertices. // The centroid of the total is the volume weighted average. function vnf_centroid(vnf) = + assert(is_vnf(vnf) && len(vnf[0])!=0 ) let( verts = vnf[0], - vol = sum([ - for(face=vnf[1], j=[1:1:len(face)-2]) let( - v0 = verts[face[0]], - v1 = verts[face[j]], - v2 = verts[face[j+1]] - ) cross(v2,v1)*v0 - ]), pos = sum([ for(face=vnf[1], j=[1:1:len(face)-2]) let( v0 = verts[face[0]], @@ -469,10 +464,11 @@ function vnf_centroid(vnf) = v2 = verts[face[j+1]], vol = cross(v2,v1)*v0 ) - (v0+v1+v2)*vol + [ vol, (v0+v1+v2)*vol ] ]) ) - pos/vol/4; + assert(!approx(pos[0],0, EPSILON), "The vnf has self-intersections.") + pos[1]/pos[0]/4; function _triangulate_planar_convex_polygons(polys) = From 5062ee1605db3827d703b3242467e8109a9aa44b Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 12:37:49 +0100 Subject: [PATCH 10/12] Minor change in is_matrix_symmetric --- arrays.scad | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/arrays.scad b/arrays.scad index ac97444..604186b 100644 --- a/arrays.scad +++ b/arrays.scad @@ -41,6 +41,7 @@ function is_homogeneous(l, depth=10) = [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) ) 0 ]; function is_homogenous(l, depth=10) = is_homogeneous(l, depth); + function _same_type(a,b, depth) = (depth==0) || @@ -50,7 +51,7 @@ function _same_type(a,b, depth) = (is_string(a) && is_string(b)) || (is_list(a) && is_list(b) && len(a)==len(b) && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); - + // Function: select() // Topics: List Handling @@ -97,7 +98,6 @@ function select(list, start, end) = // Function: slice() -// Topics: List Handling // Usage: // list = slice(list,s,e); // Description: @@ -476,7 +476,7 @@ function reverse(x) = // l9 = list_rotate([1,2,3,4,5],6); // Returns: [2,3,4,5,1] function list_rotate(list,n=1) = assert(is_list(list)||is_string(list), "Invalid list or string.") - assert(is_finite(n), "Invalid number") + assert(is_int(n), "The rotation number should be integer") let ( ll = len(list), n = ((n % ll) + ll) % ll, @@ -990,7 +990,7 @@ function _sort_vectors(arr, idxlist, _i=0) = _sort_vectors(equal, idxlist, _i+1), _sort_vectors(greater, idxlist, _i ) ); - + // sorting using compare_vals(); returns indexed list when `indexed==true` function _sort_general(arr, idx=undef, indexed=false) = (len(arr)<=1) ? arr : @@ -1332,6 +1332,8 @@ function permutations(l,n=2) = // pairs = zip(a,b); // triples = zip(a,b,c); // quads = zip([LIST1,LIST2,LIST3,LIST4]); +// Topics: List Handling, Iteration +// See Also: zip_long() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1357,6 +1359,8 @@ function zip(a,b,c) = // pairs = zip_long(a,b); // triples = zip_long(a,b,c); // quads = zip_long([LIST1,LIST2,LIST3,LIST4]); +// Topics: List Handling, Iteration +// See Also: zip() // Description: // Zips together two or more lists into a single list. For example, if you have two // lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]]. @@ -1526,7 +1530,6 @@ function subindex(M, idx) = // [[4,2], 91, false], // [6, [3,4], undef]]; // submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]] - function submatrix(M,idx1,idx2) = [for(i=idx1) [for(j=idx2) M[i][j] ] ]; @@ -1629,7 +1632,6 @@ function block_matrix(M) = assert(badrows==[], "Inconsistent or invalid input") bigM; - // Function: diagonal_matrix() // Usage: // mat = diagonal_matrix(diag, ); @@ -1855,7 +1857,7 @@ function transpose(arr, reverse=false) = // A = matrix to test // eps = epsilon for comparing equality. Default: 1e-12 function is_matrix_symmetric(A,eps=1e-12) = - approx(A,transpose(A)); + approx(A,transpose(A), eps); // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap From dab6407b0f0bcf610e81ffd87fba0e36cdad065c Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 13:01:06 +0100 Subject: [PATCH 11/12] Minor tweak --- geometry.scad | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/geometry.scad b/geometry.scad index d0e98df..26fa9ad 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1772,7 +1772,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. @@ -1784,9 +1784,9 @@ 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? undef : + plane==[]? [] : let( - n = plane_normal(plane), + n = plane_normal(plane), total = sum([ for(i=[1:1:len(poly)-2]) let( From 41616872fe0185627ce014b2cc791642c824d726 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 11 Apr 2021 23:45:21 +0100 Subject: [PATCH 12/12] Correction in is_convex_polygon --- geometry.scad | 11 ++++++----- tests/test_geometry.scad | 1 + 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/geometry.scad b/geometry.scad index 26fa9ad..6995f7f 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1304,7 +1304,7 @@ function _pointlist_greatest_distance(points,plane) = normal = point3d(plane), pt_nrm = points*normal ) - abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3])) / norm(normal); + abs(max( max(pt_nrm) - plane[3], -min(pt_nrm) + plane[3])) / norm(normal); // Function: points_on_plane() @@ -1786,7 +1786,7 @@ function polygon_area(poly, signed=false) = : let( plane = plane_from_polygon(poly) ) plane==[]? [] : let( - n = plane_normal(plane), + n = plane_normal(plane), total = sum([ for(i=[1:1:len(poly)-2]) let( @@ -1808,25 +1808,26 @@ function polygon_area(poly, signed=false) = // 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) = +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(max(crosses)) && !approx(min(crosses)), "The points are collinear" ) + ? 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(maxc-minc), "The points are collinear" ) + assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" ) minc>=0 || maxc<=0; diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 69d4705..ccb4ac7 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -847,6 +847,7 @@ module test_is_convex_polygon() { assert(is_convex_polygon(circle(r=50,$fn=1000))); assert(is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50))))); assert(!is_convex_polygon([[1,1],[0,0],[-1,1],[-1,-1],[1,-1]])); + assert(!is_convex_polygon([for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]])); // spiral } *test_is_convex_polygon();