From c5da41ed8ca653a15bf7ff7bf6de4b0b25d4c5b0 Mon Sep 17 00:00:00 2001
From: RonaldoCMP <rcmpersiano@gmail.com>
Date: Mon, 1 Nov 2021 04:42:02 +0000
Subject: [PATCH 1/8] Correct bugs in polygon_triangulation and
 _cleave_connected_region(

---
 geometry.scad            | 213 +++++++++++++++++++++++++++------------
 tests/test_all.scad      |  32 ++++++
 tests/test_geometry.scad |   8 +-
 vectors.scad             |   9 +-
 vnf.scad                 | 125 ++++++++++++++++++++++-
 5 files changed, 307 insertions(+), 80 deletions(-)
 create mode 100644 tests/test_all.scad

diff --git a/geometry.scad b/geometry.scad
index 37326e1..c09604b 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -38,14 +38,10 @@ function _is_point_on_line(point, line, bounded=false, eps=EPSILON) =
         t  = v0*v1/(v1*v1),
         bounded = force_list(bounded,2)
     ) 
-    abs(cross(v0,v1))<eps*norm(v1) 
+    abs(cross(v0,v1))<=eps*norm(v1) 
     && (!bounded[0] || t>=-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
 ///_dist2line works for any dimension
@@ -61,7 +57,67 @@ function _valid_line(line,dim,eps=EPSILON) =
 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:
+///   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), 
+          _tri_class([tri[1],tri[2],point],eps), 
+          _tri_class([tri[2],tri[0],point],eps) );
+        
+
+/// Internal Function: _point_left_of_line2d()
 /// Usage:
 ///   pt = point_left_of_line2d(point, line);
 /// Topics: Geometry, Points, Lines
@@ -72,10 +128,11 @@ function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps
 /// Arguments:
 ///   point = The point to check position of.
 ///   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." )
-    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()
 // Usage:
@@ -1648,18 +1705,19 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
 //   .
 //   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 
+//   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 
 //   (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, 
-//   it returns an empty list for 2d polygons and issues an error for 3d polygons. 
+//   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. 
 //   .
-//   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
+//   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 
 //   issue an error for this case. 
 // 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`.
 //   eps = A maximum tolerance in geometrical tests. Default: EPSILON
 // Example(2D,NoAxes):
@@ -1686,7 +1744,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
 // 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
+//   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); 
@@ -1703,97 +1761,118 @@ function polygon_triangulate(poly, ind, eps=EPSILON) =
     assert(is_undef(ind) 
            || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
            "Improper or out of bounds list of indices")
+    (! is_undef(ind) ) && len(ind) == 0 ? [] :
     let( ind = is_undef(ind) ? count(len(poly)) : ind )
     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
         assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 2*eps,
                 "The polygon vertices are collinear.") 
         [ind]
       : 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( 
                 ind = deduplicate_indexed(poly, ind, eps) 
             )
             len(ind)<3 ? [] :
             let(
                 pts = select(poly,ind),
-                nrm = polygon_normal(pts)
+                nrm = -polygon_normal(pts)
             )
             assert( 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.") 
             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])
+                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( cw = is_polygon_clockwise(select(poly, ind)) )
-            cw 
-              ? [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ]
-              : _triangulate( poly, ind, eps );
+          : is_polygon_clockwise(select(poly, ind)) 
+              ? _triangulate( poly, ind,  eps )
+              : [for(tri=_triangulate( poly, reverse(ind), eps )) 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 (tecnically) the ones whose interior
+// is homeomoph to the interior of a simple polygon
+function _triangulate(poly, ind,  eps=EPSILON, tris=[]) =
     len(ind)==3 
-    ?   _is_degenerate(select(poly,ind),eps)  
-        ?   tris // last 3 pts perform a degenerate triangle, ignore it
+    ?   _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, 
-                "The polygon has self-intersections or its vertices are collinear or non coplanar.") 
-        is_list(ear) // degenerate ear
-        ?   _triangulate(poly, select(ind,ear[0]+2, ear[0]), eps, tris) // discard it
+            "The polygon has twists or all its vertices are collinear or non coplanar.") 
+        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
         :   let(
                 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,  eps, concat(tris,[ear_tri]));
 
 
 // a returned ear will be:
-// 1. a CCW (non-degenerate) triangle, made of subsequent vertices, without other 
-//    points inside except possibly at its vertices
+// 1. a CW non-reflex triangle, made of subsequent poly vertices, without any other 
+//    poly points inside except possibly at its own 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) =
-    _i>=len(ind) ? undef : // poly has no ears
+function _get_ear(poly, ind,  eps, _i=0) =
+    let( lind = len(ind) )
+    lind==3 ? 0 :
     let( // the _i-th ear candidate
         p0 = poly[ind[_i]],
-        p1 = poly[ind[(_i+1)%len(ind)]],
-        p2 = poly[ind[(_i+2)%len(ind)]]
+        p1 = poly[ind[(_i+1)%lind]],
+        p2 = poly[ind[(_i+2)%lind]]
     )
-    // 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),
-        q      = [(p0-p2).y, (p2-p0).x],  // orthogonal to ray [p0,p2] pointing right
-        r      = [(p2-p1).y, (p1-p2).x],  // orthogonal to ray [p2,p1] pointing right
-        s      = [(p1-p0).y, (p0-p1).x],  // orthogonal to ray [p1,p0] pointing right
-        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 ]                       
+    // if vertex p1 is a convex candidate to be an ear,
+    // check if the triangle [p0,p1,p2] contains any other point
+    // except possibly p0 and p2
+    // exclude the ear candidate central vertex p1 from the verts to check 
+    _tri_class([p0,p1,p2],eps) > 0  
+    &&  _none_inside(select(ind,_i+2, _i),poly,p0,p1,p2,eps) ? _i : // found an ear
+    // otherwise check the next ear candidate 
+    _i<lind-1 ?  _get_ear(poly, ind,  eps, _i=_i+1) :
+    // poly has no ears, look for wiskers
+    let( 
+        wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ]
     )
-    inside==[] ? _i : // found an ear
-    // check the next ear candidate
-    _get_ear(poly, ind, eps, _i=_i+1);
-
-
-// true for some specific kinds of degeneracy
-function _is_degenerate(tri,eps) =
-       norm(tri[0]-tri[1])<eps || norm(tri[1]-tri[2])<eps || norm(tri[2]-tri[0])<eps ;
-
-
-function _is_cw2(a,b,c,eps=EPSILON) = cross(a-c,b-c)<eps*norm(a-c)*norm(b-c);
+    wiskers==[] ? undef : [wiskers[0]];
+    
+    
 
+// returns false ASA it finds some reflex vertex of poly[idxs[.]] 
+// inside the triangle different from p0 and p2
+// 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)]
+    )
+    // check if vert prevent [p0,p1,p2] to be an ear
+    // this conditions might have a simpler expression
+    _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 && 
+            _tri_class([p1,p2,vert],eps)>0 && 
+            _tri_class([p2,p0,vert],eps)>=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) ) 
+              ) 
+        )
+    ?   false
+    :    _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1);
 
 
 // Function: is_polygon_clockwise()
diff --git a/tests/test_all.scad b/tests/test_all.scad
new file mode 100644
index 0000000..e51434a
--- /dev/null
+++ b/tests/test_all.scad
@@ -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>
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 7eff2bc..1b649d1 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -84,14 +84,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] ];
 		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));
+		tris0 = (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]])));
+		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]])));
+		assert(approx(tris2,( [[3, 4, 5], [1, 2, 3]])));
 		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(){
diff --git a/vectors.scad b/vectors.scad
index 7343941..687f0d4 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -491,12 +491,11 @@ function _bt_tree(points, ind, leafsize=25) =
         bounds = pointlist_bounds(select(points,ind)),
         coord  = max_index(bounds[1]-bounds[0]), 
         projc  = [for(i=ind) points[i][coord] ],
-        pmc    = mean(projc), 
-        pivot  = min_index([for(p=projc) abs(p-pmc)]),
+        meanpr = mean(projc), 
+        pivot  = min_index([for(p=projc) abs(p-meanpr)]),
         radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]),
-        median = median(projc),
-        Lind   = [for(i=idx(ind)) if(projc[i]<=median && i!=pivot) ind[i] ],
-        Rind   = [for(i=idx(ind)) if(projc[i] >median && i!=pivot) ind[i] ]
+        Lind   = [for(i=idx(ind)) if(projc[i]<=meanpr && i!=pivot) ind[i] ],
+        Rind   = [for(i=idx(ind)) if(projc[i] >meanpr && i!=pivot) ind[i] ]
       )
     [ ind[pivot], radius, _bt_tree(points, Lind, leafsize), _bt_tree(points, Rind, leafsize) ];
 
diff --git a/vnf.scad b/vnf.scad
index dafad7f..b2325c1 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -318,14 +318,13 @@ function vnf_merge(vnfs, cleanup=false, eps=EPSILON) =
     cleanup? _vnf_cleanup(verts,faces,eps) : [verts,faces];
 
 
-
 function _vnf_cleanup(verts,faces,eps) = 
     let(
         dedup  = vector_search(verts,eps,verts),                 // collect vertex duplicates
         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 
         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 = 
             [ for(face=faces) 
                 let(
@@ -388,7 +387,7 @@ function _join_paths_at_vertices(path1,path2,v1,v2) =
 // 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.  
-function _cleave_connected_region(region) =
+function _old_cleave_connected_region(region) =
     len(region)==0? [] :
     len(region)<=1? clockwise_polygon(region[0]) :
     let(
@@ -411,8 +410,126 @@ function _cleave_connected_region(region) =
         ]
     )
     assert(len(orgn)<len(region))
-    _cleave_connected_region(orgn);
+    _old_cleave_connected_region(orgn);
 
+/// 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
+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) ) ],     //
+        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 )
+    ) 
+    _polyHoles(outer, holes, extremes, eps, 0);
+
+
+// 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
+    )
+    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] ]
+    )
+    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 = 
+            [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) 
+                    && (edge[1].y<=pt.y) 
+                    && ( norm(edge[1]-pt)<eps // accepts touching vertices
+                         || _tri_class([pt, edge[0], edge[1]], eps)>0 ) )
+                    [ 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]
+                    ]
+             ]
+    )
+    assert(crxs!=[], "Error: check input polygon restrictions") 
+    let( 
+        // the intersection point nearest to pt
+        minX    = min([for(p=crxs) p[1].x]),
+        crxcand = [for(crx=crxs) if(crx[1].x < minX+eps) crx ],
+        nearest = min_index([for(crx=crxcand) outer[crx[0]].y]),
+        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
+    )
+    // 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
+        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()

From 8c4022fa3c98e19fa5c5f3d660c14e55448bd24d Mon Sep 17 00:00:00 2001
From: RonaldoCMP <rcmpersiano@gmail.com>
Date: Wed, 3 Nov 2021 20:57:35 +0000
Subject: [PATCH 2/8] 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)<len(poly) ),
+    assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
            "Improper or out of bounds list of indices")
-    (! is_undef(ind) ) && len(ind) == 0 ? [] :
     let( ind = is_undef(ind) ? count(len(poly)) : ind )
+    len(ind) <=2 ? [] :
     len(ind) == 3 
       ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : 
         // non zero area
-        assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 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<lind-1 ?  _get_ear(poly, ind,  eps, _i=_i+1) :
     // poly has no ears, look for wiskers
-    let( 
-        wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ]
-    )
+    let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
     wiskers==[] ? undef : [wiskers[0]];
     
     
@@ -1867,8 +1873,8 @@ function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
             _tri_class([p2,p0,vert],eps)>=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)<eps // accepts touching vertices
-                         || _tri_class([pt, edge[0], edge[1]], eps)>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)<eps ? edge[1] :
@@ -487,22 +487,25 @@ function _bridge(pt, outer,eps) =
                     ]
              ]
     )
-    assert(crxs!=[], "Error: check input polygon restrictions") 
+//    assert(crxs!=[], "Error: check input polygon restrictions") 
+    crxs == [] ? undef :
     let( 
-        // the intersection point nearest to pt
+        // 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 = min_index([for(crx=crxcand) outer[crx[0]].y]),
+        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
     )
-    // 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 :
+    // 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 7ed34f592a46fb4e5f791d55cbab645850bf1daa Mon Sep 17 00:00:00 2001
From: Garth Minette <revarbat@gmail.com>
Date: Wed, 3 Nov 2021 18:52:51 -0700
Subject: [PATCH 3/8] Speed optimizations for approx()

---
 comparisons.scad | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/comparisons.scad b/comparisons.scad
index bd2a25c..3367de1 100644
--- a/comparisons.scad
+++ b/comparisons.scad
@@ -12,7 +12,8 @@
 // Usage:
 //   test = approx(a, b, [eps])
 // 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:
 //   a = First value.
 //   b = Second value.
@@ -21,12 +22,16 @@
 //   test1 = approx(-0.3333333333,-1/3);  // Returns: true
 //   test2 = approx(0.3333333333,1/3);    // Returns: true
 //   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
+//   test6 = approx([0,0,sin(45)],[0,0,sqrt(2)/2]);  // Returns: true
 function approx(a,b,eps=EPSILON) = 
-    (a==b && is_bool(a) == is_bool(b)) ||
-    (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]);
+    a == b? is_bool(a) == is_bool(b) :
+    is_num(a) && is_num(b)? abs(a-b) <= eps :
+    is_list(a) && is_list(b) && len(a) == len(b)? (
+        is_num(a.x) && is_num(b.x)? norm(a-b) <= eps :
+        [] == [for (i=idx(a)) if (!approx(a[i],b[i],eps=eps)) 1]
+    ) : false;
 
 
 // Function: all_zero()

From 36968a463b638c34c9d5f2a9a48e9cf499515907 Mon Sep 17 00:00:00 2001
From: RonaldoCMP <rcmpersiano@gmail.com>
Date: Thu, 4 Nov 2021 07:56:34 +0000
Subject: [PATCH 4/8] 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)<len(poly) ),
+    assert(is_undef(ind) 
+           || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
            "Improper or out of bounds list of indices")
+    (! is_undef(ind) ) && len(ind) == 0 ? [] :
     let( ind = is_undef(ind) ? count(len(poly)) : ind )
-    len(ind) <=2 ? [] :
     len(ind) == 3 
       ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : 
         // non zero area
-        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 :
+        assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 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<lind-1 ?  _get_ear(poly, ind,  eps, _i=_i+1) :
     // poly has no ears, look for wiskers
-    let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
+    let( 
+        wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ]
+    )
     wiskers==[] ? undef : [wiskers[0]];
     
     
@@ -1873,8 +1867,8 @@ function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
             _tri_class([p2,p0,vert],eps)>=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)<eps // accepts touching vertices
+                         || _tri_class([pt, edge[0], edge[1]], eps)>0 ) )
                     [ i,
                       // the point of edge with ordinate pt.y
                       abs(pt.y-edge[1].y)<eps ? edge[1] :
@@ -487,25 +487,22 @@ function _bridge(pt, outer,eps) =
                     ]
              ]
     )
-//    assert(crxs!=[], "Error: check input polygon restrictions") 
-    crxs == [] ? undef :
+    assert(crxs!=[], "Error: check input polygon restrictions") 
     let( 
-        // the intersection point of the nearest edge to pt with minimum slope
+        // the intersection point nearest to pt
         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
+        crxcand = [for(crx=crxs) if(crx[1].x < minX+eps) crx ],
+        nearest = min_index([for(crx=crxcand) outer[crx[0]].y]),
         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
     )
-    // 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
+    // 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 <rcmpersiano@gmail.com>
Date: Thu, 4 Nov 2021 12:09:29 +0000
Subject: [PATCH 5/8] 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)<len(poly) ),
+   assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
            "Improper or out of bounds list of indices")
-    (! is_undef(ind) ) && len(ind) == 0 ? [] :
     let( ind = is_undef(ind) ? count(len(poly)) : ind )
+    len(ind) <=2 ? [] :
     len(ind) == 3 
       ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : 
         // non zero area
-        assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 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<lind-1 ?  _get_ear(poly, ind,  eps, _i=_i+1) :
     // poly has no ears, look for wiskers
-    let( 
-        wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ]
-    )
+    let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
     wiskers==[] ? undef : [wiskers[0]];
     
     
@@ -1873,12 +1875,12 @@ function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
             _tri_class([p2,p0,vert],eps)>=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)<eps // accepts touching vertices
-                         || _tri_class([pt, edge[0], edge[1]], eps)>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)<eps ? edge[1] :
@@ -487,22 +487,21 @@ function _bridge(pt, outer,eps) =
                     ]
              ]
     )
-    assert(crxs!=[], "Error: check input polygon restrictions") 
+    crxs == [] ? undef :
     let( 
-        // the intersection point nearest to pt
+        // 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 = min_index([for(crx=crxcand) outer[crx[0]].y]),
+        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
     )
-    // 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 :
+    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 <rcmpersiano@gmail.com>
Date: Thu, 4 Nov 2021 12:12:05 +0000
Subject: [PATCH 6/8] 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 <rcmpersiano@gmail.com>
Date: Thu, 4 Nov 2021 19:27:14 +0000
Subject: [PATCH 7/8] 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)<len(poly) ),
+    assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
            "Improper or out of bounds list of indices")
     let( ind = is_undef(ind) ? count(len(poly)) : ind )
     len(ind) <=2 ? [] :
diff --git a/vnf.scad b/vnf.scad
index 3351bd9..0e7684f 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -384,34 +384,6 @@ function _join_paths_at_vertices(path1,path2,v1,v2) =
     ];                      
 
         
-// 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.  
-function _old_cleave_connected_region(region) =
-    len(region)==0? [] :
-    len(region)<=1? clockwise_polygon(region[0]) :
-    let(
-        dists = [
-            for (i=[1:1:len(region)-1])
-            _path_path_closest_vertices(region[0],region[i])
-        ],
-        idxi = min_index(column(dists,0)),
-        newoline = _join_paths_at_vertices(
-            region[0], region[idxi+1],
-            dists[idxi][1], dists[idxi][2]
-        )
-    ) len(region)==2? clockwise_polygon(newoline) :
-    let(
-        orgn = [
-            newoline,
-            for (i=idx(region))
-                if (i>0 && i!=idxi+1)
-                    region[i]
-        ]
-    )
-    assert(len(orgn)<len(region))
-    _old_cleave_connected_region(orgn);
-
 /// Internal Function: _cleave_connected_region(region, eps)
 /// Description:
 ///   Given a region that is connected and has its outer border in region[0],

From 2b7d0d935dba87d3420578fee35f82702145c2ef Mon Sep 17 00:00:00 2001
From: Garth Minette <revarbat@gmail.com>
Date: Thu, 4 Nov 2021 21:13:50 -0700
Subject: [PATCH 8/8] Fixes for approx() with mixed types.

---
 comparisons.scad | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/comparisons.scad b/comparisons.scad
index 3367de1..906aa41 100644
--- a/comparisons.scad
+++ b/comparisons.scad
@@ -29,8 +29,14 @@ function approx(a,b,eps=EPSILON) =
     a == b? is_bool(a) == is_bool(b) :
     is_num(a) && is_num(b)? abs(a-b) <= eps :
     is_list(a) && is_list(b) && len(a) == len(b)? (
-        is_num(a.x) && is_num(b.x)? norm(a-b) <= eps :
-        [] == [for (i=idx(a)) if (!approx(a[i],b[i],eps=eps)) 1]
+        [] == [
+            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;
 
 
@@ -783,3 +789,4 @@ function list_smallest(list, k) =
     let( bigger  = [for(li=list) if(li>v) li ] )
     concat(smaller, equal, list_smallest(bigger, k-len(smaller) -len(equal)));
 
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap