diff --git a/geometry.scad b/geometry.scad index 1e0df7b..6c916d9 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1603,17 +1603,25 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // Description: // 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 -// given, the it is used as an index list into `poly` to define the polygon. In that case, -// `poly` may have a length greater than `ind`. Otherwise, all points in `poly` +// given, it is used as an index list into `poly` to define the polygon. In that case, +// `poly` may have a length greater than `ind`. When `ind` is undefined, all points in `poly` // are considered as vertices of the polygon. // . -// The function may issue an error if it finds that the polygon is not simple -// (self-intersecting) or its vertices are collinear. It can work for 3d non-planar polygons -// if they are close enough to planar but may otherwise issue an error for this case. -// . // For 2d polygons, the output triangles will have the same winding (CW or CCW) of // the input polygon. For 3d polygons, the triangle windings will induce a normal // vector with the same direction of the polygon normal. +// . +// 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 +// simple polygons with the same winding. These polygons may have "touching" vertices +// (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 +// no self-crossing. See examples bellow. If all polygon edges are contact edges, returns an empty list. +// . +// Self-crossing polygons have no consistent winding and usually produce an error but +// when an error is not issued 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 +// issue an error for this case. // Arguments: // poly = Array of vertices for the polygon. // ind = A list indexing the vertices of the polygon in `poly`. @@ -1621,7 +1629,28 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // Example(2D,NoAxes): // poly = star(id=10, od=15,n=11); // tris = polygon_triangulate(poly); -// for(tri=tris) stroke(select(poly,tri), width=.2, closed=true); +// color("lightblue") for(tri=tris) polygon(select(poly,tri)); +// color("magenta") up(2) stroke(poly,.25,closed=true); +// color("black") up(3) vnf_debug([poly,[]],faces=false,size=1); +// Example(2D,NoAxes): a polygon with a hole and one "contact" edge +// poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ]; +// tris = polygon_triangulate(poly); +// color("lightblue") for(tri=tris) polygon(select(poly,tri)); +// color("magenta") up(2) stroke(poly,.25,closed=true); +// color("black") up(3) vnf_debug([poly,[]],faces=false,size=1); +// Example(2D,NoAxes): a polygon with "touching" vertices and no holes +// poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ]; +// tris = polygon_triangulate(poly); +// color("lightblue") for(tri=tris) polygon(select(poly,tri)); +// color("magenta") up(2) stroke(poly,.25,closed=true); +// color("black") up(3) vnf_debug([poly,[]],faces=false,size=1); +// Example(2D,NoAxes): a polygon with "contact" edges and no holes +// poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3], +// [7,7], [7,3], [3,3] ]; +// tris = polygon_triangulate(poly); // see from the top +// color("lightblue") for(tri=tris) polygon(select(poly,tri)); +// color("magenta") up(2) stroke(poly,.25,closed=true); +// color("black") up(3) vnf_debug([poly,[]],faces=false,size=1); // Example(3D): // include // vnf = regular_polyhedron_info(name="dodecahedron",side=5,info="vnf"); @@ -1630,82 +1659,43 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // color("blue") // vnf_wireframe(vnf_tri, width=.15); function polygon_triangulate(poly, ind, eps=EPSILON) = - assert(is_path(poly), "Polygon `poly` should be a list of 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) || (is_vector(ind) && min(ind)>=0 && max(ind) 2*eps, + "The polygon vertices are collinear.") + [ind] + : len(poly[ind[0]]) == 3 + ? // represents the polygon projection on its plane as a 2d polygon + let( + ind = deduplicate_indexed(poly, ind, eps) + ) + len(ind)<3 ? [] : + let( + pts = select(poly,ind), + nrm = polygon_normal(pts) + ) + // here, instead of an error, it might return [] or undef + assert( nrm!=undef, + "The polygon has self-intersections or its vertices are collinear or non coplanar.") + let( + imax = max_index([for(p=pts) norm(p-pts[0]) ]), + v1 = unit( pts[imax] - pts[0] ), + v2 = cross(v1,nrm), + prpts = pts*transpose([v1,v2]) + ) + [for(tri=_triangulate(prpts, count(len(ind)), eps)) select(ind,tri) ] + : let( cw = is_polygon_clockwise(select(poly, ind)) ) + cw + ? [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ] + : _triangulate( poly, ind, eps ); - -// requires ccw 2d polygons -// returns ccw triangles -function _old_triangulate(poly, ind, eps=EPSILON, tris=[]) = - len(ind)==3 ? concat(tris,[ind]) : - let( ear = _get_ear(poly,ind,eps) ) - assert( ear!=undef, - "The polygon has self-intersections or its vertices are collinear or non coplanar.") - let( - ear_tri = select(ind,ear,ear+2), - indr = select(ind,ear+2, ear) // indices of the remaining points - ) - _triangulate(poly, indr, eps, concat(tris,[ear_tri])); - -// search a valid ear from the remaining polygon -function _old_get_ear(poly, ind, eps, _i=0) = - _i>=len(ind) ? undef : // poly has no ears - let( // the _i-th ear candidate - p0 = poly[ind[_i]], - p1 = poly[ind[(_i+1)%len(ind)]], - p2 = poly[ind[(_i+2)%len(ind)]] - ) - // if it is not a convex vertex, try the next one - _is_cw2(p0,p1,p2,eps) ? _get_ear(poly,ind,eps, _i=_i+1) : - let( // vertex p1 is convex; check if the triangle contains any other point - to_tst = select(ind,_i+3, _i-1), - pt2tst = select(poly,to_tst), // points other than p0, p1 and p2 - q = [(p0-p2).y, (p2-p0).x], // orthogonal to ray [p0,p2] pointing right - q0 = q*p0, - atleft = [for(p=pt2tst) if(p*q<=q0) p ] - ) - atleft==[] ? _i : // no point inside -> an ear - let( - q = [(p2-p1).y, (p1-p2).x], // orthogonal to ray [p1,p2] pointing right - q0 = q*p2, - atleft = [for(p=atleft) if(p*q<=q0) p ] - ) - atleft==[] ? _i : // no point inside -> an ear - let( - q = [(p1-p0).y, (p0-p1).x], // orthogonal to ray [p1,p0] pointing right - q0 = q*p1, - atleft = [for(p=atleft) if(p*q<=q0) p ] - ) - atleft==[] ? _i : // no point inside -> an ear - // check the next ear candidate - _get_ear(poly, ind, eps, _i=_i+1); - function _triangulate(poly, ind, eps=EPSILON, tris=[]) = len(ind)==3 ? _is_degenerate(select(poly,ind),eps) @@ -1714,18 +1704,18 @@ function _triangulate(poly, ind, eps=EPSILON, tris=[]) = : let( ear = _get_ear(poly,ind,eps) ) assert( ear!=undef, "The polygon has self-intersections or its vertices are collinear or non coplanar.") - ear<0 // degenerate ear - ? let( indr = select(ind,-ear+1, -ear-1) ) // discard it - _triangulate(poly, indr, eps, tris) + is_list(ear) // degenerate ear + ? _triangulate(poly, select(ind,ear[0]+2, ear[0]), eps, tris) // discard it : let( ear_tri = select(ind,ear,ear+2), - indr = select(ind,ear+2, ear) // indices of the remaining points + indr = select(ind,ear+2, ear) // remaining point indices ) _triangulate(poly, indr, eps, concat(tris,[ear_tri])); // a returned ear will be: -// 1. a CCW triangle without points inside except possibly at its vertices +// 1. a CCW (non-degenerate) triangle, made of subsequent vertices, without other +// points inside except possibly at its vertices // 2. or a degenerate triangle where two vertices are coincident // the returned ear is specified by the index of `ind` of its first vertex function _get_ear(poly, ind, eps, _i=0) = @@ -1735,29 +1725,25 @@ function _get_ear(poly, ind, eps, _i=0) = p1 = poly[ind[(_i+1)%len(ind)]], p2 = poly[ind[(_i+2)%len(ind)]] ) - // if it is a degenerate triangle, return it (codified) - _is_degenerate([p0,p1,p2],eps) ? -(_i+1) : - // if it is not a convex vertex, try the next one + // degenerate triangles are returned codified + _is_degenerate([p0,p1,p2],eps) ? [_i] : + // if it is not a convex vertex, check the next one _is_cw2(p0,p1,p2,eps) ? _get_ear(poly,ind,eps, _i=_i+1) : let( // vertex p1 is convex // check if the triangle contains any other point // except possibly its own vertices to_tst = select(ind,_i+3, _i-1), - pt2tst = select(poly,to_tst), // points other than p0, p1 and p2 q = [(p0-p2).y, (p2-p0).x], // orthogonal to ray [p0,p2] pointing right - q0 = q*p0, r = [(p2-p1).y, (p1-p2).x], // orthogonal to ray [p2,p1] pointing right - r0 = r*p2, s = [(p1-p0).y, (p0-p1).x], // orthogonal to ray [p1,p0] pointing right - s0 = s*p1, - inside = [for(p=pt2tst) - if( p*q<=q0 && p*r<=r0 && p*s<=s0 ) // p is in the triangle - if( norm(p-p0)>eps // and doesn't coincide with - && norm(p-p1)>eps // any of its vertices + 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 && norm(p-p2)>eps ) p ] ) - inside==[] ? _i : // no point inside -> an ear + inside==[] ? _i : // found an ear // check the next ear candidate _get_ear(poly, ind, eps, _i=_i+1); diff --git a/math.scad b/math.scad index 4f46258..360989d 100644 --- a/math.scad +++ b/math.scad @@ -689,17 +689,12 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1); // cumsum([[1,2,3], [3,4,5], [5,6,7]]); // returns [[1,2,3], [4,6,8], [9,12,15]] function cumsum(v) = assert(is_consistent(v), "The input is not consistent." ) - _cumsum(v,_i=0,_acc=[]); + len(v)<=1 ? v : + _cumsum(v,_i=1,_acc=[v[0]]); function _cumsum(v,_i=0,_acc=[]) = - _i==len(v) ? _acc : - _cumsum( - v, _i+1, - concat( - _acc, - [_i==0 ? v[_i] : last(_acc) + v[_i]] - ) - ); + _i>=len(v) ? _acc : + _cumsum( v, _i+1, [ each _acc, _acc[len(_acc)-1] + v[_i] ] ); // Function: sum_of_sines() diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 60359b6..d7405df 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -48,6 +48,7 @@ test_reindex_polygon(); test_align_polygon(); test_polygon_centroid(); test_point_in_polygon(); +test_polygon_triangulate(); test_is_polygon_clockwise(); test_clockwise_polygon(); test_ccw_polygon(); @@ -78,6 +79,21 @@ function info_str(list,i=0,string=chr(10)) = : info_str(list,i+1,str(string,str(list[i][0],_valstr(list[i][1]),chr(10)))); +module test_polygon_triangulate() { + poly0 = [ [0,0,1], [10,0,2], [10,10,0] ]; + 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] ]; + 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)); + assert(approx(tris0, [[0, 1, 2]])); + 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]]))); + tris2 = (polygon_triangulate(poly2)); + assert(approx(tris2,([[0, 1, 2], [3, 4, 5]]))); + 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]]))); +} + module test__normalize_plane(){ plane = rands(-5,5,4,seed=333)+[10,0,0,0]; plane2 = _normalize_plane(plane);