diff --git a/arrays.scad b/arrays.scad index ac97444b..604186b3 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 e24c744e..6d2d44f6 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. +// that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane and norm([A,B,C])=1. +// If `fast` is false and the points in the list are collinear or not coplanar, then [] 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 +// fast = If true, don't verify the point coplanarity. 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)); @@ -914,17 +956,17 @@ 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==[] ? [] : plane + : let( + covmix = _covariance_evec_eval(points), + pm = covmix[0], + evec = covmix[1], + eval0 = covmix[2], + plane = [ each evec, pm*evec] ) + !fast && _pointlist_greatest_distance(points,plane)>eps*eval0 ? undef : + plane ; // Function: plane_from_polygon() @@ -934,7 +976,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 [] 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 +991,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 +1294,17 @@ 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]]) ) + _pointlist_greatest_distance(points,plane) < eps; + + +// the maximum distance from points to the plane +function _pointlist_greatest_distance(points,plane) = + let( + normal = point3d(plane), + pt_nrm = points*normal + ) + abs(max( max(pt_nrm) - plane[3], -min(pt_nrm) + plane[3])) / norm(normal); // Function: points_on_plane() @@ -1270,9 +1320,7 @@ function points_on_plane(points, plane, eps=EPSILON) = assert( _valid_plane(plane), "Invalid plane." ) assert( is_matrix(points,undef,3) && len(points)>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() @@ -1599,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") @@ -1665,11 +1712,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)=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; @@ -2008,7 +2056,7 @@ function point_in_polygon(point, poly, nonzero=true, eps=EPSILON) = // poly = The list of 2D path points for the perimeter of the polygon. function polygon_is_clockwise(poly) = assert(is_path(poly,dim=2), "Input should be a 2d path") - polygon_area(poly, signed=true)<0; + polygon_area(poly, signed=true)<-EPSILON; // Function: clockwise_polygon() @@ -2052,19 +2100,16 @@ 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. 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==[] ? [] : + point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ; function _split_polygon_at_x(poly, x) = diff --git a/math.scad b/math.scad index 43b13479..1f6cb920 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); diff --git a/paths.scad b/paths.scad index 74f999e3..11d05dda 100644 --- a/paths.scad +++ b/paths.scad @@ -454,120 +454,120 @@ function path_torsion(path, closed=false) = // path2 = path_chamfer_and_rounding(path, closed=true, chamfer=chamfs, rounding=rs); // stroke(path2, closed=true); function path_chamfer_and_rounding(path, closed, chamfer, rounding) = - 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); + 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; diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 6fe00c12..ccb4ac74 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -197,8 +197,8 @@ module test_plane_from_polygon(){ poly1 = [ rands(-1,1,3), rands(-1,1,3)+[2,0,0], rands(-1,1,3)+[0,2,2] ]; poly2 = concat(poly1, [sum(poly1)/3] ); info = info_str([["poly1 = ",poly1],["poly2 = ",poly2]]); - assert_std(plane_from_polygon(poly1),plane3pt(poly1[0],poly1[1],poly1[2]),info); - assert_std(plane_from_polygon(poly2),plane3pt(poly1[0],poly1[1],poly1[2]),info); + assert_approx(plane_from_polygon(poly1),plane3pt(poly1[0],poly1[1],poly1[2]),info); + assert_approx(plane_from_polygon(poly2),plane3pt(poly1[0],poly1[1],poly1[2]),info); } *test_plane_from_polygon(); @@ -208,8 +208,7 @@ module test_plane_from_normal(){ displ = normal*point; info = info_str([["normal = ",normal],["point = ",point],["displ = ",displ]]); assert_approx(plane_from_normal(normal,point)*[each point,-1],0,info); - assert_std(plane_from_normal(normal,point),normalize_plane([each normal,displ]),info); - assert_std(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]); + assert_approx(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]); } *test_plane_from_normal(); @@ -680,23 +679,23 @@ module test_triangle_area() { module test_plane3pt() { - assert_std(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]); - assert_std(plane3pt([2,0,20], [2,10,10], [2,0,0]), [1,0,0,2]); - assert_std(plane3pt([0,0,0], [10,0,10], [0,0,20]), [0,1,0,0]); - assert_std(plane3pt([0,2,0], [10,2,10], [0,2,20]), [0,1,0,2]); - assert_std(plane3pt([0,0,0], [10,10,0], [20,0,0]), [0,0,1,0]); - assert_std(plane3pt([0,0,2], [10,10,2], [20,0,2]), [0,0,1,2]); + assert_approx(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]); + assert_approx(plane3pt([2,0,20], [2,10,10], [2,0,0]), [1,0,0,2]); + assert_approx(plane3pt([0,0,0], [10,0,10], [0,0,20]), [0,1,0,0]); + assert_approx(plane3pt([0,2,0], [10,2,10], [0,2,20]), [0,1,0,2]); + assert_approx(plane3pt([0,0,0], [10,10,0], [20,0,0]), [0,0,1,0]); + assert_approx(plane3pt([0,0,2], [10,10,2], [20,0,2]), [0,0,1,2]); } *test_plane3pt(); module test_plane3pt_indexed() { pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ]; s13 = sqrt(1/3); - assert_std(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]); - assert_std(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]); - assert_std(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]); - assert_std(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]); - assert_std(plane3pt_indexed(pts, 0,2,1), [0,0,1,0]); + assert_approx(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]); + assert_approx(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]); + assert_approx(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]); + assert_approx(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]); + assert_approx(plane3pt_indexed(pts, 0,2,1), [0,0,1,0]); assert_approx(plane3pt_indexed(pts, 0,1,2), [0,0,-1,0]); assert_approx(plane3pt_indexed(pts, 3,2,1), [s13,s13,s13,10*s13]); assert_approx(plane3pt_indexed(pts, 1,2,3), [-s13,-s13,-s13,-10*s13]); @@ -709,18 +708,18 @@ module test_plane_from_points() { assert_std(plane_from_points([[0,0,0], [10,0,10], [0,0,20], [5,0,7]]), [0,1,0,0]); assert_std(plane_from_points([[0,2,0], [10,2,10], [0,2,20], [4,2,3]]), [0,1,0,2]); assert_std(plane_from_points([[0,0,0], [10,10,0], [20,0,0], [8,3,0]]), [0,0,1,0]); - assert_std(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]), [0,0,1,2]); + assert_std(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]), [0,0,1,2]); } *test_plane_from_points(); module test_plane_normal() { - assert_std(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]); - assert_std(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])), [1,0,0]); - assert_std(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])), [0,1,0]); - assert_std(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])), [0,1,0]); - assert_std(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])), [0,0,1]); - assert_std(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])), [0,0,1]); + assert_approx(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]); + assert_approx(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])), [1,0,0]); + assert_approx(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])), [0,1,0]); + assert_approx(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])), [0,1,0]); + assert_approx(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])), [0,0,1]); + assert_approx(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])), [0,0,1]); } *test_plane_normal(); @@ -780,7 +779,7 @@ module test_coplanar() { assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true); assert(coplanar([ [0,0,0],[1,0,1],[1,1,1], [0,1,2] ]) == false); assert(coplanar([ [0,0,0],[1,0,1],[1,1,2], [0,1,1] ]) == true); -} + } *test_coplanar(); @@ -836,7 +835,9 @@ module test_cleanup_path() { module test_polygon_area() { assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4)); assert(approx(polygon_area(circle(r=50,$fn=1000),signed=true), -PI*50*50, eps=0.1)); - assert(approx(polygon_area(rot([13,27,75],p=path3d(circle(r=50,$fn=1000),fill=23)),signed=true), PI*50*50, eps=0.1)); + assert(approx(polygon_area(rot([13,27,75], + p=path3d(circle(r=50,$fn=1000),fill=23)), + signed=true), -PI*50*50, eps=0.1)); } *test_polygon_area(); @@ -846,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(); diff --git a/tests/test_vnf.scad b/tests/test_vnf.scad index bb171d48..b83775e2 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 9c05e0f9..7d1c0212 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) =