1//////////////////////////////////////////////////////////////////////
   2// LibFile: paths.scad
   3//   A `path` is a list of points of the same dimensions, usually 2D or 3D, that can
   4//   be connected together to form a sequence of line segments or a polygon.
   5//   A `region` is a list of paths that represent polygons, and the functions
   6//   in this file work on paths and also 1-regions, which are regions
   7//   that include exactly one path.  When you pass a 1-region to a function, the default
   8//   value for `closed` is always `true` because regions represent polygons.  
   9//   Capabilities include computing length of paths, computing
  10//   path tangents and normals, resampling of paths, and cutting paths up into smaller paths.  
  11// Includes:
  12//   include <BOSL2/std.scad>
  13// FileGroup: Advanced Modeling
  14// FileSummary: Operations on paths: length, resampling, tangents, splitting into subpaths
  15// FileFootnotes: STD=Included in std.scad
  16//////////////////////////////////////////////////////////////////////
  17
  18// Section: Utility Functions
  19
  20// Function: is_path()
  21// Synopsis: Returns True if 'list' is a path.
  22// Topics: Paths
  23// See Also: is_region(), is_vnf()
  24// Usage:
  25//   is_path(list, [dim], [fast])
  26// Description:
  27//   Returns true if `list` is a path.  A path is a list of two or more numeric vectors (AKA points).
  28//   All vectors must of the same size, and may only contain numbers that are not inf or nan.
  29//   By default the vectors in a path must be 2d or 3d.  Set the `dim` parameter to specify a list
  30//   of allowed dimensions, or set it to `undef` to allow any dimension.  (Note that this function
  31//   returns `false` on 1-regions.)  
  32// Example:
  33//   bool1 = is_path([[3,4],[5,6]]);    // Returns true
  34//   bool2 = is_path([[3,4]]);          // Returns false
  35//   bool3 = is_path([[3,4],[4,5]],2);  // Returns true
  36//   bool4 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  37//   bool5 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  38//   bool6 = is_path([[3,4,5],undef,[4,5,6]]);  // Returns false
  39//   bool7 = is_path([[3,5],[undef,undef],[4,5]]);  // Returns false
  40//   bool8 = is_path([[3,4],[5,6],[5,3]]);     // Returns true
  41//   bool9 = is_path([3,4,5,6,7,8]);           // Returns false
  42//   bool10 = is_path([[3,4],[5,6]], dim=[2,3]);// Returns true
  43//   bool11 = is_path([[3,4],[5,6]], dim=[1,3]);// Returns false
  44//   bool12 = is_path([[3,4],"hello"], fast=true); // Returns true
  45//   bool13 = is_path([[3,4],[3,4,5]]);            // Returns false
  46//   bool14 = is_path([[1,2,3,4],[2,3,4,5]]);      // Returns false
  47//   bool15 = is_path([[1,2,3,4],[2,3,4,5]],undef);// Returns true
  48// Arguments:
  49//   list = list to check
  50//   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
  51//   fast = set to true for fast check that only looks at first entry.  Default: false
  52function is_path(list, dim=[2,3], fast=false) =
  53    fast
  54    ?   is_list(list) && is_vector(list[0]) 
  55    :   is_matrix(list) 
  56        && len(list)>1 
  57        && len(list[0])>0
  58        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));
  59
  60// Function: is_1region()
  61// Synopsis: Returns true if path is a region with one component.
  62// Topics: Paths, Regions
  63// See Also: force_path()
  64// Usage:
  65//   bool = is_1region(path, [name])
  66// Description:
  67//   If `path` is a region with one component (a 1-region) then return true.  If path is a region with more components
  68//   then display an error message about the parameter `name` requiring a path or a single component region.  If the input
  69//   is not a region then return false.  This function helps path functions accept 1-regions.
  70// Arguments:
  71//   path = input to process
  72//   name = name of parameter to use in error message.  Default: "path"
  73function is_1region(path, name="path") = 
  74     !is_region(path)? false
  75    :assert(len(path)==1,str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  76     true;
  77
  78
  79// Function: force_path()
  80// Synopsis: Checks that path is a region with one component.
  81// SynTags: Path
  82// Topics: Paths, Regions
  83// See Also: is_1region()
  84// Usage:
  85//   outpath = force_path(path, [name])
  86// Description:
  87//   If `path` is a region with one component (a 1-region) then returns that component as a path.
  88//   If path is a region with more components then displays an error message about the parameter
  89//   `name` requiring a path or a single component region.  If the input is not a region then
  90//   returns the input without any checks.  This function helps path functions accept 1-regions.
  91// Arguments:
  92//   path = input to process
  93//   name = name of parameter to use in error message.  Default: "path"
  94function force_path(path, name="path") =
  95   is_region(path) ?
  96       assert(len(path)==1, str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  97       path[0]
  98   : path;
  99
 100
 101/// Internal Function: _path_select()
 102/// Usage:
 103///   _path_select(path,s1,u1,s2,u2,[closed]):
 104/// Description:
 105///   Returns a portion of a path, from between the `u1` part of segment `s1`, to the `u2` part of
 106///   segment `s2`.  Both `u1` and `u2` are values between 0.0 and 1.0, inclusive, where 0 is the start
 107///   of the segment, and 1 is the end.  Both `s1` and `s2` are integers, where 0 is the first segment.
 108/// Arguments:
 109///   path = The path to get a section of.
 110///   s1 = The number of the starting segment.
 111///   u1 = The proportion along the starting segment, between 0.0 and 1.0, inclusive.
 112///   s2 = The number of the ending segment.
 113///   u2 = The proportion along the ending segment, between 0.0 and 1.0, inclusive.
 114///   closed = If true, treat path as a closed polygon.
 115function _path_select(path, s1, u1, s2, u2, closed=false) =
 116    let(
 117        lp = len(path),
 118        l = lp-(closed?0:1),
 119        u1 = s1<0? 0 : s1>l? 1 : u1,
 120        u2 = s2<0? 0 : s2>l? 1 : u2,
 121        s1 = constrain(s1,0,l),
 122        s2 = constrain(s2,0,l),
 123        pathout = concat(
 124            (s1<l && u1<1)? [lerp(path[s1],path[(s1+1)%lp],u1)] : [],
 125            [for (i=[s1+1:1:s2]) path[i]],
 126            (s2<l && u2>0)? [lerp(path[s2],path[(s2+1)%lp],u2)] : []
 127        )
 128    ) pathout;
 129
 130
 131
 132// Function: path_merge_collinear()
 133// Synopsis: Removes unnecessary points from a path.
 134// SynTags: Path
 135// Topics: Paths, Regions
 136// Description:
 137//   Takes a path and removes unnecessary sequential collinear points.
 138// Usage:
 139//   path_merge_collinear(path, [eps])
 140// Arguments:
 141//   path = A path of any dimension or a 1-region
 142//   closed = treat as closed polygon.  Default: false
 143//   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
 144function path_merge_collinear(path, closed, eps=EPSILON) =
 145    is_1region(path) ? path_merge_collinear(path[0], default(closed,true), eps) :
 146    let(closed=default(closed,false))
 147    assert(is_bool(closed))
 148    assert( is_path(path), "Invalid path in path_merge_collinear." )
 149    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "Invalid tolerance." )    
 150    len(path)<=2 ? path :
 151    let(
 152        indices = [
 153            0,
 154            for (i=[1:1:len(path)-(closed?1:2)]) 
 155                if (!is_collinear(path[i-1], path[i], select(path,i+1), eps=eps)) i, 
 156            if (!closed) len(path)-1 
 157        ]
 158    ) [for (i=indices) path[i]];
 159
 160
 161
 162// Section: Path length calculation
 163
 164
 165// Function: path_length()
 166// Synopsis: Returns the path length.
 167// Topics: Paths
 168// See Also: path_segment_lengths(), path_length_fractions()
 169// Usage:
 170//   path_length(path,[closed])
 171// Description:
 172//   Returns the length of the path.
 173// Arguments:
 174//   path = Path of any dimension or 1-region. 
 175//   closed = true if the path is closed.  Default: false
 176// Example:
 177//   path = [[0,0], [5,35], [60,-25], [80,0]];
 178//   echo(path_length(path));
 179function path_length(path,closed) =
 180    is_1region(path) ? path_length(path[0], default(closed,true)) :
 181    assert(is_path(path), "Invalid path in path_length")
 182    let(closed=default(closed,false))
 183    assert(is_bool(closed))
 184    len(path)<2? 0 :
 185    sum([for (i = [0:1:len(path)-2]) norm(path[i+1]-path[i])])+(closed?norm(path[len(path)-1]-path[0]):0);
 186
 187
 188// Function: path_segment_lengths()
 189// Synopsis: Returns a list of the lengths of segments in a path.
 190// Topics: Paths
 191// See Also: path_length(), path_length_fractions()
 192// Usage:
 193//   path_segment_lengths(path,[closed])
 194// Description:
 195//   Returns list of the length of each segment in a path
 196// Arguments:
 197//   path = path in any dimension or 1-region
 198//   closed = true if the path is closed.  Default: false
 199function path_segment_lengths(path, closed) =
 200    is_1region(path) ? path_segment_lengths(path[0], default(closed,true)) :
 201    let(closed=default(closed,false))
 202    assert(is_path(path),"Invalid path in path_segment_lengths.")
 203    assert(is_bool(closed))
 204    [
 205        for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
 206        if (closed) norm(path[0]-last(path))
 207    ]; 
 208
 209
 210// Function: path_length_fractions()
 211// Synopsis: Returns the fractional distance of each point along the length of a path.
 212// Topics: Paths
 213// See Also: path_length(), path_segment_lengths()
 214// Usage:
 215//   fracs = path_length_fractions(path, [closed]);
 216// Description:
 217//    Returns the distance fraction of each point in the path along the path, so the first
 218//    point is zero and the final point is 1.  If the path is closed the length of the output
 219//    will have one extra point because of the final connecting segment that connects the last
 220//    point of the path to the first point.
 221// Arguments:
 222//    path = path in any dimension or a 1-region
 223//    closed = set to true if path is closed.  Default: false
 224function path_length_fractions(path, closed) =
 225    is_1region(path) ? path_length_fractions(path[0], default(closed,true)):
 226    let(closed=default(closed, false))
 227    assert(is_path(path))
 228    assert(is_bool(closed))
 229    let(
 230        lengths = [
 231            0,
 232            each path_segment_lengths(path,closed)
 233        ],
 234        partial_len = cumsum(lengths),
 235        total_len = last(partial_len)
 236    )
 237    partial_len / total_len;
 238
 239
 240
 241/// Internal Function: _path_self_intersections()
 242/// Usage:
 243///   isects = _path_self_intersections(path, [closed], [eps]);
 244/// Description:
 245///   Locates all self intersection points of the given path.  Returns a list of intersections, where
 246///   each intersection is a list like [POINT, SEGNUM1, PROPORTION1, SEGNUM2, PROPORTION2] where
 247///   POINT is the coordinates of the intersection point, SEGNUMs are the integer indices of the
 248///   intersecting segments along the path, and the PROPORTIONS are the 0.0 to 1.0 proportions
 249///   of how far along those segments they intersect at.  A proportion of 0.0 indicates the start
 250///   of the segment, and a proportion of 1.0 indicates the end of the segment.
 251///   .
 252///   Note that this function does not return self-intersecting segments, only the points
 253///   where non-parallel segments intersect.  
 254/// Arguments:
 255///   path = The path to find self intersections of.
 256///   closed = If true, treat path like a closed polygon.  Default: true
 257///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
 258/// Example(2D):
 259///   path = [
 260///       [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100]
 261///   ];
 262///   isects = _path_self_intersections(path, closed=true);
 263///   // isects == [[[-33.3333, 0], 0, 0.666667, 4, 0.333333], [[33.3333, 0], 1, 0.333333, 3, 0.666667]]
 264///   stroke(path, closed=true, width=1);
 265///   for (isect=isects) translate(isect[0]) color("blue") sphere(d=10);
 266function _path_self_intersections(path, closed=true, eps=EPSILON) =
 267    let(
 268        path = closed ? list_wrap(path,eps=eps) : path,
 269        plen = len(path)
 270    )
 271    [ for (i = [0:1:plen-3]) let(
 272          a1 = path[i],
 273          a2 = path[i+1], 
 274          seg_normal = unit([-(a2-a1).y, (a2-a1).x],[0,0]),
 275          vals = path*seg_normal,
 276          ref  = a1*seg_normal,
 277            // The value of vals[j]-ref is positive if vertex j is one one side of the
 278            // line [a1,a2] and negative on the other side. Only a segment with opposite
 279            // signs at its two vertices can have an intersection with segment
 280            // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than
 281            // eps and the sign of vals[j]-ref otherwise.  
 282          signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) 
 283                        abs(vals[j]-ref) <  eps ? 0 : sign(vals[j]-ref) ] 
 284        )
 285        if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2]
 286        for(j=[i+2:1:plen-(i==0 && closed? 3: 2)])
 287            if( signals[j-i-2]*signals[j-i-1]<=0 ) let( // segm [b1,b2] intersects line [a1,a2]
 288                b1 = path[j],
 289                b2 = path[j+1],
 290                isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
 291            )
 292            if (isect 
 293                && isect[1]>=-eps
 294                && isect[1]<= 1+eps
 295                && isect[2]>= -eps 
 296                && isect[2]<= 1+eps)
 297                [isect[0], i, isect[1], j, isect[2]]
 298    ];
 299
 300// Section: Resampling - changing the number of points in a path
 301
 302
 303// Input `data` is a list that sums to an integer. 
 304// Returns rounded version of input data so that every 
 305// entry is rounded to an integer and the sum is the same as
 306// that of the input.  Works by rounding an entry in the list
 307// and passing the rounding error forward to the next entry.
 308// This will generally distribute the error in a uniform manner. 
 309function _sum_preserving_round(data, index=0) =
 310    index == len(data)-1 ? list_set(data, len(data)-1, round(data[len(data)-1])) :
 311    let(
 312        newval = round(data[index]),
 313        error = newval - data[index]
 314    ) _sum_preserving_round(
 315        list_set(data, [index,index+1], [newval, data[index+1]-error]),
 316        index+1
 317    );
 318
 319
 320// Function: subdivide_path()
 321// Synopsis: Subdivides a path to produce a more finely sampled path.
 322// SynTags: Path
 323// Topics: Paths, Path Subdivision
 324// See Also: subdivide_and_slice(), resample_path(), jittered_poly()
 325// Usage:
 326//   newpath = subdivide_path(path, n|refine=|maxlen=, [method=], [closed=], [exact=]);
 327// Description:
 328//   Takes a path as input (closed or open) and subdivides the path to produce a more
 329//   finely sampled path.  You control the subdivision process by using the `maxlen` arg
 330//   to specify a maximum segment length, or by specifying `n` or `refine`, which request
 331//   a certain point count in the output.
 332//   .
 333//   You can specify the point count using the `n` option, where
 334//   you give the number of points you want in the output, or you can use
 335//   the `refine` option, where you specify a resampling factor.  If `refine=3` then
 336//   the number of points would increase by a factor of three, so a four point square would
 337//   have 12 points after subdivision.  With point-count subdivision, the new points can be distributed
 338//   proportional to length (`method="length"`), which is the default, or they can be divided up evenly among all the path segments
 339//   (`method="segment"`).  If the extra points don't fit evenly on the path then the
 340//   algorithm attempts to distribute them as uniformly as possible, but the result may be uneven.
 341//   The `exact` option, which is true by default, requires that the final point count is
 342//   exactly as requested.  For example, if you subdivide a four point square and request `n=13` then one edge will have
 343//   an extra point compared to the others.  
 344//   If you set `exact=false` then the
 345//   algorithm will favor uniformity and the output path may have a different number of
 346//   points than you requested, but the sampling will be uniform.   In our example of the
 347//   square with `n=13`, you will get only 12 points output, with the same number of points on each edge.
 348//   .
 349//   The points are always distributed uniformly on each segment.  The `method="length"` option does
 350//   means that the number of points on a segment is based on its length, but the points are still
 351//   distributed uniformly on each segment, independent of the other segments.  
 352//   With the `"segment"` method you can also give `n` as a vector of counts.  This 
 353//   specifies the desired point count on each segment: with vector valued `n` the `subdivide_path`
 354//   function places `n[i]-1` points on segment `i`.  The reason for the -1 is to avoid
 355//   double counting the endpoints, which are shared by pairs of segments, so that for
 356//   a closed polygon the total number of points will be sum(n).  Note that with an open
 357//   path there is an extra point at the end, so the number of points will be sum(n)+1.
 358//   .
 359//   If you use the `maxlen` option then you specify the maximum length segment allowed in the output.
 360//   Each segment is subdivided into the largest number of segments meeting your requirement.  As above,
 361//   the sampling is uniform on each segment, independent of the other segments.  With the `maxlen` option
 362//   you cannot specify `method` or `exact`.    
 363// Arguments:
 364//   path = path in any dimension or a 1-region
 365//   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.
 366//   ---
 367//   refine = increase total number of points by this factor (Specify only one of n, refine and maxlen)
 368//   maxlen = maximum length segment in the output (Specify only one of n, refine and maxlen)
 369//   closed = set to false if the path is open.  Default: True
 370//   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
 371//   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"`
 372// Example(2D):
 373//   mypath = subdivide_path(square([2,2],center=true), 12);
 374//   move_copies(mypath)circle(r=.1,$fn=32);
 375// Example(2D):
 376//   mypath = subdivide_path(square([8,2],center=true), 12);
 377//   move_copies(mypath)circle(r=.2,$fn=32);
 378// Example(2D):
 379//   mypath = subdivide_path(square([8,2],center=true), 12, method="segment");
 380//   move_copies(mypath)circle(r=.2,$fn=32);
 381// Example(2D):
 382//   mypath = subdivide_path(square([2,2],center=true), 17, closed=false);
 383//   move_copies(mypath)circle(r=.1,$fn=32);
 384// Example(2D): Specifying different numbers of points on each segment
 385//   mypath = subdivide_path(hexagon(side=2), [2,3,4,5,6,7], method="segment");
 386//   move_copies(mypath)circle(r=.1,$fn=32);
 387// Example(2D): Requested point total is 14 but 15 points output due to extra end point
 388//   mypath = subdivide_path(pentagon(side=2), [3,4,3,4], method="segment", closed=false);
 389//   move_copies(mypath)circle(r=.1,$fn=32);
 390// Example(2D): Since 17 is not divisible by 5, a completely uniform distribution is not possible. 
 391//   mypath = subdivide_path(pentagon(side=2), 17);
 392//   move_copies(mypath)circle(r=.1,$fn=32);
 393// Example(2D): With `exact=false` a uniform distribution, but only 15 points
 394//   mypath = subdivide_path(pentagon(side=2), 17, exact=false);
 395//   move_copies(mypath)circle(r=.1,$fn=32);
 396// Example(2D): With `exact=false` you can also get extra points, here 20 instead of requested 18
 397//   mypath = subdivide_path(pentagon(side=2), 18, exact=false);
 398//   move_copies(mypath)circle(r=.1,$fn=32);
 399// Example(2D): Using refine in this example multiplies the point count by 3 by adding 2 points to each edge
 400//   mypath = subdivide_path(pentagon(side=2), refine=3);
 401//   move_copies(mypath)circle(r=.1,$fn=32);
 402// Example(2D): But note that 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.  
 403//   mypath = subdivide_path(square([8,2],center=true), refine=3);
 404//   move_copies(mypath)circle(r=.2,$fn=32);
 405// Example(2D): In this example with maxlen, every side gets a different number of new points
 406//   path = [[0,0],[0,4],[10,6],[10,0]];
 407//   spath = subdivide_path(path, maxlen=2, closed=true);
 408//   move_copies(spath) circle(r=.25,$fn=12);
 409// Example(FlatSpin,VPD=15,VPT=[0,0,1.5]): Three-dimensional paths also work
 410//   mypath = subdivide_path([[0,0,0],[2,0,1],[2,3,2]], 12);
 411//   move_copies(mypath)sphere(r=.1,$fn=32);
 412function subdivide_path(path, n, refine, maxlen, closed=true, exact, method) =
 413    let(path = force_path(path))
 414    assert(is_path(path))
 415    assert(num_defined([n,refine,maxlen]),"Must give exactly one of n, refine, and maxlen")
 416    refine==1 || n==len(path) ? path :
 417    is_def(maxlen) ?
 418        assert(is_undef(method), "Cannot give method with maxlen")
 419        assert(is_undef(exact), "Cannot give exact with maxlen")
 420        [
 421         for (p=pair(path,closed))
 422           let(steps = ceil(norm(p[1]-p[0])/maxlen))
 423           each lerpn(p[0], p[1], steps, false),
 424         if (!closed) last(path)
 425        ]               
 426    :
 427    let(
 428        exact = default(exact, true),
 429        method = default(method, "length")
 430    )
 431    assert(method=="length" || method=="segment")
 432    let(
 433        n = !is_undef(n)? n :
 434            !is_undef(refine)? len(path) * refine :
 435            undef
 436    )
 437    assert((is_num(n) && n>0) || is_vector(n),"Parameter n to subdivide_path must be postive number or vector")
 438    let(
 439        count = len(path) - (closed?0:1), 
 440        add_guess = method=="segment"?
 441                       (
 442                          is_list(n)
 443                          ? assert(len(n)==count,"Vector parameter n to subdivide_path has the wrong length")
 444                            add_scalar(n,-1)
 445                          : repeat((n-len(path)) / count, count)
 446                       )
 447                  : // method=="length"
 448                    assert(is_num(n),"Parameter n to subdivide path must be a number when method=\"length\"")
 449                    let(
 450                        path_lens = path_segment_lengths(path,closed),
 451                        add_density = (n - len(path)) / sum(path_lens)
 452                    )
 453                    path_lens * add_density,
 454        add = exact? _sum_preserving_round(add_guess)
 455                   : [for (val=add_guess) round(val)]
 456    )
 457    [
 458        for (i=[0:1:count-1]) 
 459           each lerpn(path[i],select(path,i+1), 1+add[i],endpoint=false),
 460        if (!closed) last(path)
 461    ];
 462
 463
 464
 465
 466// Function: resample_path()
 467// Synopsis: Returns an equidistant set of points along a path.
 468// SynTags: Path
 469// Topics: Paths
 470// See Also: subdivide_path()
 471// Usage:
 472//   newpath = resample_path(path, n|spacing=, [closed=]);
 473// Description:
 474//   Compute a uniform resampling of the input path.  If you specify `n` then the output path will have n
 475//   points spaced uniformly (by linear interpolation along the input path segments).  The only points of the
 476//   input path that are guaranteed to appear in the output path are the starting and ending points, and any
 477//   points that have an angular deflection of at least the number of degrees given in `keep_corners`.
 478//   If you specify `spacing` then the length you give will be rounded to the nearest spacing that gives
 479//   a uniform sampling of the path and the resulting uniformly sampled path is returned.
 480//   Note that because this function operates on a discrete input path the quality of the output depends on
 481//   the sampling of the input.  If you want very accurate output, use a lot of points for the input.
 482// Arguments:
 483//   path = path in any dimension or a 1-region
 484//   n = Number of points in output
 485//   ---
 486//   spacing = Approximate spacing desired
 487//   keep_corners = If given a scalar, path vertices with deflection angle greater than this are preserved in the output.
 488//   closed = set to true if path is closed.  Default: true
 489// Example(2D):  Subsampling lots of points from a smooth curve
 490//   path = xscale(2,circle($fn=250, r=10));
 491//   sampled = resample_path(path, 16);
 492//   stroke(path);
 493//   color("red")move_copies(sampled) circle($fn=16);
 494// Example(2D): Specified spacing is rounded to make a uniform sampling
 495//   path = xscale(2,circle($fn=250, r=10));
 496//   sampled = resample_path(path, spacing=17);
 497//   stroke(path);
 498//   color("red")move_copies(sampled) circle($fn=16);
 499// Example(2D): Notice that the corners are excluded.
 500//   path = square(20);
 501//   sampled = resample_path(path, spacing=6);
 502//   stroke(path,closed=true);
 503//   color("red")move_copies(sampled) circle($fn=16);
 504// Example(2D): Forcing preservation of corners.
 505//   path = square(20);
 506//   sampled = resample_path(path, spacing=6, keep_corners=90);
 507//   stroke(path,closed=true);
 508//   color("red")move_copies(sampled) circle($fn=16);
 509// Example(2D): Closed set to false
 510//   path = square(20);
 511//   sampled = resample_path(path, spacing=6,closed=false);
 512//   stroke(path);
 513//   color("red")move_copies(sampled) circle($fn=16);
 514
 515function resample_path(path, n, spacing, keep_corners, closed=true) =
 516    let(path = force_path(path))
 517    assert(is_path(path))
 518    assert(num_defined([n,spacing])==1,"Must define exactly one of n and spacing")
 519    assert(n==undef || (is_integer(n) && n>0))
 520    assert(spacing==undef || (is_finite(spacing) && spacing>0))
 521    assert(is_bool(closed))
 522    let(
 523        corners = is_undef(keep_corners)
 524          ? [0, len(path)-(closed?0:1)]
 525          : [
 526                0,
 527                for (i = [1:1:len(path)-(closed?1:2)])
 528                    let( ang = abs(modang(vector_angle(select(path,i-1,i+1))-180)) )
 529                    if (ang >= keep_corners) i,
 530                len(path)-(closed?0:1),
 531            ],
 532        pcnt = len(path),
 533        plen = path_length(path, closed=closed),
 534        subpaths = [ for (p = pair(corners)) [for(i = [p.x:1:p.y]) path[i%pcnt]] ],
 535        n = is_undef(n)? undef : closed? n+1 : n
 536    )
 537    assert(n==undef || n >= len(corners), "There are nore than `n=` corners whose angle is greater than `keep_corners=`.")
 538    let(
 539        lens = [for (subpath = subpaths) path_length(subpath)],
 540        part_ns = is_undef(n)
 541          ? [for (i=idx(subpaths)) max(1,round(lens[i]/spacing)-1)]
 542          : let(
 543                ccnt = len(corners),
 544                parts = [for (l=lens) (n-ccnt) * l/plen]
 545            )
 546            _sum_preserving_round(parts),
 547        out = [
 548            for (i = idx(subpaths))
 549                let(
 550                    subpath = subpaths[i],
 551                    splen = lens[i],
 552                    pn = part_ns[i] + 1,
 553                    distlist = lerpn(0, splen, pn, false),
 554                    cuts = path_cut_points(subpath, distlist, closed=false)
 555                )
 556                each column(cuts,0),
 557            if (!closed) last(path)
 558        ]
 559    ) out;
 560
 561
 562// Section: Path Geometry
 563
 564// Function: is_path_simple()
 565// Synopsis: Returns true if a path has no self intersections.
 566// Topics: Paths
 567// See Also: is_path()
 568// Usage:
 569//   bool = is_path_simple(path, [closed], [eps]);
 570// Description:
 571//   Returns true if the given 2D path is simple, meaning that it has no self-intersections.
 572//   Repeated points are not considered self-intersections: a path with such points can
 573//   still be simple.  
 574//   If closed is set to true then treat the path as a polygon.
 575// Arguments:
 576//   path = 2D path or 1-region
 577//   closed = set to true to treat path as a polygon.  Default: false
 578//   eps = Epsilon error value used for determine if points coincide.  Default: `EPSILON` (1e-9)
 579function is_path_simple(path, closed, eps=EPSILON) =
 580    is_1region(path) ? is_path_simple(path[0], default(closed,true), eps) :
 581    let(closed=default(closed,false))
 582    assert(is_path(path, 2),"Must give a 2D path")
 583    assert(is_bool(closed))
 584    // check for path reversals
 585    [for(i=[0:1:len(path)-(closed?2:3)])
 586         let(v1=path[i+1]-path[i],
 587             v2=select(path,i+2)-path[i+1],
 588             normv1 = norm(v1),
 589             normv2 = norm(v2)
 590             )
 591         if (approx(v1*v2/normv1/normv2,-1)) 1
 592    ]  == [] 
 593    &&
 594    _path_self_intersections(path,closed=closed,eps=eps) == [];
 595
 596
 597// Function: path_closest_point()
 598// Synopsis: Returns the closest place on a path to a given point.
 599// Topics: Paths
 600// See Also: point_line_distance(), line_closest_point()
 601// Usage:
 602//   index_pt = path_closest_point(path, pt);
 603// Description:
 604//   Finds the closest path segment, and point on that segment to the given point.
 605//   Returns `[SEGNUM, POINT]`
 606// Arguments:
 607//   path = Path of any dimension or a 1-region.
 608//   pt = The point to find the closest point to.
 609//   closed = If true, the path is considered closed.
 610// Example(2D):
 611//   path = circle(d=100,$fn=6);
 612//   pt = [20,10];
 613//   closest = path_closest_point(path, pt);
 614//   stroke(path, closed=true);
 615//   color("blue") translate(pt) circle(d=3, $fn=12);
 616//   color("red") translate(closest[1]) circle(d=3, $fn=12);
 617function path_closest_point(path, pt, closed=true) =
 618    let(path = force_path(path))
 619    assert(is_path(path), "Input must be a path")
 620    assert(is_vector(pt, len(path[0])), "Input pt must be a compatible vector")
 621    assert(is_bool(closed))
 622    let(
 623        pts = [for (seg=pair(path,closed)) line_closest_point(seg,pt,SEGMENT)],
 624        dists = [for (p=pts) norm(p-pt)],
 625        min_seg = min_index(dists)
 626    ) [min_seg, pts[min_seg]];
 627
 628
 629// Function: path_tangents()
 630// Synopsis: Returns tangent vectors for each point along a path.
 631// Topics: Paths
 632// See Also: path_normals()
 633// Usage:
 634//   tangs = path_tangents(path, [closed], [uniform]);
 635// Description:
 636//   Compute the tangent vector to the input path.  The derivative approximation is described in deriv().
 637//   The returns vectors will be normalized to length 1.  If any derivatives are zero then
 638//   the function fails with an error.  If you set `uniform` to false then the sampling is
 639//   assumed to be non-uniform and the derivative is computed with adjustments to produce corrected
 640//   values.
 641// Arguments:
 642//   path = path of any dimension or a 1-region
 643//   closed = set to true of the path is closed.  Default: false
 644//   uniform = set to false to correct for non-uniform sampling.  Default: true
 645// Example(2D): A shape with non-uniform sampling gives distorted derivatives that may be undesirable.  Note that derivatives tilt towards the long edges of the rectangle.  
 646//   rect = square([10,3]);
 647//   tangents = path_tangents(rect,closed=true);
 648//   stroke(rect,closed=true, width=0.25);
 649//   color("purple")
 650//       for(i=[0:len(tangents)-1])
 651//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 652// Example(2D): Setting uniform to false corrects the distorted derivatives for this example:
 653//   rect = square([10,3]);
 654//   tangents = path_tangents(rect,closed=true,uniform=false);
 655//   stroke(rect,closed=true, width=0.25);
 656//   color("purple")
 657//       for(i=[0:len(tangents)-1])
 658//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 659function path_tangents(path, closed, uniform=true) =
 660    is_1region(path) ? path_tangents(path[0], default(closed,true), uniform) :
 661    let(closed=default(closed,false))
 662    assert(is_bool(closed))
 663    assert(is_path(path))
 664    !uniform ? [for(t=deriv(path,closed=closed, h=path_segment_lengths(path,closed))) unit(t)]
 665             : [for(t=deriv(path,closed=closed)) unit(t)];
 666
 667
 668// Function: path_normals()
 669// Synopsis: Returns normal vectors for each point along a path.
 670// Topics: Paths
 671// See Also: path_tangents()
 672// Usage:
 673//   norms = path_normals(path, [tangents], [closed]);
 674// Description:
 675//   Compute the normal vector to the input path.  This vector is perpendicular to the
 676//   path tangent and lies in the plane of the curve.  For 3d paths we define the plane of the curve
 677//   at path point i to be the plane defined by point i and its two neighbors.  At the endpoints of open paths
 678//   we use the three end points.  For 3d paths the computed normal is the one lying in this plane that points
 679//   towards the center of curvature at that path point.  For 2d paths, which lie in the xy plane, the normal
 680//   is the path pointing to the right of the direction the path is traveling.  If points are collinear then
 681//   a 3d path has no center of curvature, and hence the 
 682//   normal is not uniquely defined.  In this case the function issues an error.
 683//   For 2d paths the plane is always defined so the normal fails to exist only
 684//   when the derivative is zero (in the case of repeated points).
 685// Arguments:
 686//   path = 2D or 3D path or a 1-region
 687//   tangents = path tangents optionally supplied
 688//   closed = if true path is treated as a polygon.  Default: false
 689function path_normals(path, tangents, closed) =
 690    is_1region(path) ? path_normals(path[0], tangents, default(closed,true)) :
 691    let(closed=default(closed,false))
 692    assert(is_path(path,[2,3]))
 693    assert(is_bool(closed))
 694    let(
 695         tangents = default(tangents, path_tangents(path,closed)),
 696         dim=len(path[0])
 697    )
 698    assert(is_path(tangents) && len(tangents[0])==dim,"Dimensions of path and tangents must match")
 699    [
 700     for(i=idx(path))
 701         let(
 702             pts = i==0 ? (closed? select(path,-1,1) : select(path,0,2))
 703                 : i==len(path)-1 ? (closed? select(path,i-1,i+1) : select(path,i-2,i))
 704                 : select(path,i-1,i+1)
 705        )
 706        dim == 2 ? [tangents[i].y,-tangents[i].x]
 707                 : let( v=cross(cross(pts[1]-pts[0], pts[2]-pts[0]),tangents[i]))
 708                   assert(norm(v)>EPSILON, "3D path contains collinear points")
 709                   unit(v)
 710    ];
 711
 712
 713// Function: path_curvature()
 714// Synopsis: Returns the estimated numerical curvature of the path.
 715// Topics: Paths
 716// See Also: path_tangents(), path_normals(), path_torsion()
 717// Usage:
 718//   curvs = path_curvature(path, [closed]);
 719// Description:
 720//   Numerically estimate the curvature of the path (in any dimension).
 721// Arguments:
 722//   path = path in any dimension or a 1-region
 723//   closed = if true then treat the path as a polygon.  Default: false
 724function path_curvature(path, closed) =
 725    is_1region(path) ? path_curvature(path[0], default(closed,true)) :
 726    let(closed=default(closed,false))
 727    assert(is_bool(closed))
 728    assert(is_path(path))
 729    let( 
 730        d1 = deriv(path, closed=closed),
 731        d2 = deriv2(path, closed=closed)
 732    ) [
 733        for(i=idx(path))
 734        sqrt(
 735            sqr(norm(d1[i])*norm(d2[i])) -
 736            sqr(d1[i]*d2[i])
 737        ) / pow(norm(d1[i]),3)
 738    ];
 739
 740
 741// Function: path_torsion()
 742// Synopsis: Returns the estimated numerical torsion of the path.
 743// Topics: Paths
 744// See Also: path_tangents(), path_normals(), path_curvature()
 745// Usage:
 746//   torsions = path_torsion(path, [closed]);
 747// Description:
 748//   Numerically estimate the torsion of a 3d path.
 749// Arguments:
 750//   path = 3D path
 751//   closed = if true then treat path as a polygon.  Default: false
 752function path_torsion(path, closed=false) =
 753    assert(is_path(path,3), "Input path must be a 3d path")
 754    assert(is_bool(closed))
 755    let(
 756        d1 = deriv(path,closed=closed),
 757        d2 = deriv2(path,closed=closed),
 758        d3 = deriv3(path,closed=closed)
 759    ) [
 760        for (i=idx(path)) let(
 761            crossterm = cross(d1[i],d2[i])
 762        ) crossterm * d3[i] / sqr(norm(crossterm))
 763    ];
 764
 765
 766// Section: Breaking paths up into subpaths
 767
 768
 769
 770// Function: path_cut()
 771// Synopsis: Cuts a path into subpaths at various points.
 772// SynTags: PathList
 773// Topics: Paths, Path Subdivision
 774// See Also: split_path_at_self_crossings(), path_cut_points()
 775// Usage:
 776//   path_list = path_cut(path, cutdist, [closed]);
 777// Description:
 778//   Given a list of distances in `cutdist`, cut the path into
 779//   subpaths at those lengths, returning a list of paths.
 780//   If the input path is closed then the final path will include the
 781//   original starting point.  The list of cut distances must be
 782//   in ascending order and should not include the endpoints: 0 
 783//   or len(path).  If you repeat a distance you will get an
 784//   empty list in that position in the output.  If you give an
 785//   empty cutdist array you will get the input path as output
 786//   (without the final vertex doubled in the case of a closed path).
 787// Arguments:
 788//   path = path of any dimension or a 1-region
 789//   cutdist = Distance or list of distances where path is cut
 790//   closed = If true, treat the path as a closed polygon.  Default: false
 791// Example(2D,NoAxes):
 792//   path = circle(d=100);
 793//   segs = path_cut(path, [50, 200], closed=true);
 794//   rainbow(segs) stroke($item, endcaps="butt", width=3);
 795function path_cut(path,cutdist,closed) =
 796  is_num(cutdist) ? path_cut(path,[cutdist],closed) :
 797  is_1region(path) ? path_cut(path[0], cutdist, default(closed,true)):
 798  let(closed=default(closed,false))
 799  assert(is_bool(closed))
 800  assert(is_vector(cutdist))
 801  assert(last(cutdist)<path_length(path,closed=closed),"Cut distances must be smaller than the path length")
 802  assert(cutdist[0]>0, "Cut distances must be strictly positive")
 803  let(
 804      cutlist = path_cut_points(path,cutdist,closed=closed)
 805  )
 806  _path_cut_getpaths(path, cutlist, closed);
 807
 808
 809function _path_cut_getpaths(path, cutlist, closed) =
 810  let(
 811      cuts = len(cutlist)
 812  )
 813  [
 814      [ each list_head(path,cutlist[0][1]-1),
 815        if (!approx(cutlist[0][0], path[cutlist[0][1]-1])) cutlist[0][0]
 816      ],
 817      for(i=[0:1:cuts-2])
 818          cutlist[i][0]==cutlist[i+1][0] && cutlist[i][1]==cutlist[i+1][1] ? []
 819          :
 820          [ if (!approx(cutlist[i][0], select(path,cutlist[i][1]))) cutlist[i][0],
 821            each slice(path, cutlist[i][1], cutlist[i+1][1]-1),
 822            if (!approx(cutlist[i+1][0], select(path,cutlist[i+1][1]-1))) cutlist[i+1][0],
 823          ],
 824      [
 825        if (!approx(cutlist[cuts-1][0], select(path,cutlist[cuts-1][1]))) cutlist[cuts-1][0],
 826        each select(path,cutlist[cuts-1][1],closed ? 0 : -1)
 827      ]
 828  ];
 829
 830
 831
 832// Function: path_cut_points()
 833// Synopsis: Returns a list of cut points at a list of distances from the first point in a path.
 834// Topics: Paths, Path Subdivision
 835// See Also: path_cut(), split_path_at_self_crossings()
 836// Usage:
 837//   cuts = path_cut_points(path, cutdist, [closed=], [direction=]);
 838//
 839// Description:
 840//   Cuts a path at a list of distances from the first point in the path.  Returns a list of the cut
 841//   points and indices of the next point in the path after that point.  So for example, a return
 842//   value entry of [[2,3], 5] means that the cut point was [2,3] and the next point on the path after
 843//   this point is path[5].  If the path is too short then path_cut_points returns undef.  If you set
 844//   `direction` to true then `path_cut_points` will also return the tangent vector to the path and a normal
 845//   vector to the path.  It tries to find a normal vector that is coplanar to the path near the cut
 846//   point.  If this fails it will return a normal vector parallel to the xy plane.  The output with
 847//   direction vectors will be `[point, next_index, tangent, normal]`.
 848//   .
 849//   If you give the very last point of the path as a cut point then the returned index will be
 850//   one larger than the last index (so it will not be a valid index).  If you use the closed
 851//   option then the returned index will be equal to the path length for cuts along the closing
 852//   path segment, and if you give a point equal to the path length you will get an
 853//   index of len(path)+1 for the index.  
 854//
 855// Arguments:
 856//   path = path to cut
 857//   cutdist = distances where the path should be cut (a list) or a scalar single distance
 858//   ---
 859//   closed = set to true if the curve is closed.  Default: false
 860//   direction = set to true to return direction vectors.  Default: false
 861//
 862// Example(NORENDER):
 863//   square=[[0,0],[1,0],[1,1],[0,1]];
 864//   path_cut_points(square, [.5,1.5,2.5]);   // Returns [[[0.5, 0], 1], [[1, 0.5], 2], [[0.5, 1], 3]]
 865//   path_cut_points(square, [0,1,2,3]);      // Returns [[[0, 0], 1], [[1, 0], 2], [[1, 1], 3], [[0, 1], 4]]
 866//   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]]
 867//   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]
 868function path_cut_points(path, cutdist, closed=false, direction=false) =
 869    let(long_enough = len(path) >= (closed ? 3 : 2))
 870    assert(long_enough,len(path)<2 ? "Two points needed to define a path" : "Closed path must include three points")
 871    is_num(cutdist) ? path_cut_points(path, [cutdist],closed, direction)[0] :
 872    assert(is_vector(cutdist))
 873    assert(is_increasing(cutdist), "Cut distances must be an increasing list")
 874    let(cuts = path_cut_points_recurse(path,cutdist,closed))
 875    !direction
 876       ? cuts
 877       : let(
 878             dir = _path_cuts_dir(path, cuts, closed),
 879             normals = _path_cuts_normals(path, cuts, dir, closed)
 880         )
 881         hstack(cuts, list_to_matrix(dir,1), list_to_matrix(normals,1));
 882
 883// Main recursive path cut function
 884function path_cut_points_recurse(path, dists, closed=false, pind=0, dtotal=0, dind=0, result=[]) =
 885    dind == len(dists) ? result :
 886    let(
 887        lastpt = len(result)==0? [] : last(result)[0],       // location of last cut point
 888        dpartial = len(result)==0? 0 : norm(lastpt-select(path,pind)),  // remaining length in segment
 889        nextpoint = dists[dind] < dpartial+dtotal  // Do we have enough length left on the current segment?
 890           ? [lerp(lastpt,select(path,pind),(dists[dind]-dtotal)/dpartial),pind] 
 891           : _path_cut_single(path, dists[dind]-dtotal-dpartial, closed, pind)
 892    ) 
 893    path_cut_points_recurse(path, dists, closed, nextpoint[1], dists[dind],dind+1, concat(result, [nextpoint]));
 894
 895
 896// Search for a single cut point in the path
 897function _path_cut_single(path, dist, closed=false, ind=0, eps=1e-7) =
 898    // If we get to the very end of the path (ind is last point or wraparound for closed case) then
 899    // check if we are within epsilon of the final path point.  If not we're out of path, so we fail
 900    ind==len(path)-(closed?0:1) ?
 901       assert(dist<eps,"Path is too short for specified cut distance")
 902       [select(path,ind),ind+1]
 903    :let(d = norm(path[ind]-select(path,ind+1))) d > dist ?
 904        [lerp(path[ind],select(path,ind+1),dist/d), ind+1] :
 905        _path_cut_single(path, dist-d,closed, ind+1, eps);
 906
 907// Find normal directions to the path, coplanar to local part of the path
 908// Or return a vector parallel to the x-y plane if the above fails
 909function _path_cuts_normals(path, cuts, dirs, closed=false) =
 910    [for(i=[0:len(cuts)-1])
 911        len(path[0])==2? [-dirs[i].y, dirs[i].x]
 912          : 
 913            let(
 914                plane = len(path)<3 ? undef :
 915                let(start = max(min(cuts[i][1],len(path)-1),2)) _path_plane(path, start, start-2)
 916            )
 917            plane==undef?
 918                ( dirs[i].x==0 && dirs[i].y==0 ? [1,0,0]  // If it's z direction return x vector
 919                                               : unit([-dirs[i].y, dirs[i].x,0])) // otherwise perpendicular to projection
 920                : unit(cross(dirs[i],cross(plane[0],plane[1])))
 921    ];
 922
 923// Scan from the specified point (ind) to find a noncoplanar triple to use
 924// to define the plane of the path.
 925function _path_plane(path, ind, i,closed) =
 926    i<(closed?-1:0) ? undef :
 927    !is_collinear(path[ind],path[ind-1], select(path,i))?
 928        [select(path,i)-path[ind-1],path[ind]-path[ind-1]] :
 929        _path_plane(path, ind, i-1);
 930
 931// Find the direction of the path at the cut points
 932function _path_cuts_dir(path, cuts, closed=false, eps=1e-2) =
 933    [for(ind=[0:len(cuts)-1])
 934        let(
 935            zeros = path[0]*0,
 936            nextind = cuts[ind][1],
 937            nextpath = unit(select(path, nextind+1)-select(path, nextind),zeros),
 938            thispath = unit(select(path, nextind) - select(path,nextind-1),zeros),
 939            lastpath = unit(select(path,nextind-1) - select(path, nextind-2),zeros),
 940            nextdir =
 941                nextind==len(path) && !closed? lastpath :
 942                (nextind<=len(path)-2 || closed) && approx(cuts[ind][0], path[nextind],eps)
 943                   ? unit(nextpath+thispath)
 944              : (nextind>1 || closed) && approx(cuts[ind][0],select(path,nextind-1),eps)
 945                   ? unit(thispath+lastpath)
 946              :  thispath
 947        ) nextdir
 948    ];
 949
 950
 951// internal function
 952// converts pathcut output form to a [segment, u]
 953// form list that works withi path_select
 954function _cut_to_seg_u_form(pathcut, path, closed) =
 955  let(lastind = len(path) - (closed?0:1))
 956  [for(entry=pathcut)
 957    entry[1] > lastind ? [lastind,0] :
 958    let(
 959        a = path[entry[1]-1],
 960        b = path[entry[1]],
 961        c = entry[0],
 962        i = max_index(v_abs(b-a)),
 963        factor = (c[i]-a[i])/(b[i]-a[i])
 964    )
 965    [entry[1]-1,factor]
 966  ];
 967
 968
 969
 970// Function: split_path_at_self_crossings()
 971// Synopsis: Split a 2D path wherever it crosses itself.
 972// SynTags: PathList
 973// Topics: Paths, Path Subdivision
 974// See Also: path_cut(), path_cut_points()
 975// Usage:
 976//   paths = split_path_at_self_crossings(path, [closed], [eps]);
 977// Description:
 978//   Splits a 2D path into sub-paths wherever the original path crosses itself.
 979//   Splits may occur mid-segment, so new vertices will be created at the intersection points.
 980//   Returns a list of the resulting subpaths.  
 981// Arguments:
 982//   path = A 2D path or a 1-region.
 983//   closed = If true, treat path as a closed polygon.  Default: true
 984//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 985// Example(2D,NoAxes):
 986//   path = [ [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100] ];
 987//   paths = split_path_at_self_crossings(path);
 988//   rainbow(paths) stroke($item, closed=false, width=3);
 989function split_path_at_self_crossings(path, closed=true, eps=EPSILON) =
 990    let(path = force_path(path))
 991    assert(is_path(path,2), "Must give a 2D path")
 992    assert(is_bool(closed))
 993    let(
 994        path = list_unwrap(path, eps=eps),
 995        isects = deduplicate(
 996            eps=eps,
 997            concat(
 998                [[0, 0]],
 999                sort([
1000                    for (
1001                        a = _path_self_intersections(path, closed=closed, eps=eps),
1002                        ss = [ [a[1],a[2]], [a[3],a[4]] ]
1003                    ) if (ss[0] != undef) ss
1004                ]),
1005                [[len(path)-(closed?1:2), 1]]
1006            )
1007        )
1008    ) [
1009        for (p = pair(isects))
1010            let(
1011                s1 = p[0][0],
1012                u1 = p[0][1],
1013                s2 = p[1][0],
1014                u2 = p[1][1],
1015                section = _path_select(path, s1, u1, s2, u2, closed=closed),
1016                outpath = deduplicate(eps=eps, section)
1017            )
1018            if (len(outpath)>1) outpath
1019    ];
1020
1021
1022function _tag_self_crossing_subpaths(path, nonzero, closed=true, eps=EPSILON) =
1023    let(
1024        subpaths = split_path_at_self_crossings(
1025            path, closed=true, eps=eps
1026        )
1027    ) [
1028        for (subpath = subpaths) let(
1029            seg = select(subpath,0,1),
1030            mp = mean(seg),
1031            n = line_normal(seg) / 2048,
1032            p1 = mp + n,
1033            p2 = mp - n,
1034            p1in = point_in_polygon(p1, path, nonzero=nonzero) >= 0,
1035            p2in = point_in_polygon(p2, path, nonzero=nonzero) >= 0,
1036            tag = (p1in && p2in)? "I" : "O"
1037        ) [tag, subpath]
1038    ];
1039
1040
1041// Function: polygon_parts()
1042// Synopsis: Parses a self-intersecting polygon into a list of non-intersecting polygons.
1043// SynTags: PathList
1044// Topics: Paths, Polygons
1045// See Also: split_path_at_self_crossings(), path_cut(), path_cut_points()
1046// Usage:
1047//   splitpolys = polygon_parts(poly, [nonzero], [eps]);
1048// Description:
1049//   Given a possibly self-intersecting 2d polygon, constructs a representation of the original polygon as a list of
1050//   non-intersecting simple polygons.  If nonzero is set to true then it uses the nonzero method for defining polygon membership.
1051//   For simple cases, such as the pentagram, this will produce the outer perimeter of a self-intersecting polygon.  
1052// Arguments:
1053//   poly = a 2D polygon or 1-region
1054//   nonzero = If true use the nonzero method for checking if a point is in a polygon.  Otherwise use the even-odd method.  Default: false
1055//   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1056// Example(2D,NoAxes):  This cross-crossing polygon breaks up into its 3 components (regardless of the value of nonzero).
1057//   poly = [
1058//       [-100,100], [0,-50], [100,100],
1059//       [100,-100], [0,50], [-100,-100]
1060//   ];
1061//   splitpolys = polygon_parts(poly);
1062//   rainbow(splitpolys) stroke($item, closed=true, width=3);
1063// Example(2D,NoAxes): With nonzero=false you get even-odd mode which matches OpenSCAD, so the pentagram breaks apart into its five points.
1064//   pentagram = turtle(["move",100,"left",144], repeat=4);
1065//   left(100)polygon(pentagram);
1066//   rainbow(polygon_parts(pentagram,nonzero=false))
1067//     stroke($item,closed=true,width=2.5);
1068// 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.
1069//   pentagram = turtle(["move",100,"left",144], repeat=4);
1070//   outside = polygon_parts(pentagram,nonzero=true);
1071//   left(100)region(outside);
1072//   rainbow(outside)
1073//     stroke($item,closed=true,width=2.5);
1074// Example(2D,NoAxes): 
1075//   N=12;
1076//   ang=360/N;
1077//   sr=10;
1078//   poly = turtle(["angle", 90+ang/2,
1079//                  "move", sr, "left",
1080//                  "move", 2*sr*sin(ang/2), "left",
1081//                  "repeat", 4,
1082//                     ["move", 2*sr, "left",
1083//                      "move", 2*sr*sin(ang/2), "left"],
1084//                  "move", sr]);
1085//   stroke(poly, width=.3);
1086//   right(20)rainbow(polygon_parts(poly)) polygon($item);
1087// Example(2D,NoAxes): overlapping poly segments disappear
1088//   poly = [[0,0], [10,0], [10,10], [0,10],[0,20], [20,10],[10,10], [0,10],[0,0]];
1089//   stroke(poly,width=0.3);
1090//   right(22)stroke(polygon_parts(poly)[0], width=0.3, closed=true);
1091// Example(2D,NoAxes): Poly segments disappear outside as well
1092//   poly = turtle(["repeat", 3, ["move", 17, "left", "move", 10, "left", "move", 7, "left", "move", 10, "left"]]);
1093//   back(2)stroke(poly,width=.5);
1094//   fwd(12)rainbow(polygon_parts(poly)) stroke($item, closed=true, width=0.5);
1095// Example(2D,NoAxes):  This shape has six components
1096//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 17, "left"]]);
1097//   polygon(poly);
1098//   right(22)rainbow(polygon_parts(poly)) polygon($item);
1099// Example(2D,NoAxes): When the loops of the shape overlap then nonzero gives a different result than the even-odd method.
1100//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 10, "left"]]);
1101//   polygon(poly);
1102//   right(27)rainbow(polygon_parts(poly)) polygon($item);
1103//   move([16,-14])rainbow(polygon_parts(poly,nonzero=true)) polygon($item);
1104function polygon_parts(poly, nonzero=false, eps=EPSILON) =
1105    let(poly = force_path(poly))
1106    assert(is_path(poly,2), "Must give 2D polygon")
1107    assert(is_bool(nonzero))    
1108    let(
1109        poly = list_unwrap(poly, eps=eps),
1110        tagged = _tag_self_crossing_subpaths(poly, nonzero=nonzero, closed=true, eps=eps),
1111        kept = [for (sub = tagged) if(sub[0] == "O") sub[1]],
1112        outregion = _assemble_path_fragments(kept, eps=eps)
1113    ) outregion;
1114
1115
1116function _extreme_angle_fragment(seg, fragments, rightmost=true, eps=EPSILON) =
1117    !fragments? [undef, []] :
1118    let(
1119        delta = seg[1] - seg[0],
1120        segang = atan2(delta.y,delta.x),
1121        frags = [
1122            for (i = idx(fragments)) let(
1123                fragment = fragments[i],
1124                fwdmatch = approx(seg[1], fragment[0], eps=eps),
1125                bakmatch =  approx(seg[1], last(fragment), eps=eps)
1126            ) [
1127                fwdmatch,
1128                bakmatch,
1129                bakmatch? reverse(fragment) : fragment
1130            ]
1131        ],
1132        angs = [
1133            for (frag = frags)
1134                (frag[0] || frag[1])? let(
1135                    delta2 = frag[2][1] - frag[2][0],
1136                    segang2 = atan2(delta2.y, delta2.x)
1137                ) modang(segang2 - segang) : (
1138                    rightmost? 999 : -999
1139                )
1140        ],
1141        fi = rightmost? min_index(angs) : max_index(angs)
1142    ) abs(angs[fi]) > 360? [undef, fragments] : let(
1143        remainder = [for (i=idx(fragments)) if (i!=fi) fragments[i]],
1144        frag = frags[fi],
1145        foundfrag = frag[2]
1146    ) [foundfrag, remainder];
1147
1148
1149/// Internal Function: _assemble_a_path_from_fragments()
1150/// Usage:
1151///   _assemble_a_path_from_fragments(subpaths);
1152/// Description:
1153///   Given a list of paths, assembles them together into one complete closed polygon path, and
1154///   remainder fragments.  Returns [PATH, FRAGMENTS] where FRAGMENTS is the list of remaining
1155///   unused path fragments.
1156/// Arguments:
1157///   fragments = List of paths to be assembled into complete polygons.
1158///   rightmost = If true, assemble paths using rightmost turns. Leftmost if false.
1159///   startfrag = The fragment to start with.  Default: 0
1160///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1161function _assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0, eps=EPSILON) =
1162    len(fragments)==0? [[],[]] :
1163    len(fragments)==1? [fragments[0],[]] :
1164    let(
1165        path = fragments[startfrag],
1166        newfrags = [for (i=idx(fragments)) if (i!=startfrag) fragments[i]]
1167    ) are_ends_equal(path, eps=eps)? (
1168        // starting fragment is already closed
1169        [path, newfrags]
1170    ) : let(
1171        // Find rightmost/leftmost continuation fragment
1172        seg = select(path,-2,-1),
1173        extrema = _extreme_angle_fragment(seg=seg, fragments=newfrags, rightmost=rightmost, eps=eps),
1174        foundfrag = extrema[0],
1175        remainder = extrema[1]
1176    ) is_undef(foundfrag)? (
1177        // No remaining fragments connect!  INCOMPLETE PATH!
1178        // Treat it as complete.
1179        [path, remainder]
1180    ) : are_ends_equal(foundfrag, eps=eps)? (
1181        // Found fragment is already closed
1182        [foundfrag, concat([path], remainder)]
1183    ) : let(
1184        fragend = last(foundfrag),
1185        hits = [for (i = idx(path,e=-2)) if(approx(path[i],fragend,eps=eps)) i]
1186    ) hits? (
1187        let(
1188            // Found fragment intersects with initial path
1189            hitidx = last(hits),
1190            newpath = list_head(path,hitidx),
1191            newfrags = concat(len(newpath)>1? [newpath] : [], remainder),
1192            outpath = concat(slice(path,hitidx,-2), foundfrag)
1193        )
1194        [outpath, newfrags]
1195    ) : let(
1196        // Path still incomplete.  Continue building it.
1197        newpath = concat(path, list_tail(foundfrag)),
1198        newfrags = concat([newpath], remainder)
1199    )
1200    _assemble_a_path_from_fragments(
1201        fragments=newfrags,
1202        rightmost=rightmost,
1203        eps=eps
1204    );
1205
1206
1207/// Internal Function: _assemble_path_fragments()
1208/// Usage:
1209///   _assemble_path_fragments(subpaths);
1210/// Description:
1211///   Given a list of paths, assembles them together into complete closed polygon paths if it can.
1212///   Polygons with area < eps will be discarded and not returned.  
1213/// Arguments:
1214///   fragments = List of paths to be assembled into complete polygons.
1215///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1216function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
1217    len(fragments)==0? _finished :
1218    let(
1219        minxidx = min_index([
1220            for (frag=fragments) min(column(frag,0))
1221        ]),
1222        result_l = _assemble_a_path_from_fragments(
1223            fragments=fragments,
1224            startfrag=minxidx,
1225            rightmost=false,
1226            eps=eps
1227        ),
1228        result_r = _assemble_a_path_from_fragments(
1229            fragments=fragments,
1230            startfrag=minxidx,
1231            rightmost=true,
1232            eps=eps
1233        ),
1234        l_area = abs(polygon_area(result_l[0])),
1235        r_area = abs(polygon_area(result_r[0])),
1236        result = l_area < r_area? result_l : result_r,
1237        newpath = list_unwrap(result[0]),
1238        remainder = result[1],
1239        finished = min(l_area,r_area)<eps ? _finished : concat(_finished, [newpath])
1240    ) _assemble_path_fragments(
1241        fragments=remainder,
1242        eps=eps,
1243        _finished=finished
1244    );
1245
1246
1247
1248// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap