From 8c4022fa3c98e19fa5c5f3d660c14e55448bd24d Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 3 Nov 2021 20:57:35 +0000 Subject: [PATCH 1/5] Corrections in _cleave_connected_region and polygon_triangulate and some minor tweaks, --- geometry.scad | 86 +++++++++++++++++++++++++++------------------------ lists.scad | 25 ++++++++++++--- paths.scad | 6 +--- regions.scad | 6 ++-- vnf.scad | 75 ++++++++++++++++++++++++-------------------- 5 files changed, 112 insertions(+), 86 deletions(-) diff --git a/geometry.scad b/geometry.scad index c09604b..e983f7a 100644 --- a/geometry.scad +++ b/geometry.scad @@ -104,15 +104,16 @@ function _tri_class(tri, eps=EPSILON) = /// class = _pt_in_tri(point, tri); /// Topics: Geometry, Points, Triangles /// Description: -/// Return 1 if point is inside the triangle interion. -/// Return =0 if point is on the triangle border. -/// Return -1 if point is outside the triangle. +/// 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), + 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) ); @@ -1695,7 +1696,7 @@ 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, 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` // are considered as vertices of the polygon. // . @@ -1704,47 +1705,50 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // 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 +// A polygon is non-twisted iff it is simple or it has a partition in // simple polygons with the same winding such that the intersection of any two partitions is -// made of full edges of both partitions. These polygons may have "touching" vertices +// 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 // (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 (polygons with -// zero area), 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. // . -// Twisted polygons have no consistent winding and when input to this function usually produce -// an error but when an error is not issued the outputs are not correct triangulations. The function +// Twisted polygons have no consistent winding and when input to this function usually reports +// 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 -// issue an error for this case. +// report an error for this case. // Arguments: // poly = Array of the polygon vertices. // ind = A list indexing the vertices of the polygon in `poly`. +// error = If false, reports triangulation errors returning `undef` otherwise report them by an assert error. Default: true // 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); // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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] ]; // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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] ]; // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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], // [7,7], [7,3], [3,3] ]; -// tris = polygon_triangulate(poly); // see from above +// tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); @@ -1756,18 +1760,18 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ]; // color("blue") // 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_undef(ind) - || (is_vector(ind) && min(ind)>=0 && max(ind)=0 && max(ind) 2*eps, - "The polygon vertices are collinear.") + let( degen = norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) < 2*eps ) + assert( ! error || ! degen, "The polygon vertices are collinear.") + degen ? undef : [ind] : len(poly[ind[0]]) == 3 ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane @@ -1779,44 +1783,47 @@ function polygon_triangulate(poly, ind, eps=EPSILON) = pts = select(poly,ind), nrm = -polygon_normal(pts) ) - assert( nrm!=undef, + assert( ! error || (nrm != undef), "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.") + nrm == undef ? undef : 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]) // 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) ) + tris == undef ? undef : + [for(tri=tris) select(ind,tri) ] : is_polygon_clockwise(select(poly, ind)) - ? _triangulate( poly, ind, eps ) - : [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ]; + ? _triangulate( poly, ind, error, eps ) + : let( tris = _triangulate( poly, reverse(ind), error, eps ) ) + tris == undef ? undef : + [for(tri=tris) reverse(tri) ]; // 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 (tecnically) the ones whose interior -// is homeomoph to the interior of a simple polygon -function _triangulate(poly, ind, eps=EPSILON, tris=[]) = +// is homeomoph to the interior of a circle +function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) = +//echo(eps=eps) len(ind)==3 ? _degenerate_tri(select(poly,ind),eps) ? tris // if last 3 pts perform a degenerate triangle, ignore it : concat(tris,[ind]) // otherwise, include it : let( ear = _get_ear(poly,ind,eps) ) -/* -let( x= [if(is_undef(ear)) echo(ind=ind) 0] ) -is_undef(ear) ? tris : -*/ - assert( ear!=undef, + assert( ! error || (ear != undef), "The polygon has twists or all its vertices are collinear or non coplanar.") + ear == undef ? undef : is_list(ear) // is it a degenerate ear ? ? len(ind) <= 4 ? tris : - _triangulate(poly, select(ind,ear[0]+3, ear[0]), eps, tris) // discard it + _triangulate(poly, select(ind,ear[0]+3, ear[0]), error, eps, tris) // discard it : let( ear_tri = select(ind,ear,ear+2), 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: @@ -1825,6 +1832,7 @@ is_undef(ear) ? tris : // 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) = +//let( x = ! is_num(eps) ? echo(ind=ind,eps=eps) 0:0 ) let( lind = len(ind) ) lind==3 ? 0 : let( // the _i-th ear candidate @@ -1841,9 +1849,7 @@ function _get_ear(poly, ind, eps, _i=0) = // otherwise check the next ear candidate _i=0 ) // 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,next_vert],eps) ) + && _is_at_left(p0,[prev_vert,p1],eps) && _is_at_left(p2,[p1,prev_vert],eps) + && _is_at_left(p2,[p1,next_vert],eps) && _is_at_left(p0,[next_vert,p1],eps) ) ) ? false diff --git a/lists.scad b/lists.scad index 07356e9..70b098f 100644 --- a/lists.scad +++ b/lists.scad @@ -230,14 +230,15 @@ function select(list, start, end) = : end==undef ? is_num(start) ? 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 ] ] : assert(is_finite(start), "When `end` is given, `start` parameter should be a number.") assert(is_finite(end), "Invalid end parameter.") let( s = (start%l+l)%l, e = (end%l+l)%l ) (s <= e) - ? [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:e]) list[i] ] + : [ for (i = [s:1:l-1]) list[i], + for (i = [0:1:e]) list[i] ] ; // Function: slice() @@ -952,7 +953,7 @@ function enumerate(l,idx=undef) = // Example: // l = ["A","B","C","D"]; // echo([for (p=pair(l)) str(p.y,p.x)]); // Outputs: ["BA", "CB", "DC"] -function pair(list, wrap=false) = +function _old_pair(list, wrap=false) = assert(is_list(list)||is_string(list), "Invalid input." ) assert(is_bool(wrap)) let( @@ -961,6 +962,12 @@ function pair(list, wrap=false) = ? [for (i=[0:1:ll-1]) [list[i], list[(i+1) % ll]]] : [for (i=[0:1:ll-2]) [list[i], list[i+1]]]; +function pair(list, wrap=false) = + assert(is_list(list)||is_string(list), "Invalid input." ) + assert(is_bool(wrap)) + let( ll = len(list)-1 ) + [for (i=[0:1:ll-1]) [list[i], list[i+1]], if(wrap) [list[ll], list[0]] ]; + // Function: triplet() // Usage: @@ -981,7 +988,7 @@ function pair(list, wrap=false) = // translate(b) rot(from=FWD,to=v) anchor_arrow2d(); // } // stroke(path); -function triplet(list, wrap=false) = +function _old_triplet(list, wrap=false) = assert(is_list(list)||is_string(list), "Invalid input." ) assert(is_bool(wrap)) let( @@ -991,6 +998,14 @@ function triplet(list, wrap=false) = : [for (i=[0:1:ll-3]) [ list[i], list[i+1], list[i+2] ]]; +function triplet(list, wrap=false) = + assert(is_list(list)||is_string(list), "Invalid input." ) + assert(is_bool(wrap)) + let( ll = len(list) ) + [ for (i=[0:1:ll-3]) [ list[i], list[i+1], list[i+2] ], + each if(wrap) [ [ list[ll-2], list[ll-1], list[0] ], [ list[ll-1], list[0], list[1] ] ] ]; + + // Function: combinations() // Usage: // list = combinations(l, [n]); diff --git a/paths.scad b/paths.scad index fe7fd86..7a59432 100644 --- a/paths.scad +++ b/paths.scad @@ -274,9 +274,7 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = // 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 // 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 - : vals[j]-ref < -eps ? -1 - : 0] + signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ] ) 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)]) @@ -286,10 +284,8 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) ) if (isect -// && isect[1]> (i==0 && !closed? -eps: 0) // Apparently too strict && isect[1]>=-eps && isect[1]<= 1+eps -// && isect[2]> 0 && isect[2]>= -eps && isect[2]<= 1+eps) [isect[0], i, isect[1], j, isect[2]] diff --git a/regions.scad b/regions.scad index 662d860..5f6d442 100644 --- a/regions.scad +++ b/regions.scad @@ -216,9 +216,9 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru for(p2=idx(region2)) let( 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]) if(signs[j]!=signs[j+1]) let( // exclude non-crossing and collinear segments @@ -260,7 +260,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru // where region1 intersections region2. Split region2 similarly with respect to region1. // The return is a pair of results of the form [split1, split2] where split1=[frags1,frags2,...] // and frags1 is a list of path pieces (in order) from the first path of the region. -// You can pass a single path in for either region, but the output will be a singleton list, as ify +// You can pass a single path in for either region, but the output will be a singleton list, as if // you passed in a singleton region. // Arguments: // region1 = first region diff --git a/vnf.scad b/vnf.scad index b2325c1..dc4dfd6 100644 --- a/vnf.scad +++ b/vnf.scad @@ -414,15 +414,15 @@ function _old_cleave_connected_region(region) = /// Internal Function: _cleave_connected_region(region, eps) /// Description: -/// Given a region that is connected and has its outer border in region[0], -/// produces a polygon with the same points that has overlapping connected paths -/// to join internal holes to the outer border. Output is a single path. -/// It expect that region[0] be a simple closed CW path and that each hole, -/// region[i] for i>0, be a simple closed CCW path. -/// The paths are also supposed to be disjoint except for common vertices and -/// common edges but no crossing. -/// This function implements an extension of the algorithm discussed in: -/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf +/// Given a region that is connected and has its outer border in region[0], +/// produces a overlapping connected path to join internal holes to +/// the outer border without adding points. Output is a single non-simple polygon. +/// It expect that all region paths be simple closed paths with arbitrary widings. +/// The input region paths are also supposed to be disjoint except for common +/// vertices and common edges but with no crossings. It may return `undef` if +/// these conditions are not met. +/// This function implements an extension of the algorithm discussed in: +/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf function _cleave_connected_region(region, eps=EPSILON) = len(region)==1 ? region[0] : let( @@ -439,7 +439,7 @@ function _cleave_connected_region(region, eps=EPSILON) = // 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) +// see: _cleave_connected_region(region, eps) function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = let( extr = extremes[n], // @@ -447,17 +447,18 @@ function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = 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 ) - assert(brdg!=undef, "Error: check input polygon restrictions") +// assert(brdg!=undef, "Error: check input polygon restrictions") + 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] ] + ? [ 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); @@ -472,13 +473,12 @@ function _bridge(pt, outer,eps) = let( l = len(outer), crxs = - [for( i=idx(outer) ) - let( edge = select(outer,i,i+1) ) + let( edges = pair(outer,wrap=true) ) + [for( i = idx(edges), edge=[edges[i]] ) // consider just descending outer edges at right of pt crossing ordinate pt.y - if( (edge[0].y> pt.y) - && (edge[1].y<=pt.y) - && ( norm(edge[1]-pt)0 ) ) + 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) error - assert( ! approx(pt,isect,eps) || approx(pt,vert0,eps) || approx(pt,vert1,eps), - "There is a forbidden self_intersection" ) - norm(pt-vert0) < eps ? proj[0] : // if pt touches an outer vertex, return its index - norm(pt-vert1) < eps ? (proj[0]+1)%l : + // as vert0.y > pt.y then pt!=vert0 + norm(pt-vert1) < eps ? (proj[0]+1)%l : // if pt touches an outer vertex, return its index + +// assert( ! approx(pt,isect,eps) , // if pt touches the middle of an outer edge -> error +// "A hole vertex is in the middle of an edge" ) + 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 @@ -555,7 +558,10 @@ function vnf_from_region(region, transform, reverse=false) = regions = region_parts(force_region(region)), vnfs = [ for (rgn = regions) let( - cleaved = path3d(_cleave_connected_region(rgn)), + cleaved = path3d(_cleave_connected_region(rgn)) + ) + assert( cleaved, "The region is invalid") + let( 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]] @@ -668,8 +674,11 @@ function vnf_triangulate(vnf) = let( verts = vnf[0], faces = [for (face=vnf[1]) each len(face)==3 ? [face] : - polygon_triangulate(verts, face)] - ) [verts, faces]; + polygon_triangulate(verts, face)], + invalid = [for(face=faces) if(face==undef) 1 ] + ) + assert( invalid==[], "Some `vnf` face cannot be triangulated.") + [verts, faces]; From 36968a463b638c34c9d5f2a9a48e9cf499515907 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Thu, 4 Nov 2021 07:56:34 +0000 Subject: [PATCH 2/5] Revert "Corrections in _cleave_connected_region and polygon_triangulate and some minor tweaks," This reverts commit 8c4022fa3c98e19fa5c5f3d660c14e55448bd24d. --- geometry.scad | 86 ++++++++++++++++++++++++--------------------------- lists.scad | 25 +++------------ paths.scad | 6 +++- regions.scad | 6 ++-- vnf.scad | 75 ++++++++++++++++++++------------------------ 5 files changed, 86 insertions(+), 112 deletions(-) diff --git a/geometry.scad b/geometry.scad index e983f7a..c09604b 100644 --- a/geometry.scad +++ b/geometry.scad @@ -104,16 +104,15 @@ function _tri_class(tri, eps=EPSILON) = /// 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. +/// Return 1 if point is inside the triangle interion. +/// 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), + 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) ); @@ -1696,7 +1695,7 @@ 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, it is used as an index list into `poly` to define the polygon vertices. In that case, +// 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. // . @@ -1705,50 +1704,47 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // 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 it has a partition in +// A polygon is non-twisted iff it is simple or there is a partition of it in // 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 +// made of full edges of both partitions. 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 (polygons with -// 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. +// zero area), it returns an empty list for 2d polygons and issues an error for 3d polygons. // . -// Twisted polygons have no consistent winding and when input to this function usually reports -// an error but when an error is not reported the outputs are not correct triangulations. The function +// Twisted polygons have no consistent winding and when input to this function 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 -// report an error for this case. +// issue an error for this case. // Arguments: // poly = Array of the polygon vertices. // ind = A list indexing the vertices of the polygon in `poly`. -// error = If false, reports triangulation errors returning `undef` otherwise report them by an assert error. Default: true // eps = A maximum tolerance in geometrical tests. Default: EPSILON -// Example(2D,NoAxes): a simple polygon; see from above +// Example(2D,NoAxes): // poly = star(id=10, od=15,n=11); // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); -// Example(2D,NoAxes): a polygon with a hole and one "contact" edge; see from above +// 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("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); -// Example(2D,NoAxes): a polygon with "touching" vertices and no holes; see from above +// 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("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // color("black") up(3) vnf_debug([path3d(poly),[]],faces=false,size=1); -// Example(2D,NoAxes): a polygon with "contact" edges and no holes; see from above +// 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); +// tris = polygon_triangulate(poly); // see from above // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); @@ -1760,18 +1756,18 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ]; // color("blue") // vnf_wireframe(vnf_tri, width=.15); -function polygon_triangulate(poly, ind, error=true, eps=EPSILON) = +function polygon_triangulate(poly, ind, 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_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)=0 && max(ind) 2*eps, + "The polygon vertices are collinear.") [ind] : len(poly[ind[0]]) == 3 ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane @@ -1783,47 +1779,44 @@ function polygon_triangulate(poly, ind, error=true, eps=EPSILON) = pts = select(poly,ind), nrm = -polygon_normal(pts) ) - assert( ! error || (nrm != undef), + assert( nrm!=undef, "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.") - nrm == undef ? undef : 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]) // the 2d projection of pts on the polygon plane ) - let( tris = _triangulate(prpts, count(len(ind)), error, eps) ) - tris == undef ? undef : - [for(tri=tris) select(ind,tri) ] + [for(tri=_triangulate(prpts, count(len(ind)), eps)) select(ind,tri) ] : is_polygon_clockwise(select(poly, ind)) - ? _triangulate( poly, ind, error, eps ) - : let( tris = _triangulate( poly, reverse(ind), error, eps ) ) - tris == undef ? undef : - [for(tri=tris) reverse(tri) ]; + ? _triangulate( poly, ind, eps ) + : [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ]; // 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 (tecnically) the ones whose interior -// is homeomoph to the interior of a circle -function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) = -//echo(eps=eps) +// is homeomoph to the interior of a simple polygon +function _triangulate(poly, ind, eps=EPSILON, tris=[]) = len(ind)==3 ? _degenerate_tri(select(poly,ind),eps) ? tris // if last 3 pts perform a degenerate triangle, ignore it : concat(tris,[ind]) // otherwise, include it : let( ear = _get_ear(poly,ind,eps) ) - assert( ! error || (ear != undef), +/* +let( x= [if(is_undef(ear)) echo(ind=ind) 0] ) +is_undef(ear) ? tris : +*/ + assert( ear!=undef, "The polygon has twists or all its vertices are collinear or non coplanar.") - ear == undef ? undef : 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 + _triangulate(poly, select(ind,ear[0]+3, ear[0]), eps, tris) // discard it : let( ear_tri = select(ind,ear,ear+2), indr = select(ind,ear+2, ear) // indices of the remaining path ) - _triangulate(poly, indr, error, eps, concat(tris,[ear_tri])); + _triangulate(poly, indr, eps, concat(tris,[ear_tri])); // a returned ear will be: @@ -1832,7 +1825,6 @@ function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) = // 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) = -//let( x = ! is_num(eps) ? echo(ind=ind,eps=eps) 0:0 ) let( lind = len(ind) ) lind==3 ? 0 : let( // the _i-th ear candidate @@ -1849,7 +1841,9 @@ function _get_ear(poly, ind, eps, _i=0) = // otherwise check the next ear candidate _i=0 ) // 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) - && _is_at_left(p2,[p1,next_vert],eps) && _is_at_left(p0,[next_vert,p1],eps) + && ( _is_at_left(p0,[prev_vert,p1],eps) + && _is_at_left(p2,[p1,next_vert],eps) ) ) ) ? false diff --git a/lists.scad b/lists.scad index 70b098f..07356e9 100644 --- a/lists.scad +++ b/lists.scad @@ -230,15 +230,14 @@ function select(list, start, end) = : end==undef ? is_num(start) ? list[ (start%l+l)%l ] - : assert( start==[] || is_vector(start) || is_range(start), "Invalid start parameter") + : assert( is_list(start) || is_range(start), "Invalid start parameter") [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(end), "Invalid end parameter.") let( s = (start%l+l)%l, e = (end%l+l)%l ) (s <= e) - ? [ for (i = [s:1:e]) list[i] ] - : [ for (i = [s:1:l-1]) list[i], - for (i = [0: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]]) ; // Function: slice() @@ -953,7 +952,7 @@ function enumerate(l,idx=undef) = // Example: // l = ["A","B","C","D"]; // echo([for (p=pair(l)) str(p.y,p.x)]); // Outputs: ["BA", "CB", "DC"] -function _old_pair(list, wrap=false) = +function pair(list, wrap=false) = assert(is_list(list)||is_string(list), "Invalid input." ) assert(is_bool(wrap)) let( @@ -962,12 +961,6 @@ function _old_pair(list, wrap=false) = ? [for (i=[0:1:ll-1]) [list[i], list[(i+1) % ll]]] : [for (i=[0:1:ll-2]) [list[i], list[i+1]]]; -function pair(list, wrap=false) = - assert(is_list(list)||is_string(list), "Invalid input." ) - assert(is_bool(wrap)) - let( ll = len(list)-1 ) - [for (i=[0:1:ll-1]) [list[i], list[i+1]], if(wrap) [list[ll], list[0]] ]; - // Function: triplet() // Usage: @@ -988,7 +981,7 @@ function pair(list, wrap=false) = // translate(b) rot(from=FWD,to=v) anchor_arrow2d(); // } // stroke(path); -function _old_triplet(list, wrap=false) = +function triplet(list, wrap=false) = assert(is_list(list)||is_string(list), "Invalid input." ) assert(is_bool(wrap)) let( @@ -998,14 +991,6 @@ function _old_triplet(list, wrap=false) = : [for (i=[0:1:ll-3]) [ list[i], list[i+1], list[i+2] ]]; -function triplet(list, wrap=false) = - assert(is_list(list)||is_string(list), "Invalid input." ) - assert(is_bool(wrap)) - let( ll = len(list) ) - [ for (i=[0:1:ll-3]) [ list[i], list[i+1], list[i+2] ], - each if(wrap) [ [ list[ll-2], list[ll-1], list[0] ], [ list[ll-1], list[0], list[1] ] ] ]; - - // Function: combinations() // Usage: // list = combinations(l, [n]); diff --git a/paths.scad b/paths.scad index 7a59432..fe7fd86 100644 --- a/paths.scad +++ b/paths.scad @@ -274,7 +274,9 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = // 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 // eps and the sign of vals[j]-ref otherwise. - signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ] + signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) vals[j]-ref > eps ? 1 + : vals[j]-ref < -eps ? -1 + : 0] ) 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)]) @@ -284,8 +286,10 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) ) if (isect +// && isect[1]> (i==0 && !closed? -eps: 0) // Apparently too strict && isect[1]>=-eps && isect[1]<= 1+eps +// && isect[2]> 0 && isect[2]>= -eps && isect[2]<= 1+eps) [isect[0], i, isect[1], j, isect[2]] diff --git a/regions.scad b/regions.scad index 5f6d442..662d860 100644 --- a/regions.scad +++ b/regions.scad @@ -216,9 +216,9 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru for(p2=idx(region2)) let( poly = closed2?close_path(region2[p2]):region2[p2], - signs = [for(v=poly*seg_normal) abs(v-ref) < eps ? 0 : sign(v-ref) ] + signs = [for(v=poly*seg_normal) v-ref> eps ? 1 : v-ref<-eps ? -1 : 0] ) - if(max(signs)>=0 && min(signs)<=0) // some edge intersects line [a1,a2] + if(max(signs)>=0 && min(signs)<=0) // some edge edge intersects line [a1,a2] for(j=[0:1:len(poly)-2]) if(signs[j]!=signs[j+1]) let( // exclude non-crossing and collinear segments @@ -260,7 +260,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru // where region1 intersections region2. Split region2 similarly with respect to region1. // The return is a pair of results of the form [split1, split2] where split1=[frags1,frags2,...] // and frags1 is a list of path pieces (in order) from the first path of the region. -// You can pass a single path in for either region, but the output will be a singleton list, as if +// You can pass a single path in for either region, but the output will be a singleton list, as ify // you passed in a singleton region. // Arguments: // region1 = first region diff --git a/vnf.scad b/vnf.scad index dc4dfd6..b2325c1 100644 --- a/vnf.scad +++ b/vnf.scad @@ -414,15 +414,15 @@ function _old_cleave_connected_region(region) = /// Internal Function: _cleave_connected_region(region, eps) /// Description: -/// Given a region that is connected and has its outer border in region[0], -/// produces a overlapping connected path to join internal holes to -/// the outer border without adding points. Output is a single non-simple polygon. -/// It expect that all region paths be simple closed paths with arbitrary widings. -/// The input region paths are also supposed to be disjoint except for common -/// vertices and common edges but with no crossings. It may return `undef` if -/// these conditions are not met. -/// This function implements an extension of the algorithm discussed in: -/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf +/// Given a region that is connected and has its outer border in region[0], +/// produces a polygon with the same points that has overlapping connected paths +/// to join internal holes to the outer border. Output is a single path. +/// It expect that region[0] be a simple closed CW path and that each hole, +/// region[i] for i>0, be a simple closed CCW path. +/// The paths are also supposed to be disjoint except for common vertices and +/// common edges but no crossing. +/// This function implements an extension of the algorithm discussed in: +/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf function _cleave_connected_region(region, eps=EPSILON) = len(region)==1 ? region[0] : let( @@ -439,7 +439,7 @@ function _cleave_connected_region(region, eps=EPSILON) = // 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) +// see _cleave_connected_region(region, eps) function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = let( extr = extremes[n], // @@ -447,18 +447,17 @@ function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = 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 ) -// assert(brdg!=undef, "Error: check input polygon restrictions") - brdg == undef ? undef : + assert(brdg!=undef, "Error: check input polygon restrictions") 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] ] + ? [ 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); @@ -473,12 +472,13 @@ function _bridge(pt, outer,eps) = let( l = len(outer), crxs = - let( edges = pair(outer,wrap=true) ) - [for( i = idx(edges), edge=[edges[i]] ) + [for( i=idx(outer) ) + let( edge = select(outer,i,i+1) ) // 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) ) + if( (edge[0].y> pt.y) + && (edge[1].y<=pt.y) + && ( norm(edge[1]-pt)0 ) ) [ i, // the point of edge with ordinate pt.y abs(pt.y-edge[1].y) pt.y then pt!=vert0 - norm(pt-vert1) < eps ? (proj[0]+1)%l : // if pt touches an outer vertex, return its index - -// assert( ! approx(pt,isect,eps) , // if pt touches the middle of an outer edge -> error -// "A hole vertex is in the middle of an edge" ) - norm(pt-isect) < eps ? undef : // if pt touches the middle of an outer edge -> error + // if pt touches the middle of an outer edge -> error + assert( ! approx(pt,isect,eps) || approx(pt,vert0,eps) || approx(pt,vert1,eps), + "There is a forbidden self_intersection" ) + norm(pt-vert0) < eps ? proj[0] : // if pt touches an outer vertex, return its index + norm(pt-vert1) < eps ? (proj[0]+1)%l : let( // the edge [vert0, vert1] necessarily satisfies vert0.y > vert1.y // indices of candidates to an outer bridge point @@ -558,10 +555,7 @@ function vnf_from_region(region, transform, reverse=false) = regions = region_parts(force_region(region)), vnfs = [ for (rgn = regions) let( - cleaved = path3d(_cleave_connected_region(rgn)) - ) - assert( cleaved, "The region is invalid") - let( + cleaved = path3d(_cleave_connected_region(rgn)), 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]] @@ -674,11 +668,8 @@ function vnf_triangulate(vnf) = let( verts = vnf[0], faces = [for (face=vnf[1]) each len(face)==3 ? [face] : - polygon_triangulate(verts, face)], - invalid = [for(face=faces) if(face==undef) 1 ] - ) - assert( invalid==[], "Some `vnf` face cannot be triangulated.") - [verts, faces]; + polygon_triangulate(verts, face)] + ) [verts, faces]; From 6bd1dd918f231181bcc670462958a3e817fb7dcd Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Thu, 4 Nov 2021 12:09:29 +0000 Subject: [PATCH 3/5] Correction in _cleave_connected_region and polygon_triangulate and some few tweaks --- comparisons.scad | 4 +-- geometry.scad | 88 +++++++++++++++++++++++---------------------- lists.scad | 7 ++-- paths.scad | 5 ++- regions.scad | 6 ++-- vnf.scad | 93 +++++++++++++++++++++++++----------------------- 6 files changed, 105 insertions(+), 98 deletions(-) diff --git a/comparisons.scad b/comparisons.scad index bd2a25c..300b468 100644 --- a/comparisons.scad +++ b/comparisons.scad @@ -290,12 +290,12 @@ function compare_lists(a, b) = // idx = find_approx(val, list, [start=], [eps=]); // indices = find_approx(val, list, all=true, [start=], [eps=]); // 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: // 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. // --- -// 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. // eps = The maximum allowed floating point rounding error for numeric comparisons. function find_approx(val, list, start=0, all=false, eps=EPSILON) = diff --git a/geometry.scad b/geometry.scad index c94d5ca..77930bb 100644 --- a/geometry.scad +++ b/geometry.scad @@ -104,15 +104,16 @@ function _tri_class(tri, eps=EPSILON) = /// class = _pt_in_tri(point, tri); /// Topics: Geometry, Points, Triangles /// Description: -/// Return 1 if point is inside the triangle interion. -/// Return =0 if point is on the triangle border. -/// Return -1 if point is outside the triangle. +// 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), + 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) ); @@ -1701,7 +1702,7 @@ 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, 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` // are considered as vertices of the polygon. // . @@ -1710,47 +1711,49 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // 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 +// A polygon is non-twisted iff it is simple or it has a partition in // simple polygons with the same winding such that the intersection of any two partitions is -// made of full edges of both partitions. These polygons may have "touching" vertices +// 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 // (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 (polygons with -// zero area), 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. // . -// Twisted polygons have no consistent winding and when input to this function usually produce -// an error but when an error is not issued the outputs are not correct triangulations. The function +// Twisted polygons have no consistent winding and when input to this function usually reports +// 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 -// issue an error for this case. +// report an error for this case. // Arguments: // poly = Array of the polygon vertices. // ind = A list indexing the vertices of the polygon in `poly`. // 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); // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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] ]; // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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] ]; // tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); // 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], // [7,7], [7,3], [3,3] ]; -// tris = polygon_triangulate(poly); // see from above +// tris = polygon_triangulate(poly); // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); @@ -1762,19 +1765,18 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ]; // color("blue") // 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_undef(ind) - || (is_vector(ind) && min(ind)>=0 && max(ind)=0 && max(ind) 2*eps, - "The polygon vertices are collinear.") - [ind] + let( degen = norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) < 2*eps ) + assert( ! error || ! degen, "The polygon vertices are collinear.") + degen ? undef : [ind] : len(poly[ind[0]]) == 3 ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane let( @@ -1785,44 +1787,46 @@ function polygon_triangulate(poly, ind, eps=EPSILON) = pts = select(poly,ind), nrm = -polygon_normal(pts) ) - assert( nrm!=undef, + assert( ! error || (nrm != undef), "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.") + nrm == undef ? undef : 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]) // 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) ) + tris == undef ? undef : + [for(tri=tris) select(ind,tri) ] : is_polygon_clockwise(select(poly, ind)) - ? _triangulate( poly, ind, eps ) - : [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ]; + ? _triangulate( poly, ind, error, eps ) + : let( tris = _triangulate( poly, reverse(ind), error, eps ) ) + tris == undef ? undef : + [for(tri=tris) reverse(tri) ]; // 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 (tecnically) the ones whose interior -// is homeomoph to the interior of a simple polygon -function _triangulate(poly, ind, eps=EPSILON, tris=[]) = +// the polygons accepted by this function are those decomposable in simple +// CW polygons. +function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) = len(ind)==3 ? _degenerate_tri(select(poly,ind),eps) ? tris // if last 3 pts perform a degenerate triangle, ignore it : concat(tris,[ind]) // otherwise, include it : let( ear = _get_ear(poly,ind,eps) ) -/* -let( x= [if(is_undef(ear)) echo(ind=ind) 0] ) -is_undef(ear) ? tris : -*/ - assert( ear!=undef, + assert( ! error || (ear != undef), "The polygon has twists or all its vertices are collinear or non coplanar.") + ear == undef ? undef : is_list(ear) // is it a degenerate ear ? ? len(ind) <= 4 ? tris : - _triangulate(poly, select(ind,ear[0]+3, ear[0]), eps, tris) // discard it + _triangulate(poly, select(ind,ear[0]+3, ear[0]), error, eps, tris) // discard it : let( ear_tri = select(ind,ear,ear+2), 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: @@ -1847,9 +1851,7 @@ function _get_ear(poly, ind, eps, _i=0) = // otherwise check the next ear candidate _i=0 ) // 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,next_vert],eps) ) + && _is_at_left(p0,[prev_vert,p1],eps) && _is_at_left(p2,[p1,prev_vert],eps) + && _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); + : _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1); // Function: is_polygon_clockwise() diff --git a/lists.scad b/lists.scad index 9da1972..0720545 100644 --- a/lists.scad +++ b/lists.scad @@ -230,14 +230,15 @@ function select(list, start, end) = : end==undef ? is_num(start) ? 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 ] ] : assert(is_finite(start), "When `end` is given, `start` parameter should be a number.") assert(is_finite(end), "Invalid end parameter.") let( s = (start%l+l)%l, e = (end%l+l)%l ) (s <= e) - ? [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:e]) list[i] ] + : [ for (i = [s:1:l-1]) list[i], + for (i = [0:1:e]) list[i] ] ; // Function: slice() diff --git a/paths.scad b/paths.scad index 15349dd..1c76251 100644 --- a/paths.scad +++ b/paths.scad @@ -278,9 +278,8 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = // 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 // 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 - : vals[j]-ref < -eps ? -1 - : 0] + signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) + abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ] ) 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)]) diff --git a/regions.scad b/regions.scad index 36549b3..05bdd57 100644 --- a/regions.scad +++ b/regions.scad @@ -282,9 +282,9 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru for(p2=idx(region2)) let( 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]) if(signs[j]!=signs[j+1]) let( // exclude non-crossing and collinear segments @@ -329,7 +329,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru // where region1 intersections region2. Split region2 similarly with respect to region1. // The return is a pair of results of the form [split1, split2] where split1=[frags1,frags2,...] // and frags1 is a list of path pieces (in order) from the first path of the region. -// You can pass a single path in for either region, but the output will be a singleton list, as ify +// You can pass a single path in for either region, but the output will be a singleton list, as if // you passed in a singleton region. // Arguments: // region1 = first region diff --git a/vnf.scad b/vnf.scad index 3457a45..3351bd9 100644 --- a/vnf.scad +++ b/vnf.scad @@ -414,22 +414,22 @@ function _old_cleave_connected_region(region) = /// Internal Function: _cleave_connected_region(region, eps) /// Description: -/// Given a region that is connected and has its outer border in region[0], -/// produces a polygon with the same points that has overlapping connected paths -/// to join internal holes to the outer border. Output is a single path. -/// It expect that region[0] be a simple closed CW path and that each hole, -/// region[i] for i>0, be a simple closed CCW path. -/// The paths are also supposed to be disjoint except for common vertices and -/// common edges but no crossing. -/// This function implements an extension of the algorithm discussed in: -/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf +/// Given a region that is connected and has its outer border in region[0], +/// produces a overlapping connected path to join internal holes to +/// the outer border without adding points. Output is a single non-simple polygon. +/// Requirements: +/// It expects that all region paths be simple closed paths, with region[0] CW and +/// the other paths CCW and encircled by region[0]. The input region paths are also +/// supposed to be disjoint except for common vertices and common edges but with +/// no crossings. It may return `undef` if these conditions are not met. +/// This function implements an extension of the algorithm discussed in: +/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf function _cleave_connected_region(region, eps=EPSILON) = len(region)==1 ? region[0] : let( - outer = deduplicate(clockwise_polygon(region[0])), // - holes = [for(i=[1:1:len(region)-1]) // possibly unneeded - let(poly=region[i]) // - deduplicate( ccw_polygon(poly) ) ], // + outer = deduplicate(region[0]), // + holes = [for(i=[1:1:len(region)-1]) // deduplication possibly unneeded + deduplicate( region[i] ) ], // extridx = [for(li=holes) max_index(column(li,0)) ], // the right extreme vertex for each hole sorted by decreasing x values extremes = sort( [for(i=idx(holes)) [ i, extridx[i], -holes[i][extridx[i]].x] ], idx=2 ) @@ -439,7 +439,7 @@ function _cleave_connected_region(region, eps=EPSILON) = // 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) +// see: _cleave_connected_region(region, eps) function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = let( extr = extremes[n], // @@ -447,17 +447,17 @@ function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = 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 ) - assert(brdg!=undef, "Error: check input polygon restrictions") + 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] ] + ? [ 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); @@ -472,13 +472,13 @@ function _bridge(pt, outer,eps) = let( l = len(outer), crxs = - [for( i=idx(outer) ) - let( edge = select(outer,i,i+1) ) + 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) - && (edge[1].y<=pt.y) - && ( norm(edge[1]-pt)0 ) ) + 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) error - assert( ! approx(pt,isect,eps) || approx(pt,vert0,eps) || approx(pt,vert1,eps), - "There is a forbidden self_intersection" ) - norm(pt-vert0) < eps ? proj[0] : // if pt touches an outer vertex, return its index - norm(pt-vert1) < eps ? (proj[0]+1)%l : + 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 @@ -553,13 +552,15 @@ function _bridge(pt, outer,eps) = function vnf_from_region(region, transform, reverse=false) = let ( regions = region_parts(force_region(region)), - vnfs = [ - for (rgn = regions) let( - cleaved = path3d(_cleave_connected_region(rgn)), - 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]] - ], + vnfs = + [ for (rgn = regions) + let( cleaved = path3d(_cleave_connected_region(rgn)) ) + assert( cleaved, "The region is invalid") + let( + 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) ) vnf_triangulate(outvnf); @@ -667,9 +668,13 @@ function _link_indicator(l,imin,imax) = function vnf_triangulate(vnf) = let( verts = vnf[0], - faces = [for (face=vnf[1]) each len(face)==3 ? [face] : - polygon_triangulate(verts, face)] - ) [verts, faces]; + faces = [for (face=vnf[1]) + each (len(face)==3 ? [face] : + let( tris = polygon_triangulate(verts, face) ) + assert( tris!=undef, "Some `vnf` face cannot be triangulated.") + tris ) ] + ) + [verts, faces]; From e60f5fd932c6ecaf237155cc56a19fbd700c8ed3 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Thu, 4 Nov 2021 12:12:05 +0000 Subject: [PATCH 4/5] Update paths.scad --- paths.scad | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/paths.scad b/paths.scad index 1c76251..31b55b6 100644 --- a/paths.scad +++ b/paths.scad @@ -279,7 +279,7 @@ function _path_self_intersections(path, closed=true, eps=EPSILON) = // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than // eps and the sign of vals[j]-ref otherwise. signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) - abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ] + abs(vals[j]-ref) < eps ? 0 : sign(vals[j]-ref) ] ) 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)]) From abab41a2dbe37b383f3b293bf17dc4f77394956f Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Thu, 4 Nov 2021 19:27:14 +0000 Subject: [PATCH 5/5] Doc update of polygon_triangulate --- geometry.scad | 5 +++-- vnf.scad | 28 ---------------------------- 2 files changed, 3 insertions(+), 30 deletions(-) diff --git a/geometry.scad b/geometry.scad index 77930bb..107cf03 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1698,7 +1698,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // Function: polygon_triangulate() // Usage: -// triangles = polygon_triangulate(poly, [ind], [eps]) +// triangles = polygon_triangulate(poly, [ind], [error], [eps]) // 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 @@ -1728,6 +1728,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // Arguments: // poly = Array of the polygon vertices. // 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 // Example(2D,NoAxes): a simple polygon; see from above // poly = star(id=10, od=15,n=11); @@ -1767,7 +1768,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // vnf_wireframe(vnf_tri, width=.15); 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_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)=0 && max(ind)0 && i!=idxi+1) - region[i] - ] - ) - assert(len(orgn)