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

Merge pull request #495 from RonaldoCMP/master

Browse source 
Review of geometry.scad for speed and minor other changes
This commit is contained in:
Revar Desmera 2021-04-12 00:33:47 -07:00 committed by GitHub
commit b00547b834
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 264 additions and 219 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
arrays.scad
View file

@ -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 ]; [] == [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 is_homogenous(l, depth=10) = is_homogeneous(l, depth);
function _same_type(a,b, depth) = function _same_type(a,b, depth) =
(depth==0) || (depth==0) ||
@ -50,7 +51,7 @@ function _same_type(a,b, depth) =
(is_string(a) && is_string(b)) || (is_string(a) && is_string(b)) ||
(is_list(a) && is_list(b) && len(a)==len(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] ); && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] );
// Function: select() // Function: select()
// Topics: List Handling // Topics: List Handling
@ -97,7 +98,6 @@ function select(list, start, end) =
// Function: slice() // Function: slice()
// Topics: List Handling
// Usage: // Usage:
// list = slice(list,s,e); // list = slice(list,s,e);
// Description: // Description:
@ -476,7 +476,7 @@ function reverse(x) =
// l9 = list_rotate([1,2,3,4,5],6); // Returns: [2,3,4,5,1] // l9 = list_rotate([1,2,3,4,5],6); // Returns: [2,3,4,5,1]
function list_rotate(list,n=1) = function list_rotate(list,n=1) =
assert(is_list(list)||is_string(list), "Invalid list or string.") 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 ( let (
ll = len(list), ll = len(list),
n = ((n % ll) + ll) % ll, n = ((n % ll) + ll) % ll,
@ -990,7 +990,7 @@ function _sort_vectors(arr, idxlist, _i=0) =
_sort_vectors(equal, idxlist, _i+1), _sort_vectors(equal, idxlist, _i+1),
_sort_vectors(greater, idxlist, _i ) ); _sort_vectors(greater, idxlist, _i ) );
// sorting using compare_vals(); returns indexed list when `indexed==true` // sorting using compare_vals(); returns indexed list when `indexed==true`
function _sort_general(arr, idx=undef, indexed=false) = function _sort_general(arr, idx=undef, indexed=false) =
(len(arr)<=1) ? arr : (len(arr)<=1) ? arr :
@ -1332,6 +1332,8 @@ function permutations(l,n=2) =
// pairs = zip(a,b); // pairs = zip(a,b);
// triples = zip(a,b,c); // triples = zip(a,b,c);
// quads = zip([LIST1,LIST2,LIST3,LIST4]); // quads = zip([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip_long()
// Description: // Description:
// Zips together two or more lists into a single list. For example, if you have two // 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]]. // 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); // pairs = zip_long(a,b);
// triples = zip_long(a,b,c); // triples = zip_long(a,b,c);
// quads = zip_long([LIST1,LIST2,LIST3,LIST4]); // quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip()
// Description: // Description:
// Zips together two or more lists into a single list. For example, if you have two // 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]]. // 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], // [[4,2], 91, false],
// [6, [3,4], undef]]; // [6, [3,4], undef]];
// submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]] // submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]]
function submatrix(M,idx1,idx2) = function submatrix(M,idx1,idx2) =
[for(i=idx1) [for(j=idx2) M[i][j] ] ]; [for(i=idx1) [for(j=idx2) M[i][j] ] ];
@ -1629,7 +1632,6 @@ function block_matrix(M) =
assert(badrows==[], "Inconsistent or invalid input") assert(badrows==[], "Inconsistent or invalid input")
bigM; bigM;
// Function: diagonal_matrix() // Function: diagonal_matrix()
// Usage: // Usage:
// mat = diagonal_matrix(diag, <offdiag>); // mat = diagonal_matrix(diag, <offdiag>);
@ -1855,7 +1857,7 @@ function transpose(arr, reverse=false) =
// A = matrix to test // A = matrix to test
// eps = epsilon for comparing equality. Default: 1e-12 // eps = epsilon for comparing equality. Default: 1e-12
function is_matrix_symmetric(A,eps=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 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

11
common.scad
View file

@ -205,7 +205,8 @@ function is_func(x) = version_num()>20210000 && is_function(x);
// Description: // Description:
// Tests whether input is a list of entries which all have the same list structure // 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 // 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: // Arguments:
// list = list to check // list = list to check
// pattern = optional pattern required to match // 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. // 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. // recursive = If true, sublists are checked recursively for defined values. The first sublist that has a defined item is returned.
// Examples: // 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) = function first_defined(v,recursive=false,_i=0) =
_i<len(v) && ( _i<len(v) && (
is_undef(v[_i]) || ( is_undef(v[_i]) || (
@ -605,15 +606,15 @@ function segs(r) =
// Module: no_children() // Module: no_children()
// Topics: Error Checking
// Usage: // Usage:
// no_children($children); // no_children($children);
// Topics: Error Checking
// See Also: no_function(), no_module()
// Description: // Description:
// Assert that the calling module does not support children. Prints an error message to this effect and fails if children are present, // 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. // as indicated by its argument.
// Arguments: // Arguments:
// $children = number of children the module has. // $children = number of children the module has.
// See Also: no_function(), no_module()
// Example: // Example:
// module foo() { // module foo() {
// no_children($children); // no_children($children);
@ -676,7 +677,7 @@ function _valstr(x) =
// expected = The value that was expected. // expected = The value that was expected.
// info = Extra info to print out to make the error clearer. // info = Extra info to print out to make the error clearer.
// Example: // 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) { module assert_approx(got, expected, info) {
no_children($children); no_children($children);
if (!approx(got, expected)) { if (!approx(got, expected)) {

167
geometry.scad
View file

@ -453,7 +453,7 @@ function segment_closest_point(seg,pt) =
// Usage: // Usage:
// line_from_points(points, [fast], [eps]); // line_from_points(points, [fast], [eps]);
// Description: // 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 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. // if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
// Arguments: // Arguments:
@ -888,9 +888,51 @@ function plane3pt_indexed(points, i1, i2, i3) =
// Example: // Example:
// plane_from_normal([0,0,1], [2,2,2]); // Returns the xy plane passing through the point (2,2,2) // 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]) = function plane_from_normal(normal, pt=[0,0,0]) =
assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),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); 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() // Function: plane_from_points()
@ -898,12 +940,12 @@ function plane_from_normal(normal, pt=[0,0,0]) =
// plane_from_points(points, <fast>, <eps>); // plane_from_points(points, <fast>, <eps>);
// Description: // Description:
// Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane, // 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. // 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 `undef` is returned. // If `fast` is false and the points in the list are collinear or not coplanar, then [] is returned.
// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned. // If `fast` is true, the polygon coplanarity check is skipped and a best fitted plane is returned.
// Arguments: // Arguments:
// points = The list of points to find the plane of. // 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) // eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
// Example(3D): // Example(3D):
// xyzpath = rot(45, v=[-0.3,1,0], p=path3d(star(n=6,id=70,d=100), 70)); // 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) = function plane_from_points(points, fast=false, eps=EPSILON) =
assert( is_path(points,dim=3), "Improper 3d point list." ) assert( is_path(points,dim=3), "Improper 3d point list." )
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
let( len(points) == 3
indices = noncollinear_triple(points,error=false) ? let( plane = plane3pt(points[0],points[1],points[2]) )
) plane==[] ? [] : plane
indices==[] ? undef : : let(
let( covmix = _covariance_evec_eval(points),
p1 = points[indices[0]], pm = covmix[0],
p2 = points[indices[1]], evec = covmix[1],
p3 = points[indices[2]], eval0 = covmix[2],
plane = plane3pt(p1,p2,p3) plane = [ each evec, pm*evec] )
) !fast && _pointlist_greatest_distance(points,plane)>eps*eval0 ? undef :
fast || points_on_plane(points,plane,eps=eps) ? plane : undef; plane ;
// Function: plane_from_polygon() // 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. // 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. // 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 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: // Arguments:
// poly = The planar 3D polygon to find the plane of. // 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 // 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) = function plane_from_polygon(poly, fast=false, eps=EPSILON) =
assert( is_path(poly,dim=3), "Invalid polygon." ) assert( is_path(poly,dim=3), "Invalid polygon." )
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
let( len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) :
poly = deduplicate(poly), let( triple = sort(noncollinear_triple(poly,error=false)) )
n = polygon_normal(poly), triple==[] ? [] :
plane = [n.x, n.y, n.z, n*poly[0]] let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]]))
) fast? plane: points_on_plane(poly, plane, eps=eps)? plane: [];
fast? plane: coplanar(poly,eps=eps)? plane: [];
// Function: plane_normal() // Function: plane_normal()
// Usage: // Usage:
// plane_normal(plane); // plane_normal(plane);
@ -1252,9 +1294,17 @@ function coplanar(points, eps=EPSILON) =
len(points)<=2 ? false len(points)<=2 ? false
: let( ip = noncollinear_triple(points,error=false,eps=eps) ) : let( ip = noncollinear_triple(points,error=false,eps=eps) )
ip == [] ? false : ip == [] ? false :
let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]) )
normal = point3d(plane) ) _pointlist_greatest_distance(points,plane) < eps;
max( points*normal ) - plane[3]< eps*norm(normal);
// 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() // Function: points_on_plane()
@ -1270,9 +1320,7 @@ function points_on_plane(points, plane, eps=EPSILON) =
assert( _valid_plane(plane), "Invalid plane." ) 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_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." ) assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
let( normal = point3d(plane), _pointlist_greatest_distance(points,plane) < eps;
pt_nrm = points*normal )
abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3]))< eps*norm(normal);
// Function: in_front_of_plane() // Function: in_front_of_plane()
@ -1599,20 +1647,19 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
// Function: circle_line_intersection() // Function: circle_line_intersection()
// Usage: // Usage:
// isect = circle_line_intersection(c,r,line,<bounded>,<eps>); // isect = circle_line_intersection(c,<r|d>,<line>,<bounded>,<eps>);
// isect = circle_line_intersection(c,d,line,<bounded>,<eps>);
// Description: // Description:
// Find intersection points between a 2d circle and a line, ray or segment specified by two points. // Find intersection points between a 2d circle and a line, ray or segment specified by two points.
// By default the line is unbounded. // By default the line is unbounded.
// Arguments: // Arguments:
// c = center of circle // c = center of circle
// r = radius of circle // r = radius of circle
// ---
// d = diameter of circle
// line = two points defining the unbounded line // 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 // 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 // eps = epsilon used for identifying the case with one solution. Default: 1e-9
// --- function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) =
// d = diameter of circle
function circle_line_intersection(c,r,line,d,bounded=false,eps=EPSILON) =
let(r=get_radius(r=r,d=d,dflt=undef)) let(r=get_radius(r=r,d=d,dflt=undef))
assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.") assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.")
assert(is_vector(c,2), "Circle center must be a 2-vector") assert(is_vector(c,2), "Circle center must be a 2-vector")
@ -1665,11 +1712,11 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
n = (pb-pa)/nrm, n = (pb-pa)/nrm,
distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)] 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.") ? assert(!error, "Cannot find three noncollinear points in pointlist.")
[] []
: [0,b,max_index(distlist)]; : [0,b,max_index(distlist)];
// Function: pointlist_bounds() // Function: pointlist_bounds()
// Usage: // Usage:
@ -1725,7 +1772,7 @@ function furthest_point(pt, points) =
// area = polygon_area(poly); // area = polygon_area(poly);
// Description: // Description:
// Given a 2D or 3D planar polygon, returns the area of that polygon. // 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. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon.
// Arguments: // Arguments:
// poly = Polygon to compute the area of. // poly = Polygon to compute the area of.
@ -1736,19 +1783,19 @@ function polygon_area(poly, signed=false) =
len(poly[0])==2 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 ) ? 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) signed ? total : abs(total)
: let( plane = plane_from_points(poly) ) : let( plane = plane_from_polygon(poly) )
plane==undef? undef : plane==[]? [] :
let( let(
n = plane_normal(plane), n = plane_normal(plane),
total = sum([ total = sum([
for(i=[1:1:len(poly)-2]) for(i=[1:1:len(poly)-2])
let( let(
v1 = poly[i] - poly[0], v1 = poly[i] - poly[0],
v2 = poly[i+1] - poly[0] v2 = poly[i+1] - poly[0]
) )
cross(v1,v2) * n cross(v1,v2)
])/2 ])* n/2
) )
signed ? total : abs(total); signed ? total : abs(total);
@ -1761,25 +1808,26 @@ function polygon_area(poly, signed=false) =
// If the points are collinear an error is generated. // If the points are collinear an error is generated.
// Arguments: // Arguments:
// poly = Polygon to check. // poly = Polygon to check.
// eps = Tolerance for the collinearity test. Default: EPSILON.
// Example: // Example:
// is_convex_polygon(circle(d=50)); // Returns: true // is_convex_polygon(circle(d=50)); // Returns: true
// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true // is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
// Example: // Example:
// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]]; // spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
// is_convex_polygon(spiral); // Returns: false // is_convex_polygon(spiral); // Returns: false
function is_convex_polygon(poly) = function is_convex_polygon(poly,eps=EPSILON) =
assert(is_path(poly), "The input should be a 2D or 3D polygon." ) assert(is_path(poly), "The input should be a 2D or 3D polygon." )
let( lp = len(poly), let( lp = len(poly),
p0 = poly[0] ) p0 = poly[0] )
assert( lp>=3 , "A polygon must have at least 3 points" ) 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]) ] ) 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 len(p0)==2
? assert( !approx(max(crosses)) && !approx(min(crosses)), "The points are collinear" ) ? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" )
min(crosses) >=0 || max(crosses)<=0 min(crosses) >=0 || max(crosses)<=0
: let( prod = crosses*sum(crosses), : let( prod = crosses*sum(crosses),
minc = min(prod), minc = min(prod),
maxc = max(prod) ) maxc = max(prod) )
assert( !approx(maxc-minc), "The points are collinear" ) assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" )
minc>=0 || maxc<=0; minc>=0 || maxc<=0;
@ -2008,7 +2056,7 @@ function point_in_polygon(point, poly, nonzero=true, eps=EPSILON) =
// poly = The list of 2D path points for the perimeter of the polygon. // poly = The list of 2D path points for the perimeter of the polygon.
function polygon_is_clockwise(poly) = function polygon_is_clockwise(poly) =
assert(is_path(poly,dim=2), "Input should be a 2d path") 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() // Function: clockwise_polygon()
@ -2052,19 +2100,16 @@ function reverse_polygon(poly) =
// n = polygon_normal(poly); // n = polygon_normal(poly);
// Description: // Description:
// Given a 3D planar polygon, returns a unit-length normal vector for the // 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: // Arguments:
// poly = The list of 3D path points for the perimeter of the polygon. // poly = The list of 3D path points for the perimeter of the polygon.
function polygon_normal(poly) = function polygon_normal(poly) =
assert(is_path(poly,dim=3), "Invalid 3D polygon." ) assert(is_path(poly,dim=3), "Invalid 3D polygon." )
let( len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) :
poly = cleanup_path(poly), let( triple = sort(noncollinear_triple(poly,error=false)) )
p0 = poly[0], triple==[] ? [] :
n = sum([ point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ;
for (i=[1:1:len(poly)-2])
cross(poly[i+1]-p0, poly[i]-p0)
])
) unit(n,undef);
function _split_polygon_at_x(poly, x) = function _split_polygon_at_x(poly, x) =

4
math.scad
View file

@ -985,10 +985,8 @@ function determinant(M) =
function is_matrix(A,m,n,square=false) = function is_matrix(A,m,n,square=false) =
is_list(A) is_list(A)
&& (( is_undef(m) && len(A) ) || len(A)==m) && (( 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])) && (!square || len(A) == len(A[0]))
&& is_vector(A[0]) && is_vector(A[0],n)
&& is_consistent(A); && is_consistent(A);

216
paths.scad
View file

@ -454,120 +454,120 @@ function path_torsion(path, closed=false) =
// path2 = path_chamfer_and_rounding(path, closed=true, chamfer=chamfs, rounding=rs); // path2 = path_chamfer_and_rounding(path, closed=true, chamfer=chamfs, rounding=rs);
// stroke(path2, closed=true); // stroke(path2, closed=true);
function path_chamfer_and_rounding(path, closed, chamfer, rounding) = function path_chamfer_and_rounding(path, closed, chamfer, rounding) =
let ( let (
path = deduplicate(path,closed=true), path = deduplicate(path,closed=true),
lp = len(path), lp = len(path),
chamfer = is_undef(chamfer)? repeat(0,lp) : chamfer = is_undef(chamfer)? repeat(0,lp) :
is_vector(chamfer)? list_pad(chamfer,lp,0) : is_vector(chamfer)? list_pad(chamfer,lp,0) :
is_num(chamfer)? repeat(chamfer,lp) : is_num(chamfer)? repeat(chamfer,lp) :
assert(false, "Bad chamfer value."), assert(false, "Bad chamfer value."),
rounding = is_undef(rounding)? repeat(0,lp) : rounding = is_undef(rounding)? repeat(0,lp) :
is_vector(rounding)? list_pad(rounding,lp,0) : is_vector(rounding)? list_pad(rounding,lp,0) :
is_num(rounding)? repeat(rounding,lp) : is_num(rounding)? repeat(rounding,lp) :
assert(false, "Bad rounding value."), assert(false, "Bad rounding value."),
corner_paths = [ corner_paths = [
for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let( for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let(
p1 = select(path,i-1), p1 = select(path,i-1),
p2 = select(path,i), p2 = select(path,i),
p3 = select(path,i+1) p3 = select(path,i+1)
) )
chamfer[i] > 0? _corner_chamfer_path(p1, p2, p3, side=chamfer[i]) : chamfer[i] > 0? _corner_chamfer_path(p1, p2, p3, side=chamfer[i]) :
rounding[i] > 0? _corner_roundover_path(p1, p2, p3, r=rounding[i]) : rounding[i] > 0? _corner_roundover_path(p1, p2, p3, r=rounding[i]) :
[p2] [p2]
], ],
out = [ out = [
if (!closed) path[0], if (!closed) path[0],
for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let( for (i=(closed? [0:1:lp-1] : [1:1:lp-2])) let(
p1 = select(path,i-1), p1 = select(path,i-1),
p2 = select(path,i), p2 = select(path,i),
crn1 = select(corner_paths,i-1), crn1 = select(corner_paths,i-1),
crn2 = corner_paths[i], crn2 = corner_paths[i],
l1 = norm(last(crn1)-p1), l1 = norm(last(crn1)-p1),
l2 = norm(crn2[0]-p2), l2 = norm(crn2[0]-p2),
needed = l1 + l2, needed = l1 + l2,
seglen = norm(p2-p1), seglen = norm(p2-p1),
check = assert(seglen >= needed, str("Path segment ",i," is too short to fulfill rounding/chamfering for the adjacent corners.")) check = assert(seglen >= needed, str("Path segment ",i," is too short to fulfill rounding/chamfering for the adjacent corners."))
) each crn2, ) each crn2,
if (!closed) last(path) if (!closed) last(path)
] ]
) deduplicate(out); ) deduplicate(out);
function _corner_chamfer_path(p1, p2, p3, dist1, dist2, side, angle) = function _corner_chamfer_path(p1, p2, p3, dist1, dist2, side, angle) =
let( let(
v1 = unit(p1 - p2), v1 = unit(p1 - p2),
v2 = unit(p3 - p2), v2 = unit(p3 - p2),
n = vector_axis(v1,v2), n = vector_axis(v1,v2),
ang = vector_angle(v1,v2), ang = vector_angle(v1,v2),
path = (is_num(dist1) && is_undef(dist2) && is_undef(side))? ( path = (is_num(dist1) && is_undef(dist2) && is_undef(side))? (
// dist1 & optional angle // dist1 & optional angle
assert(dist1 > 0) assert(dist1 > 0)
let(angle = default(angle,(180-ang)/2)) let(angle = default(angle,(180-ang)/2))
assert(is_num(angle)) assert(is_num(angle))
assert(angle > 0 && angle < 180) assert(angle > 0 && angle < 180)
let( let(
pta = p2 + dist1*v1, pta = p2 + dist1*v1,
a3 = 180 - angle - ang a3 = 180 - angle - ang
) assert(a3>0, "Angle too extreme.") ) assert(a3>0, "Angle too extreme.")
let( let(
side = sin(angle) * dist1/sin(a3), side = sin(angle) * dist1/sin(a3),
ptb = p2 + side*v2 ptb = p2 + side*v2
) [pta, ptb] ) [pta, ptb]
) : (is_undef(dist1) && is_num(dist2) && is_undef(side))? ( ) : (is_undef(dist1) && is_num(dist2) && is_undef(side))? (
// dist2 & optional angle // dist2 & optional angle
assert(dist2 > 0) assert(dist2 > 0)
let(angle = default(angle,(180-ang)/2)) let(angle = default(angle,(180-ang)/2))
assert(is_num(angle)) assert(is_num(angle))
assert(angle > 0 && angle < 180) assert(angle > 0 && angle < 180)
let( let(
ptb = p2 + dist2*v2, ptb = p2 + dist2*v2,
a3 = 180 - angle - ang a3 = 180 - angle - ang
) assert(a3>0, "Angle too extreme.") ) assert(a3>0, "Angle too extreme.")
let( let(
side = sin(angle) * dist2/sin(a3), side = sin(angle) * dist2/sin(a3),
pta = p2 + side*v1 pta = p2 + side*v1
) [pta, ptb] ) [pta, ptb]
) : (is_undef(dist1) && is_undef(dist2) && is_num(side))? ( ) : (is_undef(dist1) && is_undef(dist2) && is_num(side))? (
// side & optional angle // side & optional angle
assert(side > 0) assert(side > 0)
let(angle = default(angle,(180-ang)/2)) let(angle = default(angle,(180-ang)/2))
assert(is_num(angle)) assert(is_num(angle))
assert(angle > 0 && angle < 180) assert(angle > 0 && angle < 180)
let( let(
a3 = 180 - angle - ang a3 = 180 - angle - ang
) assert(a3>0, "Angle too extreme.") ) assert(a3>0, "Angle too extreme.")
let( let(
dist1 = sin(a3) * side/sin(ang), dist1 = sin(a3) * side/sin(ang),
dist2 = sin(angle) * side/sin(ang), dist2 = sin(angle) * side/sin(ang),
pta = p2 + dist1*v1, pta = p2 + dist1*v1,
ptb = p2 + dist2*v2 ptb = p2 + dist2*v2
) [pta, ptb] ) [pta, ptb]
) : (is_num(dist1) && is_num(dist2) && is_undef(side) && is_undef(side))? ( ) : (is_num(dist1) && is_num(dist2) && is_undef(side) && is_undef(side))? (
// dist1 & dist2 // dist1 & dist2
assert(dist1 > 0) assert(dist1 > 0)
assert(dist2 > 0) assert(dist2 > 0)
let( let(
pta = p2 + dist1*v1, pta = p2 + dist1*v1,
ptb = p2 + dist2*v2 ptb = p2 + dist2*v2
) [pta, ptb] ) [pta, ptb]
) : ( ) : (
assert(false,"Bad arguments.") assert(false,"Bad arguments.")
) )
) path; ) path;
function _corner_roundover_path(p1, p2, p3, r, d) = function _corner_roundover_path(p1, p2, p3, r, d) =
let( let(
r = get_radius(r=r,d=d,dflt=undef), r = get_radius(r=r,d=d,dflt=undef),
res = circle_2tangents(p1, p2, p3, r=r, tangents=true), res = circle_2tangents(p1, p2, p3, r=r, tangents=true),
cp = res[0], cp = res[0],
n = res[1], n = res[1],
tp1 = res[2], tp1 = res[2],
ang = res[4]+res[5], ang = res[4]+res[5],
steps = floor(segs(r)*ang/360+0.5), steps = floor(segs(r)*ang/360+0.5),
step = ang / steps, step = ang / steps,
path = [for (i=[0:1:steps]) move(cp, p=rot(a=-i*step, v=n, p=tp1-cp))] path = [for (i=[0:1:steps]) move(cp, p=rot(a=-i*step, v=n, p=tp1-cp))]
) path; ) path;

50
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] ]; 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] ); poly2 = concat(poly1, [sum(poly1)/3] );
info = info_str([["poly1 = ",poly1],["poly2 = ",poly2]]); info = info_str([["poly1 = ",poly1],["poly2 = ",poly2]]);
assert_std(plane_from_polygon(poly1),plane3pt(poly1[0],poly1[1],poly1[2]),info); assert_approx(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(poly2),plane3pt(poly1[0],poly1[1],poly1[2]),info);
} }
*test_plane_from_polygon(); *test_plane_from_polygon();
@ -208,8 +208,7 @@ module test_plane_from_normal(){
displ = normal*point; displ = normal*point;
info = info_str([["normal = ",normal],["point = ",point],["displ = ",displ]]); info = info_str([["normal = ",normal],["point = ",point],["displ = ",displ]]);
assert_approx(plane_from_normal(normal,point)*[each point,-1],0,info); 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_approx(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]);
assert_std(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]);
} }
*test_plane_from_normal(); *test_plane_from_normal();
@ -680,23 +679,23 @@ module test_triangle_area() {
module test_plane3pt() { module test_plane3pt() {
assert_std(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]); assert_approx(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_approx(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_approx(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_approx(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_approx(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,2], [10,10,2], [20,0,2]), [0,0,1,2]);
} }
*test_plane3pt(); *test_plane3pt();
module test_plane3pt_indexed() { module test_plane3pt_indexed() {
pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ]; pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ];
s13 = sqrt(1/3); s13 = sqrt(1/3);
assert_std(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]); assert_approx(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]);
assert_std(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]); assert_approx(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]);
assert_std(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]); assert_approx(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]);
assert_std(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]); assert_approx(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,2,1), [0,0,1,0]);
assert_approx(plane3pt_indexed(pts, 0,1,2), [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, 3,2,1), [s13,s13,s13,10*s13]);
assert_approx(plane3pt_indexed(pts, 1,2,3), [-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,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,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,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(); *test_plane_from_points();
module test_plane_normal() { module test_plane_normal() {
assert_std(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]); assert_approx(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_approx(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_approx(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_approx(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_approx(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,2], [10,10,2], [20,0,2])), [0,0,1]);
} }
*test_plane_normal(); *test_plane_normal();
@ -780,7 +779,7 @@ module test_coplanar() {
assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true); 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,1], [0,1,2] ]) == false);
assert(coplanar([ [0,0,0],[1,0,1],[1,1,2], [0,1,1] ]) == true); assert(coplanar([ [0,0,0],[1,0,1],[1,1,2], [0,1,1] ]) == true);
} }
*test_coplanar(); *test_coplanar();
@ -836,7 +835,9 @@ module test_cleanup_path() {
module test_polygon_area() { module test_polygon_area() {
assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4)); 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(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(); *test_polygon_area();
@ -846,6 +847,7 @@ module test_is_convex_polygon() {
assert(is_convex_polygon(circle(r=50,$fn=1000))); assert(is_convex_polygon(circle(r=50,$fn=1000)));
assert(is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50))))); 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([[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(); *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(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=CENTER, $fn=36)), [0,0,0]);
assert_approx(vnf_centroid(sphere(d=100, anchor=BOT, $fn=36)), [0,0,50]); 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(); test_vnf_centroid();

16
vnf.scad
View file

@ -450,18 +450,13 @@ function vnf_volume(vnf) =
// Returns the centroid of the given manifold VNF. The VNF must describe a valid polyhedron with consistent face direction and // 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. // 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. // The centroid of the total is the volume weighted average.
function vnf_centroid(vnf) = function vnf_centroid(vnf) =
assert(is_vnf(vnf) && len(vnf[0])!=0 )
let( let(
verts = vnf[0], 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([ pos = sum([
for(face=vnf[1], j=[1:1:len(face)-2]) let( for(face=vnf[1], j=[1:1:len(face)-2]) let(
v0 = verts[face[0]], v0 = verts[face[0]],
@ -469,10 +464,11 @@ function vnf_centroid(vnf) =
v2 = verts[face[j+1]], v2 = verts[face[j+1]],
vol = cross(v2,v1)*v0 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) = function _triangulate_planar_convex_polygons(polys) =