1//////////////////////////////////////////////////////////////////////
   2// LibFile: regions.scad
   3//   This file provides 2D boolean set operations on polygons, where you can
   4//   compute, for example, the intersection or union of the shape defined by point lists, producing
   5//   a new point list.  Of course, such operations may produce shapes with multiple
   6//   components.  To handle that, we use "regions" which are lists of paths representing the polygons.
   7//   In addition to set operations, you can calculate offsets, determine whether a point is in a
   8//   region and you can decompose a region into parts.  
   9// Includes:
  10//   include <BOSL2/std.scad>
  11// FileGroup: Advanced Modeling
  12// FileSummary: Offsets and boolean geometry of 2D paths and regions.
  13// FileFootnotes: STD=Included in std.scad
  14//////////////////////////////////////////////////////////////////////
  15
  16
  17// CommonCode:
  18//   include <BOSL2/rounding.scad>
  19
  20
  21// Section: Regions
  22//   A region is a list of polygons meeting these conditions:
  23//   .
  24//   - Every polygon on the list is simple, meaning it does not intersect itself
  25//   - Two polygons on the list do not cross each other
  26//   - A vertex of one polygon never meets the edge of another one except at a vertex
  27//   .
  28//   Note that this means vertex-vertex touching between two polygons is acceptable
  29//   to define a region.  Note, however, that regions with vertex-vertex contact usually
  30//   cannot be rendered with CGAL.  See {{is_valid_region()}} for examples of valid regions and
  31//   lists of polygons that are not regions.  Note that {{is_region_simple()}} will identify
  32//   regions with no polygon intersections at all, which should render successfully witih CGAL.  
  33//   .
  34//   The actual geometry of the region is defined by XORing together
  35//   all of the polygons in the list.  This may sound obscure, but it simply means that nested
  36//   boundaries make rings in the obvious fashion, and non-nested shapes simply union together.
  37//   Checking that a list of polygons is a valid region, meaning that it satisfies all of the conditions
  38//   above, can be a time consuming test, so it is not done automatically.  It is your responsibility to ensure that your regions are
  39//   compliant.  You can construct regions by making a suitable list of polygons, or by using
  40//   set operation function such as union() or difference(), which all acccept polygons, as
  41//   well as regions, as their inputs.  And if you must you can clean up an ill-formed region using make_region(),
  42//   which will break up self-intersecting polygons and polygons that cross each other.  
  43
  44
  45// Function: is_region()
  46// Usage:
  47//   bool = is_region(x);
  48// Description:
  49//   Returns true if the given item looks like a region.  A region is a list of non-crossing simple polygons.  This test just checks
  50//   that the argument is a list whose first entry is a path.  
  51function is_region(x) = is_list(x) && is_path(x.x);
  52
  53
  54// Function: is_valid_region()
  55// Usage:
  56//   bool = is_valid_region(region, [eps]);
  57// Description:
  58//   Returns true if the input is a valid region, meaning that it is a list of simple polygons whose segments do not cross each other.
  59//   This test can be time consuming with regions that contain many points.
  60//   It differs from `is_region()` which simply checks that the object is a list whose first entry is a path
  61//   because it searches all the list polygons for any self-intersections or intersections with each other.  
  62//   Will also return true if given a single simple polygon.  Use {{make_region()}} to convert sets of self-intersecting polygons into
  63//   a region.  
  64// Arguments:
  65//   region = region to check
  66//   eps = tolerance for geometric comparisons.  Default: `EPSILON` = 1e-9
  67// Example(2D,NoAxes):  In all of the examples each polygon in the region appears in a different color.  Two non-intersecting squares make a valid region.
  68//   region = [square(10), right(11,square(8))];
  69//   rainbow(region)stroke($item, width=.2,closed=true);
  70//   back(11)text(is_valid_region(region) ? "region" : "non-region", size=2);
  71// Example(2D,NoAxes):  Nested squares form a region
  72//   region = [for(i=[3:2:10]) square(i,center=true)];
  73//   rainbow(region)stroke($item, width=.2,closed=true);
  74//   back(6)text(is_valid_region(region) ? "region" : "non-region", size=2,halign="center");
  75// Example(2D,NoAxes):  Also a region:
  76//   region= [square(10,center=true), square(5,center=true), right(10,square(7))];
  77//   rainbow(region)stroke($item, width=.2,closed=true);
  78//   back(8)text(is_valid_region(region) ? "region" : "non-region", size=2);
  79// Example(2D,NoAxes):  The squares cross each other, so not a region
  80//   object = [square(10), move([8,8], square(8))];
  81//   rainbow(object)stroke($item, width=.2,closed=true);
  82//   back(17)text(is_valid_region(object) ? "region" : "non-region", size=2);
  83// Example(2D,NoAxes): A union is one way to fix the above example and get a region.  (Note that union is run here on two simple polygons, which are valid regions themselves and hence acceptable inputs to union.
  84//   region = union([square(10), move([8,8], square(8))]);
  85//   rainbow(region)stroke($item, width=.25,closed=true);
  86//   back(12)text(is_valid_region(region) ? "region" : "non-region", size=2);
  87// Example(2D,NoAxes):  Not a region due to a self-intersecting (non-simple) hourglass polygon
  88//   object = [move([-2,-2],square(14)), [[0,0],[10,0],[0,10],[10,10]]];
  89//   rainbow(object)stroke($item, width=.2,closed=true);
  90//   move([-1.5,13])text(is_valid_region(object) ? "region" : "non-region", size=2);
  91// Example(2D,NoAxes):  Breaking hourglass in half fixes it.  Now it's a region:
  92//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,5]], [[5,5],[0,10],[10,10]]];
  93//   rainbow(region)stroke($item, width=.2,closed=true);
  94// Example(2D,NoAxes):  A single polygon corner touches an edge, so not a region:
  95//   object = [[[-10,0], [-10,10], [20,10], [20,-20], [-10,-20],
  96//              [-10,-10], [0,0], [10,-10], [10,0]]];
  97//   rainbow(object)stroke($item, width=.3,closed=true);
  98//   move([-4,12])text(is_valid_region(object) ? "region" : "non-region", size=3);
  99// Example(2D,NoAxes):  Corners touch in the same polygon, so the polygon is not simple and the object is not a region.
 100//   object = [[[0,0],[10,0],[10,10],[-10,10],[-10,0],[0,0],[-5,5],[5,5]]];
 101//   rainbow(object)stroke($item, width=.3,closed=true);
 102//   move([-10,12])text(is_valid_region(object) ? "region" : "non-region", size=3);
 103// Example(2D,NoAxes):  The shape above as a valid region with two polygons:
 104//   region = [  [[0,0],[10,0],[10,10],[-10,10],[-10,0]],
 105//               [[0,0],[5,5],[-5,5]]  ];
 106//   rainbow(region)stroke($item, width=.3,closed=true);
 107//   move([-5.5,12])text(is_valid_region(region) ? "region" : "non-region", size=3);
 108// Example(2D,NoAxes):  As with the "broken" hourglass, Touching at corners is OK.  This is a region.
 109//   region = [square(10), move([10,10], square(8))];
 110//   rainbow(region)stroke($item, width=.25,closed=true);
 111//   back(12)text(is_valid_region(region) ? "region" : "non-region", size=2);
 112// Example(2D,NoAxes): These two squares share part of an edge, hence not a region
 113//   object = [square(10), move([10,2], square(7))];
 114//   stroke(object[0], width=0.2,closed=true);
 115//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 116//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 117// Example(2D,NoAxes): These two squares share a full edge, hence not a region
 118//   object = [square(10), right(10, square(10))];
 119//   stroke(object[0], width=0.2,closed=true);
 120//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 121//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 122// Example(2D,NoAxes): Sharing on edge on the inside, also not a regionn
 123//   object = [square(10), [[0,0], [2,2],[2,8],[0,10]]];
 124//   stroke(object[0], width=0.2,closed=true);
 125//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 126//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 127// Example(2D,NoAxes): Crossing at vertices is also bad
 128//   object = [square(10), [[10,0],[0,10],[8,13],[13,8]]];
 129//   rainbow(object)stroke($item, width=.2,closed=true);
 130//   back(14)text(is_valid_region(object) ? "region" : "non-region", size=2);
 131// Example(2D,NoAxes): One polygon touches another in the middle of an edge
 132//   object = [square(10), [[10,5],[15,0],[15,10]]];
 133//   rainbow(object)stroke($item, width=.2,closed=true);
 134//   back(11)text(is_valid_region(object) ? "region" : "non-region", size=2);
 135// Example(2D,NoAxes): The polygon touches the side, but the side has a vertex at the contact point so this is a region
 136//   poly1 = [ each square(30,center=true), [15,0]];
 137//   poly2 = right(10,circle(5,$fn=4));
 138//   poly3 = left(0,circle(5,$fn=4));
 139//   poly4 = move([0,-8],square([10,3]));
 140//   region = [poly1,poly2,poly3,poly4];
 141//   rainbow(region)stroke($item, width=.25,closed=true);
 142//   move([-5,16.5])text(is_valid_region(region) ? "region" : "non-region", size=3);
 143//   color("black")move_copies(region[0]) circle(r=.4);
 144// Example(2D,NoAxes): The polygon touches the side, but not at a vertex so this is not a region
 145//   poly1 = fwd(4,[ each square(30,center=true), [15,0]]);
 146//   poly2 = right(10,circle(5,$fn=4));
 147//   poly3 = left(0,circle(5,$fn=4));
 148//   poly4 = move([0,-8],square([10,3]));
 149//   object = [poly1,poly2,poly3,poly4];
 150//   rainbow(object)stroke($item, width=.25,closed=true);
 151//   move([-9,12.5])text(is_valid_region(object) ? "region" : "non-region", size=3);
 152//   color("black")move_copies(object[0]) circle(r=.4);
 153// Example(2D,NoAxes): The inner polygon touches the middle of the edges, so not a region
 154//   poly1 = square(20,center=true);
 155//   poly2 = circle(10,$fn=8);
 156//   object=[poly1,poly2];
 157//   rainbow(object)stroke($item, width=.25,closed=true);
 158//   move([-10,11.4])text(is_valid_region(object) ? "region" : "non-region", size=3);
 159// Example(2D,NoAxes): The above shape made into a region using {{difference()}} now has four components that touch at corners
 160//   poly1 = square(20,center=true);
 161//   poly2 = circle(10,$fn=8);
 162//   region = difference(poly1,poly2);
 163//   rainbow(region)stroke($item, width=.25,closed=true);
 164//   move([-5,11.4])text(is_valid_region(region) ? "region" : "non-region", size=3);
 165function is_valid_region(region, eps=EPSILON) =
 166   let(region=force_region(region))
 167   assert(is_region(region), "Input is not a region")
 168   // no short paths
 169   [for(p=region) if (len(p)<3) 1] == []
 170   &&
 171   // all paths are simple
 172   [for(p=region) if (!is_path_simple(p,closed=true,eps=eps)) 1] == []
 173   &&
 174   // paths do not cross each other
 175   [for(i=[0:1:len(region)-2])
 176            if (_polygon_crosses_region(list_tail(region,i+1),region[i], eps=eps)) 1] == []
 177   &&
 178   // one path doesn't touch another in the middle of an edge
 179   [for(i=idx(region), j=idx(region))
 180       if (i!=j) for(v=region[i], edge=pair(region[j],wrap=true))
 181           let(
 182               v1 = edge[1]-edge[0],
 183               v0 = v - edge[0],
 184               t = v0*v1/(v1*v1)
 185           )
 186           if (abs(cross(v0,v1))<eps*norm(v1) && t>eps && t<1-eps) 1
 187   ]==[];
 188
 189
 190
 191// internal function:
 192// returns true if the polygon crosses the region so that part of the 
 193// polygon is inside the region and part is outside.  
 194function _polygon_crosses_region(region, poly, eps=EPSILON) =
 195    let(  
 196        subpaths = flatten(split_region_at_region_crossings(region,[poly],eps=eps)[1])
 197    )
 198    [for(path=subpaths)
 199      let(isect=
 200         [for (subpath = subpaths)
 201          let(
 202                midpt = mean([subpath[0], subpath[1]]),
 203                rel = point_in_region(midpt,region,eps=eps)
 204          )
 205          rel
 206         ])
 207       if (!all_equal(isect) || isect[0]==0) 1 ] != [];
 208
 209
 210// Function: is_region_simple()
 211// Usage:
 212//   bool = is_region_simple(region, [eps]);
 213// Description:
 214//   We extend the notion of the simple path to regions: a simple region is entirely
 215//   non-self-intersecting, meaning that it is formed from a list of simple polygons that
 216//   don't intersect each other at all&mdash;not even with corner contact points.
 217//   Regions with corner contact are valid but may fail CGAL.  Simple regions
 218//   should not create problems with CGAL.  
 219// Arguments:
 220//   region = region to check
 221//   eps = tolerance for geometric comparisons.  Default: `EPSILON` = 1e-9
 222// Example(2D,NoAxes):  Corner contact means it's not simple
 223//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,5]], [[5,5],[0,10],[10,10]]];
 224//   rainbow(region)stroke($item, width=.2,closed=true);
 225//   move([-1,13])text(is_region_simple(region) ? "simple" : "not-simple", size=2);
 226// Example(2D,NoAxes):  Moving apart the triangles makes it simple:
 227//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,4.5]], [[5,5.5],[0,10],[10,10]]];
 228//   rainbow(region)stroke($item, width=.2,closed=true);
 229//   move([1,13])text(is_region_simple(region) ? "simple" : "not-simple", size=2);
 230function is_region_simple(region, eps=EPSILON) =
 231   let(region=force_region(region))
 232   assert(is_region(region), "Input is not a region")
 233   [for(p=region) if (!is_path_simple(p,closed=true,eps=eps)) 1] == []
 234   &&
 235   [for(i=[0:1:len(region)-2])
 236       if (_region_region_intersections([region[i]], list_tail(region,i+1), eps=eps)[0][0] != []) 1
 237   ] ==[];
 238  
 239  
 240// Function: make_region()
 241// Usage:
 242//   region = make_region(polys, [nonzero], [eps]);
 243// Description:
 244//   Takes a list of polygons that may intersect themselves or cross each other 
 245//   and converts it into a properly defined region without
 246//   these defects.
 247// Arguments:
 248//   polys = list of polygons to use
 249//   nonzero = set to true to use nonzero rule for polygon membership.  Default: false
 250//   eps = Epsilon for geometric comparisons.  Default: `EPSILON` (1e-9)
 251// Example(2D,NoAxes):  The pentagram is self-intersecting, so it is not a region.  Here it becomes five triangles:
 252//   pentagram = turtle(["move",100,"left",144], repeat=4);
 253//   region = make_region(pentagram);
 254//   rainbow(region)stroke($item, width=1,closed=true);
 255// Example(2D,NoAxes):  Alternatively with the nonzero option you can get the perimeter:
 256//   pentagram = turtle(["move",100,"left",144], repeat=4);
 257//   region = make_region(pentagram,nonzero=true);
 258//   rainbow(region)stroke($item, width=1,closed=true);
 259// Example(2D,NoAxes):  Two crossing squares become two L-shaped components
 260//   region = make_region([square(10), move([5,5],square(8))]);
 261//   rainbow(region)stroke($item, width=.3,closed=true);
 262
 263function make_region(polys,nonzero=false,eps=EPSILON) =
 264     let(polys=force_region(polys))
 265     assert(is_region(polys), "Input is not a region")
 266     exclusive_or(
 267                  [for(poly=polys) each polygon_parts(poly,nonzero,eps)],
 268                  eps=eps);
 269
 270// Function: force_region()
 271// Usage:
 272//   region = force_region(poly)
 273// Description:
 274//   If the input is a polygon then return it as a region.  Otherwise return it unaltered.
 275// Arguments:
 276//   poly = polygon to turn into a region
 277function force_region(poly) = is_path(poly) ? [poly] : poly;
 278
 279
 280// Section: Turning a region into geometry
 281
 282// Module: region()
 283// Usage:
 284//   region(r, [anchor], [spin=], [cp=], [atype=]) [ATTACHMENTS];
 285// Description:
 286//   Creates the 2D polygons described by the given region or list of polygons.  This module works on
 287//   arbitrary lists of polygons that cross each other and hence do not define a valid region.  The
 288//   displayed result is the exclusive-or of the polygons listed in the input. 
 289// Arguments:
 290//   r = region to create as geometry
 291//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `"origin"`
 292//   ---
 293//   spin = Rotate this many degrees after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
 294//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 2D point.  Default: "centroid"
 295//   atype = Set to "hull" or "intersect" to select anchor type.  Default: "hull"
 296// Example(2D): Displaying a region
 297//   region([circle(d=50), square(25,center=true)]);
 298// Example(2D): Displaying a list of polygons that intersect each other, which is not a region
 299//   rgn = concat(
 300//       [for (d=[50:-10:10]) circle(d=d-5)],
 301//       [square([60,10], center=true)]
 302//   );
 303//   region(rgn);
 304module region(r, anchor="origin", spin=0, cp="centroid", atype="hull")
 305{
 306    assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"");
 307    r = force_region(r);
 308    dummy=assert(is_region(r), "Input is not a region");
 309    points = flatten(r);
 310    lengths = [for(path=r) len(path)];
 311    starts = [0,each cumsum(lengths)];
 312    paths = [for(i=idx(r)) count(s=starts[i], n=lengths[i])];
 313    attachable(anchor, spin, two_d=true, region=r, extent=atype=="hull", cp=cp){
 314      polygon(points=points, paths=paths);
 315      children();
 316    }
 317}
 318
 319
 320
 321// Section: Gometrical calculations with regions
 322
 323// Function: point_in_region()
 324// Usage:
 325//   check = point_in_region(point, region, [eps]);
 326// Description:
 327//   Tests if a point is inside, outside, or on the border of a region.  
 328//   Returns -1 if the point is outside the region.
 329//   Returns 0 if the point is on the boundary.
 330//   Returns 1 if the point lies inside the region.
 331// Arguments:
 332//   point = The point to test.
 333//   region = The region to test against, as a list of polygon paths.
 334//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 335// Example(2D,Med):  Red points are in the region.
 336//   region = [for(i=[2:4:10]) hexagon(r=i)];
 337//   color("#ff7") region(region);
 338//   for(x=[-10:10], y=[-10:10])
 339//     if (point_in_region([x,y], region)>=0)
 340//       move([x,y]) color("red") circle(0.15, $fn=12);
 341//     else
 342//       move([x,y]) color("#ddf") circle(0.1, $fn=12);
 343function point_in_region(point, region, eps=EPSILON) =
 344    let(region=force_region(region))
 345    assert(is_region(region), "Region given to point_in_region is not a region")
 346    assert(is_vector(point,2), "Point must be a 2D point in point_in_region")
 347    _point_in_region(point, region, eps);
 348
 349function _point_in_region(point, region, eps=EPSILON, i=0, cnt=0) =
 350      i >= len(region) ? ((cnt%2==1)? 1 : -1)
 351    : let(
 352          pip = point_in_polygon(point, region[i], eps=eps)
 353      )
 354      pip == 0 ? 0
 355   : _point_in_region(point, region, eps=eps, i=i+1, cnt = cnt + (pip>0? 1 : 0));
 356
 357
 358// Function: region_area()
 359// Usage:
 360//   area = region_area(region);
 361// Description:
 362//   Computes the area of the specified valid region. (If the region is invalid and has self intersections
 363//   the result is meaningless.)
 364// Arguments:
 365//   region = region whose area to compute
 366// Examples:
 367//   area = region_area([square(10), right(20,square(8))]);  // Returns 164
 368function region_area(region) =
 369  assert(is_region(region), "Input must be a region")
 370  let(
 371      parts = region_parts(region)
 372  )
 373  -sum([for(R=parts, poly=R) polygon_area(poly,signed=true)]);
 374
 375
 376
 377function _clockwise_region(r) = [for(p=r) clockwise_polygon(p)];
 378
 379// Function: are_regions_equal()
 380// Usage:
 381//    b = are_regions_equal(region1, region2, [either_winding])
 382// Description:
 383//    Returns true if the components of region1 and region2 are the same polygons (in any order). 
 384// Arguments:
 385//    region1 = first region
 386//    region2 = second region
 387//    either_winding = if true then two shapes test equal if they wind in opposite directions.  Default: false
 388function are_regions_equal(region1, region2, either_winding=false) =
 389    let(
 390        region1=force_region(region1),
 391        region2=force_region(region2)
 392    )
 393    assert(is_region(region1) && is_region(region2), "One of the inputs is not a region")
 394    len(region1) != len(region2)? false :
 395    __are_regions_equal(either_winding?_clockwise_region(region1):region1,
 396                        either_winding?_clockwise_region(region2):region2,
 397                        0);
 398
 399function __are_regions_equal(region1, region2, i) =
 400    i >= len(region1)? true :
 401    !_is_polygon_in_list(region1[i], region2)? false :
 402    __are_regions_equal(region1, region2, i+1);
 403
 404
 405/// Internal Function: _region_region_intersections()
 406/// Usage:
 407///    risect = _region_region_intersections(region1, region2, [closed1], [closed2], [eps]
 408/// Description:
 409///    Returns a pair of sorted lists such that risect[0] is a list of intersection
 410///    points for every path in region1, and similarly risect[1] is a list of intersection
 411///    points for the paths in region2.  For each path the intersection list is
 412///    a sorted list of the form [PATHIND, SEGMENT, U].  You can specify that the paths in either
 413///    region be regarded as open paths if desired.  Default is to treat them as
 414///    regions and hence the paths as closed polygons.
 415///    .
 416///    Included as intersection points are points where region1 touches itself at a vertex or
 417///    region2 touches itself at a vertex.  (The paths are assumed to have no self crossings.
 418///    Self crossings of the paths in the regions are not returned.)
 419function _region_region_intersections(region1, region2, closed1=true,closed2=true, eps=EPSILON) =
 420   let(
 421       intersections =   [
 422           for(p1=idx(region1))
 423              let(
 424                  path = closed1?close_path(region1[p1]):region1[p1]
 425              )
 426              for(i = [0:1:len(path)-2])
 427                  let(
 428                      a1 = path[i],
 429                      a2 = path[i+1],
 430                      nrm = norm(a1-a2)
 431                  )
 432                  if( nrm>eps )  // ignore zero-length path edges
 433                       let( 
 434                           seg_normal = [-(a2-a1).y, (a2-a1).x]/nrm,
 435                           ref = a1*seg_normal
 436                       )
 437                           // `signs[j]` is the sign of the signed distance from
 438                           // poly vertex j to the line [a1,a2] where near zero
 439                           // distances are snapped to zero;  poly edges 
 440                           //  with equal signs at its vertices cannot intersect
 441                           // the path edge [a1,a2] or they are collinear and 
 442                           // further tests can be discarded.
 443                       for(p2=idx(region2))
 444                           let(
 445                               poly  = closed2?close_path(region2[p2]):region2[p2],
 446                               signs = [for(v=poly*seg_normal) abs(v-ref) < eps ? 0 : sign(v-ref) ]
 447                           ) 
 448                           if(max(signs)>=0 && min(signs)<=0) // some edge intersects line [a1,a2]
 449                               for(j=[0:1:len(poly)-2]) 
 450                                   if(signs[j]!=signs[j+1])
 451                                        let( // exclude non-crossing and collinear segments
 452                                            b1 = poly[j],
 453                                            b2 = poly[j+1],
 454                                            isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
 455                                        )
 456                                        if (isect 
 457                                            && isect[1]>= -eps 
 458                                            && isect[1]<= 1+eps 
 459                                            && isect[2]>= -eps
 460                                            && isect[2]<= 1+eps)       
 461                                         [[p1,i,isect[1]], [p2,j,isect[2]]]
 462         ],
 463         regions=[region1,region2],
 464         // Create a flattened index list corresponding to the points in region1 and region2
 465         // that gives each point as an intersection point
 466         ptind = [for(i=[0:1])   
 467                    [for(p=idx(regions[i]))
 468                       for(j=idx(regions[i][p])) [p,j,0]]],
 469         points = [for(i=[0:1]) flatten(regions[i])],
 470         // Corner points are those points where the region touches itself, hence duplicate
 471         // points in the region's point set
 472         cornerpts = [for(i=[0:1])
 473                         [for(k=vector_search(points[i],eps,points[i]))
 474                             each if (len(k)>1) select(ptind[i],k)]],
 475         risect = [for(i=[0:1]) concat(column(intersections,i), cornerpts[i])],
 476         counts = [count(len(region1)), count(len(region2))],
 477         pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
 478       )
 479       [for(i=[0:1]) [for(j=counts[i]) _sort_vectors(select(risect[i],pathind[i][j]))]];
 480         
 481
 482// Section: Breaking up regions into subregions
 483
 484
 485// Function: split_region_at_region_crossings()
 486// Usage:
 487//   split_region = split_region_at_region_crossings(region1, region2, [closed1], [closed2], [eps])
 488// Description:
 489//   Splits region1 at the places where polygons in region1 touches each other at corners and at locations
 490//   where region1 intersections region2.  Split region2 similarly with respect to region1.
 491//   The return is a pair of results of the form [split1, split2] where split1=[frags1,frags2,...]
 492//   and frags1 is a list of paths that when placed end to end (in the given order), give the first polygon of region1.
 493//   Each path in the list is either entirely inside or entirely outside region2.  
 494//   Then frags2 is the decomposition of the second polygon into path pieces, and so on.  Finally split2 is
 495//   the same list, but for the polygons in region2.  
 496//   You can pass a single polygon in for either region, but the output will be a singleton list, as if
 497//   you passed in a singleton region.  If you set the closed parameters to false then the region components
 498//   will be treated as open paths instead of polygons.  
 499// Arguments:
 500//   region1 = first region
 501//   region2 = second region
 502//   closed1 = if false then treat region1 as list of open paths.  Default: true
 503//   closed2 = if false then treat region2 as list of open paths.  Default: true
 504//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 505// Example(2D): 
 506//   path = square(50,center=false);
 507//   region = [circle(d=80), circle(d=40)];
 508//   paths = split_region_at_region_crossings(path, region);
 509//   color("#aaa") region(region);
 510//   rainbow(paths[0][0]) stroke($item, width=2);
 511//   right(110){
 512//     color("#aaa") region([path]);
 513//     rainbow(flatten(paths[1])) stroke($item, width=2);
 514//   }
 515function split_region_at_region_crossings(region1, region2, closed1=true, closed2=true, eps=EPSILON) = 
 516    let(
 517        region1=force_region(region1),
 518        region2=force_region(region2)
 519    )
 520    assert(is_region(region1) && is_region(region2),"One of the inputs is not a region")
 521    let(
 522        xings = _region_region_intersections(region1, region2, closed1, closed2, eps),
 523        regions = [region1,region2],
 524        closed = [closed1,closed2]
 525    )
 526    [for(i=[0:1])
 527      [for(p=idx(xings[i]))
 528        let(
 529            crossings = deduplicate([
 530                                     [p,0,0],
 531                                     each xings[i][p],
 532                                     [p,len(regions[i][p])-(closed[i]?1:2), 1],
 533                                    ],eps=eps),
 534            subpaths = [
 535                for (frag = pair(crossings)) 
 536                    deduplicate(
 537                        _path_select(regions[i][p], frag[0][1], frag[0][2], frag[1][1], frag[1][2], closed=closed[i]),
 538                        eps=eps
 539                    )
 540            ]
 541        )
 542        [for(s=subpaths) if (len(s)>1) s]
 543       ]
 544    ];
 545                
 546                
 547
 548// Function: region_parts()
 549// Usage:
 550//   rgns = region_parts(region);
 551// Description:
 552//   Divides a region into a list of connected regions.  Each connected region has exactly one clockwise outside boundary
 553//   and zero or more counter-clockwise outlines defining internal holes.  Note that behavior is undefined on invalid regions whose
 554//   components cross each other.
 555// Example(2D,NoAxes):
 556//   R = [for(i=[1:7]) square(i,center=true)];
 557//   region_list = region_parts(R);
 558//   rainbow(region_list) region($item);
 559// Example(2D,NoAxes):
 560//   R = [back(7,square(3,center=true)),
 561//        square([20,10],center=true),
 562//        left(5,square(8,center=true)),
 563//        for(i=[4:2:8])
 564//          right(5,square(i,center=true))];
 565//   region_list = region_parts(R);
 566//   rainbow(region_list) region($item);
 567function region_parts(region) =
 568   let(
 569       region = force_region(region)
 570   )
 571   assert(is_region(region), "Input is not a region")
 572   let(
 573       inside = [for(i=idx(region))
 574                    let(pt = mean([region[i][0], region[i][1]]))
 575                    [for(j=idx(region))  i==j ? 0
 576                                       : point_in_polygon(pt,region[j]) >=0 ? 1 : 0]
 577                ],
 578       level = inside*repeat(1,len(region))
 579   )
 580   [ for(i=idx(region))
 581      if(level[i]%2==0)
 582         let(
 583             possible_children = search([level[i]+1],level,0)[0],
 584             keep=search([1], select(inside,possible_children), 0, i)[0]
 585         )
 586         [
 587           clockwise_polygon(region[i]),
 588           for(good=keep)
 589              ccw_polygon(region[possible_children[good]])
 590         ]
 591    ];
 592
 593
 594
 595
 596// Section: Offset and 2D Boolean Set Operations
 597
 598
 599function _offset_chamfer(center, points, delta) =
 600    let(
 601        dist = sign(delta)*norm(center-line_intersection(select(points,[0,2]), [center, points[1]])),
 602        endline = _shift_segment(select(points,[0,2]), delta-dist)
 603    ) [
 604        line_intersection(endline, select(points,[0,1])),
 605        line_intersection(endline, select(points,[1,2]))
 606    ];
 607
 608
 609function _shift_segment(segment, d) =
 610    assert(!approx(segment[0],segment[1]),"Path has repeated points")
 611    move(d*line_normal(segment),segment);
 612
 613
 614// Extend to segments to their intersection point.  First check if the segments already have a point in common,
 615// which can happen if two colinear segments are input to the path variant of `offset()`
 616function _segment_extension(s1,s2) =
 617    norm(s1[1]-s2[0])<1e-6 ? s1[1] : line_intersection(s1,s2,LINE,LINE);
 618
 619
 620function _makefaces(direction, startind, good, pointcount, closed) =
 621    let(
 622        lenlist = list_bset(good, pointcount),
 623        numfirst = len(lenlist),
 624        numsecond = sum(lenlist),
 625        prelim_faces = _makefaces_recurse(startind, startind+len(lenlist), numfirst, numsecond, lenlist, closed)
 626    )
 627    direction? [for(entry=prelim_faces) reverse(entry)] : prelim_faces;
 628
 629
 630function _makefaces_recurse(startind1, startind2, numfirst, numsecond, lenlist, closed, firstind=0, secondind=0, faces=[]) =
 631    // We are done if *both* firstind and secondind reach their max value, which is the last point if !closed or one past
 632    // the last point if closed (wrapping around).  If you don't check both you can leave a triangular gap in the output.
 633    ((firstind == numfirst - (closed?0:1)) && (secondind == numsecond - (closed?0:1)))? faces :
 634    _makefaces_recurse(
 635        startind1, startind2, numfirst, numsecond, lenlist, closed, firstind+1, secondind+lenlist[firstind],
 636        lenlist[firstind]==0? (
 637            // point in original path has been deleted in offset path, so it has no match.  We therefore
 638            // make a triangular face using the current point from the offset (second) path
 639            // (The current point in the second path can be equal to numsecond if firstind is the last point)
 640            concat(faces,[[secondind%numsecond+startind2, firstind+startind1, (firstind+1)%numfirst+startind1]])
 641            // in this case a point or points exist in the offset path corresponding to the original path
 642        ) : (
 643            concat(faces,
 644                // First generate triangular faces for all of the extra points (if there are any---loop may be empty)
 645                [for(i=[0:1:lenlist[firstind]-2]) [firstind+startind1, secondind+i+1+startind2, secondind+i+startind2]],
 646                // Finish (unconditionally) with a quadrilateral face
 647                [
 648                    [
 649                        firstind+startind1,
 650                        (firstind+1)%numfirst+startind1,
 651                        (secondind+lenlist[firstind])%numsecond+startind2,
 652                        (secondind+lenlist[firstind]-1)%numsecond+startind2
 653                    ]
 654                ]
 655            )
 656        )
 657    );
 658
 659
 660// Determine which of the shifted segments are good
 661function _good_segments(path, d, shiftsegs, closed, quality) =
 662    let(
 663        maxind = len(path)-(closed ? 1 : 2),
 664        pathseg = [for(i=[0:maxind]) select(path,i+1)-path[i]],
 665        pathseg_len =  [for(seg=pathseg) norm(seg)],
 666        pathseg_unit = [for(i=[0:maxind]) pathseg[i]/pathseg_len[i]],
 667        // Order matters because as soon as a valid point is found, the test stops
 668        // This order works better for circular paths because they succeed in the center
 669        alpha = concat([for(i=[1:1:quality]) i/(quality+1)],[0,1])
 670    ) [
 671        for (i=[0:len(shiftsegs)-1])
 672            (i>maxind)? true :
 673            _segment_good(path,pathseg_unit,pathseg_len, d - 1e-7, shiftsegs[i], alpha)
 674    ];
 675
 676
 677// Determine if a segment is good (approximately)
 678// Input is the path, the path segments normalized to unit length, the length of each path segment
 679// the distance threshold, the segment to test, and the locations on the segment to test (normalized to [0,1])
 680// The last parameter, index, gives the current alpha index.
 681//
 682// A segment is good if any part of it is farther than distance d from the path.  The test is expensive, so
 683// we want to quit as soon as we find a point with distance > d, hence the recursive code structure.
 684//
 685// This test is approximate because it only samples the points listed in alpha.  Listing more points
 686// will make the test more accurate, but slower.
 687function _segment_good(path,pathseg_unit,pathseg_len, d, seg,alpha ,index=0) =
 688    index == len(alpha) ? false :
 689    _point_dist(path,pathseg_unit,pathseg_len, alpha[index]*seg[0]+(1-alpha[index])*seg[1]) > d ? true :
 690    _segment_good(path,pathseg_unit,pathseg_len,d,seg,alpha,index+1);
 691
 692
 693// Input is the path, the path segments normalized to unit length, the length of each path segment
 694// and a test point.  Computes the (minimum) distance from the path to the point, taking into
 695// account that the minimal distance may be anywhere along a path segment, not just at the ends.
 696function _point_dist(path,pathseg_unit,pathseg_len,pt) =
 697    min([
 698        for(i=[0:len(pathseg_unit)-1]) let(
 699            v = pt-path[i],
 700            projection = v*pathseg_unit[i],
 701            segdist = projection < 0? norm(pt-path[i]) :
 702                projection > pathseg_len[i]? norm(pt-select(path,i+1)) :
 703                norm(v-projection*pathseg_unit[i])
 704        ) segdist
 705    ]);
 706
 707
 708// Function: offset()
 709// Usage:
 710//   offsetpath = offset(path, [r=|delta=], [chamfer=], [closed=], [check_valid=], [quality=], [same_length=])
 711//   path_faces = offset(path, return_faces=true, [r=|delta=], [chamfer=], [closed=], [check_valid=], [quality=], [firstface_index=], [flip_faces=])
 712// Description:
 713//   Takes a 2D input path, polygon or region and returns a path offset by the specified amount.  As with the built-in
 714//   offset() module, you can use `r` to specify rounded offset and `delta` to specify offset with
 715//   corners.  If you used `delta` you can set `chamfer` to true to get chamfers.
 716//   For paths and polygons positive offsets make the polygons larger.  For paths, 
 717//   positive offsets shift the path to the left, relative to the direction of the path.  Note
 718//   that the path must not include any 180 degree turns, where the path reverses direction.  
 719//   .
 720//   When offsets shrink the path, segments cross and become invalid.  By default `offset()` checks
 721//   for this situation.  To test validity the code checks that segments have distance larger than (r
 722//   or delta) from the input path.  This check takes O(N^2) time and may mistakenly eliminate
 723//   segments you wanted included in various situations, so you can disable it if you wish by setting
 724//   check_valid=false.  Another situation is that the test is not sufficiently thorough and some
 725//   segments persist that should be eliminated.  In this case, increase `quality` to 2 or 3.  (This
 726//   increases the number of samples on the segment that are checked.)  Run time will increase.  In
 727//   some situations you may be able to decrease run time by setting quality to 0, which causes only
 728//   segment ends to be checked.
 729//   .
 730//   When invalid segments are eliminated, the path length decreases.  If you use chamfering or rounding, then
 731//   the chamfers and roundings can increase the length of the output path.  Hence points in the output may be 
 732//   difficult to associate with the input.  If you want to maintain alignment between the points you
 733//   can use the `same_length` option.  This option requires that you use `delta=` with `chamfer=false` to ensure
 734//   that no points are added.  When points collapse to a single point in the offset, the output includes
 735//   that point repeated to preserve the correct length.  
 736//   .
 737//   Another way to obtain alignment information is to use the return_faces option, which can
 738//   provide alignment information for all offset parameters: it returns a face list which lists faces between
 739//   the original path and the offset path where the vertices are ordered with the original path
 740//   first, starting at `firstface_index` and the offset path vertices appearing afterwords.  The
 741//   direction of the faces can be flipped using `flip_faces`.  When you request faces the return
 742//   value is a list: [offset_path, face_list].
 743// Arguments:
 744//   path = the path to process.  A list of 2d points.
 745//   ---
 746//   r = offset radius.  Distance to offset.  Will round over corners.
 747//   delta = offset distance.  Distance to offset with pointed corners.
 748//   chamfer = chamfer corners when you specify `delta`.  Default: false
 749//   closed = if true path is treate as a polygon. Default: False.
 750//   check_valid = perform segment validity check.  Default: True.
 751//   quality = validity check quality parameter, a small integer.  Default: 1.
 752//   same_length = return a path with the same length as the input.  Only compatible with `delta=`.  Default: false
 753//   return_faces = return face list.  Default: False.
 754//   firstface_index = starting index for face list.  Default: 0.
 755//   flip_faces = flip face direction.  Default: false
 756// Example(2D,NoAxes):
 757//   star = star(5, r=100, ir=30);
 758//   #stroke(closed=true, star, width=3);
 759//   stroke(closed=true, width=3, offset(star, delta=10, closed=true));
 760// Example(2D,NoAxes):
 761//   star = star(5, r=100, ir=30);
 762//   #stroke(closed=true, star, width=3);
 763//   stroke(closed=true, width=3,
 764//          offset(star, delta=10, chamfer=true, closed=true));
 765// Example(2D,NoAxes):
 766//   star = star(5, r=100, ir=30);
 767//   #stroke(closed=true, star, width=3);
 768//   stroke(closed=true, width=3,
 769//          offset(star, r=10, closed=true));
 770// Example(2D,NoAxes):
 771//   star = star(7, r=120, ir=50);
 772//   #stroke(closed=true, width=3, star);
 773//   stroke(closed=true, width=3,
 774//          offset(star, delta=-15, closed=true));
 775// Example(2D,NoAxes):
 776//   star = star(7, r=120, ir=50);
 777//   #stroke(closed=true, width=3, star);
 778//   stroke(closed=true, width=3,
 779//          offset(star, delta=-15, chamfer=true, closed=true));
 780// Example(2D,NoAxes):
 781//   star = star(7, r=120, ir=50);
 782//   #stroke(closed=true, width=3, star);
 783//   stroke(closed=true, width=3,
 784//          offset(star, r=-15, closed=true, $fn=20));
 785// Example(2D,NoAxes):  This case needs `quality=2` for success
 786//   test = [[0,0],[10,0],[10,7],[0,7], [-1,-3]];
 787//   polygon(offset(test,r=-1.9, closed=true, quality=2));
 788//   //polygon(offset(test,r=-1.9, closed=true, quality=1));  // Fails with erroneous 180 deg path error
 789//   %down(.1)polygon(test);
 790// Example(2D,NoAxes): This case fails if `check_valid=true` when delta is large enough because segments are too close to the opposite side of the curve.  
 791//   star = star(5, r=22, ir=13);
 792//   stroke(star,width=.3,closed=true);                                                           
 793//   color("green")
 794//     stroke(offset(star, delta=-9, closed=true),width=.3,closed=true); // Works with check_valid=true (the default)
 795//   color("red")
 796//     stroke(offset(star, delta=-10, closed=true, check_valid=false),   // Fails if check_valid=true 
 797//            width=.3,closed=true); 
 798// Example(2D): But if you use rounding with offset then you need `check_valid=true` when `r` is big enough.  It works without the validity check as long as the offset shape retains a some of the straight edges at the star tip, but once the shape shrinks smaller than that, it fails.  There is no simple way to get a correct result for the case with `r=10`, because as in the previous example, it will fail if you turn on validity checks.  
 799//   star = star(5, r=22, ir=13);
 800//   color("green")
 801//     stroke(offset(star, r=-8, closed=true,check_valid=false), width=.1, closed=true);
 802//   color("red")
 803//     stroke(offset(star, r=-10, closed=true,check_valid=false), width=.1, closed=true);
 804// Example(2D,NoAxes): The extra triangles in this example show that the validity check cannot be skipped
 805//   ellipse = scale([20,4], p=circle(r=1,$fn=64));
 806//   stroke(ellipse, closed=true, width=0.3);
 807//   stroke(offset(ellipse, r=-3, check_valid=false, closed=true),
 808//          width=0.3, closed=true);
 809// Example(2D,NoAxes): The triangles are removed by the validity check
 810//   ellipse = scale([20,4], p=circle(r=1,$fn=64));
 811//   stroke(ellipse, closed=true, width=0.3);
 812//   stroke(offset(ellipse, r=-3, check_valid=true, closed=true),
 813//          width=0.3, closed=true);
 814// Example(2D): Open path.  The path moves from left to right and the positive offset shifts to the left of the initial red path.
 815//   sinpath = 2*[for(theta=[-180:5:180]) [theta/4,45*sin(theta)]];
 816//   #stroke(sinpath, width=2);
 817//   stroke(offset(sinpath, r=17.5),width=2);
 818// Example(2D,NoAxes): Region
 819//   rgn = difference(circle(d=100),
 820//                    union(square([20,40], center=true),
 821//                          square([40,20], center=true)));
 822//   #linear_extrude(height=1.1) stroke(rgn, width=1);
 823//   region(offset(rgn, r=-5));
 824// Example(2D,NoAxes): Using `same_length=true` to align the original curve to the offset.  Note that lots of points map to the corner at the top.
 825//   closed=false;
 826//   path = [for(angle=[0:5:180]) 10*[angle/100,2*sin(angle)]];
 827//   opath = offset(path, delta=-3,same_length=true,closed=closed);
 828//   stroke(path,closed=closed,width=.3);
 829//   stroke(opath,closed=closed,width=.3);
 830//   color("red") for(i=idx(path)) stroke([path[i],opath[i]],width=.3);
 831
 832function offset(
 833    path, r=undef, delta=undef, chamfer=false,
 834    closed=false, check_valid=true,
 835    quality=1, return_faces=false, firstface_index=0,
 836    flip_faces=false, same_length=false
 837) =
 838    assert(!(same_length && return_faces), "Cannot combine return_faces with same_length")
 839    is_region(path)?
 840        assert(!return_faces, "return_faces not supported for regions.")
 841        let(
 842            ofsregs = [for(R=region_parts(path))
 843                difference([for(i=idx(R)) offset(R[i], r=u_mul(i>0?-1:1,r), delta=u_mul(i>0?-1:1,delta),
 844                                      chamfer=chamfer, check_valid=check_valid, quality=quality,closed=true)])]
 845        )
 846        union(ofsregs)
 847    :
 848    let(rcount = num_defined([r,delta]))
 849    assert(rcount==1,"Must define exactly one of 'delta' and 'r'")
 850    assert(!same_length || (is_def(delta) && !chamfer), "Must specify delta, with chamfer=false, when same_length=true")
 851    assert(is_path(path), "Input must be a path or region")
 852    let(
 853        chamfer = is_def(r) ? false : chamfer,
 854        quality = max(0,round(quality)),
 855        flip_dir = closed && !is_polygon_clockwise(path)? -1 : 1,
 856        d = flip_dir * (is_def(r) ? r : delta),
 857//        shiftsegs = [for(i=[0:len(path)-1]) _shift_segment(select(path,i,i+1), d)],
 858        shiftsegs = [for(i=[0:len(path)-2]) _shift_segment([path[i],path[i+1]], d),
 859                     if (closed) _shift_segment([last(path),path[0]],d)
 860                     else [path[0],path[1]]  // dummy segment, not used
 861                    ],
 862        // good segments are ones where no point on the segment is less than distance d from any point on the path
 863        good = check_valid ? _good_segments(path, abs(d), shiftsegs, closed, quality) : repeat(true,len(shiftsegs)),
 864        goodsegs = bselect(shiftsegs, good),
 865        goodpath = bselect(path,good)
 866    )
 867    assert(len(goodsegs)-(!closed && select(good,-1)?1:0)>0,"Offset of path is degenerate")
 868    let(
 869        // Extend the shifted segments to their intersection points
 870        sharpcorners = [for(i=[0:len(goodsegs)-1]) _segment_extension(select(goodsegs,i-1), select(goodsegs,i))],
 871        // If some segments are parallel then the extended segments are undefined.  This case is not handled
 872        // Note if !closed the last corner doesn't matter, so exclude it
 873        parallelcheck =
 874            (len(sharpcorners)==2 && !closed) ||
 875            all_defined(closed? sharpcorners : select(sharpcorners, 1,-2))
 876    )
 877    assert(parallelcheck, "Path contains a segment that reverses direction (180 deg turn)")
 878    let(
 879        // This is a boolean array that indicates whether a corner is an outside or inside corner
 880        // For outside corners, the newcorner is an extension (angle 0), for inside corners, it turns backward
 881        // If either side turns back it is an inside corner---must check both.
 882        // Outside corners can get rounded (if r is specified and there is space to round them)
 883        outsidecorner = len(sharpcorners)==2 ? [false,false]
 884           :
 885            [for(i=[0:len(goodsegs)-1])
 886                let(prevseg=select(goodsegs,i-1))
 887                (i==0 || i==len(goodsegs)-1) && !closed ? false  // In open case first entry is bogus
 888               :  
 889                (goodsegs[i][1]-goodsegs[i][0]) * (goodsegs[i][0]-sharpcorners[i]) > 0
 890                 && (prevseg[1]-prevseg[0]) * (sharpcorners[i]-prevseg[1]) > 0
 891            ],
 892        steps = is_def(delta) ? [] : [
 893            for(i=[0:len(goodsegs)-1])  
 894                r==0 ? 0
 895                // if path is open but first and last entries match value is not used, but
 896                // computation below gives error, so special case handle it
 897              : i==len(goodsegs)-1 && !closed && approx(goodpath[i],goodsegs[i][0]) ? 0 
 898                // floor is important here to ensure we don't generate extra segments when nearly straight paths expand outward
 899              : 1+floor(segs(r)*vector_angle(   
 900                                             select(goodsegs,i-1)[1]-goodpath[i],
 901                                             goodsegs[i][0]-goodpath[i])
 902                        /360)
 903        ],
 904        // If rounding is true then newcorners replaces sharpcorners with rounded arcs where needed
 905        // Otherwise it's the same as sharpcorners
 906        // If rounding is on then newcorners[i] will be the point list that replaces goodpath[i] and newcorners later
 907        // gets flattened.  If rounding is off then we set it to [sharpcorners] so we can later flatten it and get
 908        // plain sharpcorners back.
 909        newcorners = is_def(delta) && !chamfer ? [sharpcorners]
 910            : [for(i=[0:len(goodsegs)-1])
 911                  (!chamfer && steps[i] <=1)  // Don't round if steps is smaller than 2
 912                  || !outsidecorner[i]        // Don't round inside corners
 913                  || (!closed && (i==0 || i==len(goodsegs)-1))  // Don't round ends of an open path
 914                ? [sharpcorners[i]]
 915                : chamfer ? _offset_chamfer(
 916                                  goodpath[i], [
 917                                      select(goodsegs,i-1)[1],
 918                                      sharpcorners[i],
 919                                      goodsegs[i][0]
 920                                  ], d
 921                              )     
 922                : // rounded case
 923                  arc(cp=goodpath[i],
 924                      points=[
 925                          select(goodsegs,i-1)[1],
 926                          goodsegs[i][0]
 927                      ],
 928                      n=steps[i])
 929              ],
 930        pointcount = (is_def(delta) && !chamfer)?
 931            repeat(1,len(sharpcorners)) :
 932            [for(i=[0:len(goodsegs)-1]) len(newcorners[i])],
 933        start = [goodsegs[0][0]],
 934        end = [goodsegs[len(goodsegs)-2][1]],
 935        edges =  closed?
 936            flatten(newcorners) :
 937            concat(start,slice(flatten(newcorners),1,-2),end),
 938        faces = !return_faces? [] :
 939            _makefaces(
 940                flip_faces, firstface_index, good,
 941                pointcount, closed
 942            ),
 943        final_edges = same_length ? select(edges, [0,each list_head (cumsum([for(g=good) g?1:0]))])
 944                                  : edges
 945    ) return_faces? [edges,faces] : final_edges;
 946
 947
 948
 949/// Internal Function: _filter_region_parts()
 950///
 951/// splits region1 into subpaths where either it touches itself or crosses region2.  Classifies all of the
 952/// subpaths as described below and keeps the ones listed in keep1.  A similar process is performed for region2.
 953/// All of the kept subpaths are assembled into polygons and returned as a lst.
 954/// .
 955/// The four types of subpath from the region are defined relative to the second region:
 956///    "O" - the subpath is outside the second region
 957///    "I" - the subpath is in the second region's interior
 958///    "S" - the subpath is on the 2nd region's border and the two regions interiors are on the same side of the subpath
 959///    "U" - the subpath is on the 2nd region's border and the two regions meet at the subpath from opposite sides
 960/// You specify which type of subpaths to keep with a string of the desired types such as "OS".  
 961function _filter_region_parts(region1, region2, keep, eps=EPSILON) = 
 962    // We have to compute common vertices between paths in the region because
 963    // they can be places where the path must be cut, even though they aren't
 964    // found my the split_path function.  
 965    let(
 966        subpaths = split_region_at_region_crossings(region1,region2,eps=eps),
 967        regions=[force_region(region1),
 968                 force_region(region2)]
 969    )        
 970    _assemble_path_fragments(
 971        [for(i=[0:1])
 972           let(
 973               keepS = search("S",keep[i])!=[],
 974               keepU = search("U",keep[i])!=[],        
 975               keepoutside = search("O",keep[i]) !=[],
 976               keepinside = search("I",keep[i]) !=[],
 977               all_subpaths = flatten(subpaths[i])
 978           )
 979           for (subpath = all_subpaths)
 980               let(
 981                   midpt = mean([subpath[0], subpath[1]]),
 982                   rel = point_in_region(midpt,regions[1-i],eps=eps),
 983                   keepthis = rel<0 ? keepoutside
 984                            : rel>0 ? keepinside
 985                            : !(keepS || keepU) ? false
 986                            : let(
 987                                  sidept = midpt + 0.01*line_normal(subpath[0],subpath[1]),
 988                                  rel1 = point_in_region(sidept,regions[0],eps=eps)>0,
 989                                  rel2 = point_in_region(sidept,regions[1],eps=eps)>0
 990                              )
 991                              rel1==rel2 ? keepS : keepU
 992               )
 993               if (keepthis) subpath
 994        ]
 995    );
 996
 997
 998function _list_three(a,b,c) =
 999   is_undef(b) ? a : 
1000   [
1001     a,
1002     if (is_def(b)) b,
1003     if (is_def(c)) c
1004   ];
1005
1006
1007
1008// Function&Module: union()
1009// Usage:
1010//   union() CHILDREN;
1011//   region = union(regions);
1012//   region = union(REGION1,REGION2);
1013//   region = union(REGION1,REGION2,REGION3);
1014// Description:
1015//   When called as a function and given a list of regions or 2D polygons,
1016//   returns the union of all given regions and polygons.  Result is a single region.
1017//   When called as the built-in module, makes the union of the given children.
1018// Arguments:
1019//   regions = List of regions to union.
1020// Example(2D):
1021//   shape1 = move([-8,-8,0], p=circle(d=50));
1022//   shape2 = move([ 8, 8,0], p=circle(d=50));
1023//   color("green") region(union(shape1,shape2));
1024//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1025function union(regions=[],b=undef,c=undef,eps=EPSILON) =
1026    let(regions=_list_three(regions,b,c))
1027    len(regions)==0? [] :
1028    len(regions)==1? regions[0] :
1029    let(regions=[for (r=regions) is_path(r)? [r] : r])
1030    union([
1031           _filter_region_parts(regions[0],regions[1],["OS", "O"], eps=eps),           
1032           for (i=[2:1:len(regions)-1]) regions[i]
1033          ],
1034          eps=eps
1035    );
1036
1037
1038// Function&Module: difference()
1039// Usage:
1040//   difference() CHILDREN;
1041//   region = difference(regions);
1042//   region = difference(REGION1,REGION2);
1043//   region = difference(REGION1,REGION2,REGION3);
1044// Description:
1045//   When called as a function, and given a list of regions or 2D polygons, 
1046//   takes the first region or polygon and differences away all other regions/polygons from it.  The resulting
1047//   region is returned.
1048//   When called as the built-in module, makes the set difference of the given children.
1049// Arguments:
1050//   regions = List of regions or polygons to difference.
1051// Example(2D):
1052//   shape1 = move([-8,-8,0], p=circle(d=50));
1053//   shape2 = move([ 8, 8,0], p=circle(d=50));
1054//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1055//   color("green") region(difference(shape1,shape2));
1056function difference(regions=[],b=undef,c=undef,eps=EPSILON) =
1057     let(regions = _list_three(regions,b,c))
1058     len(regions)==0? []
1059   : len(regions)==1? regions[0]
1060   : regions[0]==[] ? []
1061   : let(regions=[for (r=regions) is_path(r)? [r] : r])
1062     difference([
1063                 _filter_region_parts(regions[0],regions[1],["OU", "I"], eps=eps),                
1064                 for (i=[2:1:len(regions)-1]) regions[i]
1065                ],
1066                eps=eps
1067     );
1068
1069
1070// Function&Module: intersection()
1071// Usage:
1072//   intersection() CHILDREN;
1073//   region = intersection(regions);
1074//   region = intersection(REGION1,REGION2);
1075//   region = intersection(REGION1,REGION2,REGION3);
1076// Description:
1077//   When called as a function, and given a list of regions or polygons returns the
1078//   intersection of all given regions.  Result is a single region.
1079//   When called as the built-in module, makes the intersection of all the given children.
1080// Arguments:
1081//   regions = List of regions to intersect.
1082// Example(2D):
1083//   shape1 = move([-8,-8,0], p=circle(d=50));
1084//   shape2 = move([ 8, 8,0], p=circle(d=50));
1085//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1086//   color("green") region(intersection(shape1,shape2));
1087function intersection(regions=[],b=undef,c=undef,eps=EPSILON) =
1088     let(regions = _list_three(regions,b,c))
1089     len(regions)==0 ? []
1090   : len(regions)==1? regions[0]
1091   : regions[0]==[] || regions[1]==[] ? []   
1092   : intersection([
1093                   _filter_region_parts(regions[0],regions[1],["IS","I"],eps=eps),                       
1094                   for (i=[2:1:len(regions)-1]) regions[i]
1095                  ],
1096                  eps=eps
1097     );
1098
1099
1100
1101// Function&Module: exclusive_or()
1102// Usage:
1103//   exclusive_or() CHILDREN;
1104//   region = exclusive_or(regions);
1105//   region = exclusive_or(REGION1,REGION2);
1106//   region = exclusive_or(REGION1,REGION2,REGION3);
1107// Description:
1108//   When called as a function and given a list of regions or 2D polygons, 
1109//   returns the exclusive_or of all given regions.  Result is a single region.
1110//   When called as a module, performs a boolean exclusive-or of up to 10 children.  Note that when
1111//   the input regions cross each other the exclusive-or operator will produce shapes that
1112//   meet at corners (non-simple regions), which do not render in CGAL.  
1113// Arguments:
1114//   regions = List of regions or polygons to exclusive_or
1115// Example(2D): As Function.  A linear_sweep of this shape fails to render in CGAL.  
1116//   shape1 = move([-8,-8,0], p=circle(d=50));
1117//   shape2 = move([ 8, 8,0], p=circle(d=50));
1118//   for (shape = [shape1,shape2])
1119//       color("red") stroke(shape, width=0.5, closed=true);
1120//   color("green") region(exclusive_or(shape1,shape2));
1121// Example(2D): As Module.  A linear_extrude() of the resulting geometry fails to render in CGAL.  
1122//   exclusive_or() {
1123//       square(40,center=false);
1124//       circle(d=40);
1125//   }
1126function exclusive_or(regions=[],b=undef,c=undef,eps=EPSILON) =
1127     let(regions = _list_three(regions,b,c))
1128     len(regions)==0? []
1129   : len(regions)==1? force_region(regions[0])
1130   : regions[0]==[] ? exclusive_or(list_tail(regions))
1131   : regions[1]==[] ? exclusive_or(list_remove(regions,1))
1132   : exclusive_or([
1133                   _filter_region_parts(regions[0],regions[1],["IO","IO"],eps=eps),                  
1134                   for (i=[2:1:len(regions)-1]) regions[i]
1135                  ],
1136                  eps=eps
1137     );
1138
1139
1140module exclusive_or() {
1141    if ($children==1) {
1142        children();
1143    } else if ($children==2) {
1144        difference() {
1145            children(0);
1146            children(1);
1147        }
1148        difference() {
1149            children(1);
1150            children(0);
1151        }
1152    } else if ($children==3) {
1153        exclusive_or() {
1154            exclusive_or() {
1155                children(0);
1156                children(1);
1157            }
1158            children(2);
1159        }
1160    } else if ($children==4) {
1161        exclusive_or() {
1162            exclusive_or() {
1163                children(0);
1164                children(1);
1165            }
1166            exclusive_or() {
1167                children(2);
1168                children(3);
1169            }
1170        }
1171    } else if ($children==5) {
1172        exclusive_or() {
1173            exclusive_or() {
1174                children(0);
1175                children(1);
1176                children(2);
1177                children(3);
1178            }
1179            children(4);
1180        }
1181    } else if ($children==6) {
1182        exclusive_or() {
1183            exclusive_or() {
1184                children(0);
1185                children(1);
1186                children(2);
1187                children(3);
1188            }
1189            children(4);
1190            children(5);
1191        }
1192    } else if ($children==7) {
1193        exclusive_or() {
1194            exclusive_or() {
1195                children(0);
1196                children(1);
1197                children(2);
1198                children(3);
1199            }
1200            children(4);
1201            children(5);
1202            children(6);
1203        }
1204    } else if ($children==8) {
1205        exclusive_or() {
1206            exclusive_or() {
1207                children(0);
1208                children(1);
1209                children(2);
1210                children(3);
1211            }
1212            exclusive_or() {
1213                children(4);
1214                children(5);
1215                children(6);
1216                children(7);
1217            }
1218        }
1219    } else if ($children==9) {
1220        exclusive_or() {
1221            exclusive_or() {
1222                children(0);
1223                children(1);
1224                children(2);
1225                children(3);
1226            }
1227            exclusive_or() {
1228                children(4);
1229                children(5);
1230                children(6);
1231                children(7);
1232            }
1233            children(8);
1234        }
1235    } else if ($children==10) {
1236        exclusive_or() {
1237            exclusive_or() {
1238                children(0);
1239                children(1);
1240                children(2);
1241                children(3);
1242            }
1243            exclusive_or() {
1244                children(4);
1245                children(5);
1246                children(6);
1247                children(7);
1248            }
1249            children(8);
1250            children(9);
1251        }
1252    } else {
1253        assert($children<=10, "exclusive_or() can only handle up to 10 children.");
1254    }
1255}
1256
1257
1258// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap