1//////////////////////////////////////////////////////////////////////
   2// LibFile: vnf.scad
   3//   The Vertices'N'Faces structure (VNF) holds the data used by polyhedron() to construct objects: a vertex
   4//   list and a list of faces.  This library makes it easier to construct polyhedra by providing
   5//   functions to construct, merge, and modify VNF data, while avoiding common pitfalls such as
   6//   reversed faces.  It can find faults in your polyhedrons.  Note that this file is for low level manipulation
   7//   of lists of vertices and faces: it can perform some simple transformations on VNF structures
   8//   but cannot perform boolean operations on the polyhedrons represented by VNFs.
   9// Includes:
  10//   include <BOSL2/std.scad>
  11// FileGroup: Advanced Modeling
  12// FileSummary: Vertices 'n' Faces structure.  Makes polyhedron() easier to use.
  13// FileFootnotes: STD=Included in std.scad
  14//////////////////////////////////////////////////////////////////////
  15
  16
  17// Section: Creating Polyhedrons with VNF Structures
  18//   VNF stands for "Vertices'N'Faces".  VNF structures are 2-item lists, `[VERTICES,FACES]` where the
  19//   first item is a list of vertex points, and the second is a list of face indices into the vertex
  20//   list.  Each VNF is self contained, with face indices referring only to its own vertex list.
  21//   You can construct a `polyhedron()` in parts by describing each part in a self-contained VNF, then
  22//   merge the various VNFs to get the completed polyhedron vertex list and faces.
  23
  24/// Constant: EMPTY_VNF
  25/// Description:
  26///   The empty VNF data structure.  Equal to `[[],[]]`.
  27EMPTY_VNF = [[],[]];  // The standard empty VNF with no vertices or faces.
  28
  29
  30// Function: vnf_vertex_array()
  31// Synopsis: Returns a VNF structure from a rectangular vertex list.
  32// SynTags: VNF
  33// Topics: VNF Generators, Lists
  34// See Also: vnf_tri_array(), vnf_join(), vnf_from_polygons(), vnf_from_region()
  35// Usage:
  36//   vnf = vnf_vertex_array(points, [caps=], [cap1=], [cap2=], [style=], [reverse=], [col_wrap=], [row_wrap=], [triangulate=]);
  37// Description:
  38//   Creates a VNF structure from a rectangular vertex list, creating edges that connect the adjacent vertices in the vertex list
  39//   and creating the faces defined by those edges.  You can optionally create the edges and faces to wrap the last column
  40//   back to the first column, or wrap the last row to the first.  Endcaps can be added to either
  41//   the first and/or last rows.  The style parameter determines how the quadrilaterals are divided into
  42//   triangles.  The default style is an arbitrary, systematic subdivision in the same direction.  The "alt" style
  43//   is the uniform subdivision in the other (alternate) direction.  The "flip1" style is an arbitrary division which alternates the
  44//   direction for any adjacent pair of quadrilaterals.  The "flip2" style is the alternating division that is the opposite of "flip1". 
  45//   The "min_edge" style picks the shorter edge to
  46//   subdivide for each quadrilateral, so the division may not be uniform across the shape.  The "quincunx" style
  47//   adds a vertex in the center of each quadrilateral and creates four triangles, and the "convex" and "concave" styles
  48//   choose the locally convex/concave subdivision.  The "min_area" option creates the triangulation with the minimal area.  Degenerate faces
  49//   are not included in the output, but if this results in unused vertices they will still appear in the output.
  50// Arguments:
  51//   points = A list of vertices to divide into columns and rows.
  52//   ---
  53//   caps = If true, add endcap faces to the first AND last rows.
  54//   cap1 = If true, add an endcap face to the first row.
  55//   cap2 = If true, add an endcap face to the last row.
  56//   col_wrap = If true, add faces to connect the last column to the first.
  57//   row_wrap = If true, add faces to connect the last row to the first.
  58//   reverse = If true, reverse all face normals.
  59//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", "flip1", "flip2",  "min_edge", "min_area", "quincunx", "convex" and "concave".
  60//   triangulate = If true, triangulates endcaps to resolve possible CGAL issues.  This can be an expensive operation if the endcaps are complex.  Default: false
  61// Example(3D):
  62//   vnf = vnf_vertex_array(
  63//       points=[
  64//           for (h = [0:5:180-EPSILON]) [
  65//               for (t = [0:5:360-EPSILON])
  66//                   cylindrical_to_xyz(100 + 12 * cos((h/2 + t)*6), t, h)
  67//           ]
  68//       ],
  69//       col_wrap=true, caps=true, reverse=true, style="alt"
  70//   );
  71//   vnf_polyhedron(vnf);
  72// Example(3D): Both `col_wrap` and `row_wrap` are true to make a torus.
  73//   vnf = vnf_vertex_array(
  74//       points=[
  75//           for (a=[0:5:360-EPSILON])
  76//               apply(
  77//                   zrot(a) * right(30) * xrot(90),
  78//                   path3d(circle(d=20))
  79//               )
  80//       ],
  81//       col_wrap=true, row_wrap=true, reverse=true
  82//   );
  83//   vnf_polyhedron(vnf);
  84// Example(3D): Möbius Strip.  Note that `row_wrap` is not used, and the first and last profile copies are the same.
  85//   vnf = vnf_vertex_array(
  86//       points=[
  87//           for (a=[0:5:360]) apply(
  88//               zrot(a) * right(30) * xrot(90) * zrot(a/2+60),
  89//               path3d(square([1,10], center=true))
  90//           )
  91//       ],
  92//       col_wrap=true, reverse=true
  93//   );
  94//   vnf_polyhedron(vnf);
  95// Example(3D): Assembling a Polyhedron from Multiple Parts
  96//   wall_points = [
  97//       for (a = [-90:2:90]) apply(
  98//           up(a) * scale([1-0.1*cos(a*6),1-0.1*cos((a+90)*6),1]),
  99//           path3d(circle(d=100))
 100//       )
 101//   ];
 102//   cap = [
 103//       for (a = [0:0.01:1+EPSILON]) apply(
 104//           up(90-5*sin(a*360*2)) * scale([a,a,1]),
 105//           wall_points[0]
 106//       )
 107//   ];
 108//   cap1 = [for (p=cap) down(90, p=zscale(-1, p=p))];
 109//   cap2 = [for (p=cap) up(90, p=p)];
 110//   vnf1 = vnf_vertex_array(points=wall_points, col_wrap=true);
 111//   vnf2 = vnf_vertex_array(points=cap1, col_wrap=true);
 112//   vnf3 = vnf_vertex_array(points=cap2, col_wrap=true, reverse=true);
 113//   vnf_polyhedron([vnf1, vnf2, vnf3]);
 114// Example(3D): Building a Multi-Stage Cylindrical Ramp
 115//   include <BOSL2/rounding.scad>
 116//   major_r = 50;
 117//   groove_profile = [
 118//       [-10,0], each arc(points=[[-7,0],[0,-3],[7,0]]), [10,0]
 119//   ];
 120//   ramp_profile = [ [-10,25], [90,25], [180,5], [190,5] ];
 121//   rgroove = apply(right(major_r) * xrot(90), path3d(groove_profile));
 122//   rprofile = round_corners(ramp_profile, radius=20, closed=false, $fn=72);
 123//   vnf = vnf_vertex_array([
 124//       for (a = [ramp_profile[0].x : 1 : last(ramp_profile).x]) let(
 125//           z = lookup(a,rprofile),
 126//           m = zrot(a) * up(z)
 127//       )
 128//       apply(m, [ [rgroove[0].x,0,-z], each rgroove, [last(rgroove).x,0,-z] ])
 129//   ], caps=true, col_wrap=true, reverse=true);
 130//   vnf_polyhedron(vnf, convexity=8);
 131function vnf_vertex_array(
 132    points,
 133    caps, cap1, cap2,
 134    col_wrap=false,
 135    row_wrap=false,
 136    reverse=false,
 137    style="default",
 138    triangulate = false
 139) =
 140    assert(!(any([caps,cap1,cap2]) && !col_wrap), "col_wrap must be true if caps are requested")
 141    assert(!(any([caps,cap1,cap2]) && row_wrap), "Cannot combine caps with row_wrap")
 142    assert(in_list(style,["default","alt","quincunx", "convex","concave", "min_edge","min_area","flip1","flip2"]))
 143    assert(is_matrix(points[0], n=3),"Point array has the wrong shape or points are not 3d")
 144    assert(is_consistent(points), "Non-rectangular or invalid point array")
 145    assert(is_bool(triangulate))
 146    let(
 147        pts = flatten(points),
 148        pcnt = len(pts),
 149        rows = len(points),
 150        cols = len(points[0])
 151    )
 152    rows<=1 || cols<=1 ? EMPTY_VNF :
 153    let(
 154        cap1 = first_defined([cap1,caps,false]),
 155        cap2 = first_defined([cap2,caps,false]),
 156        colcnt = cols - (col_wrap?0:1),
 157        rowcnt = rows - (row_wrap?0:1),
 158        verts = [
 159            each pts,
 160            if (style=="quincunx")
 161                for (r = [0:1:rowcnt-1], c = [0:1:colcnt-1])
 162                   let(
 163                       i1 = ((r+0)%rows)*cols + ((c+0)%cols),
 164                       i2 = ((r+1)%rows)*cols + ((c+0)%cols),
 165                       i3 = ((r+1)%rows)*cols + ((c+1)%cols),
 166                       i4 = ((r+0)%rows)*cols + ((c+1)%cols)
 167                   )
 168                   mean([pts[i1], pts[i2], pts[i3], pts[i4]])
 169        ],
 170        allfaces = [
 171            if (cap1) count(cols,reverse=!reverse),
 172            if (cap2) count(cols,(rows-1)*cols, reverse=reverse),
 173            for (r = [0:1:rowcnt-1], c=[0:1:colcnt-1])
 174               each
 175               let(
 176                   i1 = ((r+0)%rows)*cols + ((c+0)%cols),
 177                   i2 = ((r+1)%rows)*cols + ((c+0)%cols),
 178                   i3 = ((r+1)%rows)*cols + ((c+1)%cols),
 179                   i4 = ((r+0)%rows)*cols + ((c+1)%cols),
 180                   faces =
 181                        style=="quincunx"?
 182                          let(i5 = pcnt + r*colcnt + c)
 183                          [[i1,i5,i2],[i2,i5,i3],[i3,i5,i4],[i4,i5,i1]]
 184                      : style=="alt" || (style=="flip1" && ((r+c)%2==0)) || (style=="flip2" && ((r+c)%2==1)) || (style=="random" && rands(0,1,1)[0]<.5)?
 185                          [[i1,i4,i2],[i2,i4,i3]]
 186                      : style=="min_area"?
 187                          let(
 188                               area42 = norm(cross(pts[i2]-pts[i1], pts[i4]-pts[i1]))+norm(cross(pts[i4]-pts[i3], pts[i2]-pts[i3])),
 189                               area13 = norm(cross(pts[i1]-pts[i4], pts[i3]-pts[i4]))+norm(cross(pts[i3]-pts[i2], pts[i1]-pts[i2])),
 190                               minarea_edge = area42 < area13 + EPSILON
 191                                 ? [[i1,i4,i2],[i2,i4,i3]]
 192                                 : [[i1,i3,i2],[i1,i4,i3]]
 193                          )
 194                          minarea_edge
 195                      : style=="min_edge"?
 196                          let(
 197                               d42=norm(pts[i4]-pts[i2]),
 198                               d13=norm(pts[i1]-pts[i3]),
 199                               shortedge = d42<d13+EPSILON
 200                                 ? [[i1,i4,i2],[i2,i4,i3]]
 201                                 : [[i1,i3,i2],[i1,i4,i3]]
 202                          )
 203                          shortedge
 204                      : style=="convex"?
 205                          let(   // Find normal for 3 of the points.  Is the other point above or below?
 206                              n = (reverse?-1:1)*cross(pts[i2]-pts[i1],pts[i3]-pts[i1]),
 207                              convexfaces = n==0
 208                                ? [[i1,i4,i3]]
 209                                : n*pts[i4] > n*pts[i1]
 210                                    ? [[i1,i4,i2],[i2,i4,i3]]
 211                                    : [[i1,i3,i2],[i1,i4,i3]]
 212                          )
 213                          convexfaces
 214                      : style=="concave"?
 215                          let(   // Find normal for 3 of the points.  Is the other point above or below?
 216                              n = (reverse?-1:1)*cross(pts[i2]-pts[i1],pts[i3]-pts[i1]),
 217                              concavefaces = n==0
 218                                ? [[i1,i4,i3]]
 219                                : n*pts[i4] <= n*pts[i1]
 220                                    ? [[i1,i4,i2],[i2,i4,i3]]
 221                                    : [[i1,i3,i2],[i1,i4,i3]]
 222                          )
 223                          concavefaces
 224                      : [[i1,i3,i2],[i1,i4,i3]],
 225                   // remove degenerate faces
 226                   culled_faces= [for(face=faces)
 227                       if (norm(cross(verts[face[1]]-verts[face[0]],
 228                                      verts[face[2]]-verts[face[0]]))>EPSILON)
 229                           face
 230                   ],
 231                   rfaces = reverse? [for (face=culled_faces) reverse(face)] : culled_faces
 232               )
 233               rfaces,
 234        ],
 235        vnf = [verts, allfaces]
 236    ) triangulate? vnf_triangulate(vnf) : vnf;
 237
 238
 239// Function: vnf_tri_array()
 240// Synopsis: Returns a VNF from an array of points.
 241// SynTags: VNF
 242// Topics: VNF Generators, Lists
 243// See Also: vnf_vertex_array(), vnf_join(), vnf_from_polygons(), vnf_from_region()
 244// Usage:
 245//   vnf = vnf_tri_array(points, [row_wrap], [reverse])
 246// Description:
 247//   Produces a VNF from an array of points where each row length can differ from the adjacent rows by up to 2 in length.  This enables
 248//   the construction of triangular VNF patches.  The resulting VNF can be wrapped along the rows by setting `row_wrap` to true.
 249//   You cannot wrap columns: if you need to do that you'll need to merge two VNF arrays that share edges.  Degenerate faces
 250//   are not included in the output, but if this results in unused vertices they will still appear in the output.
 251// Arguments:
 252//   points = List of point lists for each row
 253//   row_wrap = If true then add faces connecting the first row and last row.  These rows must differ by at most 2 in length.
 254//   reverse = Set this to reverse the direction of the faces
 255// Example(3D,NoAxes): Each row has one more point than the preceeding one.
 256//   pts = [for(y=[1:1:10]) [for(x=[0:y-1]) [x,y,y]]];
 257//   vnf = vnf_tri_array(pts);
 258//   vnf_wireframe(vnf,width=0.1);
 259//   color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9);
 260// Example(3D,NoAxes): Each row has two more points than the preceeding one.
 261//   pts = [for(y=[0:2:10]) [for(x=[-y/2:y/2]) [x,y,y]]];
 262//   vnf = vnf_tri_array(pts);
 263//   vnf_wireframe(vnf,width=0.1);
 264//   color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9);
 265// Example(3D): Merging two VNFs to construct a cone with one point length change between rows.
 266//   pts1 = [for(z=[0:10]) path3d(arc(3+z,r=z/2+1, angle=[0,180]),10-z)];
 267//   pts2 = [for(z=[0:10]) path3d(arc(3+z,r=z/2+1, angle=[180,360]),10-z)];
 268//   vnf = vnf_join([vnf_tri_array(pts1),
 269//                     vnf_tri_array(pts2)]);
 270//   color("green")vnf_wireframe(vnf,width=0.1);
 271//   vnf_polyhedron(vnf);
 272// Example(3D): Cone with length change two between rows
 273//   pts1 = [for(z=[0:1:10]) path3d(arc(3+2*z,r=z/2+1, angle=[0,180]),10-z)];
 274//   pts2 = [for(z=[0:1:10]) path3d(arc(3+2*z,r=z/2+1, angle=[180,360]),10-z)];
 275//   vnf = vnf_join([vnf_tri_array(pts1),
 276//                    vnf_tri_array(pts2)]);
 277//   color("green")vnf_wireframe(vnf,width=0.1);
 278//   vnf_polyhedron(vnf);
 279// Example(3D,NoAxes): Point count can change irregularly
 280//   lens = [10,9,7,5,6,8,8,10];
 281//   pts = [for(y=idx(lens)) lerpn([-lens[y],y,y],[lens[y],y,y],lens[y])];
 282//   vnf = vnf_tri_array(pts);
 283//   vnf_wireframe(vnf,width=0.1);
 284//   color("red")move_copies(flatten(pts)) sphere(r=.15,$fn=9);
 285function vnf_tri_array(points, row_wrap=false, reverse=false) =
 286    let(
 287       lens = [for(row=points) len(row)],
 288       rowstarts = [0,each cumsum(lens)],
 289       faces =
 290          [for(i=[0:1:len(points) - 1 - (row_wrap ? 0 : 1)]) each
 291            let(
 292                rowstart = rowstarts[i],
 293                nextrow = select(rowstarts,i+1),
 294                delta = select(lens,i+1)-lens[i]
 295            )
 296            delta == 0 ?
 297              [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow],
 298               for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1] : [j+rowstart+1, j+nextrow+1, j+nextrow]] :
 299            delta == 1 ?
 300              [for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1],
 301               for(j=[0:1:lens[i]-1]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow]] :
 302            delta == -1 ?
 303              [for(j=[0:1:lens[i]-3]) reverse ? [j+rowstart+1, j+nextrow, j+nextrow+1]: [j+rowstart+1, j+nextrow+1, j+nextrow],
 304               for(j=[0:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow] : [j+rowstart, j+rowstart+1, j+nextrow]] :
 305            let(count = floor((lens[i]-1)/2))
 306            delta == 2 ?
 307              [
 308               for(j=[0:1:count-1]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+1] : [j+rowstart, j+rowstart+1, j+nextrow+1],       // top triangles left
 309               for(j=[count:1:lens[i]-2]) reverse ? [j+rowstart+1, j+rowstart, j+nextrow+2] : [j+rowstart, j+rowstart+1, j+nextrow+2], // top triangles right
 310               for(j=[0:1:count]) reverse ? [j+rowstart, j+nextrow, j+nextrow+1] : [j+rowstart, j+nextrow+1, j+nextrow],                        // bot triangles left
 311               for(j=[count+1:1:select(lens,i+1)-2]) reverse ? [j+rowstart-1, j+nextrow, j+nextrow+1] : [j+rowstart-1, j+nextrow+1, j+nextrow], // bot triangles right
 312              ] :
 313            delta == -2 ?
 314              [
 315               for(j=[0:1:count-2]) reverse ? [j+nextrow, j+nextrow+1, j+rowstart+1] : [j+nextrow, j+rowstart+1, j+nextrow+1],
 316               for(j=[count-1:1:lens[i]-4]) reverse ? [j+nextrow,j+nextrow+1,j+rowstart+2] : [j+nextrow,j+rowstart+2, j+nextrow+1],
 317               for(j=[0:1:count-1]) reverse ? [j+nextrow, j+rowstart+1, j+rowstart] : [j+nextrow, j+rowstart, j+rowstart+1],
 318               for(j=[count:1:select(lens,i+1)]) reverse ? [ j+nextrow-1, j+rowstart+1, j+rowstart]: [ j+nextrow-1, j+rowstart, j+rowstart+1],
 319              ] :
 320            assert(false,str("Unsupported row length difference of ",delta, " between row ",i," and ",(i+1)%len(points)))
 321          ],
 322       verts = flatten(points),
 323       culled_faces=
 324           [for(face=faces)
 325               if (norm(verts[face[0]]-verts[face[1]])>EPSILON &&
 326                   norm(verts[face[1]]-verts[face[2]])>EPSILON &&
 327                   norm(verts[face[2]]-verts[face[0]])>EPSILON)
 328                   face
 329           ]
 330    )
 331    [flatten(points), culled_faces];
 332
 333
 334
 335// Function: vnf_join()
 336// Synopsis: Returns a single VNF structure from a list of VNF structures.
 337// SynTags: VNF
 338// Topics: VNF Generators, Lists
 339// See Also: vnf_tri_array(), vnf_vertex_array(), vnf_from_polygons(), vnf_from_region()
 340// Usage:
 341//   vnf = vnf_join([VNF, VNF, VNF, ...]);
 342// Description:
 343//   Given a list of VNF structures, merges them all into a single VNF structure.
 344//   Combines all the points of the input VNFs and labels the faces appropriately.
 345//   All the points in the input VNFs will appear in the output, even if they are
 346//   duplicates of each other.  It is valid to repeat points in a VNF, but if you
 347//   with to remove the duplicates that will occur along joined edges, use {{vnf_merge_points()}}.
 348//   .
 349//   Note that this is a tool for manipulating polyhedron data.  It is for
 350//   building up a full polyhedron from partial polyhedra.
 351//   It is *not* a union operator for VNFs.  The VNFs to be joined must not intersect each other,
 352//   except at edges, or the result will be an invalid polyhedron.  Similarly the
 353//   result must not have any other illegal polyhedron characteristics, such as creating
 354//   more than two faces sharing the same edge.
 355//   If you want a valid result it is your responsibility to ensure that the polyhedron
 356//   has no holes, no intersecting faces or edges, and obeys all the requirements
 357//   that CGAL expects.
 358//   .
 359//   For example, if you combine two pyramids to try to make an octahedron, the result will
 360//   be invalid because of the two internal faces created by the pyramid bases.  A valid
 361//   use would be to build a cube missing one face and a pyramid missing its base and
 362//   then join them into a cube with a point.
 363// Arguments:
 364//   vnfs = a list of the VNFs to joint into one VNF.
 365// Example(3D,VPR=[60,0,26],VPD=55,VPT=[5.6,-5.3,9.8]): Here is a VNF where the top face is missing.  It is not a valid polyhedron like this, but we can use it as a building block to make a polyhedron.
 366//   bottom = vnf_vertex_array([path3d(rect(8)), path3d(rect(5),4)],col_wrap=true,cap1=true);
 367//   vnf_polyhedron(bottom);
 368// Example(3D,VPR=[60,0,26],VPD=55,VPT=[5.6,-5.3,9.8]): Here is a VNF that also has a missing face.
 369//   triangle = yrot(-90,path3d(regular_ngon(n=3,side=5,anchor=LEFT)));
 370//   top = up(4,vnf_vertex_array([list_set(right(2.5,triangle),0,[0,0,7]),
 371//                               right(6,triangle)
 372//                               ], col_wrap=true, cap2=true));
 373//   vnf_polyhedron(zrot(90,top));
 374// Example(3D,VPR=[60,0,26],VPD=55,VPT=[5.6,-5.3,9.8]): Using vnf_join combines the two VNFs into a single VNF.  Note that they share an edge.  But the result still isn't closed, so it is not yet a valid polyhedron.
 375//   bottom = vnf_vertex_array([path3d(rect(8)), path3d(rect(5),4)],col_wrap=true,cap1=true);
 376//   triangle = yrot(-90,path3d(regular_ngon(n=3,side=5,anchor=LEFT)));
 377//   top = up(4,vnf_vertex_array([list_set(right(2.5,triangle),0,[0,0,7]),
 378//                                right(6,triangle)
 379//                               ], col_wrap=true, cap2=true));
 380//   full = vnf_join([bottom,zrot(90,top)]);
 381//   vnf_polyhedron(full);
 382// Example(3D,VPR=[60,0,26],VPD=55,VPT=[5.6,-5.3,9.8]): If we add enough pieces, and the pieces are all consistent with each other, then we can arrive at a valid polyhedron like this one.  To be valid you need to meet all the CGAL requirements: every edge has exactly two faces, all faces are in clockwise order, no intersections of edges.
 383//   bottom = vnf_vertex_array([path3d(rect(8)), path3d(rect(5),4)],col_wrap=true,cap1=true);
 384//   triangle = yrot(-90,path3d(regular_ngon(n=3,side=5,anchor=LEFT)));
 385//   top = up(4,vnf_vertex_array([list_set(right(2.5,triangle),0,[0,0,7]),
 386//                                right(6,triangle)
 387//                               ], col_wrap=true, cap2=true));
 388//   full = vnf_join([bottom,
 389//                     for(theta=[0:90:359]) zrot(theta,top)
 390//                    ]);
 391//   vnf_polyhedron(full);
 392// Example(3D): The vnf_join function is not a union operator for polyhedra.  If any faces intersect, like they do in this example where we combine the faces of two cubes, the result is invalid and will give rise to CGAL errors when you add more objects into the model.
 393//   cube1 = cube(5);
 394//   cube2 = move([2,2,2],cube1);
 395//   badvnf = vnf_join([cube1,cube2]);
 396//   vnf_polyhedron(badvnf);
 397//   right(2.5)up(3)color("red")
 398//         text3d("Invalid",size=1,anchor=CENTER,
 399//         orient=FRONT,h=.1);
 400function vnf_join(vnfs) =
 401    assert(is_vnf_list(vnfs) , "Input must be a list of VNFs")
 402    len(vnfs)==1 ? vnfs[0]
 403    :
 404    let (
 405        offs  = cumsum([ 0, for (vnf = vnfs) len(vnf[0]) ]),
 406        verts = [for (vnf=vnfs) each vnf[0]],
 407        faces =
 408            [ for (i = idx(vnfs))
 409                let( faces = vnfs[i][1] )
 410                for (face = faces)
 411                    if ( len(face) >= 3 )
 412                        [ for (j = face)
 413                            assert( j>=0 && j<len(vnfs[i][0]),
 414                                    str("VNF number ", i, " has a face indexing an nonexistent vertex") )
 415                            offs[i] + j ]
 416            ]
 417    )
 418    [verts,faces];
 419
 420
 421
 422// Function: vnf_from_polygons()
 423// Synopsis: Returns a VNF from a list of 3D polygons.
 424// SynTags: VNF
 425// Topics: VNF Generators, Lists
 426// See Also: vnf_tri_array(), vnf_join(), vnf_vertex_array(), vnf_from_region()
 427// Usage:
 428//   vnf = vnf_from_polygons(polygons, [eps]);
 429// Description:
 430//   Given a list of 3D polygons, produces a VNF containing those polygons.
 431//   It is up to the caller to make sure that the points are in the correct order to make the face
 432//   normals point outwards.  No checking for duplicate vertices is done.  If you want to
 433//   remove duplicate vertices use {{vnf_merge_points()}}.  Polygons with zero area are discarded from the face list by default.
 434//   If you give non-coplanar faces an error is displayed.  These checks increase run time by about 2x for triangular polygons, but
 435//   about 10x for pentagons; the checks can be disabled by setting fast=true.  
 436// Arguments:
 437//   polygons = The list of 3D polygons to turn into a VNF
 438//   fast = Set to true to skip area and coplanarity checks for increased speed.  Default: false
 439//   eps = Polygons with area small than this are discarded.  Default: EPSILON
 440function vnf_from_polygons(polygons,fast=false,eps=EPSILON) =
 441   assert(is_list(polygons) && is_path(polygons[0]),"Input should be a list of polygons")
 442   let(
 443       offs = cumsum([0, for(p=polygons) len(p)]),
 444       faces = [for(i=idx(polygons))
 445                  let(
 446                      area=fast ? 1 : polygon_area(polygons[i]),
 447                      dummy=assert(is_def(area) || is_collinear(polygons[i],eps=eps),str("Polygon ", i, " is not coplanar"))
 448                  )
 449                  if (is_def(area) && area > eps)
 450                    [for (j=idx(polygons[i])) offs[i]+j]
 451               ]
 452   )
 453   [flatten(polygons), faces];
 454
 455
 456
 457
 458function _path_path_closest_vertices(path1,path2) =
 459    let(
 460        dists = [for (i=idx(path1)) let(j=closest_point(path1[i],path2)) [j,norm(path2[j]-path1[i])]],
 461        i1 = min_index(column(dists,1)),
 462        i2 = dists[i1][0]
 463    ) [dists[i1][1], i1, i2];
 464
 465
 466function _join_paths_at_vertices(path1,path2,v1,v2) =
 467    let(
 468        repeat_start = !approx(path1[v1],path2[v2]),
 469        path1 = clockwise_polygon(list_rotate(path1,v1)),
 470        path2 = ccw_polygon(list_rotate(path2,v2))
 471    )
 472    [
 473        each path1,
 474        if (repeat_start) path1[0],
 475        each path2,
 476        if (repeat_start) path2[0],
 477    ];
 478
 479
 480/// Internal Function: _cleave_connected_region(region, eps)
 481/// Description:
 482///   Given a region that is connected and has its outer border in region[0],
 483///   produces a overlapping connected path to join internal holes to
 484///   the outer border without adding points. Output is a single non-simple polygon.
 485/// Requirements:
 486///   It expects that all region paths be simple closed paths, with region[0] CW and
 487///   the other paths CCW and encircled by region[0]. The input region paths are also
 488///   supposed to be disjoint except for common vertices and common edges but with
 489///   no crossings. It may return `undef` if these conditions are not met.
 490///   This function implements an extension of the algorithm discussed in:
 491///   https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
 492function _cleave_connected_region(region, eps=EPSILON) =
 493    len(region)==1 ? region[0] :
 494    let(
 495        outer   = deduplicate(region[0]),             //
 496        holes   = [for(i=[1:1:len(region)-1])         // deduplication possibly unneeded
 497                      deduplicate( region[i] ) ],     //
 498        extridx = [for(li=holes) max_index(column(li,0)) ],
 499        // the right extreme vertex for each hole sorted by decreasing x values
 500        extremes = sort( [for(i=idx(holes)) [ i, extridx[i], -holes[i][extridx[i]].x] ], idx=2 )
 501    )
 502    _polyHoles(outer, holes, extremes, eps, 0);
 503
 504
 505// connect the hole paths one at a time to the outer path.
 506// 'extremes' is the list of the right extreme vertex of each hole sorted by decreasing abscissas
 507// see: _cleave_connected_region(region, eps)
 508function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) =
 509    let(
 510        extr = extremes[n],    //
 511        hole = holes[extr[0]], // hole path to bridge to the outer path
 512        ipt  = extr[1],        // index of the hole point with maximum abscissa
 513        brdg = _bridge(hole[ipt], outer, eps)  // the index of a point in outer to bridge hole[ipt] to
 514    )
 515    brdg == undef ? undef :
 516    let(
 517        l  = len(outer),
 518        lh = len(hole),
 519        // the new outer polygon bridging the hole to the old outer
 520        npoly =
 521            approx(outer[brdg], hole[ipt], eps)
 522            ?   [ for(i=[brdg:  1: brdg+l])   outer[i%l] ,
 523                  for(i=[ipt+1: 1: ipt+lh-1]) hole[i%lh] ]
 524            :   [ for(i=[brdg:  1: brdg+l])   outer[i%l] ,
 525                  for(i=[ipt:   1: ipt+lh])   hole[i%lh] ]
 526    )
 527    n==len(holes)-1 ?  npoly :
 528    _polyHoles(npoly, holes, extremes, eps, n+1);
 529
 530// find a point in outer to be connected to pt in the interior of outer
 531// by a segment that not cross or touch any non adjacente edge of outer.
 532// return the index of a vertex in the outer path where the bridge should end
 533// see _polyHoles(outer, holes, extremes, eps)
 534function _bridge(pt, outer,eps) =
 535    // find the intersection of a ray from pt to the right
 536    // with the boundary of the outer cycle
 537    let(
 538        l    = len(outer),
 539        crxs =
 540            let( edges = pair(outer,wrap=true) )
 541            [for( i = idx(edges) )
 542                let( edge = edges[i] )
 543                // consider just descending outer edges at right of pt crossing ordinate pt.y
 544                if(    (edge[0].y >  pt.y) //+eps)
 545                    && (edge[1].y <= pt.y)
 546                    && _is_at_left(pt, [edge[1], edge[0]], eps) )
 547                    [ i,
 548                      // the point of edge with ordinate pt.y
 549                      abs(pt.y-edge[1].y)<eps ? edge[1] :
 550                      let( u = (pt-edge[1]).y / (edge[0]-edge[1]).y )
 551                      (1-u)*edge[1] + u*edge[0]
 552                    ]
 553             ]
 554    )
 555    crxs == [] ? undef :
 556    let(
 557        // the intersection point of the nearest edge to pt with minimum slope
 558        minX    = min([for(p=crxs) p[1].x]),
 559        crxcand = [for(crx=crxs) if(crx[1].x < minX+eps) crx ], // nearest edges
 560        nearest = min_index([for(crx=crxcand)
 561                                (outer[crx[0]].x - pt.x) / (outer[crx[0]].y - pt.y) ]), // minimum slope
 562        proj    = crxcand[nearest],
 563        vert0   = outer[proj[0]],    // the two vertices of the nearest crossing edge
 564        vert1   = outer[(proj[0]+1)%l],
 565        isect   = proj[1]            // the intersection point
 566    )
 567    norm(pt-vert1) < eps ? (proj[0]+1)%l : // if pt touches an outer vertex, return its index
 568    // as vert0.y > pt.y then pt!=vert0
 569    norm(pt-isect) < eps ? undef :         // if pt touches the middle of an outer edge -> error
 570    let(
 571        // the edge [vert0, vert1] necessarily satisfies vert0.y > vert1.y
 572        // indices of candidates to an outer bridge point
 573        cand  =
 574            (vert0.x > pt.x)
 575            ?   [ proj[0],
 576                  // select reflex vertices inside of the triangle [pt, vert0, isect]
 577                  for(i=idx(outer))
 578                      if( _tri_class(select(outer,i-1,i+1),eps) <= 0
 579                          && _pt_in_tri(outer[i], [pt, vert0, isect], eps)>=0 )
 580                        i
 581                ]
 582            :   [ (proj[0]+1)%l,
 583                  // select reflex vertices inside of the triangle [pt, isect, vert1]
 584                  for(i=idx(outer))
 585                      if( _tri_class(select(outer,i-1,i+1),eps) <= 0
 586                          &&  _pt_in_tri(outer[i], [pt, isect, vert1], eps)>=0 )
 587                        i
 588                ],
 589        // choose the candidate outer[i] such that the line [pt, outer[i]] has minimum slope
 590        // among those with minimum slope choose the nearest to pt
 591        slopes  = [for(i=cand) 1-abs(outer[i].x-pt.x)/norm(outer[i]-pt) ],
 592        min_slp = min(slopes),
 593        cand2   = [for(i=idx(cand)) if(slopes[i]<=min_slp+eps) cand[i] ],
 594        nearest = min_index([for(i=cand2) norm(pt-outer[i]) ])
 595    )
 596    cand2[nearest];
 597
 598
 599// Function: vnf_from_region()
 600// Synopsis: Returns a 3D VNF given a 2D region.
 601// SynTags: VNF
 602// Topics: VNF Generators, Lists
 603// See Also: vnf_vertex_array(), vnf_tri_array(), vnf_join(), vnf_from_polygons()
 604// Usage:
 605//   vnf = vnf_from_region(region, [transform], [reverse]);
 606// Description:
 607//   Given a (two-dimensional) region, applies the given transformation matrix to it and makes a (three-dimensional) triangulated VNF of
 608//   faces for that region, reversed if desired.
 609// Arguments:
 610//   region = The region to convert to a vnf.
 611//   transform = If given, a transformation matrix to apply to the faces generated from the region.  Default: No transformation applied.
 612//   reverse = If true, reverse the normals of the faces generated from the region.  An untransformed region will have face normals pointing `UP`.  Default: false
 613// Example(3D):
 614//   region = [square([20,10],center=true),
 615//             right(5,square(4,center=true)),
 616//             left(5,square(6,center=true))];
 617//   vnf = vnf_from_region(region);
 618//   color("gray")down(.125)
 619//        linear_extrude(height=.125)region(region);
 620//   vnf_wireframe(vnf,width=.25);
 621function vnf_from_region(region, transform, reverse=false) =
 622    let (
 623        region = [for (path = region) deduplicate(path, closed=true)],
 624        regions = region_parts(force_region(region)),
 625        vnfs =
 626            [
 627                for (rgn = regions)
 628                let(
 629                    cleaved = path3d(_cleave_connected_region(rgn))
 630                )
 631                assert( cleaved, "The region is invalid")
 632                let(
 633                    face = is_undef(transform)? cleaved : apply(transform,cleaved),
 634                    faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i]
 635                ) [face, [faceidxs]]
 636            ],
 637        outvnf = vnf_join(vnfs)
 638    )
 639    vnf_triangulate(outvnf);
 640
 641
 642
 643// Section: VNF Testing and Access
 644
 645
 646// Function: is_vnf()
 647// Synopsis: Returns true given a VNF-like structure.
 648// Topics: VNF Manipulation
 649// See Also: is_vnf_list(), vnf_vertices(), vnf_faces()
 650// Usage:
 651//   bool = is_vnf(x);
 652// Description:
 653//   Returns true if the given value looks like a VNF structure.
 654function is_vnf(x) =
 655    is_list(x) &&
 656    len(x)==2 &&
 657    is_list(x[0]) &&
 658    is_list(x[1]) &&
 659    (x[0]==[] || (len(x[0])>=3 && is_vector(x[0][0],3))) &&
 660    (x[1]==[] || is_vector(x[1][0]));
 661
 662
 663// Function: is_vnf_list()
 664// Synopsis: Returns true given a list of VNF-like structures.
 665// Topics: VNF Manipulation
 666// See Also: is_vnf(), vnf_vertices(), vnf_faces()
 667//
 668// Description: Returns true if the given value looks passingly like a list of VNF structures.
 669function is_vnf_list(x) = is_list(x) && all([for (v=x) is_vnf(v)]);
 670
 671
 672// Function: vnf_vertices()
 673// Synopsis: Returns the list of vertex points from a VNF.
 674// Topics: VNF Manipulation
 675// See Also: is_vnf(), is_vnf_list(), vnf_faces()
 676// Description: Given a VNF structure, returns the list of vertex points.
 677function vnf_vertices(vnf) = vnf[0];
 678
 679
 680// Function: vnf_faces()
 681// Synopsis: Returns the list of faces from a VNF.
 682// Topics: VNF Manipulation
 683// See Also: is_vnf(), is_vnf_list(), vnf_vertices()
 684// Description: Given a VNF structure, returns the list of faces, where each face is a list of indices into the VNF vertex list.
 685function vnf_faces(vnf) = vnf[1];
 686
 687
 688
 689// Section: Altering the VNF Internals
 690
 691
 692// Function: vnf_reverse_faces()
 693// Synopsis: Reverses the faces of a VNF.
 694// SynTags: VNF
 695// Topics: VNF Manipulation
 696// See Also: vnf_quantize(), vnf_merge_points(), vnf_drop_unused_points(), vnf_triangulate(), vnf_slice(), vnf_unify_faces() 
 697// Usage:
 698//   rvnf = vnf_reverse_faces(vnf);
 699// Description:
 700//   Reverses the orientation of all the faces in the given VNF.
 701function vnf_reverse_faces(vnf) =
 702    [vnf[0], [for (face=vnf[1]) reverse(face)]];
 703
 704
 705// Function: vnf_quantize()
 706// Synopsis: Quantizes the vertex coordinates of a VNF.
 707// SynTags: VNF
 708// Topics: VNF Manipulation
 709// See Also: vnf_reverse_faces(), vnf_merge_points(), vnf_drop_unused_points(), vnf_triangulate(), vnf_slice() 
 710// Usage:
 711//   vnf2 = vnf_quantize(vnf,[q]);
 712// Description:
 713//   Quantizes the vertex coordinates of the VNF to the given quanta `q`.
 714// Arguments:
 715//   vnf = The VNF to quantize.
 716//   q = The quanta to quantize the VNF coordinates to.
 717function vnf_quantize(vnf,q=pow(2,-12)) =
 718    [[for (pt = vnf[0]) quant(pt,q)], vnf[1]];
 719
 720
 721
 722// Function: vnf_merge_points()
 723// Synopsis: Consolidates duplicate vertices of a VNF.
 724// SynTags: VNF
 725// Topics: VNF Manipulation
 726// See Also: vnf_reverse_faces(), vnf_quantize(), vnf_drop_unused_points(), vnf_triangulate(), vnf_slice(), vnf_unify_faces() 
 727// Usage:
 728//   new_vnf = vnf_merge_points(vnf, [eps]);
 729// Description:
 730//   Given a VNF, consolidates all duplicate vertices with a tolerance `eps`, relabeling the faces as necessary,
 731//   and eliminating any face with fewer than 3 vertices.  Unreferenced vertices of the input VNF are not dropped.
 732//   To remove such vertices uses {{vnf_drop_unused_points()}}.
 733// Arguments:
 734//   vnf = a VNF to consolidate
 735//   eps = the tolerance in finding duplicates. Default: EPSILON
 736function vnf_merge_points(vnf,eps=EPSILON) =
 737    let(
 738        verts = vnf[0],
 739        dedup  = vector_search(verts,eps,verts),                 // collect vertex duplicates
 740        map    = [for(i=idx(verts)) min(dedup[i]) ],             // remap duplic vertices
 741        offset = cumsum([for(i=idx(verts)) map[i]==i ? 0 : 1 ]), // remaping face vertex offsets
 742        map2   = list(idx(verts))-offset,                        // map old vertex indices to new indices
 743        nverts = [for(i=idx(verts)) if(map[i]==i) verts[i] ],    // this doesn't eliminate unreferenced vertices
 744        nfaces =
 745            [ for(face=vnf[1])
 746                let(
 747                    nface = [ for(vi=face) map2[map[vi]] ],
 748                    dface = [for (i=idx(nface))
 749                                if( nface[i]!=nface[(i+1)%len(nface)])
 750                                    nface[i] ]
 751                )
 752                if(len(dface) >= 3) dface
 753            ]
 754    )
 755    [nverts, nfaces];
 756
 757
 758// Function: vnf_drop_unused_points()
 759// Synopsis: Removes unreferenced vertices from a VNF.
 760// SynTags: VNF
 761// Topics: VNF Manipulation
 762// See Also: vnf_reverse_faces(), vnf_quantize(), vnf_merge_points(), vnf_triangulate(), vnf_slice(), vnf_unify_faces() 
 763// Usage:
 764//   clean_vnf = vnf_drop_unused_points(vnf);
 765// Description:
 766//   Remove all unreferenced vertices from a VNF.  Note that in most
 767//   cases unreferenced vertices cause no harm, and this function may
 768//   be slow on large VNFs.
 769function vnf_drop_unused_points(vnf) =
 770    let(
 771        flat = flatten(vnf[1]),
 772        ind  = _link_indicator(flat,0,len(vnf[0])-1),
 773        verts = [for(i=idx(vnf[0])) if(ind[i]==1) vnf[0][i] ],
 774        map   = cumsum(ind)
 775    )
 776    [ verts, [for(face=vnf[1]) [for(v=face) map[v]-1 ] ] ];
 777
 778function _link_indicator(l,imin,imax) =
 779    len(l) == 0  ? repeat(imax-imin+1,0) :
 780    imax-imin<100 || len(l)<400 ? [for(si=search(list([imin:1:imax]),l,1)) si!=[] ? 1: 0 ] :
 781    let(
 782        pivot   = floor((imax+imin)/2),
 783        lesser  = [ for(li=l) if( li< pivot) li ],
 784        greater = [ for(li=l) if( li> pivot) li ]
 785    )
 786    concat( _link_indicator(lesser ,imin,pivot-1),
 787            search(pivot,l,1) ? 1 : 0 ,
 788            _link_indicator(greater,pivot+1,imax) ) ;
 789
 790// Function: vnf_triangulate()
 791// Synopsis: Triangulates the faces of a VNF.
 792// SynTags: VNF
 793// Topics: VNF Manipulation
 794// See Also: vnf_reverse_faces(), vnf_quantize(), vnf_merge_points(), vnf_drop_unused_points(), vnf_slice(), vnf_unify_faces() 
 795// Usage:
 796//   vnf2 = vnf_triangulate(vnf);
 797// Description:
 798//   Triangulates faces in the VNF that have more than 3 vertices.
 799// Arguments:
 800//   vnf = VNF to triangulate
 801// Example(3D):
 802//   include <BOSL2/polyhedra.scad>
 803//   vnf = zrot(33,regular_polyhedron_info("vnf", "dodecahedron", side=12));
 804//   vnf_polyhedron(vnf);
 805//   triangulated = vnf_triangulate(vnf);
 806//   color("red")vnf_wireframe(triangulated,width=.3);
 807function vnf_triangulate(vnf) =
 808    let(
 809        verts = vnf[0],
 810        faces = [for (face=vnf[1])
 811                    each (len(face)==3 ? [face] :
 812                    let( tris = polygon_triangulate(verts, face) )
 813                    assert( tris!=undef, "Some `vnf` face cannot be triangulated.")
 814                    tris ) ]
 815    )
 816    [verts, faces];
 817
 818
 819
 820// Function: vnf_unify_faces()
 821// Synopsis: Remove triangulation from VNF, returning a copy with full faces
 822// SynTags: VNF
 823// Topics: VNF Manipulation
 824// See Also: vnf_reverse_faces(), vnf_quantize(), vnf_merge_points(), vnf_triangulate(), vnf_slice() 
 825// Usage:
 826//   newvnf = vnf_unify_faces(vnf);
 827// Description:
 828//   When a VNF has been triangulated, the polygons that form the true faces have been chopped up into
 829//   triangles.  This can create problems for algorithms that operate on the VNF itself, where you might
 830//   want to be able to identify the true faces.  This function merges together the triangles that
 831//   form those true faces, turning a VNF where each true face is represented by a single entry
 832//   in the faces list of the VNF.  This function requires that the true faces have no internal vertices.
 833//   This will always be true for a triangulated VNF, but might fail for a VNF with some other
 834//   face partition.  If internal vertices are present, the output will include backtracking paths from
 835//   the boundary to all of those vertices.
 836// Arguments:
 837//   vnf = vnf whose faces you want to unify
 838// Example(3D,Med,NoAxes): Original prism on the left is triangulated.  On the right, the result of unifying the faces.
 839//   $fn=16;
 840//   poly = linear_sweep(hexagon(side=10),h=35);
 841//   vnf = vnf_unify_faces(poly);
 842//   vnf_wireframe(poly);
 843//   color([0,1,1,.70])vnf_polyhedron(poly);
 844//   right(25){
 845//     vnf_wireframe(vnf);
 846//     color([0,1,1,.70])vnf_polyhedron(vnf);
 847//   }
 848
 849function vnf_unify_faces(vnf) =
 850   let(
 851       faces = vnf[1],
 852       edges =  [for(i=idx(faces), edge=pair(faces[i],wrap=true))
 853                   [[min(edge),max(edge)],i]],
 854       normals = [for(face=faces) polygon_normal(select(vnf[0],face))],
 855       facelist = count(faces), //[for(i=[1:1:len(faces)-1]) i],
 856       newfaces = _detri_combine_faces(edges,faces,normals,facelist,0)
 857   )
 858   [vnf[0],newfaces];
 859
 860
 861function _detri_combine_faces(edgelist,faces,normals,facelist,curface) =
 862    curface==len(faces)? select(faces,facelist)
 863  : !in_list(curface,facelist) ? _detri_combine_faces(edgelist,faces,normals,facelist,curface+1)
 864  :
 865    let(
 866        thisface=faces[curface],
 867        neighbors = [for(i=idx(thisface))
 868                       let(
 869                           edgepair = search([sort(select(thisface,i,i+1))],edgelist,0)[0],
 870                           choices = select(edgelist,edgepair),
 871                           good_choice=[for(choice=choices)
 872                              if (choice[1]!=curface && in_list(choice[1],facelist) && normals[choice[1]]*normals[curface]>1-EPSILON)
 873                                choice],
 874                           d=assert(len(good_choice)<=1)
 875                       )
 876                       len(good_choice)==1 ? good_choice[0][1] : -1
 877                    ],
 878        // Check for duplicates in the neighbor list so we don't add them twice
 879        dups = search([for(n=neighbors) if (n>=0) n], neighbors,0),
 880        goodind = column(dups,0),
 881        newface = [for(i=idx(thisface))
 882                    each
 883                     !in_list(i,goodind) ? [thisface[i]]
 884                    :
 885                     let(
 886                         ind = search(select(thisface,i,i+1), faces[neighbors[i]])
 887                     )
 888                     select(faces[neighbors[i]],ind[0],ind[1]-1)
 889                   ],
 890        usedfaces = [for(n=neighbors) if (n>=0) n],
 891        faces = list_set(faces,curface,newface),
 892        facelist = list_remove_values(facelist,usedfaces)
 893     )
 894     _detri_combine_faces(edgelist,faces,normals,facelist,len(usedfaces)==0?curface+1:curface);
 895
 896
 897
 898function _vnf_sort_vertices(vnf, idx=[2,1,0]) =
 899    let(
 900        verts = vnf[0],
 901        faces = vnf[1],
 902        vidx = sortidx(verts, idx=idx),
 903        rvidx = sortidx(vidx),
 904        sorted_vnf = [
 905            [ for (i = vidx) verts[i] ],
 906            [ for (face = faces) [ for (i = face) rvidx[i] ] ],
 907        ]
 908    ) sorted_vnf;
 909
 910
 911
 912
 913// Function: vnf_slice()
 914// Synopsis: Slice the faces of a VNF along an axis.
 915// SynTags: VNF
 916// Topics: VNF Manipulation
 917// See Also: vnf_reverse_faces(), vnf_quantize(), vnf_merge_points(), vnf_drop_unused_points(), vnf_triangulate()
 918// Usage:
 919//   sliced = vnf_slice(vnf, dir, cuts);
 920// Description:
 921//   Slice the faces of a VNF along a specified axis direction at a given list of cut points.
 922//   The cut points can appear in any order.  You can use this to refine the faces of a VNF before
 923//   applying a nonlinear transformation to its vertex set.
 924// Arguments:
 925//   vnf = VNF to slice
 926//   dir = normal direction to the slices, either "X", "Y" or "Z"
 927//   cuts = X, Y or Z values where cuts occur
 928// Example(3D):
 929//   include <BOSL2/polyhedra.scad>
 930//   vnf = regular_polyhedron_info("vnf", "dodecahedron", side=12);
 931//   vnf_polyhedron(vnf);
 932//   sliced = vnf_slice(vnf, "X", [-6,-1,10]);
 933//   color("red")vnf_wireframe(sliced,width=.3);
 934function vnf_slice(vnf,dir,cuts) =
 935    let(
 936        //  Code below seems to be unnecessary
 937        //cuts = [for (cut=cuts) _shift_cut_plane(vnf,dir,cut)],
 938        vert = vnf[0],
 939        faces = [for(face=vnf[1]) select(vert,face)],
 940        poly_list = _slice_3dpolygons(faces, dir, cuts)
 941    )
 942    vnf_merge_points(vnf_from_polygons(poly_list));
 943
 944
 945function _shift_cut_plane(vnf,dir,cut,off=0.001) =
 946    let(
 947        I = ident(3),
 948        dir_ind = ord(dir)-ord("X"),
 949        verts = vnf[0],
 950        on_cut = [for (x = verts * I[dir_ind]) if(approx(x,cut,eps=1e-4)) 1] != []
 951    ) !on_cut? cut :
 952    _shift_cut_plane(vnf,dir,cut+off);
 953
 954
 955function _split_polygon_at_x(poly, x) =
 956    let(
 957        xs = column(poly,0)
 958    ) (min(xs) >= x || max(xs) <= x)? [poly] :
 959    let(
 960        poly2 = [
 961            for (p = pair(poly,true)) each [
 962                p[0],
 963                if(
 964                    (p[0].x < x && p[1].x > x) ||
 965                    (p[1].x < x && p[0].x > x)
 966                ) let(
 967                    u = (x - p[0].x) / (p[1].x - p[0].x)
 968                ) [
 969                    x,  // Important for later exact match tests
 970                    u*(p[1].y-p[0].y)+p[0].y
 971                ]
 972            ]
 973        ],
 974        out1 = [for (p = poly2) if(p.x <= x) p],
 975        out2 = [for (p = poly2) if(p.x >= x) p],
 976        out3 = [
 977            if (len(out1)>=3) each split_path_at_self_crossings(out1),
 978            if (len(out2)>=3) each split_path_at_self_crossings(out2),
 979        ],
 980        out = [for (p=out3) if (len(p) > 2) list_unwrap(p)]
 981    ) out;
 982
 983
 984function _split_2dpolygons_at_each_x(polys, xs, _i=0) =
 985    _i>=len(xs)? polys :
 986    _split_2dpolygons_at_each_x(
 987        [
 988            for (poly = polys)
 989            each _split_polygon_at_x(poly, xs[_i])
 990        ], xs, _i=_i+1
 991    );
 992
 993/// Internal Function: _slice_3dpolygons()
 994/// Usage:
 995///   splitpolys = _slice_3dpolygons(polys, dir, cuts);
 996/// Topics: Geometry, Polygons, Intersections
 997/// Description:
 998///   Given a list of 3D polygons, a choice of X, Y, or Z, and a cut list, `cuts`, splits all of the polygons where they cross
 999///   X/Y/Z at any value given in cuts.
1000/// Arguments:
1001///   polys = A list of 3D polygons to split.
1002///   dir_ind = slice direction, 0=X, 1=Y, or 2=Z
1003///   cuts = A list of scalar values for locating the cuts
1004function _slice_3dpolygons(polys, dir, cuts) =
1005    assert( [for (poly=polys) if (!is_path(poly,3)) 1] == [], "Expects list of 3D paths.")
1006    assert( is_vector(cuts), "The split list must be a vector.")
1007    assert( in_list(dir, ["X", "Y", "Z"]))
1008    let(
1009        I = ident(3),
1010        dir_ind = ord(dir)-ord("X")
1011    )
1012    flatten([
1013        for (poly = polys)
1014            if (polygon_area(poly)>EPSILON)   // Discard zero area polygons
1015            let( 
1016                 plane = plane_from_polygon(poly,1e-4))
1017            assert(plane,"Found non-coplanar face.")
1018            let(
1019                normal = point3d(plane),
1020                pnormal = normal - (normal*I[dir_ind])*I[dir_ind]
1021            )
1022            approx(pnormal,[0,0,0]) ? [poly]     // Polygons parallel to cut plane just pass through
1023          : let(
1024                pind = max_index(v_abs(pnormal)),  // project along this direction
1025                otherind = 3-pind-dir_ind,         // keep dir_ind and this direction
1026                keep = [I[dir_ind], I[otherind]],  // dir ind becomes the x dir
1027                poly2d = poly*transpose(keep),     // project to 2d, putting selected direction in the X position
1028                poly_list = [for(p=_split_2dpolygons_at_each_x([poly2d], cuts))
1029                                let(
1030                                    a = p*keep,    // unproject, but pind dimension data is missing
1031                                    ofs = outer_product((repeat(plane[3], len(a))-a*normal)/plane[pind],I[pind])
1032                                 )
1033                                 a+ofs]    // ofs computes the missing pind dimension data and adds it back in
1034            )
1035            poly_list
1036    ]);
1037
1038
1039
1040
1041
1042// Section: Turning a VNF into geometry
1043
1044
1045// Module: vnf_polyhedron()
1046// Synopsis: Returns a polyhedron from a VNF or list of VNFs.
1047// SynTags: Geom
1048// Topics: VNF Manipulation
1049// See Also: vnf_wireframe()
1050// Usage:
1051//   vnf_polyhedron(vnf) [ATTACHMENTS];
1052//   vnf_polyhedron([VNF, VNF, VNF, ...]) [ATTACHMENTS];
1053// Description:
1054//   Given a VNF structure, or a list of VNF structures, creates a polyhedron from them.
1055// Arguments:
1056//   vnf = A VNF structure, or list of VNF structures.
1057//   convexity = Max number of times a line could intersect a wall of the shape.
1058//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
1059//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `"origin"`
1060//   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
1061//   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#subsection-orient).  Default: `UP`
1062//   atype = Select "hull" or "intersect" anchor type.  Default: "hull"
1063// Anchor Types:
1064//   "hull" = Anchors to the virtual convex hull of the shape.
1065//   "intersect" = Anchors to the surface of the shape.
1066// Named Anchors:
1067//   "origin" = Anchor at the origin, oriented UP.
1068module vnf_polyhedron(vnf, convexity=2, cp="centroid", anchor="origin", spin=0, orient=UP, atype="hull") {
1069    vnf = is_vnf_list(vnf)? vnf_join(vnf) : vnf;
1070    assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"");
1071    attachable(anchor,spin,orient, vnf=vnf, extent=atype=="hull", cp=cp) {
1072        polyhedron(vnf[0], vnf[1], convexity=convexity);
1073        children();
1074    }
1075}
1076
1077
1078// Module: vnf_wireframe()
1079// Synopsis: Creates a wireframe model from a VNF.
1080// SynTags: VNF
1081// Topics: VNF Manipulation
1082// See Also: vnf_polyhedron()
1083// Usage:
1084//   vnf_wireframe(vnf, [width]);
1085// Description:
1086//   Given a VNF, creates a wire frame ball-and-stick model of the polyhedron with a cylinder for
1087//   each edge and a sphere at each vertex.  The width parameter specifies the width of the sticks
1088//   that form the wire frame and the diameter of the balls.
1089// Arguments:
1090//   vnf = A vnf structure
1091//   width = width of the cylinders forming the wire frame.  Default: 1
1092// Example:
1093//   $fn=32;
1094//   ball = sphere(r=20, $fn=6);
1095//   vnf_wireframe(ball,width=1);
1096// Example:
1097//   include <BOSL2/polyhedra.scad>
1098//   $fn=32;
1099//   cube_oct = regular_polyhedron_info("vnf",
1100//                      name="cuboctahedron", or=20);
1101//   vnf_wireframe(cube_oct);
1102// Example: The spheres at the vertex are imperfect at aligning with the cylinders, so especially at low $fn things look prety ugly.  This is normal.
1103//   include <BOSL2/polyhedra.scad>
1104//   $fn=8;
1105//   octahedron = regular_polyhedron_info("vnf",
1106//                         name="octahedron", or=20);
1107//   vnf_wireframe(octahedron,width=5);
1108module vnf_wireframe(vnf, width=1)
1109{
1110  no_children($children);
1111  vertex = vnf[0];
1112  edges = unique([for (face=vnf[1], i=idx(face))
1113                    sort([face[i], select(face,i+1)])
1114                 ]);
1115  attachable()
1116  {
1117    union(){
1118      for (e=edges) extrude_from_to(vertex[e[0]],vertex[e[1]]) circle(d=width);
1119      // Identify vertices actually used and draw them
1120      vertused = search(count(len(vertex)), flatten(edges), 1);
1121      for(i=idx(vertex)) if(vertused[i]!=[]) move(vertex[i]) sphere(d=width);
1122    }
1123    union();
1124  }
1125}
1126
1127
1128// Section: Operations on VNFs
1129
1130// Function: vnf_volume()
1131// Synopsis: Returns the volume of a VNF.
1132// Topics: VNF Manipulation
1133// See Also: vnf_area(), vnf_halfspace(), vnf_bend() 
1134// Usage:
1135//   vol = vnf_volume(vnf);
1136// Description:
1137//   Returns the volume enclosed by the given manifold VNF.   The VNF must describe a valid polyhedron with consistent face direction and
1138//   no holes; otherwise the results are undefined.  Returns a positive volume if face direction is clockwise and a negative volume
1139//   if face direction is counter-clockwise.
1140
1141// Divide the polyhedron into tetrahedra with the origin as one vertex and sum up the signed volume.
1142function vnf_volume(vnf) =
1143    let(verts = vnf[0])
1144    sum([
1145         for(face=vnf[1], j=[1:1:len(face)-2])
1146             cross(verts[face[j+1]], verts[face[j]]) * verts[face[0]]
1147    ])/6;
1148
1149
1150// Function: vnf_area()
1151// Synopsis: Returns the surface area of a VNF.
1152// Topics: VNF Manipulation
1153// See Also: vnf_volume(), vnf_halfspace(), vnf_bend() 
1154// Usage:
1155//   area = vnf_area(vnf);
1156// Description:
1157//   Returns the surface area in any VNF by adding up the area of all its faces.  The VNF need not be a manifold.
1158function vnf_area(vnf) =
1159    let(verts=vnf[0])
1160    sum([for(face=vnf[1]) polygon_area(select(verts,face))]);
1161
1162
1163/// Internal Function: _vnf_centroid()
1164/// Usage:
1165///   vol = _vnf_centroid(vnf);
1166/// Description:
1167///   Returns the centroid of the given manifold VNF.  The VNF must describe a valid polyhedron with consistent face direction and
1168///   no holes; otherwise the results are undefined.
1169
1170/// Divide the solid up into tetrahedra with the origin as one vertex.
1171/// The centroid of a tetrahedron is the average of its vertices.
1172/// The centroid of the total is the volume weighted average.
1173function _vnf_centroid(vnf,eps=EPSILON) =
1174    assert(is_vnf(vnf) && len(vnf[0])!=0 && len(vnf[1])!=0,"Invalid or empty VNF given to centroid")
1175    let(
1176        verts = vnf[0],
1177        pos = sum([
1178            for(face=vnf[1], j=[1:1:len(face)-2]) let(
1179                v0  = verts[face[0]],
1180                v1  = verts[face[j]],
1181                v2  = verts[face[j+1]],
1182                vol = cross(v2,v1)*v0
1183            )
1184            [ vol, (v0+v1+v2)*vol ]
1185        ])
1186    )
1187    assert(!approx(pos[0],0, eps), "The vnf has self-intersections.")
1188    pos[1]/pos[0]/4;
1189
1190
1191// Function: vnf_halfspace()
1192// Synopsis: Returns the intersection of the vnf with a half space.
1193// SynTags: VNF
1194// Topics: VNF Manipulation
1195// See Also: vnf_volume(), vnf_area(), vnf_bend() 
1196// Usage:
1197//   newvnf = vnf_halfspace(plane, vnf, [closed], [boundary]);
1198// Description:
1199//   Returns the intersection of the vnf with a half space.  The half space is defined by
1200//   plane = [A,B,C,D], taking the side where the normal [A,B,C] points: Ax+By+Cz≥D.
1201//   If closed is set to false then the cut face is not included in the vnf.  This could
1202//   allow further extension of the vnf by join with other vnfs using {{vnf_join()}}.
1203//   Note that if your given VNF has holes (missing faces) or is not a complete polyhedron
1204//   then closed=true is may produce invalid results when it tries to construct closing faces
1205//   on the cut plane.  Set closed=false for such inputs.
1206//   .
1207//   If you set boundary to true then the return will be the pair [vnf,boundary] where vnf is the
1208//   vnf as usual (with closed=false) and boundary is a list giving each connected component of the cut
1209//   boundary surface.  Each entry in boundary is a list of index values that index into the vnf vertex list (vnf[0]).
1210//   This makes it possible to construct mating shapes, e.g. with {{skin()}} or {{vnf_vertex_array()}} that
1211//   can be combined using {{vnf_join()}} to make a valid polyhedron.
1212//   .
1213//   Note that the input to vnf_halfspace() does not need to be a closed, manifold polyhedron.
1214//   Because it adds the faces on the cut surface, you can use vnf_halfspace() to cap off an open shape if you
1215//   slice through a region that excludes all of the gaps in the input VNF.  
1216// Arguments:
1217//   plane = plane defining the boundary of the half space
1218//   vnf = VNF to cut
1219//   closed = if false do not return the cut face(s) in the returned VNF.  Default: true
1220//   boundary = if true return a pair [vnf,boundary] where boundary is a list of paths on the cut boundary indexed into the VNF vertex list.  If boundary is true, then closed is set to false.  Default: false
1221// Example(3D):
1222//   vnf = cube(10,center=true);
1223//   cutvnf = vnf_halfspace([-1,1,-1,0], vnf);
1224//   vnf_polyhedron(cutvnf);
1225// Example(3D):  Cut face has 2 components
1226//   vnf = path_sweep(circle(r=4, $fn=16),
1227//                    circle(r=20, $fn=64),closed=true);
1228//   cutvnf = vnf_halfspace([-1,1,-4,0], vnf);
1229//   vnf_polyhedron(cutvnf);
1230// Example(3D): Cut face is not simply connected
1231//   vnf = path_sweep(circle(r=4, $fn=16),
1232//                    circle(r=20, $fn=64),closed=true);
1233//   cutvnf = vnf_halfspace([0,0.7,-4,0], vnf);
1234//   vnf_polyhedron(cutvnf);
1235// Example(3D): Cut object has multiple components
1236//   function knot(a,b,t) =   // rolling knot
1237//        [ a * cos (3 * t) / (1 - b* sin (2 *t)),
1238//          a * sin( 3 * t) / (1 - b* sin (2 *t)),
1239//        1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))];
1240//   a = 0.8; b = sqrt (1 - a * a);
1241//   ksteps = 400;
1242//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
1243//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1244//   knot=path_sweep(ushape, knot_path, closed=true, method="incremental");
1245//   cut_knot = vnf_halfspace([1,0,0,0], knot);
1246//   vnf_polyhedron(cut_knot);
1247// Example(VPR=[80,0,15]): Cut a sphere with an arbitrary plane
1248//   vnf1=sphere(r=50, style="icosa", $fn=16);
1249//   vnf2=vnf_halfspace([.8,1,-1.5,0], vnf1);
1250//   vnf_polyhedron(vnf2);
1251// Example(VPR=[80,0,15]): Cut it again, but with closed=false to leave an open boundary. 
1252//   vnf1=sphere(r=50, style="icosa", $fn=16);
1253//   vnf2=vnf_halfspace([.8,1,-1.5,0], vnf1);
1254//   vnf3=vnf_halfspace([0,0,-1,0], vnf2, closed=false);
1255//   vnf_polyhedron(vnf3);
1256// Example(VPR=[80,0,15]): Use {vnf_join()} to combine with a mating vnf, in this case a reflection of the part we made. 
1257//   vnf1=sphere(r=50, style="icosa", $fn=16);
1258//   vnf2=vnf_halfspace([.8,1,-1.5,0], vnf1);
1259//   vnf3=vnf_halfspace([0,0,-1,0], vnf2, closed=false);
1260//   vnf4=vnf_join([vnf3, zflip(vnf3,1)]);
1261//   vnf_polyhedron(vnf4);
1262// Example: When the input VNF is a surface with a boundary, if you use the default setting closed=true, then vnf_halfspace() tries to construct closing faces from the edges created by the cut.  These faces may be invalid, for example if the cut points are collinear.  In this example the constructed face is a valid face.
1263//   include <BOSL2/beziers.scad>
1264//   patch=[
1265//          [[10,-10,0],[1,-1,0],[-1,-1,0],[-10,-10,0]],
1266//          [[10,-10,20],[1,-1,20],[-1,-1,20],[-10,-10,20]]
1267//         ];
1268//   vnf=bezier_vnf(patch);
1269//   vnfcut = vnf_halfspace([-.8,0,-1,-14],vnf);
1270//   vnf_polyhedron(vnfcut);
1271// Example: Setting closed to false eliminates this (possibly invalid) face:
1272//   include <BOSL2/beziers.scad>
1273//   patch=[
1274//          [[10,-10,0],[1,-1,0],[-1,-1,0],[-10,-10,0]],
1275//          [[10,-10,20],[1,-1,20],[-1,-1,20],[-10,-10,20]]
1276//         ];
1277//   vnf=bezier_vnf(patch);
1278//   vnfcut = vnf_halfspace([-.8,0,-1,-14],vnf,closed=false);
1279//   vnf_polyhedron(vnfcut);
1280// Example: Here is a VNF that has holes, so it is not a valid manifold. 
1281//   outside = linear_sweep(circle(r=30), h=100, caps=false);
1282//   inside = yrot(7,linear_sweep(circle(r=10), h=120, caps=false));
1283//   open_vnf=vnf_join([outside, vnf_reverse_faces(inside)]);
1284//   vnf_polyhedron(open_vnf);
1285// Example: By cutting it at each end we can create closing faces, resulting in a valid manifold without holes.
1286//   outside = linear_sweep(circle(r=30), h=100, caps=false);
1287//   inside = yrot(11,linear_sweep(circle(r=10), h=120, caps=false));
1288//   open_vnf=vnf_join([outside, vnf_reverse_faces(inside)]);
1289//   vnf = vnf_halfspace([0,0,1,5], vnf_halfspace([0,.7,-1,-75], open_vnf));
1290//   vnf_polyhedron(vnf);
1291// Example: If boundary=true then the return is a list with the VNF and boundary data.  
1292//   vnf = path_sweep(circle(r=4, $fn=16),
1293//                    circle(r=20, $fn=64),closed=true);
1294//   cut_bnd = vnf_halfspace([-1,1,-4,0], vnf, boundary=true);*/
1295//   cutvnf = cut_bnd[0];
1296//   boundary = [for(b=cut_bnd[1]) select(cutvnf[0],b)];
1297//   vnf_polyhedron(cutvnf);
1298//   stroke(boundary,color="red");
1299function vnf_halfspace(plane, vnf, closed=true, boundary=false) =
1300    assert(_valid_plane(plane), "Invalid plane")
1301    assert(is_vnf(vnf), "Invalid vnf")
1302    let(
1303         inside = [for(x=vnf[0]) plane*[each x,-1] >= -EPSILON ? 1 : 0],
1304         vertexmap = [0,each cumsum(inside)],
1305         faces_edges_vertices = _vnfcut(plane, vnf[0],vertexmap,inside, vnf[1], last(vertexmap)),
1306         newvert = concat(bselect(vnf[0],inside), faces_edges_vertices[2])
1307    )
1308    closed==false && !boundary ? [newvert, faces_edges_vertices[0]]
1309  : let(
1310        allpaths = _assemble_paths(newvert, faces_edges_vertices[1]),
1311        newpaths = [for(p=allpaths) if (len(p)>=3) p
1312                                    else assert(approx(p[0],p[1]),"Orphan edge found when assembling cut edges.")
1313           ]
1314    )
1315    boundary ? [[newvert, faces_edges_vertices[0]], newpaths]
1316  : len(newpaths)<=1 ? [newvert, concat(faces_edges_vertices[0], newpaths)]
1317  : let(
1318           M = project_plane(plane),
1319           faceregion = [for(path=newpaths) path2d(apply(M,select(newvert,path)))],
1320           facevnf = vnf_from_region(faceregion,transform=rot_inverse(M),reverse=true)
1321      )
1322      vnf_join([[newvert, faces_edges_vertices[0]], facevnf]);
1323
1324function _assemble_paths(vertices, edges, paths=[],i=0) =
1325     i==len(edges) ? paths :
1326     norm(vertices[edges[i][0]]-vertices[edges[i][1]])<EPSILON ? _assemble_paths(vertices,edges,paths,i+1) :
1327     let(    // Find paths that connects on left side and right side of the edges (if one exists)
1328         left = [for(j=idx(paths)) if (approx(vertices[last(paths[j])],vertices[edges[i][0]])) j],
1329         right = [for(j=idx(paths)) if (approx(vertices[edges[i][1]],vertices[paths[j][0]])) j]
1330     )
1331     assert(len(left)<=1 && len(right)<=1)
1332     let(
1333          keep_path = list_remove(paths,concat(left,right)),
1334          update_path = left==[] && right==[] ? edges[i]
1335                      : left==[] ? concat([edges[i][0]],paths[right[0]])
1336                      : right==[] ?  concat(paths[left[0]],[edges[i][1]])
1337                      : left != right ? concat(paths[left[0]], paths[right[0]])
1338                      : paths[left[0]]
1339     )
1340     _assemble_paths(vertices, edges, concat(keep_path, [update_path]), i+1);
1341
1342
1343function _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount, newfaces=[], newedges=[], newvertices=[], i=0) =
1344   i==len(faces) ? [newfaces, newedges, newvertices] :
1345   let(
1346        pts_inside = select(inside,faces[i])
1347   )
1348   all(pts_inside) ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount,
1349                             concat(newfaces, [select(vertexmap,faces[i])]), newedges, newvertices, i+1):
1350   !any(pts_inside) ? _vnfcut(plane, vertices, vertexmap,inside, faces, vertcount, newfaces, newedges, newvertices, i+1):
1351   let(
1352        first = search([[1,0]],pair(pts_inside,wrap=true),0)[0],
1353        second = search([[0,1]],pair(pts_inside,wrap=true),0)[0]
1354   )
1355   assert(len(first)==1 && len(second)==1, "Found concave face in VNF.  Run vnf_triangulate first to ensure convex faces.")
1356   let(
1357        newface = [each select(vertexmap,select(faces[i],second[0]+1,first[0])),vertcount, vertcount+1],
1358        newvert = [plane_line_intersection(plane, select(vertices,select(faces[i],first[0],first[0]+1)),eps=0),
1359                   plane_line_intersection(plane, select(vertices,select(faces[i],second[0],second[0]+1)),eps=0)]
1360   )
1361   true //!approx(newvert[0],newvert[1])
1362       ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount+2,
1363                 concat(newfaces, [newface]), concat(newedges,[[vertcount+1,vertcount]]),concat(newvertices,newvert),i+1)
1364   :len(newface)>3
1365       ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount+1,
1366                 concat(newfaces, [list_head(newface)]), newedges,concat(newvertices,[newvert[0]]),i+1)
1367   :
1368   _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount,newfaces, newedges, newvert, i+1);
1369
1370
1371
1372
1373function _triangulate_planar_convex_polygons(polys) =
1374    polys==[]? [] :
1375    let(
1376        tris = [for (poly=polys) if (len(poly)==3) poly],
1377        bigs = [for (poly=polys) if (len(poly)>3) poly],
1378        newtris = [for (poly=bigs) select(poly,-2,0)],
1379        newbigs = [for (poly=bigs) select(poly,0,-2)],
1380        newtris2 = _triangulate_planar_convex_polygons(newbigs),
1381        outtris = concat(tris, newtris, newtris2)
1382    ) outtris;
1383
1384//**
1385// this function may produce degenerate triangles:
1386//    _triangulate_planar_convex_polygons([ [for(i=[0:1]) [i,i],
1387//                                           [1,-1], [-1,-1],
1388//                                           for(i=[-1:0]) [i,i] ] ] )
1389//    == [[[-1, -1], [ 0,  0], [0,  0]]
1390//        [[-1, -1], [-1, -1], [0,  0]]
1391//        [[ 1, -1], [-1, -1], [0,  0]]
1392//        [[ 0,  0], [ 1,  1], [1, -1]] ]
1393//
1394
1395// Function: vnf_bend()
1396// Synopsis: Bends a VNF around an axis.
1397// SynTags: VNF
1398// Topics: VNF Manipulation
1399// See Also: vnf_volume(), vnf_area(), vnf_halfspace() 
1400// Usage:
1401//   bentvnf = vnf_bend(vnf,r|d=,[axis=]);
1402// Description:
1403//   Bend a VNF around the X, Y or Z axis, splitting up faces as necessary.  Returns the bent
1404//   VNF.  For bending around the Z axis the input VNF must not cross the Y=0 plane.  For bending
1405//   around the X or Y axes the VNF must not cross the Z=0 plane.  Note that if you wrap a VNF all the way around
1406//   it may intersect itself, which produces an invalid polyhedron.  It is your responsibility to
1407//   avoid this situation.  The 1:1
1408//   radius is where the curved length of the bent VNF matches the length of the original VNF.  If the
1409//   `r` or `d` arguments are given, then they will specify the 1:1 radius or diameter.  If they are
1410//   not given, then the 1:1 radius will be defined by the distance of the furthest vertex in the
1411//   original VNF from the Z=0 plane.  You can adjust the granularity of the bend using the standard
1412//   `$fa`, `$fs`, and `$fn` variables.
1413// Arguments:
1414//   vnf = The original VNF to bend.
1415//   r = If given, the radius where the size of the original shape is the same as in the original.
1416//   ---
1417//   d = If given, the diameter where the size of the original shape is the same as in the original.
1418//   axis = The axis to wrap around.  "X", "Y", or "Z".  Default: "Z"
1419// Example(3D):
1420//   vnf0 = cube([100,40,10], center=true);
1421//   vnf1 = up(50, p=vnf0);
1422//   vnf2 = down(50, p=vnf0);
1423//   bent1 = vnf_bend(vnf1, axis="Y");
1424//   bent2 = vnf_bend(vnf2, axis="Y");
1425//   vnf_polyhedron([bent1,bent2]);
1426// Example(3D):
1427//   vnf0 = linear_sweep(star(n=5,step=2,d=100), height=10);
1428//   vnf1 = up(50, p=vnf0);
1429//   vnf2 = down(50, p=vnf0);
1430//   bent1 = vnf_bend(vnf1, axis="Y");
1431//   bent2 = vnf_bend(vnf2, axis="Y");
1432//   vnf_polyhedron([bent1,bent2]);
1433// Example(3D):
1434//   rgn = union(rect([100,20]),
1435//               rect([20,100]));
1436//   vnf0 = linear_sweep(zrot(45,p=rgn), height=10);
1437//   vnf1 = up(50, p=vnf0);
1438//   vnf2 = down(50, p=vnf0);
1439//   bent1 = vnf_bend(vnf1, axis="Y");
1440//   bent2 = vnf_bend(vnf2, axis="Y");
1441//   vnf_polyhedron([bent1,bent2]);
1442// Example(3D): Bending Around X Axis.
1443//   rgnr = union(
1444//       rect([20,100]),
1445//       back(50, p=trapezoid(w1=40, w2=0, h=20, anchor=FRONT))
1446//   );
1447//   vnf0 = xrot(00,p=linear_sweep(rgnr, height=10));
1448//   vnf1 = up(50, p=vnf0);
1449//   #vnf_polyhedron(vnf1);
1450//   bent1 = vnf_bend(vnf1, axis="X");
1451//   vnf_polyhedron([bent1]);
1452// Example(3D): Bending Around Y Axis.
1453//   rgn = union(
1454//       rect([20,100]),
1455//       back(50, p=trapezoid(w1=40, w2=0, h=20, anchor=FRONT))
1456//   );
1457//   rgnr = zrot(-90, p=rgn);
1458//   vnf0 = xrot(00,p=linear_sweep(rgnr, height=10));
1459//   vnf1 = up(50, p=vnf0);
1460//   #vnf_polyhedron(vnf1);
1461//   bent1 = vnf_bend(vnf1, axis="Y");
1462//   vnf_polyhedron([bent1]);
1463// Example(3D): Bending Around Z Axis.
1464//   rgn = union(
1465//       rect([20,100]),
1466//       back(50, p=trapezoid(w1=40, w2=0, h=20, anchor=FRONT))
1467//   );
1468//   rgnr = zrot(90, p=rgn);
1469//   vnf0 = xrot(90,p=linear_sweep(rgnr, height=10));
1470//   vnf1 = fwd(50, p=vnf0);
1471//   #vnf_polyhedron(vnf1);
1472//   bent1 = vnf_bend(vnf1, axis="Z");
1473//   vnf_polyhedron([bent1]);
1474// Example(3D): Bending more than once around the cylinder
1475//   $fn=32;
1476//   vnf = apply(fwd(5)*yrot(30),cube([100,2,5],center=true));
1477//   bent = vnf_bend(vnf, axis="Z");
1478//   vnf_polyhedron(bent);
1479function vnf_bend(vnf,r,d,axis="Z") =
1480    let(
1481        chk_axis = assert(in_list(axis,["X","Y","Z"])),
1482        verts = vnf[0],
1483        bounds = pointlist_bounds(verts),
1484        bmin = bounds[0],
1485        bmax = bounds[1],
1486        dflt = axis=="Z"?
1487            max(abs(bmax.y), abs(bmin.y)) :
1488            max(abs(bmax.z), abs(bmin.z)),
1489        r = get_radius(r=r,d=d,dflt=dflt),
1490        extent = axis=="X" ? [bmin.y, bmax.y] : [bmin.x, bmax.x]
1491    )
1492    let(
1493        span_chk = axis=="Z"?
1494            assert(bmin.y > 0 || bmax.y < 0, "Entire shape MUST be completely in front of or behind y=0.") :
1495            assert(bmin.z > 0 || bmax.z < 0, "Entire shape MUST be completely above or below z=0."),
1496        steps = 1+ceil(segs(r) * (extent[1]-extent[0])/(2*PI*r)),
1497        step = (extent[1]-extent[0]) / steps,
1498        bend_at = [for(i = [1:1:steps-1]) i*step+extent[0]],
1499        slicedir = axis=="X"? "Y" : "X",   // slice in y dir for X axis case, and x dir otherwise
1500        sliced = vnf_slice(vnf, slicedir, bend_at),
1501        coord = axis=="X" ? [0,sign(bmax.z),0] : axis=="Y" ? [sign(bmax.z),0,0] : [sign(bmax.y),0,0],
1502        new_vert = [for(p=sliced[0])
1503                       let(a=coord*p*180/(PI*r))
1504                       axis=="X"? [p.x, p.z*sin(a), p.z*cos(a)] :
1505                       axis=="Y"? [p.z*sin(a), p.y, p.z*cos(a)] :
1506                       [p.y*sin(a), p.y*cos(a), p.z]]
1507   ) [new_vert,sliced[1]];
1508
1509
1510
1511// Function&Module: vnf_hull()
1512// Synopsis: Compute convex hull of VNF or 3d path
1513// Usage:
1514//    vnf_hull = hull_vnf(vnf);
1515//    hull_vnf(vnf,[fast]);
1516// Description:
1517//   Given a VNF or a list of 3d points, compute the convex hull
1518//   and return it as a VNF.  This differs from {{hull()}} and {{hull3d_faces()}} which
1519//   return just the face list referenced to the input point list.  Note that the point
1520//   list that is returned will contain all the points that are actually used in the input
1521//   VNF, which may be many more points than are needed to represent the convex hull.
1522//   This is not usually a problem, but you can run the somewhat slow {{vnf_drop_unused_points()}}
1523//   function to fix this if necessary.  
1524// Arguments:
1525//   region = region or path listing points to compute the hull from.
1526//   fast = (module only) if input is a point list (not a VNF) use a fasterer cheat that may handle more points, but could emit warnings.  Ignored if input is a VNF.  Default: false.  
1527// Example(3D,Big,NoAxes,VPR=[55,0,25],VPT=[9.47096,-4.50217,8.45727],VPD=60.2654): Input is a VNF
1528//   ellipse = xscale(2, p=circle($fn=48, r=3));
1529//   pentagon = subdivide_path(pentagon(r=1), 20);
1530//   vnf=path_sweep(pentagon, path3d(ellipse),
1531//                  closed=true, twist=360*2);
1532//   vnfhull = vnf_hull(vnf);
1533//   vnf_polyhedron(vnf);
1534//   move([10,10])
1535//     vnf_polyhedron(vnfhull);
1536// Example(3D,Med,NoAxes,VPR=[70.4,0,110.4],VPT=[5.97456,1.26459,18.0317],VPD=126): Input is a point list
1537//   h=helix(l=40, turns=1, r=8);
1538//   color("red")move_copies(h)
1539//     sphere(r=0.5,$fn=12);
1540//   vnf_polyhedron(vnf_hull(h));
1541function vnf_hull(vnf) =
1542  assert(is_vnf(vnf) || is_path(vnf,3),"Input must be a VNF or a 3d path")
1543  let(
1544      pts = is_vnf(vnf) ? select(vnf[0],unique(flatten(vnf[1])))
1545                        : vnf,
1546      faces = hull3d_faces(pts)
1547  )
1548  [pts, faces];
1549
1550module vnf_hull(vnf, fast=false)
1551{
1552  if (is_vnf(vnf)) hull()vnf_polyhedron(vnf);
1553  else hull_points(vnf, fast);
1554}  
1555    
1556
1557// Section: Debugging Polyhedrons
1558
1559/// Internal Module: _show_vertices()
1560/// Usage:
1561///   _show_vertices(vertices, [size], [filter=])
1562/// Description:
1563///   Draws all the vertices in an array, at their 3D position, numbered by their
1564///   position in the vertex array.  Also draws any children of this module with
1565///   transparency.
1566/// Arguments:
1567///   vertices = Array of point vertices.
1568///   size = The size of the text used to label the vertices.  Default: 1
1569/// Example:
1570///   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
1571///   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
1572///   _show_vertices(vertices=verts, size=2) {
1573///       polyhedron(points=verts, faces=faces);
1574///   }
1575module _show_vertices(vertices, size=1, filter) {
1576    color("blue") {
1577        dups = vector_search(vertices, EPSILON, vertices);
1578        for (ind = dups) {
1579            if (is_undef(filter) || any(ind, filter)) {
1580                numstr = str_join([for(i=ind) str(i)],",");
1581                v = vertices[ind[0]];
1582                translate(v) {
1583                    rot($vpr) back(size/8){
1584                       linear_extrude(height=size/10, center=true, convexity=10) {
1585                          text(text=numstr, size=size, halign="center");
1586                       }
1587                    }
1588                    sphere(size/10);
1589                }
1590            }
1591        }
1592    }
1593}
1594
1595
1596/// Internal Module: _show_faces()
1597/// Usage:
1598///   _show_faces(vertices, faces, [size=], [filter=]);
1599/// Description:
1600///   Draws all the vertices at their 3D position, numbered in blue by their
1601///   position in the vertex array.  Each face will have their face number drawn
1602///   in red, aligned with the center of face.  All children of this module are drawn
1603///   with transparency.
1604/// Arguments:
1605///   vertices = Array of point vertices.
1606///   faces = Array of faces by vertex numbers.
1607///   size = The size of the text used to label the faces and vertices.  Default: 1
1608/// Example(EdgesMed):
1609///   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
1610///   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
1611///   _show_faces(vertices=verts, faces=faces, size=2) {
1612///       polyhedron(points=verts, faces=faces);
1613///   }
1614module _show_faces(vertices, faces, size=1, filter) {
1615    vlen = len(vertices);
1616    color("red") {
1617        for (i = [0:1:len(faces)-1]) {
1618            face = faces[i];
1619            if (face[0] < 0 || face[1] < 0 || face[2] < 0 || face[0] >= vlen || face[1] >= vlen || face[2] >= vlen) {
1620                echo(str("INVALID FACE: indices of face ",i," are out of bounds [0,",vlen-1,"]: face=",face));
1621            }
1622            else if (is_undef(filter) || any(face,filter)) {
1623                verts = select(vertices,face);
1624                normal = polygon_normal(verts);
1625                if (is_undef(normal))
1626                    echo(str("DEGENERATE FACE: face ",i," has no normal vector, face=", face));
1627                else {
1628                    axis = vector_axis(normal, DOWN);
1629                    ang = vector_angle(normal, DOWN);
1630                    theta = atan2(normal[1], normal[0]);
1631                    translate(mean(verts)) 
1632                      rotate(a=(180-ang), v=axis)
1633                      zrot(theta+90)
1634                      linear_extrude(height=size/10, center=true, convexity=10) {
1635                                text(text=str(i), size=size, halign="center");
1636                                text(text=str("_"), size=size, halign="center");
1637                      }
1638                }
1639            }
1640        }
1641    }        
1642}
1643
1644
1645
1646// Module: debug_vnf()
1647// Synopsis: A replacement for `vnf_polyhedron()` to help with debugging.
1648// SynTags: VNF
1649// Topics: VNF Manipulation, Debugging
1650// See Also: vnf_validate()
1651// Usage:
1652//   debug_vnf(vnfs, [faces=], [vertices=], [opacity=], [size=], [convexity=], [filter=]);
1653// Description:
1654//   A drop-in module to replace `vnf_polyhedron()` to help debug vertices and faces.
1655//   Draws all the vertices at their 3D position, numbered in blue by their
1656//   position in the vertex array.  Each face will have its face number drawn
1657//   in red, aligned with the center of face.  All given faces are drawn with
1658//   transparency. All children of this module are drawn with transparency.
1659//   Works best with Thrown-Together preview mode, to see reversed faces.
1660//   You can set opacity to 0 if you want to supress the display of the polyhedron faces.
1661//   .
1662//   The vertex numbers are shown rotated to face you.  As you rotate your polyhedron you
1663//   can rerun the preview to display them oriented for viewing from a different viewpoint.
1664// Topics: Polyhedra, Debugging
1665// Arguments:
1666//   vnf = VNF to display
1667//   ---
1668//   faces = if true display face numbers.  Default: true
1669//   vertices = if true display vertex numbers.  Default: true
1670//   opacity = Opacity of the polyhedron faces.  Default: 0.5
1671//   convexity = The max number of walls a ray can pass through the given polygon paths.
1672//   size = The size of the text used to label the faces and vertices.  Default: 1
1673//   filter = If given a function literal of signature `function(i)`, will only show labels for vertices and faces that have a vertex index that gets a true result from that function.  Default: no filter.
1674// Example(EdgesMed):
1675//   verts = [for (z=[-10,10], a=[0:120:359.9]) [10*cos(a),10*sin(a),z]];
1676//   faces = [[0,1,2], [5,4,3], [0,3,4], [0,4,1], [1,4,5], [1,5,2], [2,5,3], [2,3,0]];
1677//   debug_vnf([verts,faces], size=2);
1678module debug_vnf(vnf, faces=true, vertices=true, opacity=0.5, size=1, convexity=6, filter ) {
1679    no_children($children);
1680    if (faces)
1681      _show_faces(vertices=vnf[0], faces=vnf[1], size=size, filter=filter);
1682    if (vertices)
1683      _show_vertices(vertices=vnf[0], size=size, filter=filter);
1684    if (opacity > 0)
1685      color([0.2, 1.0, 0, opacity])
1686        vnf_polyhedron(vnf,convexity=convexity);
1687}
1688
1689
1690// Module: vnf_validate()
1691// Synopsis: Echos non-manifold VNF errors to the console.
1692// SynTags: VNF
1693// Topics: VNF Manipulation, Debugging
1694// See Also: debug_vnf()
1695// 
1696// Usage: 
1697//   vnf_validate(vnf, [size], [show_warns=], [check_isects=], [opacity=], [adjacent=], [label_verts=], [label_faces=], [wireframe=]);
1698// Description:
1699//   When called as a module, echoes the non-manifold errors to the console, and color hilites the
1700//   bad edges and vertices, overlaid on a transparent gray polyhedron of the VNF.
1701//   .
1702//   Currently checks for these problems:
1703//   .
1704//   Type    | Color    | Code         | Message
1705//   ------- | -------- | ------------ | ---------------------------------
1706//   WARNING | Yellow   | BIG_FACE     | Face has more than 3 vertices, and may confuse CGAL.
1707//   WARNING | Blue     | NULL_FACE    | Face has zero area.
1708//   ERROR   | Cyan     | NONPLANAR    | Face vertices are not coplanar.
1709//   ERROR   | Brown    | DUP_FACE     | Multiple instances of the same face.
1710//   ERROR   | Orange   | MULTCONN     | Multiply Connected Geometry. Too many faces attached at Edge.
1711//   ERROR   | Violet   | REVERSAL     | Faces reverse across edge.
1712//   ERROR   | Red      | T_JUNCTION   | Vertex is mid-edge on another Face.
1713//   ERROR   | Brown    | FACE_ISECT   | Faces intersect.
1714//   ERROR   | Magenta  | HOLE_EDGE    | Edge bounds Hole.
1715//   .
1716//   Still to implement:
1717//   - Overlapping coplanar faces.
1718// Arguments:
1719//   vnf = The VNF to validate.
1720//   size = The width of the lines and diameter of points used to highlight edges and vertices.  Module only.  Default: 1
1721//   ---
1722//   show_warns = If true show warnings for non-triangular faces.  Default: true
1723//   check_isects = If true, performs slow checks for intersecting faces.  Default: false
1724//   opacity = The opacity level to show the polyhedron itself with.    Default: 0.67
1725//   label_verts = If true, shows labels at each vertex that show the vertex number.    Default: false
1726//   label_faces = If true, shows labels at the center of each face that show the face number.    Default: false
1727//   wireframe = If true, shows edges more clearly so you can see them in Thrown Together mode.    Default: false
1728//   adjacent = If true, only display faces adjacent to a vertex listed in the errors.    Default: false
1729// Example(3D,Edges): BIG_FACE Warnings; Faces with More Than 3 Vertices.  CGAL often will fail to accept that a face is planar after a rotation, if it has more than 3 vertices.
1730//   vnf = skin([
1731//       path3d(regular_ngon(n=3, d=100),0),
1732//       path3d(regular_ngon(n=5, d=100),100)
1733//   ], slices=0, caps=true, method="tangent");
1734//   vnf_validate(vnf);
1735// Example(3D,Edges): NONPLANAR Errors; Face Vertices are Not Coplanar
1736//   a = [  0,  0,-50];
1737//   b = [-50,-50, 50];
1738//   c = [-50, 50, 50];
1739//   d = [ 50, 50, 60];
1740//   e = [ 50,-50, 50];
1741//   vnf = vnf_from_polygons([
1742//       [a, b, e], [a, c, b], [a, d, c], [a, e, d], [b, c, d, e]
1743//   ],fast=true);
1744//   vnf_validate(vnf);
1745// Example(3D,Edges): MULTCONN Errors; More Than Two Faces Attached to the Same Edge.  This confuses CGAL, and can lead to failed renders.
1746//   vnf = vnf_triangulate(linear_sweep(union(square(50), square(50,anchor=BACK+RIGHT)), height=50));
1747//   vnf_validate(vnf);
1748// Example(3D,Edges): REVERSAL Errors; Faces Reversed Across Edge
1749//   vnf1 = skin([
1750//       path3d(square(100,center=true),0),
1751//       path3d(square(100,center=true),100),
1752//   ], slices=0, caps=false);
1753//   vnf = vnf_join([vnf1, vnf_from_polygons([
1754//       [[-50,-50,  0], [ 50, 50,  0], [-50, 50,  0]],
1755//       [[-50,-50,  0], [ 50,-50,  0], [ 50, 50,  0]],
1756//       [[-50,-50,100], [-50, 50,100], [ 50, 50,100]],
1757//       [[-50,-50,100], [ 50,-50,100], [ 50, 50,100]],
1758//   ])]);
1759//   vnf_validate(vnf);
1760// Example(3D,Edges): T_JUNCTION Errors; Vertex is Mid-Edge on Another Face.
1761//   vnf = [
1762//       [
1763//           each path3d(square(100,center=true),0),
1764//           each path3d(square(100,center=true),100),
1765//           [0,-50,100],
1766//       ], [
1767//          [0,2,1], [0,3,2], [0,8,4], [0,1,8], [1,5,8],
1768//          [0,4,3], [4,7,3], [1,2,5], [2,6,5], [3,7,6],
1769//          [3,6,2], [4,5,6], [4,6,7],
1770//       ]
1771//   ];
1772//   vnf_validate(vnf);
1773// Example(3D,Edges): FACE_ISECT Errors; Faces Intersect
1774//   vnf = vnf_join([
1775//       linear_sweep(square(100,center=true), height=100),
1776//       move([75,35,30],p=linear_sweep(square(100,center=true), height=100))
1777//   ]);
1778//   vnf_validate(vnf,size=2,check_isects=true);
1779// Example(3D,Edges): HOLE_EDGE Errors; Edges Adjacent to Holes.
1780//   vnf = skin([
1781//       path3d(regular_ngon(n=4, d=100),0),
1782//       path3d(regular_ngon(n=5, d=100),100)
1783//   ], slices=0, caps=false);
1784//   vnf_validate(vnf,size=2);
1785
1786
1787//   Returns a list of non-manifold errors with the given VNF.
1788//   Each error has the format `[ERR_OR_WARN,CODE,MESG,POINTS,COLOR]`.
1789function _vnf_validate(vnf, show_warns=true, check_isects=false) =
1790    assert(is_vnf(vnf), "Invalid VNF")
1791    let(
1792        varr = vnf[0],
1793        faces = vnf[1],
1794        lvarr = len(varr),
1795        edges = sort([
1796            for (face=faces, edge=pair(face,true))
1797            edge[0]<edge[1]? edge : [edge[1],edge[0]]
1798        ]),
1799        dfaces = [
1800            for (face=faces) let(
1801                face=deduplicate_indexed(varr,face,closed=true)
1802            ) if(len(face)>=3)
1803            face
1804        ],
1805        face_areas = [
1806            for (face = faces)
1807            len(face) < 3? 0 :
1808            polygon_area([for (k=face) varr[k]])
1809        ],
1810        edgecnts = unique_count(edges),
1811        uniq_edges = edgecnts[0],
1812        issues = []
1813    )
1814    let(
1815        big_faces = !show_warns? [] : [
1816            for (face = faces)
1817            if (len(face) > 3)
1818            _vnf_validate_err("BIG_FACE", face)
1819        ],
1820        null_faces = !show_warns? [] : [
1821            for (i = idx(faces)) let(
1822                face = faces[i],
1823                area = face_areas[i]
1824            )
1825            if (is_num(area) && abs(area) < EPSILON)
1826            _vnf_validate_err("NULL_FACE", face)
1827        ],
1828        issues = concat(big_faces, null_faces)
1829    )
1830    let(
1831        bad_indices = [
1832            for (face = faces, idx = face)
1833            if (idx < 0 || idx >= lvarr)
1834            _vnf_validate_err("BAD_INDEX", [idx])
1835        ],
1836        issues = concat(issues, bad_indices)
1837    ) bad_indices? issues :
1838    let(
1839        repeated_faces = [
1840            for (i=idx(dfaces), j=idx(dfaces))
1841            if (i!=j) let(
1842                face1 = dfaces[i],
1843                face2 = dfaces[j]
1844            ) if (min(face1) == min(face2)) let(
1845                min1 = min_index(face1),
1846                min2 = min_index(face2)
1847            ) if (min1 == min2) let(
1848                sface1 = list_rotate(face1,min1),
1849                sface2 = list_rotate(face2,min2)
1850            ) if (sface1 == sface2)
1851            _vnf_validate_err("DUP_FACE", sface1)
1852        ],
1853        issues = concat(issues, repeated_faces)
1854    ) repeated_faces? issues :
1855    let(
1856        multconn_edges = unique([
1857            for (i = idx(uniq_edges))
1858            if (edgecnts[1][i]>2)
1859            _vnf_validate_err("MULTCONN", uniq_edges[i])
1860        ]),
1861        issues = concat(issues, multconn_edges)
1862    ) multconn_edges? issues :
1863    let(
1864        reversals = unique([
1865            for(i = idx(dfaces), j = idx(dfaces)) if(i != j)
1866            for(edge1 = pair(faces[i],true))
1867            for(edge2 = pair(faces[j],true))
1868            if(edge1 == edge2)  // Valid adjacent faces will never have the same vertex ordering.
1869            if(_edge_not_reported(edge1, varr, multconn_edges))
1870            _vnf_validate_err("REVERSAL", edge1)
1871        ]),
1872        issues = concat(issues, reversals)
1873    ) reversals? issues :
1874    let(
1875        t_juncts = unique([
1876            for (v=idx(varr), edge=uniq_edges) let(
1877                ia = edge[0],
1878                ib = v,
1879                ic = edge[1]
1880            )
1881            if (ia!=ib && ib!=ic && ia!=ic) let(
1882                a = varr[ia],
1883                b = varr[ib],
1884                c = varr[ic]
1885            )
1886            if (!approx(a,b) && !approx(b,c) && !approx(a,c)) let(
1887                pt = line_closest_point([a,c],b,SEGMENT)
1888            )
1889            if (approx(pt,b))
1890            _vnf_validate_err("T_JUNCTION", [ib])
1891        ]),
1892        issues = concat(issues, t_juncts)
1893    ) t_juncts? issues :
1894    let(
1895        isect_faces = !check_isects? [] : unique([
1896            for (i = [0:1:len(faces)-2])
1897              let(
1898                  f1 = faces[i],
1899                  poly1   = select(varr, faces[i]),
1900                  plane1  = plane3pt(poly1[0], poly1[1], poly1[2]),
1901                  normal1 = [plane1[0], plane1[1], plane1[2]]
1902              )
1903              for (j = [i+1:1:len(faces)-1])
1904                let(
1905                  f2 = faces[j],
1906                  poly2 = select(varr, f2),
1907                  val = poly2 * normal1
1908                )
1909                // The next test skips f2 if it lies entirely on one side of the plane of poly1
1910                if( min(val)<=plane1[3] && max(val)>=plane1[3] )
1911                  let(
1912                      plane2  = plane_from_polygon(poly2),
1913                      normal2 = [plane2[0], plane2[1], plane2[2]],
1914                      val = poly1 * normal2
1915                  )
1916                  // Skip if f1 lies entirely on one side of the plane defined by poly2
1917                  if( min(val)<=plane2[3] && max(val)>=plane2[3] )
1918                    let(
1919                        shared_edges = [
1920                                        for (edge1 = pair(f1, true), edge2 = pair(f2, true))
1921                                           if (edge1 == [edge2[1], edge2[0]]) 1
1922                                       ]
1923                    )
1924                    if (shared_edges==[])
1925                       let(
1926                           line = plane_intersection(plane1, plane2)
1927                       )
1928                       if (is_def(line))
1929                          let(
1930                              isects = polygon_line_intersection(poly1, line)
1931                          )
1932                          if (is_def(isects))
1933                            for (isect = isects)
1934                              if (len(isect) > 1)
1935                                let(
1936                                    isects2 = polygon_line_intersection(poly2, isect, bounded=true)
1937                                )
1938                                if (is_def(isects2))
1939                                  for (seg = isects2) 
1940                                    if (len(seg)>1 && seg[0] != seg[1]) _vnf_validate_err("FACE_ISECT", seg)
1941        ]),
1942        issues = concat(issues, isect_faces)
1943    ) isect_faces? issues :
1944    let(
1945        hole_edges = unique([
1946            for (i=idx(uniq_edges))
1947            if (edgecnts[1][i]<2)
1948            if (_pts_not_reported(uniq_edges[i], varr, t_juncts))
1949            if (_pts_not_reported(uniq_edges[i], varr, isect_faces))
1950            _vnf_validate_err("HOLE_EDGE", uniq_edges[i])
1951        ]),
1952        issues = concat(issues, hole_edges)
1953    ) hole_edges? issues :
1954    let(
1955        nonplanars = unique([
1956            for (i = idx(faces))
1957               if (is_undef(face_areas[i])) 
1958                  _vnf_validate_err("NONPLANAR", faces[i])
1959        ]),
1960        issues = concat(issues, nonplanars)
1961    ) issues;
1962
1963
1964_vnf_validate_errs = [
1965    ["BIG_FACE",    "WARNING", "cyan",    "Face has more than 3 vertices, and may confuse CGAL"],
1966    ["NULL_FACE",   "WARNING", "blue",    "Face has zero area."],
1967    ["BAD_INDEX",   "ERROR",   "cyan",    "Invalid face vertex index."],
1968    ["NONPLANAR",   "ERROR",   "yellow",  "Face vertices are not coplanar"],
1969    ["DUP_FACE",    "ERROR",   "brown",   "Multiple instances of the same face."],
1970    ["MULTCONN",    "ERROR",   "orange",  "Multiply Connected Geometry. Too many faces attached at Edge"],
1971    ["REVERSAL",    "ERROR",   "violet",  "Faces Reverse Across Edge"],
1972    ["T_JUNCTION",  "ERROR",   "magenta", "Vertex is mid-edge on another Face"],
1973    ["FACE_ISECT",  "ERROR",   "brown",   "Faces intersect"],
1974    ["HOLE_EDGE",   "ERROR",   "red",     "Edge bounds Hole"]
1975];
1976
1977
1978function _vnf_validate_err(name, extra) =
1979    let(
1980        info = [for (x = _vnf_validate_errs) if (x[0] == name) x][0]
1981    ) concat(info, [extra]);
1982
1983
1984function _pts_not_reported(pts, varr, reports) =
1985    [
1986        for (i = pts, report = reports, pt = report[3])
1987        if (varr[i] == pt) 1
1988    ] == [];
1989
1990
1991function _edge_not_reported(edge, varr, reports) =
1992    let(
1993        edge = sort([for (i=edge) varr[i]])
1994    ) [
1995        for (report = reports) let(
1996            pts = sort(report[3])
1997        ) if (len(pts)==2 && edge == pts) 1
1998    ] == [];
1999
2000
2001module vnf_validate(vnf, size=1, show_warns=true, check_isects=false, opacity=0.67, adjacent=false, label_verts=false, label_faces=false, wireframe=false) {
2002    no_children($children);
2003    vcount = len(vnf[0]);
2004    fcount = len(vnf[1]);
2005    vnf = vnf_merge_points(vnf);
2006    faults = _vnf_validate(
2007        vnf, show_warns=show_warns,
2008        check_isects=check_isects
2009    );
2010    verts = vnf[0];
2011    vnf_changed = len(verts)!=vcount || len(vnf[1])!=fcount;
2012    if (!faults) {
2013        echo("VNF appears valid.");
2014    }
2015    if (vnf_changed) echo("VNF changed when merging points; unable to display indices");
2016    for (fault = faults) {
2017        err = fault[0];
2018        typ = fault[1];
2019        clr = fault[2];
2020        msg = fault[3];
2021        idxs = fault[4];
2022        pts = err=="FACE_ISECT" ? idxs : [for (i=idxs) if(is_finite(i) && i>=0 && i<len(verts)) verts[i]];
2023        if (vnf_changed || err=="FACE_ISECT")
2024          echo(str(typ, " ", err, " (", clr ,"): ", msg, " at ", pts));
2025        else
2026          echo(str(typ, " ", err, " (", clr ,"): ", msg, " at ", pts, " indices: ", idxs));
2027        color(clr) {
2028            if (is_vector(pts[0])) {
2029                if (len(pts)==2) {
2030                    stroke(pts, width=size, endcaps="butt", $fn=8);
2031                } else if (len(pts)>2) {
2032                    stroke(pts, width=size, closed=true, $fn=8);
2033                    polyhedron(pts,[[for (i=idx(pts)) i]]);
2034                } else {
2035                    move_copies(pts) sphere(d=size*3, $fn=18);
2036                }
2037            }
2038        }
2039    }
2040    badverts = unique([for (fault=faults) each fault[4]]);
2041    badverts2 = unique([for (j=idx(verts), i=badverts) if (i!=j && verts[i]==verts[j]) j]);
2042    all_badverts = unique(concat(badverts, badverts2));
2043    adjacent = !faults? false : adjacent;
2044    filter_fn = !adjacent? undef : function(i) in_list(i,all_badverts);
2045    adj_vnf = !adjacent? vnf : [
2046        verts, [for (face=vnf[1]) if (any(face,filter_fn)) face]
2047    ];
2048    if (wireframe) {
2049        vnf_wireframe(adj_vnf, width=size*0.25);
2050    }
2051    if (label_verts) {
2052        debug_vnf(adj_vnf, size=size*3, opacity=0, faces=false, vertices=true, filter=filter_fn);
2053    }
2054    if (label_faces) {
2055        debug_vnf(vnf, size=size*3, opacity=0, faces=true, vertices=false, filter=filter_fn);
2056    }
2057    if (opacity > 0) {
2058        color([0.5,1,0.5,opacity]) vnf_polyhedron(adj_vnf);
2059    }
2060}
2061
2062
2063
2064// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap