1//////////////////////////////////////////////////////////////////////
   2// LibFile: skin.scad
   3//   This file provides functions and modules that construct shapes from a list of cross sections.
   4//   In the case of skin() you specify each cross sectional shape yourself, and the number of
   5//   points can vary.  The various forms of sweep use a fixed shape, which may follow a path, or
   6//   be transformed in other ways to produce the list of cross sections.  In all cases it is the
   7//   user's responsibility to avoid creating a self-intersecting shape, which will produce
   8//   cryptic CGAL errors.  This file was inspired by list-comprehension-demos skin():
   9//   - https://github.com/openscad/list-comprehension-demos/blob/master/skin.scad
  10// Includes:
  11//   include <BOSL2/std.scad>
  12// FileGroup: Advanced Modeling
  13// FileSummary: Construct 3D shapes from 2D cross sections of the desired shape.
  14// FileFootnotes: STD=Included in std.scad
  15//////////////////////////////////////////////////////////////////////
  16
  17
  18// Section: Skin and sweep
  19
  20// Function&Module: skin()
  21// Synopsis: Connect a sequence of arbitrary polygons into a 3D object. 
  22// SynTags: VNF, Geom
  23// Topics: Extrusion, Skin
  24// See Also: vnf_vertex_array(), sweep(), linear_sweep(), rotate_sweep(), spiral_sweep(), path_sweep(), offset_sweep()
  25// Usage: As module:
  26//   skin(profiles, slices, [z=], [refine=], [method=], [sampling=], [caps=], [closed=], [style=], [convexity=], [anchor=],[cp=],[spin=],[orient=],[atype=]) [ATTACHMENTS];
  27// Usage: As function:
  28//   vnf = skin(profiles, slices, [z=], [refine=], [method=], [sampling=], [caps=], [closed=], [style=], [anchor=],[cp=],[spin=],[orient=],[atype=]);
  29// Description:
  30//   Given a list of two or more path `profiles` in 3d space, produces faces to skin a surface between
  31//   the profiles.  Optionally the first and last profiles can have endcaps, or the first and last profiles
  32//   can be connected together.  Each profile should be roughly planar, but some variation is allowed.
  33//   Each profile must rotate in the same clockwise direction.  If called as a function, returns a
  34//   [VNF structure](vnf.scad) `[VERTICES, FACES]`.  If called as a module, creates a polyhedron
  35//    of the skinned profiles.
  36//   .
  37//   The profiles can be specified either as a list of 3d curves or they can be specified as
  38//   2d curves with heights given in the `z` parameter.  It is your responsibility to ensure
  39//   that the resulting polyhedron is free from self-intersections, which would make it invalid
  40//   and can result in cryptic CGAL errors upon rendering with a second object present, even though the polyhedron appears
  41//   OK during preview or when rendered by itself.
  42//   .
  43//   For this operation to be well-defined, the profiles must all have the same vertex count and
  44//   we must assume that profiles are aligned so that vertex `i` links to vertex `i` on all polygons.
  45//   Many interesting cases do not comply with this restriction.  Two basic methods can handle
  46//   these cases: either subdivide edges (insert additional points along edges)
  47//   or duplicate vertcies (insert edges of length 0) so that both polygons have
  48//   the same number of points.
  49//   Duplicating vertices allows two distinct points in one polygon to connect to a single point
  50//   in the other one, creating
  51//   triangular faces.  You can adjust non-matching polygons yourself
  52//   either by resampling them using {{subdivide_path()}} or by duplicating vertices using
  53//   `repeat_entries`.  It is OK to pass a polygon that has the same vertex repeated, such as
  54//   a square with 5 points (two of which are identical), so that it can match up to a pentagon.
  55//   Such a combination would create a triangular face at the location of the duplicated vertex.
  56//   Alternatively, `skin` provides methods (described below) for inserting additional vertices
  57//   automatically to make incompatible paths match.
  58//   .
  59//   In order for skinned surfaces to look good it is usually necessary to use a fine sampling of
  60//   points on all of the profiles, and a large number of extra interpolated slices between the
  61//   profiles that you specify.  It is generally best if the triangles forming your polyhedron
  62//   are approximately equilateral.  The `slices` parameter specifies the number of slices to insert
  63//   between each pair of profiles, either a scalar to insert the same number everywhere, or a vector
  64//   to insert a different number between each pair.
  65//   .
  66//   Resampling may occur, depending on the `method` parameter, to make profiles compatible.
  67//   To force (possibly additional) resampling of the profiles to increase the point density you can set `refine=N`, which
  68//   will multiply the number of points on your profile by `N`.  You can choose between two resampling
  69//   schemes using the `sampling` option, which you can set to `"length"` or `"segment"`.
  70//   The length resampling method resamples proportional to length.
  71//   The segment method divides each segment of a profile into the same number of points.
  72//   This means that if you refine a profile with the "segment" method you will get N points
  73//   on each edge, but if you refine a profile with the "length" method you will get new points
  74//   distributed around the profile based on length, so small segments will get fewer new points than longer ones.
  75//   A uniform division may be impossible, in which case the code computes an approximation, which may result
  76//   in arbitrary distribution of extra points.  See {{subdivide_path()}} for more details.
  77//   Note that when dealing with continuous curves it is always better to adjust the
  78//   sampling in your code to generate the desired sampling rather than using the `refine` argument.
  79//   .
  80//   You can choose from five methods for specifying alignment for incommensurate profiles.
  81//   The available methods are `"distance"`, `"fast_distance"`, `"tangent"`, `"direct"` and `"reindex"`.
  82//   It is useful to distinguish between continuous curves like a circle and discrete profiles
  83//   like a hexagon or star, because the algorithms' suitability depend on this distinction.
  84//   .
  85//   The default method for aligning profiles is `method="direct"`.
  86//   If you simply supply a list of compatible profiles it will link them up
  87//   exactly as you have provided them.  You may find that profiles you want to connect define the
  88//   right shapes but the point lists don't start from points that you want aligned in your skinned
  89//   polyhedron.  You can correct this yourself using `reindex_polygon`, or you can use the "reindex"
  90//   method which will look for the index choice that will minimize the length of all of the edges
  91//   in the polyhedron&mdash;it will produce the least twisted possible result.  This algorithm has quadratic
  92//   run time so it can be slow with very large profiles.
  93//   .
  94//   When the profiles are incommensurate, the "direct" and "reindex" resample them to match.  As noted above,
  95//   for continuous input curves, it is better to generate your curves directly at the desired sample size,
  96//   but for mapping between a discrete profile like a hexagon and a circle, the hexagon must be resampled
  97//   to match the circle.  When you use "direct" or "reindex" the default `sampling` value is
  98//   of `sampling="length"` to approximate a uniform length sampling of the profile.  This will generally
  99//   produce the natural result for connecting two continuously sampled profiles or a continuous
 100//   profile and a polygonal one.  However depending on your particular case,
 101//   `sampling="segment"` may produce a more pleasing result.  These two approaches differ only when
 102//   the segments of your input profiles have unequal length.
 103//   .
 104//   The "distance", "fast_distance" and "tangent" methods work by duplicating vertices to create
 105//   triangular faces.  In the skined object created by two polygons, every vertex of a polygon must
 106//   have an edge that connects to some vertex on the other one.  If you connect two squares this can be
 107//   accomplished with four edges, but if you want to connect a square to a pentagon you must add a
 108//   fifth edge for the "extra" vertex on the pentagon.  You must now decide which vertex on the square to
 109//   connect the "extra" edge to.  How do you decide where to put that fifth edge?  The "distance" method answers this
 110//   question by using an optimization: it minimizes the total length of all the edges connecting
 111//   the two polygons.   This algorithm generally produces a good result when both profiles are discrete ones with
 112//   a small number of vertices.  It is computationally intensive (O(N^3)) and may be
 113//   slow on large inputs.  The resulting surfaces generally have curved faces, so be
 114//   sure to select a sufficiently large value for `slices` and `refine`.  Note that for
 115//   this method, `sampling` must be set to `"segment"`, and hence this is the default setting.
 116//   Using sampling by length would ignore the repeated vertices and ruin the alignment.
 117//   The "fast_distance" method restricts the optimization by assuming that an edge should connect
 118//   vertex 0 of the two polygons.  This reduces the run time to O(N^2) and makes
 119//   the method usable on profiles with more points if you take care to index the inputs to match.
 120//   .
 121//   The `"tangent"` method generally produces good results when
 122//   connecting a discrete polygon to a convex, finely sampled curve.  Given a polygon and a curve, consider one edge
 123//   on the polygon.  Find a plane passing through the edge that is tangent to the curve.  The endpoints of the edge and
 124//   the point of tangency define a triangular face in the output polyhedron.  If you work your way around the polygon
 125//   edges, you can establish a series of triangular faces in this way, with edges linking the polygon to the curve.
 126//   You can then complete the edge assignment by connecting all the edges in between the triangular faces together,
 127//   with many edges meeting at each polygon vertex.  The result is an alternation of flat triangular faces with conical
 128//   curves joining them.  Another way to think about it is that it splits the points on the curve up into groups and
 129//   connects all the points in one group to the same vertex on the polygon.
 130//   .
 131//   The "tangent" method may fail if the curved profile is non-convex, or doesn't have enough points to distinguish
 132//   all of the tangent points from each other.    The algorithm treats whichever input profile has fewer points as the polygon
 133//   and the other one as the curve.  Using `refine` with this method will have little effect on the model, so
 134//   you should do it only for agreement with other profiles, and these models are linear, so extra slices also
 135//   have no effect.  For best efficiency set `refine=1` and `slices=0`.  As with the "distance" method, refinement
 136//   must be done using the "segment" sampling scheme to preserve alignment across duplicated points.
 137//   Note that the "tangent" method produces similar results to the "distance" method on curved inputs.  If this
 138//   method fails due to concavity, "fast_distance" may be a good option.
 139//   .
 140//   It is possible to specify `method` and `refine` as arrays, but it is important to observe
 141//   matching rules when you do this.  If a pair of profiles is connected using "tangent" or "distance"
 142//   then the `refine` values for those two profiles must be equal.  If a profile is connected by
 143//   a vertex duplicating method on one side and a resampling method on the other side, then
 144//   `refine` must be set so that the resulting number of vertices matches the number that is
 145//   used for the resampled profiles.  The best way to avoid confusion is to ensure that the
 146//   profiles connected by "direct" or "reindex" all have the same number of points and at the
 147//   transition, the refined number of points matches.
 148//   .
 149// Arguments:
 150//   profiles = list of 2d or 3d profiles to be skinned.  (If 2d must also give `z`.)
 151//   slices = scalar or vector number of slices to insert between each pair of profiles.  Set to zero to use only the profiles you provided.  Recommend starting with a value around 10.
 152//   ---
 153//   refine = resample profiles to this number of points per edge.  Can be a list to give a refinement for each profile.  Recommend using a value above 10 when using the "distance" or "fast_distance" methods.  Default: 1.
 154//   sampling = sampling method to use with "direct" and "reindex" methods.  Can be "length" or "segment".  Ignored if any profile pair uses either the "distance", "fast_distance", or "tangent" methods.  Default: "length".
 155//   closed = set to true to connect first and last profile (to make a torus).  Default: false
 156//   caps = true to create endcap faces when closed is false.  Can be a length 2 boolean array.  Default is true if closed is false.
 157//   method = method for connecting profiles, one of "distance", "fast_distance", "tangent", "direct" or "reindex".  Default: "direct".
 158//   z = array of height values for each profile if the profiles are 2d
 159//   convexity = convexity setting for use with polyhedron.  (module only) Default: 10
 160//   anchor = Translate so anchor point is at the origin.  Default: "origin"
 161//   spin = Rotate this many degrees around Z axis after anchor.  Default: 0
 162//   orient = Vector to rotate top towards after spin
 163//   atype = Select "hull" or "intersect" anchor types. Default: "hull"
 164//   cp = Centerpoint for determining "intersect" anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
 165//   style = vnf_vertex_array style.  Default: "min_edge"
 166// Named Anchors:
 167//   "origin" = The native position of the shape.  
 168// Anchor Types:
 169//   "hull" = Anchors to the virtual convex hull of the shape.
 170//   "intersect" = Anchors to the surface of the shape.
 171// Example:
 172//   skin([octagon(4), circle($fn=70,r=2)], z=[0,3], slices=10);
 173// Example: Rotating the pentagon place the zero index at different locations, giving a twist
 174//   skin([rot(90,p=pentagon(4)), circle($fn=80,r=2)], z=[0,3], slices=10);
 175// Example: You can untwist it with the "reindex" method
 176//   skin([rot(90,p=pentagon(4)), circle($fn=80,r=2)], z=[0,3], slices=10, method="reindex");
 177// Example: Offsetting the starting edge connects to circles in an interesting way:
 178//   circ = circle($fn=80, r=3);
 179//   skin([circ, rot(110,p=circ)], z=[0,5], slices=20);
 180// Example(FlatSpin,VPD=20):
 181//   skin([ yrot(37,p=path3d(circle($fn=128, r=4))), path3d(square(3),3)], method="reindex",slices=10);
 182// Example(FlatSpin,VPD=16): Ellipses connected with twist
 183//   ellipse = xscale(2.5,p=circle($fn=80));
 184//   skin([ellipse, rot(45,p=ellipse)], z=[0,1.5], slices=10);
 185// Example(FlatSpin,VPD=16): Ellipses connected without a twist.  (Note ellipses stay in the same position: just the connecting edges are different.)
 186//   ellipse = xscale(2.5,p=circle($fn=80));
 187//   skin([ellipse, rot(45,p=ellipse)], z=[0,1.5], slices=10, method="reindex");
 188// Example(FlatSpin,VPD=500):
 189//   $fn=24;
 190//   skin([
 191//         yrot(0, p=yscale(2,p=path3d(circle(d=75)))),
 192//         [[40,0,100], [35,-15,100], [20,-30,100],[0,-40,100],[-40,0,100],[0,40,100],[20,30,100], [35,15,100]]
 193//   ],slices=10);
 194// Example(FlatSpin,VPD=600):
 195//   $fn=48;
 196//   skin([
 197//       for (b=[0,90]) [
 198//           for (a=[360:-360/$fn:0.01])
 199//               point3d(polar_to_xy((100+50*cos((a+b)*2))/2,a),b/90*100)
 200//       ]
 201//   ], slices=20);
 202// Example: Vaccum connector example from list-comprehension-demos
 203//   include <BOSL2/rounding.scad>
 204//   $fn=32;
 205//   base = round_corners(square([2,4],center=true), radius=0.5);
 206//   skin([
 207//       path3d(base,0),
 208//       path3d(base,2),
 209//       path3d(circle(r=0.5),3),
 210//       path3d(circle(r=0.5),4),
 211//       for(i=[0:2]) each [path3d(circle(r=0.6), i+4),
 212//                          path3d(circle(r=0.5), i+5)]
 213//   ],slices=0);
 214// Example: Vaccum nozzle example from list-comprehension-demos, using "length" sampling (the default)
 215//   xrot(90)down(1.5)
 216//   difference() {
 217//       skin(
 218//           [square([2,.2],center=true),
 219//            circle($fn=64,r=0.5)], z=[0,3],
 220//           slices=40,sampling="length",method="reindex");
 221//       skin(
 222//           [square([1.9,.1],center=true),
 223//            circle($fn=64,r=0.45)], z=[-.01,3.01],
 224//           slices=40,sampling="length",method="reindex");
 225//   }
 226// Example: Same thing with "segment" sampling
 227//   xrot(90)down(1.5)
 228//   difference() {
 229//       skin(
 230//           [square([2,.2],center=true),
 231//            circle($fn=64,r=0.5)], z=[0,3],
 232//           slices=40,sampling="segment",method="reindex");
 233//       skin(
 234//           [square([1.9,.1],center=true),
 235//            circle($fn=64,r=0.45)], z=[-.01,3.01],
 236//           slices=40,sampling="segment",method="reindex");
 237//   }
 238// Example: Forma Candle Holder (from list-comprehension-demos)
 239//   r = 50;
 240//   height = 140;
 241//   layers = 10;
 242//   wallthickness = 5;
 243//   holeradius = r - wallthickness;
 244//   difference() {
 245//       skin([for (i=[0:layers-1]) zrot(-30*i,p=path3d(hexagon(ir=r),i*height/layers))],slices=0);
 246//       up(height/layers) cylinder(r=holeradius, h=height);
 247//   }
 248// Example(FlatSpin,VPD=300): A box that is octagonal on the outside and circular on the inside
 249//   height = 45;
 250//   sub_base = octagon(d=71, rounding=2, $fn=128);
 251//   base = octagon(d=75, rounding=2, $fn=128);
 252//   interior = regular_ngon(n=len(base), d=60);
 253//   right_half()
 254//     skin([ sub_base, base, base, sub_base, interior], z=[0,2,height, height, 2], slices=0, refine=1, method="reindex");
 255// Example: Connecting a pentagon and circle with the "tangent" method produces large triangular faces and cone shaped corners.
 256//   skin([pentagon(4), circle($fn=80,r=2)], z=[0,3], slices=10, method="tangent");
 257// Example: rounding corners of a square.  Note that `$fn` makes the number of points constant, and avoiding the `rounding=0` case keeps everything simple.  In this case, the connections between profiles are linear, so there is no benefit to setting `slices` bigger than zero.
 258//   shapes = [for(i=[.01:.045:2])zrot(-i*180/2,cp=[-8,0,0],p=xrot(90,p=path3d(regular_ngon(n=4, side=4, rounding=i, $fn=64))))];
 259//   rotate(180) skin( shapes, slices=0);
 260// Example: Here's a simplified version of the above, with `i=0` included.  That first layer doesn't look good.
 261//   shapes = [for(i=[0:.2:1]) path3d(regular_ngon(n=4, side=4, rounding=i, $fn=32),i*5)];
 262//   skin(shapes, slices=0);
 263// Example: You can fix it by specifying "tangent" for the first method, but you still need "direct" for the rest.
 264//   shapes = [for(i=[0:.2:1]) path3d(regular_ngon(n=4, side=4, rounding=i, $fn=32),i*5)];
 265//   skin(shapes, slices=0, method=concat(["tangent"],repeat("direct",len(shapes)-2)));
 266// Example(FlatSpin,VPD=35): Connecting square to pentagon using "direct" method.
 267//   skin([regular_ngon(n=4, r=4), regular_ngon(n=5,r=5)], z=[0,4], refine=10, slices=10);
 268// Example(FlatSpin,VPD=35): Connecting square to shifted pentagon using "direct" method.
 269//   skin([regular_ngon(n=4, r=4), right(4,p=regular_ngon(n=5,r=5))], z=[0,4], refine=10, slices=10);
 270// Example(FlatSpin,VPD=185): In this example reindexing does not fix the orientation of the triangle because it happens in 3d within skin(), so we have to reverse the triangle manually
 271//   ellipse = yscale(3,circle(r=10, $fn=32));
 272//   tri = move([-50/3,-9],[[0,0], [50,0], [0,27]]);
 273//   skin([ellipse, reverse(tri)], z=[0,20], slices=20, method="reindex");
 274// Example(FlatSpin,VPD=185): You can get a nicer transition by rotating the polygons for better alignment.  You have to resample yourself before calling `align_polygon`. The orientation is fixed so we do not need to reverse.
 275//   ellipse = yscale(3,circle(r=10, $fn=32));
 276//   tri = move([-50/3,-9],
 277//              subdivide_path([[0,0], [50,0], [0,27]], 32));
 278//   aligned = align_polygon(ellipse,tri, [0:5:180]);
 279//   skin([ellipse, aligned], z=[0,20], slices=20);
 280// Example(FlatSpin,VPD=35): The "distance" method is a completely different approach.
 281//   skin([regular_ngon(n=4, r=4), regular_ngon(n=5,r=5)], z=[0,4], refine=10, slices=10, method="distance");
 282// Example(FlatSpin,VPD=35,VPT=[0,0,4]): Connecting pentagon to heptagon inserts two triangular faces on each side
 283//   small = path3d(circle(r=3, $fn=5));
 284//   big = up(2,p=yrot( 0,p=path3d(circle(r=3, $fn=7), 6)));
 285//   skin([small,big],method="distance", slices=10, refine=10);
 286// Example(FlatSpin,VPD=35,VPT=[0,0,4]): But just a slight rotation of the top profile moves the two triangles to one end
 287//   small = path3d(circle(r=3, $fn=5));
 288//   big = up(2,p=yrot(14,p=path3d(circle(r=3, $fn=7), 6)));
 289//   skin([small,big],method="distance", slices=10, refine=10);
 290// Example(FlatSpin,VPD=32,VPT=[1.2,4.3,2]): Another "distance" example:
 291//   off = [0,2];
 292//   shape = turtle(["right",45,"move", "left",45,"move", "left",45, "move", "jump", [.5+sqrt(2)/2,8]]);
 293//   rshape = rot(180,cp=centroid(shape)+off, p=shape);
 294//   skin([shape,rshape],z=[0,4], method="distance",slices=10,refine=15);
 295// Example(FlatSpin,VPD=32,VPT=[1.2,4.3,2]): Slightly shifting the profile changes the optimal linkage
 296//   off = [0,1];
 297//   shape = turtle(["right",45,"move", "left",45,"move", "left",45, "move", "jump", [.5+sqrt(2)/2,8]]);
 298//   rshape = rot(180,cp=centroid(shape)+off, p=shape);
 299//   skin([shape,rshape],z=[0,4], method="distance",slices=10,refine=15);
 300// Example(FlatSpin,VPD=444,VPT=[0,0,50]): This optimal solution doesn't look terrible:
 301//   prof1 = path3d([[-50,-50], [-50,50], [50,50], [25,25], [50,0], [25,-25], [50,-50]]);
 302//   prof2 = path3d(regular_ngon(n=7, r=50),100);
 303//   skin([prof1, prof2], method="distance", slices=10, refine=10);
 304// Example(FlatSpin,VPD=444,VPT=[0,0,50]): But this one looks better.  The "distance" method doesn't find it because it uses two more edges, so it clearly has a higher total edge distance.  We force it by doubling the first two vertices of one of the profiles.
 305//   prof1 = path3d([[-50,-50], [-50,50], [50,50], [25,25], [50,0], [25,-25], [50,-50]]);
 306//   prof2 = path3d(regular_ngon(n=7, r=50),100);
 307//   skin([repeat_entries(prof1,[2,2,1,1,1,1,1]),
 308//         prof2],
 309//        method="distance", slices=10, refine=10);
 310// Example(FlatSpin,VPD=80,VPT=[0,0,7]): The "distance" method will often produces results similar to the "tangent" method if you use it with a polygon and a curve, but the results can also look like this:
 311//   skin([path3d(circle($fn=128, r=10)), xrot(39, p=path3d(square([8,10]),10))],  method="distance", slices=0);
 312// Example(FlatSpin,VPD=80,VPT=[0,0,7]): Using the "tangent" method produces:
 313//   skin([path3d(circle($fn=128, r=10)), xrot(39, p=path3d(square([8,10]),10))],  method="tangent", slices=0);
 314// Example(FlatSpin,VPD=74): Torus using hexagons and pentagons, where `closed=true`
 315//   hex = right(7,p=path3d(hexagon(r=3)));
 316//   pent = right(7,p=path3d(pentagon(r=3)));
 317//   N=5;
 318//   skin(
 319//        [for(i=[0:2*N-1]) yrot(360*i/2/N, p=(i%2==0 ? hex : pent))],
 320//        refine=1,slices=0,method="distance",closed=true);
 321// Example: A smooth morph is achieved when you can calculate all the slices yourself.  Since you provide all the slices, set `slices=0`.
 322//   skin([for(n=[.1:.02:.5])
 323//            yrot(n*60-.5*60,p=path3d(supershape(step=360/128,m1=5,n1=n, n2=1.7),5-10*n))],
 324//        slices=0);
 325// Example: Another smooth supershape morph:
 326//   skin([for(alpha=[-.2:.05:1.5])
 327//            path3d(supershape(step=360/256,m1=7, n1=lerp(2,3,alpha),
 328//                              n2=lerp(8,4,alpha), n3=lerp(4,17,alpha)),alpha*5)],
 329//        slices=0);
 330// Example: Several polygons connected using "distance"
 331//   skin([regular_ngon(n=4, r=3),
 332//         regular_ngon(n=6, r=3),
 333//         regular_ngon(n=9, r=4),
 334//         rot(17,p=regular_ngon(n=6, r=3)),
 335//         rot(37,p=regular_ngon(n=4, r=3))],
 336//        z=[0,2,4,6,9], method="distance", slices=10, refine=10);
 337// Example(FlatSpin,VPD=935,VPT=[75,0,123]): Vertex count of the polygon changes at every profile
 338//   skin([
 339//       for (ang = [0:10:90])
 340//       rot([0,ang,0], cp=[200,0,0], p=path3d(circle(d=100,$fn=12-(ang/10))))
 341//   ],method="distance",slices=10,refine=10);
 342// Example: Möbius Strip.  This is a tricky model because when you work your way around to the connection, the direction of the profiles is flipped, so how can the proper geometry be created?  The trick is to duplicate the first profile and turn the caps off.  The model closes up and forms a valid polyhedron.
 343//   skin([
 344//     for (ang = [0:5:360])
 345//     rot([0,ang,0], cp=[100,0,0], p=rot(ang/2, p=path3d(square([1,30],center=true))))
 346//   ], caps=false, slices=0, refine=20);
 347// Example: This model of two scutoids packed together is based on https://www.thingiverse.com/thing:3024272 by mathgrrl
 348//   sidelen = 10;  // Side length of scutoid
 349//   height = 25;   // Height of scutoid
 350//   angle = -15;   // Angle (twists the entire form)
 351//   push = -5;     // Push (translates the base away from the top)
 352//   flare = 1;     // Flare (the two pieces will be different unless this is 1)
 353//   midpoint = .5; // Height of the extra vertex (as a fraction of total height); the two pieces will be different unless this is .5)
 354//   pushvec = rot(angle/2,p=push*RIGHT);  // Push direction is the average of the top and bottom mating edges
 355//   pent = path3d(apply(move(pushvec)*rot(angle),pentagon(side=sidelen,align_side=RIGHT,anchor="side0")));
 356//   hex = path3d(hexagon(side=flare*sidelen, align_side=RIGHT, anchor="side0"),height);
 357//   pentmate = path3d(pentagon(side=flare*sidelen,align_side=LEFT,anchor="side0"),height);
 358//             // Native index would require mapping first and last vertices together, which is not allowed, so shift
 359//   hexmate = list_rotate(
 360//                           path3d(apply(move(pushvec)*rot(angle),hexagon(side=sidelen,align_side=LEFT,anchor="side0"))),
 361//                           -1);
 362//   join_vertex = lerp(
 363//                       mean(select(hex,1,2)),     // midpoint of "extra" hex edge
 364//                       mean(select(hexmate,0,1)), // midpoint of "extra" hexmate edge
 365//                       midpoint);
 366//   augpent = repeat_entries(pent, [1,2,1,1,1]);         // Vertex 1 will split at the top forming a triangular face with the hexagon
 367//   augpent_mate = repeat_entries(pentmate,[2,1,1,1,1]); // For mating pentagon it is vertex 0 that splits
 368//              // Middle is the interpolation between top and bottom except for the join vertex, which is doubled because it splits
 369//   middle = list_set(lerp(augpent,hex,midpoint),[1,2],[join_vertex,join_vertex]);
 370//   middle_mate = list_set(lerp(hexmate,augpent_mate,midpoint), [0,1], [join_vertex,join_vertex]);
 371//   skin([augpent,middle,hex],  slices=10, refine=10, sampling="segment");
 372//   color("green")skin([augpent_mate,middle_mate,hexmate],  slices=10,refine=10, sampling="segment");
 373// Example: If you create a self-intersecting polyhedron the result is invalid.  In some cases self-intersection may be obvous.  Here is a more subtle example.
 374//   skin([
 375//          for (a = [0:30:180]) let(
 376//              pos  = [-60*sin(a),     0, a    ],
 377//              pos2 = [-60*sin(a+0.1), 0, a+0.1]
 378//          ) move(pos,
 379//              p=rot(from=UP, to=pos2-pos,
 380//                  p=path3d(circle(d=150))
 381//              )
 382//          )
 383//      ],refine=1,slices=0);
 384//      color("red") {
 385//          zrot(25) fwd(130) xrot(75) {
 386//              linear_extrude(height=0.1) {
 387//                  ydistribute(25) {
 388//                      text(text="BAD POLYHEDRONS!", size=20, halign="center", valign="center");
 389//                      text(text="CREASES MAKE", size=20, halign="center", valign="center");
 390//                  }
 391//              }
 392//          }
 393//          up(160) zrot(25) fwd(130) xrot(75) {
 394//              stroke(zrot(30, p=yscale(0.5, p=circle(d=120))),width=10,closed=true);
 395//          }
 396//      }
 397module skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=false, z, style="min_edge", convexity=10,
 398            anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull")
 399{
 400    vnf = skin(profiles, slices, refine, method, sampling, caps, closed, z, style=style);
 401    vnf_polyhedron(vnf,convexity=convexity,spin=spin,anchor=anchor,orient=orient,atype=atype,cp=cp)
 402        children();
 403}
 404
 405
 406function skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=false, z, style="min_edge",
 407              anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull") =
 408  assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"")
 409  assert(is_def(slices),"The slices argument must be specified.")
 410  assert(is_list(profiles) && len(profiles)>1, "Must provide at least two profiles")
 411  let(
 412       profiles = [for(p=profiles) if (is_region(p) && len(p)==1) p[0] else p]
 413  )
 414  let( bad = [for(i=idx(profiles)) if (!(is_path(profiles[i]) && len(profiles[i])>2)) i])
 415  assert(len(bad)==0, str("Profiles ",bad," are not a paths or have length less than 3"))
 416  let(
 417    profcount = len(profiles) - (closed?0:1),
 418    legal_methods = ["direct","reindex","distance","fast_distance","tangent"],
 419    caps = is_def(caps) ? caps :
 420           closed ? false : true,
 421    capsOK = is_bool(caps) || is_bool_list(caps,2),
 422    fullcaps = is_bool(caps) ? [caps,caps] : caps,
 423    refine = is_list(refine) ? refine : repeat(refine, len(profiles)),
 424    slices = is_list(slices) ? slices : repeat(slices, profcount),
 425    refineOK = [for(i=idx(refine)) if (refine[i]<=0 || !is_integer(refine[i])) i],
 426    slicesOK = [for(i=idx(slices)) if (!is_integer(slices[i]) || slices[i]<0) i],
 427    maxsize = max_length(profiles),
 428    methodok = is_list(method) || in_list(method, legal_methods),
 429    methodlistok = is_list(method) ? [for(i=idx(method)) if (!in_list(method[i], legal_methods)) i] : [],
 430    method = is_string(method) ? repeat(method, profcount) : method,
 431    // Define to be zero where a resampling method is used and 1 where a vertex duplicator is used
 432    RESAMPLING = 0,
 433    DUPLICATOR = 1,
 434    method_type = [for(m = method) m=="direct" || m=="reindex" ? 0 : 1],
 435    sampling = is_def(sampling) ? sampling :
 436               in_list(DUPLICATOR,method_type) ? "segment" : "length"
 437  )
 438  assert(len(refine)==len(profiles), "refine list is the wrong length")
 439  assert(len(slices)==profcount, str("slices list must have length ",profcount))
 440  assert(slicesOK==[],str("slices must be nonnegative integers"))
 441  assert(refineOK==[],str("refine must be postive integer"))
 442  assert(methodok,str("method must be one of ",legal_methods,". Got ",method))
 443  assert(methodlistok==[], str("method list contains invalid method at ",methodlistok))
 444  assert(len(method) == profcount,"Method list is the wrong length")
 445  assert(in_list(sampling,["length","segment"]), "sampling must be set to \"length\" or \"segment\"")
 446  assert(sampling=="segment" || (!in_list("distance",method) && !in_list("fast_distance",method) && !in_list("tangent",method)), "sampling is set to \"length\" which is only allowed with methods \"direct\" and \"reindex\"")
 447  assert(capsOK, "caps must be boolean or a list of two booleans")
 448  assert(!closed || !caps, "Cannot make closed shape with caps")
 449  let(
 450    profile_dim=list_shape(profiles,2),
 451    profiles_zcheck = (profile_dim != 2) || (profile_dim==2 && is_list(z) && len(z)==len(profiles)),
 452    profiles_ok = (profile_dim==2 && is_list(z) && len(z)==len(profiles)) || profile_dim==3
 453  )
 454  assert(profiles_zcheck, "z parameter is invalid or has the wrong length.")
 455  assert(profiles_ok,"Profiles must all be 3d or must all be 2d, with matching length z parameter.")
 456  assert(is_undef(z) || profile_dim==2, "Do not specify z with 3d profiles")
 457  assert(profile_dim==3 || len(z)==len(profiles),"Length of z does not match length of profiles.")
 458  let(
 459    // Adjoin Z coordinates to 2d profiles
 460    profiles = profile_dim==3 ? profiles :
 461               [for(i=idx(profiles)) path3d(profiles[i], z[i])],
 462    // True length (not counting repeated vertices) of profiles after refinement
 463    refined_len = [for(i=idx(profiles)) refine[i]*len(profiles[i])],
 464    // Define this to be 1 if a profile is used on either side by a resampling method, zero otherwise.
 465    profile_resampled = [for(i=idx(profiles))
 466      1-(
 467           i==0 ?  method_type[0] * (closed? last(method_type) : 1) :
 468           i==len(profiles)-1 ? last(method_type) * (closed ? select(method_type,-2) : 1) :
 469         method_type[i] * method_type[i-1])],
 470    parts = search(1,[1,for(i=[0:1:len(profile_resampled)-2]) profile_resampled[i]!=profile_resampled[i+1] ? 1 : 0],0),
 471    plen = [for(i=idx(parts)) (i== len(parts)-1? len(refined_len) : parts[i+1]) - parts[i]],
 472    max_list = [for(i=idx(parts)) each repeat(max(select(refined_len, parts[i], parts[i]+plen[i]-1)), plen[i])],
 473    transition_profiles = [for(i=[(closed?0:1):1:profcount-1]) if (select(method_type,i-1) != method_type[i]) i],
 474    badind = [for(tranprof=transition_profiles) if (refined_len[tranprof] != max_list[tranprof]) tranprof]
 475  )
 476  assert(badind==[],str("Profile length mismatch at method transition at indices ",badind," in skin()"))
 477  let(
 478    full_list =    // If there are no duplicators then use more efficient where the whole input is treated together
 479      !in_list(DUPLICATOR,method_type) ?
 480         let(
 481             resampled = [for(i=idx(profiles)) subdivide_path(profiles[i], max_list[i], method=sampling)],
 482             fixedprof = [for(i=idx(profiles))
 483                             i==0 || method[i-1]=="direct" ? resampled[i]
 484                                                         : reindex_polygon(resampled[i-1],resampled[i])],
 485             sliced = slice_profiles(fixedprof, slices, closed)
 486            )
 487            [!closed ? sliced : concat(sliced,[sliced[0]])]
 488      :  // There are duplicators, so use approach where each pair is treated separately
 489      [for(i=[0:profcount-1])
 490        let(
 491          pair =
 492            method[i]=="distance" ? _skin_distance_match(profiles[i],select(profiles,i+1)) :
 493            method[i]=="fast_distance" ? _skin_aligned_distance_match(profiles[i], select(profiles,i+1)) :
 494            method[i]=="tangent" ? _skin_tangent_match(profiles[i],select(profiles,i+1)) :
 495            /*method[i]=="reindex" || method[i]=="direct" ?*/
 496               let( p1 = subdivide_path(profiles[i],max_list[i], method=sampling),
 497                    p2 = subdivide_path(select(profiles,i+1),max_list[i], method=sampling)
 498               ) (method[i]=="direct" ? [p1,p2] : [p1, reindex_polygon(p1, p2)]),
 499            nsamples =  method_type[i]==RESAMPLING ? len(pair[0]) :
 500               assert(refine[i]==select(refine,i+1),str("Refine value mismatch at indices ",[i,(i+1)%len(refine)],
 501                                                        ".  Method ",method[i]," requires equal values"))
 502               refine[i] * len(pair[0])
 503          )
 504          subdivide_and_slice(pair,slices[i], nsamples, method=sampling)],
 505      vnf=vnf_join(
 506          [for(i=idx(full_list))
 507              vnf_vertex_array(full_list[i], cap1=i==0 && fullcaps[0], cap2=i==len(full_list)-1 && fullcaps[1],
 508                               col_wrap=true, style=style)])
 509  )
 510  reorient(anchor,spin,orient,vnf=vnf,p=vnf,extent=atype=="hull",cp=cp);
 511
 512
 513
 514// Function&Module: linear_sweep()
 515// Synopsis: Create a linear extrusion from a path, with optional texturing. 
 516// SynTags: VNF, Geom
 517// Topics: Extrusion, Textures, Sweep
 518// See Also: rotate_sweep(), sweep(), spiral_sweep(), path_sweep(), offset_sweep()
 519// Usage: As Module
 520//   linear_sweep(region, [height], [center=], [slices=], [twist=], [scale=], [style=], [caps=], [convexity=]) [ATTACHMENTS];
 521// Usage: With Texturing
 522//   linear_sweep(region, [height], [center=], texture=, [tex_size=]|[tex_reps=], [tex_depth=], [style=], [tex_samples=], ...) [ATTACHMENTS];
 523// Usage: As Function
 524//   vnf = linear_sweep(region, [height], [center=], [slices=], [twist=], [scale=], [style=], [caps=]);
 525//   vnf = linear_sweep(region, [height], [center=], texture=, [tex_size=]|[tex_reps=], [tex_depth=], [style=], [tex_samples=], ...);
 526// Description:
 527//   If called as a module, creates a polyhedron that is the linear extrusion of the given 2D region or polygon.
 528//   If called as a function, returns a VNF that can be used to generate a polyhedron of the linear extrusion
 529//   of the given 2D region or polygon.  The benefit of using this, over using `linear_extrude region(rgn)` is
 530//   that it supports `anchor`, `spin`, `orient` and attachments.  You can also make more refined
 531//   twisted extrusions by using `maxseg` to subsample flat faces.
 532//   .
 533//   Anchoring for linear_sweep is based on the anchors for the swept region rather than from the polyhedron that is created.  This can produce more
 534//   predictable anchors for LEFT, RIGHT, FWD and BACK in many cases, but the anchors may only
 535//   be aproximately correct for twisted objects, and corner anchors may point in unexpected directions in some cases.
 536//   If you need anchors directly computed from the surface you can pass the vnf from linear_sweep
 537//   to {{vnf_polyhedron()}}, which will compute anchors directly from the full VNF.  
 538// Arguments:
 539//   region = The 2D [Region](regions.scad) or polygon that is to be extruded.
 540//   h / height / l / length = The height to extrude the region.  Default: 1
 541//   center = If true, the created polyhedron will be vertically centered.  If false, it will be extruded upwards from the XY plane.  Default: `false`
 542//   ---
 543//   twist = The number of degrees to rotate the top of the shape, clockwise around the Z axis, relative to the bottom.  Default: 0
 544//   scale = The amount to scale the top of the shape, in the X and Y directions, relative to the size of the bottom.  Default: 1
 545//   shift = The amount to shift the top of the shape, in the X and Y directions, relative to the position of the bottom.  Default: [0,0]
 546//   slices = The number of slices to divide the shape into along the Z axis, to allow refinement of detail, especially when working with a twist.  Default: `twist/5`
 547//   maxseg = If given, then any long segments of the region will be subdivided to be shorter than this length.  This can refine twisting flat faces a lot.  Default: `undef` (no subsampling)
 548//   texture = A texture name string, or a rectangular array of scalar height values (0.0 to 1.0), or a VNF tile that defines the texture to apply to vertical surfaces.  See {{texture()}} for what named textures are supported.
 549//   tex_size = An optional 2D target size for the textures.  Actual texture sizes will be scaled somewhat to evenly fit the available surface. Default: `[5,5]`
 550//   tex_reps = If given instead of tex_size, a 2-vector giving the number of texture tile repetitions in the horizontal and vertical directions on the extrusion.
 551//   tex_inset = If numeric, lowers the texture into the surface by the specified proportion, e.g. 0.5 would lower it half way into the surface.  If `true`, insets by exactly its full depth.  Default: `false`
 552//   tex_rot = Rotate texture by specified angle, which must be a multiple of 90 degrees.  Default: 0
 553//   tex_depth = Specify texture depth; if negative, invert the texture.  Default: 1.
 554//   tex_samples = Minimum number of "bend points" to have in VNF texture tiles.  Default: 8
 555//   style = The style to use when triangulating the surface of the object.  Valid values are `"default"`, `"alt"`, or `"quincunx"`.
 556//   caps = If false do not create end caps.  Can be a boolean vector.  Default: true
 557//   convexity = Max number of surfaces any single ray could pass through.  Module use only.
 558//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determines the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: `"centroid"`
 559//   atype = Set to "hull" or "intersect" to select anchor type.  Default: "hull"
 560//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `"origin"`
 561//   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
 562//   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#subsection-orient).  Default: `UP`
 563// Anchor Types:
 564//   "hull" = Anchors to the virtual convex hull of the shape.
 565//   "intersect" = Anchors to the surface of the shape.
 566//   "bbox" = Anchors to the bounding box of the extruded shape.
 567// Named Anchors:
 568//   "origin" = Centers the extruded shape vertically only, but keeps the original path positions in the X and Y.  Oriented UP.
 569//   "original_base" = Keeps the original path positions in the X and Y, but at the bottom of the extrusion.  Oriented UP.
 570// Example: Extruding a Compound Region.
 571//   rgn1 = [for (d=[10:10:60]) circle(d=d,$fn=8)];
 572//   rgn2 = [square(30,center=false)];
 573//   rgn3 = [for (size=[10:10:20]) move([15,15],p=square(size=size, center=true))];
 574//   mrgn = union(rgn1,rgn2);
 575//   orgn = difference(mrgn,rgn3);
 576//   linear_sweep(orgn,height=20,convexity=16);
 577// Example: With Twist, Scale, Shift, Slices and Maxseg.
 578//   rgn1 = [for (d=[10:10:60]) circle(d=d,$fn=8)];
 579//   rgn2 = [square(30,center=false)];
 580//   rgn3 = [
 581//       for (size=[10:10:20])
 582//       apply(
 583//          move([15,15]),
 584//          square(size=size, center=true)
 585//       )
 586//   ];
 587//   mrgn = union(rgn1,rgn2);
 588//   orgn = difference(mrgn,rgn3);
 589//   linear_sweep(
 590//       orgn, height=50, maxseg=2, slices=40,
 591//       twist=90, scale=0.5, shift=[10,5],
 592//       convexity=16
 593//   );
 594// Example: Anchors on an Extruded Region
 595//   rgn1 = [for (d=[10:10:60]) circle(d=d,$fn=8)];
 596//   rgn2 = [square(30,center=false)];
 597//   rgn3 = [
 598//       for (size=[10:10:20])
 599//       apply(
 600//           move([15,15]),
 601//           rect(size=size)
 602//       )
 603//   ];
 604//   mrgn = union(rgn1,rgn2);
 605//   orgn = difference(mrgn,rgn3);
 606//   linear_sweep(orgn,height=20,convexity=16)
 607//       show_anchors();
 608// Example: "diamonds" texture.
 609//   path = glued_circles(r=15, spread=40, tangent=45);
 610//   linear_sweep(
 611//       path, texture="diamonds", tex_size=[5,10],
 612//       h=40, style="concave");
 613// Example: "pyramids" texture.
 614//   linear_sweep(
 615//       rect(50), texture="pyramids", tex_size=[10,10],
 616//       h=40, style="convex");
 617// Example: "bricks_vnf" texture.
 618//   path = glued_circles(r=15, spread=40, tangent=45);
 619//   linear_sweep(
 620//       path, texture="bricks_vnf", tex_size=[10,10],
 621//       tex_depth=0.25, h=40);
 622// Example: User defined heightfield texture.
 623//   path = ellipse(r=[20,10]);
 624//   texture = [for (i=[0:9])
 625//       [for (j=[0:9])
 626//           1/max(0.5,norm([i,j]-[5,5])) ]];
 627//   linear_sweep(
 628//       path, texture=texture, tex_size=[5,5],
 629//       h=40, style="min_edge", anchor=BOT);
 630// Example: User defined VNF tile texture.
 631//   path = ellipse(r=[20,10]);
 632//   tex = let(n=16,m=0.25) [
 633//        [
 634//            each resample_path(path3d(square(1)),n),
 635//            each move([0.5,0.5],
 636//                p=path3d(circle(d=0.5,$fn=n),m)),
 637//            [1/2,1/2,0],
 638//        ], [
 639//            for (i=[0:1:n-1]) each [
 640//                [i,(i+1)%n,(i+3)%n+n],
 641//                [i,(i+3)%n+n,(i+2)%n+n],
 642//                [2*n,n+i,n+(i+1)%n],
 643//            ]
 644//        ]
 645//   ];
 646//   linear_sweep(path, texture=tex, tex_size=[5,5], h=40);
 647// Example: Textured with twist and scale.
 648//   linear_sweep(regular_ngon(n=3, d=50),
 649//       texture="rough", h=100, tex_depth=2,
 650//       tex_size=[20,20], style="min_edge",
 651//       convexity=10, scale=0.2, twist=120);
 652// Example: As Function
 653//   path = glued_circles(r=15, spread=40, tangent=45);
 654//   vnf = linear_sweep(
 655//       path, h=40, texture="trunc_pyramids", tex_size=[5,5],
 656//       tex_depth=1, style="convex");
 657//   vnf_polyhedron(vnf, convexity=10);
 658// Example: VNF tile that has no top/bottom edges and produces a disconnected result
 659//   shape = skin([rect(2/5),
 660//                 rect(2/3),
 661//                 rect(2/5)],
 662//                z=[0,1/2,1],
 663//                slices=0,
 664//                caps=false);
 665//   tile = move([0,1/2,2/3],yrot(90,shape));
 666//   linear_sweep(circle(20), texture=tile,
 667//                tex_size=[10,10],tex_depth=5,
 668//                h=40,convexity=4);
 669// Example: The same tile from above, turned 90 degrees, creates problems at the ends, because the end cap is not a connected polygon.  When the ends are disconnected you may find that some parts of the end cap are missing and spurious polygons included.  
 670//  shape = skin([rect(2/5),
 671//                rect(2/3),
 672//                rect(2/5)],
 673//               z=[0,1/2,1],
 674//               slices=0,
 675//               caps=false);
 676//  tile = move([1/2,1,2/3],xrot(90,shape));
 677//  linear_sweep(circle(20), texture=tile,
 678//               tex_size=[30,20],tex_depth=15,
 679//               h=40,convexity=4);
 680// Example: This example shows some endcap polygons missing and a spurious triangle
 681//   shape = skin([rect(2/5),
 682//                 rect(2/3),
 683//                 rect(2/5)],
 684//                z=[0,1/2,1],
 685//                slices=0,
 686//                caps=false);
 687//   tile = xscale(.5,move([1/2,1,2/3],xrot(90,shape)));
 688//   doubletile = vnf_join([tile, right(.5,tile)]);
 689//   linear_sweep(circle(20), texture=doubletile,
 690//                tex_size=[45,45],tex_depth=15, h=40);
 691// Example: You can fix ends for disconnected cases using {{top_half()}} and {{bottom_half()}}
 692//   shape = skin([rect(2/5),
 693//                 rect(2/3),
 694//                 rect(2/5)],
 695//                z=[0,1/2,1],
 696//                slices=0,
 697//                caps=false);
 698//   tile = move([1/2,1,2/3],xrot(90,shape));
 699//   vnf_polyhedron(
 700//     top_half(
 701//       bottom_half(
 702//         linear_sweep(circle(20), texture=tile,
 703//                     tex_size=[30,20],tex_depth=15,
 704//                     h=40.2,caps=false),
 705//       z=20),
 706//     z=-20)); 
 707
 708module linear_sweep(
 709    region, height, center,
 710    twist=0, scale=1, shift=[0,0],
 711    slices, maxseg, style="default", convexity, caps=true, 
 712    texture, tex_size=[5,5], tex_reps, tex_counts,
 713    tex_inset=false, tex_rot=0,
 714    tex_depth, tex_scale, tex_samples,
 715    cp, atype="hull", h,l,length,
 716    anchor, spin=0, orient=UP
 717) {
 718    h = one_defined([h, height,l,length],"h,height,l,length",dflt=1);
 719    region = force_region(region);
 720    check = assert(is_region(region),"Input is not a region");
 721    anchor = center==true? "origin" :
 722        center == false? "original_base" :
 723        default(anchor, "original_base");
 724    vnf = linear_sweep(
 725        region, height=h, style=style, caps=caps, 
 726        twist=twist, scale=scale, shift=shift,
 727        texture=texture,
 728        tex_size=tex_size,
 729        tex_reps=tex_reps,
 730        tex_counts=tex_counts,
 731        tex_inset=tex_inset,
 732        tex_rot=tex_rot,
 733        tex_depth=tex_depth,
 734        tex_samples=tex_samples,
 735        slices=slices,
 736        maxseg=maxseg,
 737        anchor="origin"
 738    );
 739    anchors = [
 740        named_anchor("original_base", [0,0,-h/2], UP)
 741    ];
 742    cp = default(cp, "centroid");
 743    geom = atype=="hull"?  attach_geom(cp=cp, region=region, h=h, extent=true, shift=shift, scale=scale, twist=twist, anchors=anchors) :
 744        atype=="intersect"?  attach_geom(cp=cp, region=region, h=h, extent=false, shift=shift, scale=scale, twist=twist, anchors=anchors) :
 745        atype=="bbox"?
 746            let(
 747                bounds = pointlist_bounds(flatten(region)),
 748                size = bounds[1] - bounds[0],
 749                midpt = (bounds[0] + bounds[1])/2
 750            )
 751            attach_geom(cp=[0,0,0], size=point3d(size,h), offset=point3d(midpt), shift=shift, scale=scale, twist=twist, anchors=anchors) :
 752        assert(in_list(atype, ["hull","intersect","bbox"]), "Anchor type must be \"hull\", \"intersect\", or \"bbox\".");
 753    attachable(anchor,spin,orient, geom=geom) {
 754        vnf_polyhedron(vnf, convexity=convexity);
 755        children();
 756    }
 757}
 758
 759
 760function linear_sweep(
 761    region, height, center,
 762    twist=0, scale=1, shift=[0,0],
 763    slices, maxseg, style="default", caps=true, 
 764    cp, atype="hull", h,
 765    texture, tex_size=[5,5], tex_reps, tex_counts,
 766    tex_inset=false, tex_rot=0,
 767    tex_scale, tex_depth, tex_samples, h, l, length, 
 768    anchor, spin=0, orient=UP
 769) =
 770    assert(num_defined([tex_reps,tex_counts])<2, "In linear_sweep() the 'tex_counts' parameter has been replaced by 'tex_reps'.  You cannot give both.")
 771    assert(num_defined([tex_scale,tex_depth])<2, "In linear_sweep() the 'tex_scale' parameter has been replaced by 'tex_depth'.  You cannot give both.")
 772    let(
 773        region = force_region(region),
 774        tex_reps = is_def(tex_counts)? echo("In linear_sweep() the 'tex_counts' parameter is deprecated and has been replaced by 'tex_reps'")tex_counts
 775                 : tex_reps,
 776        tex_depth = is_def(tex_scale)? echo("In linear_sweep() the 'tex_scale' parameter is deprecated and has been replaced by 'tex_depth'")tex_scale
 777                  : default(tex_depth,1)
 778    )
 779    assert(is_region(region), "Input is not a region or polygon.")
 780    assert(is_num(scale) || is_vector(scale))
 781    assert(is_vector(shift, 2), str(shift))
 782    assert(is_bool(caps) || is_bool_list(caps,2), "caps must be boolean or a list of two booleans")
 783    let(
 784        h = one_defined([h, height,l,length],"h,height,l,length",dflt=1)
 785    )
 786    !is_undef(texture)? _textured_linear_sweep(
 787        region, h=h, caps=caps, 
 788        texture=texture, tex_size=tex_size,
 789        counts=tex_reps, inset=tex_inset,
 790        rot=tex_rot, tex_scale=tex_depth,
 791        twist=twist, scale=scale, shift=shift,
 792        style=style, samples=tex_samples,
 793        anchor=anchor, spin=spin, orient=orient
 794    ) :
 795    let(
 796        caps = is_bool(caps) ? [caps,caps] : caps, 
 797        anchor = center==true? "origin" :
 798            center == false? "original_base" :
 799            default(anchor, "original_base"),
 800        regions = region_parts(region),
 801        slices = default(slices, max(1,ceil(abs(twist)/5))),
 802        scale = is_num(scale)? [scale,scale] : point2d(scale),
 803        topmat = move(shift) * scale(scale) * rot(-twist),
 804        trgns = [
 805            for (rgn = regions) [
 806                for (path = rgn) let(
 807                    p = list_unwrap(path),
 808                    path = is_undef(maxseg)? p : [
 809                        for (seg = pair(p,true)) each
 810                        let( steps = ceil(norm(seg.y - seg.x) / maxseg) )
 811                        lerpn(seg.x, seg.y, steps, false)
 812                    ]
 813                ) apply(topmat, path)
 814            ]
 815        ],
 816        vnf = vnf_join([
 817            for (rgn = regions)
 818            for (pathnum = idx(rgn)) let(
 819                p = list_unwrap(rgn[pathnum]),
 820                path = is_undef(maxseg)? p : [
 821                    for (seg=pair(p,true)) each
 822                    let(steps=ceil(norm(seg.y-seg.x)/maxseg))
 823                    lerpn(seg.x, seg.y, steps, false)
 824                ],
 825                verts = [
 826                    for (i=[0:1:slices]) let(
 827                        u = i / slices,
 828                        scl = lerp([1,1], scale, u),
 829                        ang = lerp(0, -twist, u),
 830                        off = lerp([0,0,-h/2], point3d(shift,h/2), u),
 831                        m = move(off) * scale(scl) * rot(ang)
 832                    ) apply(m, path3d(path))
 833                ]
 834            ) vnf_vertex_array(verts, caps=false, col_wrap=true, style=style),
 835            if (caps[0]) for (rgn = regions) vnf_from_region(rgn, down(h/2), reverse=true),
 836            if (caps[1]) for (rgn = trgns) vnf_from_region(rgn, up(h/2), reverse=false)
 837        ]),
 838        anchors = [
 839            named_anchor("original_base", [0,0,-h/2], UP)
 840        ],
 841        cp = default(cp, "centroid"),
 842        geom = atype=="hull"?  attach_geom(cp=cp, region=region, h=h, extent=true, shift=shift, scale=scale, twist=twist, anchors=anchors) :
 843            atype=="intersect"?  attach_geom(cp=cp, region=region, h=h, extent=false, shift=shift, scale=scale, twist=twist, anchors=anchors) :
 844            atype=="bbox"?
 845                let(
 846                    bounds = pointlist_bounds(flatten(region)),
 847                    size = bounds[1] - bounds[0],
 848                    midpt = (bounds[0] + bounds[1])/2
 849                )
 850                attach_geom(cp=[0,0,0], size=point3d(size,h), offset=point3d(midpt), shift=shift, scale=scale, twist=twist, anchors=anchors) :
 851            assert(in_list(atype, ["hull","intersect","bbox"]), "Anchor type must be \"hull\", \"intersect\", or \"bbox\".")
 852    ) reorient(anchor,spin,orient, geom=geom, p=vnf);
 853
 854
 855// Function&Module: rotate_sweep()
 856// Synopsis: Create a surface of revolution from a path with optional texturing. 
 857// SynTags: VNF, Geom
 858// Topics: Extrusion, Sweep, Revolution, Textures
 859// See Also: linear_sweep(), sweep(), spiral_sweep(), path_sweep(), offset_sweep()
 860// Usage: As Function
 861//   vnf = rotate_sweep(shape, [angle], ...);
 862// Usage: As Module
 863//   rotate_sweep(shape, [angle], ...) [ATTACHMENTS];
 864// Usage: With Texturing
 865//   rotate_sweep(shape, texture=, [tex_size=]|[tex_reps=], [tex_depth=], [tex_samples=], [tex_rot=], [tex_inset=], ...) [ATTACHMENTS];
 866// Description:
 867//   Takes a polygon or [region](regions.scad) and sweeps it in a rotation around the Z axis, with optional texturing.
 868//   When called as a function, returns a [VNF](vnf.scad).
 869//   When called as a module, creates the sweep as geometry.
 870// Arguments:
 871//   shape = The polygon or [region](regions.scad) to sweep around the Z axis.
 872//   angle = If given, specifies the number of degrees to sweep the shape around the Z axis, counterclockwise from the X+ axis.  Default: 360 (full rotation)
 873//   ---
 874//   texture = A texture name string, or a rectangular array of scalar height values (0.0 to 1.0), or a VNF tile that defines the texture to apply to vertical surfaces.  See {{texture()}} for what named textures are supported.
 875//   tex_size = An optional 2D target size for the textures.  Actual texture sizes will be scaled somewhat to evenly fit the available surface. Default: `[5,5]`
 876//   tex_reps = If given instead of tex_size, a 2-vector giving the number of texture tile repetitions in the direction perpendicular to extrusion and in the direction parallel to extrusion.  
 877//   tex_inset = If numeric, lowers the texture into the surface by the specified proportion, e.g. 0.5 would lower it half way into the surface.  If `true`, insets by exactly its full depth.  Default: `false`
 878//   tex_rot = Rotate texture by specified angle, which must be a multiple of 90 degrees.  Default: 0
 879//   tex_depth = Specify texture depth; if negative, invert the texture.  Default: 1.
 880//   tex_samples = Minimum number of "bend points" to have in VNF texture tiles.  Default: 8
 881//   style = {{vnf_vertex_array()}} style.  Default: "min_edge"
 882//   closed = If false, and shape is given as a path, then the revolved path will be sealed to the axis of rotation with untextured caps.  Default: `true`
 883//   convexity = (Module only) Convexity setting for use with polyhedron.  Default: 10
 884//   cp = Centerpoint for determining "intersect" anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
 885//   atype = Select "hull" or "intersect" anchor types.  Default: "hull"
 886//   anchor = Translate so anchor point is at the origin. Default: "origin"
 887//   spin = Rotate this many degrees around Z axis after anchor. Default: 0
 888//   orient = Vector to rotate top towards after spin  (module only)
 889// Named Anchors:
 890//   "origin" = The native position of the shape.  
 891// Anchor Types:
 892//   "hull" = Anchors to the virtual convex hull of the shape.
 893//   "intersect" = Anchors to the surface of the shape.
 894// Example:
 895//   rgn = [
 896//       for (a = [0, 120, 240]) let(
 897//           cp = polar_to_xy(15, a) + [30,0]
 898//       ) each [
 899//           move(cp, p=circle(r=10)),
 900//           move(cp, p=hexagon(d=15)),
 901//       ]
 902//   ];
 903//   rotate_sweep(rgn, angle=240);
 904// Example:
 905//   rgn = right(30, p=union([for (a = [0, 90]) rot(a, p=rect([15,5]))]));
 906//   rotate_sweep(rgn);
 907// Example:
 908//   path = right(50, p=circle(d=40));
 909//   rotate_sweep(path, texture="bricks_vnf", tex_size=[10,10], tex_depth=0.5, style="concave");
 910// Example:
 911//   tex = [
 912//       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 913//       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 914//       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
 915//       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
 916//       [0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1],
 917//       [0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1],
 918//       [0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1],
 919//       [0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1],
 920//       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
 921//       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
 922//       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 923//       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 924//   ];
 925//   path = arc(cp=[0,0], r=40, start=60, angle=-120);
 926//   rotate_sweep(
 927//       path, closed=false,
 928//       texture=tex, tex_size=[20,20],
 929//       tex_depth=1, style="concave");
 930// Example:
 931//   include <BOSL2/beziers.scad>
 932//   bezpath = [
 933//       [15, 30], [10,15],
 934//       [10,  0], [20, 10], [30,12],
 935//       [30,-12], [20,-10], [10, 0],
 936//       [10,-15], [15,-30]
 937//   ];
 938//   path = bezpath_curve(bezpath, splinesteps=32);
 939//   rotate_sweep(
 940//       path, closed=false,
 941//       texture="diamonds", tex_size=[10,10],
 942//       tex_depth=1, style="concave");
 943// Example:
 944//   path = [
 945//       [20, 30], [20, 20],
 946//       each arc(r=20, corner=[[20,20],[10,0],[20,-20]]),
 947//       [20,-20], [20,-30],
 948//   ];
 949//   vnf = rotate_sweep(
 950//       path, closed=false,
 951//       texture="trunc_pyramids",
 952//       tex_size=[5,5], tex_depth=1,
 953//       style="convex");
 954//   vnf_polyhedron(vnf, convexity=10);
 955// Example:
 956//   rgn = [
 957//       right(40, p=circle(d=50)),
 958//       right(40, p=circle(d=40,$fn=6)),
 959//   ];
 960//   rotate_sweep(
 961//       rgn, texture="diamonds",
 962//       tex_size=[10,10], tex_depth=1,
 963//       angle=240, style="concave");
 964
 965function rotate_sweep(
 966    shape, angle=360,
 967    texture, tex_size=[5,5], tex_counts, tex_reps, 
 968    tex_inset=false, tex_rot=0,
 969    tex_scale, tex_depth, tex_samples,
 970    tex_taper, shift=[0,0], closed=true,
 971    style="min_edge", cp="centroid",
 972    atype="hull", anchor="origin",
 973    spin=0, orient=UP,
 974    _tex_inhibit_y_slicing=false
 975) =
 976    assert(num_defined([tex_reps,tex_counts])<2, "In rotate_sweep() the 'tex_counts' parameters has been replaced by 'tex_reps'.  You cannot give both.")
 977    assert(num_defined([tex_scale,tex_depth])<2, "In linear_sweep() the 'tex_scale' parameter has been replaced by 'tex_depth'.  You cannot give both.")
 978    let( region = force_region(shape),
 979         tex_reps = is_def(tex_counts)? echo("In rotate_sweep() the 'tex_counts' parameter is deprecated and has been replaced by 'tex_reps'")tex_counts
 980                  : tex_reps,
 981        tex_depth = is_def(tex_scale)? echo("In rotate_sweep() the 'tex_scale' parameter is deprecated and has been replaced by 'tex_depth'")tex_scale
 982                  : default(tex_depth,1)
 983    )
 984    assert(is_region(region), "Input is not a region or polygon.")
 985    let(
 986        bounds = pointlist_bounds(flatten(region)),
 987        min_x = bounds[0].x,
 988        max_x = bounds[1].x,
 989        min_y = bounds[0].y,
 990        max_y = bounds[1].y,
 991        h = max_y - min_y
 992    )
 993    assert(min_x>=0, "Input region must exist entirely in the X+ half-plane.")
 994    !is_undef(texture)? _textured_revolution(
 995        shape,
 996        texture=texture,
 997        tex_size=tex_size,
 998        counts=tex_reps,
 999        tex_scale=tex_depth,
1000        inset=tex_inset,
1001        rot=tex_rot,
1002        samples=tex_samples,
1003        inhibit_y_slicing=_tex_inhibit_y_slicing,
1004        taper=tex_taper,
1005        shift=shift,
1006        closed=closed,
1007        angle=angle,
1008        style=style
1009    ) :
1010    let(
1011        steps = ceil(segs(max_x) * angle / 360) + (angle<360? 1 : 0),
1012        skmat = down(min_y) * skew(sxz=shift.x/h, syz=shift.y/h) * up(min_y),
1013        transforms = [
1014            if (angle==360) for (i=[0:1:steps-1]) skmat * rot([90,0,360-i*360/steps]),
1015            if (angle<360) for (i=[0:1:steps-1]) skmat * rot([90,0,angle-i*angle/(steps-1)]),
1016        ],
1017        vnf = sweep(
1018            region, transforms,
1019            closed=angle==360,
1020            caps=angle!=360,
1021            style=style, cp=cp,
1022            atype=atype, anchor=anchor,
1023            spin=spin, orient=orient
1024        )
1025    ) vnf;
1026
1027
1028module rotate_sweep(
1029    shape, angle=360,
1030    texture, tex_size=[5,5], tex_counts, tex_reps,
1031    tex_inset=false, tex_rot=0,
1032    tex_scale, tex_depth, tex_samples,
1033    tex_taper, shift=[0,0],
1034    style="min_edge",
1035    closed=true,
1036    cp="centroid",
1037    convexity=10,
1038    atype="hull",
1039    anchor="origin",
1040    spin=0,
1041    orient=UP,
1042    _tex_inhibit_y_slicing=false
1043) {
1044    dummy =
1045       assert(num_defined([tex_reps,tex_counts])<2, "In rotate_sweep() the 'tex_counts' parameters has been replaced by 'tex_reps'.  You cannot give both.")
1046       assert(num_defined([tex_scale,tex_depth])<2, "In rotate_sweep() the 'tex_scale' parameter has been replaced by 'tex_depth'.  You cannot give both.");
1047    tex_reps = is_def(tex_counts)? echo("In rotate_sweep() the 'tex_counts' parameter is deprecated and has been replaced by 'tex_reps'")tex_counts
1048             : tex_reps;
1049    tex_depth = is_def(tex_scale)? echo("In rotate_sweep() the 'tex_scale' parameter is deprecated and has been replaced by 'tex_depth'")tex_scale
1050              : default(tex_depth,1);
1051    region = force_region(shape);
1052    check = assert(is_region(region), "Input is not a region or polygon.");
1053    bounds = pointlist_bounds(flatten(region));
1054    min_x = bounds[0].x;
1055    max_x = bounds[1].x;
1056    min_y = bounds[0].y;
1057    max_y = bounds[1].y;
1058    h = max_y - min_y;
1059    check2 = assert(min_x>=0, "Input region must exist entirely in the X+ half-plane.");
1060    if (!is_undef(texture)) {
1061        _textured_revolution(
1062            shape,
1063            texture=texture,
1064            tex_size=tex_size,
1065            counts=tex_reps,
1066            tex_scale=tex_depth,
1067            inset=tex_inset,
1068            rot=tex_rot,
1069            samples=tex_samples,
1070            taper=tex_taper,
1071            shift=shift,
1072            closed=closed,
1073            inhibit_y_slicing=_tex_inhibit_y_slicing,
1074            angle=angle,
1075            style=style,
1076            atype=atype, anchor=anchor,
1077            spin=spin, orient=orient
1078        ) children();
1079    } else {
1080        steps = ceil(segs(max_x) * angle / 360) + (angle<360? 1 : 0);
1081        skmat = down(min_y) * skew(sxz=shift.x/h, syz=shift.y/h) * up(min_y);
1082        transforms = [
1083            if (angle==360) for (i=[0:1:steps-1]) skmat * rot([90,0,360-i*360/steps]),
1084            if (angle<360) for (i=[0:1:steps-1]) skmat * rot([90,0,angle-i*angle/(steps-1)]),
1085        ];
1086        sweep(
1087            region, transforms,
1088            closed=angle==360,
1089            caps=angle!=360,
1090            style=style, cp=cp,
1091            convexity=convexity,
1092            atype=atype, anchor=anchor,
1093            spin=spin, orient=orient
1094        ) children();
1095    }
1096}
1097
1098
1099// Function&Module: spiral_sweep()
1100// Synopsis: Sweep a path along a helix.
1101// SynTags: VNF, Geom
1102// Topics: Extrusion, Sweep, Spiral
1103// See Also: thread_helix(), linear_sweep(), rotate_sweep(), sweep(), path_sweep(), offset_sweep()
1104// Usage: As Module
1105//   spiral_sweep(poly, h, r|d=, turns, [taper=], [center=], [taper1=], [taper2=], [internal=], ...)[ATTACHMENTS];
1106//   spiral_sweep(poly, h, r1=|d1=, r2=|d2=, turns, [taper=], [center=], [taper1=], [taper2=], [internal=], ...)[ATTACHMENTS];
1107// Usage: As Function
1108//   vnf = spiral_sweep(poly, h, r|d=, turns, ...);
1109//   vnf = spiral_sweep(poly, h, r1=|d1=, r1=|d2=, turns, ...);
1110// Description:
1111//   Takes a closed 2D polygon path, centered on the XY plane, and sweeps/extrudes it along a 3D spiral path
1112//   of a given radius, height and degrees of rotation.  The origin in the profile traces out the helix of the specified radius.
1113//   If turns is positive the path will be right-handed;  if turns is negative the path will be left-handed.
1114//   Such an extrusion can be used to make screw threads.  
1115//   .
1116//   The lead_in options specify a lead-in section where the ends of the spiral scale down to avoid a sharp cut face at the ends.
1117//   You can specify the length of this scaling directly with the lead_in parameters or as an angle using the lead_in_ang parameters.
1118//   If you give a positive value, the extrusion is lengthenend by the specified distance or angle; if you give a negative
1119//   value then the scaled end is included in the extrusion length specified by `turns`.  If the value is zero then no scaled ends
1120//   are produced.  The shape of the scaled ends can be controlled with the lead_in_shape parameter.  Supported options are "sqrt", "linear"
1121//   "smooth" and "cut".  
1122//   .
1123//   The inside argument changes how the extrusion lead-in sections are formed.  If it is true then they scale
1124//   towards the outside, like would be needed for internal threading.  If internal is fale then the lead-in sections scale
1125//   towards the inside, like would be appropriate for external threads.  
1126// Arguments:
1127//   poly = Array of points of a polygon path, to be extruded.
1128//   h = height of the spiral extrusion path
1129//   r = Radius of the spiral extrusion path
1130//   turns = number of revolutions to include in the spiral
1131//   ---
1132//   d = Diameter of the spiral extrusion path.
1133//   d1/r1 = Bottom inside diameter or radius of spiral to extrude along.
1134//   d2/r2 = Top inside diameter or radius of spiral to extrude along.
1135//   lead_in = Specify linear length of the lead-in scaled section of the spiral.  Default: 0
1136//   lead_in1 = Specify linear length of the lead-in scaled section of the spiral at the bottom
1137//   lead_in2 = Specify linear length of the lead-in scaled section of the spiral at the top
1138//   lead_in_ang = Specify angular  length of the lead-in scaled section of the spiral
1139//   lead_in_ang1 = Specify angular length of the lead-in scaled section of the spiral at the bottom
1140//   lead_in_ang2 = Specify angular length of the lead-in scaled section of the spiral at the top
1141//   lead_in_shape = Specify the shape of the thread lead in by giving a text string or function.  Default: "sqrt"
1142//   lead_in_shape1 = Specify the shape of the thread lead-in at the bottom by giving a text string or function.  
1143//   lead_in_shape2 = Specify the shape of the thread lead-in at the top by giving a text string or function.
1144//   lead_in_sample = Factor to increase sample rate in the lead-in section.  Default: 10
1145//   internal = if true make internal threads.  The only effect this has is to change how the extrusion lead-in section are formed. When true, the extrusion scales towards the outside; when false, it scales towards the inside.  Default: false
1146//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `CENTER`
1147//   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
1148//   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#subsection-orient).  Default: `UP`
1149// Example:
1150//   poly = [[-10,0], [-3,-5], [3,-5], [10,0], [0,-30]];
1151//   spiral_sweep(poly, h=200, r=50, turns=3, $fn=36);
1152_leadin_ogive=function (x,L) 
1153     let( minscale = .05,
1154          r=(L^2+(1-minscale^2))/2/(1-minscale),
1155          scale = sqrt(r^2-(L*(1-x))^2) -(r-1)
1156     )
1157     x>1 ? [1,1]
1158   : x<0 ? [lerp(minscale,1,.25),0] 
1159   : [lerp(scale,1,.25),scale];     
1160
1161_leadin_cut = function(x,L) x>0 ? [1,1] : [1,0];
1162
1163_leadin_sqrt = function(x,L)
1164     let(end=0.05)   // Smallest scale at the end
1165     x>1 ? [1,1]
1166   : x<0 ? [lerp(end,1,.25),0]   
1167   : let(  
1168          s = sqrt(x + end^2 * (1-x))
1169     )
1170     [lerp(s,1,.25),s];    // thread width scale, thread height scale
1171
1172_leadin_linear = function(x,L)
1173     let(minscale=.1)
1174     x>1 ? [1,1]
1175   : x<0 ? [lerp(minscale,1,.25),0]
1176   : let(scale = lerp(minscale,1,x))
1177     [lerp(scale,1,.25),scale];
1178
1179_lead_in_table = [
1180     ["default", _leadin_sqrt],
1181     ["sqrt", _leadin_sqrt],
1182     ["cut", _leadin_cut],
1183     ["smooth", _leadin_ogive],
1184     ["linear", _leadin_linear]
1185];
1186
1187      
1188function _ss_polygon_r(N,theta) =
1189        let( alpha = 360/N )
1190        cos(alpha/2)/(cos(posmod(theta,alpha)-alpha/2));
1191function spiral_sweep(poly, h, r, turns=1, taper, r1, r2, d, d1, d2, internal=false,
1192                      lead_in_shape,lead_in_shape1, lead_in_shape2,
1193                      lead_in, lead_in1, lead_in2,
1194                      lead_in_ang, lead_in_ang1, lead_in_ang2,
1195                      height,l,length,
1196                      lead_in_sample = 10,
1197                      anchor=CENTER, spin=0, orient=UP) =
1198    assert(is_num(turns) && turns != 0, "turns must be a nonzero number")
1199    assert(all_positive([h]), "Spiral height must be a positive number")
1200    let(
1201        dir = sign(turns),
1202        r1 = get_radius(r1=r1, r=r, d1=d1, d=d),
1203        r2 = get_radius(r1=r2, r=r, d1=d2, d=d),
1204        bounds = pointlist_bounds(poly),
1205        yctr = (bounds[0].y+bounds[1].y)/2,
1206        xmin = bounds[0].x,
1207        xmax = bounds[1].x,
1208        poly = path3d(clockwise_polygon(poly)),
1209        sides = segs(max(r1,r2)),
1210        ang_step = 360/sides,
1211        turns = abs(turns),
1212        lead_in1 = first_defined([lead_in1, lead_in]),
1213        lead_in2 = first_defined([lead_in1, lead_in]),        
1214        lead_in_ang1 =
1215                      let(
1216                           user_ang = first_defined([lead_in_ang1,lead_in_ang])
1217                      )
1218                      assert(is_undef(user_ang) || is_undef(lead_in1), "Cannot define lead_in/lead_in1 by both length and angle")
1219                      is_def(user_ang) ? user_ang : default(lead_in1,0)*360/(2*PI*r1),
1220        lead_in_ang2 =
1221                      let(
1222                           user_ang = first_defined([lead_in_ang2,lead_in_ang])
1223                      )
1224                      assert(is_undef(user_ang) || is_undef(lead_in2), "Cannot define lead_in/lead_in2 by both length and angle")
1225                      is_def(user_ang) ? user_ang : default(lead_in2,0)*360/(2*PI*r2),
1226        minang = -max(0,lead_in_ang1),
1227        maxang = 360*turns + max(0,lead_in_ang2),
1228        cut_ang1 = minang+abs(lead_in_ang1),
1229        cut_ang2 = maxang-abs(lead_in_ang1),        
1230        lead_in_shape1 = first_defined([lead_in_shape1, lead_in_shape, "default"]),
1231        lead_in_shape2 = first_defined([lead_in_shape2, lead_in_shape, "default"]),             
1232        lead_in_func1 = is_func(lead_in_shape1) ? lead_in_shape1
1233                      : assert(is_string(lead_in_shape1),"lead_in_shape/lead_in_shape1 must be a function or string")
1234                        let(ind = search([lead_in_shape1], _lead_in_table,0)[0])
1235                        assert(ind!=[],str("Unknown lead_in_shape, \"",lead_in_shape1,"\""))
1236                        _lead_in_table[ind[0]][1],
1237        lead_in_func2 = is_func(lead_in_shape2) ? lead_in_shape2
1238                      : assert(is_string(lead_in_shape2),"lead_in_shape/lead_in_shape2 must be a function or string")
1239                        let(ind = search([lead_in_shape2], _lead_in_table,0)[0])
1240                        assert(ind!=[],str("Unknown lead_in_shape, \"",lead_in_shape2,"\""))
1241                        _lead_in_table[ind[0]][1]
1242    )
1243    assert( cut_ang1<cut_ang2, "Tapers are too long to fit")
1244    assert( all_positive([r1,r2]), "Diameter/radius must be positive")
1245    let(
1246  
1247        // This complicated sampling scheme is designed to ensure that faceting always starts at angle zero
1248        // for alignment with cylinders, and there is always a facet boundary at the $fn specified locations, 
1249        // regardless of what kind of subsampling occurs for tapers.
1250        orig_anglist = [
1251            if (minang<0) minang,
1252            each reverse([for(ang = [-ang_step:-ang_step:minang+EPSILON]) ang]),
1253            for(ang = [0:ang_step:maxang-EPSILON]) ang,
1254            maxang
1255        ],
1256        anglist = [
1257           for(a=orig_anglist) if (a<cut_ang1-EPSILON) a,
1258           cut_ang1,
1259           for(a=orig_anglist) if (a>cut_ang1+EPSILON && a<cut_ang2-EPSILON) a,
1260           cut_ang2,
1261           for(a=orig_anglist) if (a>cut_ang2+EPSILON) a
1262        ],
1263        interp_ang = [
1264                      for(i=idx(anglist,e=-2)) 
1265                          each lerpn(anglist[i],anglist[i+1],
1266                                         (lead_in_ang1!=0 && anglist[i+1]<=cut_ang1) || (lead_in_ang2!=0 && anglist[i]>=cut_ang2)
1267                                            ? ceil((anglist[i+1]-anglist[i])/ang_step*lead_in_sample)
1268                                            : 1,
1269                                     endpoint=false),
1270                      last(anglist)
1271                     ],
1272        skewmat = affine3d_skew_xz(xa=atan2(r2-r1,h)),
1273        points = [
1274            for (a = interp_ang) let (
1275                hsc = a<cut_ang1 ? lead_in_func1((a-minang)/abs(lead_in_ang1),abs(lead_in_ang1)*2*PI*r1/360)
1276                    : a>cut_ang2 ? lead_in_func2((maxang-a)/abs(lead_in_ang2),abs(lead_in_ang2)*2*PI*r2/360)
1277                    : [1,1],
1278                u = a/(360*turns), 
1279                r = lerp(r1,r2,u),
1280                mat = affine3d_zrot(dir*a)
1281                    * affine3d_translate([_ss_polygon_r(sides,dir*a)*r, 0, h * (u-0.5)])
1282                    * affine3d_xrot(90)
1283                    * skewmat
1284                    * scale([hsc.y,hsc.x,1], cp=[internal ? xmax : xmin, yctr, 0]),
1285                pts = apply(mat, poly)
1286            ) pts
1287        ],
1288        vnf = vnf_vertex_array(
1289            points, col_wrap=true, caps=true, reverse=dir>0,
1290        //    style=higbee1>0 || higbee2>0 ? "quincunx" : "alt"
1291            style="convex"
1292        ),
1293        vnf2 = vnf_triangulate(vnf)
1294    )
1295    reorient(anchor,spin,orient, vnf=vnf2, r1=r1, r2=r2, l=h, p=vnf2);
1296
1297
1298
1299module spiral_sweep(poly, h, r, turns=1, taper, r1, r2, d, d1, d2, internal=false,
1300                    lead_in_shape,lead_in_shape1, lead_in_shape2,
1301                    lead_in, lead_in1, lead_in2,
1302                    lead_in_ang, lead_in_ang1, lead_in_ang2,
1303                    height,l,length,
1304                    lead_in_sample=10,
1305                    anchor=CENTER, spin=0, orient=UP)
1306{
1307    vnf = spiral_sweep(poly=poly, h=h, r=r, turns=turns, r1=r1, r2=r2, d=d, d1=d1, d2=d2, internal=internal,
1308                       lead_in_shape=lead_in_shape,lead_in_shape1=lead_in_shape1, lead_in_shape2=lead_in_shape2,
1309                       lead_in=lead_in, lead_in1=lead_in1, lead_in2=lead_in2,
1310                       lead_in_ang=lead_in_ang, lead_in_ang1=lead_in_ang1, lead_in_ang2=lead_in_ang2,
1311                       height=height,l=length,length=length,
1312                       lead_in_sample=lead_in_sample);
1313    h = one_defined([h,height,length,l],"h,height,length,l");
1314    r1 = get_radius(r1=r1, r=r, d1=d1, d=d);
1315    r2 = get_radius(r1=r2, r=r, d1=d2, d=d);
1316    lead_in1 = u_mul(first_defined([lead_in1,lead_in]),1/(2*PI*r1));
1317    lead_in2 = u_mul(first_defined([lead_in2,lead_in]),1/(2*PI*r2));
1318    lead_in_ang1 = first_defined([lead_in_ang1,lead_in_ang]);
1319    lead_in_ang2 = first_defined([lead_in_ang2,lead_in_ang]);
1320    extra_turns = max(0,first_defined([lead_in1,lead_in_ang1,0]))+max(0,first_defined([lead_in2,lead_in_ang2,0]));
1321    attachable(anchor,spin,orient, r1=r1, r2=r2, l=h) {
1322        vnf_polyhedron(vnf, convexity=ceil(2*(abs(turns)+extra_turns)));
1323        children();
1324    }
1325}
1326
1327
1328
1329// Function&Module: path_sweep()
1330// Synopsis: Sweep a 2d polygon path along a 2d or 3d path. 
1331// SynTags: VNF, Geom
1332// Topics: Extrusion, Sweep, Paths
1333// See Also: sweep_attach(), linear_sweep(), rotate_sweep(), sweep(), spiral_sweep(), path_sweep2d(), offset_sweep()
1334// Usage: As module
1335//   path_sweep(shape, path, [method], [normal=], [closed=], [twist=], [twist_by_length=], [symmetry=], [scale=], [scale_by_length=], [last_normal=], [tangent=], [uniform=], [relaxed=], [caps=], [style=], [convexity=], [anchor=], [cp=], [spin=], [orient=], [atype=]) [ATTACHMENTS];
1336// Usage: As function
1337//   vnf = path_sweep(shape, path, [method], [normal=], [closed=], [twist=], [twist_by_length=], [symmetry=], [scale=], [scale_by_length=], [last_normal=], [tangent=], [uniform=], [relaxed=], [caps=], [style=], [transforms=], [anchor=], [cp=], [spin=], [orient=], [atype=]);
1338// Description:
1339//   Takes as input `shape`, a 2D polygon path (list of points), and `path`, a 2d or 3d path (also a list of points)
1340//   and constructs a polyhedron by sweeping the shape along the path. When run as a module returns the polyhedron geometry.
1341//   When run as a function returns a VNF by default or if you set `transforms=true` then it returns a list of transformations suitable as input to `sweep`.
1342//   .
1343//   The sweeping process places one copy of the shape for each point in the path.  The origin in `shape` is translated to
1344//   the point in `path`.  The normal vector of the shape, which points in the Z direction, is aligned with the tangent
1345//   vector for the path, so this process is constructing a shape whose normal cross sections are equal to your specified shape.
1346//   If you do not supply a list of tangent vectors then an approximate tangent vector is computed
1347//   based on the path points you supply using {{path_tangents()}}.
1348// Figure(3D,Big,VPR=[70,0,345],VPD=20,VPT=[5.5,10.8,-2.7],NoScales): This example shows how the shape, in this case the quadrilateral defined by `[[0, 0], [0, 1], [0.25, 1], [1, 0]]`, appears as the cross section of the swept polyhedron.  The blue line shows the path.  The normal vector to the shape is shown in black; it is based at the origin and points upwards in the Z direction.  The sweep aligns this normal vector with the blue path tangent, which in this case, flips the shape around.  Note that for a 2D path like this one, the Y direction in the shape is mapped to the Z direction in the sweep.
1349//   tri= [[0, 0], [0, 1], [.25,1], [1, 0]];
1350//   path = arc(r=5,n=81,angle=[-20,65]);
1351//   % path_sweep(tri,path);
1352//   T = path_sweep(tri,path,transforms=true);
1353//   color("red")for(i=[0:20:80]) stroke(apply(T[i],path3d(tri)),width=.1,closed=true);
1354//   color("blue")stroke(path3d(arc(r=5,n=101,angle=[-20,80])),width=.1,endcap2="arrow2");
1355//   color("red")stroke([path3d(tri)],width=.1);
1356//   stroke([CENTER,UP], width=.07,endcap2="arrow2",color="black");
1357// Continues:
1358//   In the figure you can see that the swept polyhedron, shown in transparent gray, has the quadrilateral as its cross
1359//   section.  The quadrilateral is positioned perpendicular to the path, which is shown in blue, so that the normal
1360//   vector for the quadrilateral is parallel to the tangent vector for the path.  The origin for the shape is the point
1361//   which follows the path.  For a 2D path, the Y axis of the shape is mapped to the Z axis and in this case,
1362//   pointing the quadrilateral's normal vector (in black) along the tangent line of
1363//   the path, which is going in the direction of the blue arrow, requires that the quadrilateral be "turned around".  If we
1364//   reverse the order of points in the path we get a different result:
1365// Figure(3D,Big,VPR=[70,0,20],VPD=20,VPT=[1.25,9.25,-2.65],NoScales): The same sweep operation with the path traveling in the opposite direction.  Note that in order to line up the normal correctly, the shape is reversed compared to Figure 1, so the resulting sweep looks quite different.
1366//   tri= [[0, 0], [0, 1], [.25,1], [1, 0]];
1367//   path = reverse(arc(r=5,n=81,angle=[-20,65]));
1368//   % path_sweep(tri,path);
1369//   T = path_sweep(tri,path,transforms=true);
1370//   color("red")for(i=[0:20:80]) stroke(apply(T[i],path3d(tri)),width=.1,closed=true);
1371//   color("blue")stroke(reverse(path3d(arc(r=5,n=101,angle=[-20-15,65]))),width=.1,endcap2="arrow2");
1372//   color("red")stroke([path3d(tri)],width=.1);
1373//   stroke([CENTER,UP], width=.07,endcap2="arrow2",color="black");
1374// Continues:
1375//   If your shape is too large for the curves in the path you can create a situation where the shapes cross each
1376//   other.  This results in an invalid polyhedron, which may appear OK when previewed or rendered alone, but will give rise
1377//   to cryptic CGAL errors when rendered with a second object in your model.  You may be able to use {{path_sweep2d()}}
1378//   to produce a valid model in cases like this.  You can debug models like this using the `profiles=true` option which will show all
1379//   the cross sections in your polyhedron.  If any of them intersect, the polyhedron will be invalid.
1380// Figure(3D,Big,VPR=[47,0,325],VPD=23,VPT=[6.8,4,-3.8],NoScales): We have scaled the path to an ellipse and show a large triangle as the shape.  The triangle is sometimes bigger than the local radius of the path, leading to an invalid polyhedron, which you can identify because the red lines cross in the middle.
1381//   tri= scale([4.5,2.5],[[0, 0], [0, 1], [1, 0]]);
1382//   path = xscale(1.5,arc(r=5,n=81,angle=[-70,70]));
1383//   % path_sweep(tri,path);
1384//   T = path_sweep(tri,path,transforms=true);
1385//   color("red")for(i=[0:20:80]) stroke(apply(T[i],path3d(tri)),width=.1,closed=true);
1386//   color("blue")stroke(path3d(xscale(1.5,arc(r=5,n=81,angle=[-70,80]))),width=.1,endcap2="arrow2");
1387// Continues:
1388//   During the sweep operation the shape's normal vector aligns with the tangent vector of the path.  Note that
1389//   this leaves an ambiguity about how the shape is rotated as it sweeps along the path.
1390//   For 2D paths, this ambiguity is resolved by aligning the Y axis of the shape to the Z axis of the swept polyhedron.
1391//   You can force the  shape to twist as it sweeps along the path using the `twist` parameter, which specifies the total
1392//   number of degrees to twist along the whole swept polyhedron.  This produces a result like the one shown below.
1393// Figure(3D,Big,VPR=[66,0,14],VPD=20,VPT=[3.4,4.5,-0.8]): The shape twists as we sweep.  Note that it still aligns the origin in the shape with the path, and still aligns the normal vector with the path tangent vector.
1394//   tri= [[0, 0], [0, 1], [.25,1],[1, 0]];
1395//   path = arc(r=5,n=81,angle=[-20,65]);
1396//   % path_sweep(tri,path,twist=-60);
1397//   T = path_sweep(tri,path,transforms=true,twist=-60);
1398//   color("red")for(i=[0:20:80]) stroke(apply(T[i],path3d(tri)),width=.1,closed=true);
1399//   color("blue")stroke(path3d(arc(r=5,n=101,angle=[-20,80])),width=.1,endcap2="arrow2");
1400// Continues:
1401//   The `twist` argument adds the specified number of degrees of twist into the model, and it may be positive or
1402//   negative.  When `closed=true` the starting shape and ending shape must match to avoid a sudden extreme twist at the
1403//   joint.  By default `twist` is therefore required to be a multiple of 360.  However, if your shape has rotational
1404//   symmetry, this requirement is overly strict.  You can specify the symmetry using the `symmetry` argument, and then
1405//   you can choose smaller twists consistent with the specified symmetry.  The symmetry argument gives the number of
1406//   rotations that map the shape exactly onto itself, so a pentagon has 5-fold symmetry.  This argument is only valid
1407//   for closed sweeps.  When you specify symmetry, the twist must be a multiple of 360/symmetry.
1408//   .
1409//   The twist is normally spread uniformly along your shape based on the path length.  If you set `twist_by_length` to
1410//   false then the twist will be uniform based on the point count of your path.  Twisted shapes will produce twisted
1411//   faces, so if you want them to look good you should use lots of points on your path and also lots of points on the
1412//   shape.  If your shape is a simple polygon, use {{subdivide_path()}} to increase
1413//   the number of points.
1414//   .
1415//   As noted above, the sweep process has an ambiguity regarding the twist.  For 2D paths it is easy to resolve this
1416//   ambiguity by aligning the Y axis in the shape to the Z axis in the swept polyhedron.  When the path is
1417//   three-dimensional, things become more complex.  It is no longer possible to use a simple alignment rule like the
1418//   one we use in 2D.  You may find that the shape rotates unexpectedly around its axis as it traverses the path.  The
1419//   `method` parameter allows you to specify how the shapes are aligned, resulting in different twist in the resulting
1420//   polyhedron.  You can choose from three different methods for selecting the rotation of your shape.  None of these
1421//   methods will produce good, or even valid, results on all inputs, so it is important to select a suitable method.
1422//   .
1423//   The three methods you can choose using the `method` parameter are:
1424//   .
1425//   The "incremental" method (the default) works by adjusting the shape at each step by the minimal rotation that makes the shape normal to the tangent
1426//   at the next point.  This method is robust in that it always produces a valid result for well-behaved paths with sufficiently high
1427//   sampling.  Unfortunately, it can produce a large amount of undesirable twist.  When constructing a closed shape this algorithm in
1428//   its basic form provides no guarantee that the start and end shapes match up.  To prevent a sudden twist at the last segment,
1429//   the method calculates the required twist for a good match and distributes it over the whole model (as if you had specified a
1430//   twist amount).  If you specify `symmetry` this may allow the algorithm to choose a smaller twist for this alignment.
1431//   To start the algorithm, we need an initial condition.  This is supplied by
1432//   using the `normal` argument to give a direction to align the Y axis of your shape.  By default the normal points UP if the path
1433//   makes an angle of 45 deg or less with the xy plane and it points BACK if the path makes a higher angle with the XY plane.  You
1434//   can also supply `last_normal` which provides an ending orientation constraint.  Be aware that the curve may still exhibit
1435//   twisting in the middle.  This method is the default because it is the most robust, not because it generally produces the best result.
1436//   .
1437//   The "natural" method works by computing the Frenet frame at each point on the path.  This is defined by the tangent to the curve and
1438//   the normal which lies in the plane defined by the curve at each point.  This normal points in the direction of curvature of the curve.
1439//   The result is a very well behaved set of shape positions without any unexpected twisting&mdash;as long as the curvature never falls to zero.  At a
1440//   point of zero curvature (a flat point), the curve does not define a plane and the natural normal is not defined.  Furthermore, even if
1441//   you skip over this troublesome point so the normal is defined, it can change direction abruptly when the curvature is zero, leading to
1442//   a nasty twist and an invalid model.  A simple example is a circular arc joined to another arc that curves the other direction.  Note
1443//   that the X axis of the shape is aligned with the normal from the Frenet frame.
1444//   .
1445//   The "manual" method allows you to specify your desired normal either globally with a single vector, or locally with
1446//   a list of normal vectors for every path point.  The normal you supply is projected to be orthogonal to the tangent to the
1447//   path and the Y direction of your shape will be aligned with the projected normal.  (Note this is different from the "natural" method.)
1448//   Careless choice of a normal may result in a twist in the shape, or an error if your normal is parallel to the path tangent.
1449//   If you set `relax=true` then the condition that the cross sections are orthogonal to the path is relaxed and the swept object
1450//   uses the actual specified normal.  In this case, the tangent is projected to be orthogonal to your supplied normal to define
1451//   the cross section orientation.  Specifying a list of normal vectors gives you complete control over the orientation of your
1452//   cross sections and can be useful if you want to position your model to be on the surface of some solid.
1453//   .
1454//   You can also apply scaling to the profile along the path.  You can give a list of scalar scale factors or a list of 2-vector scale. 
1455//   In the latter scale the x and y scales of the profile are scaled separately before the profile is placed onto the path.  For non-closed
1456//   paths you can also give a single scale value or a 2-vector which is treated as the final scale.  The intermediate sections
1457//   are then scaled by linear interpolation either relative to length (if scale_by_length is true) or by point count otherwise.  
1458//   .
1459//   You can use set `transforms` to true to return a list of transformation matrices instead of the swept shape.  In this case, you can
1460//   often omit shape entirely.  The exception is when `closed=true` and you are using the "incremental" method.  In this case, `path_sweep`
1461//   uses the shape to correct for twist when the shape closes on itself, so you must include a valid shape.
1462// Arguments:
1463//   shape = A 2D polygon path or region describing the shape to be swept.
1464//   path = 2D or 3D path giving the path to sweep over
1465//   method = one of "incremental", "natural" or "manual".  Default: "incremental"
1466//   ---
1467//   normal = normal vector for initializing the incremental method, or for setting normals with method="manual".  Default: UP if the path makes an angle lower than 45 degrees to the xy plane, BACK otherwise.
1468//   closed = path is a closed loop.  Default: false
1469//   twist = amount of twist to add in degrees.  For closed sweeps must be a multiple of 360/symmetry.  Default: 0
1470//   twist_by_length = if true then interpolate twist based on the path length of the path. If false interoplate based on point count.  Default: true
1471//   symmetry = symmetry of the shape when closed=true.  Allows the shape to join with a 360/symmetry rotation instead of a full 360 rotation.  Default: 1
1472//   scale = Amount to scale the profiles.  If you give a scalar the scale starts at 1 and ends at your specified value. The same is true for a 2-vector, but x and y are scaled separately.   You can also give a vector of values, one for each path point, and you can give a list of 2-vectors that give the x and y scales of your profile for every point on the path (a Nx2 matrix for a path of length N.  Default: 1 (no scaling)
1473//   scale_by_length = if true then interpolate scale based on the path length of the path. If false interoplate based on point count.  Default: true
1474//   last_normal = normal to last point in the path for the "incremental" method.  Constrains the orientation of the last cross section if you supply it.
1475//   uniform = if set to false then compute tangents using the uniform=false argument, which may give better results when your path is non-uniformly sampled.  This argument is passed to {{path_tangents()}}.  Default: true
1476//   tangent = a list of tangent vectors in case you need more accuracy (particularly at the end points of your curve)
1477//   relaxed = set to true with the "manual" method to relax the orthogonality requirement of cross sections to the path tangent.  Default: false
1478//   caps = Can be a boolean or vector of two booleans.  Set to false to disable caps at the two ends.  Default: true
1479//   style = vnf_vertex_array style.  Default: "min_edge"
1480//   profiles = if true then display all the cross section profiles instead of the solid shape.  Can help debug a sweep.  (module only) Default: false
1481//   width = the width of lines used for profile display.  (module only) Default: 1
1482//   transforms = set to true to return transforms instead of a VNF.  These transforms can be manipulated and passed to sweep().  (function only)  Default: false.
1483//   convexity = convexity parameter for polyhedron().  (module only)  Default: 10
1484//   anchor = Translate so anchor point is at the origin. Default: "origin"
1485//   spin = Rotate this many degrees around Z axis after anchor. Default: 0
1486//   orient = Vector to rotate top towards after spin
1487//   atype  = Select "hull" or "intersect" anchor types.  Default: "hull"
1488//   cp = Centerpoint for determining "intersect" anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
1489// Side Effects:
1490//   `$sweep_path` is set to the path thd defining the swept object
1491//   `$sweep_shape` is set to the shape being swept
1492//   `$sweep_closed` is true if the sweep is closed and false otherwise
1493//   `$sweep_transforms` is set to the array of transformation matrices that define the swept object.
1494//   `$sweep_scales` is set to the array of scales that were applied at each point to create the swept object.
1495//   `$sweep_twist` set to a scalar value giving the total twist across the path sweep object.
1496// Anchor Types:
1497//   "hull" = Anchors to the virtual convex hull of the shape.
1498//   "intersect" = Anchors to the surface of the shape.
1499// Named Anchors:
1500//   "origin" = The native position of the shape
1501//   "start" = When `closed==false`, the origin point of the shape, on the starting face of the object
1502//   "end" = When `closed==false`, the origin point of the shape, on the ending face of the object
1503//   "start-centroid" = When `closed==false`, the centroid of the shape, on the starting face of the object
1504//   "end-centroid" = When `closed==false`, the centroid of the shape, on the ending face of the object
1505// Example(NoScales): A simple sweep of a square along a sine wave:
1506//   path = [for(theta=[-180:5:180]) [theta/10, 10*sin(theta)]];
1507//   sq = square(6,center=true);
1508//   path_sweep(sq,path);
1509// Example(NoScales): If the square is not centered, then we get a different result because the shape is in a different place relative to the origin:
1510//   path = [for(theta=[-180:5:180]) [theta/10, 10*sin(theta)]];
1511//   sq = square(6);
1512//   path_sweep(sq,path);
1513// Example(Med,VPR=[34,0,8],NoScales): It may not be obvious, but the polyhedron in the previous example is invalid.  It will eventually give CGAL errors when you combine it with other shapes.  To see this, set profiles to true and look at the left side.  The profiles cross each other and intersect.  Any time this happens, your polyhedron is invalid, even if it seems to be working at first.  Another observation from the profile display is that we have more profiles than needed over a lot of the shape, so if the model is slow, using fewer profiles in the flat portion of the curve might speed up the calculation.
1514//   path = [for(theta=[-180:5:180]) [theta/10, 10*sin(theta)]];
1515//   sq = square(6);
1516//   path_sweep(sq,path,profiles=true,width=.1,$fn=8);
1517// Example(2D): We'll use this shape in several examples
1518//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1519//   polygon(ushape);
1520// Example(NoScales): Sweep along a clockwise elliptical arc, using default "incremental" method.
1521//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1522//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[180,00], r=30));  // Clockwise
1523//   path_sweep(ushape, path3d(elliptic_arc));
1524// Example(NoScales): Sweep along a counter-clockwise elliptical arc.  Note that the orientation of the shape flips.
1525//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1526//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[0,180], r=30));   // Counter-clockwise
1527//   path_sweep(ushape, path3d(elliptic_arc));
1528// Example(NoScales): Sweep along a clockwise elliptical arc, using "natural" method, which lines up the X axis of the shape with the direction of curvature.  This means the X axis will point inward, so a counterclockwise arc gives:
1529//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1530//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[0,180], r=30));  // Counter-clockwise
1531//   path_sweep(ushape, elliptic_arc, method="natural");
1532// Example(NoScales): Sweep along a clockwise elliptical arc, using "natural" method.  If the curve is clockwise then the shape flips upside-down to align the X axis.
1533//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1534//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[180,0], r=30));  // Clockwise
1535//   path_sweep(ushape, path3d(elliptic_arc), method="natural");
1536// Example(NoScales): Sweep along a clockwise elliptical arc, using "manual" method.  You can orient the shape in a direction you choose (subject to the constraint that the profiles remain normal to the path):
1537//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1538//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[180,0], r=30));  // Clockwise
1539//   path_sweep(ushape, path3d(elliptic_arc), method="manual", normal=UP+RIGHT);
1540// Example(NoScales): Here we changed the ellipse to be more pointy, and with the same results as above we get a shape with an irregularity in the middle where it maintains the specified direction around the point of the ellipse.  If the ellipse were more pointy, this would result in a bad polyhedron:
1541//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1542//   elliptic_arc = yscale(2, p=arc($fn=64,angle=[180,0], r=30));  // Clockwise
1543//   path_sweep(ushape, path3d(elliptic_arc), method="manual", normal=UP+RIGHT);
1544// Example(NoScales): It is easy to produce an invalid shape when your path has a smaller radius of curvature than the width of your shape.  The exact threshold where the shape becomes invalid depends on the density of points on your path.  The error may not be immediately obvious, as the swept shape appears fine when alone in your model, but adding a cube to the model reveals the problem.  In this case the pentagon is turned so its longest direction points inward to create the singularity.
1545//   qpath = [for(x=[-3:.01:3]) [x,x*x/1.8,0]];
1546//   // Prints 0.9, but we use pentagon with radius of 1.0 > 0.9
1547//   echo(radius_of_curvature = 1/max(path_curvature(qpath)));
1548//   path_sweep(apply(rot(90),pentagon(r=1)), qpath, normal=BACK, method="manual");
1549//   cube(0.5);    // Adding a small cube forces a CGAL computation which reveals
1550//                 // the error by displaying nothing or giving a cryptic message
1551// Example(NoScales): Using the `relax` option we allow the profiles to deviate from orthogonality to the path.  This eliminates the crease that broke the previous example because the sections are all parallel to each other.
1552//   qpath = [for(x=[-3:.01:3]) [x,x*x/1.8,0]];
1553//   path_sweep(apply(rot(90),pentagon(r=1)), qpath, normal=BACK, method="manual", relaxed=true);
1554//   cube(0.5);    // Adding a small cube is not a problem with this valid model
1555// Example(Med,VPR=[16,0,100],VPT=[0.05,0.6,0.6],VPD=25,NoScales): Using the `profiles=true` option can help debug bad polyhedra such as this one.  If any of the profiles intersect or cross each other, the polyhedron will be invalid.  In this case, you can see these intersections in the middle of the shape, which may give insight into how to fix your shape.   The profiles may also help you identify cases with a valid polyhedron where you have more profiles than needed to adequately define the shape.
1556//   tri= scale([4.5,2.5],[[0, 0], [0, 1], [1, 0]]);
1557//   path = left(4,xscale(1.5,arc(r=5,n=25,angle=[-70,70])));
1558//   path_sweep(tri,path,profiles=true,width=.1);
1559// Example(NoScales):  This 3d arc produces a result that twists to an undefined angle.  By default the incremental method sets the starting normal to UP, but the ending normal is unconstrained.
1560//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1561//   arc = yrot(37, p=path3d(arc($fn=64, r=30, angle=[0,180])));
1562//   path_sweep(ushape, arc, method="incremental");
1563// Example(NoScales): You can constrain the last normal as well.  Here we point it right, which produces a nice result.
1564//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1565//   arc = yrot(37, p=path3d(arc($fn=64, r=30, angle=[0,180])));
1566//   path_sweep(ushape, arc, method="incremental", last_normal=RIGHT);
1567// Example(NoScales): Here we constrain the last normal to UP.  Be aware that the behavior in the middle is unconstrained.
1568//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1569//   arc = yrot(37, p=path3d(arc($fn=64, r=30, angle=[0,180])));
1570//   path_sweep(ushape, arc, method="incremental", last_normal=UP);
1571// Example(NoScales): The "natural" method produces a very different result
1572//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1573//   arc = yrot(37, p=path3d(arc($fn=64, r=30, angle=[0,180])));
1574//   path_sweep(ushape, arc, method="natural");
1575// Example(NoScales): When the path starts at an angle of more that 45 deg to the xy plane the initial normal for "incremental" is BACK.  This produces the effect of the shape rising up out of the xy plane.  (Using UP for a vertical path is invalid, hence the need for a split in the defaults.)
1576//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1577//   arc = xrot(75, p=path3d(arc($fn=64, r=30, angle=[0,180])));
1578//   path_sweep(ushape, arc, method="incremental");
1579// Example(NoScales): Adding twist
1580//   // Counter-clockwise
1581//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[0,180], r=3));
1582//   path_sweep(pentagon(r=1), path3d(elliptic_arc), twist=72);
1583// Example(NoScales): Closed shape
1584//   ellipse = xscale(2, p=circle($fn=64, r=3));
1585//   path_sweep(pentagon(r=1), path3d(ellipse), closed=true);
1586// Example(NoScales): Closed shape with added twist
1587//   ellipse = xscale(2, p=circle($fn=64, r=3));
1588//   // Looks better with finer sampling
1589//   pentagon = subdivide_path(pentagon(r=1), 30);
1590//   path_sweep(pentagon, path3d(ellipse),
1591//              closed=true, twist=360);
1592// Example(NoScales): The last example was a lot of twist.  In order to use less twist you have to tell `path_sweep` that your shape has symmetry, in this case 5-fold.  Mobius strip with pentagon cross section:
1593//   ellipse = xscale(2, p=circle($fn=64, r=3));
1594//   // Looks better with finer sampling
1595//   pentagon = subdivide_path(pentagon(r=1), 30);
1596//   path_sweep(pentagon, path3d(ellipse), closed=true,
1597//              symmetry = 5, twist=2*360/5);
1598// Example(Med,NoScales): A helical path reveals the big problem with the "incremental" method: it can introduce unexpected and extreme twisting.  (Note helix example came from list-comprehension-demos)
1599//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1600//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1601//                         200 * (1 - t)];
1602//   helix_steps = 200;
1603//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1604//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1605//   path_sweep(ushape, helix);
1606// Example(Med,NoScales): You can constrain both ends, but still the twist remains:
1607//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1608//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1609//                         200 * (1 - t)];
1610//   helix_steps = 200;
1611//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1612//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1613//   path_sweep(ushape, helix, normal=UP, last_normal=UP);
1614// Example(Med,NoScales): Even if you manually guess the amount of twist and remove it, the result twists one way and then the other:
1615//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1616//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1617//                         200 * (1 - t)];
1618//   helix_steps = 200;
1619//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1620//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1621//   path_sweep(ushape, helix, normal=UP, last_normal=UP, twist=360);
1622// Example(Med,NoScales): To get a good result you must use a different method.
1623//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1624//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1625//                         200 * (1 - t)];
1626//   helix_steps = 200;
1627//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1628//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1629//   path_sweep(ushape, helix, method="natural");
1630// Example(Med,NoScales): Note that it may look like the shape above is flat, but the profiles are very slightly tilted due to the nonzero torsion of the curve.  If you want as flat as possible, specify it so with the "manual" method:
1631//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1632//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1633//                         200 * (1 - t)];
1634//   helix_steps = 200;
1635//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1636//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1637//   path_sweep(ushape, helix, method="manual", normal=UP);
1638// Example(Med,NoScales): What if you want to angle the shape inward?  This requires a different normal at every point in the path:
1639//   function helix(t) = [(t / 1.5 + 0.5) * 30 * cos(6 * 360 * t),
1640//                        (t / 1.5 + 0.5) * 30 * sin(6 * 360 * t),
1641//                         200 * (1 - t)];
1642//   helix_steps = 200;
1643//   helix = [for (i=[0:helix_steps]) helix(i/helix_steps)];
1644//   normals = [for(i=[0:helix_steps]) [-cos(6*360*i/helix_steps), -sin(6*360*i/helix_steps), 2.5]];
1645//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1646//   path_sweep(ushape, helix, method="manual", normal=normals);
1647// Example(NoScales): When using "manual" it is important to choose a normal that works for the whole path, producing a consistent result.  Here we have specified an upward normal, and indeed the shape is pointed up everywhere, but two abrupt transitional twists render the model invalid.
1648//   yzcircle = yrot(90,p=path3d(circle($fn=64, r=30)));
1649//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1650//   path_sweep(ushape, yzcircle, method="manual", normal=UP, closed=true);
1651// Example(NoScales): The "natural" method will introduce twists when the curvature changes direction.  A warning is displayed.
1652//   arc1 = path3d(arc(angle=90, r=30));
1653//   arc2 = xrot(-90, cp=[0,30],p=path3d(arc(angle=[90,180], r=30)));
1654//   two_arcs = path_merge_collinear(concat(arc1,arc2));
1655//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1656//   path_sweep(ushape, two_arcs, method="natural");
1657// Example(NoScales): The only simple way to get a good result is the "incremental" method:
1658//   arc1 = path3d(arc(angle=90, r=30));
1659//   arc2 = xrot(-90, cp=[0,30],p=path3d(arc(angle=[90,180], r=30)));
1660//   arc3 = apply( translate([-30,60,30])*yrot(90), path3d(arc(angle=[270,180], r=30)));
1661//   three_arcs = path_merge_collinear(concat(arc1,arc2,arc3));
1662//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1663//   path_sweep(ushape, three_arcs, method="incremental");
1664// Example(Med,NoScales): knot example from list-comprehension-demos, "incremental" method
1665//   function knot(a,b,t) =   // rolling knot
1666//        [ a * cos (3 * t) / (1 - b* sin (2 *t)),
1667//          a * sin( 3 * t) / (1 - b* sin (2 *t)),
1668//        1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))];
1669//   a = 0.8; b = sqrt (1 - a * a);
1670//   ksteps = 400;
1671//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
1672//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1673//   path_sweep(ushape, knot_path, closed=true, method="incremental");
1674// Example(Med,NoScales): knot example from list-comprehension-demos, "natural" method.  Which one do you like better?
1675//   function knot(a,b,t) =   // rolling knot
1676//        [ a * cos (3 * t) / (1 - b* sin (2 *t)),
1677//          a * sin( 3 * t) / (1 - b* sin (2 *t)),
1678//        1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))];
1679//   a = 0.8; b = sqrt (1 - a * a);
1680//   ksteps = 400;
1681//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
1682//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1683//   path_sweep(ushape, knot_path, closed=true, method="natural");
1684// Example(Med,NoScales): knot with twist.  Note if you twist it the other direction the center section untwists because of the natural twist there.  Also compare to the "incremental" method which has less twist in the center.
1685//   function knot(a,b,t) =   // rolling knot
1686//        [ a * cos (3 * t) / (1 - b* sin (2 *t)),
1687//          a * sin( 3 * t) / (1 - b* sin (2 *t)),
1688//        1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))];
1689//   a = 0.8; b = sqrt (1 - a * a);
1690//   ksteps = 400;
1691//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
1692//   path_sweep(subdivide_path(pentagon(r=12),30), knot_path, closed=true,
1693//              twist=-360*8, symmetry=5, method="natural");
1694// Example(Med,NoScales): twisted knot with twist distributed by path sample points instead of by length using `twist_by_length=false`
1695//   function knot(a,b,t) =   // rolling knot
1696//           [ a * cos (3 * t) / (1 - b* sin (2 *t)),
1697//             a * sin( 3 * t) / (1 - b* sin (2 *t)),
1698//           1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))];
1699//   a = 0.8; b = sqrt (1 - a * a);
1700//   ksteps = 400;
1701//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
1702//   path_sweep(subdivide_path(pentagon(r=12),30), knot_path, closed=true,
1703//              twist=-360*8, symmetry=5, method="natural", twist_by_length=false);
1704// Example(Big,NoScales): This torus knot example comes from list-comprehension-demos.  The knot lies on the surface of a torus.  When we use the "natural" method the swept figure is angled compared to the surface of the torus because the curve doesn't follow geodesics of the torus.
1705//   function knot(phi,R,r,p,q) =
1706//       [ (r * cos(q * phi) + R) * cos(p * phi),
1707//         (r * cos(q * phi) + R) * sin(p * phi),
1708//          r * sin(q * phi) ];
1709//   ushape = 3*[[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1710//   points = 50;       // points per loop
1711//   R = 400; r = 150;  // Torus size
1712//   p = 2;  q = 5;     // Knot parameters
1713//   %torus(r_maj=R,r_min=r);
1714//   k = max(p,q) / gcd(p,q) * points;
1715//   knot_path   = [ for (i=[0:k-1]) knot(360*i/k/gcd(p,q),R,r,p,q) ];
1716//   path_sweep(rot(90,p=ushape),knot_path,  method="natural", closed=true);
1717// Example(Big,NoScales): By computing the normal to the torus at the path we can orient the path to lie on the surface of the torus:
1718//   function knot(phi,R,r,p,q) =
1719//       [ (r * cos(q * phi) + R) * cos(p * phi),
1720//         (r * cos(q * phi) + R) * sin(p * phi),
1721//          r * sin(q * phi) ];
1722//   function knot_normal(phi,R,r,p,q) =
1723//       knot(phi,R,r,p,q)
1724//           - R*unit(knot(phi,R,r,p,q)
1725//               - [0,0, knot(phi,R,r,p,q)[2]]) ;
1726//   ushape = 3*[[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
1727//   points = 50;       // points per loop
1728//   R = 400; r = 150;  // Torus size
1729//   p = 2;  q = 5;     // Knot parameters
1730//   %torus(r_maj=R,r_min=r);
1731//   k = max(p,q) / gcd(p,q) * points;
1732//   knot_path   = [ for (i=[0:k-1]) knot(360*i/k/gcd(p,q),R,r,p,q) ];
1733//   normals = [ for (i=[0:k-1]) knot_normal(360*i/k/gcd(p,q),R,r,p,q) ];
1734//   path_sweep(ushape,knot_path,normal=normals, method="manual", closed=true);
1735// Example(NoScales): You can request the transformations and manipulate them before passing them on to sweep.  Here we construct a tube that changes scale by first generating the transforms and then applying the scale factor and connecting the inside and outside.  Note that the wall thickness varies because it is produced by scaling.
1736//   shape = star(n=5, r=10, ir=5);
1737//   rpath = arc(25, points=[[29,6,-4], [3,4,6], [1,1,7]]);
1738//   trans = path_sweep(shape, rpath, transforms=true);
1739//   outside = [for(i=[0:len(trans)-1]) trans[i]*scale(lerp(1,1.5,i/(len(trans)-1)))];
1740//   inside = [for(i=[len(trans)-1:-1:0]) trans[i]*scale(lerp(1.1,1.4,i/(len(trans)-1)))];
1741//   sweep(shape, concat(outside,inside),closed=true);
1742// Example(NoScales): An easier way to scale your model is to use the scale parameter.
1743//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[0,180], r=3));
1744//   path_sweep(pentagon(r=1), path3d(elliptic_arc), scale=2);
1745// Example(NoScales): Scaling only in the y direction of the profile (z direction in the model in this case)
1746//   elliptic_arc = xscale(2, p=arc($fn=64,angle=[0,180], r=3));
1747//   path_sweep(rect(2), path3d(elliptic_arc), scale=[1,2]);
1748// Example(NoScales): Specifying scale at every point for a closed path
1749//   N=64;
1750//   path = circle(r=5, $fn=64);
1751//   theta = lerpn(0,360,N,endpoint=false);
1752//   scale = [for(t=theta) sin(6*t)/5+1];
1753//   path_sweep(rect(2), path3d(path), closed=true, scale=scale);
1754// Example(Med,NoScales): Using path_sweep on a region
1755//   rgn1 = [for (d=[10:10:60]) circle(d=d,$fn=8)];
1756//   rgn2 = [square(30,center=false)];
1757//   rgn3 = [for (size=[10:10:20]) move([15,15],p=square(size=size, center=true))];
1758//   mrgn = union(rgn1,rgn2);
1759//   orgn = difference(mrgn,rgn3);
1760//   path_sweep(orgn,arc(r=40,angle=180));
1761// Example(Med,NoScales): A region with a twist
1762//   region = [for(i=pentagon(5)) move(i,p=circle(r=2,$fn=25))];
1763//   path_sweep(region,
1764//              circle(r=16,$fn=75),closed=true,
1765//              twist=360/5*2,symmetry=5);
1766// Example(Med,NoScales): Cutting a cylinder with a curved path.  Note that in this case, the incremental method produces just a slight twist but the natural method produces an extreme twist.  But manual specification produces no twist, as desired:
1767//   $fn=90;
1768//   r=8;
1769//   thickness=1;
1770//   len=21;
1771//   curve = [for(theta=[0:4:359])
1772//              [r*cos(theta), r*sin(theta), 10+sin(6*theta)]];
1773//   difference(){
1774//     cylinder(r=r, h=len);
1775//     down(.5)cylinder(r=r-thickness, h=len+1);
1776//     path_sweep(left(.05,square([1.1,1])), curve, closed=true,
1777//                method="manual", normal=UP);
1778//   }
1779// Example(Med,NoScales,VPR=[78.1,0,43.2],VPT=[2.18042,-0.485127,1.90371],VPD=74.4017): The "start" and "end" anchors are located at the origin point of the swept shape.
1780//   shape = back_half(right_half(star(n=5,id=5,od=10)),y=-1);
1781//   path = arc(angle=[0,180],d=30);
1782//   path_sweep(shape,path,method="natural"){
1783//     attach(["start","end"]) anchor_arrow(s=5);
1784//   }
1785// Example(Med,NoScales,VPR=[78.1,0,43.2],VPT=[2.18042,-0.485127,1.90371],VPD=74.4017): The "start" and "end" anchors are located at the origin point of the swept shape.
1786//   shape = back_half(right_half(star(n=5,id=5,od=10)),y=-1);
1787//   path = arc(angle=[0,180],d=30);
1788//   path_sweep(shape,path,method="natural"){
1789//     attach(["start-centroid","end-centroid"]) anchor_arrow(s=5);
1790//   }
1791// Example(Med,NoScales,VPR=[78.1,0,43.2],VPT=[2.18042,-0.485127,1.90371],VPD=74.4017): Note that the "start" anchors are backwards compared to the direction of the sweep, so you have to attach the TOP to align the shape with its ends.  
1792//   shape = back_half(right_half(star(n=5,id=5,od=10)),y=-1)[0];
1793//   path = arc(angle=[0,180],d=30);
1794//   path_sweep(shape,path,method="natural",scale=[1,1.5])
1795//     recolor("red"){
1796//       attach("start",TOP) stroke([path3d(shape)],width=.5);
1797//       attach("end") stroke([path3d(yscale(1.5,shape))],width=.5);       
1798//     }
1799
1800module path_sweep(shape, path, method="incremental", normal, closed, twist=0, twist_by_length=true, scale=1, scale_by_length=true,
1801                    symmetry=1, last_normal, tangent, uniform=true, relaxed=false, caps, style="min_edge", convexity=10,
1802                    anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull",profiles=false,width=1)
1803{
1804    dummy = assert(is_region(shape) || is_path(shape,2), "shape must be a 2D path or region")
1805            assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"");
1806    trans_scale = path_sweep(shape, path, method, normal, closed, twist, twist_by_length, scale, scale_by_length,
1807                            symmetry, last_normal, tangent, uniform, relaxed, caps, style, transforms=true,_return_scales=true);
1808    caps = is_def(caps) ? caps :
1809           closed ? false : true;
1810    fullcaps = is_bool(caps) ? [caps,caps] : caps;
1811    transforms = trans_scale[0];
1812    scales = trans_scale[1];
1813    firstscale = is_num(scales[0]) ? 1/scales[0] : [1/scales[0].x, 1/scales[0].y];
1814    lastscale = is_num(last(scales)) ? 1/last(scales) : [1/last(scales).x, 1/last(scales).y];
1815    vnf = sweep(is_path(shape)?clockwise_polygon(shape):shape, transforms, closed=false, caps=fullcaps,style=style);
1816    shapecent = point3d(centroid(shape));
1817    $sweep_transforms = transforms;
1818    $sweep_scales = scales;
1819    $sweep_shape = shape;
1820    $sweep_path = path;
1821    $sweep_closed = closed;
1822    $sweep_twist = twist;
1823    anchors = closed ? []
1824            :
1825              [
1826                named_anchor("start", rot=transforms[0]*scale(firstscale), flip=true), 
1827                named_anchor("end", rot=last(transforms)*scale(lastscale)),
1828                named_anchor("start-centroid", rot=transforms[0]*move(shapecent)*scale(firstscale), flip=true),
1829                named_anchor("end-centroid", rot=last(transforms)*move(shapecent)*scale(lastscale))
1830    ];
1831    if (profiles){
1832        rshape = is_path(shape) ? [path3d(shape)]
1833                                : [for(s=shape) path3d(s)];
1834        attachable(anchor,spin,orient, vnf=vnf, extent=atype=="hull", cp=cp, anchors=anchors) {
1835            for(T=transforms) stroke([for(part=rshape)apply(T,part)],width=width);
1836            children();
1837        }
1838    }
1839    else
1840      attachable(anchor,spin,orient,vnf=vnf,extent=atype=="hull", cp=cp,anchors=anchors){
1841          vnf_polyhedron(vnf,convexity=convexity);
1842          children();
1843      }
1844}
1845
1846
1847function path_sweep(shape, path, method="incremental", normal, closed, twist=0, twist_by_length=true, scale=1, scale_by_length=true, 
1848                    symmetry=1, last_normal, tangent, uniform=true, relaxed=false, caps, style="min_edge", transforms=false,
1849                    anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull",_return_scales=false) =
1850  is_1region(path) ? path_sweep(shape=shape,path=path[0], method=method, normal=normal, closed=default(closed,true), 
1851                                twist=twist, scale=scale, scale_by_length=scale_by_length, twist_by_length=twist_by_length, symmetry=symmetry, last_normal=last_normal,
1852                                tangent=tangent, uniform=uniform, relaxed=relaxed, caps=caps, style=style, transforms=transforms,
1853                                anchor=anchor, cp=cp, spin=spin, orient=orient, atype=atype, _return_scales=_return_scales) :
1854  let(closed=default(closed,false))
1855  assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"")
1856  assert(!closed || twist % (360/symmetry)==0, str("For a closed sweep, twist must be a multiple of 360/symmetry = ",360/symmetry))
1857  assert(closed || symmetry==1, "symmetry must be 1 when closed is false")
1858  assert(is_integer(symmetry) && symmetry>0, "symmetry must be a positive integer")
1859  let(path = force_path(path))
1860  assert(is_path(path,[2,3]), "input path is not a 2D or 3D path")
1861  assert(!closed || !approx(path[0],last(path)), "Closed path includes start point at the end")
1862  assert((is_region(shape) || is_path(shape,2)) || (transforms && !(closed && method=="incremental")),"shape must be a 2d path or region")
1863  let(
1864    path = path3d(path),
1865    caps = is_def(caps) ? caps :
1866           closed ? false : true,
1867    capsOK = is_bool(caps) || is_bool_list(caps,2),
1868    fullcaps = is_bool(caps) ? [caps,caps] : caps,
1869    normalOK = is_undef(normal) || (method!="natural" && is_vector(normal,3))
1870                                || (method=="manual" && same_shape(normal,path)),
1871    scaleOK = scale==1 || ((is_num(scale) || is_vector(scale,2)) && !closed) || is_vector(scale,len(path)) || is_matrix(scale,len(path),2)
1872    
1873  )
1874  assert(normalOK,  method=="natural" ? "Cannot specify normal with the \"natural\" method"
1875                  : method=="incremental" ? "Normal with \"incremental\" method must be a 3-vector"
1876                  : str("Incompatible normal given.  Must be a 3-vector or a list of ",len(path)," 3-vectors"))
1877  assert(capsOK, "caps must be boolean or a list of two booleans")
1878  assert(!closed || !caps, "Cannot make closed shape with caps")
1879  assert(is_undef(normal) || (is_vector(normal) && len(normal)==3) || (is_path(normal) && len(normal)==len(path) && len(normal[0])==3), "Invalid normal specified")
1880  assert(is_undef(tangent) || (is_path(tangent) && len(tangent)==len(path) && len(tangent[0])==3), "Invalid tangent specified")
1881  assert(scaleOK,str("Incompatible or invalid scale",closed?" for closed path":"",": must be ", closed?"":"a scalar, a 2-vector, ",
1882                     "a vector of length ",len(path)," or a ",len(path),"x2 matrix of scales"))
1883  let(
1884    scale = !(is_num(scale) || is_vector(scale,2)) ? scale
1885          : let(s=is_num(scale) ? [scale,scale] : scale)
1886            !scale_by_length ? lerpn([1,1],s,len(path))
1887          : lerp([1,1],s, path_length_fractions(path,false)),
1888    scale_list = [for(s=scale) scale(s),if (closed) scale(scale[0])],
1889    tangents = is_undef(tangent) ? path_tangents(path,uniform=uniform,closed=closed) : [for(t=tangent) unit(t)],
1890    normal = is_path(normal) ? [for(n=normal) unit(n)] :
1891             is_def(normal) ? unit(normal) :
1892             method =="incremental" && abs(tangents[0].z) > 1/sqrt(2) ? BACK : UP,
1893    normals = is_path(normal) ? normal : repeat(normal,len(path)),
1894    tpathfrac = twist_by_length ? path_length_fractions(path, closed) : [for(i=[0:1:len(path)]) i / (len(path)-(closed?0:1))],
1895    spathfrac = scale_by_length ? path_length_fractions(path, closed) : [for(i=[0:1:len(path)]) i / (len(path)-(closed?0:1))],    
1896    L = len(path),
1897    unscaled_transform_list =
1898        method=="old_incremental" ?
1899          let(rotations =
1900                 [for( i  = 0,
1901                       ynormal = normal - (normal * tangents[0])*tangents[0],
1902                       rotation = frame_map(y=ynormal, z=tangents[0])
1903                         ;
1904                       i < len(tangents) + (closed?1:0)
1905                         ;
1906                       rotation = i<len(tangents)-1+(closed?1:0)? rot(from=tangents[i],to=tangents[(i+1)%L])*rotation : undef,
1907                       i=i+1
1908                      )
1909                   rotation],
1910              // The mismatch is the inverse of the last transform times the first one for the closed case, or the inverse of the
1911              // desired final transform times the realized final transform in the open case.  Note that when closed==true the last transform
1912              // is a actually looped around and applies to the first point position, so if we got back exactly where we started
1913              // then it will be the identity, but we might have accumulated some twist which will show up as a rotation around the
1914              // X axis.  Similarly, in the closed==false case the desired and actual transformations can only differ in the twist,
1915              // so we can need to calculate the twist angle so we can apply a correction, which we distribute uniformly over the whole path.
1916              reference_rot = closed ? rotations[0] :
1917                           is_undef(last_normal) ? last(rotations) :
1918                             let(
1919                                 last_tangent = last(tangents),
1920                                 lastynormal = last_normal - (last_normal * last_tangent) * last_tangent
1921                             )
1922                           frame_map(y=lastynormal, z=last_tangent),
1923              mismatch = transpose(last(rotations)) * reference_rot,
1924              correction_twist = atan2(mismatch[1][0], mismatch[0][0]),
1925              // Spread out this extra twist over the whole sweep so that it doesn't occur
1926              // abruptly as an artifact at the last step.
1927              twistfix = correction_twist%(360/symmetry),
1928              adjusted_final = !closed ? undef :
1929                            translate(path[0]) * rotations[0] * zrot(-correction_twist+correction_twist%(360/symmetry)-twist)
1930          )  [for(i=idx(path)) translate(path[i]) * rotations[i] * zrot((twistfix-twist)*tpathfrac[i]), if(closed) adjusted_final] 
1931      : method=="incremental" ?   // Implements Rotation Minimizing Frame from "Computation of Rotation Minimizing Frames"
1932                                  // by Wenping Yang, Bert Büttler, Dayue Zheng, Yang Liu, 2008
1933                                  // http://doi.acm.org/10.1145/1330511.1330513
1934          let(rotations =         // https://www.microsoft.com/en-us/research/wp-content/uploads/2016/12/Computation-of-rotation-minimizing-frames.pdf
1935                 [for( i  = 0,
1936                       ynormal = normal - (normal * tangents[0])*tangents[0],
1937                       rotation = frame_map(y=ynormal, z=tangents[0]),
1938                       r=ynormal
1939                         ;
1940                       i < len(tangents) + (closed?1:0)
1941                         ;
1942                       v1 = path[(i+1)%L]-path[i%L],
1943                       c1 = v1*v1,
1944                       rL = r - 2*(v1*r)/c1 * v1,
1945                       tL = tangents[i%L] - 2*(v1*tangents[i%L])/c1 * v1,
1946                       v2 = tangents[(i+1)%L]-tL,
1947                       c2 = v2*v2,
1948                       r = rL - (2/c2)*(v2*rL)*v2,
1949                       rotation = i<len(tangents)-1+(closed?1:0)? frame_map(y=r,z=tangents[(i+1)%L]) : undef,
1950                       i=i+1
1951                      )
1952                   rotation],
1953              // The mismatch is the inverse of the last transform times the first one for the closed case, or the inverse of the
1954              // desired final transform times the realized final transform in the open case.  Note that when closed==true the last transform
1955              // is a actually looped around and applies to the first point position, so if we got back exactly where we started
1956              // then it will be the identity, but we might have accumulated some twist which will show up as a rotation around the
1957              // X axis.  Similarly, in the closed==false case the desired and actual transformations can only differ in the twist,
1958              // so we can need to calculate the twist angle so we can apply a correction, which we distribute uniformly over the whole path.
1959              reference_rot = closed ? rotations[0] :
1960                           is_undef(last_normal) ? last(rotations) :
1961                             let(
1962                                 last_tangent = last(tangents),
1963                                 lastynormal = last_normal - (last_normal * last_tangent) * last_tangent
1964                             )
1965                           frame_map(y=lastynormal, z=last_tangent),
1966              mismatch = transpose(last(rotations)) * reference_rot,
1967              correction_twist = atan2(mismatch[1][0], mismatch[0][0]),
1968              // Spread out this extra twist over the whole sweep so that it doesn't occur
1969              // abruptly as an artifact at the last step.
1970              twistfix = correction_twist%(360/symmetry),
1971              adjusted_final = !closed ? undef :
1972                            translate(path[0]) * rotations[0] * zrot(-correction_twist+correction_twist%(360/symmetry)-twist)
1973          )  [for(i=idx(path)) translate(path[i]) * rotations[i] * zrot((twistfix-twist)*tpathfrac[i]), if(closed) adjusted_final] 
1974      : method=="manual" ?
1975              [for(i=[0:L-(closed?0:1)]) let(
1976                       ynormal = relaxed ? normals[i%L] : normals[i%L] - (normals[i%L] * tangents[i%L])*tangents[i%L],
1977                       znormal = relaxed ? tangents[i%L] - (normals[i%L] * tangents[i%L])*normals[i%L] : tangents[i%L],
1978                       rotation = frame_map(y=ynormal, z=znormal)
1979                   )
1980                   assert(approx(ynormal*znormal,0),str("Supplied normal is parallel to the path tangent at point ",i))
1981                   translate(path[i%L])*rotation*zrot(-twist*tpathfrac[i])
1982              ]
1983      : method=="natural" ?   // map x axis of shape to the path normal, which points in direction of curvature
1984              let (pathnormal = path_normals(path, tangents, closed))
1985              assert(all_defined(pathnormal),"Natural normal vanishes on your curve, select a different method")
1986              let( testnormals = [for(i=[0:len(pathnormal)-1-(closed?1:2)]) pathnormal[i]*select(pathnormal,i+2)],
1987                   a=[for(i=idx(testnormals)) testnormals[i]<.5 ? echo(str("Big change at index ",i," pn=",pathnormal[i]," pn2= ",select(pathnormal,i+2))):0],
1988                   dummy = min(testnormals) < .5 ? echo("WARNING: ***** Abrupt change in normal direction.  Consider a different method in path_sweep() *****") :0
1989                 )
1990              [for(i=[0:L-(closed?0:1)]) let(
1991                       rotation = frame_map(x=pathnormal[i%L], z=tangents[i%L])
1992                   )
1993                   translate(path[i%L])*rotation*zrot(-twist*tpathfrac[i])
1994                 ] 
1995      : assert(false,"Unknown method or no method given"), // unknown method
1996    transform_list = v_mul(unscaled_transform_list, scale_list),
1997    ends_match = !closed ? true
1998                 : let( rshape = is_path(shape) ? [path3d(shape)]
1999                                                : [for(s=shape) path3d(s)]
2000                   )
2001                   are_regions_equal(apply(transform_list[0], rshape),
2002                                     apply(transform_list[L], rshape)),
2003    dummy = ends_match ? 0 : echo("WARNING: ***** The points do not match when closing the model in path_sweep() *****")
2004  )
2005  transforms && _return_scales
2006             ? [transform_list,scale]
2007: transforms ? transform_list
2008             : sweep(is_path(shape)?clockwise_polygon(shape):shape, transform_list, closed=false, caps=fullcaps,style=style,
2009                       anchor=anchor,cp=cp,spin=spin,orient=orient,atype=atype);
2010
2011
2012// Function&Module: path_sweep2d()
2013// Synopsis: Sweep a 2d polygon path along a 2d path allowing self-intersection. 
2014// SynTags: VNF, Geom
2015// Topics: Extrusion, Sweep, Paths
2016// See Also: linear_sweep(), rotate_sweep(), sweep(), spiral_sweep(), path_sweep(), offset_sweep()
2017// Usage: as module
2018//   path_sweep2d(shape, path, [closed], [caps], [quality], [style], [convexity=], [anchor=], [spin=], [orient=], [atype=], [cp=]) [ATTACHMENTS];
2019// Usage: as function
2020//   vnf = path_sweep2d(shape, path, [closed], [caps], [quality], [style], [anchor=], [spin=], [orient=], [atype=], [cp=]);
2021// Description:
2022//   Takes an input 2D polygon (the shape) and a 2d path, and constructs a polyhedron by sweeping the shape along the path.
2023//   When run as a module returns the polyhedron geometry.  When run as a function returns a VNF.
2024//   .
2025//   See {{path_sweep()}} for more details on how the sweep operation works and for introductory examples.
2026//   This 2d version is different because local self-intersections (creases in the output) are allowed and do not produce CGAL errors.
2027//   This is accomplished by using offset() calculations, which are more expensive than simply copying the shape along
2028//   the path, so if you do not have local self-intersections, use {{path_sweep()}} instead.  If xmax is the largest x value (in absolute value)
2029//   of the shape, then path_sweep2d() will work as long as the offset of `path` exists at `delta=xmax`.  If the offset vanishes, as in the
2030//   case of a circle offset by more than its radius, then you will get an error about a degenerate offset.
2031//   Note that global self-intersections will still give rise to CGAL errors.  You should be able to handle these by partitioning your model.  The y axis of the
2032//   shape is mapped to the z axis in the swept polyhedron, and no twisting can occur.
2033//   The quality parameter is passed to offset to determine the offset quality.
2034// Arguments:
2035//   shape = a 2D polygon describing the shape to be swept
2036//   path = a 2D path giving the path to sweep over
2037//   closed = path is a closed loop.  Default: false
2038//   caps = true to create endcap faces when closed is false.  Can be a length 2 boolean array.  Default is true if closed is false.
2039//   quality = quality of offset used in calculation.  Default: 1
2040//   style = vnf_vertex_array style.  Default: "min_edge"
2041//   ---
2042//   convexity = convexity parameter for polyhedron (module only)  Default: 10
2043//   anchor = Translate so anchor point is at the origin.  Default: "origin"
2044//   spin = Rotate this many degrees around Z axis after anchor.  Default: 0
2045//   orient = Vector to rotate top towards after spin
2046//   atype = Select "hull" or "intersect" anchor types.  Default: "hull"
2047//   cp = Centerpoint for determining "intersect" anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
2048// Named Anchors:
2049//   "origin" = The native position of the shape.  
2050// Anchor Types:
2051//   "hull" = Anchors to the virtual convex hull of the shape.
2052//   "intersect" = Anchors to the surface of the shape.
2053// Example: Sine wave example with self-intersections at each peak.  This would fail with path_sweep().
2054//   sinewave = [for(i=[-30:10:360*2+30]) [i/40,3*sin(i)]];
2055//   path_sweep2d(circle(r=3,$fn=15), sinewave);
2056// Example: The ends can look weird if they are in a place where self intersection occurs.  This is a natural result of how offset behaves at ends of a path.
2057//   coswave = [for(i=[0:10:360*1.5]) [i/40,3*cos(i)]];
2058//   zrot(-20)
2059//     path_sweep2d( circle(r=3,$fn=15), coswave);
2060// Example: This closed path example works ok as long as the hole in the center remains open.
2061//   ellipse = yscale(3,p=circle(r=3,$fn=120));
2062//   path_sweep2d(circle(r=2.5,$fn=32), reverse(ellipse), closed=true);
2063// Example: When the hole is closed a global intersection renders the model invalid.  You can fix this by taking the union of the two (valid) halves.
2064//   ellipse = yscale(3,p=circle(r=3,$fn=120));
2065//   L = len(ellipse);
2066//   path_sweep2d(circle(r=3.25, $fn=32), select(ellipse,floor(L*.2),ceil(L*.8)),closed=false);
2067//   path_sweep2d(circle(r=3.25, $fn=32), select(ellipse,floor(L*.7),ceil(L*.3)),closed=false);
2068
2069function path_sweep2d(shape, path, closed=false, caps, quality=1, style="min_edge",
2070                      anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull") =
2071   let(
2072        caps = is_def(caps) ? caps
2073             : closed ? false : true,
2074        capsOK = is_bool(caps) || is_bool_list(caps,2),
2075        fullcaps = is_bool(caps) ? [caps,caps] : caps,
2076        shape = force_path(shape,"shape"),
2077        path = force_path(path)
2078   )
2079   assert(is_path(shape,2), "shape must be a 2D path")
2080   assert(is_path(path,2), "path must be a 2D path")
2081   assert(capsOK, "caps must be boolean or a list of two booleans")
2082   assert(!closed || !caps, "Cannot make closed shape with caps")
2083   let(
2084        profile = ccw_polygon(shape),
2085        flip = closed && is_polygon_clockwise(path) ? -1 : 1,
2086        path = flip ? reverse(path) : path,
2087        proflist= transpose(
2088                     [for(pt = profile)
2089                        let(
2090                            ofs = offset(path, delta=-flip*pt.x, return_faces=true,closed=closed, quality=quality),
2091                            map = column(_ofs_vmap(ofs,closed=closed),1)
2092                        )
2093                        select(path3d(ofs[0],pt.y),map)
2094                      ]
2095                  ),
2096        vnf = vnf_vertex_array([
2097                         each proflist,
2098                         if (closed) proflist[0]
2099                        ],cap1=fullcaps[0],cap2=fullcaps[1],col_wrap=true,style=style)
2100   )
2101   reorient(anchor,spin,orient,vnf=vnf,p=vnf,extent=atype=="hull",cp=cp);
2102
2103
2104module path_sweep2d(profile, path, closed=false, caps, quality=1, style="min_edge", convexity=10,
2105                    anchor="origin", cp="centroid", spin=0, orient=UP, atype="hull")
2106{
2107   vnf = path_sweep2d(profile, path, closed, caps, quality, style);
2108   vnf_polyhedron(vnf,convexity=convexity,anchor=anchor, spin=spin, orient=orient, atype=atype, cp=cp)
2109        children();
2110}
2111
2112// Extract vertex mapping from offset face list.  The output of this function
2113// is a list of pairs [i,j] where i is an index into the parent curve and j is
2114// an index into the offset curve.  It would probably make sense to rewrite
2115// offset() to return this instead of the face list and have offset_sweep
2116// use this input to assemble the faces it needs.
2117
2118function _ofs_vmap(ofs,closed=false) =
2119    let(   // Caclulate length of the first (parent) curve
2120        firstlen = max(flatten(ofs[1]))+1-len(ofs[0])
2121    )
2122    [
2123     for(entry=ofs[1]) _ofs_face_edge(entry,firstlen),
2124     if (!closed) _ofs_face_edge(last(ofs[1]),firstlen,second=true)
2125    ];
2126
2127
2128// Extract first (default) or second edge that connects the parent curve to its offset.  The first input
2129// face is a list of 3 or 4 vertices as indices into the two curves where the parent curve vertices are
2130// numbered from 0 to firstlen-1 and the offset from firstlen and up.  The firstlen pararameter is used
2131// to determine which curve the vertices belong to and to remove the offset so that the return gives
2132// the index into each curve with a 0 base.
2133function _ofs_face_edge(face,firstlen,second=false) =
2134   let(
2135       itry = min_index(face),
2136       i = select(face,itry-1)<firstlen ? itry-1:itry,
2137       edge1 = select(face,[i,i-1]),
2138       edge2 = select(face,i+1)<firstlen ? select(face,[i+1,i+2])
2139                                         : select(face,[i,i+1])
2140   )
2141   (second ? edge2 : edge1)-[0,firstlen];
2142
2143
2144
2145// Function&Module: sweep()
2146// Synopsis: Construct a 3d object from arbitrary transformations of a 2d polygon path.
2147// SynTags: VNF, Geom
2148// Topics: Extrusion, Sweep, Paths
2149// See Also: sweep_attach(), linear_sweep(), rotate_sweep(), spiral_sweep(), path_sweep(), path_sweep2d(), offset_sweep()
2150// Usage: As Module
2151//   sweep(shape, transforms, [closed], [caps], [style], [convexity=], [anchor=], [spin=], [orient=], [atype=]) [ATTACHMENTS];
2152// Usage: As Function
2153//   vnf = sweep(shape, transforms, [closed], [caps], [style], [anchor=], [spin=], [orient=], [atype=]);
2154// Description:
2155//   The input `shape` must be a non-self-intersecting 2D polygon or region, and `transforms`
2156//   is a list of 4x4 transformation matrices.  The sweep algorithm applies each transformation in sequence
2157//   to the shape input and links the resulting polygons together to form a polyhedron.
2158//   If `closed=true` then the first and last transformation are linked together.
2159//   The `caps` parameter controls whether the ends of the shape are closed.
2160//   As a function, returns the VNF for the polyhedron.  As a module, computes the polyhedron.
2161//   .
2162//   Note that this is a very powerful, general framework for producing polyhedra.  It is important
2163//   to ensure that your resulting polyhedron does not include any self-intersections, or it will
2164//   be invalid and will generate CGAL errors.  If you get such errors, most likely you have an
2165//   overlooked self-intersection.  Note also that the errors will not occur when your shape is alone
2166//   in your model, but will arise if you add a second object to the model.  This may mislead you into
2167//   thinking the second object caused a problem.  Even adding a simple cube to the model will reveal the problem.
2168// Arguments:
2169//   shape = 2d path or region, describing the shape to be swept.
2170//   transforms = list of 4x4 matrices to apply
2171//   closed = set to true to form a closed (torus) model.  Default: false
2172//   caps = true to create endcap faces when closed is false.  Can be a singe boolean to specify endcaps at both ends, or a length 2 boolean array.  Default is true if closed is false.
2173//   style = vnf_vertex_array style.  Default: "min_edge"
2174//   ---
2175//   convexity = convexity setting for use with polyhedron.  (module only) Default: 10
2176//   cp = Centerpoint for determining "intersect" anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
2177//   atype = Select "hull" or "intersect" anchor types.  Default: "hull"
2178//   anchor = Translate so anchor point is at the origin. Default: "origin"
2179//   spin = Rotate this many degrees around Z axis after anchor. Default: 0
2180//   orient = Vector to rotate top towards after spin  (module only)
2181// Named Anchors:
2182//   "origin" = The native position of the shape.  
2183// Anchor Types:
2184//   "hull" = Anchors to the virtual convex hull of the shape.
2185//   "intersect" = Anchors to the surface of the shape.
2186// Example(VPR=[45,0,74],VPD=175,VPT=[-3.8,12.4,19]): A bent object that also changes shape along its length.
2187//   radius = 75;
2188//   angle = 40;
2189//   shape = circle(r=5,$fn=32);
2190//   T = [for(i=[0:25]) xrot(-angle*i/25,cp=[0,radius,0])*scale([1+i/25, 2-i/25,1])];
2191//   sweep(shape,T);
2192// Example: This is the "sweep-drop" example from list-comprehension-demos.
2193//   function drop(t) = 100 * 0.5 * (1 - cos(180 * t)) * sin(180 * t) + 1;
2194//   function path(t) = [0, 0, 80 + 80 * cos(180 * t)];
2195//   function rotate(t) = 180 * pow((1 - t), 3);
2196//   step = 0.01;
2197//   path_transforms = [for (t=[0:step:1-step]) translate(path(t)) * zrot(rotate(t)) * scale([drop(t), drop(t), 1])];
2198//   sweep(circle(1, $fn=12), path_transforms);
2199// Example: Another example from list-comprehension-demos
2200//   function f(x) = 3 - 2.5 * x;
2201//   function r(x) = 2 * 180 * x * x * x;
2202//   pathstep = 1;
2203//   height = 100;
2204//   shape_points = subdivide_path(square(10),40,closed=true);
2205//   path_transforms = [for (i=[0:pathstep:height]) let(t=i/height) up(i) * scale([f(t),f(t),i]) * zrot(r(t))];
2206//   sweep(shape_points, path_transforms);
2207// Example: Twisted container.  Note that this technique doesn't create a fixed container wall thickness.
2208//   shape = subdivide_path(square(30,center=true), 40, closed=true);
2209//   outside = [for(i=[0:24]) up(i)*rot(i)*scale(1.25*i/24+1)];
2210//   inside = [for(i=[24:-1:2]) up(i)*rot(i)*scale(1.2*i/24+1)];
2211//   sweep(shape, concat(outside,inside));
2212
2213function sweep(shape, transforms, closed=false, caps, style="min_edge",
2214               anchor="origin", cp="centroid", spin=0, orient=UP, atype="hull") =
2215    assert(is_consistent(transforms, ident(4)), "Input transforms must be a list of numeric 4x4 matrices in sweep")
2216    assert(is_path(shape,2) || is_region(shape), "Input shape must be a 2d path or a region.")
2217    let(
2218        caps = is_def(caps) ? caps :
2219            closed ? false : true,
2220        capsOK = is_bool(caps) || is_bool_list(caps,2),
2221        fullcaps = is_bool(caps) ? [caps,caps] : caps
2222    )
2223    assert(len(transforms)>=2, "transformation must be length 2 or more")
2224    assert(capsOK, "caps must be boolean or a list of two booleans")
2225    assert(!closed || !caps, "Cannot make closed shape with caps")
2226    is_region(shape)? let(
2227        regions = region_parts(shape),
2228        rtrans = reverse(transforms),
2229        vnfs = [
2230            for (rgn=regions) each [
2231                for (path=rgn)
2232                    sweep(path, transforms, closed=closed, caps=false, style=style),
2233                if (fullcaps[0]) vnf_from_region(rgn, transform=transforms[0], reverse=true),
2234                if (fullcaps[1]) vnf_from_region(rgn, transform=last(transforms)),
2235            ],
2236        ],
2237        vnf = vnf_join(vnfs)
2238    ) vnf :
2239    assert(len(shape)>=3, "shape must be a path of at least 3 non-colinear points")
2240    vnf_vertex_array([for(i=[0:len(transforms)-(closed?0:1)]) apply(transforms[i%len(transforms)],path3d(shape))],
2241                     cap1=fullcaps[0],cap2=fullcaps[1],col_wrap=true,style=style);
2242
2243
2244module sweep(shape, transforms, closed=false, caps, style="min_edge", convexity=10,
2245             anchor="origin",cp="centroid",spin=0, orient=UP, atype="hull")
2246{
2247    $sweep_transforms=transforms;
2248    $sweep_shape=shape;
2249    $sweep_closed=closed;
2250    vnf = sweep(shape, transforms, closed, caps, style);
2251    vnf_polyhedron(vnf, convexity=convexity, anchor=anchor, spin=spin, orient=orient, atype=atype, cp=cp)
2252        children();
2253}
2254
2255
2256
2257// Section: Attaching children to sweeps
2258
2259
2260// Module: sweep_attach()
2261// Synopsis: Attach children to sides of a path_sweep parent object
2262// SynTags: Geom
2263// Topics: Extrusion, Sweep, Paths
2264// See Also: path_sweep()
2265// Usage:
2266//   path_sweep(...) { sweep_attach(parent, [child], [frac], [idx=], [len=], [spin=], [overlap=], [atype=]) CHILDREN; }
2267//   sweep(...) { sweep_attach(parent, [child], [frac], [idx=], [len=], [spin=], [overlap=], [atype=]) CHILDREN; }
2268// Description:
2269//   Attaches children to the sides of a {{path_sweep()}} or {{sweep()}} object.  You supply a position along the path,
2270//   either by path fraction, length, or index.  In the case of `sweep()` objects the path is defined as the path traced out
2271//   by the origin of the shape under the transformation list.  Objects are attached with their UP direction aligned with
2272//   the anchor for the profile and their BACK direction pointing in the direction of the sweep.
2273//   .
2274//   Like {{attach()}} this module has a parent-child anchor mode where you specify the child anchor and it is
2275//   aligned with the anchor on the sweep.  As with {{attach()}}, the child `anchor` and `orient` parameters are ignored.
2276//   Alternative you can use parent anchor mode where give only the parent anchor and the child appears at its
2277//   child-specified (default) anchor point.  The spin parameter spins the child around the attachment anchor axis.  
2278//   .
2279//   For a path_sweep() with no scaling, if you give a location or index that is exactly at one of the sections the normal will be in the plane
2280//   of the section.  In the general case if you give a location in between sections the normal will be normal to the facet.  If you
2281//   give a location at a section in the general case the normal will be the average of the normals of the two adjacent facets.  
2282//   For twisted or other complicated sweeps the normals may not be accurate.  If you need accurate normals for such shapes, you must
2283//   use the anchors for the VNF swept shape directly---it is a tradeoff between easy specification of the anchor location on the
2284//   swept object, which may be very difficult with direct anchors, and accuracy of the normal.
2285//   .
2286//   For closed sweeps the index will wrap around and can be positive or negative.  For sweeps that are not closed the index must
2287//   be positive and no longer than the length of the path.  In some cases for closed path_sweeps the shape can be a mobius strip
2288//   and it may take more than one cycle to return to the starting point.  The extra twist will be properly handled in this case.
2289//   If you construct a mobius strip using the generic {{sweep()}} then information about the amount of twist is not available
2290//   to `sweep_attach()` so it will not be handled automatically.  
2291//   .
2292//   The anchor you give acts as a 2D anchor to the path or region used by the sweep, in the XY plane as that shape appears
2293//   before it is transformed to form the swept object.  As with {{region()}}, you can control the anchor using `cp` and `atype`, 
2294//   and you can check the anchors by using the same anchors with {{region()}} in a two dimensional test case.
2295//   .
2296//   Note that {{path_sweep2d()}} does not support `sweep_attach()` because it doesn't compute the transform list, which is
2297//   the input used to calculate the attachment transform.  
2298// Arguments:
2299//   anchor = 2d anchor to the shape used in the path_sweep parent
2300//   frac = position along the path_sweep path as a fraction of total length
2301//   ---
2302//   idx = index into the path_sweep path (use instead of frac)
2303//   len = absolute length along the path_sweep path (use instead of frac)
2304//   spin = spin the child this amount around the anchor axis.  Default: 0
2305//   overlap = Amount to lower the shape into the parent.  Default: 0
2306//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 2D point.  Default: "centroid"
2307//   atype = Set to "hull" or "intersect" to select anchor type.  Default: "hull"
2308// Anchor Types:
2309//   "hull" = Anchors to the virtual convex hull of the region.
2310//   "intersect" = Anchors to the outer edge of the region.
2311// Example(Med,NoAxes,VPT=[4.75027,0.805639,-3.50893],VPR=[66.9,0,219.6],VPD=213.382): This example shows several children positioned at different places on the parent.  The blue cone is positioned using its TOP anchor and is sunk into the parent with overlay.  The three orange cubes show how they rotate to follow the local sweep direction.  
2312//   function a(h) = arc(points=[[-20,0],[0,h],[20,0]],n=24);
2313//   shape = concat(
2314//                  a(2), // bottom
2315//                  back(6,reverse(a(4))) // top
2316//           );
2317//   path = xrot(90,path3d(arc(points=[[-40,0],[0,5],[40,-20]],n=36)));
2318//   path_sweep(shape,path) {
2319//       sweep_attach(BACK,BOT,0.2) recolor("red") cyl(d1=5,d2=0,h=8,$fn=12);
2320//       sweep_attach(BACK,TOP,0.5,overlap=3) recolor("blue") cyl(d1=5,d2=0,h=8,$fn=12);
2321//       sweep_attach(RIGHT,BOT,idx=15) recolor("orange") cuboid([3,3,5]);
2322//       sweep_attach(RIGHT,BOT,idx=1) recolor("orange") cuboid([3,3,5]);
2323//       sweep_attach(RIGHT,BOT,idx=32) recolor("orange") cuboid([3,3,5]);        
2324//   }
2325// Example(VPT=[20.7561,8.89872,0.901718],VPR=[32.6,0,338.8],VPD=66.9616,NoAxes): In this example with scaling the objects' normals are not in the plane of the path_sweep sections.  
2326//   shape = hexagon(r=4);
2327//   path = xscale(2,arc(r=15, angle=[0,75],n=10));
2328//   path_sweep(shape,path,scale=3)
2329//   {  
2330//      sweep_attach(RIGHT,BOT,0)
2331//         color_this("red")cuboid([1,1,4]);
2332//      sweep_attach(RIGHT,BOT,0.5)
2333//         color_this("blue")cuboid([1,1,4]);
2334//      sweep_attach(BACK,BOT,1/3)
2335//         color_this("lightblue")prismoid(3,1,3);
2336//   }
2337// Example(Med): This pentagonal torus is a mobius strip.  It takes five times around to return to your starting point.  Here the red box has gone 4.4 times around.  
2338//   ellipse = xscale(2, p=circle($fn=64, r=3));
2339//   pentagon = subdivide_path(pentagon(r=1), 30);
2340//   path_sweep(pentagon, path3d(ellipse),
2341//              closed=true, twist=360*2/5,symmetry=5)
2342//     sweep_attach(RIGHT,BOT,4.4) color("red") cuboid([.25,.25,3]);
2343// Example(VPT=[17.1585,9.05454,50.69],VPR=[67.6,0,64.9],VPD=292.705,NoAxes): Example using {{sweep()}}
2344//   function f(x) = 3 - 2.5 * x;
2345//   function r(x) = 2 * 180 * x * x * x;
2346//   pathstep = 1;
2347//   height = 100;
2348//   shape_points = subdivide_path(square(10),40,closed=true);
2349//   path_transforms = [for (i=[0:pathstep:height]) let(t=i/height) up(i) * scale([f(t),f(t),i]) * zrot(r(t))];
2350//   sweep(shape_points, path_transforms){
2351//     sweep_attach(RIGHT,BOT,idx=33)
2352//           color_this("red")cuboid([5,5,5]);
2353//     sweep_attach(FWD,BOT,idx=65)
2354//           color_this("red")cuboid([5,5,5]);
2355//   }
2356
2357module sweep_attach(parent, child, frac, idx, pathlen, spin=0, overlap=0, atype="hull", cp="centroid")
2358{
2359   $attach_to=child;
2360   req_children($children);
2361   dummy =  assert(!is_undef($sweep_transforms), "sweep_attach() must be used as a child of sweep() or path_sweep()")
2362            assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"")
2363            assert(num_defined([idx,frac,pathlen])==1, "Must define exactly one of idx, frac and pathlen")
2364            assert(is_undef(idx) || is_finite(idx), "idx must be a number")
2365            assert(is_undef(frac) || is_finite(frac), "frac must be a number");
2366   parmset = is_def(frac) ? "frac"
2367           : is_def(pathlen) ? "pathlen"
2368           : "idx";
2369   path = !is_undef($sweep_path) ? $sweep_path
2370        : [for(T=$sweep_transforms) apply(T,CTR)];
2371   seglen = path_segment_lengths(path,closed=$sweep_closed);
2372   pathcum = [0, each cumsum(seglen)];
2373   totlen = last(pathcum);
2374   pathtable = [for(i=idx(pathcum)) [pathcum[i],i]];
2375   i = _force_int(is_def(idx) ? idx
2376                :let(
2377                      pathlen = is_def(pathlen) ? pathlen : frac*totlen
2378                  )
2379                  lookup(posmod(pathlen,totlen),pathtable)+len($sweep_transforms)*floor(pathlen/totlen) //floor(abs(pathlen)/totlen)*sign(pathlen)
2380   );
2381   twist = is_undef($sweep_twist) ? ident(4)
2382         : let(
2383                L = len($sweep_transforms),
2384                absturn = floor(abs(i)/L),
2385                turns = floor(i/L) //sign(i)*absturn-1
2386           )
2387           zrot(-turns*$sweep_twist);
2388   geom = attach_geom(region=force_region($sweep_shape), two_d=true, extent=atype=="hull", cp=cp);
2389   anchor_data = _find_anchor(parent, geom);
2390   anchor_pos = point3d(anchor_data[1]);
2391   anchor_dir = point3d(anchor_data[2]);
2392   length = len($sweep_transforms);
2393   nextind = is_int(i) ? i>=length-1 && !$sweep_closed ? assert(i==length-1,str(parmset," is too large for the path")) undef
2394                       : i+1
2395          : $sweep_closed ?  posmod(ceil(i),length)
2396          : assert(i<length-1,str(parmset," is too large for the path")) ceil(i);
2397   prevind = is_int(i) ? i<=0 && !$sweep_closed ? assert(i==0,str(parmset," must be nonnegative")) undef
2398                       : i-1 
2399           : $sweep_closed ? floor(i)
2400           : assert(i>0,str(parmset, " must be nonnegative")) floor(i);
2401   uniform = is_undef($sweep_scales) ? false
2402           : let( 
2403                   slist = [if (is_def(prevind)) select($sweep_scales,prevind),
2404                            select($sweep_scales,i),
2405                            if (is_def(nextind)) select($sweep_scales,nextind)]
2406             )
2407             all_equal(slist);
2408   if (is_int(i) && uniform){      // Unscaled integer case: just use the profile transformation
2409       multmatrix(select($sweep_transforms,i)*twist)
2410         translate(anchor_pos)
2411         yrot(spin)
2412           frame_map(z=point3d(anchor_dir),y=UP) down(overlap) children();
2413   }
2414   else if (is_int(i) && all_defined([nextind,prevind])) {      // Scaled integer case, must average two adjacent facets
2415       frac1 = 0.1*min(seglen[i-1],seglen[i])/seglen[i-1];   // But can't average two facets at ends so exclude that case    
2416       frac2 = 0.1*min(seglen[i-1],seglen[i])/seglen[i];       
2417       dirsprev = _find_ps_dir(frac1,prevind,i,twist,anchor_pos,anchor_dir); 
2418       dirsnext = _find_ps_dir(frac2,i,nextind,twist,anchor_pos,anchor_dir);
2419       pos = apply($sweep_transforms[i]*twist, anchor_pos);
2420       mixdir = dirsprev[2]+dirsnext[2];   // Normal direction
2421       ydir=cross(cross(mixdir, dirsprev[1]+dirsnext[1]),mixdir);  // y direction perpendicular to mixdir
2422       translate(pos)
2423         rotate(v=mixdir,a=spin)
2424         frame_map(y=ydir, z=mixdir)
2425           down(overlap)
2426           children();
2427  }
2428  else {                       // Non-integer case or scaled integer at the ends: compute directions from single facet
2429    interp = is_undef(prevind)?0
2430           : is_undef(nextind)?1
2431           : i-floor(i);
2432    dirs = _find_ps_dir(interp,first_defined([prevind,i]),first_defined([nextind,i]),twist,anchor_pos,anchor_dir);
2433    translate(dirs[0])
2434        rotate(v=dirs[2],a=spin)
2435        frame_map(y=dirs[1], z=dirs[2])
2436        down(overlap) children();
2437  }
2438}     
2439
2440function _force_int(x) = approx(round(x),x) ? round(x) : x;
2441
2442// This function finds the normal to a facet on the path sweep
2443// prevind and nextind are the indices into the path, frac is the
2444// interpolation value bewteen them.
2445// anchor_pos and anchor_dir are the anchor data for the 2d shape
2446// Return is [position, ydirection, zdirection], where zdirection
2447// is normal to the facet.  Note that frac is only needed because
2448// of the possibility of twist.  
2449
2450function _find_ps_dir(frac,prevind,nextind,twist,anchor_pos,anchor_dir) =
2451  let(
2452      length = len($sweep_transforms),
2453      prevpos = apply(select($sweep_transforms,prevind)*twist,anchor_pos),
2454      nextpos = apply(select($sweep_transforms,nextind)*twist,anchor_pos),
2455      curpos = lerp(prevpos,nextpos,frac),
2456
2457      prevposdir = apply(select($sweep_transforms,prevind)*twist,anchor_pos+anchor_dir),
2458      nextposdir = apply(select($sweep_transforms,nextind)*twist,anchor_pos+anchor_dir),
2459      curposdir = lerp(prevposdir, nextposdir, frac),
2460      dir = curposdir-curpos,
2461      
2462      normal_plane = plane_from_normal(nextpos-prevpos,curpos),
2463      other_plane = plane3pt(nextpos, prevpos, curposdir),
2464      normal=plane_intersection(normal_plane, other_plane),
2465      ndir = unit(normal[1]-normal[0]),
2466      flip = sign(ndir*dir)
2467  )
2468  [curpos, nextpos-prevpos, flip*ndir];
2469
2470
2471
2472
2473
2474
2475// Section: Functions for resampling and slicing profile lists
2476
2477// Function: subdivide_and_slice()
2478// Synopsis: Resample list of paths to have the same point count and interpolate additional paths. 
2479// SynTags: PathList
2480// Topics: Paths, Path Subdivision
2481// See Also: slice_profiles()
2482// Usage:
2483//   newprof = subdivide_and_slice(profiles, slices, [numpoints], [method], [closed]);
2484// Description:
2485//   Subdivides the input profiles to have length `numpoints` where `numpoints` must be at least as
2486//   big as the largest input profile.  By default `numpoints` is set equal to the length of the
2487//   largest profile.  You can set `numpoints="lcm"` to sample to the least common multiple of all
2488//   curves, which will avoid sampling artifacts but may produce a huge output.  After subdivision,
2489//   profiles are sliced.
2490// Arguments:
2491//   profiles = profiles to operate on
2492//   slices = number of slices to insert between each pair of profiles.  May be a vector
2493//   numpoints = number of points after sampling.
2494//   method = method used for calling {{subdivide_path()}}, either `"length"` or `"segment"`.  Default: `"length"`
2495//   closed = the first and last profile are connected.  Default: false
2496function subdivide_and_slice(profiles, slices, numpoints, method="length", closed=false) =
2497  let(
2498    maxsize = max_length(profiles),
2499    numpoints = is_undef(numpoints) ? maxsize :
2500                numpoints == "lcm" ? lcmlist([for(p=profiles) len(p)]) :
2501                is_num(numpoints) ? round(numpoints) : undef
2502  )
2503  assert(is_def(numpoints), "Parameter numpoints must be \"max\", \"lcm\" or a positive number")
2504  assert(numpoints>=maxsize, "Number of points requested is smaller than largest profile")
2505  let(fixpoly = [for(poly=profiles) subdivide_path(poly, numpoints,method=method)])
2506  slice_profiles(fixpoly, slices, closed);
2507
2508
2509
2510// Function: slice_profiles()
2511// Synopsis: Linearly interpolates between path profiles.
2512// SynTags: PathList
2513// Topics: Paths, Path Subdivision
2514// See Also: subdivide_and_slice()
2515// Usage:
2516//   profs = slice_profiles(profiles, slices, [closed]);
2517// Description:
2518//   Given an input list of profiles, linearly interpolate between each pair to produce a
2519//   more finely sampled list.  The parameters `slices` specifies the number of slices to
2520//   be inserted between each pair of profiles and can be a number or a list.
2521// Arguments:
2522//   profiles = list of paths to operate on.  They must be lists of the same shape and length.
2523//   slices = number of slices to insert between each pair, or a list to vary the number inserted.
2524//   closed = set to true if last profile connects to first one.  Default: false
2525function slice_profiles(profiles,slices,closed=false) =
2526  assert(is_num(slices) || is_list(slices))
2527  let(listok = !is_list(slices) || len(slices)==len(profiles)-(closed?0:1))
2528  assert(listok, "Input slices to slice_profiles is a list with the wrong length")
2529  let(
2530    count = is_num(slices) ? repeat(slices,len(profiles)-(closed?0:1)) : slices,
2531    slicelist = [for (i=[0:len(profiles)-(closed?1:2)])
2532      each lerpn(profiles[i], select(profiles,i+1), count[i]+1, false)
2533    ]
2534  )
2535  concat(slicelist, closed?[]:[profiles[len(profiles)-1]]);
2536
2537
2538
2539function _closest_angle(alpha,beta) =
2540    is_vector(beta) ? [for(entry=beta) _closest_angle(alpha,entry)]
2541  : beta-alpha > 180 ? beta - ceil((beta-alpha-180)/360) * 360
2542  : beta-alpha < -180 ? beta + ceil((alpha-beta-180)/360) * 360
2543  : beta;
2544
2545
2546// Smooth data with N point moving average.  If angle=true handles data as angles.
2547// If closed=true assumes last point is adjacent to the first one.
2548// If closed=false pads data with left/right value (probably wrong behavior...should do linear interp)
2549function _smooth(data,len,closed=false,angle=false) =
2550  let(  halfwidth = floor(len/2),
2551        result = closed ? [for(i=idx(data))
2552                           let(
2553                             window = angle ? _closest_angle(data[i],select(data,i-halfwidth,i+halfwidth))
2554                                            : select(data,i-halfwidth,i+halfwidth)
2555                           )
2556                           mean(window)]
2557               : [for(i=idx(data))
2558                   let(
2559                       window = select(data,max(i-halfwidth,0),min(i+halfwidth,len(data)-1)),
2560                       left = i-halfwidth<0,
2561                       pad = left ? data[0] : last(data)
2562                   )
2563                   sum(window)+pad*(len-len(window))] / len
2564   )
2565   result;
2566
2567
2568// Function: rot_resample()
2569// Synopsis: Resample a list of rotation operators. 
2570// SynTags: MatList
2571// Topics: Matrices, Interpolation, Rotation
2572// See Also: subdivide_and_slice(), slice_profiles()
2573// Usage:
2574//   rlist = rot_resample(rotlist, n, [method=], [twist=], [scale=], [smoothlen=], [long=], [turns=], [closed=])
2575// Description:
2576//   Takes as input a list of rotation matrices in 3d.  Produces as output a resampled
2577//   list of rotation operators (4x4 matrixes) suitable for use with sweep().  You can optionally apply twist to
2578//   the output with the twist parameter, which is either a scalar to apply a uniform
2579//   overall twist, or a vector to apply twist non-uniformly.  Similarly you can apply
2580//   scaling either overall or with a vector.  The smoothlen parameter applies smoothing
2581//   to the twist and scaling to prevent abrupt changes.  This is done by a moving average
2582//   of the smoothing or scaling values.  The default of 1 means no smoothing.  The long parameter causes
2583//   the interpolation to be done the "long" way around the rotation instead of the short way.
2584//   Note that the rotation matrix cannot distinguish which way you rotate, only the place you
2585//   end after rotation.  Another ambiguity arises if your rotation is more than 360 degrees.
2586//   You can add turns with the turns parameter, so giving turns=1 will add 360 degrees to the
2587//   rotation so it completes one full turn plus the additional rotation given my the transform.
2588//   You can give long as a scalar or as a vector.  Finally if closed is true then the
2589//   resampling will connect back to the beginning.
2590//   .
2591//   The default is to resample based on the length of the arc defined by each rotation operator.  This produces
2592//   uniform sampling over all of the transformations.  It requires that each rotation has nonzero length.
2593//   In this case n specifies the total number of samples.  If you set method to "count" then you get
2594//   n samples for each transform.  You can set n to a vector to vary the samples at each step.
2595// Arguments:
2596//   rotlist = list of rotation operators in 3d to resample
2597//   n = Number of rotations to produce as output when method is "length" or number for each transformation if method is "count".  Can be a vector when method is "count"
2598//   ---
2599//   method = sampling method, either "length" or "count"
2600//   twist = scalar or vector giving twist to add overall or at each rotation.  Default: none
2601//   scale = scalar or vector giving scale factor to add overall or at each rotation.  Default: none
2602//   smoothlen = amount of smoothing to apply to scaling and twist.  Should be an odd integer.  Default: 1
2603//   long = resample the "long way" around the rotation, a boolean or list of booleans.  Default: false
2604//   turns = add extra turns.  If a scalar adds the turns to every rotation, or give a vector.  Default: 0
2605//   closed = if true then the rotation list is treated as closed.  Default: false
2606// Example(3D): Resampling the arc from a compound rotation with translations thrown in.
2607//   tran = rot_resample([ident(4), back(5)*up(4)*xrot(-10)*zrot(-20)*yrot(117,cp=[10,0,0])], n=25);
2608//   sweep(circle(r=1,$fn=3), tran);
2609// Example(3D): Applying a scale factor
2610//   tran = rot_resample([ident(4), back(5)*up(4)*xrot(-10)*zrot(-20)*yrot(117,cp=[10,0,0])], n=25, scale=2);
2611//   sweep(circle(r=1,$fn=3), tran);
2612// Example(3D): Applying twist
2613//   tran = rot_resample([ident(4), back(5)*up(4)*xrot(-10)*zrot(-20)*yrot(117,cp=[10,0,0])], n=25, twist=60);
2614//   sweep(circle(r=1,$fn=3), tran);
2615// Example(3D): Going the long way
2616//   tran = rot_resample([ident(4), back(5)*up(4)*xrot(-10)*zrot(-20)*yrot(117,cp=[10,0,0])], n=25, long=true);
2617//   sweep(circle(r=1,$fn=3), tran);
2618// Example(3D): Getting transformations from turtle3d
2619//   include<BOSL2/turtle3d.scad>
2620//   tran=turtle3d(["arcsteps",1,"up", 10, "arczrot", 10,170],transforms=true);
2621//   sweep(circle(r=1,$fn=3),rot_resample(tran, n=40));
2622// Example(3D): If you specify a larger angle in turtle you need to use the long argument
2623//   include<BOSL2/turtle3d.scad>
2624//   tran=turtle3d(["arcsteps",1,"up", 10, "arczrot", 10,270],transforms=true);
2625//   sweep(circle(r=1,$fn=3),rot_resample(tran, n=40,long=true));
2626// Example(3D): And if the angle is over 360 you need to add turns to get the right result.  Note long is false when the remaining angle after subtracting full turns is below 180:
2627//   include<BOSL2/turtle3d.scad>
2628//   tran=turtle3d(["arcsteps",1,"up", 10, "arczrot", 10,90+360],transforms=true);
2629//   sweep(circle(r=1,$fn=3),rot_resample(tran, n=40,long=false,turns=1));
2630// Example(3D): Here the remaining angle is 270, so long must be set to true
2631//   include<BOSL2/turtle3d.scad>
2632//   tran=turtle3d(["arcsteps",1,"up", 10, "arczrot", 10,270+360],transforms=true);
2633//   sweep(circle(r=1,$fn=3),rot_resample(tran, n=40,long=true,turns=1));
2634// Example(3D): Note the visible line at the scale transition
2635//   include<BOSL2/turtle3d.scad>
2636//   tran = turtle3d(["arcsteps",1,"arcup", 10, 90, "arcdown", 10, 90], transforms=true);
2637//   rtran = rot_resample(tran,200,scale=[1,6]);
2638//   sweep(circle(1,$fn=32),rtran);
2639// Example(3D): Observe how using a large smoothlen value eases that transition
2640//   include<BOSL2/turtle3d.scad>
2641//   tran = turtle3d(["arcsteps",1,"arcup", 10, 90, "arcdown", 10, 90], transforms=true);
2642//   rtran = rot_resample(tran,200,scale=[1,6],smoothlen=17);
2643//   sweep(circle(1,$fn=32),rtran);
2644// Example(3D): A similar issues can arise with twist, where a "line" is visible at the transition
2645//   include<BOSL2/turtle3d.scad>
2646//   tran = turtle3d(["arcsteps", 1, "arcup", 10, 90, "move", 10], transforms=true,state=[1,-.5,0]);
2647//   rtran = rot_resample(tran,100,twist=[0,60],smoothlen=1);
2648//   sweep(subdivide_path(rect([3,3]),40),rtran);
2649// Example(3D): Here's the smoothed twist transition
2650//   include<BOSL2/turtle3d.scad>
2651//   tran = turtle3d(["arcsteps", 1, "arcup", 10, 90, "move", 10], transforms=true,state=[1,-.5,0]);
2652//   rtran = rot_resample(tran,100,twist=[0,60],smoothlen=17);
2653//   sweep(subdivide_path(rect([3,3]),40),rtran);
2654// Example(3D): Toothed belt based on a list-comprehension-demos example.  This version has a smoothed twist transition.  Try changing smoothlen to 1 to see the more abrupt transition that occurs without smoothing.
2655//   include<BOSL2/turtle3d.scad>
2656//   r_small = 19;       // radius of small curve
2657//   r_large = 46;       // radius of large curve
2658//   flat_length = 100;  // length of flat belt section
2659//   teeth=42;           // number of teeth
2660//   belt_width = 12;
2661//   tooth_height = 9;
2662//   belt_thickness = 3;
2663//   angle = 180 - 2*atan((r_large-r_small)/flat_length);
2664//   beltprofile = path3d(subdivide_path(
2665//                   square([belt_width, belt_thickness],anchor=FWD),
2666//                   20));
2667//   beltrots =
2668//     turtle3d(["arcsteps",1,
2669//               "move", flat_length,
2670//               "arcleft", r_small, angle,
2671//               "move", flat_length,
2672//     // Closing path will be interpolated
2673//     //        "arcleft", r_large, 360-angle
2674//              ],transforms=true);
2675//   beltpath = rot_resample(beltrots,teeth*4,
2676//                           twist=[180,0,-180,0],
2677//                           long=[false,false,false,true],
2678//                           smoothlen=15,closed=true);
2679//   belt = [for(i=idx(beltpath))
2680//             let(tooth = floor((i+$t*4)/2)%2)
2681//             apply(beltpath[i]*
2682//                     yscale(tooth
2683//                            ? tooth_height/belt_thickness
2684//                            : 1),
2685//                   beltprofile)
2686//          ];
2687//   skin(belt,slices=0,closed=true);
2688function rot_resample(rotlist,n,twist,scale,smoothlen=1,long=false,turns=0,closed=false,method="length") =
2689    assert(is_int(smoothlen) && smoothlen>0 && smoothlen%2==1, "smoothlen must be a positive odd integer")
2690    assert(method=="length" || method=="count")
2691    let(tcount = len(rotlist) + (closed?0:-1))
2692    assert(method=="count" || is_int(n), "n must be an integer when method is \"length\"")
2693    assert(is_int(n) || is_vector(n,tcount), str("n must be scalar or vector with length ",tcount))
2694    let(
2695          count = method=="length" ? (closed ? n+1 : n)
2696                                   : (is_vector(n) ? sum(n) : tcount*n)+1  //(closed?0:1)
2697    )
2698    assert(is_bool(long) || len(long)==tcount,str("Input long must be a scalar or have length ",tcount))
2699    let(
2700        long = force_list(long,tcount),
2701        turns = force_list(turns,tcount),
2702        T = [for(i=[0:1:tcount-1]) rot_inverse(rotlist[i])*select(rotlist,i+1)],
2703        parms = [for(i=idx(T))
2704                    let(tparm = rot_decode(T[i],long[i]))
2705                    [tparm[0]+turns[i]*360,tparm[1],tparm[2],tparm[3]]
2706                ],
2707        radius = [for(i=idx(parms)) norm(parms[i][2])],
2708        length = [for(i=idx(parms)) norm([norm(parms[i][3]), parms[i][0]/360*2*PI*radius[i]])]
2709    )
2710    assert(method=="count" || all_positive(length),
2711           "Rotation list includes a repeated entry or a rotation around the origin, not allowed when method=\"length\"")
2712    let(
2713        cumlen = [0, each cumsum(length)],
2714        totlen = last(cumlen),
2715        stepsize = totlen/(count-1),
2716        samples = method=="count"
2717                  ? let( n = force_list(n,tcount))
2718                    [for(N=n) lerpn(0,1,N,endpoint=false)]
2719                  :[for(i=idx(parms))
2720                    let(
2721                        remainder = cumlen[i] % stepsize,
2722                        offset = remainder==0 ? 0
2723                                              : stepsize-remainder,
2724                        num = ceil((length[i]-offset)/stepsize)
2725                    )
2726                    count(num,offset,stepsize)/length[i]],
2727         twist = first_defined([twist,0]),
2728         scale = first_defined([scale,1]),
2729         needlast = !approx(last(last(samples)),1),
2730         sampletwist = is_num(twist) ? lerpn(0,twist,count)
2731                     : let(
2732                          cumtwist = [0,each cumsum(twist)]
2733                      )
2734                      [for(i=idx(parms)) each lerp(cumtwist[i],cumtwist[i+1],samples[i]),
2735                      if (needlast) last(cumtwist)
2736                      ],
2737         samplescale = is_num(scale) ? lerp(1,scale,lerpn(0,1,count))
2738                     : let(
2739                          cumscale = [1,each cumprod(scale)]
2740                      )
2741                      [for(i=idx(parms)) each lerp(cumscale[i],cumscale[i+1],samples[i]),
2742                       if (needlast) last(cumscale)],
2743         smoothtwist = _smooth(closed?select(sampletwist,0,-2):sampletwist,smoothlen,closed=closed,angle=true),
2744         smoothscale = _smooth(samplescale,smoothlen,closed=closed),
2745         interpolated = [
2746           for(i=idx(parms))
2747             each [for(u=samples[i]) rotlist[i] * move(u*parms[i][3]) * rot(a=u*parms[i][0],v=parms[i][1],cp=parms[i][2])],
2748           if (needlast) last(rotlist)
2749         ]
2750     )
2751     [for(i=idx(interpolated,e=closed?-2:-1)) interpolated[i]*zrot(smoothtwist[i])*scale(smoothscale[i])];
2752
2753
2754
2755
2756
2757//////////////////////////////////////////////////////////////////
2758//
2759// Minimum Distance Mapping using Dynamic Programming
2760//
2761// Given inputs of a two polygons, computes a mapping between their vertices that minimizes the sum the sum of
2762// the distances between every matched pair of vertices.  The algorithm uses dynamic programming to calculate
2763// the optimal mapping under the assumption that poly1[0] <-> poly2[0].  We then rotate through all the
2764// possible indexings of the longer polygon.  The theoretical run time is quadratic in the longer polygon and
2765// linear in the shorter one.
2766//
2767// The top level function, _skin_distance_match(), cycles through all the of the indexings of the larger
2768// polygon, computes the optimal value for each indexing, and chooses the overall best result.  It uses
2769// _dp_extract_map() to thread back through the dynamic programming array to determine the actual mapping, and
2770// then converts the result to an index repetition count list, which is passed to repeat_entries().
2771//
2772// The function _dp_distance_array builds up the rows of the dynamic programming matrix with reference
2773// to the previous rows, where `tdist` holds the total distance for a given mapping, and `map`
2774// holds the information about which path was optimal for each position.
2775//
2776// The function _dp_distance_row constructs each row of the dynamic programming matrix in the usual
2777// way where entries fill in based on the three entries above and to the left.  Note that we duplicate
2778// entry zero so account for wrap-around at the ends, and we initialize the distance to zero to avoid
2779// double counting the length of the 0-0 pair.
2780//
2781// This function builds up the dynamic programming distance array where each entry in the
2782// array gives the optimal distance for aligning the corresponding subparts of the two inputs.
2783// When the array is fully populated, the bottom right corner gives the minimum distance
2784// for matching the full input lists.  The `map` array contains a the three key values for the three
2785// directions, where _MAP_DIAG means you map the next vertex of `big` to the next vertex of `small`,
2786// _MAP_LEFT means you map the next vertex of `big` to the current vertex of `small`, and _MAP_UP
2787// means you map the next vertex of `small` to the current vertex of `big`.
2788//
2789// Return value is [min_distance, map], where map is the array that is used to extract the actual
2790// vertex map.
2791
2792_MAP_DIAG = 0;
2793_MAP_LEFT = 1;
2794_MAP_UP = 2;
2795
2796/*
2797function _dp_distance_array(small, big, abort_thresh=1/0, small_ind=0, tdist=[], map=[]) =
2798   small_ind == len(small)+1 ? [tdist[len(tdist)-1][len(big)-1], map] :
2799   let( newrow = _dp_distance_row(small, big, small_ind, tdist) )
2800   min(newrow[0]) > abort_thresh ? [tdist[len(tdist)-1][len(big)-1],map] :
2801   _dp_distance_array(small, big, abort_thresh, small_ind+1, concat(tdist, [newrow[0]]), concat(map, [newrow[1]]));
2802*/
2803
2804
2805function _dp_distance_array(small, big, abort_thresh=1/0) =
2806   [for(
2807        small_ind = 0,
2808        tdist = [],
2809        map = []
2810           ;
2811        small_ind<=len(small)+1
2812           ;
2813        newrow =small_ind==len(small)+1 ? [0,0,0] :  // dummy end case
2814                           _dp_distance_row(small,big,small_ind,tdist),
2815        tdist = concat(tdist, [newrow[0]]),
2816        map = concat(map, [newrow[1]]),
2817        small_ind = min(newrow[0])>abort_thresh ? len(small)+1 : small_ind+1
2818       )
2819     if (small_ind==len(small)+1) each [tdist[len(tdist)-1][len(big)], map]];
2820                                     //[tdist,map]];
2821
2822
2823function _dp_distance_row(small, big, small_ind, tdist) =
2824                    // Top left corner is zero because it gets counted at the end in bottom right corner
2825   small_ind == 0 ? [cumsum([0,for(i=[1:len(big)]) norm(big[i%len(big)]-small[0])]), repeat(_MAP_LEFT,len(big)+1)] :
2826   [for(big_ind=1,
2827       newrow=[ norm(big[0] - small[small_ind%len(small)]) + tdist[small_ind-1][0] ],
2828       newmap = [_MAP_UP]
2829         ;
2830       big_ind<=len(big)+1
2831         ;
2832       costs = big_ind == len(big)+1 ? [0] :    // handle extra iteration
2833                             [tdist[small_ind-1][big_ind-1],  // diag
2834                              newrow[big_ind-1],              // left
2835                              tdist[small_ind-1][big_ind]],   // up
2836       newrow = concat(newrow, [min(costs)+norm(big[big_ind%len(big)]-small[small_ind%len(small)])]),
2837       newmap = concat(newmap, [min_index(costs)]),
2838       big_ind = big_ind+1
2839   ) if (big_ind==len(big)+1) each [newrow,newmap]];
2840
2841
2842function _dp_extract_map(map) =
2843      [for(
2844           i=len(map)-1,
2845           j=len(map[0])-1,
2846           smallmap=[],
2847           bigmap = []
2848              ;
2849           j >= 0
2850              ;
2851           advance_i = map[i][j]==_MAP_UP || map[i][j]==_MAP_DIAG,
2852           advance_j = map[i][j]==_MAP_LEFT || map[i][j]==_MAP_DIAG,
2853           i = i - (advance_i ? 1 : 0),
2854           j = j - (advance_j ? 1 : 0),
2855           bigmap = concat( [j%(len(map[0])-1)] ,  bigmap),
2856           smallmap = concat( [i%(len(map)-1)]  , smallmap)
2857          )
2858        if (i==0 && j==0) each [smallmap,bigmap]];
2859
2860
2861/// Internal Function: _skin_distance_match(poly1,poly2)
2862/// Usage:
2863///   polys = _skin_distance_match(poly1,poly2);
2864/// Description:
2865///   Find a way of associating the vertices of poly1 and vertices of poly2
2866///   that minimizes the sum of the length of the edges that connect the two polygons.
2867///   Polygons can be in 2d or 3d.  The algorithm has cubic run time, so it can be
2868///   slow if you pass large polygons.  The output is a pair of polygons with vertices
2869///   duplicated as appropriate to be used as input to `skin()`.
2870/// Arguments:
2871///   poly1 = first polygon to match
2872///   poly2 = second polygon to match
2873function _skin_distance_match(poly1,poly2) =
2874   let(
2875      swap = len(poly1)>len(poly2),
2876      big = swap ? poly1 : poly2,
2877      small = swap ? poly2 : poly1,
2878      map_poly = [ for(
2879              i=0,
2880              bestcost = 1/0,
2881              bestmap = -1,
2882              bestpoly = -1
2883              ;
2884              i<=len(big)
2885              ;
2886              shifted = list_rotate(big,i),
2887              result =_dp_distance_array(small, shifted, abort_thresh = bestcost),
2888              bestmap = result[0]<bestcost ? result[1] : bestmap,
2889              bestpoly = result[0]<bestcost ? shifted : bestpoly,
2890              best_i = result[0]<bestcost ? i : best_i,
2891              bestcost = min(result[0], bestcost),
2892              i=i+1
2893              )
2894              if (i==len(big)) each [bestmap,bestpoly,best_i]],
2895      map = _dp_extract_map(map_poly[0]),
2896      smallmap = map[0],
2897      bigmap = map[1],
2898      // These shifts are needed to handle the case when points from both ends of one curve map to a single point on the other
2899      bigshift =  len(bigmap) - max(max_index(bigmap,all=true))-1,
2900      smallshift = len(smallmap) - max(max_index(smallmap,all=true))-1,
2901      newsmall = list_rotate(repeat_entries(small,unique_count(smallmap)[1]),smallshift),
2902      newbig = list_rotate(repeat_entries(map_poly[1],unique_count(bigmap)[1]),bigshift)
2903      )
2904      swap ? [newbig, newsmall] : [newsmall,newbig];
2905
2906
2907// This function associates vertices but with the assumption that index 0 is associated between the
2908// two inputs.  This gives only quadratic run time.  As above, output is pair of polygons with
2909// vertices duplicated as suited to use as input to skin().
2910
2911function _skin_aligned_distance_match(poly1, poly2) =
2912    let(
2913      result = _dp_distance_array(poly1, poly2, abort_thresh=1/0),
2914      map = _dp_extract_map(result[1]),
2915      shift0 = len(map[0]) - max(max_index(map[0],all=true))-1,
2916      shift1 = len(map[1]) - max(max_index(map[1],all=true))-1,
2917      new0 = list_rotate(repeat_entries(poly1,unique_count(map[0])[1]),shift0),
2918      new1 = list_rotate(repeat_entries(poly2,unique_count(map[1])[1]),shift1)
2919  )
2920  [new0,new1];
2921
2922
2923//////////////////////////////////////////////////////////////////////////////////////////////////////////////
2924/// Internal Function: _skin_tangent_match()
2925/// Usage:
2926///   x = _skin_tangent_match(poly1, poly2)
2927/// Description:
2928///   Finds a mapping of the vertices of the larger polygon onto the smaller one.  Whichever input is the
2929///   shorter path is the polygon, and the longer input is the curve.  For every edge of the polygon, the algorithm seeks a plane that contains that
2930///   edge and is tangent to the curve.  There will be more than one such point.  To choose one, the algorithm centers the polygon and curve on their centroids
2931///   and chooses the closer tangent point.  The algorithm works its way around the polygon, computing a series of tangent points and then maps all of the
2932///   points on the curve between two tangent points into one vertex of the polygon.  This algorithm can fail if the curve has too few points or if it is concave.
2933/// Arguments:
2934///   poly1 = input polygon
2935///   poly2 = input polygon
2936function _skin_tangent_match(poly1, poly2) =
2937    let(
2938        swap = len(poly1)>len(poly2),
2939        big = swap ? poly1 : poly2,
2940        small = swap ? poly2 : poly1,
2941        curve_offset = centroid(small)-centroid(big),
2942        cutpts = [for(i=[0:len(small)-1]) _find_one_tangent(big, select(small,i,i+1),curve_offset=curve_offset)],
2943        shift = last(cutpts)+1,
2944        newbig = list_rotate(big, shift),
2945        repeat_counts = [for(i=[0:len(small)-1]) posmod(cutpts[i]-select(cutpts,i-1),len(big))],
2946        newsmall = repeat_entries(small,repeat_counts)
2947    )
2948    assert(len(newsmall)==len(newbig), "Tangent alignment failed, probably because of insufficient points or a concave curve")
2949    swap ? [newbig, newsmall] : [newsmall, newbig];
2950
2951
2952function _find_one_tangent(curve, edge, curve_offset=[0,0,0], closed=true) =
2953    let(
2954        angles = [
2955            for (i = [0:len(curve)-(closed?1:2)])
2956            let(
2957                plane = plane3pt( edge[0], edge[1], curve[i]),
2958                tangent = [curve[i], select(curve,i+1)]
2959            ) plane_line_angle(plane,tangent)
2960        ],
2961        zero_cross = [
2962            for (i = [0:len(curve)-(closed?1:2)])
2963            if (sign(angles[i]) != sign(select(angles,i+1)))
2964            i
2965        ],
2966        d = [
2967            for (i = zero_cross)
2968            point_line_distance(curve[i]+curve_offset, edge)
2969        ]
2970    ) zero_cross[min_index(d)];
2971
2972
2973// Function: associate_vertices()
2974// Synopsis: Create vertex association to control how {{skin()}} links vertices. 
2975// SynTags: PathList
2976// Topics: Extrusion, Skinning, Paths
2977// See Also: skin()
2978// Usage:
2979//   newpoly = associate_vertices(polygons, split);
2980// Description:
2981//   Takes as input a list of polygons and duplicates specified vertices in each polygon in the list through the series so
2982//   that the input can be passed to `skin()`.  This allows you to decide how the vertices are linked up rather than accepting
2983//   the automatically computed minimal distance linkage.  However, the number of vertices in the polygons must not decrease in the list.
2984//   The output is a list of polygons that all have the same number of vertices with some duplicates.  You specify the vertex splitting
2985//   using the `split` which is a list where each entry corresponds to a polygon: split[i] is a value or list specifying which vertices in polygon i to split.
2986//   Give the empty list if you don't want a split for a particular polygon.  If you list a vertex once then it will be split and mapped to
2987//   two vertices in the next polygon.  If you list it N times then N copies will be created to map to N+1 vertices in the next polygon.
2988//   You must ensure that each mapping produces the correct number of vertices to exactly map onto every vertex of the next polygon.
2989//   Note that if you split (only) vertex i of a polygon that means it will map to vertices i and i+1 of the next polygon.  Vertex 0 will always
2990//   map to vertex 0 and the last vertices will always map to each other, so if you want something different than that you'll need to reindex
2991//   your polygons.
2992// Arguments:
2993//   polygons = list of polygons to split
2994//   split = list of lists of split vertices
2995// Example(FlatSpin,VPD=17,VPT=[0,0,2]):  If you skin together a square and hexagon using the optimal distance method you get two triangular faces on opposite sides:
2996//   sq = regular_ngon(4,side=2);
2997//   hex = apply(rot(15),hexagon(side=2));
2998//   skin([sq,hex], slices=10, refine=10, method="distance", z=[0,4]);
2999// Example(FlatSpin,VPD=17,VPT=[0,0,2]):  Using associate_vertices you can change the location of the triangular faces.  Here they are connect to two adjacent vertices of the square:
3000//   sq = regular_ngon(4,side=2);
3001//   hex = apply(rot(15),hexagon(side=2));
3002//   skin(associate_vertices([sq,hex],[[1,2]]), slices=10, refine=10, sampling="segment", z=[0,4]);
3003// Example(FlatSpin,VPD=17,VPT=[0,0,2]): Here the two triangular faces connect to a single vertex on the square.  Note that we had to rotate the hexagon to line them up because the vertices match counting forward, so in this case vertex 0 of the square matches to vertices 0, 1, and 2 of the hexagon.
3004//   sq = regular_ngon(4,side=2);
3005//   hex = apply(rot(60),hexagon(side=2));
3006//   skin(associate_vertices([sq,hex],[[0,0]]), slices=10, refine=10, sampling="segment", z=[0,4]);
3007// Example(3D): This example shows several polygons, with only a single vertex split at each step:
3008//   sq = regular_ngon(4,side=2);
3009//   pent = pentagon(side=2);
3010//   hex = hexagon(side=2);
3011//   sep = regular_ngon(7,side=2);
3012//   profiles = associate_vertices([sq,pent,hex,sep], [1,3,4]);
3013//   skin(profiles ,slices=10, refine=10, method="distance", z=[0,2,4,6]);
3014// Example(3D): The polygons cannot shrink, so if you want to have decreasing polygons you'll need to concatenate multiple results.  Note that it is perfectly ok to duplicate a profile as shown here, where the pentagon is duplicated:
3015//   sq = regular_ngon(4,side=2);
3016//   pent = pentagon(side=2);
3017//   grow = associate_vertices([sq,pent], [1]);
3018//   shrink = associate_vertices([sq,pent], [2]);
3019//   skin(concat(grow, reverse(shrink)), slices=10, refine=10, method="distance", z=[0,2,2,4]);
3020function associate_vertices(polygons, split, curpoly=0) =
3021   curpoly==len(polygons)-1 ? polygons :
3022   let(
3023      polylen = len(polygons[curpoly]),
3024      cursplit = force_list(split[curpoly])
3025   )
3026    assert(len(split)==len(polygons)-1,str(split,"Split list length mismatch: it has length ", len(split)," but must have length ",len(polygons)-1))
3027    assert(polylen<=len(polygons[curpoly+1]),str("Polygon ",curpoly," has more vertices than the next one."))
3028    assert(len(cursplit)+polylen == len(polygons[curpoly+1]),
3029           str("Polygon ", curpoly, " has ", polylen, " vertices.  Next polygon has ", len(polygons[curpoly+1]),
3030                  " vertices.  Split list has length ", len(cursplit), " but must have length ", len(polygons[curpoly+1])-polylen))
3031    assert(len(cursplit) == 0 || max(cursplit)<polylen && min(curpoly)>=0,
3032           str("Split ",cursplit," at polygon ",curpoly," has invalid vertices.  Must be in [0:",polylen-1,"]"))
3033    len(cursplit)==0 ? associate_vertices(polygons,split,curpoly+1) :
3034    let(
3035      splitindex = sort(concat(count(polylen), cursplit)),
3036      newpoly = [for(i=[0:len(polygons)-1]) i<=curpoly ? select(polygons[i],splitindex) : polygons[i]]
3037    )
3038   associate_vertices(newpoly, split, curpoly+1);
3039
3040
3041
3042// DefineHeader(Table;Headers=Texture Name|Type|Description): Texture Values
3043
3044// Section: Texturing
3045//   Some operations are able to add texture to the objects they create.  A texture can be any regularly repeated variation in the height of the surface.
3046//   To define a texture you need to specify how the height should vary over a rectangular block that will be repeated to tile the object.  Because textures
3047//   are based on rectangular tiling, this means adding textures to curved shapes may result in distortion of the basic texture unit.  For example, if you
3048//   texture a cone, the scale of the texture will be larger at the wide end of the cone and smaller at the narrower end of the cone.
3049//   .
3050//   You can specify a texture using two methods: a height field or a VNF.  For each method you also must specify the scale of the texture, which
3051//   gives the size of the rectangular unit in your object that will correspond to one texture tile.  Note that this scale does not preserve
3052//   aspect ratio: you can stretch the texture as desired.  
3053// Subsection: Height Field Texture Maps
3054//   The simplest way to specify a texture map is to give a 2d array of
3055//   height values which specify the height of the texture on a grid.
3056//   Values in the height field should range from 0 to 1.  A zero height
3057//   in the height field corresponds to the height of the surface and 1
3058//   the highest point in the texture above the surface being textured.
3059// Figure(2D,Big,NoScales,VPT=[6.21418,0.242814,0],VPD=28.8248,VPR=[0,0,0]): Here is a 2d texture described by a "grid" that just contains a single row.  Such a texture can be used to create ribbing. The texture is `[[0, 1, 1, 0]]`, and the fixture shows three repetitions of the basic texture unit.
3060//   ftex1 = [0,1,1,0,0];
3061//   stroke( transpose([count(5),ftex1]), dots=true, dots_width=3,width=.05);
3062//   right(4)stroke( transpose([count(5),ftex1]), dots=true, width=.05,dots_color="red",color="blue",dots_width=3);
3063//   right(8)stroke( transpose([count(5),ftex1]), dots=true, dots_width=3,width=.05);
3064//   stroke([[4,-.3],[8,-.3]],width=.05,endcaps="arrow2",color="black");
3065//   move([6,-.4])color("black")text("Texture Size", size=0.3,anchor=BACK);
3066// Continues:
3067//   Line segments connect the dots within the texture and also the dots between adjacent texture tiles.
3068//   The size of the texture (specified with `tex_size`) includes the segment that connects the tile to the next one.
3069//   Note that the grid is always uniformly spaced.
3070//   By default textures are created with unit depth, meaning that the top surface
3071//   of the texture is 1 unit above the surface being textured, assuming that the texture
3072//   is correctly designed to span the range from 0 to 1.  The `tex_depth` parameter can adjust
3073//   this dimension of a texture without changing anything else, and setting `tex_depth` negative
3074//   will invert a texture.
3075// Figure(2D,Big,NoScales,VPR=[0,0,0],VPT=[6.86022,-1.91238,0],VPD=28.8248):
3076//   ftex1 = [0,1,1,0,0];
3077//   left(0)color(.6*[1,1,1])rect([12,1],anchor=BACK+LEFT);
3078//   stroke( transpose([count(5),ftex1]), dots=true, dots_width=3,width=.05);
3079//   polygon( transpose([count(5),ftex1]));
3080//   right(4){stroke( transpose([count(5),ftex1]), dots=true, width=.05,dots_width=3);
3081//        polygon( transpose([count(5),ftex1]));
3082//        }
3083//   right(8){stroke( transpose([count(5),ftex1]), dots=true, dots_width=3,width=.05);
3084//             polygon( transpose([count(5),ftex1]));
3085//        }
3086//   stroke([[12.25,0],[12.25,1]],width=.05,endcaps="arrow2",color="black");
3087//   move([12.35,.5])color("black")text("Depth=1", size=0.3,anchor=LEFT);
3088//   fwd(4){
3089//   left(0)color(.6*[1,1,1])rect([12,1],anchor=BACK+LEFT);
3090//   stroke( transpose([count(5),2*ftex1]), dots=true, dots_width=3,width=.05);
3091//   polygon( transpose([count(5),2*ftex1]));
3092//   right(4){stroke( transpose([count(5),2*ftex1]), dots=true, width=.05,dots_width=3);
3093//        polygon( transpose([count(5),2*ftex1]));
3094//        }
3095//   right(8){stroke( transpose([count(5),2*ftex1]), dots=true, dots_width=3,width=.05);
3096//             polygon( transpose([count(5),2*ftex1]));
3097//        }
3098//   stroke([[12.25,0],[12.25,2]],width=.05,endcaps="arrow2",color="black");
3099//   move([12.35,1])color("black")text("Depth=2", size=0.3,anchor=LEFT);
3100//   }
3101// Continues:
3102//   If you want to keep the texture the same size but make the slope
3103//   steeper you need to add more points to make the uniform grid fine enough
3104//   to represent the slope you want.  This means that creating sharp edges
3105//   can require a large number of points, resulting in longer run times.
3106//   When using the built-in textures you can control the number of points
3107//   using the `n=` argument to {{texture()}}.  
3108// Figure(2D,Big,NoScales,VPT=[6.21418,0.242814,0],VPD=28.8248,VPR=[0,0,0]):  
3109//   ftex2 = xscale(4/11,transpose([count(12),[0,1,1,1,1,1,1,1,1,1,0,0]]));
3110//   stroke( ftex2, dots=true, dots_width=3,width=.05);
3111//   right(4)stroke( ftex2, dots=true, width=.05,dots_color="red",color="blue",dots_width=3);
3112//   right(8)stroke( ftex2, dots=true, dots_width=3,width=.05);
3113//   stroke([[4,-.3],[8,-.3]],width=.05,endcaps="arrow2",color="black");
3114//   move([6,-.4])color("black")text("Texture Size", size=0.3,anchor=BACK);
3115// Continues:
3116//   A more serious limitation of height field textures is that some shapes, such as hexagons or circles, cannot be accurately represented because
3117//   their points don't fall on any grid.  Trying to create such shapes is difficult and will require many points to approximate the
3118//   true point positions for the desired shape.  This will make the texture slow to compute.  
3119//   Another serious limitation is more subtle.  In the 2D examples above, it is obvious how to connect the
3120//   dots together.  But in 3D example we need to triangulate the points on a grid, and this triangulation is not unique.
3121//   The `style` argument lets you specify how the points are triangulated using the styles supported by {{vnf_vertex_array()}}.
3122//   In the example below we have expanded the 2D example into 3D:
3123//   ```openscad
3124//       [[0,0,0,0],
3125//        [0,1,1,0],
3126//        [0,1,1,0],
3127//        [0,0,0,0]]
3128//   ```
3129//   and we show the 3D triangulations produced by the different styles:
3130// Figure(3D,Big,NoAxes,VPR=[45.5,0,18.2],VPT=[2.3442,-6.25815,3.91529],VPD=35.5861):
3131//   tex = [
3132//          [0,0,0,0,0],
3133//          [0,1,1,0,0],
3134//          [0,1,1,0,0],
3135//          [0,0,0,0,0],
3136//          [0,0,0,0,0]       
3137//         ];
3138//   hm = [for(i=[0:4]) [for(j=[0:4]) [i,-j,tex[i][j]]]];      
3139//   types = ["quincunx", "convex", "concave","min_area", "default","alt","min_edge"]; 
3140//   grid_copies(spacing=5, n=[4,2]){
3141//     let(s = types[$row*4+$col]){
3142//       if (is_def(s)){
3143//       vnf_polyhedron(vnf_vertex_array(hm,style=s));
3144//       if ($row==1)
3145//         back(.8)right(2)rotate($vpr)color("black")text(s,size=.5,anchor=CENTER);
3146//       else
3147//         fwd(4.7)right(2)rotate($vpr)color("black")text(s,size=.5,anchor=CENTER);
3148//       }
3149//     }
3150//   }  
3151// Continues:
3152//   Note that of the seven available styles, five produce a different result.  There may exist some concave shape where none of the styles
3153//   produce the right result everywhere on the shape.  If this happens it would be another limitation of height field textures.  (If you have an
3154//   example of such a texture and shape please let us know!)
3155// Subsection: VNF Textures
3156//   VNF textures overcome all of the limitations of height field textures, but with two costs.  They can be more difficult to construct than
3157//   a simple array of height values, and they are significantly slower to compute for a tile with the same number of points.  Note, however, for
3158//   textures that don't neatly lie on a grid, a VNF tile will be more efficient than a finely sampled height field.  With VNF textures you can create
3159//   textures that have disconnected components, or concavities that cannot be expressed with a single valued height map.  However, you can also
3160//   create invalid textures that fail to close at the ends, so care is required to ensure that your resulting shape is valid.  
3161//   .
3162//   A VNF texture is defined by defining the texture tile with a VNF whose projection onto the XY plane is contained in the unit square [0,1] x [0,1] so
3163//   that the VNF can be tiled.   The VNF is tiled without a gap, matching the edges, so the vertices along corresponding edges must match to make a
3164//   consistent triangulation possible.  The VNF cannot have any X or Y values outside the interval [0,1].  If you want a valid polyhedron
3165//   that OpenSCAD will render then you need to take care with edges of the tiles that correspond to endcap faces in the textured object.
3166//   So for example, in a linear sweep, the top and bottom edges of tiles end abruptly to form the end cap of the object.  You can make a valid object
3167//   in two ways.  One way is to create a tile with a single, complete edge along Y=0, and of course a corresponding edges along Y=1.  The second way
3168//   to make a valid object is to have no points at all on the Y=0 line, and of course none on Y=1.  In this case, the resulting texture produces
3169//   a collection of disconnected objects.  Note that the Z coordinates of your tile can be anything, but for the dimensional settings on textures
3170//   to work intuitively, you should construct your tile so that Z ranges from 0 to 1.
3171// Figure(3D): This is the "hexgrid" VNF tile, which creates a hexagonal grid texture, something which doesn't work well with a height field because the edges of the hexagon don't align with the grid.  Note how the tile ranges between 0 and 1 in both X, Y and Z.  In fact, to get a proper aspect ratio in your final texture you need to use the `tex_size` parameter to introduct a sqrt(3) scale factor.  
3172//   tex = texture("hex_grid");
3173//   vnf_polyhedron(tex);
3174// Figure(3D): This is an example of a tile that has no edges at the top or bottom, so it creates disconnected rings.  See {{linear_sweep()}} for examples showing this tile in use.
3175//   shape = skin([
3176//                 rect(2/5),
3177//                 rect(2/3),
3178//                 rect(2/5)
3179//                ],
3180//                z=[0,1/2,1],
3181//                slices=0,
3182//                caps=false);
3183//   tile = move([0,1/2,2/3],yrot(90,shape));
3184//   vnf_polyhedron(tile);
3185// Continues:
3186//   A VNF texture provides a flat structure.  In order to apply this structure to a cylinder or other curved object, the VNF must be sliced
3187//   and "folded" so it can follow the curve.  This folding is controlled by the `tex_samples` parameter to {{cyl()}}, {{linear_sweep()}},
3188//   and {{rotate_sweep()}}.  Note that you specify it when you **use** the texture, not when you create it.  This differs from height
3189//   fields, where the analogous parameter is the `n=` parameter of the {{texture()}} function.  When `tex_samples` is too small, only the
3190//   points given in the VNF will follow the surface, resulting in a blocky look and geometrical artifacts.  
3191// Figure(3D,Med,NoAxes): On the left the `tex_samples` value is small and the texture is blocky.  On the right, the default value of 8 allows a reasonable fit to the cylinder. 
3192//   xdistribute(spacing=5){
3193//      cyl(d=10/PI, h=5, chamfer=0,
3194//         texture=texture("bricks_vnf"), tex_samples=1, tex_reps=[6,3], tex_depth=.2);
3195//      cyl(d=10/PI, h=5, chamfer=0,
3196//         texture=texture("bricks_vnf"), tex_samples=8, tex_reps=[6,3], tex_depth=.2);
3197//   }
3198// Continues:
3199//   Note that when the VNF is sliced,
3200//   extra points can be introduced in the interior of faces leading to unexpected irregularities in the textures, which appear
3201//   as extra triangles.  These artifacts can be minimized by making the VNF texture's faces as large as possible rather than using
3202//   a triangulated VNF, but depending on the specific VNF texture, it may be impossible to entirely eliminate them.
3203// Figure(3D,Big,NoAxes,VPR=[140.9,0,345.7],VPT=[9.48289,-0.88709,5.7837],VPD=39.5401): The left shows a normal bricks_vnf texture.  The right shows a texture that was first passed through {{vnf_triangulate()}}.  Note the extra triangle artifacts visible at the ends on the brick faces.
3204//   tex = texture("bricks_vnf");
3205//   cyl(d=10,h=15,texture=tex, tex_reps=[4,2],tex_samples=5,rounding=2);
3206//   up(7)fwd(-3)right(15)cyl(d=10,h=15,texture=vnf_triangulate(tex), tex_reps=[4,2],tex_samples=5,rounding=2);
3207
3208
3209// Function: texture()
3210// Topics: Textures, Knurling
3211// Synopsis: Produce a standard texture. 
3212// Topics: Extrusion, Textures
3213// See Also: linear_sweep(), rotate_sweep(), heightfield(), cylindrical_heightfield()
3214// Usage:
3215//   tx = texture(tex, [n=], [inset=], [gap=], [roughness=]);
3216// Description:
3217//   Given a texture name, returns a texture.  Textures can come in two varieties:
3218//   - Heightfield textures which are 2D arrays of scalars.  These are usually faster to render, but can be less precise and prone to triangulation errors.  The table below gives the recommended style for the best triangulation.  If results are still incorrect, switch to the similar VNF tile by adding the "_vnf" suffix.
3219//   - VNF Tile textures, which are VNFs that cover the unit square [0,0] x [1,1].  These tend to be slower to render, but allow greater flexibility and precision for shapes that don't align with a grid.
3220//   .
3221//   In the descriptions below, imagine the textures positioned on the XY plane, so "horizontal" refers to the "sideways" dimensions of the texture and
3222//   "up" and "down" refer to the depth dimension, perpendicular to the surface being textured.  If a texture is placed on a cylinder the "depth" will become the radial direction and the "horizontal"
3223//   direction will be the vertical and tangential directions on the cylindrical surface.  All horizontal dimensions for VNF textures are relative to the unit square
3224//   on which the textures are defined, so a value of 0.25 for a gap or border will refer to 1/4 of the texture's full length and/or width.  All supported textures appear below in the examples.  
3225// Arguments:
3226//   tex = The name of the texture to get.
3227//   ---
3228//   n = The number of samples to use for defining a heightfield texture.  Depending on the texture, result will be either n x n or 1 x n.  Not allowed for VNF textures.  See the `tex_samples` argument to {{cyl()}}, {{linear_sweep()}} and {{rotate_sweep()}} for controlling the sampling of VNF textures.
3229//   border = The size of a border region on some VNF tile textures.  Generally between 0 and 0.5.
3230//   gap = The gap between logically distinct parts of some VNF tiles.  (ie: gap between bricks, gap between truncated ribs, etc.)
3231//   roughness = The amount of roughness used on the surface of some heightfield textures.  Generally between 0 and 0.5.
3232// Example(3D): **"bricks"** (Heightfield) = A brick-wall pattern.  Giving `n=` sets the number of heightfield samples to `n x n`.  Default: 24.  Giving `roughness=` adds a level of height randomization to add roughness to the texture.  Default: 0.05.  Use `style="convex"`.
3233//   tex = texture("bricks");
3234//   linear_sweep(
3235//       rect(30), texture=tex, h=30,
3236//       tex_size=[10,10]
3237//   );
3238// Example(3D): **"bricks_vnf"** (VNF) = VNF version of "bricks".  Giving `gap=` sets the "mortar" gap between adjacent bricks, default 0.05.  Giving `border=` specifies that the top face of the brick is smaller than the bottom of the brick by `border` on each of the four sides.  If `gap` is zero then a `border` value close to 0.5 will cause bricks to come to a sharp pointed edge, with just a tiny flat top surface.  Note that `gap+border` must be strictly smaller than 0.5.   Default is `border=0.05`.  
3239//   tex = texture("bricks_vnf");
3240//   linear_sweep(
3241//       rect(30), texture=tex, h=30,
3242//       tex_size=[10,10]
3243//   );
3244// Example(3D): "bricks_vnf" texture with large border. 
3245//   tex = texture("bricks_vnf",border=0.25);
3246//   linear_sweep(
3247//       rect(30), texture=tex, h=30,
3248//       tex_size=[10,10]
3249//   );
3250// Example(3D,VPR=[84.4,0,4.7],VPT=[2.44496,6.53317,14.6135],VPD = 126): **"checkers"** (VNF) = A pattern of alternating checkerboard squares.  Giving `border=` specifies that the top face of the checker surface is smaller than the bottom by `border` on each of the four sides.  As `border` approaches 0.5 the tops come to sharp corners.  You must set `border` strictly between 0 and 0.5.  Default: 0.05.
3251//   tex = texture("checkers");
3252//   linear_sweep(
3253//       rect(30), texture=tex, h=30,
3254//       tex_size=[10,10]
3255//   );
3256// Example(3D,VPR=[84.4,0,4.7],VPT=[2.44496,6.53317,14.6135],VPD = 126): "checkers" texture with large border.  
3257//   tex = texture("checkers",border=0.25);
3258//   linear_sweep(
3259//       rect(30), texture=tex, h=30,
3260//       tex_size=[10,10]
3261//   );
3262// Example(3D): **"cones"** (VNF) = Raised conical spikes.  Specify `$fn` to set the number of segments on the cone (will be rounded to a multiple of 4).  The default is `$fn=16`.  Note that `$fa` and `$fs` are ignored, since the scale of the texture is unknown at the time of definition.  Giving `border=` specifies the horizontal border width between the edge of the tile and the base of the cone.  The `border` value must be nonnegative and smaller than 0.5.  Default: 0.
3263//   tex = texture("cones", $fn=16);
3264//   linear_sweep(
3265//       rect(30), texture=tex, h=30, tex_depth=3,
3266//       tex_size=[10,10]
3267//   );
3268// Example(3D): **"cubes"** (VNF) = Corner-cubes texture.  This texture needs to be scaled in vertically by sqrt(3) to have its correct aspect
3269//   tex = texture("cubes");
3270//   linear_sweep(
3271//       rect(30), texture=tex, h=30,
3272//       tex_size=[10,10]
3273//   );
3274// Example(3D): "cubes" texture at the correct scale.  
3275//   tex = texture("cubes");
3276//   linear_sweep(
3277//       rect(30), texture=tex, h=20*sqrt(3), tex_depth=3,
3278//       tex_size=[10,10*sqrt(3)]
3279//   );
3280// Example(3D): **"diamonds"** (Heightfield) = Four-sided pyramid with the corners of the base aligned with the axes.  Compare to "pyramids".  Useful for knurling.  Giving `n=` sets the number of heightfield samples to `n x n`. Default: 2.  Use `style="concave"` for pointed bumps, or `style="default"` or `style="alt"` for a diagonal ribs.  
3281//   tex = texture("diamonds");
3282//   linear_sweep(
3283//       rect(30), texture=tex, h=30,
3284//       tex_size=[10,10], style="concave"
3285//   );
3286// Example(3D): "diamonds" texture can give diagonal ribbing with "default" style. 
3287//   tex = texture("diamonds");
3288//   linear_sweep(
3289//       rect(30), texture=tex, h=30,
3290//       tex_size=[10,10], style="default"
3291//   );
3292// Example(3D): "diamonds" texture gives diagonal ribbing the other direction with "alt" style.  
3293//   tex = texture("diamonds");
3294//   linear_sweep(
3295//       rect(30), texture=tex, h=30,
3296//       tex_size=[10,10], style="alt"
3297//   );
3298// Example(3D): **"diamonds_vnf"** (VNF) = VNF version of "diamonds".
3299//   tex = texture("diamonds_vnf");
3300//   linear_sweep(
3301//       rect(30), texture=tex, h=30,
3302//       tex_size=[10,10]
3303//   );
3304// Example(3D): **"dimples"** (VNF) = Round divots.  Specify `$fn` to set the number of segments on the dimples (will be rounded to a multiple of 4).  The default is `$fn=16`.  Note that `$fa` and `$fs` are ignored, since the scale of the texture is unknown at the time of definition.  Giving `border=` specifies the horizontal width of the flat border region between the tile edges and the edge of the dimple.  Must be nonnegative and strictly less than 0.5.  Default: 0.05.  
3305//   tex = texture("dimples", $fn=16);
3306//   linear_sweep(
3307//       rect(30), texture=tex, h=30, 
3308//       tex_size=[10,10]
3309//   );
3310// Example(3D): **"dots"** (VNF) = Raised round bumps.  Specify `$fn` to set the number of segments on the dots (will be rounded to a multiple of 4).  The default is `$fn=16`.  Note that `$fa` and `$fs` are ignored, since the scale of the texture is unknown at the time of definition.  Giving `border=` specifies the horizontal width of the flat border region between the tile edge and the edge of the dots.  Must be nonnegative and strictly less than 0.5.  Default: 0.05.
3311//   tex = texture("dots", $fn=16);
3312//   linear_sweep(
3313//       rect(30), texture=tex, h=30,
3314//       tex_size=[10,10]
3315//   );
3316// Example(3D): **"hex_grid"** (VNF) = A hexagonal grid defined by V-grove borders.  Giving `border=` specifies that the top face of the hexagon is smaller than the bottom by `border` on the left and right sides.  This means the V-groove top width for grooves running parallel to the Y axis will be double the border value.  If the texture is scaled in the Y direction by sqrt(3) then the groove will be uniform on all six sides of the hexagon.  Border must be strictly between 0 and 0.5, default: 0.1.
3317//   tex = texture("hex_grid");
3318//   linear_sweep(
3319//       rect(30), texture=tex, h=30,
3320//       tex_size=[10,10]
3321//   );
3322// Example(3D): "hex_grid" texture with large border
3323//   tex = texture("hex_grid", border=0.4);
3324//   linear_sweep(
3325//       rect(30), texture=tex, h=30,
3326//       tex_size=[10,10]
3327//   );
3328// Example(3D): "hex_grid" scaled in Y by sqrt(3) so hexagons are regular and grooves are all the same width.  Note height of cube is also scaled so tile fits without being automatically adjusted to fit, ruining our choice of scale.
3329//   tex = texture("hex_grid",border=.07);
3330//   linear_sweep(
3331//       rect(30), texture=tex, h=quantup(30,10*sqrt(3)),
3332//       tex_size=[10,10*sqrt(3)], tex_depth=3
3333//   );
3334// Example(3D): "hex_grid" texture, with approximate scaling because 17 is close to sqrt(3) times 10.
3335//   tex = texture("hex_grid");
3336//   linear_sweep(
3337//       rect(30), texture=tex, h=34,
3338//       tex_size=[10,17]
3339//   );
3340// Example(3D): **"hills"** (Heightfield) = Wavy sine-wave hills and valleys,  Giving `n=` sets the number of heightfield samples to `n` x `n`.  Default: 12.  Set `style="quincunx"`.
3341//   tex = texture("hills");
3342//   linear_sweep(
3343//       rect(30), texture=tex, h=30,
3344//       tex_size=[10,10], style="quincunx"
3345//   );
3346// Example(3D): **"pyramids"** (Heightfield) = Four-sided pyramid with the edges of the base aligned with the axess.  Compare to "diamonds". Useful for knurling.  Giving `n=` sets the number of heightfield samples to `n` by `n`. Default: 2. Set style to "convex".  Note that style="concave" or style="min_edge" produce mini-diamonds with flat squares in between.
3347//   tex = texture("pyramids");
3348//   linear_sweep(
3349//       rect(30), texture=tex, h=30,
3350//       tex_size=[10,10], style="convex"
3351//   );
3352// Example(3D): "pyramids" texture, with "concave" produces a mini-diamond texture.  Note that "min_edge" also gives this result.
3353//   tex = texture("pyramids");
3354//   linear_sweep(
3355//       rect(30), texture=tex, h=30,
3356//       tex_size=[10,10], style="concave"
3357//   );
3358// Example(3D): **"pyramids_vnf"** (VNF) = VNF version of "pyramids".
3359//   tex = texture("pyramids_vnf");
3360//   linear_sweep(
3361//       rect(30), texture=tex, h=30,
3362//       tex_size=[10,10]
3363//   );
3364// Example(3D): **"ribs"** (Heightfield) = Vertically aligned triangular ribs.  Giving `n=` sets the number of heightfield samples to `n` by 1.  Default: 2.  The choice of style does not matter.
3365//   tex = texture("ribs");
3366//   linear_sweep(
3367//       rect(30), texture=tex, h=30, tex_depth=3,
3368//       tex_size=[10,10], style="concave"
3369//   );
3370// Example(3D): **"rough"** (Heightfield) = A pseudo-randomized rough texture.  Giving `n=` sets the number of heightfield samples to `n` by `n`.  Default: 32.  The `roughness=` parameter specifies the height of the random texture.  Default: 0.2.
3371//   tex = texture("rough");
3372//   linear_sweep(
3373//       rect(30), texture=tex, h=30,
3374//       tex_size=[10,10], style="min_edge"
3375//   );
3376// Example(3D): **"tri_grid"** (VNF) = A triangular grid defined by V-groove borders  Giving `border=` specifies that the top face of the triangular surface is smaller than the bottom by `border` along the horizontal edges (parallel to the X axis).  This means the V-groove top width of the grooves parallel to the X axis will be double the border value.  (The other grooves are wider.) If the tile is scaled in the Y direction by sqrt(3) then the groove will be uniform on the three sides of the triangle.  The border must be strictly between 0 and 1/6, default: 0.05.
3377//   tex = texture("tri_grid");
3378//   linear_sweep(
3379//       rect(30), texture=tex, h=30,
3380//       tex_size=[10,10]
3381//   );
3382// Example(3D): "tri_grid" texture with large border.  (Max border for tri_grid is 1/6.)  
3383//   tex = texture("tri_grid",border=.12);
3384//   linear_sweep(
3385//       rect(30), texture=tex, h=30,
3386//       tex_size=[10,10]
3387//   );
3388// Example(3D): "tri_grid" texture scaled in Y by sqrt(3) so triangles are equilateral and grooves are all the same width.  Note we have to ensure the height evenly fits the scaled texture tiles.
3389//   tex = texture("tri_grid",border=.04);
3390//   linear_sweep(
3391//       rect(30), texture=tex, h=quantup(30,10*sqrt(3)),
3392//       tex_size=[10,10*sqrt(3)], tex_depth=3
3393//   );
3394// Example(3D): "tri_grid" texture.  Here scale makes Y approximately sqrt(3) larger than X so triangles are close to equilateral.
3395//   tex = texture("tri_grid");
3396//   linear_sweep(
3397//       rect(30), texture=tex, h=34,
3398//       tex_size=[10,17]
3399//   );
3400// Example(3D): **"trunc_diamonds"** (VNF) = Truncated diamonds, four-sided pyramids with the base corners aligned with the axes and the top cut off.  Or you can interpret it as V-groove lines at 45º angles.  Giving `border=` specifies that the width and height of the top surface of the diamond are smaller by `border` at the left, right, top and bottom.  The border is measured in the **horizontal** direction.  This means the V-groove width will be sqrt(2) times the border value.  The border must be strictly between 0 and sqrt(2)/4, which is about 0.35.  Default: 0.1.
3401//   tex = texture("trunc_diamonds");
3402//   linear_sweep(
3403//       rect(30), texture=tex, h=30,
3404//       tex_size=[10,10]
3405//   );
3406// Example(3D): "trunc_diamonds" texture with large border. 
3407//   tex = texture("trunc_diamonds",border=.25);
3408//   linear_sweep(
3409//       rect(30), texture=tex, h=30,
3410//       tex_size=[10,10]
3411//   );
3412// Example(3D): **"trunc_pyramids"** (Heightfield) = Truncated pyramids, four sided pyramids with the base edges aligned to the axes and the top cut off.  Giving `n=` sets the number of heightfield samples to `n` by `n`.  Default: 6.  Set `style="convex"`.
3413//   tex = texture("trunc_pyramids");
3414//   linear_sweep(
3415//       rect(30), texture=tex, h=30,
3416//       tex_size=[10,10], style="convex"
3417//   );
3418// Example(3D): **"trunc_pyramids_vnf"** (VNF) = Truncated pyramids, four sided pyramids with the base edges aligned to the axes and the top cut off.  You can also regard this as a grid of V-grooves.  Giving `border=` specifies that the top face is smaller than the top by `border` on all four sides.  This means the V-groove top width will be double the border value.  The border must be strictly between 0 and 0.5.  Default: 0.1.
3419//   tex = texture("trunc_pyramids_vnf");
3420//   linear_sweep(
3421//       rect(30), texture=tex, h=30,
3422//       tex_size=[10,10]
3423//   );
3424// Example(3D): "trunc_pyramids_vnf" texture with large border
3425//   tex = texture("trunc_pyramids_vnf", border=.4);
3426//   linear_sweep(
3427//       rect(30), texture=tex, h=30,
3428//       tex_size=[10,10]
3429//   );
3430// Example(3D): **"trunc_ribs"** (Heightfield) = Truncated ribs.  Vertically aligned triangular ribs with the tops cut off, and with rib separation equal to the width of the flat tops.  Giving `n=` sets the number of heightfield samples to `n` by `1`.  Default: 4.  The style does not matter.
3431//   tex = texture("trunc_ribs");
3432//   linear_sweep(
3433//       rect(30), h=30, texture=tex,
3434//       tex_depth=3, tex_size=[10,10],
3435//       style="concave"
3436//   );
3437// Example(3D): **"trunc_ribs_vnf"** (VNF) = Vertically aligned triangular ribs with the tops cut off.  Giving `gap=` sets the bottom gap between ribs.  Giving `border=` specifies that the top rib face is smaller than its base by `border` on both the left and right sides.  The gap measures the flat part between ribs and the border the width of the sloping portion. In order to fit, gap+2*border must be less than 1.  (This is because the gap is counted once but the border counts on both sides.)  Defaults: gap=1/4, border=1/4.
3438//   tex = texture("trunc_ribs_vnf", gap=0.25, border=1/6);
3439//   linear_sweep(
3440//       rect(30), h=30, texture=tex,
3441//       tex_depth=3, tex_size=[10,10]
3442//   );
3443// Example(3D): **"wave_ribs"** (Heightfield) = Vertically aligned wavy ribs.  Giving `n=` sets the number of heightfield samples to `n` by `1`.  Default: 8.  The style does not matter.  
3444//   tex = texture("wave_ribs");
3445//   linear_sweep(
3446//       rect(30), h=30, texture=tex, 
3447//       tex_size=[10,10], tex_depth=3, style="concave"
3448//   );
3449
3450
3451function _tex_fn_default() = 16;
3452
3453__vnf_no_n_mesg=" texture is a VNF so it does not accept n.  Set sample rate for VNF textures using the tex_samples parameter to cyl(), linear_sweep() or rotate_sweep().";
3454
3455function texture(tex, n, border, gap, roughness, inset) =
3456    assert(num_defined([border,inset])<2, "In texture() the 'inset' parameter has been replaced by 'border'.  You cannot give both parameters.")
3457    let(
3458        border = is_def(inset)?echo("In texture() the argument 'inset' has been deprecated and will be removed.  Use 'border' instead")
3459                               inset
3460                              :border
3461    )
3462    assert(is_undef(n) || all_positive([n]), "n must be a positive value if given")
3463    assert(is_undef(border) || is_finite(border), "border must be a number if given")
3464    assert(is_undef(gap) || is_finite(gap), "gap must be a number if given")
3465    assert(is_undef(roughness) || all_nonnegative([roughness]), "roughness must be a nonnegative value if given")  
3466    tex=="ribs"?
3467        assert(num_defined([gap, border, roughness])==0, "ribs texture does not accept gap, border or roughness")
3468
3469        let(
3470            n = quantup(default(n,2),2)
3471        ) [[
3472            each lerpn(1,0,n/2,endpoint=false),
3473            each lerpn(0,1,n/2,endpoint=false),
3474        ]] :
3475    tex=="trunc_ribs"?
3476        assert(num_defined([gap, border, roughness])==0, "trunc_ribs texture does not accept gap, border or roughness")
3477        let(
3478            n = quantup(default(n,4),4)
3479        ) [[
3480            each repeat(0,n/4),
3481            each lerpn(0,1,n/4,endpoint=false),
3482            each repeat(1,n/4),
3483            each lerpn(1,0,n/4,endpoint=false),
3484        ]] :
3485    tex=="trunc_ribs_vnf"?
3486        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3487        let(
3488            border = default(border,1/4)*2,
3489            gap = default(gap,1/4),
3490            f=echo(gap, border, gap+border, gap+2*border)
3491        )
3492        assert(all_nonnegative([border,gap]), "trunc_ribs_vnf texture requires gap>=0 and border>=0")
3493        assert(gap+border <= 1, "trunc_ribs_vnf texture requires that gap+2*border<=1")
3494        [
3495            [
3496               each move([0.5,0.5], p=path3d(rect([1-gap,1]),0)),
3497               each move([0.5,0.5], p=path3d(rect([1-gap-border,1]),1)),
3498               each path3d(square(1)),
3499            ], [
3500                [4,7,3,0], [1,2,6,5],
3501                if (gap+border < 1-EPSILON) [4,5,6,7],
3502                if (gap > EPSILON) each [[1,9,10,2], [0,3,11,8]],
3503            ]
3504        ] :
3505    tex=="wave_ribs"?
3506        assert(num_defined([gap, border, roughness])==0, "wave_ribs texture does not accept gap, border or roughness")  
3507        let(
3508            n = max(6,default(n,8))
3509        ) [[
3510            for(a=[0:360/n:360-EPSILON])
3511            (cos(a)+1)/2
3512        ]] :
3513    tex=="diamonds"?
3514        assert(num_defined([gap, border, roughness])==0, "diamonds texture does not accept gap, border or roughness")  
3515        let(
3516            n = quantup(default(n,2),2)
3517        ) [
3518            let(
3519                path = [
3520                    each lerpn(0,1,n/2,endpoint=false),
3521                    each lerpn(1,0,n/2,endpoint=false),
3522                ]
3523            )
3524            for (i=[0:1:n-1]) [
3525                for (j=[0:1:n-1]) min(
3526                    select(path,i+j),
3527                    select(path,i-j)
3528                )
3529            ],
3530        ] :
3531    tex=="diamonds_vnf"?
3532        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3533        assert(num_defined([gap, border, roughness])==0, "diamonds_vnf texture does not accept gap, border or roughness")
3534        [
3535            [
3536                [0,   1, 1], [1/2,   1, 0], [1,   1, 1],
3537                [0, 1/2, 0], [1/2, 1/2, 1], [1, 1/2, 0],
3538                [0,   0, 1], [1/2,   0, 0], [1,   0, 1],
3539            ], [
3540                [0,1,3], [2,5,1], [8,7,5], [6,3,7],
3541                [1,5,4], [5,7,4], [7,3,4], [4,3,1],
3542            ]
3543        ] :
3544    tex=="pyramids"?
3545        assert(num_defined([gap, border, roughness])==0, "pyramids texture does not accept gap, border or roughness")
3546        let(
3547            n = quantup(default(n,2),2)
3548        ) [
3549            for (i = [0:1:n-1]) [
3550                for (j = [0:1:n-1])
3551                1 - (max(abs(i-n/2), abs(j-n/2)) / (n/2))
3552            ]
3553        ] :
3554    tex=="pyramids_vnf"?
3555        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3556        assert(num_defined([gap, border, roughness])==0, "pyramids_vnf texture does not accept gap, border or roughness")  
3557        [
3558            [ [0,1,0], [1,1,0], [1/2,1/2,1], [0,0,0], [1,0,0] ],
3559            [ [2,0,1], [2,1,4], [2,4,3], [2,3,0] ]
3560        ] :
3561    tex=="trunc_pyramids"?
3562        assert(num_defined([gap, border, roughness])==0, "trunc_pyramids texture does not accept gap, border or roughness")  
3563        let(
3564            n = quantup(default(n,6),3)
3565        ) [
3566            for (i = [0:1:n-1]) [
3567                for (j = [0:1:n-1])
3568                (1 - (max(n/6, abs(i-n/2), abs(j-n/2)) / (n/2))) * 1.5
3569            ]
3570        ] :
3571    tex=="trunc_pyramids_vnf"?
3572        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3573        assert(num_defined([gap, roughness])==0, "trunc_pyramids_vnf texture does not accept gap, or roughness")
3574        let(
3575            border = default(border,0.1)
3576        )
3577        assert(border>0 && border<.5, "trunc_pyramids_vnf texture requires border in (0,0.5)")
3578        [
3579            [
3580                each path3d(square(1)),
3581                each move([1/2,1/2,1], p=path3d(rect(1-2*border))),
3582            ], [
3583                for (i=[0:3])
3584                    [i, (i+1)%4, (i+1)%4+4,i+4],
3585                [4,5,6,7]
3586            ]
3587        ] :
3588    tex=="hills"?
3589        assert(num_defined([gap, border, roughness])==0, "hills texture does not accept gap, border or roughness")  
3590        let(
3591            n = default(n,12)
3592        ) [
3593            for (a=[0:360/n:359.999]) [
3594                for (b=[0:360/n:359.999])
3595                (cos(a)*cos(b)+1)/2
3596            ]
3597        ] :
3598    tex=="bricks"?
3599        assert(num_defined([gap,border])==0, "bricks texture does not accept gap or border")  
3600        let(
3601            n = quantup(default(n,24),2),
3602            rough = default(roughness,0.05)
3603        ) [
3604            for (y = [0:1:n-1])
3605            rands(-rough/2, rough/2, n, seed=12345+y*678) + [
3606                for (x = [0:1:n-1])
3607                (y%(n/2) <= max(1,n/16))? 0 :
3608                let( even = floor(y/(n/2))%2? n/2 : 0 )
3609                (x+even) % n <= max(1,n/16)? 0 : 0.5
3610            ]
3611        ] :
3612    tex=="bricks_vnf"?
3613        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3614        assert(num_defined([roughness])==0, "bricks_vnf texture does not accept roughness")
3615        let(
3616            border = default(border,0.05),
3617            gap = default(gap,0.05)
3618        )
3619        assert(border>=0,"bricks_vnf texture requires nonnegative border")
3620        assert(gap>0, "bricks_vnf requires gap greater than 0")
3621        assert(gap+border<0.5, "bricks_vnf requires gap+border<0.5")
3622          [
3623            [
3624                each path3d(square(1)),
3625                each move([gap/2, gap/2, 0], p=path3d(square([1-gap, 0.5-gap]))),
3626                each move([gap/2+border/2, gap/2+border/2, 1], p=path3d(square([1-gap-border, 0.5-gap-border]))),
3627                each move([0, 0.5+gap/2, 0], p=path3d(square([0.5-gap/2, 0.5-gap]))),
3628                each move([0, 0.5+gap/2+border/2, 1], p=path3d(square([0.5-gap/2-border/2, 0.5-gap-border]))),
3629                each move([0.5+gap/2, 0.5+gap/2, 0], p=path3d(square([0.5-gap/2, 0.5-gap]))),
3630                each move([0.5+gap/2+border/2, 0.5+gap/2+border/2, 1], p=path3d(square([0.5-gap/2-border/2, 0.5-gap-border]))),
3631            ], [
3632                [0,4,7,20], [4,8,11,7], [9,8,4,5], [4,0,1,5], [10,9,5,6],
3633                [20,7,6,13,12,21] ,[2,3,23,22,15,14], [15,19,18,14], [22,23,27,26], [16,19,15,12],[13,6,5,1],
3634                [26,25,21,22], [8,9,10,11],[7,11,10,6],[17,16,12,13],[22,21,12,15],[16,17,18,19],[24,25,26,27],[25,24,20,21]
3635            ]
3636        ] :
3637    tex=="checkers"?
3638        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3639        assert(num_defined([gap, roughness])==0, "checkers texture does not accept gap, or roughness")
3640        let(
3641            border = default(border,0.05)
3642        )
3643        assert(border>0 && border<.5, "checkers texture requires border in (0,0.5)")
3644          [
3645            [
3646                each move([0,0], p=path3d(square(0.5-border),1)),
3647                each move([0,0.5], p=path3d(square(0.5-border))),
3648                each move([0.5,0], p=path3d(square(0.5-border))),
3649                each move([0.5,0.5], p=path3d(square(0.5-border),1)),
3650                [1/2-border/2,1/2-border/2,1/2], [0,1,1], [1/2-border,1,1],
3651                [1/2,1,0], [1-border,1,0], [1,0,1], [1,1/2-border,1],
3652                [1,1/2,0], [1,1-border,0], [1,1,1], [1/2-border/2,1-border/2,1/2],
3653                [1-border/2,1-border/2,1/2], [1-border/2,1/2-border/2,1/2],
3654            ], [
3655                for (i=[0:4:12]) each [[i,i+1,i+2,i+3]],
3656                [10,16,13,12,28,11],[9,0,3,16,10], [11,28,22,21,8],
3657                [4,7,26,14,13,16], [7,6,17,18,26], [5,4,16,3,2],
3658                [19,20,27,15,14,26], [20,25,27], [19,26,18],
3659                [23,28,12,15,27,24], [23,22,28], [24,27,25]
3660            ]
3661        ] :
3662    tex=="cones"?
3663        assert(is_undef(n),str("To set number of segments on cones use $fn. ", tex,__vnf_no_n_mesg))
3664        assert(num_defined([gap,roughness])==0, "cones texture does not accept gap or roughness")  
3665        let(
3666            border = default(border,0),
3667            n = $fn > 0 ? quantup($fn,4) : _tex_fn_default()
3668        )
3669        assert(border>=0 && border<0.5)
3670        [
3671            [
3672                each move([1/2,1/2], p=path3d(circle(d=1-2*border,$fn=n))),
3673                [1/2,1/2,1],
3674                each border>0 ? path3d(subdivide_path(square(1),refine=2,closed=true))
3675                              : path3d(square(1))
3676            ], [
3677                for (i=[0:1:n-1]) [i, (i+1)%n, n],
3678                if (border>0) for (i=[0:3]) [for(j=[(i+1)*n/4:-1:i*n/4]) j%n,
3679                                            (2*i+7)%8+n+1,(2*i)%8+n+1, (2*i+1)%8+n+1],
3680                if (border==0) for (i=[0:3]) [for(j=[(i+1)*n/4:-1:i*n/4]) j%n, i+n+1]
3681            ]
3682        ] :
3683    tex=="cubes"?
3684        assert(is_undef(n), str(tex,__vnf_no_n_mesg))
3685        assert(num_defined([gap, border, roughness])==0, "cubes texture does not accept gap, border or roughness")  
3686        [
3687            [
3688                [0,1,1/2], [1,1,1/2], [1/2,5/6,1], [0,4/6,0], [1,4/6,0],
3689                [1/2,3/6,1/2], [0,2/6,1], [1,2/6,1], [1/2,1/6,0], [0,0,1/2],
3690                [1,0,1/2],
3691            ], [
3692                [0,1,2], [0,2,3], [1,4,2], [2,5,3], [2,4,5],
3693                [6,3,5], [4,7,5], [7,8,5], [6,5,8], [10,8,7],
3694                [9,6,8], [10,9,8],
3695            ]
3696        ] :
3697    tex=="trunc_diamonds"?
3698        assert(is_undef(n), str(tex,__vnf_no_n_mesg))  
3699        assert(num_defined([gap, roughness])==0, "trunc_diamonds texture does not accept gap or roughness")
3700        let(
3701            border = default(border,0.1)/sqrt(2)*2
3702        )
3703        assert(border>0 && border<0.5)
3704        [
3705            [
3706                each move([1/2,1/2,0], p=path3d(circle(d=1,$fn=4))),
3707                each move([1/2,1/2,1], p=path3d(circle(d=1-border*2,$fn=4))),
3708                for (a=[0:90:359]) each move([1/2,1/2], p=zrot(-a, p=[[1/2,border,1], [border,1/2,1], [1/2,1/2,1]]))
3709            ], [
3710                for (i=[0:3]) each let(j=i*3+8) [
3711                    [i,(i+1)%4,(i+1)%4+4,i+4],
3712                    [j,j+1,j+2], [i, (i+3)%4,j+1, j], 
3713                ],
3714                [4,5,6,7],
3715            ]
3716        ] :
3717    tex=="dimples" || tex=="dots" ?
3718        assert(is_undef(n),str("To set number of segments on ",tex," use $fn. ", tex,__vnf_no_n_mesg))
3719        assert(num_defined([gap,roughness])==0, str(tex," texture does not accept gap or roughness"))
3720        let(
3721            border = default(border,0.05),
3722            n = $fn > 0 ? quantup($fn,4) : _tex_fn_default()
3723        )
3724        assert(border>=0 && border < 0.5)
3725        let(
3726            rows=ceil(n/4),
3727            r=adj_ang_to_hyp(1/2-border,45),
3728            dots = tex=="dots",
3729            cp = [1/2, 1/2, r*sin(45)*(dots?-1:1)],
3730            sc = 1 / (r - abs(cp.z)),
3731            uverts = [
3732                for (p=[0:1:rows-1], t=[0:360/n:359.999])
3733                    cp + (
3734                        dots? spherical_to_xyz(r, -t, 45-45*p/rows) :
3735                        spherical_to_xyz(r, -t, 135+45*p/rows)
3736                    ),
3737                cp + r * (dots?UP:DOWN),
3738                each border>0 ? path3d(subdivide_path(square(1),refine=2,closed=true))
3739                              : path3d(square(1)),
3740                      
3741            ],
3742            verts = zscale(sc, p=uverts),
3743            faces = [
3744                for (i=[0:1:rows-2], j=[0:1:n-1]) each [
3745                    [i*n+j, i*n+(j+1)%n, (i+1)*n+(j+1)%n,(i+1)*n+j],
3746                ],
3747                for (i=[0:1:n-1]) [(rows-1)*n+i, (rows-1)*n+(i+1)%n, rows*n],
3748                if (border>0) for (i=[0:3]) [for(j=[(i+1)*n/4:-1:i*n/4]) j%n,
3749                                            (2*i+7)%8+rows*n+1,(2*i)%8+rows*n+1, (2*i+1)%8+rows*n+1],
3750                if (border==0) for (i=[0:3]) [for(j=[(i+1)*n/4:-1:i*n/4]) j%n, i+rows*n+1]
3751            ]
3752        ) [verts, faces] :
3753    tex=="tri_grid"?
3754        assert(is_undef(n), str(tex,__vnf_no_n_mesg))  
3755        assert(num_defined([gap, roughness])==0, str(tex," texture does not accept gap or roughness"))  
3756        let(
3757            border = default(border,0.05)*sqrt(3)
3758        )
3759        assert(border>0 && border<sqrt(3)/6, "tri_grid texture requires border in (0,1/6)")
3760        let(
3761            adj = opp_ang_to_adj(border, 30),
3762            y1 = border / adj_ang_to_opp(1,60),     // i/sqrt(3)
3763            y2 = 2*y1,            // 2*i/sqrt(3)
3764            y3 = 0.5 - y1,
3765            y4 = 0.5 + y1,
3766            y5 = 1 - y2,
3767            y6 = 1 - y1
3768        )
3769        [
3770            [
3771                [0,0,0], [1,0,0],
3772                [adj,y1,1], [1-adj,y1,1],
3773                [0,y2,1], [1,y2,1],
3774                [0.5,0.5-y2,1],
3775                [0,y3,1], [0.5-adj,y3,1], [0.5+adj,y3,1], [1,y3,1],
3776                [0,0.5,0], [0.5,0.5,0], [1,0.5,0],
3777                [0,y4,1], [0.5-adj,y4,1], [0.5+adj,y4,1], [1,y4,1],
3778                [0.5,0.5+y2,1],
3779                [0,y5,1], [1,y5,1],
3780                [adj,y6,1], [1-adj,y6,1],
3781                [0,1,0], [1,1,0],
3782            ], [
3783               [0,2,3,1],
3784               [21,23,24,22],
3785               [2,6,3], [0,12,6,2], [1,3,6,12],
3786               [0,4,8,12], [4,7,8], [8,7,11,12],
3787               [1,12,9,5], [5,9,10], [10,9,12,13], 
3788               [11,14,15,12], [19,15,14], [19,23,12,15],
3789               [16,17,13,12], [16,20,17], [12,24,20,16], 
3790               [21,22,18], [12,23,21,18],
3791               [12,18,22,24],
3792            ]
3793        ] :
3794    tex=="hex_grid"?
3795        assert(is_undef(n), str(tex,__vnf_no_n_mesg))  
3796        assert(num_defined([gap, roughness])==0, str(tex," texture does not accept gap or roughness"))
3797        let(
3798            border=default(border,0.1)
3799        )
3800        assert(border>0 && border<0.5)
3801        let(
3802            diag=opp_ang_to_hyp(border,60),
3803            side=adj_ang_to_opp(1,30),
3804            hyp=adj_ang_to_hyp(0.5,30),
3805            sc = 1/3/hyp,
3806            hex=[ [1,2/6,0], [1/2,1/6,0], [0,2/6,0], [0,4/6,0], [1/2,5/6,0], [1,4/6,0] ]
3807        ) [
3808            [
3809                each hex,
3810                each move([0.5,0.5], p=yscale(sc, p=path3d(ellipse(d=1-2*border, circum=true, spin=-30,$fn=6),1))),
3811                hex[0]-[0,diag*sc,-1],
3812                for (ang=[270+60,270-60]) hex[1]+yscale(sc, p=cylindrical_to_xyz(diag,ang,1)),
3813                hex[2]-[0,diag*sc,-1],
3814                [0,0,1], [0.5-border,0,1], [0.5,0,0], [0.5+border,0,1], [1,0,1],
3815                hex[3]+[0,diag*sc,1],
3816                for (ang=[90+60,90-60]) hex[4]+yscale(sc, p=cylindrical_to_xyz(diag,ang,1)),
3817                hex[5]+[0,diag*sc,1],
3818                [0,1,1], [0.5-border,1,1], [0.5,1,0], [0.5+border,1,1], [1,1,1],
3819            ], [
3820                count(6,s=6),
3821                for (i=[0:1:5]) [i,(i+1)%6, (i+1)%6+6, i+6],
3822                [20,19,13,12], [17,16,15,14],
3823                [21,25,26,22], [23,28,29,24],
3824                [0,12,13,1], [1,14,15,2],
3825                [3,21,22,4], [4,23,24,5],
3826                [1,13,19,18], [1,18,17,14], 
3827                [4,22,26,27], [4,27,28,23],
3828            ]
3829        ] :
3830    tex=="rough"?
3831        assert(num_defined([gap,border])==0, str(tex," texture does not accept gap or border"))
3832        let(
3833            n = default(n,32),
3834            rough = default(roughness, 0.2)
3835        ) [
3836            for (y = [0:1:n-1])
3837            rands(0, rough, n, seed=123456+29*y)
3838        ] :
3839    assert(false, str("Unrecognized texture name: ", tex));
3840
3841
3842/// Function&Module: _textured_linear_sweep()
3843/// Usage: As Function
3844///   vnf = _textured_linear_sweep(region, texture, tex_size, h, ...);
3845///   vnf = _textured_linear_sweep(region, texture, counts=, h=, ...);
3846/// Usage: As Module
3847///   _textured_linear_sweep(region, texture, tex_size, h, ...) [ATTACHMENTS];
3848///   _textured_linear_sweep(region, texture, counts=, h=, ...) [ATTACHMENTS];
3849/// Topics: Sweep, Extrusion, Textures, Knurling
3850/// See Also: heightfield(), cylindrical_heightfield(), texture()
3851/// Description:
3852///   Given a [[Region|regions.scad]], creates a linear extrusion of it vertically, optionally twisted, scaled, and/or shifted,
3853///   with a given texture tiled evenly over the side surfaces.  The texture can be given in one of three ways:
3854///   - As a texture name string. (See {{texture()}} for supported named textures.)
3855///   - As a 2D array of evenly spread height values. (AKA a heightfield.)
3856///   - As a VNF texture tile.  A VNF tile exactly defines a surface from `[0,0]` to `[1,1]`, with the Z coordinates
3857///     being the height of the texture point from the surface.  VNF tiles MUST be able to tile in both X and Y
3858///     directions with no gaps, with the front and back edges aligned exactly, and the left and right edges as well.
3859///   . 
3860///   One script to convert a grayscale image to a texture heightfield array in a .scad file can be found at:
3861///   https://raw.githubusercontent.com/BelfrySCAD/BOSL2/master/scripts/img2scad.py
3862/// Arguments:
3863///   region = The [[Region|regions.scad]] to sweep/extrude.
3864///   texture = A texture name string, or a rectangular array of scalar height values (0.0 to 1.0), or a VNF tile that defines the texture to apply to vertical surfaces.  See {{texture()}} for what named textures are supported.
3865///   tex_size = An optional 2D target size for the textures.  Actual texture sizes will be scaled somewhat to evenly fit the available surface. Default: `[5,5]`
3866///   h / l = The height to extrude/sweep the path.
3867///   ---
3868///   counts = If given instead of tex_size, gives the tile repetition counts for textures over the surface length and height.
3869///   inset = If numeric, lowers the texture into the surface by that amount, before the tex_scale multiplier is applied.  If `true`, insets by exactly `1`.  Default: `false`
3870///   rot = If true, rotates the texture 90º.
3871///   tex_scale = Scaling multiplier for the texture depth.
3872///   twist = Degrees of twist for the top of the extrustion/sweep, compared to the bottom.  Default: 0
3873///   scale = Scaling multiplier for the top of the extrustion/sweep, compared to the bottom.  Default: 1
3874///   shift = [X,Y] amount to translate the top, relative to the bottom.  Default: [0,0]
3875///   style = The triangulation style used.  See {{vnf_vertex_array()}} for valid styles.  Used only with heightfield type textures. Default: `"min_edge"`
3876///   samples = Minimum number of "bend points" to have in VNF texture tiles.  Default: 8
3877///   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `CENTER`
3878///   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
3879///   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#subsection-orient).  Default: `UP`
3880/// Named Anchors:
3881///   "centroid_top" = The centroid of the top of the shape, oriented UP.
3882///   "centroid" = The centroid of the center of the shape, oriented UP.
3883///   "centroid_bot" = The centroid of the bottom of the shape, oriented DOWN.
3884
3885function _get_vnf_tile_edges(texture) =
3886    let(
3887        verts = texture[0],
3888        faces = texture[1],
3889        everts = [for (v = verts) (v.x==0 || v.y==0 || v.x==1 || v.y==1)],
3890        uc = unique_count([
3891            for (face = faces, i = idx(face))
3892            let(edge = select(face,i,i+1), i1 = min(edge), i2 = max(edge))
3893            if (everts[i1] && everts[i2])
3894            [i1, i2]
3895        ]),
3896        edges = uc[0], counts = uc[1],
3897        uedges = [for (i = idx(edges)) if (counts[i] == 1) edges[i] ]
3898    ) uedges;
3899
3900
3901function _validate_texture(texture) =
3902    is_vnf(texture)
3903      ? let( // Validate VNF tile texture
3904            bounds = pointlist_bounds(texture[0]),
3905            min_xy = point2d(bounds[0]),
3906            max_xy = point2d(bounds[1])
3907        )
3908        //assert(min_xy==[0,0] && max_xy==[1,1],"VNF tiles must span exactly from [0,0] to [1,1] in the X and Y components."))
3909        assert(all_nonnegative(concat(min_xy,[1,1]-max_xy)), "VNF tile X and Y components must be between 0 and 1.")
3910        let(
3911            verts = texture[0],
3912            uedges = _get_vnf_tile_edges(texture),
3913            edge_verts = [for (i = unique(flatten(uedges))) verts[i] ],
3914            hverts = [for(v = edge_verts) if(v.x==0 || v.x==1) v],
3915            vverts = [for(v = edge_verts) if(v.y==0 || v.y==1) v],
3916            allgoodx = all(hverts, function(v) any(hverts, function(w) approx(w,[1-v.x, v.y, v.z]))),
3917            allgoody = all(vverts, function(v) any(vverts, function(w) approx(w,[v.x, 1-v.y, v.z])))
3918        )
3919        assert(allgoodx && allgoody, "All VNF tile edge vertices must line up with a vertex on the opposite side of the tile.")
3920        true
3921      : // Validate heightfield texture.
3922        assert(is_matrix(texture), "Malformed texture.")
3923        let( tex_dim = list_shape(texture) )
3924        assert(len(tex_dim) == 2, "Heightfield texture must be a 2D square array of scalar heights.")
3925        assert(all_defined(tex_dim), "Heightfield texture must be a 2D square array of scalar heights.")
3926        true;
3927
3928
3929function _textured_linear_sweep(
3930    region, texture, tex_size=[5,5],
3931    h, counts, inset=false, rot=0,
3932    tex_scale=1, twist, scale, shift,
3933    style="min_edge", l, caps=true, 
3934    height, length, samples,
3935    anchor=CENTER, spin=0, orient=UP
3936) =
3937    assert(is_path(region,[2]) || is_region(region))
3938    assert(is_undef(samples) || is_int(samples))
3939    assert(counts==undef || is_vector(counts,2))
3940    assert(tex_size==undef || is_vector(tex_size,2))
3941    assert(is_bool(rot) || in_list(rot,[0,90,180,270]))
3942    assert(is_bool(caps) || is_bool_list(caps,2))
3943    let(
3944        caps = is_bool(caps) ? [caps,caps] : caps,
3945        regions = is_path(region,2)? [[region]] : region_parts(region),
3946        tex = is_string(texture)? texture(texture,$fn=_tex_fn_default()) : texture,
3947        dummy = assert(is_undef(samples) || is_vnf(tex), "You gave the tex_samples argument with a heightfield texture, which is not permitted.  Use the n= argument to texture() instead"),
3948        dummy2=is_bool(rot)?echo("boolean value for tex_rot is deprecated.  Use a numerical angle, one of 0, 90, 180, or 270.")0:0,
3949        texture = !rot? tex :
3950            is_vnf(tex)? zrot(is_num(rot)?rot:90, cp=[1/2,1/2], p=tex) :
3951            rot==180? reverse([for (row=tex) reverse(row)]) :
3952            rot==270? [for (row=transpose(tex)) reverse(row)] :
3953            reverse(transpose(tex)),
3954        h = first_defined([h, l, height, length, 1]),
3955        inset = is_num(inset)? inset : inset? 1 : 0,
3956        twist = default(twist, 0),
3957        shift = default(shift, [0,0]),
3958        scale = scale==undef? [1,1,1] :
3959            is_num(scale)? [scale,scale,1] : scale,
3960        samples = !is_vnf(texture)? len(texture[0]) :
3961            is_num(samples)? samples : 8,
3962        check_tex = _validate_texture(texture),
3963        sorted_tile =
3964            !is_vnf(texture)? texture :
3965            let(
3966                s = 1 / max(1, samples),
3967                vnf = samples<=1? texture :
3968                    let(
3969                        slice_us = list([s:s:1-s/2]),
3970                        vnft1 = vnf_slice(texture, "X", slice_us),
3971                        vnft = twist? vnf_slice(vnft1, "Y", slice_us) : vnft1,
3972                        zvnf = [
3973                            [
3974                                for (p=vnft[0]) [
3975                                    approx(p.x,0)? 0 : approx(p.x,1)? 1 : p.x,
3976                                    approx(p.y,0)? 0 : approx(p.y,1)? 1 : p.y,
3977                                    p.z
3978                                ]
3979                            ],
3980                            vnft[1]
3981                        ]
3982                    ) zvnf
3983            ) _vnf_sort_vertices(vnf, idx=[1,0]),
3984        vertzs = !is_vnf(sorted_tile)? undef :
3985            group_sort(sorted_tile[0], idx=1),
3986        tpath = is_vnf(sorted_tile)
3987            ? _find_vnf_tile_edge_path(sorted_tile,0)
3988            : let(
3989                  row = sorted_tile[0],
3990                  rlen = len(row)
3991              ) [for (i = [0:1:rlen]) [i/rlen, row[i%rlen]]],
3992        tmat = scale(scale) * zrot(twist) * up(h/2),
3993        pre_skew_vnf = vnf_join([
3994            for (rgn = regions) let(
3995                walls_vnf = vnf_join([
3996                    for (path = rgn) let(
3997                        path = reverse(path),
3998                        plen = path_length(path, closed=true),
3999                        counts = is_vector(counts,2)? counts :
4000                            is_vector(tex_size,2)
4001                              ? [round(plen/tex_size.x), max(1,round(h/tex_size.y)), ]
4002                              : [ceil(6*plen/h), 6],
4003                        obases = resample_path(path, n=counts.x * samples, closed=true),
4004                        onorms = path_normals(obases, closed=true),
4005                        bases = list_wrap(obases),
4006                        norms = list_wrap(onorms),
4007                        vnf = is_vnf(texture)
4008                          ? vnf_join( // VNF tile texture
4009                                let(
4010                                    row_vnf = vnf_join([
4011                                        for (i = [0:1:(scale==1?0:counts.y-1)], j = [0:1:counts.x-1]) [
4012                                            [
4013                                                for (group = vertzs)
4014                                                each [
4015                                                    for (vert = group) let(
4016                                                        u = floor((j + vert.x) * samples),
4017                                                        uu = ((j + vert.x) * samples) - u,
4018                                                        texh = tex_scale<0 ? -(1-vert.z - inset) * tex_scale
4019                                                                           : (vert.z - inset) * tex_scale,
4020                                                        base = lerp(bases[u], select(bases,u+1), uu),
4021                                                        norm = unit(lerp(norms[u], select(norms,u+1), uu)),
4022                                                        xy = base + norm * texh,
4023                                                        pt = point3d(xy,vert.y),
4024                                                        v = vert.y / counts.y,
4025                                                        vv = i / counts.y,
4026                                                        sc = lerp([1,1,1], scale, vv+v),
4027                                                        mat =
4028                                                            up((vv-0.5)*h) *
4029                                                            scale(sc) *
4030                                                            zrot(twist*(v+vv)) *
4031                                                            zscale(h/counts.y)
4032                                                    ) apply(mat, pt)
4033                                                ]
4034                                            ],
4035                                            sorted_tile[1]
4036                                        ]
4037                                    ])
4038                                ) [
4039                                    for (i = [0:1:0*(scale!=1?0:counts.y-1)])
4040                                    let(
4041                                        v = i / (scale==1?counts.y:1),
4042                                        sc = lerp([1,1,1], scale, v),
4043                                        mat =
4044                                            up((v)*h) *
4045                                            scale(sc) *
4046                                            zrot(twist*v)
4047                                    )
4048                                    apply(mat, row_vnf)
4049                                ]
4050                            )
4051                          : let( // Heightfield texture
4052                                texcnt = [len(texture[0]), len(texture)],
4053                                tile_rows = [
4054                                    for (ti = [0:1:texcnt.y-1])
4055                                    path3d([
4056                                        for (j = [0:1:counts.x])
4057                                        for (tj = [0:1:texcnt.x-1])
4058                                        if (j != counts.x || tj == 0)
4059                                        let(
4060                                            part = (j + (tj/texcnt.x)) * samples,
4061                                            u = floor(part),
4062                                            uu = part - u,
4063                                            texh = tex_scale<0 ? -(1-texture[ti][tj] - inset) * tex_scale
4064                                                               : (texture[ti][tj] - inset) * tex_scale,
4065                                            base = lerp(bases[u], select(bases,u+1), uu),
4066                                            norm = unit(lerp(norms[u], select(norms,u+1), uu)),
4067                                            xy = base + norm * texh
4068                                        ) xy
4069                                    ])
4070                                ],
4071                                tiles = [
4072                                    for (i = [0:1:counts.y], ti = [0:1:texcnt.y-1])
4073                                    if (i != counts.y || ti == 0)
4074                                    let(
4075                                        v = (i + (ti/texcnt.y)) / counts.y,
4076                                        sc = lerp([1, 1, 1], scale, v),
4077                                        mat = up((v-0.5)*h) *
4078                                              scale(sc) *
4079                                              zrot(twist*v)
4080                                    ) apply(mat, tile_rows[(texcnt.y-ti)%texcnt.y])
4081                                ]
4082                            ) vnf_vertex_array(
4083                                tiles, caps=false, style=style,
4084                                col_wrap=true, row_wrap=false,
4085                                reverse=true
4086                            )
4087                    ) vnf
4088                ]),
4089                brgn = [
4090                    for (path = rgn) let(
4091                        path = reverse(path),
4092                        plen = path_length(path, closed=true),
4093                        counts = is_vector(counts,2)? counts :
4094                            is_vector(tex_size,2)
4095                              ? [round(plen/tex_size.x), max(1,round(h/tex_size.y)), ]
4096                              : [ceil(6*plen/h), 6],
4097                        obases = resample_path(path, n=counts.x * samples, closed=true),
4098                        onorms = path_normals(obases, closed=true),
4099                        bases = list_wrap(obases),
4100                        norms = list_wrap(onorms),
4101                        nupath = [
4102                            for (j = [0:1:counts.x-1], vert = tpath) let(
4103                                part = (j + vert.x) * samples,
4104                                u = floor(part),
4105                                uu = part - u,
4106                                texh = tex_scale<0 ? -(1-vert.y - inset) * tex_scale
4107                                                   : (vert.y - inset) * tex_scale,
4108                                base = lerp(bases[u], select(bases,u+1), uu),
4109                                norm = unit(lerp(norms[u], select(norms,u+1), uu)),
4110                                xy = base + norm * texh
4111                            ) xy
4112                        ]
4113                    ) nupath
4114                ],
4115                bot_vnf = !caps[0] || brgn==[[]] ? EMPTY_VNF
4116                    : vnf_from_region(brgn, down(h/2), reverse=true),
4117                top_vnf = !caps[1] || brgn==[[]] ? EMPTY_VNF
4118                    : vnf_from_region(brgn, tmat, reverse=false)
4119            ) vnf_join([walls_vnf, bot_vnf, top_vnf])
4120        ]),
4121        skmat = down(h/2) * skew(sxz=shift.x/h, syz=shift.y/h) * up(h/2),
4122        final_vnf = apply(skmat, pre_skew_vnf),
4123        cent = centroid(region),
4124        anchors = [
4125            named_anchor("centroid_top", point3d(cent, h/2), UP),
4126            named_anchor("centroid",     point3d(cent),      UP),
4127            named_anchor("centroid_bot", point3d(cent,-h/2), DOWN)
4128        ]
4129    ) reorient(anchor,spin,orient, vnf=final_vnf, extent=true, anchors=anchors, p=final_vnf);
4130
4131
4132
4133function _find_vnf_tile_edge_path(vnf, val) =
4134    let(
4135        verts = vnf[0],
4136        fragments = [
4137            for(edge = _get_vnf_tile_edges(vnf))
4138            let(v0 = verts[edge[0]], v1 = verts[edge[1]])
4139            if (approx(v0.y, val) && approx(v1.y, val))
4140            v0.x <= v1.x? [[v0.x,v0.z], [v1.x,v1.z]] :
4141            [[v1.x,v1.z], [v0.x,v0.z]]
4142        ],
4143        sfrags = sort(fragments, idx=[0,1]),
4144        rpath = _assemble_a_path_from_fragments(sfrags)[0],
4145        opath = rpath==[]? []
4146              : rpath[0].x > last(rpath).x ? reverse(rpath)
4147              : rpath
4148    ) opath;
4149
4150
4151/// Function&Module: _textured_revolution()
4152/// Usage: As Function
4153///   vnf = _textured_revolution(shape, texture, tex_size, [tex_scale=], ...);
4154///   vnf = _textured_revolution(shape, texture, counts=, [tex_scale=], ...);
4155/// Usage: As Module
4156///   _textured_revolution(shape, texture, tex_size, [tex_scale=], ...) [ATTACHMENTS];
4157///   _textured_revolution(shape, texture, counts=, [tex_scale=], ...) [ATTACHMENTS];
4158/// Topics: Sweep, Extrusion, Textures, Knurling
4159/// See Also: heightfield(), cylindrical_heightfield(), texture()
4160/// Description:
4161///   Given a 2D region or path, fully in the X+ half-plane, revolves that shape around the Z axis (after rotating its Y+ to Z+).
4162///   This creates a solid from that surface of revolution, possibly capped top and bottom, with the sides covered in a given tiled texture.
4163///   The texture can be given in one of three ways:
4164///   - As a texture name string. (See {{texture()}} for supported named textures.)
4165///   - As a 2D array of evenly spread height values. (AKA a heightfield.)
4166///   - As a VNF texture tile.  A VNF tile exactly defines a surface from `[0,0]` to `[1,1]`, with the Z coordinates
4167///     being the height of the texture point from the surface.  VNF tiles MUST be able to tile in both X and Y
4168///     directions with no gaps, with the front and back edges aligned exactly, and the left and right edges as well.
4169///   .
4170///   One script to convert a grayscale image to a texture heightfield array in a .scad file can be found at:
4171///   https://raw.githubusercontent.com/BelfrySCAD/BOSL2/master/scripts/img2scad.py
4172/// Arguments:
4173///   shape = The path or region to sweep/extrude.
4174///   texture = A texture name string, or a rectangular array of scalar height values (0.0 to 1.0), or a VNF tile that defines the texture to apply to the revolution surface.  See {{texture()}} for what named textures are supported.
4175///   tex_size = An optional 2D target size for the textures.  Actual texture sizes will be scaled somewhat to evenly fit the available surface. Default: `[5,5]`
4176///   tex_scale = Scaling multiplier for the texture depth.
4177///   ---
4178///   inset = If numeric, lowers the texture into the surface by that amount, before the tex_scale multiplier is applied.  If `true`, insets by exactly `1`.  Default: `false`
4179///   rot = If true, rotates the texture 90º.
4180///   shift = [X,Y] amount to translate the top, relative to the bottom.  Default: [0,0]
4181///   closed = If false, and shape is given as a path, then the revolved path will be sealed to the axis of rotation with untextured caps.  Default: `true`
4182///   taper = If given, and `closed=false`, tapers the texture height to zero over the first and last given percentage of the path.  If given as a lookup table with indices between 0 and 100, uses the percentage lookup table to ramp the texture heights.  Default: `undef` (no taper)
4183///   angle = The number of degrees counter-clockwise from X+ to revolve around the Z axis.  Default: `360`
4184///   style = The triangulation style used.  See {{vnf_vertex_array()}} for valid styles.  Used only with heightfield type textures. Default: `"min_edge"`
4185///   counts = If given instead of tex_size, gives the tile repetition counts for textures over the surface length and height.
4186///   samples = Minimum number of "bend points" to have in VNF texture tiles.  Default: 8
4187///   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `CENTER`
4188///   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
4189///   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#subsection-orient).  Default: `UP`
4190/// Anchor Types:
4191///   "hull" = Anchors to the virtual convex hull of the shape.
4192///   "intersect" = Anchors to the surface of the shape.
4193
4194function _textured_revolution(
4195    shape, texture, tex_size, tex_scale=1,
4196    inset=false, rot=false, shift=[0,0],
4197    taper, closed=true, angle=360,
4198    inhibit_y_slicing=false,
4199    counts, samples,
4200    style="min_edge", atype="intersect",
4201    anchor=CENTER, spin=0, orient=UP
4202) = 
4203    assert(angle>0 && angle<=360)
4204    assert(is_path(shape,[2]) || is_region(shape))
4205    assert(is_undef(samples) || is_int(samples))
4206    assert(is_bool(closed))
4207    assert(counts==undef || is_vector(counts,2))
4208    assert(tex_size==undef || is_vector(tex_size,2))
4209    assert(is_bool(rot) || in_list(rot,[0,90,180,270]))
4210    let( taper_is_ok = is_undef(taper) || (is_finite(taper) && taper>=0 && taper<50) || is_path(taper,2) )
4211    assert(taper_is_ok, "Bad taper= value.")
4212    assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"")
4213    let(
4214        regions = !is_path(shape,2)? region_parts(shape) :
4215            closed? region_parts([shape]) :
4216            let(
4217                clpoly = [[0,shape[0].y], each shape, [0,last(shape).y]],
4218                dpoly = deduplicate(clpoly),
4219                cwpoly = is_polygon_clockwise(dpoly) ? dpoly : reverse(dpoly)
4220            )
4221            [[ select(cwpoly,1,-2) ]],
4222        checks = [
4223            for (rgn=regions, path=rgn)
4224            assert(all(path, function(pt) pt.x>=0))
4225        ]
4226    )
4227    assert(closed || is_path(shape,2))
4228    let(
4229        tex = is_string(texture)? texture(texture,$fn=_tex_fn_default()) : texture,
4230        dummy = assert(is_undef(samples) || is_vnf(tex), "You gave the tex_samples argument with a heightfield texture, which is not permitted.  Use the n= argument to texture() instead"),
4231        dummy2=is_bool(rot)?echo("boolean value for tex_rot is deprecated.  Use a numerical angle, one of 0, 90, 180, or 270.")0:0,        
4232        texture = !rot? tex :
4233            is_vnf(tex)? zrot(is_num(rot)?rot:90, cp=[1/2,1/2], p=tex) :
4234            rot==180? reverse([for (row=tex) reverse(row)]) :
4235            rot==270? [for (row=transpose(tex)) reverse(row)] :
4236            reverse(transpose(tex)),
4237        check_tex = _validate_texture(texture),
4238        inset = is_num(inset)? inset : inset? 1 : 0,
4239        samples = !is_vnf(texture)? len(texture) :
4240            is_num(samples)? samples : 8,
4241        bounds = pointlist_bounds(flatten(flatten(regions))),
4242        maxx = bounds[1].x,
4243        miny = bounds[0].y,
4244        maxy = bounds[1].y,
4245        h = maxy - miny,
4246        circumf = 2 * PI * maxx,
4247        tile = !is_vnf(texture)? texture :
4248            let(
4249                utex = samples<=1? texture :
4250                    let(
4251                        s = 1 / samples,
4252                        slices = list([s : s : 1-s/2]),
4253                        vnfx = vnf_slice(texture, "X", slices),
4254                        vnfy = inhibit_y_slicing? vnfx : vnf_slice(vnfx, "Y", slices),
4255                        vnft = vnf_triangulate(vnfy),
4256                        zvnf = [
4257                            [
4258                                for (p=vnft[0]) [
4259                                    approx(p.x,0)? 0 : approx(p.x,1)? 1 : p.x,
4260                                    approx(p.y,0)? 0 : approx(p.y,1)? 1 : p.y,
4261                                    p.z
4262                                ]
4263                            ],
4264                            vnft[1]
4265                        ]
4266                    ) zvnf
4267            ) _vnf_sort_vertices(utex, idx=[0,1]),
4268        vertzs = is_vnf(texture)? group_sort(tile[0], idx=0) : undef,
4269        bpath = is_vnf(tile)
4270            ? _find_vnf_tile_edge_path(tile,1)
4271            : let(
4272                  row = tile[0],
4273                  rlen = len(row)
4274              ) [for (i = [0:1:rlen]) [i/rlen, row[i%rlen]]],
4275        counts_x = is_vector(counts,2)? counts.x :
4276            is_vector(tex_size,2)
4277              ? max(1,round(angle/360*circumf/tex_size.x))
4278              : ceil(6*angle/360*circumf/h),
4279        taper_lup = closed || is_undef(taper)? [[-1,1],[2,1]] :
4280            is_num(taper)? [[-1,0], [0,0], [taper/100+EPSILON,1], [1-taper/100-EPSILON,1], [1,0], [2,0]] :
4281            is_path(taper,2)? let(
4282                retaper = [
4283                    for (t=taper)
4284                    assert(t[0]>=0 && t[0]<=100, "taper lookup indices must be between 0 and 100 inclusive.")
4285                    [t[0]/100, t[1]]
4286                ],
4287                taperout = [[-1,retaper[0][1]], each retaper, [2,last(retaper)[1]]]
4288            ) taperout :
4289            assert(false, "Bad taper= argument value."),
4290        full_vnf = vnf_join([
4291            for (rgn = regions) let(
4292                rgn_wall_vnf = vnf_join([
4293                    for (path = rgn) let(
4294                        plen = path_length(path, closed=closed),
4295                        counts_y = is_vector(counts,2)? counts.y :
4296                            is_vector(tex_size,2)? max(1,round(plen/tex_size.y)) : 6,
4297                        obases = resample_path(path, n=counts_y * samples + (closed?0:1), closed=closed),
4298                        onorms = path_normals(obases, closed=closed),
4299                        rbases = closed? list_wrap(obases) : obases,
4300                        rnorms = closed? list_wrap(onorms) : onorms,
4301                        bases = xrot(90, p=path3d(rbases)),
4302                        norms = xrot(90, p=path3d(rnorms)),
4303                        vnf = is_vnf(texture)
4304                          ? vnf_join([ // VNF tile texture
4305                                for (j = [0:1:counts_y-1])
4306                                [
4307                                    [
4308                                        for (group = vertzs) each [
4309                                            for (vert = group) let(
4310                                                part = (j + (1-vert.y)) * samples,
4311                                                u = floor(part),
4312                                                uu = part - u,
4313                                                base = lerp(select(bases,u), select(bases,u+1), uu),
4314                                                norm = unit(lerp(select(norms,u), select(norms,u+1), uu)),
4315                                                tex_scale = tex_scale * lookup(part/samples/counts_y, taper_lup),
4316                                                texh = tex_scale<0 ? -(1-vert.z - inset) * tex_scale * (base.x / maxx)
4317                                                                   : (vert.z - inset) * tex_scale * (base.x / maxx),
4318                                                xyz = base - norm * texh
4319                                            ) zrot(vert.x*angle/counts_x, p=xyz)
4320                                        ]
4321                                    ],
4322                                    tile[1]
4323                                ]
4324                            ])
4325                          : let( // Heightfield texture
4326                                texcnt = [len(texture[0]), len(texture)],
4327                                tiles = transpose([
4328                                    for (j = [0,1], tj = [0:1:texcnt.x-1])
4329                                    if (j == 0 || tj == 0)
4330                                    let(
4331                                        v = (j + (tj/texcnt.x)) / counts_x,
4332                                        mat = zrot(v*angle)
4333                                    ) apply(mat, [
4334                                        for (i = [0:1:counts_y-(closed?1:0)], ti = [0:1:texcnt.y-1])
4335                                        if (i != counts_y || ti == 0)
4336                                        let(
4337                                            part = (i + (ti/texcnt.y)) * samples,
4338                                            u = floor(part),
4339                                            uu = part - u,
4340                                            base = lerp(bases[u], select(bases,u+1), uu),
4341                                            norm = unit(lerp(norms[u], select(norms,u+1), uu)),
4342                                            tex_scale = tex_scale * lookup(part/samples/counts_y, taper_lup),
4343                                            texh = tex_scale<0 ? -(1-texture[ti][tj] - inset) * tex_scale * (base.x / maxx)
4344                                                               : (texture[ti][tj] - inset) * tex_scale * (base.x / maxx),
4345                                            xyz = base - norm * texh
4346                                        ) xyz
4347                                    ])
4348                                ])
4349                            ) vnf_vertex_array(
4350                                tiles, caps=false, style=style,
4351                                col_wrap=false, row_wrap=closed
4352                            )
4353                    ) vnf
4354                ]),
4355                walls_vnf = vnf_join([
4356                    for (i = [0:1:counts_x-1])
4357                    zrot(i*angle/counts_x, rgn_wall_vnf)
4358                ]),
4359                endcap_vnf = angle == 360? EMPTY_VNF :
4360                    let(
4361                        cap_rgn = [
4362                            for (path = rgn) let(
4363                                plen = path_length(path, closed=closed),
4364                                counts_y = is_vector(counts,2)? counts.y :
4365                                    is_vector(tex_size,2)? max(1,round(plen/tex_size.y)) : 6,
4366                                obases = resample_path(path, n=counts_y * samples + (closed?0:1), closed=closed),
4367                                onorms = path_normals(obases, closed=closed),
4368                                bases = closed? list_wrap(obases) : obases,
4369                                norms = closed? list_wrap(onorms) : onorms,
4370                                ppath = is_vnf(texture)
4371                                  ? [ // VNF tile texture
4372                                        for (j = [0:1:counts_y-1])
4373                                        for (group = vertzs, vert = reverse(group))
4374                                        if (approx(vert.x, 0)) let(
4375                                            part = (j + (1 - vert.y)) * samples,
4376                                            u = floor(part),
4377                                            uu = part - u,
4378                                            base = lerp(select(bases,u), select(bases,u+1), uu),
4379                                            norm = unit(lerp(select(norms,u), select(norms,u+1), uu)),
4380                                            tex_scale = tex_scale * lookup(part/samples/counts_y, taper_lup),
4381                                            texh = tex_scale<0 ? -(1-vert.z - inset) * tex_scale * (base.x / maxx)
4382                                                               : (vert.z - inset) * tex_scale * (base.x / maxx),
4383                                            xyz = base - norm * texh
4384                                        ) xyz
4385                                    ]
4386                                  : let( // Heightfield texture
4387                                        texcnt = [len(texture[0]), len(texture)]
4388                                    ) [
4389                                        for (i = [0:1:counts_y-(closed?1:0)], ti = [0:1:texcnt.y-1])
4390                                        if (i != counts_y || ti == 0)
4391                                        let(
4392                                            part = (i + (ti/texcnt.y)) * samples,
4393                                            u = floor(part),
4394                                            uu = part - u,
4395                                            base = lerp(bases[u], select(bases,u+1), uu),
4396                                            norm = unit(lerp(norms[u], select(norms,u+1), uu)),
4397                                            tex_scale = tex_scale * lookup(part/samples/counts_y, taper_lup),
4398                                            texh = tex_scale<0 ? -(1-texture[ti][0] - inset) * tex_scale * (base.x / maxx)
4399                                                               : (texture[ti][0] - inset) * tex_scale * (base.x / maxx),
4400                                            xyz = base - norm * texh
4401                                        ) xyz
4402                                    ],
4403                                path = closed? ppath : [
4404                                    [0, ppath[0].y],
4405                                    each ppath,
4406                                    [0, last(ppath).y],
4407                                ]
4408                            ) deduplicate(path, closed=closed)
4409                        ],
4410                        vnf2 = vnf_from_region(cap_rgn, xrot(90), reverse=false),
4411                        vnf3 = vnf_from_region(cap_rgn, rot([90,0,angle]), reverse=true)
4412                    ) vnf_join([vnf2, vnf3]),
4413                allcaps_vnf = closed? EMPTY_VNF :
4414                    let(
4415                        plen = path_length(rgn[0], closed=closed),
4416                        counts_y = is_vector(counts,2)? counts.y :
4417                            is_vector(tex_size,2)? max(1,round(plen/tex_size.y)) : 6,
4418                        obases = resample_path(rgn[0], n=counts_y * samples + (closed?0:1), closed=closed),
4419                        onorms = path_normals(obases, closed=closed),
4420                        rbases = closed? list_wrap(obases) : obases,
4421                        rnorms = closed? list_wrap(onorms) : onorms,
4422                        bases = xrot(90, p=path3d(rbases)),
4423                        norms = xrot(90, p=path3d(rnorms)),
4424                        caps_vnf = vnf_join([
4425                            for (j = [-1,0]) let(
4426                                base = select(bases,j),
4427                                norm = unit(select(norms,j)),
4428                                ppath = [
4429                                    for (vert = bpath) let(
4430                                        uang = vert.x / counts_x,
4431                                        tex_scale = tex_scale * lookup([0,1][j+1], taper_lup),
4432                                        texh = tex_scale<0 ? -(1-vert.y - inset) * tex_scale * (base.x / maxx)
4433                                                           : (vert.y - inset) * tex_scale * (base.x / maxx),
4434                                        xyz = base - norm * texh
4435                                    ) zrot(angle*uang, p=xyz)
4436                                ],
4437                                pplen = len(ppath),
4438                                zed = j<0? max(column(ppath,2)) :
4439                                    min(column(ppath,2)),
4440                                slice_vnf = [
4441                                    [
4442                                        each ppath,
4443                                        [0, 0, zed],
4444                                    ], [
4445                                        for (i = [0:1:pplen-2])
4446                                            j<0? [pplen, i, (i+1)%pplen] :
4447                                            [pplen, (i+1)%pplen, i]
4448                                    ]
4449                                ],
4450                                cap_vnf = vnf_join([
4451                                    for (i = [0:1:counts_x-1])
4452                                        zrot(i*angle/counts_x, p=slice_vnf)
4453                                ])
4454                            ) cap_vnf
4455                        ])
4456                    ) caps_vnf
4457            ) vnf_join([walls_vnf, endcap_vnf, allcaps_vnf])
4458        ]),
4459        skmat = down(-miny) * skew(sxz=shift.x/h, syz=shift.y/h) * up(-miny),
4460        skvnf = apply(skmat, full_vnf),
4461        geom = atype=="intersect"
4462              ? attach_geom(vnf=skvnf, extent=false)
4463              : attach_geom(vnf=skvnf, extent=true)
4464    ) reorient(anchor,spin,orient, geom=geom, p=skvnf);
4465
4466
4467module _textured_revolution(
4468    shape, texture, tex_size, tex_scale=1,
4469    inset=false, rot=false, shift=[0,0],
4470    taper, closed=true, angle=360,
4471    style="min_edge", atype="intersect",
4472    inhibit_y_slicing=false,
4473    convexity=10, counts, samples,
4474    anchor=CENTER, spin=0, orient=UP
4475) {
4476    dummy = assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"");
4477    vnf = _textured_revolution(
4478        shape, texture, tex_size=tex_size,
4479        tex_scale=tex_scale, inset=inset, rot=rot,
4480        taper=taper, closed=closed, style=style,
4481        shift=shift, angle=angle,
4482        samples=samples, counts=counts,
4483        inhibit_y_slicing=inhibit_y_slicing
4484    );
4485    geom = atype=="intersect"
4486          ? attach_geom(vnf=vnf, extent=false)
4487          : attach_geom(vnf=vnf, extent=true);
4488    attachable(anchor,spin,orient, geom=geom) {
4489        vnf_polyhedron(vnf, convexity=convexity);
4490        children();
4491    }
4492}
4493
4494
4495
4496// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap