grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2025-01-17 09:59:39 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'master' of github.com:revarbat/BOSL2 into revarbat_dev

Browse source
This commit is contained in:
Garth Minette 2021-04-12 00:42:41 -07:00
commit 6146cc3175
12 changed files with 610 additions and 373 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
arrays.scad
View file

@ -42,6 +42,7 @@ function is_homogeneous(l, depth=10) =
function is_homogenous(l, depth=10) = is_homogeneous(l, depth);
function _same_type(a,b, depth) =
(depth==0) ||
(is_undef(a) && is_undef(b)) ||
@ -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,
@ -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, <offdiag>);
@ -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

187
beziers.scad
View file

@ -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, <splinesteps>, <reverse>);
// vnf_edges = bezier_patch_degenerate(patch, <splinesteps>, <reverse>, 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;

11
common.scad
View file

@ -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<len(v) && (
is_undef(v[_i]) || (
@ -605,15 +606,15 @@ function segs(r) =
// Module: no_children()
// Topics: Error Checking
// Usage:
// no_children($children);
// Topics: Error Checking
// See Also: no_function(), no_module()
// Description:
// Assert that the calling module does not support children. Prints an error message to this effect and fails if children are present,
// as indicated by its argument.
// Arguments:
// $children = number of children the module has.
// See Also: no_function(), no_module()
// Example:
// module foo() {
// no_children($children);
@ -676,7 +677,7 @@ function _valstr(x) =
// expected = The value that was expected.
// info = Extra info to print out to make the error clearer.
// Example:
// assert_approx(1/3, 0.333333333333333, str("numer=",1,", demon=",3));
// assert_approx(1/3, 0.333333333333333, str("number=",1,", demon=",3));
module assert_approx(got, expected, info) {
no_children($children);
if (!approx(got, expected)) {

189
geometry.scad
View file

@ -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:
@ -889,21 +889,63 @@ function plane3pt_indexed(points, i1, i2, i3) =
// 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." )
"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<EPSILON)
? -sort(-[ M[0][0], M[1][1], M[2][2] ]) // diagonal matrix: eigenvals in decreasing order
: let( q = (M[0][0]+M[1][1]+M[2][2])/3,
B = (M - q*ident(3)),
dB = [B[0][0], B[1][1], B[2][2]],
p2 = dB*dB + 2*p1,
p = sqrt(p2/6),
r = det3(B/p)/2,
ph = acos(constrain(r,-1,1))/3,
e1 = q + 2*p*cos(ph),
e3 = q + 2*p*cos(ph+120),
e2 = 3*q - e1 - e3 )
[ e1, e2, e3 ];
// the i-th normalized eigenvector of a 3x3 symmetrical matrix M from its eigenvalues
// 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])<EPSILON ? I[k] : A[k]/norm(A[k]);
// finds the eigenvector corresponding to the smallest eigenvalue of the covariance matrix of a pointlist
// returns the mean of the points, the eigenvector and the greatest eigenvalue
function _covariance_evec_eval(points) =
let( pm = sum(points)/len(points), // mean point
Y = [ for(i=[0:len(points)-1]) points[i] - pm ],
M = transpose(Y)*Y , // covariance matrix
evals = _eigenvals_symm_3(M), // eigenvalues in decreasing order
evec = _eigenvec_symm_3(M,evals,i=2) )
[pm, evec, evals[0] ];
// Function: plane_from_points()
// Usage:
// plane_from_points(points, <fast>, <eps>);
// 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,12 +991,11 @@ 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()
@ -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,<bounded>,<eps>);
// isect = circle_line_intersection(c,d,line,<bounded>,<eps>);
// isect = circle_line_intersection(c,<r|d>,<line>,<bounded>,<eps>);
// 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,7 +1712,7 @@ 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)<eps
max(distlist)<eps*nrm
? assert(!error, "Cannot find three noncollinear points in pointlist.")
[]
: [0,b,max_index(distlist)];
@ -1727,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.
@ -1738,8 +1783,8 @@ 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) )
plane==undef? undef :
: let( plane = plane_from_polygon(poly) )
plane==[]? [] :
let(
n = plane_normal(plane),
total = sum([
@ -1748,8 +1793,8 @@ function polygon_area(poly, signed=false) =
v1 = poly[i] - poly[0],
v2 = poly[i+1] - poly[0]
)
cross(v1,v2) * n
])/2
cross(v1,v2)
])* n/2
)
signed ? total : abs(total);
@ -1758,27 +1803,32 @@ function polygon_area(poly, signed=false) =
// Usage:
// is_convex_polygon(poly);
// Description:
// Returns true if the given 2D polygon is convex. The result is meaningless if the polygon is not simple (self-intersecting).
// If the points are collinear the result is true.
// Returns true if the given 2D or 3D polygon is convex.
// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar.
// If the points are collinear an error is generated.
// Arguments:
// poly = Polygon to check.
// eps = Tolerance for the collinearity test. Default: EPSILON.
// Example:
// is_convex_polygon(circle(d=50)); // Returns: true
// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
// Example:
// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
// is_convex_polygon(spiral); // Returns: false
function is_convex_polygon(poly) =
assert(is_path(poly,dim=2), "The input should be a 2D polygon." )
let( l = len(poly) )
len([for( i = l-1,
c = cross(poly[(i+1)%l]-poly[i], poly[(i+2)%l]-poly[(i+1)%l]),
s = sign(c);
i>=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) =

18
math.scad
View file

@ -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,<eps>);
// 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,<eps>)

47
polyhedra.scad
View file

@ -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

132
rounding.scad
View file

@ -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) =

47
tests/test_geometry.scad
View file

@ -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]);
@ -715,12 +714,12 @@ module 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();
@ -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();

3
tests/test_vnf.scad
View file

@ -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();

2
version.scad
View file

@ -6,7 +6,7 @@
//////////////////////////////////////////////////////////////////////
BOSL_VERSION = [2,0,601];
BOSL_VERSION = [2,0,605];
// Section: BOSL Library Version Functions

97
vnf.scad
View file

@ -367,6 +367,87 @@ function vnf_vertex_array(
]);
// Function: vnf_tri_array()
// Usage:
// vnf = vnf_tri_array(points, <row_wrap>, <reverse>)
// 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) =