BOSL2/paths.scad

//////////////////////////////////////////////////////////////////////
// LibFile: paths.scad
//   A `path` is a list of points of the same dimensions, usually 2D or 3D, that can
//   be connected together to form a sequence of line segments or a polygon.
//   A `region` is a list of paths that represent polygons, and the functions
//   in this file work on paths and also 1-regions, which are regions
//   that include exactly one path.  When you pass a 1-region to a function, the default
//   value for `closed` is always `true` because regions represent polygons.  
//   Capabilities include computing length of paths, computing
//   path tangents and normals, resampling of paths, and cutting paths up into smaller paths.  
// Includes:
//   include <BOSL2/std.scad>
// FileGroup: Advanced Modeling
// FileSummary: Operations on paths: length, resampling, tangents, splitting into subpaths
// FileFootnotes: STD=Included in std.scad
//////////////////////////////////////////////////////////////////////

// Section: Utility Functions
// Definitions:
//   Point|Points = A list of numbers, also called a vector.  Usually has length 2 or 3 to represent points in the place on points in space.  
//   Pointlist|Pointlists|Point List|Point Lists = An unordered list of {{points}}.
//   Path|Paths = An ordered list of two or more {{points}} specifying a path through space.  Usually points are 2D.  
//   Polygon|Polygons = A {{path}}, usually 2D, that describes a polygon by asuming that the first and last point are connected.

// Function: is_path()
// Synopsis: Returns True if 'list' is a {{path}}.
// Topics: Paths
// See Also: is_region(), is_vnf()
// Usage:
//   is_path(list, [dim], [fast])
// Description:
//   Returns true if `list` is a {{path}}.  A path is a list of two or more numeric vectors (AKA {{points}}).
//   All vectors must of the same size, and must contain numbers that are not inf or nan.
//   By default the vectors in a path must be 2D or 3D.  Set the `dim` parameter to specify a list
//   of allowed dimensions, or set it to `undef` to allow any dimension.  (This function
//   returns `false` on 1-regions.)  
// Example:
//   bool1 = is_path([[3,4],[5,6]]);    // Returns true
//   bool2 = is_path([[3,4]]);          // Returns false
//   bool3 = is_path([[3,4],[4,5]],2);  // Returns true
//   bool4 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
//   bool5 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
//   bool6 = is_path([[3,4,5],undef,[4,5,6]]);  // Returns false
//   bool7 = is_path([[3,5],[undef,undef],[4,5]]);  // Returns false
//   bool8 = is_path([[3,4],[5,6],[5,3]]);     // Returns true
//   bool9 = is_path([3,4,5,6,7,8]);           // Returns false
//   bool10 = is_path([[3,4],[5,6]], dim=[2,3]);// Returns true
//   bool11 = is_path([[3,4],[5,6]], dim=[1,3]);// Returns false
//   bool12 = is_path([[3,4],"hello"], fast=true); // Returns true
//   bool13 = is_path([[3,4],[3,4,5]]);            // Returns false
//   bool14 = is_path([[1,2,3,4],[2,3,4,5]]);      // Returns false
//   bool15 = is_path([[1,2,3,4],[2,3,4,5]],undef);// Returns true
// Arguments:
//   list = list to check
//   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
//   fast = set to true for fast check that looks at only the first entry.  Default: false
function is_path(list, dim=[2,3], fast=false) =
    fast
    ?   is_list(list) && is_vector(list[0]) 
    :   is_matrix(list) 
        && len(list)>1 
        && len(list[0])>0
        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));

// Function: is_1region()
// Synopsis: Returns true if {{path}} is a {{region}} with one component.
// Topics: Paths, Regions
// See Also: force_path()
// Usage:
//   bool = is_1region(path, [name])
// Description:
//   If `path` is a {{region}} with one component (a single-{{path}} region, or 1-region) then returns true.  If path is a region with more components
//   then display an error message about the parameter `name` requiring a path or a single component region.  If the input
//   is not a region then return false.  This function helps path functions accept 1-regions.
// Arguments:
//   path = input to process
//   name = name of parameter to use in error message.  Default: "path"
function is_1region(path, name="path") = 
     !is_region(path)? false
    :assert(len(path)==1,str("\nParameter \"",name,"\" must be a path or singleton region, but is a multicomponent region."))
     true;


// Function: force_path()
// Synopsis: Checks that path is a region with one component.
// SynTags: Path
// Topics: Paths, Regions
// See Also: is_1region()
// Usage:
//   outpath = force_path(path, [name])
// Description:
//   If `path` is a {{region}} with one component (a single-{{path}} region, or 1-region) then returns that component as a path.
//   If `path` is a region with more components then displays an error message about the parameter
//   `name` requiring a path or a single component region.  If the input is not a region then
//   returns the input without any checks.  This function helps path functions accept 1-regions.
// Arguments:
//   path = input to process
//   name = name of parameter to use in error message.  Default: "path"
function force_path(path, name="path") =
   is_region(path) ?
       assert(len(path)==1, str("\nParameter \"",name,"\" must be a path or singleton region, but is a multicomponent region."))
       path[0]
   : path;


/// Internal Function: _path_select()
/// Usage:
///   _path_select(path,s1,u1,s2,u2,[closed]):
/// Description:
///   Returns a portion of a path, from between the `u1` part of segment `s1`, to the `u2` part of
///   segment `s2`.  Both `u1` and `u2` are values between 0.0 and 1.0, inclusive, where 0 is the start
///   of the segment, and 1 is the end.  Both `s1` and `s2` are integers, where 0 is the first segment.
/// Arguments:
///   path = The path to get a section of.
///   s1 = The number of the starting segment.
///   u1 = The proportion along the starting segment, between 0.0 and 1.0, inclusive.
///   s2 = The number of the ending segment.
///   u2 = The proportion along the ending segment, between 0.0 and 1.0, inclusive.
///   closed = If true, treat path as a closed polygon.
function _path_select(path, s1, u1, s2, u2, closed=false) =
    let(
        lp = len(path),
        l = lp-(closed?0:1),
        u1 = s1<0? 0 : s1>l? 1 : u1,
        u2 = s2<0? 0 : s2>l? 1 : u2,
        s1 = constrain(s1,0,l),
        s2 = constrain(s2,0,l),
        pathout = concat(
            (s1<l && u1<1)? [lerp(path[s1],path[(s1+1)%lp],u1)] : [],
            [for (i=[s1+1:1:s2]) path[i]],
            (s2<l && u2>0)? [lerp(path[s2],path[(s2+1)%lp],u2)] : []
        )
    ) pathout;



// Function: path_merge_collinear()
// Synopsis: Removes unnecessary points from a path.
// SynTags: Path
// Topics: Paths, Regions
// Description:
//   Takes a {{path}} and removes unnecessary sequential collinear {{points}}.  When `closed=true` either of the path
//   endpoints may be removed.  
// Usage:
//   path_merge_collinear(path, [eps])
// Arguments:
//   path = A path of any dimension or a 1-region
//   closed = treat as closed polygon.  Default: false
//   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
function path_merge_collinear(path, closed, eps=EPSILON) =
    is_1region(path) ? path_merge_collinear(path[0], default(closed,true), eps) :
    let(closed=default(closed,false))
    assert(is_bool(closed))
    assert( is_path(path), "\nInvalid path in path_merge_collinear.")
    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "\nInvalid tolerance.")
    len(path)<=2 ? path :
    let(path = deduplicate(path, closed=closed))
    [
      if(!closed) path[0],
      for(triple=triplet(path,wrap=closed))
        if (!is_collinear(triple,eps=eps)) triple[1],
      if(!closed) last(path)
    ];


// Section: Path length calculation


// Function: path_length()
// Synopsis: Returns the path length.
// Topics: Paths
// See Also: path_segment_lengths(), path_length_fractions()
// Usage:
//   path_length(path,[closed])
// Description:
//   Returns the length of the given {{path}}.
// Arguments:
//   path = Path of any dimension or 1-region. 
//   closed = true if the path is closed.  Default: false
// Example:
//   path = [[0,0], [5,35], [60,-25], [80,0]];
//   echo(path_length(path));
function path_length(path,closed) =
    is_1region(path) ? path_length(path[0], default(closed,true)) :
    assert(is_path(path), "\nInvalid path in path_length.")
    let(closed=default(closed,false))
    assert(is_bool(closed))
    len(path)<2? 0 :
    sum([for (i = [0:1:len(path)-2]) norm(path[i+1]-path[i])])+(closed?norm(path[len(path)-1]-path[0]):0);


// Function: path_segment_lengths()
// Synopsis: Returns a list of the lengths of segments in a {{path}}.
// Topics: Paths
// See Also: path_length(), path_length_fractions()
// Usage:
//   path_segment_lengths(path,[closed])
// Description:
//   Returns list of the length of each segment in a path
// Arguments:
//   path = path in any dimension or 1-region
//   closed = true if the path is closed.  Default: false
function path_segment_lengths(path, closed) =
    is_1region(path) ? path_segment_lengths(path[0], default(closed,true)) :
    let(closed=default(closed,false))
    assert(is_path(path),"\nInvalid path in path_segment_lengths.")
    assert(is_bool(closed))
    [
        for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
        if (closed) norm(path[0]-last(path))
    ]; 


// Function: path_length_fractions()
// Synopsis: Returns the fractional distance of each point along the length of a path.
// Topics: Paths
// See Also: path_length(), path_segment_lengths()
// Usage:
//   fracs = path_length_fractions(path, [closed]);
// Description:
//    Returns the distance fraction of each point in the {{path}} along the path, so the first
//    point is zero and the final point is 1.  If the path is closed, the length of the output
//    will have one extra point because of the final connecting segment that connects the last
//    point of the path to the first point.
// Arguments:
//    path = path in any dimension or a 1-region
//    closed = set to true if path is closed.  Default: false
function path_length_fractions(path, closed) =
    is_1region(path) ? path_length_fractions(path[0], default(closed,true)):
    let(closed=default(closed, false))
    assert(is_path(path))
    assert(is_bool(closed))
    let(
        lengths = [
            0,
            each path_segment_lengths(path,closed)
        ],
        partial_len = cumsum(lengths),
        total_len = last(partial_len)
    )
    partial_len / total_len;



/// Internal Function: _path_self_intersections()
/// Usage:
///   isects = _path_self_intersections(path, [closed], [eps]);
/// Description:
///   Locates all self intersection {{points}} of the given {{path}}.  Returns a list of intersections, where
///   each intersection is a list like [POINT, SEGNUM1, PROPORTION1, SEGNUM2, PROPORTION2] where
///   POINT is the coordinates of the intersection point, SEGNUMs are the integer indices of the
///   intersecting segments along the path, and the PROPORTIONS are the 0.0 to 1.0 proportions
///   of how far along those segments they intersect at.  A proportion of 0.0 indicates the start
///   of the segment, and a proportion of 1.0 indicates the end of the segment.
///   .
///   This function does not return self-intersecting segments, only the points
///   where non-parallel segments intersect.  
/// Arguments:
///   path = The path to find self intersections of.
///   closed = If true, treat path like a closed polygon.  Default: true
///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
/// Example(2D):
///   path = [
///       [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100]
///   ];
///   isects = _path_self_intersections(path, closed=true);
///   // isects == [[[-33.3333, 0], 0, 0.666667, 4, 0.333333], [[33.3333, 0], 1, 0.333333, 3, 0.666667]]
///   stroke(path, closed=true, width=1);
///   for (isect=isects) translate(isect[0]) color("blue") sphere(d=10);
function _path_self_intersections(path, closed=true, eps=EPSILON) =
    let(
        path = closed ? list_wrap(path,eps=eps) : path,
        plen = len(path)
    )
    [ for (i = [0:1:plen-3]) let(
          a1 = path[i],
          a2 = path[i+1], 
          seg_normal = unit([-(a2-a1).y, (a2-a1).x],[0,0]),
          vals = path*seg_normal,
          ref  = a1*seg_normal,
            // The value of vals[j]-ref is positive if vertex j is one one side of the
            // line [a1,a2] and negative on the other side. Only a segment with opposite
            // signs at its two vertices can have an intersection with segment
            // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than
            // eps and the sign of vals[j]-ref otherwise.  
          signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) 
                        abs(vals[j]-ref) <  eps ? 0 : sign(vals[j]-ref) ]
        )
        if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2]
        for(j=[i+2:1:plen-(i==0 && closed? 3: 2)])
            if( signals[j-i-2]*signals[j-i-1]<=0 ) let( // segm [b1,b2] intersects line [a1,a2]
                b1 = path[j],
                b2 = path[j+1],
                isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
            )
            if (isect 
                && isect[1]>=-eps
                && isect[1]<= 1+eps
                && isect[2]>= -eps 
                && isect[2]<= 1+eps)
                [isect[0], i, isect[1], j, isect[2]]
    ];



// Section: Resampling - changing the number of points in a path


// Input `data` is a list that sums to an integer. 
// Returns rounded version of input data so that every 
// entry is rounded to an integer and the sum is the same as
// that of the input.  Works by rounding an entry in the list
// and passing the rounding error forward to the next entry.
// This generally distributes the error in a uniform manner. 
function _sum_preserving_round(data, index=0) =
    index == len(data)-1 ? list_set(data, len(data)-1, round(data[len(data)-1])) :
    let(
        newval = round(data[index]),
        error = newval - data[index]
    ) _sum_preserving_round(
        list_set(data, [index,index+1], [newval, data[index+1]-error]),
        index+1
    );


// Function: subdivide_path()
// Synopsis: Subdivides a path to produce a more finely sampled path.
// SynTags: Path
// Topics: Paths, Path Subdivision
// See Also: subdivide_and_slice(), resample_path(), jittered_poly()
// Usage:
//   newpath = subdivide_path(path, n|refine=|maxlen=, [method=], [closed=], [exact=]);
// Description:
//   Takes a {{path}} as input (closed or open) and subdivides the path to produce a more
//   finely sampled path.  You control the subdivision process by using the `maxlen` arg
//   to specify a maximum segment length, or by specifying `n` or `refine`, which request
//   a certain {{point}} count in the output.
//   .
//   You can specify the point count using the `n` option, where
//   you give the number of points you want in the output, or you can use
//   the `refine` option, where you specify a resampling factor.  If `refine=3` then
//   the number of points would increase by a factor of three, so a four point square would
//   have 12 points after subdivision.  With point-count subdivision, the new points can be distributed
//   proportional to length (`method="length"`), which is the default, or they can be divided up evenly among all the path segments
//   (`method="segment"`).  If the extra points don't fit evenly on the path then the
//   algorithm attempts to distribute them as uniformly as possible, but the result may be uneven.
//   The `exact` option, which is true by default, requires that the final point count is
//   exactly as requested.  For example, if you subdivide a four point square and request `n=13` then one
//   edge will have an extra point compared to the others.  
//   If you set `exact=false` then the
//   algorithm favors uniformity and the output path may have a different number of
//   points than you requested, but the sampling is still uniform.  In our example of the
//   square with `n=13`, you get only 12 points output, with the same number of points on each edge.
//   .
//   The points are always distributed uniformly on each segment.  The `method="length"` option does
//   means that the number of points on a segment is based on its length, but the points are still
//   distributed uniformly on each segment, independent of the other segments.  
//   With the `"segment"` method you can also give `n` as a vector of counts.  This 
//   specifies the desired point count on each segment: with vector valued `n` the `subdivide_path`
//   function places `n[i]-1` points on segment `i`.  The reason for the -1 is to avoid
//   double counting the endpoints, which are shared by pairs of segments, so that for
//   a closed polygon the total number of points is sum(n).  With an open
//   path there is an extra point at the end, so the number of points is sum(n)+1.
//   .
//   If you use the `maxlen` option then you specify the maximum length segment allowed in the output.
//   Each segment is subdivided into the largest number of segments meeting your requirement.  As above,
//   the sampling is uniform on each segment, independent of the other segments.  With the `maxlen` option
//   you cannot specify `method` or `exact`.    
// Arguments:
//   path = path in any dimension or a 1-region
//   n = scalar total number of points desired or with `method="segment"` can be a vector requesting `n[i]-1` new points added to segment i.
//   ---
//   refine = increase total number of points by this factor (specify only one of n, refine and maxlen)
//   maxlen = maximum length segment in the output (specify only one of n, refine and maxlen)
//   closed = set to false if the path is open.  Default: True
//   exact = if true return exactly the requested number of points, possibly sacrificing uniformity.  If false, return uniform point sample that may not match the number of points requested.  (Not allowed with maxlen.) Default: true
//   method = One of `"length"` or `"segment"`.  If `"length"`, adds vertices in proportion to segment length, so short segments get fewer points.  If `"segment"`, add points evenly among the segments, so all segments get the same number of points.  (Not allowed with maxlen.) Default: `"length"`
// Example(2D):
//   mypath = subdivide_path(square([2,2],center=true), 12);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D):
//   mypath = subdivide_path(square([8,2],center=true), 12);
//   move_copies(mypath)circle(r=.2,$fn=32);
// Example(2D):
//   mypath = subdivide_path(square([8,2],center=true), 12, method="segment");
//   move_copies(mypath)circle(r=.2,$fn=32);
// Example(2D):
//   mypath = subdivide_path(square([2,2],center=true), 17, closed=false);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): Specifying different numbers of points on each segment
//   mypath = subdivide_path(hexagon(side=2), [2,3,4,5,6,7], method="segment");
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): Requested point total is 14 but 15 points output due to extra end point
//   mypath = subdivide_path(pentagon(side=2), [3,4,3,4], method="segment", closed=false);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): Since 17 is not divisible by 5, a completely uniform distribution is not possible. 
//   mypath = subdivide_path(pentagon(side=2), 17);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): With `exact=false` a uniform distribution, but only 15 points
//   mypath = subdivide_path(pentagon(side=2), 17, exact=false);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): With `exact=false` you can also get extra points, here 20 instead of requested 18
//   mypath = subdivide_path(pentagon(side=2), 18, exact=false);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): Using refine in this example multiplies the point count by 3 by adding 2 points to each edge.
//   mypath = subdivide_path(pentagon(side=2), refine=3);
//   move_copies(mypath)circle(r=.1,$fn=32);
// Example(2D): However, refine doesn't distribute evenly by segment unless you change the method.  With the default method set to `"length"`, the points are distributed with more on the long segments in this example using refine.  
//   mypath = subdivide_path(square([8,2],center=true), refine=3);
//   move_copies(mypath)circle(r=.2,$fn=32);
// Example(2D): In this example with maxlen, every side gets a different number of new points
//   path = [[0,0],[0,4],[10,6],[10,0]];
//   spath = subdivide_path(path, maxlen=2, closed=true);
//   move_copies(spath) circle(r=.25,$fn=12);
// Example(FlatSpin,VPD=15,VPT=[0,0,1.5]): Three-dimensional paths also work
//   mypath = subdivide_path([[0,0,0],[2,0,1],[2,3,2]], 12);
//   move_copies(mypath)sphere(r=.1,$fn=32);
function subdivide_path(path, n, refine, maxlen, closed=true, exact, method) =
    let(path = force_path(path))
    assert(is_path(path), "\nInvalid path or 1-region.")
    assert(num_defined([n,refine,maxlen]), "\nMust give exactly one of n, refine, and maxlen.")
    refine==1 || n==len(path) ? path :
    is_def(maxlen) ?
        assert(is_undef(method), "\nCannot give method with maxlen.")
        assert(is_undef(exact), "\nCannot give exact with maxlen.")
        [
         for (p=pair(path,closed))
           let(steps = ceil(norm(p[1]-p[0])/maxlen))
           each lerpn(p[0], p[1], steps, false),
         if (!closed) last(path)
        ]               
    :
    let(
        exact = default(exact, true),
        method = default(method, "length")
    )
    assert(method=="length" || method=="segment")
    let(
        n = !is_undef(n)? n :
            !is_undef(refine)? len(path) * refine :
            undef
    )
    assert((is_num(n) && n>0) || is_vector(n), "\nParameter n to subdivide_path must be postive number or vector.")
    let(
        count = len(path) - (closed?0:1), 
        add_guess = method=="segment"?
                       (
                          is_list(n)
                          ? assert(len(n)==count, "\nVector parameter n to subdivide_path has the wrong length.")
                            add_scalar(n,-1)
                          : repeat((n-len(path)) / count, count)
                       )
                  : // method=="length"
                    assert(is_num(n), "\nParameter n to subdivide path must be a number when method=\"length\".")
                    let(
                        path_lens = path_segment_lengths(path,closed),
                        add_density = (n - len(path)) / sum(path_lens)
                    )
                    path_lens * add_density,
        add = exact? _sum_preserving_round(add_guess)
                   : [for (val=add_guess) round(val)]
    )
    [
        for (i=[0:1:count-1]) 
           each lerpn(path[i],select(path,i+1), 1+add[i],endpoint=false),
        if (!closed) last(path)
    ];




// Function: resample_path()
// Synopsis: Returns an equidistant set of points along a path.
// SynTags: Path
// Topics: Paths
// See Also: simplify_path(), subdivide_path()
// Usage:
//   newpath = resample_path(path, n|spacing=, [closed=]);
// Description:
//   Compute a uniform resampling of the input {{path}}.  If you specify `n` then the output path has `n`
//   {{points}} spaced uniformly (by linear interpolation along the input path segments).  The only points of the
//   input path that are guaranteed to appear in the output path are the starting and ending points, and any
//   points that have an angular deflection of at least the number of degrees given in `keep_corners`.
//   If you specify `spacing` then the length you give is rounded to the nearest spacing that gives
//   a uniform sampling of the path and the resulting uniformly sampled path is returned.
//   .
//   Because this function operates on a discrete input path the quality of the output depends on
//   the sampling of the input.  If you want accurate output, use many points for the input.
// Arguments:
//   path = path in any dimension or a 1-region
//   n = Number of points in output
//   ---
//   spacing = Approximate spacing desired
//   keep_corners = If given a scalar, path vertices with deflection angle greater than this are preserved in the output.
//   closed = set to true if path is closed.  Default: true
// Example(2D):  Subsampling lots of points from a smooth curve
//   path = xscale(2,circle($fn=250, r=10));
//   sampled = resample_path(path, 16);
//   stroke(path);
//   color("red")move_copies(sampled) circle($fn=16);
// Example(2D): Specified spacing is rounded to make a uniform sampling
//   path = xscale(2,circle($fn=250, r=10));
//   sampled = resample_path(path, spacing=17);
//   stroke(path);
//   color("red")move_copies(sampled) circle($fn=16);
// Example(2D): Notice that the corners are excluded.
//   path = square(20);
//   sampled = resample_path(path, spacing=6);
//   stroke(path,closed=true);
//   color("red")move_copies(sampled) circle($fn=16);
// Example(2D): Forcing preservation of corners.
//   path = square(20);
//   sampled = resample_path(path, spacing=6, keep_corners=90);
//   stroke(path,closed=true);
//   color("red")move_copies(sampled) circle($fn=16);
// Example(2D): Closed set to false
//   path = square(20);
//   sampled = resample_path(path, spacing=6,closed=false);
//   stroke(path);
//   color("red")move_copies(sampled) circle($fn=16);

function resample_path(path, n, spacing, keep_corners, closed=true) =
    let(path = force_path(path))
    assert(is_path(path), "\nInvalid path or 1-region.")
    assert(num_defined([n,spacing])==1,"\nMust define exactly one of n and spacing.")
    assert(n==undef || (is_integer(n) && n>0))
    assert(spacing==undef || (is_finite(spacing) && spacing>0))
    assert(is_bool(closed))
    let(
        corners = is_undef(keep_corners)
          ? [0, len(path)-(closed?0:1)]
          : [
                0,
                for (i = [1:1:len(path)-(closed?1:2)])
                    let( ang = abs(modang(vector_angle(select(path,i-1,i+1))-180)) )
                    if (ang >= keep_corners) i,
                len(path)-(closed?0:1),
            ],
        pcnt = len(path),
        plen = path_length(path, closed=closed),
        subpaths = [ for (p = pair(corners)) [for(i = [p.x:1:p.y]) path[i%pcnt]] ],
        n = is_undef(n)? undef : closed? n+1 : n
    )
    assert(n==undef || n >= len(corners), "\nThere are nore than `n=` corners whose angle is greater than `keep_corners=`.")
    let(
        lens = [for (subpath = subpaths) path_length(subpath)],
        part_ns = is_undef(n)
          ? [for (i=idx(subpaths)) max(1,round(lens[i]/spacing)-1)]
          : let(
                ccnt = len(corners),
                parts = [for (l=lens) (n-ccnt) * l/plen]
            )
            _sum_preserving_round(parts),
        out = [
            for (i = idx(subpaths))
                let(
                    subpath = subpaths[i],
                    splen = lens[i],
                    pn = part_ns[i] + 1,
                    distlist = lerpn(0, splen, pn, false),
                    cuts = path_cut_points(subpath, distlist, closed=false)
                )
                each column(cuts,0),
            if (!closed) last(path)
        ]
    ) out;


// Function: simplify_path()
// Synopsis: Removes points from an irregular path, preserving dominant features.
// SynTags: Path
// Topics: Paths
// See Also: resample_path()
// Usage:
//   newpath = simplify_path(path, maxerr, [closed=]);
// Description:
//   This is intended for irregular paths such as coastlines, or paths having fractal self-similarity.
//   The original path is simplified by removing points that fall within a specified margin of error,
//   leaving behind those points that contribute to dominant features of the path. This operation has the
//   effect of making the point spacing somewhat more uniform. For coastlines, up to 80% reduction in path
//   length is possible with small degradation of the original shape. The input path may be 2D or 3D.
//   .
//   The `maxerr` parameter determines which points of the original path are kept. A point is kept if it
//   deviates beyond `maxerr` distance from a straight line between the last kept point and a point further
//   along the path. When a new deviation is found, that deviating point is kept and the process repeats from
//   that new kept point. The best value of `maxerr` depends on the geometry of the path and the amount of
//   reduction you want. A smaller value of `maxerr` returns more detail in the output, and a larger value
//   causes details to be lost. For paths such as coastlines, a `maxerr` value less than 1% of the maximum
//   bounding box dimension is a good starting value.
//   .
//   For unclosed paths (where `closed=false`) the endpoints of the path are preserved. When `closed=true`,
//   the path is treated as continuous and only dominant features that happen to be near the endpoints are
//   included in the output.
// Arguments:
//   path = Path in any dimension or 1-region
//   maxerr = Maximum deviation from line connecting last kept point to a further point; points beyond this deviation are kept.
//   ---
//   closed = Set to true if path is closed. If false, endpoints are retained in the output. Default: false
// Example(2D,Med,VPD=38000,VPT=[5600,6500,0]): A map of California, originally a 262-point polygon (yellow, on left), reduced to 39 points (green, on right).
//   calif = [
//   [225,12681], [199,12544], [180,12490], [221,12435], [300,12342], [310,12315], [320,12263], [350,12154],
//   [374,11968], [350,11820], [328,11707], [291,11586], [259,11553], [275,11499], [304,11420], [312,11321],
//   [273,11189], [233,11066], [200,10995], [160,10942], [104,10820], [0,10568], [25,10510], [50,10420],
//   [65,10312], [271,10108], [368,10004], [438,9909], [517,9809], [569,9741], [600,9666], [615,9600],
//   [630,9567], [649,9526], [679,9385], [670,9245], [650,9187], [635,9113], [644,8985], [673,8938], [694,8846],
//   [740,8745], [770,8678], [780,8635], [771,8528], [745,8449], [738,8403], [807,8364], [872,8298], [894,8264],
//   [1090,8076], [1270,7877], [1366,7798], [1440,7679], [1495,7596], [1543,7541], [1560,7487], [1575,7447],
//   [1576,7350], [1536,7234], [1521,7168], [1587,7184], [1761,7129], [1838,7050], [1893,7050], [1995,6995],
//   [2081,6940], [2109,7006], [2100,7045], [2100,7090], [2109,7155], [2115,7210], [2100,7269], [2124,7334],
//   [2179,7365], [2209,7391], [2242,7362], [2308,7311], [2280,7215], [2220,7164], [2210,7150], [2200,7095],
//   [2200,7040], [2234,7040], [2274,6932], [2415,6775], [2459,6691], [2483,6578], [2558,6497], [2610,6449],
//   [2598,6430], [2490,6475], [2444,6500], [2410,6515], [2406,6530], [2375,6570], [2305,6610], [2224,6638],
//   [2225,6806], [2211,6867], [2159,6913], [2109,6912], [2075,6810], [2074,6583], [2068,6521], [2104,6503],
//   [2140,6454], [2153,6417], [2184,6336], [2187,6243], [2173,6158], [2213,6065], [2250,6005], [2283,5970],
//   [2343,5928], [2370,5875], [2428,5822], [2485,5779], [2606,5782], [2728,5785], [2772,5725], [2850,5561],
//   [2839,5472], [2820,5391], [2797,5322], [2734,5321], [2676,5330], [2640,5289], [2656,5236], [2661,5205],
//   [2671,5144], [2712,5083], [2720,4973], [2738,4882], [2806,4819], [2891,4780], [2966,4737], [3004,4662],
//   [3043,4604], [3080,4542], [3128,4491], [3170,4453], [3294,4262], [3370,4150], [3384,4090], [3402,4057],
//   [3442,4029], [3602,3909], [3753,3777], [3855,3600], [3830,3521], [3900,3425], [3957,3394], [4000,3390],
//   [4045,3393], [4109,3315], [4121,3235], [4089,3125], [4074,3085], [4081,3019], [4098,2923], [4116,2848],
//   [4160,2774], [4135,2734], [4116,2697], [4100,2645], [4123,2585], [4208,2558], [4272,2478], [4367,2441],
//   [4453,2461], [4533,2485], [4635,2488], [4795,2472], [4875,2450], [4896,2425], [4933,2402], [4988,2404],
//   [5036,2409], [5124,2401], [5388,2338], [5479,2252], [5552,2199], [5610,2121], [5658,2051], [5802,1986],
//   [5858,1955], [5994,1930], [6110,1905], [6122,1880], [6174,1895], [6309,1910], [6447,1875], [6510,1797],
//   [6532,1730], [6525,1594], [6531,1521], [6615,1490], [6697,1525], [6729,1555], [6812,1543], [6901,1506],
//   [7084,1351], [7193,1305], [7250,1205], [7255,1155], [7305,1152], [7420,1087], [7522,995], [7751,608],
//   [7771,455], [7780,278], [7856,218], [7914,192], [7928,143], [7860,170], [7910,80], [7931,39], [7944,-2],
//   [8009,4], [8490,46], [10165,189], [10542,220], [10575,240], [10654,239], [10701,218], [10735,242],
//   [10776,308], [10811,372], [10834,527], [10781,598], [10685,625], [10595,763], [10588,913], [10553,965],
//   [10572,995], [10585,1062], [10617,1120], [10672,1153], [10755,1340], [10781,1436], [10811,1637],
//   [10791,1768], [10771,1807], [10824,1852], [10927,2015], [10995,2073], [11184,2212], [11204,2270],
//   [11174,2312], [11045,2430], [10931,2585], [10881,2678], [10806,2818], [10739,2936], [10670,3102],
//   [10670,3166], [4823,8540], [4804,12775], [4798,12800], [2515,12800], [232,12800]
//   ];
//   newpoly = simplify_path(calif, 120, closed=true);
//   left(4000) polygon(calif);
//   right(4000) color("lightgreen") polygon(newpoly);

function simplify_path(path, maxerr, closed=false) =
    let(path = force_path(path))
    assert(is_path(path), "\nInvalid path or 1-region.")
    assert(is_num(maxerr) && maxerr>0, "\nParameter 'maxerr' must be a positive number.")
    let(
        n = len(path),
        unclosed = _err_resample(path, maxerr, n) // get simplified path including original endpoints
    ) closed ? let( // search for new corners between the corners found on either side of the end points
        nu = len(unclosed),
        cornerpath = [
            for(i=[unclosed[nu-2]:n-1]) path[i],
            for(i=[0:unclosed[1]]) path[i]
        ],
        corner_resample = _err_resample(cornerpath, maxerr, len(cornerpath)),
        nc = len(corner_resample)
    ) [
        for(i=[1:nu-2]) path[unclosed[i]], // exclude endpoints
        if(nc>2) for(i=[1:nc-2]) cornerpath[corner_resample[i]] // insert new corners if any
    ]
    : [ for(i=unclosed) path[i] ];

/// return a resampled path based on error deviation, retaining path endpoints (i.e. assume path is not closed)
function _err_resample(path, maxerr, n, i1=0, i2=2, resultidx=[0], iter=0) =
    n <= 2 ? path :
    i2 >= n || i2-i1<2 ? concat(resultidx, [n-1]) : let(
        dists = [ for(i=[i1+1:i2-1]) let(j=i%n) point_line_distance(path[j], [path[i1], path[i2%n]]) ],
        imaxdist = max_index(dists),
        newfound = dists[imaxdist] >= maxerr,
        newidx1 = newfound ? i1+imaxdist+1 : i1,
        newidx2 = newfound ? min(newidx1+2, n) : min(i2+1,n)
    )
    _err_resample(path, maxerr, n, newidx1, newidx2, newfound ? concat(resultidx, [newidx1]) : resultidx, iter+1);



// Section: Path Geometry

// Function: is_path_simple()
// Synopsis: Returns true if a {{path}} has no self intersections.
// Topics: Paths
// See Also: is_path()
// Usage:
//   bool = is_path_simple(path, [closed], [eps]);
// Description:
//   Returns true if the given 2D {{path}} is simple, meaning that it has no self-intersections.
//   Repeated {{points}} are not considered self-intersections: a path with such points can
//   still be simple.  
//   If closed is set to true then treat the path as a polygon.
// Arguments:
//   path = 2D path or 1-region
//   closed = set to true to treat path as a polygon.  Default: false
//   eps = Epsilon error value used for determine if points coincide.  Default: `EPSILON` (1e-9)
function is_path_simple(path, closed, eps=EPSILON) =
    is_1region(path) ? is_path_simple(path[0], default(closed,true), eps) :
    let(closed=default(closed,false))
    assert(is_path(path, 2),"\nMust give a 2D path.")
    assert(is_bool(closed))
    let(
        path = deduplicate(path,closed=closed,eps=eps)
    )
    // check for path reversals
    [for(i=[0:1:len(path)-(closed?2:3)])
         let(v1=path[i+1]-path[i],
             v2=select(path,i+2)-path[i+1],
             normv1 = norm(v1),
             normv2 = norm(v2)
             )
         if (approx(v1*v2/normv1/normv2,-1)) 1
    ]  == [] 
    &&
    _path_self_intersections(path,closed=closed,eps=eps) == [];


// Function: path_closest_point()
// Synopsis: Returns the closest place on a {{path}} to a given {{point}}.
// Topics: Paths
// See Also: point_line_distance(), line_closest_point()
// Usage:
//   index_pt = path_closest_point(path, pt);
// Description:
//   Finds the closest {{path}} segment, and {{point}} on that segment to the given point.
//   Returns `[SEGNUM, POINT]`
// Arguments:
//   path = Path of any dimension or a 1-region.
//   pt = The point to find the closest point to.
//   closed = If true, the path is considered closed.
// Example(2D):
//   path = circle(d=100,$fn=6);
//   pt = [20,10];
//   closest = path_closest_point(path, pt);
//   stroke(path, closed=true);
//   color("blue") translate(pt) circle(d=3, $fn=12);
//   color("red") translate(closest[1]) circle(d=3, $fn=12);
function path_closest_point(path, pt, closed=true) =
    let(path = force_path(path))
    assert(is_path(path), "\nInput must be a path.")
    assert(is_vector(pt, len(path[0])), "\nInput pt must be a compatible vector.")
    assert(is_bool(closed))
    let(
        pts = [for (seg=pair(path,closed)) line_closest_point(seg,pt,SEGMENT)],
        dists = [for (p=pts) norm(p-pt)],
        min_seg = min_index(dists)
    ) [min_seg, pts[min_seg]];


// Function: path_tangents()
// Synopsis: Returns tangent vectors for each point along a path.
// Topics: Paths
// See Also: path_normals()
// Usage:
//   tangs = path_tangents(path, [closed], [uniform]);
// Description:
//   Compute the tangent vector to the input {{path}}.  The derivative approximation is described in deriv().
//   The returned vectors are normalized to length 1.  If any derivatives are zero then
//   the function fails with an error.  If you set `uniform` to false then the sampling is
//   assumed to be non-uniform and the derivative is computed with adjustments to produce corrected
//   values.
// Arguments:
//   path = path of any dimension or a 1-region
//   closed = set to true of the path is closed.  Default: false
//   uniform = set to false to correct for non-uniform sampling.  Default: true
// Example(2D): A shape with non-uniform sampling gives distorted derivatives that may be undesirable.  Derivatives tilt toward the long edges of the rectangle.  
//   rect = square([10,3]);
//   tangents = path_tangents(rect,closed=true);
//   stroke(rect,closed=true, width=0.25);
//   color("purple")
//       for(i=[0:len(tangents)-1])
//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
// Example(2D): Setting uniform to false corrects the distorted derivatives for this example:
//   rect = square([10,3]);
//   tangents = path_tangents(rect,closed=true,uniform=false);
//   stroke(rect,closed=true, width=0.25);
//   color("purple")
//       for(i=[0:len(tangents)-1])
//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
function path_tangents(path, closed, uniform=true) =
    is_1region(path) ? path_tangents(path[0], default(closed,true), uniform) :
    let(closed=default(closed,false))
    assert(is_bool(closed))
    assert(is_path(path))
    !uniform ? [for(t=deriv(path,closed=closed, h=path_segment_lengths(path,closed))) unit(t)]
             : [for(t=deriv(path,closed=closed)) unit(t)];


// Function: path_normals()
// Synopsis: Returns normal vectors for each point along a path.
// Topics: Paths
// See Also: path_tangents()
// Usage:
//   norms = path_normals(path, [tangents], [closed]);
// Description:
//   Compute the normal vector to the input {{path}}.  This vector is perpendicular to the
//   path tangent and lies in the plane of the curve.  For 3d paths we define the plane of the curve
//   at path {{point}} i to be the plane defined by point i and its two neighbors.  At the endpoints of open paths
//   we use the three end points.  For 3d paths the computed normal is the one lying in this plane that points
//   toward the center of curvature at that path point.  For 2D paths, which lie in the xy plane, the normal
//   is the path pointing to the right of the direction the path is traveling.  If points are collinear then
//   a 3d path has no center of curvature, and hence the 
//   normal is not uniquely defined.  In this case the function issues an error.
//   For 2D paths the plane is always defined so the normal fails to exist only
//   when the derivative is zero (in the case of repeated points).
// Arguments:
//   path = 2D or 3D path or a 1-region
//   tangents = path tangents optionally supplied
//   closed = if true path is treated as a polygon.  Default: false
function path_normals(path, tangents, closed) =
    is_1region(path) ? path_normals(path[0], tangents, default(closed,true)) :
    let(closed=default(closed,false))
    assert(is_path(path,[2,3]))
    assert(is_bool(closed))
    let(
         tangents = default(tangents, path_tangents(path,closed)),
         dim=len(path[0])
    )
    assert(is_path(tangents) && len(tangents[0])==dim,"\nDimensions of path and tangents must match.")
    [
     for(i=idx(path))
         let(
             pts = i==0 ? (closed? select(path,-1,1) : select(path,0,2))
                 : i==len(path)-1 ? (closed? select(path,i-1,i+1) : select(path,i-2,i))
                 : select(path,i-1,i+1)
        )
        dim == 2 ? [tangents[i].y,-tangents[i].x]
                 : let( v=cross(cross(pts[1]-pts[0], pts[2]-pts[0]),tangents[i]))
                   assert(norm(v)>EPSILON, "\n3D path contains collinear points.")
                   unit(v)
    ];


// Function: path_curvature()
// Synopsis: Returns the estimated numerical curvature of the {{path}}.
// Topics: Paths
// See Also: path_tangents(), path_normals(), path_torsion()
// Usage:
//   curvs = path_curvature(path, [closed]);
// Description:
//   Numerically estimate the curvature of the {{path}} (in any dimension).
// Arguments:
//   path = path in any dimension or a 1-region
//   closed = if true then treat the path as a polygon.  Default: false
function path_curvature(path, closed) =
    is_1region(path) ? path_curvature(path[0], default(closed,true)) :
    let(closed=default(closed,false))
    assert(is_bool(closed))
    assert(is_path(path))
    let( 
        d1 = deriv(path, closed=closed),
        d2 = deriv2(path, closed=closed)
    ) [
        for(i=idx(path))
        sqrt(
            sqr(norm(d1[i])*norm(d2[i])) -
            sqr(d1[i]*d2[i])
        ) / pow(norm(d1[i]),3)
    ];


// Function: path_torsion()
// Synopsis: Returns the estimated numerical torsion of the {{path}}.
// Topics: Paths
// See Also: path_tangents(), path_normals(), path_curvature()
// Usage:
//   torsions = path_torsion(path, [closed]);
// Description:
//   Numerically estimate the torsion of a 3d {{path}}.
// Arguments:
//   path = 3D path
//   closed = if true then treat path as a polygon.  Default: false
function path_torsion(path, closed=false) =
    assert(is_path(path,3), "\nInput path must be a 3d path.")
    assert(is_bool(closed))
    let(
        d1 = deriv(path,closed=closed),
        d2 = deriv2(path,closed=closed),
        d3 = deriv3(path,closed=closed)
    ) [
        for (i=idx(path)) let(
            crossterm = cross(d1[i],d2[i])
        ) crossterm * d3[i] / sqr(norm(crossterm))
    ];


// Function: surface_normals()
// Synopsis: Estimates the normals to a surface defined by a {{point}} array
// Topics: Math, Geometry
// See Also: path_tangents(), path_normals()
// Usage:
//   normals = surface_normals(surf, [col_wrap=], [row_wrap=]);
// Description:
//   Numerically estimate the normals to a surface defined by a 2D array of 3d {{points}}, which can
//   also be regarded as an array of {{paths}} (all of the same length).  
// Arguments:
//   surf = surface in 3d defined by a 2D array of points
//   ---
//   row_wrap = if true then wrap path in the row direction (first index)
//   col_wrap = if true then wrap path in the column direction (second index)

function surface_normals(surf, col_wrap=false, row_wrap=false) =
  let(
      rowderivs = [for(y=[0:1:len(surf)-1])  path_tangents(surf[y],closed=col_wrap)],
      colderivs = [for(x=[0:1:len(surf[0])-1]) path_tangents(column(surf,x), closed=row_wrap)]
  )
  [for(y=[0:1:len(surf)-1])
     [for(x=[0:1:len(surf[0])-1])
         cross(colderivs[x][y],rowderivs[y][x])]];



// Section: Breaking paths up into subpaths



// Function: path_cut()
// Synopsis: Cuts a {{path}} into subpaths at various {{points}}.
// SynTags: PathList
// Topics: Paths, Path Subdivision
// See Also: split_path_at_self_crossings(), path_cut_points()
// Usage:
//   path_list = path_cut(path, cutdist, [closed]);
// Description:
//   Given a list of distances in `cutdist`, cut the {{path}} into
//   subpaths at those lengths, returning a list of paths.
//   If the input path is closed then the final path includes the
//   original starting {{point}}.  The list of cut distances must be
//   in ascending order and should not include the endpoints: 0 
//   or `len(path)`.  If you repeat a distance, you get an
//   empty list in that position in the output.  If you give an
//   empty cutdist array, you get the input path as output
//   (without the final vertex doubled in the case of a closed path).
// Arguments:
//   path = path of any dimension or a 1-region
//   cutdist = Distance or list of distances where path is cut
//   closed = If true, treat the path as a closed polygon.  Default: false
// Example(2D,NoAxes):
//   path = circle(d=100);
//   segs = path_cut(path, [50, 200], closed=true);
//   rainbow(segs) stroke($item, endcaps="butt", width=3);
function path_cut(path,cutdist,closed) =
  is_num(cutdist) ? path_cut(path,[cutdist],closed) :
  is_1region(path) ? path_cut(path[0], cutdist, default(closed,true)):
  let(closed=default(closed,false))
  assert(is_bool(closed))
  assert(is_vector(cutdist))
  assert(last(cutdist)<path_length(path,closed=closed)-EPSILON,"\nCut distances must be smaller than the path length.")
  assert(cutdist[0]>EPSILON, "\nCut distances must be strictly positive.")
  let(
      cutlist = path_cut_points(path,cutdist,closed=closed)
  )
  _path_cut_getpaths(path, cutlist, closed);


function _path_cut_getpaths(path, cutlist, closed) =
  let(
      cuts = len(cutlist)
  )
  [
      [ each list_head(path,cutlist[0][1]-1),
        if (!approx(cutlist[0][0], path[cutlist[0][1]-1])) cutlist[0][0]
      ],
      for(i=[0:1:cuts-2])
          cutlist[i][0]==cutlist[i+1][0] && cutlist[i][1]==cutlist[i+1][1] ? []
          :
          [ if (!approx(cutlist[i][0], select(path,cutlist[i][1]))) cutlist[i][0],
            each slice(path, cutlist[i][1], cutlist[i+1][1]-1),
            if (!approx(cutlist[i+1][0], select(path,cutlist[i+1][1]-1))) cutlist[i+1][0],
          ],
      [
        if (!approx(cutlist[cuts-1][0], select(path,cutlist[cuts-1][1]))) cutlist[cuts-1][0],
        each select(path,cutlist[cuts-1][1],closed ? 0 : -1)
      ]
  ];



// Function: path_cut_points()
// Synopsis: Returns a list of cut {{points}} at a list of distances from the first point in a {{path}}.
// Topics: Paths, Path Subdivision
// See Also: path_cut(), split_path_at_self_crossings()
// Usage:
//   cuts = path_cut_points(path, cutdist, [closed=], [direction=]);
//
// Description:
//   Cuts a {{path}} at a list of distances from the first {{point}} in the path.  Returns a list of the cut
//   points and indices of the next point in the path after that point.  So for example, a return
//   value entry of [[2,3], 5] means that the cut point was [2,3] and the next point on the path after
//   this point is path[5].  If the path is too short then path_cut_points returns undef.  If you set
//   `direction` to true then `path_cut_points` also returns the tangent vector to the path and a normal
//   vector to the path.  It tries to find a normal vector that is coplanar to the path near the cut
//   point.  If this fails, it returns a normal vector parallel to the xy plane.  The output with
//   direction vectors are in the form `[point, next_index, tangent, normal]`.
//   .
//   If you give the last point of the path as a cut point, then the returned index is
//   one larger than the last index (so it would not be a valid index).  If you use the closed
//   option then the returned index is equal to the path length for cuts along the closing
//   path segment, and if you give a point equal to the path length you get an
//   index of len(path)+1 for the index.  
//
// Arguments:
//   path = path to cut
//   cutdist = distances where the path should be cut (a list) or a scalar single distance
//   ---
//   closed = set to true if the curve is closed.  Default: false
//   direction = set to true to return direction vectors.  Default: false
//
// Example(NORENDER):
//   square=[[0,0],[1,0],[1,1],[0,1]];
//   path_cut_points(square, [.5,1.5,2.5]);   // Returns [[[0.5, 0], 1], [[1, 0.5], 2], [[0.5, 1], 3]]
//   path_cut_points(square, [0,1,2,3]);      // Returns [[[0, 0], 1], [[1, 0], 2], [[1, 1], 3], [[0, 1], 4]]
//   path_cut_points(square, [0,0.8,1.6,2.4,3.2], closed=true);  // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], [[0, 0.8], 4]]
//   path_cut_points(square, [0,0.8,1.6,2.4,3.2]);               // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], undef]
function path_cut_points(path, cutdist, closed=false, direction=false) =
    let(long_enough = len(path) >= (closed ? 3 : 2))
    assert(long_enough,len(path)<2 ? "\nTwo points needed to define a path." : "\nClosed path must include three points.")
    is_num(cutdist) ? path_cut_points(path, [cutdist],closed, direction)[0] :
    assert(is_vector(cutdist))
    assert(is_increasing(cutdist), "\nCut distances must be an increasing list.")
    let(cuts = path_cut_points_recurse(path,cutdist,closed))
    !direction
       ? cuts
       : let(
             dir = _path_cuts_dir(path, cuts, closed),
             normals = _path_cuts_normals(path, cuts, dir, closed)
         )
         hstack(cuts, list_to_matrix(dir,1), list_to_matrix(normals,1));

// Main recursive path cut function
function path_cut_points_recurse(path, dists, closed=false, pind=0, dtotal=0, dind=0, result=[]) =
    dind == len(dists) ? result :
    let(
        lastpt = len(result)==0? [] : last(result)[0],       // location of last cut point
        dpartial = len(result)==0? 0 : norm(lastpt-select(path,pind)),  // remaining length in segment
        nextpoint = dists[dind] < dpartial+dtotal  // Do we have enough length left on the current segment?
           ? [lerp(lastpt,select(path,pind),(dists[dind]-dtotal)/dpartial),pind] 
           : _path_cut_single(path, dists[dind]-dtotal-dpartial, closed, pind)
    ) 
    path_cut_points_recurse(path, dists, closed, nextpoint[1], dists[dind],dind+1, concat(result, [nextpoint]));


// Search for a single cut point in the path
function _path_cut_single(path, dist, closed=false, ind=0, eps=1e-7) =
    // If we get to the end of the path (ind is last point or wraparound for closed case) then
    // check if we are within epsilon of the final path point.  If not we're out of path, so we fail
    ind==len(path)-(closed?0:1) ?
       assert(dist<eps,"\nPath is too short for specified cut distance.")
       [select(path,ind),ind+1]
    :let(d = norm(path[ind]-select(path,ind+1))) d > dist ?
        [lerp(path[ind],select(path,ind+1),dist/d), ind+1] :
        _path_cut_single(path, dist-d,closed, ind+1, eps);

// Find normal directions to the path, coplanar to local part of the path
// Or return a vector parallel to the x-y plane if the above fails
function _path_cuts_normals(path, cuts, dirs, closed=false) =
    [for(i=[0:len(cuts)-1])
        len(path[0])==2? [-dirs[i].y, dirs[i].x]
          : 
            let(
                plane = len(path)<3 ? undef :
                let(start = max(min(cuts[i][1],len(path)-1),2)) _path_plane(path, start, start-2)
            )
            plane==undef?
                ( dirs[i].x==0 && dirs[i].y==0 ? [1,0,0]  // If it's z direction return x vector
                                               : unit([-dirs[i].y, dirs[i].x,0])) // otherwise perpendicular to projection
                : unit(cross(dirs[i],cross(plane[0],plane[1])))
    ];

// Scan from the specified point (ind) to find a noncoplanar triple to use
// to define the plane of the path.
function _path_plane(path, ind, i,closed) =
    i<(closed?-1:0) ? undef :
    !is_collinear(path[ind],path[ind-1], select(path,i))?
        [select(path,i)-path[ind-1],path[ind]-path[ind-1]] :
        _path_plane(path, ind, i-1);

// Find the direction of the path at the cut points
function _path_cuts_dir(path, cuts, closed=false, eps=1e-2) =
    [for(ind=[0:len(cuts)-1])
        let(
            zeros = path[0]*0,
            nextind = cuts[ind][1],
            nextpath = unit(select(path, nextind+1)-select(path, nextind),zeros),
            thispath = unit(select(path, nextind) - select(path,nextind-1),zeros),
            lastpath = unit(select(path,nextind-1) - select(path, nextind-2),zeros),
            nextdir =
                nextind==len(path) && !closed? lastpath :
                (nextind<=len(path)-2 || closed) && approx(cuts[ind][0], path[nextind],eps)
                   ? unit(nextpath+thispath)
              : (nextind>1 || closed) && approx(cuts[ind][0],select(path,nextind-1),eps)
                   ? unit(thispath+lastpath)
              :  thispath
        ) nextdir
    ];


// internal function
// converts pathcut output form to a [segment, u]
// form list that works withi path_select
function _cut_to_seg_u_form(pathcut, path, closed) =
  let(lastind = len(path) - (closed?0:1))
  [for(entry=pathcut)
    entry[1] > lastind ? [lastind,0] :
    let(
        a = path[entry[1]-1],
        b = path[entry[1]],
        c = entry[0],
        i = max_index(v_abs(b-a)),
        factor = (c[i]-a[i])/(b[i]-a[i])
    )
    [entry[1]-1,factor]
  ];



// Function: split_path_at_self_crossings()
// Synopsis: Split a 2D {{path}} wherever it crosses itself.
// SynTags: PathList
// Topics: Paths, Path Subdivision
// See Also: path_cut(), path_cut_points()
// Usage:
//   paths = split_path_at_self_crossings(path, [closed], [eps]);
// Description:
//   Splits a 2D {{path}} into sub-paths wherever the original path crosses itself.
//   Splits may occur mid-segment, so new vertices are created at the intersection points.
//   Returns a list of the resulting subpaths.  
// Arguments:
//   path = A 2D path or a 1-region.
//   closed = If true, treat path as a closed polygon.  Default: true
//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
// Example(2D,NoAxes):
//   path = [ [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100] ];
//   paths = split_path_at_self_crossings(path);
//   rainbow(paths) stroke($item, closed=false, width=3);
function split_path_at_self_crossings(path, closed=true, eps=EPSILON) =
    let(path = force_path(path))
    assert(is_path(path,2), "\nMust give a 2D path.")
    assert(is_bool(closed))
    let(
        path = list_unwrap(path, eps=eps),
        isects = deduplicate(
            eps=eps,
            concat(
                [[0, 0]],
                sort([
                    for (
                        a = _path_self_intersections(path, closed=closed, eps=eps),
                        ss = [ [a[1],a[2]], [a[3],a[4]] ]
                    ) if (ss[0] != undef) ss
                ]),
                [[len(path)-(closed?1:2), 1]]
            )
        )
    ) [
        for (p = pair(isects))
            let(
                s1 = p[0][0],
                u1 = p[0][1],
                s2 = p[1][0],
                u2 = p[1][1],
                section = _path_select(path, s1, u1, s2, u2, closed=closed),
                outpath = deduplicate(eps=eps, section)
            )
            if (len(outpath)>1) outpath
    ];


function _tag_self_crossing_subpaths(path, nonzero, closed=true, eps=EPSILON) =
    let(
        subpaths = split_path_at_self_crossings(
            path, closed=true, eps=eps
        )
    ) [
        for (subpath = subpaths) let(
            seg = select(subpath,0,1),
            mp = mean(seg),
            n = line_normal(seg) / 2048,
            p1 = mp + n,
            p2 = mp - n,
            p1in = point_in_polygon(p1, path, nonzero=nonzero) >= 0,
            p2in = point_in_polygon(p2, path, nonzero=nonzero) >= 0,
            tag = (p1in && p2in)? "I" : "O"
        ) [tag, subpath]
    ];


// Function: polygon_parts()
// Synopsis: Parses a self-intersecting polygon into a list of non-intersecting {{polygons}}.
// SynTags: PathList
// Topics: Paths, Polygons
// See Also: split_path_at_self_crossings(), path_cut(), path_cut_points()
// Usage:
//   splitpolys = polygon_parts(poly, [nonzero], [eps]);
// Description:
//   Given a possibly self-intersecting 2D {{polygon}}, constructs a representation of the original polygon as a list of
//   non-intersecting simple polygons.  If nonzero is set to true then it uses the nonzero method for defining polygon membership.
//   For simple cases, such as the pentagram, this produces the outer perimeter of a self-intersecting polygon.  
// Arguments:
//   poly = a 2D polygon or 1-region
//   nonzero = If true use the nonzero method for checking if a point is in a polygon.  Otherwise use the even-odd method.  Default: false
//   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
// Example(2D,NoAxes):  This cross-crossing polygon breaks up into its 3 components (regardless of the value of nonzero).
//   poly = [
//       [-100,100], [0,-50], [100,100],
//       [100,-100], [0,50], [-100,-100]
//   ];
//   splitpolys = polygon_parts(poly);
//   rainbow(splitpolys) stroke($item, closed=true, width=3);
// Example(2D,NoAxes): With nonzero=false you get even-odd mode that matches OpenSCAD, so the pentagram breaks apart into its five points.
//   pentagram = turtle(["move",100,"left",144], repeat=4);
//   left(100)polygon(pentagram);
//   rainbow(polygon_parts(pentagram,nonzero=false))
//     stroke($item,closed=true,width=2.5);
// Example(2D,NoAxes): With nonzero=true you get only the outer perimeter.  You can use this to create the polygon using the nonzero method, which is not supported by OpenSCAD.
//   pentagram = turtle(["move",100,"left",144], repeat=4);
//   outside = polygon_parts(pentagram,nonzero=true);
//   left(100)region(outside);
//   rainbow(outside)
//     stroke($item,closed=true,width=2.5);
// Example(2D,NoAxes): 
//   N=12;
//   ang=360/N;
//   sr=10;
//   poly = turtle(["angle", 90+ang/2,
//                  "move", sr, "left",
//                  "move", 2*sr*sin(ang/2), "left",
//                  "repeat", 4,
//                     ["move", 2*sr, "left",
//                      "move", 2*sr*sin(ang/2), "left"],
//                  "move", sr]);
//   stroke(poly, width=.3);
//   right(20)rainbow(polygon_parts(poly)) polygon($item);
// Example(2D,NoAxes): overlapping poly segments disappear
//   poly = [[0,0], [10,0], [10,10], [0,10],[0,20], [20,10],[10,10], [0,10],[0,0]];
//   stroke(poly,width=0.3);
//   right(22)stroke(polygon_parts(poly)[0], width=0.3, closed=true);
// Example(2D,NoAxes): Poly segments disappear outside as well
//   poly = turtle(["repeat", 3, ["move", 17, "left", "move", 10, "left", "move", 7, "left", "move", 10, "left"]]);
//   back(2)stroke(poly,width=.5);
//   fwd(12)rainbow(polygon_parts(poly)) stroke($item, closed=true, width=0.5);
// Example(2D,NoAxes):  This shape has six components
//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 17, "left"]]);
//   polygon(poly);
//   right(22)rainbow(polygon_parts(poly)) polygon($item);
// Example(2D,NoAxes): When the loops of the shape overlap then nonzero gives a different result than the even-odd method.
//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 10, "left"]]);
//   polygon(poly);
//   right(27)rainbow(polygon_parts(poly)) polygon($item);
//   move([16,-14])rainbow(polygon_parts(poly,nonzero=true)) polygon($item);
function polygon_parts(poly, nonzero=false, eps=EPSILON) =
    let(poly = force_path(poly))
    assert(is_path(poly,2), "\nMust give 2D polygon.")
    assert(is_bool(nonzero))    
    let(
        poly = list_unwrap(poly, eps=eps),
        tagged = _tag_self_crossing_subpaths(poly, nonzero=nonzero, closed=true, eps=eps),
        kept = [for (sub = tagged) if(sub[0] == "O") sub[1]],
        outregion = _assemble_path_fragments(kept, eps=eps)
    ) outregion;


function _extreme_angle_fragment(seg, fragments, rightmost=true, eps=EPSILON) =
    !fragments? [undef, []] :
    let(
        delta = seg[1] - seg[0],
        segang = atan2(delta.y,delta.x),
        frags = [
            for (i = idx(fragments)) let(
                fragment = fragments[i],
                fwdmatch = approx(seg[1], fragment[0], eps=eps),
                bakmatch =  approx(seg[1], last(fragment), eps=eps)
            ) [
                fwdmatch,
                bakmatch,
                bakmatch? reverse(fragment) : fragment
            ]
        ],
        angs = [
            for (frag = frags)
                (frag[0] || frag[1])? let(
                    delta2 = frag[2][1] - frag[2][0],
                    segang2 = atan2(delta2.y, delta2.x)
                ) modang(segang2 - segang) : (
                    rightmost? 999 : -999
                )
        ],
        fi = rightmost? min_index(angs) : max_index(angs)
    ) abs(angs[fi]) > 360? [undef, fragments] : let(
        remainder = [for (i=idx(fragments)) if (i!=fi) fragments[i]],
        frag = frags[fi],
        foundfrag = frag[2]
    ) [foundfrag, remainder];


/// Internal Function: _assemble_a_path_from_fragments()
/// Usage:
///   _assemble_a_path_from_fragments(subpaths);
/// Description:
///   Given a list of paths, assembles them together into one complete closed polygon path, and
///   remainder fragments.  Returns [PATH, FRAGMENTS] where FRAGMENTS is the list of remaining
///   unused path fragments.
/// Arguments:
///   fragments = List of paths to be assembled into complete polygons.
///   rightmost = If true, assemble paths using rightmost turns. Leftmost if false.
///   startfrag = The fragment to start with.  Default: 0
///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
function _assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0, eps=EPSILON) =
    len(fragments)==0? [[],[]] :
    len(fragments)==1? [fragments[0],[]] :
    let(
        path = fragments[startfrag],
        newfrags = [for (i=idx(fragments)) if (i!=startfrag) fragments[i]]
    ) are_ends_equal(path, eps=eps)? (
        // starting fragment is already closed
        [path, newfrags]
    ) : let(
        // Find rightmost/leftmost continuation fragment
        seg = select(path,-2,-1),
        extrema = _extreme_angle_fragment(seg=seg, fragments=newfrags, rightmost=rightmost, eps=eps),
        foundfrag = extrema[0],
        remainder = extrema[1]
    ) is_undef(foundfrag)? (
        // No remaining fragments connect!  INCOMPLETE PATH!
        // Treat it as complete.
        [path, remainder]
    ) : are_ends_equal(foundfrag, eps=eps)? (
        // Found fragment is already closed
        [foundfrag, concat([path], remainder)]
    ) : let(
        fragend = last(foundfrag),
        hits = [for (i = idx(path,e=-2)) if(approx(path[i],fragend,eps=eps)) i]
    ) hits? (
        let(
            // Found fragment intersects with initial path
            hitidx = last(hits),
            newpath = list_head(path,hitidx),
            newfrags = concat(len(newpath)>1? [newpath] : [], remainder),
            outpath = concat(slice(path,hitidx,-2), foundfrag)
        )
        [outpath, newfrags]
    ) : let(
        // Path still incomplete.  Continue building it.
        newpath = concat(path, list_tail(foundfrag)),
        newfrags = concat([newpath], remainder)
    )
    _assemble_a_path_from_fragments(
        fragments=newfrags,
        rightmost=rightmost,
        eps=eps
    );


/// Internal Function: _assemble_path_fragments()
/// Usage:
///   _assemble_path_fragments(subpaths);
/// Description:
///   Given a list of paths, assembles them together into complete closed polygon paths if it can.
///   Polygons with area < eps are discarded and not returned.  
/// Arguments:
///   fragments = List of paths to be assembled into complete polygons.
///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
    len(fragments)==0? _finished :
    let(
        minxidx = min_index([
            for (frag=fragments) min(column(frag,0))
        ]),
        result_l = _assemble_a_path_from_fragments(
            fragments=fragments,
            startfrag=minxidx,
            rightmost=false,
            eps=eps
        ),
        result_r = _assemble_a_path_from_fragments(
            fragments=fragments,
            startfrag=minxidx,
            rightmost=true,
            eps=eps
        ),
        l_area = abs(polygon_area(result_l[0])),
        r_area = abs(polygon_area(result_r[0])),
        result = l_area < r_area? result_l : result_r,
        newpath = list_unwrap(result[0]),
        remainder = result[1],
        finished = min(l_area,r_area)<eps ? _finished : concat(_finished, [newpath])
    ) _assemble_path_fragments(
        fragments=remainder,
        eps=eps,
        _finished=finished
    );


/// Different but similar path assembly function that is much faster than
/// _assemble_path_fragments and can work in 3d, but cannot handle loops.
///
/// Takes a list of paths that are in the correct direction and assembles
/// them into a list of paths.  Returns a list of assembled paths.
/// If `closed=false` then any paths that are closed have duplicate
/// endpoints, and open paths do not have duplicate endpoints.
/// If `closed=true` then all paths are assumed closed and none of the returned
/// paths have duplicate endpoints.
///
/// It is assumed that the paths do not intersect each other.
/// Paths can be in any dimension

function _assemble_partial_paths(paths, closed=false, eps=1e-7) =
    let(
        pathlist = _assemble_partial_paths_recur(paths, eps)
        //// this eliminates crossing paths that cross only at vertices in the input paths lists
        // splitpaths =
        //     [for(path=pathlist) each
        //        let(
        //            searchlist = vector_search(path,eps,path),
        //            duplist = [for(i=idx(searchlist)) if (len(searchlist[i])>1) i]
        //        )
        //        duplist==[] ? [path]
        //       :               
        //        let(
        //            fragments = [for(i=idx(duplist)) select(path, duplist[i], select(duplist,i+1))]
        //        )
        //        len(fragments)==1 ? fragments
        //                          : _assemble_path_fragments(fragments)
        //     ]
    )
    closed ? [for(path=pathlist) list_unwrap(path)] : pathlist;


function _assemble_partial_paths_recur(edges, eps, paths=[], i=0) =
    i==len(edges) ? paths :
    norm(edges[i][0]-last(edges[i]))<eps ? _assemble_partial_paths_recur(edges, eps, paths,i+1) :
    let(    // Find paths that connects on left side and right side of the edges (if one exists)
        
        left = [for(j=idx(paths)) if (approx(last(paths[j]),edges[i][0],eps)) j],
        right = [for(j=idx(paths)) if (approx(last(edges[i]),paths[j][0],eps)) j]
    )
    let(
        keep_path = list_remove(paths,[if (len(left)>0) left[0],if (len(right)>0) right[0]]),
        update_path =  left==[] && right==[] ? edges[i]
                    : left==[] ? concat(list_head(edges[i]),paths[right[0]])
                    : right==[] ?  concat(paths[left[0]],slice(edges[i],1,-1))
                    : left[0] != right[0] ? concat(paths[left[0]],slice(edges[i],1,-2), paths[right[0]])
                    : concat(paths[left[0]], slice(edges[i],1,-1)) // last arg -2 removes duplicate endpoints but this is handled in passthrough function
    )
    _assemble_partial_paths_recur(edges, eps, concat(keep_path, [update_path]), i+1);



// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap