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/beziers.scad b/beziers.scad index bd6a8f0..7cda75b 100644 --- a/beziers.scad +++ b/beziers.scad @@ -1222,6 +1222,193 @@ function bezier_patch(patch, splinesteps=16, vnf=EMPTY_VNF, style="default") = ) vnf; + +// Function: bezier_patch_degenerate() +// Usage: +// vnf = bezier_patch_degenerate(patch, , ); +// vnf_edges = bezier_patch_degenerate(patch, , , return_edges=true); +// Description: +// Returns a VNF for a degenerate rectangular bezier patch where some of the corners of the patch are +// equal. If the resulting patch has no faces then returns an empty VNF. Note that due to the degeneracy, +// the shape of the patch can be triangular even though the actual underlying patch is a rectangle. This is +// a different method for creating triangular bezier patches than the triangular patch. +// If you specify return_edges then the return is a list whose first element is the vnf and whose second +// element lists the edges in the order [left, right, top, bottom], where each list is a list of the actual +// point values, but possibly only a single point if that edge is degenerate. +// The method checks for various types of degeneracy and uses a triangular or partly triangular array of sample points. +// See examples below for the types of degeneracy detected and how the patch is sampled for those cases. +// Note that splinesteps is the same for both directions of the patch, so it cannot be an array. +// Arguments: +// patch = Patch to process +// splinesteps = Number of segments to produce on each side. Default: 16 +// reverse = reverse direction of faces. Default: false +// return_edges = if true return the points on the four edges: [left, right, top, bottom]. Default: false +// Example: This quartic patch is degenerate at one corner, where a row of control points are equal. Processing this degenerate patch normally produces excess triangles near the degenerate point. +// splinesteps=8; +// patch=[ +// repeat([-12.5, 12.5, 15],5), +// [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]], +// [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]], +// [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]], +// [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]] +// ]; +// vnf_wireframe((bezier_patch(patch, splinesteps)),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// Example: With bezier_patch_degenerate the degenerate point does not have excess triangles. The top half of the patch decreases the number of sampled points by 2 for each row. +// splinesteps=8; +// patch=[ +// repeat([-12.5, 12.5, 15],5), +// [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]], +// [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]], +// [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]], +// [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]] +// ]; +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// Example: With splinesteps odd you get one "odd" row where the point count decreases by 1 instead of 2. You may prefer even values for splinesteps to avoid this. +// splinesteps=7; +// patch=[ +// repeat([-12.5, 12.5, 15],5), +// [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]], +// [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]], +// [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]], +// [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]] +// ]; +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// Example: A more extreme degeneracy occurs when the top half of a patch is degenerate to a line. (For odd length patches the middle row must be degenerate to trigger this style.) In this case the number of points in each row decreases by 1 for every row. It doesn't matter of splinesteps is odd or even. +// splinesteps=8; +// patch = [[[10, 0, 0], [10, -10.4, 0], [10, -20.8, 0], [1.876, -14.30, 0], [-6.24, -7.8, 0]], +// [[5, 0, 0], [5, -5.2, 0], [5, -10.4, 0], [0.938, -7.15, 0], [-3.12, -3.9, 0]], +// repeat([0,0,0],5), +// repeat([0,0,5],5), +// repeat([0,0,10],5) +// ]; +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// Example: Here is a degenerate cubic patch. +// splinesteps=8; +// patch = [ [ [-20,0,0], [-10,0,0],[0,10,0],[0,20,0] ], +// [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10]], +// [ [-10,0,20], [-5,0,20], [0,5,20], [0,10,20]], +// repeat([0,0,30],4) +// ]; +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// Example: A more extreme degenerate cubic patch, where two rows are equal. +// splinesteps=8; +// patch = [ [ [-20,0,0], [-10,0,0],[0,10,0],[0,20,0] ], +// [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10] ], +// repeat([-10,10,20],4), +// repeat([-10,10,30],4) +// ]; +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// Example: Quadratic patch degenerate at the right side: +// splinesteps=8; +// patch = [[[0, -10, 0],[10, -5, 0],[20, 0, 0]], +// [[0, 0, 0], [10, 0, 0], [20, 0, 0]], +// [[0, 0, 10], [10, 0, 5], [20, 0, 0]]]; +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +// Example: Cubic patch degenerate at both ends. In this case the point count changes by 2 at every row. +// splinesteps=8; +// patch = [ +// repeat([10,-10,0],4), +// [ [-20,0,0], [-1,0,0],[0,10,0],[0,20,0] ], +// [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10] ], +// repeat([-10,10,20],4), +// ]; +// vnf_wireframe(bezier_patch_degenerate(patch, splinesteps),d=0.1); +// color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9); +function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_edges=false) = + !return_edges ? bezier_patch_degenerate(patch, splinesteps, reverse, true)[0] : + assert(is_rectpatch(patch), "Must supply rectangular bezier patch") + assert(is_int(splinesteps) && splinesteps>=3, "splinesteps must be an integer 3 or larger") + let( + row_degen = [for(row=patch) all_equal(row)], + col_degen = [for(col=transpose(patch)) all_equal(col)], + + top_degen = row_degen[0], + bot_degen = last(row_degen), + left_degen = col_degen[0], + right_degen = last(col_degen), + samplepts = lerpn(0,1,splinesteps+1) + ) + all(row_degen) && all(col_degen) ? // fully degenerate case + [EMPTY_VNF, repeat([patch[0][0]],4)] : + all(row_degen) ? // degenerate to a line (top to bottom) + let(pts = bezier_points(subindex(patch,0), samplepts)) + [EMPTY_VNF, [pts,pts,[pts[0]],[last(pts)]]] : + all(col_degen) ? // degenerate to a line (left to right) + let(pts = bezier_points(patch[0], samplepts)) + [EMPTY_VNF, [[pts[0]], [last(pts)], pts, pts]] : + !top_degen && !bot_degen && !left_degen && !right_degen ? // non-degenerate case + let(pts = bezier_patch_points(patch, samplepts, samplepts)) + [ + vnf_vertex_array(pts, reverse=!reverse), + [subindex(pts,0), subindex(pts,len(pts)-1), pts[0], last(pts)] + ] : + top_degen && bot_degen ? + let( + rowcount = [ + each list([3:2:splinesteps]), + if (splinesteps%2==0) splinesteps+1, + each reverse(list([3:2:splinesteps])) + ], + bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)], + pts = [ + [bpatch[0][0]], + for(j=[0:splinesteps-2]) bezier_points(subindex(bpatch,j+1), lerpn(0,1,rowcount[j])), + [last(bpatch[0])] + ], + vnf = vnf_tri_array(pts, reverse=reverse) + ) [ + vnf, + [ + subindex(pts,0), + [for(row=pts) last(row)], + pts[0], + last(pts), + ] + ] : + bot_degen ? // only bottom is degenerate + let( + result = bezier_patch_degenerate(reverse(patch), splinesteps=splinesteps, reverse=!reverse, return_edges=true) + ) + [ + result[0], + [reverse(result[1][0]), reverse(result[1][1]), (result[1][3]), (result[1][2])] + ] : + top_degen ? // only top is degenerate + let( + full_degen = len(patch)>=4 && all(select(row_degen,1,ceil(len(patch)/2-1))), + rowmax = full_degen ? count(splinesteps+1) : + [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps], + bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)], + pts = [ + [bpatch[0][0]], + for(j=[1:splinesteps]) bezier_points(subindex(bpatch,j), lerpn(0,1,rowmax[j]+1)) + ], + vnf = vnf_tri_array(pts, reverse=reverse) + ) [ + vnf, + [ + subindex(pts,0), + [for(row=pts) last(row)], + pts[0], + last(pts), + ] + ] : + // must have left or right degeneracy, so transpose and recurse + let( + result = bezier_patch_degenerate(transpose(patch), splinesteps=splinesteps, reverse=!reverse, return_edges=true) + ) + [result[0], + select(result[1],[2,3,0,1]) + ]; + + function _tri_count(n) = (n*(1+n))/2; 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. +// 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); @@ -1074,6 +1116,7 @@ function distance_from_plane(plane, point) = let( plane = normalize_plane(plane) ) point3d(plane)* point - plane[3]; + // Returns [POINT, U] if line intersects plane at one point. // Returns [LINE, undef] if the line is on the plane. // Returns undef if line is parallel to, but not on the given plane. @@ -1251,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() @@ -1269,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() @@ -1596,23 +1645,21 @@ 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") @@ -1629,7 +1676,6 @@ function circle_line_intersection(c,r,line,d,bounded=false,eps=EPSILON) = let( offset = sqrt(r*r-d*d), uvec=unit(line[1]-line[0]) ) [closest-offset*uvec, closest+offset*uvec] - ) [for(p=isect) if ((!bounded[0] || (p-line[0])*(line[1]-line[0])>=0) @@ -1637,7 +1683,6 @@ function circle_line_intersection(c,r,line,d,bounded=false,eps=EPSILON) = - // Section: Pointlists @@ -1667,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)=0 && sign(c)==s; - i = i-1, - c = i<0? 0: cross(poly[(i+1)%l]-poly[i],poly[(i+2)%l]-poly[(i+1)%l]), - s = s==0 ? sign(c) : s - ) i - ])== l; +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() @@ -1854,10 +1904,10 @@ function reindex_polygon(reference, poly, return_error=false) = polygon_is_clockwise(reference) ? clockwise_polygon(poly) : ccw_polygon(poly), - I = [for(i=[0:N-1]) 1], + I = [for(i=reference) 1], val = [ for(k=[0:N-1]) -           [for(i=[0:N-1]) -              (reference[i]*poly[(i+k)%N]) ] ]*I, + [for(i=[0:N-1]) + (reference[i]*poly[(i+k)%N]) ] ]*I, optimal_poly = polygon_shift(fixpoly, max_index(val)) ) return_error? [optimal_poly, min(poly*(I*poly)-2*val)] : @@ -2006,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() @@ -2050,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 43b1347..408afe2 100644 --- a/math.scad +++ b/math.scad @@ -583,7 +583,7 @@ function lcm(a,b=[]) = function sum(v, dflt=0) = v==[]? dflt : assert(is_consistent(v), "Input to sum is non-numeric or inconsistent") - is_vector(v) || is_matrix(v) ? [for(i=[1:len(v)]) 1]*v : + is_vector(v) || is_matrix(v) ? [for(i=v) 1]*v : _sum(v,v[0]*0); function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1); @@ -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); @@ -1155,6 +1153,18 @@ function all_nonnegative(x) = false; +// Function all_equal() +// Usage: +// b = all_equal(vec,); +// Description: +// Returns true if all of the entries in vec are equal to each other, or approximately equal to each other if eps is set. +// Arguments: +// vec = vector to check +// eps = Set to tolerance for approximate equality. Default: 0 +function all_equal(vec,eps=0) = + eps==0 ? [for(v=vec) if (v!=vec[0]) v] == [] + : [for(v=vec) if (!approx(v,vec[0])) v] == []; + // Function: approx() // Usage: // b = approx(a,b,) diff --git a/paths.scad b/paths.scad index 74f999e..11d05dd 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/polyhedra.scad b/polyhedra.scad index ac0d575..0914a60 100644 --- a/polyhedra.scad +++ b/polyhedra.scad @@ -118,6 +118,7 @@ function _unique_groups(m) = [ // * `"great dodecahedron"` // * `"small stellated dodecahedron"` // * `"great stellated dodecahedron"` +// * `"small triambic icosahedron"` // // Arguments: // name = Name of polyhedron to create. @@ -147,7 +148,7 @@ function _unique_groups(m) = [ // Side Effects: // `$faceindex` - Index number of the face // `$face` - Coordinates of the face (2d if rotate_children==true, 3d if not) -// `$center` - Polyhedron center in the child coordinate system +// `$center` - Face center in the child coordinate system // // Examples: All of the available polyhedra by name in their native orientation // regular_polyhedron("tetrahedron", facedown=false); @@ -185,6 +186,7 @@ function _unique_groups(m) = [ // regular_polyhedron("great dodecahedron"); // regular_polyhedron("small stellated dodecahedron"); // regular_polyhedron("great stellated dodecahedron"); +// regular_polyhedron("small triambic icosahedron"); // Example: Third Archimedean solid // regular_polyhedron(type="archimedean", index=2); // Example(Med): Solids that have 8 or 10 vertex faces @@ -275,6 +277,10 @@ function _unique_groups(m) = [ // %sphere(r=.98); // regular_polyhedron("pentagonal hexecontahedron", or=1,facedown=false); // } +// Example: Stellate an Archimedian solid, which has mixed faces +// regular_polyhedron("truncated icosahedron",stellate=1.5,or=1); +// Example: Stellate a Catalan solid where faces are not regular +// regular_polyhedron("triakis tetrahedron",stellate=0.5,or=1); module regular_polyhedron( name=undef, index=undef, @@ -330,20 +336,19 @@ module regular_polyhedron( } } } + translate(translation) if ($children>0) { maxrange = repeat ? len(faces)-1 : $children-1; for(i=[0:1:maxrange]) { // Would like to orient so an edge (longest edge?) is parallel to x axis - facepts = move(translation, p=select(scaled_points, faces[i])); - center = mean(facepts); - rotatedface = rot(from=face_normals[i], to=[0,0,1], p=move(-center, p=facepts)); - clockwise = sortidx([for(pt=rotatedface) -atan2(pt.y,pt.x)]); - $face = rotate_children? - path2d(select(rotatedface,clockwise)) : - select(move(-center,p=facepts), clockwise); + facepts = select(scaled_points, faces[i]); + $center = -mean(facepts); + cfacepts = move($center, p=facepts); + $face = rotate_children + ? path2d(rot(from=face_normals[i], to=[0,0,1], p=cfacepts)) + : cfacepts; $faceindex = i; - $center = -translation-center; - translate(center) + translate(-$center) if (rotate_children) { rot(from=[0,0,1], to=face_normals[i]) children(i % $children); @@ -537,6 +542,7 @@ _stellated_polyhedra_ = [ ["great dodecahedron", "icosahedron", -sqrt(5/3-PHI)], ["small stellated dodecahedron", "dodecahedron", sqrt((5+2*sqrt(5))/5)], ["great stellated dodecahedron", "icosahedron", sqrt(2/3+PHI)], + ["small triambic icosahedron", "icosahedron", sqrt(3/5) - 1/sqrt(3)] ]; @@ -699,7 +705,7 @@ function regular_polyhedron_info( face_normals, radius_scale*entry[in_radius] ] : - info == "vnf" ? [move(translation,p=scaled_points), stellate ? faces : face_triangles] : + info == "vnf" ? [move(translation,p=scaled_points), faces] : info == "vertices" ? move(translation,p=scaled_points) : info == "faces" ? faces : info == "face normals" ? face_normals : @@ -770,15 +776,28 @@ function _facenormal(pts, face) = unit(cross(pts[face[2]]-pts[face[0]], pts[face // hull() function returns triangulated faces. This function identifies the vertices that belong to each face // by grouping together the face triangles that share normal vectors. The output gives the face polygon -// point indices in arbitrary order (not usable as input to a polygon call) and a normal vector. +// point indices in arbitrary order (not usable as input to a polygon call) and a normal vector. Finally +// the faces are ordered based on angle with their center (will always give a valid order for convex polygons). +// Final return is [ordered_faces, facenormals] where the first is a list of indices into the point list +// and the second is a list of vectors. function _full_faces(pts,faces) = let( normals = [for(face=faces) quant(_facenormal(pts,face),1e-12)], groups = _unique_groups(normals), faces = [for(entry=groups) unique(flatten(select(faces, entry)))], - facenormals = [for(entry=groups) normals[entry[0]]] - ) [faces, facenormals]; + facenormals = [for(entry=groups) normals[entry[0]]], + ordered_faces = [ + for(i=idx(faces)) + let( + facepts = select(pts, faces[i]), + center = mean(facepts), + rotatedface = rot(from=facenormals[i], to=[0,0,1], p=move(-center, p=facepts)), + clockwise = sortidx([for(pt=rotatedface) -atan2(pt.y,pt.x)]) + ) + select(faces[i],clockwise) + ] + ) [ordered_faces, facenormals]; // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/rounding.scad b/rounding.scad index b85a70c..d2493c6 100644 --- a/rounding.scad +++ b/rounding.scad @@ -1779,8 +1779,8 @@ function _rp_compute_patches(top, bot, rtop, rsides, ktop, ksides, concave) = // M = path3d(turtle(["left", 180, "length",3,"move", "left", "move", 3, "right", "move", "right", "move", 4, "right", "move", 3, "right", "move", 2])); // rounded_prism(M, apply(right(1)*scale(.75)*up(3),M), joint_top=0.5, joint_bot=0.2, joint_sides=[.2,1,1,0.5,1.5,.5,2], splinesteps=32); // Example: this example shows most of the different types of patches that rounded_prism creates. Note that some of the patches are close to interfering with each other across the top of the polyhedron, which would create an invalid result. -// N = apply(rot(180)*yscale(.8),turtle(["length",3,"left", "move", 2, "right", 135, "move", sqrt(2), "left", "move", sqrt(2), "right", 135, "move", 2])); -// rounded_prism(N, height=3, joint_bot=0.5, joint_top=1.25, joint_sides=[[1,1.75],0,.5,.5,2], debug=true); + N = apply(rot(180)*yscale(.8),turtle(["length",3,"left", "move", 2, "right", 135, "move", sqrt(2), "left", "move", sqrt(2), "right", 135, "move", 2])); + rounded_prism(N, height=3, joint_bot=0.5, joint_top=1.25, joint_sides=[[1,1.75],0,.5,.5,2], debug=true); // Example: This object has different scales on its different axies. Here is the largest symmetric rounding that fits. Note that the rounding is slightly smaller than the object dimensions because of roundoff error. // rounded_prism(square([100.1,30.1]), height=8.1, joint_top=4, joint_bot=4, joint_sides=15, k_sides=0.3, splinesteps=32); // Example: Using asymetric rounding enables a much more rounded form: @@ -1886,8 +1886,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b let( // Entries in the next two lists have the form [edges, vnf] where // edges is a list [leftedge, rightedge, topedge, botedge] - top_samples = [for(patch=top_patch) bezier_patch_degenerate(patch,splinesteps,reverse=true) ], - bot_samples = [for(patch=bot_patch) bezier_patch_degenerate(patch,splinesteps,reverse=false) ], + top_samples = [for(patch=top_patch) bezier_patch_degenerate(patch,splinesteps,reverse=true,return_edges=true) ], + bot_samples = [for(patch=bot_patch) bezier_patch_degenerate(patch,splinesteps,reverse=false,return_edges=true) ], leftidx=0, rightidx=1, topidx=2, @@ -1895,14 +1895,14 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b edge_points = [for(i=[0:N-1]) let( - top_edge = [ top_samples[i][0][rightidx], select(top_samples, i+1)[0][leftidx]], - bot_edge = [ select(bot_samples, i+1)[0][leftidx], bot_samples[i][0][rightidx]], - vert_edge = [ bot_samples[i][0][botidx], top_samples[i][0][botidx]] + top_edge = [ top_samples[i][1][rightidx], select(top_samples, i+1)[1][leftidx]], + bot_edge = [ select(bot_samples, i+1)[1][leftidx], bot_samples[i][1][rightidx]], + vert_edge = [ bot_samples[i][1][botidx], top_samples[i][1][botidx]] ) each [top_edge, bot_edge, vert_edge] ], faces = [ - [for(i=[0:N-1]) each top_samples[i][0][topidx]], - [for(i=[N-1:-1:0]) each reverse(bot_samples[i][0][topidx])], + [for(i=[0:N-1]) each top_samples[i][1][topidx]], + [for(i=[N-1:-1:0]) each reverse(bot_samples[i][1][topidx])], for(i=[0:N-1]) [ bot_patch[i][4][4], select(bot_patch,i+1)[4][0], @@ -1934,8 +1934,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b "Roundovers interfere with each other on bottom face: either input is self intersecting or top joint length is too large") assert(debug || (verify_vert==[] && verify_horiz==[]), "Curvature continuity failed") let( - vnf = vnf_merge([ each subindex(top_samples,1), - each subindex(bot_samples,1), + vnf = vnf_merge([ each subindex(top_samples,0), + each subindex(bot_samples,0), for(pts=edge_points) vnf_vertex_array(pts), vnf_triangulate(vnf_add_faces(EMPTY_VNF,faces)) ]) @@ -1943,116 +1943,6 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b debug ? [concat(top_patch, bot_patch), vnf] : vnf; -// This function takes a bezier patch as input and returns [edges, vnf], where -// edges = [leftedge, rightedge, topedge, bottomedge] -// gives the points along the edges of the patch, and the vnf is the patch vnf. -// It checks for various types of degeneracy and uses half or full triangular -// sampling on degenerate patches. - -function bezier_patch_degenerate(patch, splinesteps=16, reverse=false) = - assert(is_num(splinesteps), "splinesteps must be a number") - let( - top_degen = patch[0][0] == last(patch[0]), - bot_degen = last(patch)[0] == last(last(patch)), - left_degen = patch[0][0] == last(patch)[0], - right_degen = last(patch[0]) == last(last(patch)), - samplepts = count(splinesteps+1)/splinesteps - ) - top_degen && bot_degen && left_degen && right_degen ? // fully degenerate case - [repeat([patch[0][0]],4), EMPTY_VNF] : - top_degen && bot_degen ? // double degenerate case (top/bot) - let( - pts = bezier_points(subindex(patch,0), samplepts) - ) - [[pts,pts,[pts[0]],[last(pts)]], EMPTY_VNF] : - left_degen && right_degen ? // double degenerate case (sides) - let( - pts = bezier_points(patch[0], samplepts) - ) - [[[pts[0]], [last(pts)], pts, pts], EMPTY_VNF] : - !top_degen && !bot_degen ? // non-degenerate case - let( - k=echo("non-degenerate case"), - pts = bezier_patch_points(patch, samplepts, samplepts) - ) - [ - [subindex(pts,0), subindex(pts,len(pts)-1), pts[0], last(pts)], - vnf_vertex_array(pts, reverse=reverse) - ] : - bot_degen ? // only bottom is degenerate - let( - result = bezier_patch_degenerate(reverse(patch), splinesteps=splinesteps, reverse=!reverse) - ) - [ - [reverse(result[0][0]), reverse(result[0][1]), (result[0][3]), (result[0][2])], - result[1] - ] : - // at this point top_degen is true // only top is degenerate - let( - full_degen = patch[1][0] == last(patch[1]), - rowmax = full_degen ? count(splinesteps+1) : - [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps], - vbb=echo("single degenerate case"), - bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)], - pts = [ - [bpatch[0][0]], - for(j=[1:splinesteps]) bezier_points(subindex(bpatch,j), lerpn(0,1,rowmax[j]+1)) - ], - vnf = vnf_tri_array(pts, reverse=reverse) - ) [ - [ - subindex(pts,0), - [for(row=pts) last(row)], - pts[0], - last(pts), - ], - vnf - ]; - - -// This function produces a vnf with a triangulation for a list of rows -// where the number of points between rows differs by at most 2. -// It's a generalization of vnf_vertex_array. -function vnf_tri_array(points, row_wrap=false, reverse=false) = - let( - lens = [for(row=points) len(row)], - rowstarts = [0,each cumsum(lens)], - faces = - [for(i=[0:1:len(points) - 1 - (row_wrap ? 0 : 1)]) each - let( - rowstart = rowstarts[i], - nextrow = select(rowstarts,i+1), - delta = select(lens,i+1)-lens[i] - ) - delta == 0 ? - [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow], - for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1] : [j+rowstart+1, j+nextrow+1, j+nextrow]] : - delta == 1 ? - [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1], - for(j=[0:1:lens[i]-1]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow]] : - delta == -1 ? - [for(j=[0:1:lens[i]-3]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1]: [j+rowstart+1, j+nextrow+1, j+nextrow], - for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow]] : - let(count = floor((lens[i]-1)/2)) - delta == 2 ? - [ - for(j=[0:1:count-1]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1], // top triangles left - for(j=[count:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+2] : [j+rowstart, j+rowstart+1, j+nextrow+2], // top triangles right - for(j=[0:1:count]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow], // bot triangles left - for(j=[count+1:1:select(lens,i+1)-2]) reverse ? [j+rowstart-1, j+nextrow, j+nextrow+1] : [j+rowstart-1, j+nextrow+1, j+nextrow], // bot triangles right - ] : - delta == -2 ? - [ - for(j=[0:1:count-2]) reverse ? [j+nextrow, j+nextrow+1, j+rowstart+1] : [j+nextrow, j+rowstart+1, j+nextrow+1], - for(j=[count-1:1:lens[i]-4]) reverse ? [j+nextrow,j+nextrow+1,j+rowstart+2] : [j+nextrow,j+rowstart+2, j+nextrow+1], - for(j=[0:1:count-1]) reverse ? [j+nextrow, j+rowstart+1, j+rowstart] : [j+nextrow, j+rowstart, j+rowstart+1], - for(j=[count:1:select(lens,i+1)]) reverse ? [ j+nextrow-1, j+rowstart+1, j+rowstart]: [ j+nextrow-1, j+rowstart, j+rowstart+1], - ] : - assert(false,str("Unsupported row length difference of ",delta, " between row ",i," and ",i+1)) - ]) - [flatten(points), faces]; - - // Converts a 2d path to a path on a cylinder at radius r function _cyl_hole(r, path) = diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index a80894f..ccb4ac7 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(); @@ -844,7 +845,9 @@ module test_polygon_area() { module test_is_convex_polygon() { assert(is_convex_polygon([[1,1],[-1,1],[-1,-1],[1,-1]])); 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 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/version.scad b/version.scad index 89b3852..25bda53 100644 --- a/version.scad +++ b/version.scad @@ -6,7 +6,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,601]; +BOSL_VERSION = [2,0,605]; // Section: BOSL Library Version Functions diff --git a/vnf.scad b/vnf.scad index 9c05e0f..964d4e7 100644 --- a/vnf.scad +++ b/vnf.scad @@ -367,6 +367,87 @@ function vnf_vertex_array( ]); +// Function: vnf_tri_array() +// Usage: +// vnf = vnf_tri_array(points, , ) +// Description: +// Produces a vnf from an array of points where each row length can differ from the adjacent rows by up to 2 in length. This enables +// the construction of triangular VNF patches. The resulting VNF can be wrapped along the rows by setting `row_wrap` to true. +// Arguments: +// points = List of point lists for each row +// row_wrap = If true then add faces connecting the first row and last row. These rows must differ by at most 2 in length. +// reverse = Set this to reverse the direction of the faces +// Examples: Each row has one more point than the preceeding one. +// pts = [for(y=[1:1:10]) [for(x=[0:y-1]) [x,y,y]]]; +// vnf = vnf_tri_array(pts); +// vnf_wireframe(vnf,d=.1); +// color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9); +// Examples: Each row has one more point than the preceeding one. +// pts = [for(y=[0:2:10]) [for(x=[-y/2:y/2]) [x,y,y]]]; +// vnf = vnf_tri_array(pts); +// vnf_wireframe(vnf,d=.1); +// color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9); +// Example: Chaining two VNFs to construct a cone with one point length change between rows. +// pts1 = [for(z=[0:10]) path3d(arc(3+z,r=z/2+1, angle=[0,180]),10-z)]; +// pts2 = [for(z=[0:10]) path3d(arc(3+z,r=z/2+1, angle=[180,360]),10-z)]; +// vnf = vnf_tri_array(pts1, +// vnf=vnf_tri_array(pts2)); +// color("green")vnf_wireframe(vnf,d=.1); +// vnf_polyhedron(vnf); +// Example: Cone with length change two between rows +// pts1 = [for(z=[0:1:10]) path3d(arc(3+2*z,r=z/2+1, angle=[0,180]),10-z)]; +// pts2 = [for(z=[0:1:10]) path3d(arc(3+2*z,r=z/2+1, angle=[180,360]),10-z)]; +// vnf = vnf_tri_array(pts1, +// vnf=vnf_tri_array(pts2)); +// color("green")vnf_wireframe(vnf,d=.1); +// vnf_polyhedron(vnf); +// Example: Point count can change irregularly +// lens = [10,9,7,5,6,8,8,10]; +// pts = [for(y=idx(lens)) lerpn([-lens[y],y,y],[lens[y],y,y],lens[y])]; +// vnf = vnf_tri_array(pts); +// vnf_wireframe(vnf,d=.1); +// color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9); +function vnf_tri_array(points, row_wrap=false, reverse=false, vnf=EMPTY_VNF) = + let( + lens = [for(row=points) len(row)], + rowstarts = [0,each cumsum(lens)], + faces = + [for(i=[0:1:len(points) - 1 - (row_wrap ? 0 : 1)]) each + let( + rowstart = rowstarts[i], + nextrow = select(rowstarts,i+1), + delta = select(lens,i+1)-lens[i] + ) + delta == 0 ? + [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow], + for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1] : [j+rowstart+1, j+nextrow+1, j+nextrow]] : + delta == 1 ? + [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1], + for(j=[0:1:lens[i]-1]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow]] : + delta == -1 ? + [for(j=[0:1:lens[i]-3]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1]: [j+rowstart+1, j+nextrow+1, j+nextrow], + for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow]] : + let(count = floor((lens[i]-1)/2)) + delta == 2 ? + [ + for(j=[0:1:count-1]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1], // top triangles left + for(j=[count:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+2] : [j+rowstart, j+rowstart+1, j+nextrow+2], // top triangles right + for(j=[0:1:count]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow], // bot triangles left + for(j=[count+1:1:select(lens,i+1)-2]) reverse ? [j+rowstart-1, j+nextrow, j+nextrow+1] : [j+rowstart-1, j+nextrow+1, j+nextrow], // bot triangles right + ] : + delta == -2 ? + [ + for(j=[0:1:count-2]) reverse ? [j+nextrow, j+nextrow+1, j+rowstart+1] : [j+nextrow, j+rowstart+1, j+nextrow+1], + for(j=[count-1:1:lens[i]-4]) reverse ? [j+nextrow,j+nextrow+1,j+rowstart+2] : [j+nextrow,j+rowstart+2, j+nextrow+1], + for(j=[0:1:count-1]) reverse ? [j+nextrow, j+rowstart+1, j+rowstart] : [j+nextrow, j+rowstart, j+rowstart+1], + for(j=[count:1:select(lens,i+1)]) reverse ? [ j+nextrow-1, j+rowstart+1, j+rowstart]: [ j+nextrow-1, j+rowstart, j+rowstart+1], + ] : + assert(false,str("Unsupported row length difference of ",delta, " between row ",i," and ",(i+1)%len(points))) + ]) + vnf_merge(cleanup=true, [vnf, [flatten(points), faces]]); + + + // Module: vnf_polyhedron() // Usage: // vnf_polyhedron(vnf); @@ -450,18 +531,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 +545,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) =