Merge branch 'master' into master

This commit is contained in:
Revar Desmera 2021-11-06 15:19:01 -07:00 committed by GitHub
commit 41b7c210e8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 367 additions and 148 deletions

View file

@ -12,7 +12,8 @@
// Usage: // Usage:
// test = approx(a, b, [eps]) // test = approx(a, b, [eps])
// Description: // Description:
// Compares two numbers or vectors, and returns true if they are closer than `eps` to each other. // Compares two numbers, vectors, or matrices. Returns true if they are closer than `eps` to each other.
// Results are undefined if `a` and `b` are of different types, or if vectors or matrices contain non-numbers.
// Arguments: // Arguments:
// a = First value. // a = First value.
// b = Second value. // b = Second value.
@ -21,12 +22,22 @@
// test1 = approx(-0.3333333333,-1/3); // Returns: true // test1 = approx(-0.3333333333,-1/3); // Returns: true
// test2 = approx(0.3333333333,1/3); // Returns: true // test2 = approx(0.3333333333,1/3); // Returns: true
// test3 = approx(0.3333,1/3); // Returns: false // test3 = approx(0.3333,1/3); // Returns: false
// test4 = approx(0.3333,1/3,eps=1e-3); // Returns: true // test4 = approx(0.3333,1/3,eps=1e-3); // Returns: true
// test5 = approx(PI,3.1415926536); // Returns: true // test5 = approx(PI,3.1415926536); // Returns: true
// test6 = approx([0,0,sin(45)],[0,0,sqrt(2)/2]); // Returns: true
function approx(a,b,eps=EPSILON) = function approx(a,b,eps=EPSILON) =
(a==b && is_bool(a) == is_bool(b)) || a == b? is_bool(a) == is_bool(b) :
(is_num(a) && is_num(b) && abs(a-b) <= eps) || is_num(a) && is_num(b)? abs(a-b) <= eps :
(is_list(a) && is_list(b) && len(a) == len(b) && [] == [for (i=idx(a)) if (!approx(a[i],b[i],eps=eps)) 1]); is_list(a) && is_list(b) && len(a) == len(b)? (
[] == [
for (i=idx(a))
let(aa=a[i], bb=b[i])
if(
is_num(aa) && is_num(bb)? abs(aa-bb) > eps :
!approx(aa,bb,eps=eps)
) 1
]
) : false;
// Function: all_zero() // Function: all_zero()
@ -290,12 +301,12 @@ function compare_lists(a, b) =
// idx = find_approx(val, list, [start=], [eps=]); // idx = find_approx(val, list, [start=], [eps=]);
// indices = find_approx(val, list, all=true, [start=], [eps=]); // indices = find_approx(val, list, all=true, [start=], [eps=]);
// Description: // Description:
// Finds the first item in `list` that matches `val`, returning the index. // Finds the first item in `list` that matches `val`, returning the index. Returns `undef` if there is no match.
// Arguments: // Arguments:
// val = The value to search for. If given a function literal of signature `function (x)`, uses that function to check list items. Returns true for a match. // val = The value to search for. If given a function literal of signature `function (x)`, uses that function to check list items. Returns true for a match.
// list = The list to search through. // list = The list to search through.
// --- // ---
// start = The index to start searching from. // start = The index to start searching from. Default: 0
// all = If true, returns a list of all matching item indices. // all = If true, returns a list of all matching item indices.
// eps = The maximum allowed floating point rounding error for numeric comparisons. // eps = The maximum allowed floating point rounding error for numeric comparisons.
function find_approx(val, list, start=0, all=false, eps=EPSILON) = function find_approx(val, list, start=0, all=false, eps=EPSILON) =
@ -778,3 +789,4 @@ function list_smallest(list, k) =
let( bigger = [for(li=list) if(li>v) li ] ) let( bigger = [for(li=list) if(li>v) li ] )
concat(smaller, equal, list_smallest(bigger, k-len(smaller) -len(equal))); concat(smaller, equal, list_smallest(bigger, k-len(smaller) -len(equal)));
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

View file

@ -38,14 +38,10 @@ function _is_point_on_line(point, line, bounded=false, eps=EPSILON) =
t = v0*v1/(v1*v1), t = v0*v1/(v1*v1),
bounded = force_list(bounded,2) bounded = force_list(bounded,2)
) )
abs(cross(v0,v1))<eps*norm(v1) abs(cross(v0,v1))<=eps*norm(v1)
&& (!bounded[0] || t>=-eps) && (!bounded[0] || t>=-eps)
&& (!bounded[1] || t<1+eps) ; && (!bounded[1] || t<1+eps) ;
function xis_point_on_line(point, line, bounded=false, eps=EPSILON) =
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
point_line_distance(point, line, bounded)<eps;
///Internal - distance from point `d` to the line passing through the origin with unit direction n ///Internal - distance from point `d` to the line passing through the origin with unit direction n
///_dist2line works for any dimension ///_dist2line works for any dimension
@ -61,7 +57,68 @@ function _valid_line(line,dim,eps=EPSILON) =
function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
/// Internal Function: point_left_of_line2d() /// Internal Function: _is_at_left()
/// Usage:
/// pt = point_left_of_line2d(point, line);
/// Topics: Geometry, Points, Lines
/// Description:
/// Return true iff a 2d point is on or at left of the line defined by `line`.
/// Arguments:
/// pt = The 2d point to check position of.
/// line = Array of two 2d points forming the line segment to test against.
/// eps = Tolerance in the geometrical tests.
function _is_at_left(pt,line,eps=EPSILON) = _tri_class([pt,line[0],line[1]],eps) <= 0;
/// Internal Function: _degenerate_tri()
/// Usage:
/// degen = _degenerate_tri(triangle);
/// Topics: Geometry, Triangles
/// Description:
/// Return true for a specific kind of degeneracy: any two triangle vertices are equal
/// Arguments:
/// tri = A list of three 2d points
/// eps = Tolerance in the geometrical tests.
function _degenerate_tri(tri,eps) =
max(norm(tri[0]-tri[1]), norm(tri[1]-tri[2]), norm(tri[2]-tri[0])) < eps ;
/// Internal Function: _tri_class()
/// Usage:
/// class = _tri_class(triangle);
/// Topics: Geometry, Triangles
/// Description:
/// Return 1 if the triangle `tri` is CW.
/// Return 0 if the triangle `tri` has colinear vertices.
/// Return -1 if the triangle `tri` is CCW.
/// Arguments:
/// tri = A list of the three 2d vertices of a triangle.
/// eps = Tolerance in the geometrical tests.
function _tri_class(tri, eps=EPSILON) =
let( crx = cross(tri[1]-tri[2],tri[0]-tri[2]) )
abs( crx ) <= eps*norm(tri[1]-tri[2])*norm(tri[0]-tri[2]) ? 0 : sign( crx );
/// Internal Function: _pt_in_tri()
/// Usage:
/// class = _pt_in_tri(point, tri);
/// Topics: Geometry, Points, Triangles
/// Description:
// For CW triangles `tri` :
/// return 1 if point is inside the triangle interior.
/// return =0 if point is on the triangle border.
/// return -1 if point is outside the triangle.
/// Arguments:
/// point = The point to check position of.
/// tri = A list of the three 2d vertices of a triangle.
/// eps = Tolerance in the geometrical tests.
function _pt_in_tri(point, tri, eps=EPSILON) =
min( _tri_class([tri[0],tri[1],point],eps),
_tri_class([tri[1],tri[2],point],eps),
_tri_class([tri[2],tri[0],point],eps) );
/// Internal Function: _point_left_of_line2d()
/// Usage: /// Usage:
/// pt = point_left_of_line2d(point, line); /// pt = point_left_of_line2d(point, line);
/// Topics: Geometry, Points, Lines /// Topics: Geometry, Points, Lines
@ -72,10 +129,11 @@ function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps
/// Arguments: /// Arguments:
/// point = The point to check position of. /// point = The point to check position of.
/// line = Array of two points forming the line segment to test against. /// line = Array of two points forming the line segment to test against.
function _point_left_of_line2d(point, line) = function _point_left_of_line2d(point, line, eps=EPSILON) =
assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." ) assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." )
cross(line[0]-point, line[1]-line[0]); // cross(line[0]-point, line[1]-line[0]);
_tri_class([point,line[1],line[0]],eps);
// Function: is_collinear() // Function: is_collinear()
// Usage: // Usage:
@ -1640,11 +1698,11 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
// Function: polygon_triangulate() // Function: polygon_triangulate()
// Usage: // Usage:
// triangles = polygon_triangulate(poly, [ind], [eps]) // triangles = polygon_triangulate(poly, [ind], [error], [eps])
// Description: // Description:
// Given a simple polygon in 2D or 3D, triangulates it and returns a list // Given a simple polygon in 2D or 3D, triangulates it and returns a list
// of triples indexing into the polygon vertices. When the optional argument `ind` is // of triples indexing into the polygon vertices. When the optional argument `ind` is
// given, it is used as an index list into `poly` to define the polygon. In that case, // given, it is used as an index list into `poly` to define the polygon vertices. In that case,
// `poly` may have a length greater than `ind`. When `ind` is undefined, all points in `poly` // `poly` may have a length greater than `ind`. When `ind` is undefined, all points in `poly`
// are considered as vertices of the polygon. // are considered as vertices of the polygon.
// . // .
@ -1653,46 +1711,50 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
// vector with the same direction of the polygon normal. // vector with the same direction of the polygon normal.
// . // .
// The function produce correct triangulations for some non-twisted non-simple polygons. // The function produce correct triangulations for some non-twisted non-simple polygons.
// A polygon is non-twisted iff it is simple or there is a partition of it in // A polygon is non-twisted iff it is simple or it has a partition in
// simple polygons with the same winding. These polygons may have "touching" vertices // simple polygons with the same winding such that the intersection of any two partitions is
// made of full edges and/or vertices of both partitions. These polygons may have "touching" vertices
// (two vertices having the same coordinates, but distinct adjacencies) and "contact" edges // (two vertices having the same coordinates, but distinct adjacencies) and "contact" edges
// (edges whose vertex pairs have the same pairwise coordinates but are in reversed order) but has // (edges whose vertex pairs have the same pairwise coordinates but are in reversed order) but has
// no self-crossing. See examples bellow. If all polygon edges are contact edges, // no self-crossing. See examples bellow. If all polygon edges are contact edges (polygons with
// it returns an empty list for 2d polygons and issues an error for 3d polygons. // zero area), it returns an empty list for 2d polygons and reports an error for 3d polygons.
// Triangulation errors are reported either by an assert error (when `error=true`) or by returning
// `undef` (when `error=false`). Invalid arguments always produce an assert error.
// . // .
// Self-crossing polygons have no consistent winding and usually produce an error but // Twisted polygons have no consistent winding and when input to this function usually reports
// when an error is not issued the outputs are not correct triangulations. The function // an error but when an error is not reported the outputs are not correct triangulations. The function
// can work for 3d non-planar polygons if they are close enough to planar but may otherwise // can work for 3d non-planar polygons if they are close enough to planar but may otherwise
// issue an error for this case. // report an error for this case.
// Arguments: // Arguments:
// poly = Array of vertices for the polygon. // poly = Array of the polygon vertices.
// ind = A list indexing the vertices of the polygon in `poly`. // ind = A list indexing the vertices of the polygon in `poly`.
// error = If false, returns `undef` when the polygon cannot be triangulated; otherwise, issues an assert error. Default: true.
// eps = A maximum tolerance in geometrical tests. Default: EPSILON // eps = A maximum tolerance in geometrical tests. Default: EPSILON
// Example(2D,NoAxes): // Example(2D,NoAxes): a simple polygon; see from above
// poly = star(id=10, od=15,n=11); // poly = star(id=10, od=15,n=11);
// tris = polygon_triangulate(poly); // tris = polygon_triangulate(poly);
// color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("lightblue") for(tri=tris) polygon(select(poly,tri));
// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
// color("magenta") up(2) stroke(poly,.25,closed=true); // color("magenta") up(2) stroke(poly,.25,closed=true);
// color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
// Example(2D,NoAxes): a polygon with a hole and one "contact" edge // Example(2D,NoAxes): a polygon with a hole and one "contact" edge; see from above
// poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ]; // poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ];
// tris = polygon_triangulate(poly); // tris = polygon_triangulate(poly);
// color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("lightblue") for(tri=tris) polygon(select(poly,tri));
// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
// color("magenta") up(2) stroke(poly,.25,closed=true); // color("magenta") up(2) stroke(poly,.25,closed=true);
// color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
// Example(2D,NoAxes): a polygon with "touching" vertices and no holes // Example(2D,NoAxes): a polygon with "touching" vertices and no holes; see from above
// poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ]; // poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ];
// tris = polygon_triangulate(poly); // tris = polygon_triangulate(poly);
// color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("lightblue") for(tri=tris) polygon(select(poly,tri));
// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
// color("magenta") up(2) stroke(poly,.25,closed=true); // color("magenta") up(2) stroke(poly,.25,closed=true);
// color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
// Example(2D,NoAxes): a polygon with "contact" edges and no holes // Example(2D,NoAxes): a polygon with "contact" edges and no holes; see from above
// poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3], // poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3],
// [7,7], [7,3], [3,3] ]; // [7,7], [7,3], [3,3] ];
// tris = polygon_triangulate(poly); // see from the top // tris = polygon_triangulate(poly);
// color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("lightblue") for(tri=tris) polygon(select(poly,tri));
// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
// color("magenta") up(2) stroke(poly,.25,closed=true); // color("magenta") up(2) stroke(poly,.25,closed=true);
@ -1704,102 +1766,122 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
// vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ]; // vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ];
// color("blue") // color("blue")
// vnf_wireframe(vnf_tri, width=.15); // vnf_wireframe(vnf_tri, width=.15);
function polygon_triangulate(poly, ind, eps=EPSILON) = function polygon_triangulate(poly, ind, error=true, eps=EPSILON) =
assert(is_path(poly) && len(poly)>=3, "Polygon `poly` should be a list of at least three 2d or 3d points") assert(is_path(poly) && len(poly)>=3, "Polygon `poly` should be a list of at least three 2d or 3d points")
assert(is_undef(ind) assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
|| (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
"Improper or out of bounds list of indices") "Improper or out of bounds list of indices")
let( ind = is_undef(ind) ? count(len(poly)) : ind ) let( ind = is_undef(ind) ? count(len(poly)) : ind )
len(ind) <=2 ? [] :
len(ind) == 3 len(ind) == 3
? _is_degenerate([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] :
// non zero area // non zero area
assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 2*eps, let( degen = norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) < 2*eps )
"The polygon vertices are collinear.") assert( ! error || ! degen, "The polygon vertices are collinear.")
[ind] degen ? undef : [ind]
: len(poly[ind[0]]) == 3 : len(poly[ind[0]]) == 3
? // represents the polygon projection on its plane as a 2d polygon ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane
let( let(
ind = deduplicate_indexed(poly, ind, eps) ind = deduplicate_indexed(poly, ind, eps)
) )
len(ind)<3 ? [] : len(ind)<3 ? [] :
let( let(
pts = select(poly,ind), pts = select(poly,ind),
nrm = polygon_normal(pts) nrm = -polygon_normal(pts)
) )
assert( nrm!=undef, assert( ! error || (nrm != undef),
"The polygon has self-intersections or its vertices are collinear or non coplanar.") "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.")
nrm == undef ? undef :
let( let(
imax = max_index([for(p=pts) norm(p-pts[0]) ]), imax = max_index([for(p=pts) norm(p-pts[0]) ]),
v1 = unit( pts[imax] - pts[0] ), v1 = unit( pts[imax] - pts[0] ),
v2 = cross(v1,nrm), v2 = cross(v1,nrm),
prpts = pts*transpose([v1,v2]) prpts = pts*transpose([v1,v2]) // the 2d projection of pts on the polygon plane
) )
[for(tri=_triangulate(prpts, count(len(ind)), eps)) select(ind,tri) ] let( tris = _triangulate(prpts, count(len(ind)), error, eps) )
: let( cw = is_polygon_clockwise(select(poly, ind)) ) tris == undef ? undef :
cw [for(tri=tris) select(ind,tri) ]
? [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ] : is_polygon_clockwise(select(poly, ind))
: _triangulate( poly, ind, eps ); ? _triangulate( poly, ind, error, eps )
: let( tris = _triangulate( poly, reverse(ind), error, eps ) )
tris == undef ? undef :
[for(tri=tris) reverse(tri) ];
function _triangulate(poly, ind, eps=EPSILON, tris=[]) = // poly is supposed to be a 2d cw polygon
// implements a modified version of ear cut method for non-twisted polygons
// the polygons accepted by this function are those decomposable in simple
// CW polygons.
function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) =
len(ind)==3 len(ind)==3
? _is_degenerate(select(poly,ind),eps) ? _degenerate_tri(select(poly,ind),eps)
? tris // last 3 pts perform a degenerate triangle, ignore it ? tris // if last 3 pts perform a degenerate triangle, ignore it
: concat(tris,[ind]) // otherwise, include it : concat(tris,[ind]) // otherwise, include it
: let( ear = _get_ear(poly,ind,eps) ) : let( ear = _get_ear(poly,ind,eps) )
assert( ear!=undef, assert( ! error || (ear != undef),
"The polygon has self-intersections or its vertices are collinear or non coplanar.") "The polygon has twists or all its vertices are collinear or non coplanar.")
is_list(ear) // degenerate ear ear == undef ? undef :
? _triangulate(poly, select(ind,ear[0]+2, ear[0]), eps, tris) // discard it is_list(ear) // is it a degenerate ear ?
? len(ind) <= 4 ? tris :
_triangulate(poly, select(ind,ear[0]+3, ear[0]), error, eps, tris) // discard it
: let( : let(
ear_tri = select(ind,ear,ear+2), ear_tri = select(ind,ear,ear+2),
indr = select(ind,ear+2, ear) // remaining point indices indr = select(ind,ear+2, ear) // indices of the remaining path
) )
_triangulate(poly, indr, eps, concat(tris,[ear_tri])); _triangulate(poly, indr, error, eps, concat(tris,[ear_tri]));
// a returned ear will be: // a returned ear will be:
// 1. a CCW (non-degenerate) triangle, made of subsequent vertices, without other // 1. a CW non-reflex triangle, made of subsequent poly vertices, without any other
// points inside except possibly at its vertices // poly points inside except possibly at its own vertices
// 2. or a degenerate triangle where two vertices are coincident // 2. or a degenerate triangle where two vertices are coincident
// the returned ear is specified by the index of `ind` of its first vertex // the returned ear is specified by the index of `ind` of its first vertex
function _get_ear(poly, ind, eps, _i=0) = function _get_ear(poly, ind, eps, _i=0) =
_i>=len(ind) ? undef : // poly has no ears let( lind = len(ind) )
lind==3 ? 0 :
let( // the _i-th ear candidate let( // the _i-th ear candidate
p0 = poly[ind[_i]], p0 = poly[ind[_i]],
p1 = poly[ind[(_i+1)%len(ind)]], p1 = poly[ind[(_i+1)%lind]],
p2 = poly[ind[(_i+2)%len(ind)]] p2 = poly[ind[(_i+2)%lind]]
) )
// degenerate triangles are returned codified // if vertex p1 is a convex candidate to be an ear,
_is_degenerate([p0,p1,p2],eps) ? [_i] : // check if the triangle [p0,p1,p2] contains any other point
// if it is not a convex vertex, check the next one // except possibly p0 and p2
_is_cw2(p0,p1,p2,eps) ? _get_ear(poly,ind,eps, _i=_i+1) : // exclude the ear candidate central vertex p1 from the verts to check
let( // vertex p1 is convex _tri_class([p0,p1,p2],eps) > 0
// check if the triangle contains any other point && _none_inside(select(ind,_i+2, _i),poly,p0,p1,p2,eps) ? _i : // found an ear
// except possibly its own vertices // otherwise check the next ear candidate
to_tst = select(ind,_i+3, _i-1), _i<lind-1 ? _get_ear(poly, ind, eps, _i=_i+1) :
q = [(p0-p2).y, (p2-p0).x], // orthogonal to ray [p0,p2] pointing right // poly has no ears, look for wiskers
r = [(p2-p1).y, (p1-p2).x], // orthogonal to ray [p2,p1] pointing right let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
s = [(p1-p0).y, (p0-p1).x], // orthogonal to ray [p1,p0] pointing right wiskers==[] ? undef : [wiskers[0]];
inside = [for(p=select(poly,to_tst)) // for vertices other than p0, p1 and p2
if( (p-p0)*q<=0 && (p-p2)*r<=0 && (p-p1)*s<=0 // p is on the triangle
&& norm(p-p0)>eps // but not on any vertex of it
&& norm(p-p1)>eps // returns false ASA it finds some reflex vertex of poly[idxs[.]]
&& norm(p-p2)>eps ) // inside the triangle different from p0 and p2
p ] // note: to simplify the expressions it is assumed that the input polygon has no twists
function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
i>=len(idxs) ? true :
let(
vert = poly[idxs[i]],
prev_vert = poly[select(idxs,i-1)],
next_vert = poly[select(idxs,i+1)]
) )
inside==[] ? _i : // found an ear // check if vert prevent [p0,p1,p2] to be an ear
// check the next ear candidate // this conditions might have a simpler expression
_get_ear(poly, ind, eps, _i=_i+1); _tri_class([prev_vert, vert, next_vert],eps) <= 0 // reflex condition
&& ( // vert is a cw reflex poly vertex inside the triangle [p0,p1,p2]
( _tri_class([p0,p1,vert],eps)>0 &&
// true for some specific kinds of degeneracy _tri_class([p1,p2,vert],eps)>0 &&
function _is_degenerate(tri,eps) = _tri_class([p2,p0,vert],eps)>=0 )
norm(tri[0]-tri[1])<eps || norm(tri[1]-tri[2])<eps || norm(tri[2]-tri[0])<eps ; // or it is equal to p1 and some of its adjacent edges cross the open segment (p0,p2)
|| ( norm(vert-p1) < eps
&& _is_at_left(p0,[prev_vert,p1],eps) && _is_at_left(p2,[p1,prev_vert],eps)
function _is_cw2(a,b,c,eps=EPSILON) = cross(a-c,b-c)<eps*norm(a-c)*norm(b-c); && _is_at_left(p2,[p1,next_vert],eps) && _is_at_left(p0,[next_vert,p1],eps)
)
)
? false
: _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1);
// Function: is_polygon_clockwise() // Function: is_polygon_clockwise()

View file

@ -230,14 +230,15 @@ function select(list, start, end) =
: end==undef : end==undef
? is_num(start) ? is_num(start)
? list[ (start%l+l)%l ] ? list[ (start%l+l)%l ]
: assert( is_list(start) || is_range(start), "Invalid start parameter") : assert( start==[] || is_vector(start) || is_range(start), "Invalid start parameter")
[for (i=start) list[ (i%l+l)%l ] ] [for (i=start) list[ (i%l+l)%l ] ]
: assert(is_finite(start), "When `end` is given, `start` parameter should be a number.") : assert(is_finite(start), "When `end` is given, `start` parameter should be a number.")
assert(is_finite(end), "Invalid end parameter.") assert(is_finite(end), "Invalid end parameter.")
let( s = (start%l+l)%l, e = (end%l+l)%l ) let( s = (start%l+l)%l, e = (end%l+l)%l )
(s <= e) (s <= e)
? [for (i = [s:1:e]) list[i]] ? [ for (i = [s:1:e]) list[i] ]
: concat([for (i = [s:1:l-1]) list[i]], [for (i = [0:1:e]) list[i]]) ; : [ for (i = [s:1:l-1]) list[i],
for (i = [0:1:e]) list[i] ] ;
// Function: slice() // Function: slice()

View file

@ -278,9 +278,8 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) =
// signs at its two vertices can have an intersection with segment // signs at its two vertices can have an intersection with segment
// [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than
// eps and the sign of vals[j]-ref otherwise. // eps and the sign of vals[j]-ref otherwise.
signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) vals[j]-ref > eps ? 1 signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)])
: vals[j]-ref < -eps ? -1 abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ]
: 0]
) )
if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2] if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2]
for(j=[i+2:1:plen-(i==0 && closed? 3: 2)]) for(j=[i+2:1:plen-(i==0 && closed? 3: 2)])

View file

@ -427,9 +427,9 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru
for(p2=idx(region2)) for(p2=idx(region2))
let( let(
poly = closed2?close_path(region2[p2]):region2[p2], poly = closed2?close_path(region2[p2]):region2[p2],
signs = [for(v=poly*seg_normal) v-ref> eps ? 1 : v-ref<-eps ? -1 : 0] signs = [for(v=poly*seg_normal) abs(v-ref) < eps ? 0 : sign(v-ref) ]
) )
if(max(signs)>=0 && min(signs)<=0) // some edge edge intersects line [a1,a2] if(max(signs)>=0 && min(signs)<=0) // some edge intersects line [a1,a2]
for(j=[0:1:len(poly)-2]) for(j=[0:1:len(poly)-2])
if(signs[j]!=signs[j+1]) if(signs[j]!=signs[j+1])
let( // exclude non-crossing and collinear segments let( // exclude non-crossing and collinear segments

32
tests/test_all.scad Normal file
View file

@ -0,0 +1,32 @@
include <test_affine.scad>
include <test_attachments.scad>
include <test_comparisons.scad>
include <test_coords.scad>
include <test_cubetruss.scad>
include <test_distributors.scad>
include <test_drawing.scad>
include <test_edges.scad>
include <test_fnliterals.scad>
include <test_geometry.scad>
include <test_hull.scad>
include <test_linalg.scad>
include <test_linear_bearings.scad>
include <test_lists.scad>
include <test_math.scad>
include <test_mutators.scad>
include <test_paths.scad>
include <test_quaternions.scad>
include <test_regions.scad>
include <test_rounding.scad>
include <test_screw_drive.scad>
include <test_shapes2d.scad>
include <test_shapes3d.scad>
include <test_skin.scad>
include <test_strings.scad>
include <test_structs.scad>
include <test_transforms.scad>
include <test_trigonometry.scad>
include <test_utility.scad>
include <test_vectors.scad>
include <test_version.scad>
include <test_vnf.scad>

View file

@ -85,14 +85,14 @@ module test_polygon_triangulate() {
poly1 = [ [-10,0,-10], [10,0,10], [0,10,0], [-10,0,-10], [-4,4,-4], [4,4,4], [0,2,0], [-4,4,-4] ]; poly1 = [ [-10,0,-10], [10,0,10], [0,10,0], [-10,0,-10], [-4,4,-4], [4,4,4], [0,2,0], [-4,4,-4] ];
poly2 = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ]; poly2 = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ];
poly3 = [ [0,0], [10,0], [10,10], [10,13], [10,10], [0,10], [0,0], [3,3], [7,3], [7,7], [7,3], [3,3] ]; poly3 = [ [0,0], [10,0], [10,10], [10,13], [10,10], [0,10], [0,0], [3,3], [7,3], [7,7], [7,3], [3,3] ];
tris0 = sort(polygon_triangulate(poly0)); tris0 = (polygon_triangulate(poly0));
assert(approx(tris0, [[0, 1, 2]])); assert(approx(tris0, [[0, 1, 2]]));
tris1 = (polygon_triangulate(poly1)); tris1 = (polygon_triangulate(poly1));
assert(approx(tris1,( [[2, 3, 4], [6, 7, 0], [2, 4, 5], [6, 0, 1], [1, 2, 5], [5, 6, 1]]))); assert(approx(tris1,( [[2, 3, 4], [6, 7, 0], [2, 4, 5], [6, 0, 1], [1, 2, 5], [5, 6, 1]])));
tris2 = (polygon_triangulate(poly2)); tris2 = (polygon_triangulate(poly2));
assert(approx(tris2,([[0, 1, 2], [3, 4, 5]]))); assert(approx(tris2,( [[3, 4, 5], [1, 2, 3]])));
tris3 = (polygon_triangulate(poly3)); tris3 = (polygon_triangulate(poly3));
assert(approx(tris3,( [[5, 6, 7], [7, 8, 9], [10, 11, 0], [5, 7, 9], [10, 0, 1], [4, 5, 9], [9, 10, 1], [1, 4, 9]]))); assert(approx(tris3,( [[5, 6, 7], [11, 0, 1], [5, 7, 8], [10, 11, 1], [5, 8, 9], [10, 1, 2], [4, 5, 9], [9, 10, 2]])));
} }
module test__normalize_plane(){ module test__normalize_plane(){

View file

@ -491,12 +491,11 @@ function _bt_tree(points, ind, leafsize=25) =
bounds = pointlist_bounds(select(points,ind)), bounds = pointlist_bounds(select(points,ind)),
coord = max_index(bounds[1]-bounds[0]), coord = max_index(bounds[1]-bounds[0]),
projc = [for(i=ind) points[i][coord] ], projc = [for(i=ind) points[i][coord] ],
pmc = mean(projc), meanpr = mean(projc),
pivot = min_index([for(p=projc) abs(p-pmc)]), pivot = min_index([for(p=projc) abs(p-meanpr)]),
radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]), radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]),
median = median(projc), Lind = [for(i=idx(ind)) if(projc[i]<=meanpr && i!=pivot) ind[i] ],
Lind = [for(i=idx(ind)) if(projc[i]<=median && i!=pivot) ind[i] ], Rind = [for(i=idx(ind)) if(projc[i] >meanpr && i!=pivot) ind[i] ]
Rind = [for(i=idx(ind)) if(projc[i] >median && i!=pivot) ind[i] ]
) )
[ ind[pivot], radius, _bt_tree(points, Lind, leafsize), _bt_tree(points, Rind, leafsize) ]; [ ind[pivot], radius, _bt_tree(points, Lind, leafsize), _bt_tree(points, Rind, leafsize) ];

172
vnf.scad
View file

@ -318,14 +318,13 @@ function vnf_merge(vnfs, cleanup=false, eps=EPSILON) =
cleanup? _vnf_cleanup(verts,faces,eps) : [verts,faces]; cleanup? _vnf_cleanup(verts,faces,eps) : [verts,faces];
function _vnf_cleanup(verts,faces,eps) = function _vnf_cleanup(verts,faces,eps) =
let( let(
dedup = vector_search(verts,eps,verts), // collect vertex duplicates dedup = vector_search(verts,eps,verts), // collect vertex duplicates
map = [for(i=idx(verts)) min(dedup[i]) ], // remap duplic vertices map = [for(i=idx(verts)) min(dedup[i]) ], // remap duplic vertices
offset = cumsum([for(i=idx(verts)) map[i]==i ? 0 : 1 ]), // remaping face vertex offsets offset = cumsum([for(i=idx(verts)) map[i]==i ? 0 : 1 ]), // remaping face vertex offsets
map2 = list(idx(verts))-offset, // map old vertex indices to new indices map2 = list(idx(verts))-offset, // map old vertex indices to new indices
nverts = [for(i=idx(verts)) if(map[i]==i) verts[i] ], // eliminates all unreferenced vertices nverts = [for(i=idx(verts)) if(map[i]==i) verts[i] ], // this doesn't eliminate unreferenced vertices
nfaces = nfaces =
[ for(face=faces) [ for(face=faces)
let( let(
@ -385,35 +384,124 @@ function _join_paths_at_vertices(path1,path2,v1,v2) =
]; ];
// Given a region that is connected and has its outer border in region[0], /// Internal Function: _cleave_connected_region(region, eps)
// produces a polygon with the same points that has overlapping connected paths /// Description:
// to join internal holes to the outer border. Output is a single path. /// Given a region that is connected and has its outer border in region[0],
function _cleave_connected_region(region) = /// produces a overlapping connected path to join internal holes to
len(region)==0? [] : /// the outer border without adding points. Output is a single non-simple polygon.
len(region)<=1? clockwise_polygon(region[0]) : /// Requirements:
let( /// It expects that all region paths be simple closed paths, with region[0] CW and
dists = [ /// the other paths CCW and encircled by region[0]. The input region paths are also
for (i=[1:1:len(region)-1]) /// supposed to be disjoint except for common vertices and common edges but with
_path_path_closest_vertices(region[0],region[i]) /// no crossings. It may return `undef` if these conditions are not met.
], /// This function implements an extension of the algorithm discussed in:
idxi = min_index(column(dists,0)), /// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
newoline = _join_paths_at_vertices( function _cleave_connected_region(region, eps=EPSILON) =
region[0], region[idxi+1], len(region)==1 ? region[0] :
dists[idxi][1], dists[idxi][2] let(
) outer = deduplicate(region[0]), //
) len(region)==2? clockwise_polygon(newoline) : holes = [for(i=[1:1:len(region)-1]) // deduplication possibly unneeded
let( deduplicate( region[i] ) ], //
orgn = [ extridx = [for(li=holes) max_index(column(li,0)) ],
newoline, // the right extreme vertex for each hole sorted by decreasing x values
for (i=idx(region)) extremes = sort( [for(i=idx(holes)) [ i, extridx[i], -holes[i][extridx[i]].x] ], idx=2 )
if (i>0 && i!=idxi+1) )
region[i] _polyHoles(outer, holes, extremes, eps, 0);
]
)
assert(len(orgn)<len(region))
_cleave_connected_region(orgn);
// connect the hole paths one at a time to the outer path.
// 'extremes' is the list of the right extreme vertex of each hole sorted by decreasing abscissas
// see: _cleave_connected_region(region, eps)
function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) =
let(
extr = extremes[n], //
hole = holes[extr[0]], // hole path to bridge to the outer path
ipt = extr[1], // index of the hole point with maximum abscissa
brdg = _bridge(hole[ipt], outer, eps) // the index of a point in outer to bridge hole[ipt] to
)
brdg == undef ? undef :
let(
l = len(outer),
lh = len(hole),
// the new outer polygon bridging the hole to the old outer
npoly =
approx(outer[brdg], hole[ipt], eps)
? [ for(i=[brdg: 1: brdg+l]) outer[i%l] ,
for(i=[ipt+1: 1: ipt+lh-1]) hole[i%lh] ]
: [ for(i=[brdg: 1: brdg+l]) outer[i%l] ,
for(i=[ipt: 1: ipt+lh]) hole[i%lh] ]
)
n==len(holes)-1 ? npoly :
_polyHoles(npoly, holes, extremes, eps, n+1);
// find a point in outer to be connected to pt in the interior of outer
// by a segment that not cross or touch any non adjacente edge of outer.
// return the index of a vertex in the outer path where the bridge should end
// see _polyHoles(outer, holes, extremes, eps)
function _bridge(pt, outer,eps) =
// find the intersection of a ray from pt to the right
// with the boundary of the outer cycle
let(
l = len(outer),
crxs =
let( edges = pair(outer,wrap=true) )
[for( i = idx(edges) )
let( edge = edges[i] )
// consider just descending outer edges at right of pt crossing ordinate pt.y
if( (edge[0].y > pt.y+eps)
&& (edge[1].y <= pt.y)
&& _is_at_left(pt, [edge[1], edge[0]], eps) )
[ i,
// the point of edge with ordinate pt.y
abs(pt.y-edge[1].y)<eps ? edge[1] :
let( u = (pt-edge[1]).y / (edge[0]-edge[1]).y )
(1-u)*edge[1] + u*edge[0]
]
]
)
crxs == [] ? undef :
let(
// the intersection point of the nearest edge to pt with minimum slope
minX = min([for(p=crxs) p[1].x]),
crxcand = [for(crx=crxs) if(crx[1].x < minX+eps) crx ], // nearest edges
nearest = min_index([for(crx=crxcand)
(outer[crx[0]].x - pt.x) / (outer[crx[0]].y - pt.y) ]), // minimum slope
proj = crxcand[nearest],
vert0 = outer[proj[0]], // the two vertices of the nearest crossing edge
vert1 = outer[(proj[0]+1)%l],
isect = proj[1] // the intersection point
)
norm(pt-vert1) < eps ? (proj[0]+1)%l : // if pt touches an outer vertex, return its index
// as vert0.y > pt.y then pt!=vert0
norm(pt-isect) < eps ? undef : // if pt touches the middle of an outer edge -> error
let(
// the edge [vert0, vert1] necessarily satisfies vert0.y > vert1.y
// indices of candidates to an outer bridge point
cand =
(vert0.x > pt.x)
? [ proj[0],
// select reflex vertices inside of the triangle [pt, vert0, isect]
for(i=idx(outer))
if( _tri_class(select(outer,i-1,i+1),eps) <= 0
&& _pt_in_tri(outer[i], [pt, vert0, isect], eps)>=0 )
i
]
: [ (proj[0]+1)%l,
// select reflex vertices inside of the triangle [pt, isect, vert1]
for(i=idx(outer))
if( _tri_class(select(outer,i-1,i+1),eps) <= 0
&& _pt_in_tri(outer[i], [pt, isect, vert1], eps)>=0 )
i
],
// choose the candidate outer[i] such that the line [pt, outer[i]] has minimum slope
// among those with minimum slope choose the nearest to pt
slopes = [for(i=cand) 1-abs(outer[i].x-pt.x)/norm(outer[i]-pt) ],
min_slp = min(slopes),
cand2 = [for(i=idx(cand)) if(slopes[i]<=min_slp+eps) cand[i] ],
nearest = min_index([for(i=cand2) norm(pt-outer[i]) ])
)
cand2[nearest];
// Function: vnf_from_region() // Function: vnf_from_region()
// Usage: // Usage:
@ -436,13 +524,15 @@ function _cleave_connected_region(region) =
function vnf_from_region(region, transform, reverse=false) = function vnf_from_region(region, transform, reverse=false) =
let ( let (
regions = region_parts(force_region(region)), regions = region_parts(force_region(region)),
vnfs = [ vnfs =
for (rgn = regions) let( [ for (rgn = regions)
cleaved = path3d(_cleave_connected_region(rgn)), let( cleaved = path3d(_cleave_connected_region(rgn)) )
face = is_undef(transform)? cleaved : apply(transform,cleaved), assert( cleaved, "The region is invalid")
faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i] let(
) [face, [faceidxs]] face = is_undef(transform)? cleaved : apply(transform,cleaved),
], faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i]
) [face, [faceidxs]]
],
outvnf = vnf_merge(vnfs) outvnf = vnf_merge(vnfs)
) )
vnf_triangulate(outvnf); vnf_triangulate(outvnf);
@ -550,9 +640,13 @@ function _link_indicator(l,imin,imax) =
function vnf_triangulate(vnf) = function vnf_triangulate(vnf) =
let( let(
verts = vnf[0], verts = vnf[0],
faces = [for (face=vnf[1]) each len(face)==3 ? [face] : faces = [for (face=vnf[1])
polygon_triangulate(verts, face)] each (len(face)==3 ? [face] :
) [verts, faces]; let( tris = polygon_triangulate(verts, face) )
assert( tris!=undef, "Some `vnf` face cannot be triangulated.")
tris ) ]
)
[verts, faces];