diff --git a/drawing.scad b/drawing.scad
index 52a9789..60ebcbe 100644
--- a/drawing.scad
+++ b/drawing.scad
@@ -167,6 +167,7 @@ module stroke(
     trim, trim1, trim2,
     convexity=10, hull=true
 ) {
+    no_children($children);
     module setcolor(clr) {
         if (clr==undef) {
             children();
diff --git a/mutators.scad b/mutators.scad
index ebc8f68..c69fcb7 100644
--- a/mutators.scad
+++ b/mutators.scad
@@ -1037,6 +1037,8 @@ module HSV(h,s=1,v=1,a=1) color(HSV(h,s,v),a) children();
 //   list = The list of items to iterate through.
 //   stride = Consecutive colors stride around the color wheel divided into this many parts.
 //   maxhues = max number of hues to use (to prevent lots of indistinguishable hues)
+//   shuffle = if true then shuffle the hues in a random order.  Default: false
+//   seed = seed to use for shuffle
 // Side Effects:
 //   Sets the color to progressive values along the ROYGBIV spectrum for each item.
 //   Sets `$idx` to the index of the current item in `list` that we want to show.
@@ -1046,14 +1048,13 @@ module HSV(h,s=1,v=1,a=1) color(HSV(h,s,v),a) children();
 // Example(2D):
 //   rgn = [circle(d=45,$fn=3), circle(d=75,$fn=4), circle(d=50)];
 //   rainbow(rgn) stroke($item, closed=true);
-module rainbow(list, stride=1, maxhues)
+module rainbow(list, stride=1, maxhues, shuffle=false, seed)
 {
     ll = len(list);
     maxhues = first_defined([maxhues,ll]);
     huestep = 360 / maxhues;
-    hues = [for (i=[0:1:ll-1]) posmod(i*huestep+i*360/stride,360)];
-    echo(hues=hues);
-    s = [for (i=[0:1:ll-1]) [.5,.7,1][posmod(i,3)]];
+    huelist = [for (i=[0:1:ll-1]) posmod(i*huestep+i*360/stride,360)];
+    hues = shuffle ? shuffle(huelist, seed=seed) : huelist;
     for($idx=idx(list)) {
         $item = list[$idx];
         HSV(h=hues[$idx]) children();
diff --git a/paths.scad b/paths.scad
index 833e7f2..40159bc 100644
--- a/paths.scad
+++ b/paths.scad
@@ -1170,7 +1170,7 @@ function _assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0,
 ///   fragments = List of paths to be assembled into complete polygons.
 ///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
 function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
-    len(fragments)==0? _finished :
+   len(fragments)==0? _finished :
     let(
         minxidx = min_index([
             for (frag=fragments) min(subindex(frag,0))
diff --git a/regions.scad b/regions.scad
index ed51be4..4255731 100644
--- a/regions.scad
+++ b/regions.scad
@@ -74,13 +74,22 @@ function check_and_fix_path(path, valid_dim=undef, closed=false, name="path") =
 
 // Function: sanitize_region()
 // Usage:
-//   r_fixed = sanitize_region(r);
+//   r_fixed = sanitize_region(r, [nonzero], [eps]);
 // Description:
 //   Takes a malformed input region that contains self-intersecting polygons or polygons
 //   that cross each other and converts it into a properly defined region without
 //   these defects.
-function sanitize_region(r) = exclusive_or([for(poly=r) each polygon_parts(poly)]);
-
+// Arguments:
+//   r = region to sanitize
+//   nonzero = set to true to use nonzero rule for polygon membership.  Default: false
+//   eps = Epsilon for geometric comparisons.  Default: `EPSILON` (1e-9)
+// Examples:
+//   
+function sanitize_region(r,nonzero=false,eps=EPSILON) =
+     assert(is_region(r))
+     exclusive_or(
+                  [for(poly=r) each polygon_parts(poly,nonzero,eps)],
+                  eps=eps);
 
 
 // Module: region()
@@ -100,12 +109,12 @@ function sanitize_region(r) = exclusive_or([for(poly=r) each polygon_parts(poly)
 //   region(rgn);
 module region(r)
 {
+    no_children($children);
+    r = is_path(r) ? [r] : r;
     points = flatten(r);
-    paths = [
-        for (i=[0:1:len(r)-1]) let(
-            start = default(sum([for (j=[0:1:i-1]) len(r[j])]),0)
-        ) [for (k=[0:1:len(r[i])-1]) start+k]
-    ];
+    lengths = [for(path=r) len(path)];
+    starts = [0,each cumsum(lengths)];
+    paths = [for(i=idx(r)) count(s=starts[i], n=lengths[i])];
     polygon(points=points, paths=paths);
 }
 
@@ -143,12 +152,13 @@ function point_in_region(point, region, eps=EPSILON, _i=0, _cnt=0) =
 //   region = region to check
 //   eps = tolerance for geometric omparisons.  Default: `EPSILON` = 1e-9
 function is_region_simple(region, eps=EPSILON) =
-   [for(p=region) if (!is_path_simple(p,closed=true,eps)) 1] == []
-   &&
+   [for(p=region) if (!is_path_simple(p,closed=true,eps)) 1] == [];
+
+/*   &&
    [for(i=[0:1:len(region)-2])
        if (_path_region_intersections(region[i], list_tail(region,i+1), eps=eps) != []) 1
    ] ==[];
-
+*/
 
 // Function: are_regions_equal()
 // Usage:
@@ -221,10 +231,11 @@ function old_path_region_intersections(path, region, closed=true, eps=EPSILON) =
 // find the intersection points of a path and the polygons of a
 // region; only crossing intersections are caught, no collinear
 // intersection is returned.
-function _path_region_intersections(path, region, closed=true, eps=EPSILON) =
+function _path_region_intersections(path, region, closed=true, eps=EPSILON, extra=[]) =
     let( path = closed ? close_path(path,eps=eps) : path )
     _sort_vectors(
-         [for(si = [0:1:len(path)-2]) let(
+         [ each extra,
+           for(si = [0:1:len(path)-2]) let(
             a1 = path[si],
             a2 = path[si+1],
             nrm = norm(a1-a2)
@@ -278,11 +289,11 @@ function _path_region_intersections(path, region, closed=true, eps=EPSILON) =
 //   paths = split_path_at_region_crossings(path, region);
 //   color("#aaa") region(region);
 //   rainbow(paths) stroke($item, closed=false, width=2);
-function split_path_at_region_crossings(path, region, closed=true, eps=EPSILON) =
+function split_path_at_region_crossings(path, region, closed=true, eps=EPSILON, extra=[]) =
     let(
         path = deduplicate(path, eps=eps),
         region = [for (path=region) deduplicate(path, eps=eps)],
-        xings = _path_region_intersections(path, region, closed=closed, eps=eps),
+        xings = _path_region_intersections(path, region, closed=closed, eps=eps, extra=extra),
         crossings = deduplicate(
             concat([[0,0]], xings, [[len(path)-1,1]]),
             eps=eps
@@ -317,35 +328,90 @@ function split_path_at_region_crossings(path, region, closed=true, eps=EPSILON)
 //          right(5,square(i,center=true))];
 //   region_list = region_parts(R);
 //   rainbow(region_list) region($item);
-function region_parts(region) =
-    let(
-        paths = sort(idx=0, [
-            for(i = idx(region)) let(
-                cnt = sum([
-                    for (j = idx(region)) if (i!=j)
-                    let(pt = lerp(region[i][0],region[i][1],0.5))
-                    point_in_polygon(pt, region[j]) >=0 ? 1 : 0
-                ])
-            ) [cnt, region[i]]
-        ]),
-        outs = [
-            for (candout = paths) let(
-                lev = candout[0],
-                parent = candout[1]
-            ) if (lev % 2 == 0) [
-                clockwise_polygon(parent),
-                for (path = paths) if (
-                    path[0] == lev+1 &&
-                    point_in_polygon(
-                        lerp(path[1][0], path[1][1], 0.5),
-                        parent
-                    ) >= 0
-                ) ccw_polygon(path[1])
-            ]
-        ]
-    ) outs;
 
 
+function old_region_parts(region) =
+    let(
+        paths = sort(idx=0,
+           [
+            for(i = idx(region))
+              let(
+                pt = mean([region[i][0],region[i][1]]),
+                cnt = sum([for (j = idx(region))
+                              if (i!=j && point_in_polygon(pt, region[j]) >=0) 1])
+              )
+              [cnt, region[i]]
+        ]),
+        
+        outs = [
+            for (candout = paths)
+               let(
+                  lev = candout[0],
+                  parent = candout[1]
+               )
+               if (lev % 2 == 0)
+                  [
+                    clockwise_polygon(parent),
+                    for (path = paths)
+                        if (
+                            path[0] == lev+1
+                            && point_in_polygon(
+                                                lerp(path[1][0], path[1][1], 0.5)
+                                                ,parent
+                                               ) >= 0
+                        )
+                        ccw_polygon(path[1])
+                  ]
+        ]
+    )
+    outs;
+
+
+function inside(region, ins=[]) =
+  let(
+      i = len(ins)
+  )
+  i==len(region) ? ins
+  :
+  let(
+      pt=mean([region[i][0],region[i][1]])
+  )
+  i==0 ? inside(region,
+                       [[0,
+                         for(j=[1:1:len(region)-1]) point_in_polygon(pt,region[j])>=0 ? 1 : 0]])
+  : let(
+        prev = [for(j=[0:i-1]) point_in_polygon(pt,region[j])>=0 ? 1 : 0],
+        check = sum(bselect(ins,prev),repeat(0,len(region))),
+        next = [for(j=[i+1:1:len(region)-1]) check[j]>0 ? 1 : point_in_polygon(pt,region[j])>=0 ? 1 : 0]
+    )
+    inside(region, [
+                    each ins,
+                    [each prev, 0, each next]
+                   ]);
+
+
+function region_parts(region) =
+   let(
+       inside = [for(i=idx(region))
+                    let(pt = mean([region[i][0], region[i][1]]))
+                    [for(j=idx(region))  i==j ? 0
+                                       : point_in_polygon(pt,region[j]) >=0 ? 1 : 0]
+                ],
+       level = inside*repeat(1,len(region))
+   )
+   [ for(i=idx(region))
+      if(level[i]%2==0)
+         let(
+             possible_children = search([level[i]+1],level,0)[0],
+             keep=search([1], select(inside,possible_children), 0, i)[0]
+         )
+         [
+           clockwise_polygon(region[i]),
+           for(good=keep)
+              ccw_polygon(region[possible_children[good]])
+         ]
+    ];
+
 
 // Section: Region Extrusion and VNFs
 
@@ -647,7 +713,18 @@ function _point_dist(path,pathseg_unit,pathseg_len,pt) =
     ]);
 
 
-function _offset_region(
+function _offset_region(region, r, delta, chamfer, check_valid, quality,closed,return_faces,firstface_index,flip_faces) =
+    let(
+//dg=        echo(region=region),
+         reglist = [for(R=region_parts(region)) is_path(R) ? [R] : R],
+//fdsa=         echo(reglist),
+         ofsregs = [for(R=reglist)
+             [for(i=idx(R)) offset(R[i], r=u_mul(i>0?-1:1,r), delta=u_mul(i>0?-1:1,delta), chamfer=chamfer, check_valid=check_valid,
+                                    quality=quality,closed=true)]]
+    )
+    union(ofsregs);
+
+function d_offset_region(
     paths, r, delta, chamfer, closed,
     check_valid, quality,
     return_faces, firstface_index,
@@ -785,7 +862,10 @@ function offset(
     quality=1, return_faces=false, firstface_index=0,
     flip_faces=false
 ) = 
-    is_region(path)? (
+    is_region(path)? _offset_region(path,r=r,delta=delta,chamfer=chamfer,quality=quality,check_valid=check_valid)
+  
+/*
+        (
         assert(!return_faces, "return_faces not supported for regions.")
         let(
             path = [for (p=path) clockwise_polygon(p)],
@@ -805,7 +885,8 @@ function offset(
             return_faces=return_faces, firstface_index=firstface_index,
             flip_faces=flip_faces
         )
-    ) : let(rcount = num_defined([r,delta]))
+    )*/
+    : let(rcount = num_defined([r,delta]))
     assert(rcount==1,"Must define exactly one of 'delta' and 'r'")
     let(
         chamfer = is_def(r) ? false : chamfer,
@@ -904,37 +985,40 @@ function offset(
 ///    "S" - the subpath is on the region's border and the polygon and region are on the same side of the subpath
 ///    "U" - the subpath is on the region's border and the polygon and region meet at the subpath (from opposite sides)
 /// The return has the form of a list with entries [TAG, SUBPATH]
-function _tag_subpaths(path, region, eps=EPSILON) =
+function _tag_subpaths(region1, region2, eps=EPSILON, SUtags=true) =
+    // We have to compute common vertices between paths in the region because
+    // they can be places where the path must be cut, even though they aren't
+    // found my the split_path function.  
     let(
-        subpaths = split_path_at_region_crossings(path, region, eps=eps),
-        tagged = [
-            for (subpath = subpaths)
+        points = flatten(region1),
+        tree =  len(points)>0 ? vector_search_tree(points): undef
+    )
+    [for(path=region1)
+        let(
+            self_int = is_undef(tree)?[]:[for(i=idx(path)) if (len(vector_search(path[i], eps, tree))>1) [i,0]],
+            subpaths = split_path_at_region_crossings(path, region2, eps=eps, extra=self_int)
+        )
+        for (subpath = subpaths)
+            let(
+                midpt = mean([subpath[0], subpath[1]]),
+                rel = point_in_region(midpt,region2,eps=eps)
+            )
+            if(rel<0) ["O", subpath]
+            else if (rel>0) ["I", subpath]
+            else if (SUtags)
                 let(
-                    midpt = mean([subpath[0], subpath[1]]),
-                    rel = point_in_region(midpt,region,eps=eps)
-                )
-                rel<0? ["O", subpath]
-              : rel>0? ["I", subpath]
-              : let(
-                    vec = unit(subpath[1]-subpath[0]),
-                    perp = rot(90, planar=true, p=vec),
-                    sidept = midpt + perp*0.01,
-                    rel1 = point_in_polygon(sidept,path,eps=eps)>0,
-                    rel2 = point_in_region(sidept,region,eps=eps)>0
+                    sidept = midpt + 0.01*line_normal(subpath[0],subpath[1]),
+                    rel1 = point_in_region(sidept,region1,eps=eps)>0,
+                    rel2 = point_in_region(sidept,region2,eps=eps)>0
                  )
                  rel1==rel2? ["S", subpath] : ["U", subpath]
-        ]
-    ) tagged;
-
-
-function _tag_region_subpaths(region1, region2, eps=EPSILON) =
-    [for (path=region1) each _tag_subpaths(path, region2, eps=eps)];
+    ];
 
 
 function _tagged_region(region1,region2,keep1,keep2,eps=EPSILON) =
     let(
-        tagged1 = _tag_region_subpaths(region1, region2, eps=eps),
-        tagged2 = _tag_region_subpaths(region2, region1, eps=eps),
+        tagged1 = _tag_subpaths(region1, region2, eps=eps),
+        tagged2 = _tag_subpaths(region2, region1, eps=eps),
         tagged = [
                   for (tagpath = tagged1) if (in_list(tagpath[0], keep1)) tagpath[1],
                   for (tagpath = tagged2) if (in_list(tagpath[0], keep2)) tagpath[1]
@@ -963,14 +1047,14 @@ function _tagged_region(region1,region2,keep1,keep2,eps=EPSILON) =
 //   color("green") region(union(shape1,shape2));
 function union(regions=[],b=undef,c=undef,eps=EPSILON) =
     b!=undef? union(concat([regions],[b],c==undef?[]:[c]), eps=eps) :
-    len(regions)<=1? regions[0] :
-    union(
-        let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
-        concat(
-            [_tagged_region(regions[0],regions[1],["O","S"],["O"], eps=eps)],
-            [for (i=[2:1:len(regions)-1]) regions[i]]
-        ),
-        eps=eps
+    len(regions)==0? [] :
+    len(regions)==1? regions[0] :
+    let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
+    union([
+           _tagged_region(regions[0],regions[1],["O","S"],["O"], eps=eps),
+            for (i=[2:1:len(regions)-1]) regions[i]
+          ],
+          eps=eps
     );
 
 
@@ -996,13 +1080,13 @@ function difference(regions=[],b=undef,c=undef,eps=EPSILON) =
     b!=undef? difference(concat([regions],[b],c==undef?[]:[c]), eps=eps) :
     len(regions)==0? [] : 
     len(regions)==1? regions[0] :
-    difference(
-        let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
-        concat(
-            [_tagged_region(regions[0],regions[1],["O","U"],["I"], eps=eps)],
-            [for (i=[2:1:len(regions)-1]) regions[i]]
-        ),
-        eps=eps
+    regions[0]==[] ? [] : 
+    let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
+    difference([
+                _tagged_region(regions[0],regions[1],["O","U"],["I"], eps=eps),
+                for (i=[2:1:len(regions)-1]) regions[i]
+               ],
+               eps=eps
     );
 
 
@@ -1024,9 +1108,11 @@ function difference(regions=[],b=undef,c=undef,eps=EPSILON) =
 //   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
 //   color("green") region(intersection(shape1,shape2));
 function intersection(regions=[],b=undef,c=undef,eps=EPSILON) =
+//    echo(regions=regions)
      b!=undef? intersection(concat([regions],[b],c==undef?[]:[c]),eps=eps)
    : len(regions)==0 ? []
    : len(regions)==1? regions[0]
+   : regions[0]==[] || regions[1]==[] ? []   
    : let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
          intersection([
                        _tagged_region(regions[0],regions[1],["I","S"],["I"],eps=eps),
@@ -1063,17 +1149,14 @@ function intersection(regions=[],b=undef,c=undef,eps=EPSILON) =
 //   }
 function exclusive_or(regions=[],b=undef,c=undef,eps=EPSILON) =
     b!=undef? exclusive_or([regions, b, if(is_def(c)) c],eps=eps) :
-    len(regions)<=1? regions[0] :
-    exclusive_or(
-        let(regions=[for (r=regions) is_path(r)? [r] : r])
-        concat(
-            [union([
-                difference([regions[0],regions[1]], eps=eps),
-                difference([regions[1],regions[0]], eps=eps)
-            ], eps=eps)],
-            [for (i=[2:1:len(regions)-1]) regions[i]]
-        ),
-        eps=eps
+    len(regions)==0? [] :
+    len(regions)==1? regions[0] :
+    let(regions=[for (r=regions) is_path(r)? [r] : r])
+    exclusive_or([
+                  _tagged_region(regions[0],regions[1],["I","O"],["I","O"],eps=eps),
+                  for (i=[2:1:len(regions)-1]) regions[i]
+                 ],
+                 eps=eps
     );
 
 
diff --git a/tests/test_regions.scad b/tests/test_regions.scad
index 3369cd3..61e2eeb 100644
--- a/tests/test_regions.scad
+++ b/tests/test_regions.scad
@@ -10,3 +10,33 @@ module test_is_region() {
     assert(!is_region("foo"));
 }
 test_is_region();
+
+
+
+module test_union() {
+  R1 = [square(10,center=true), square(9,center=true)];
+  R2 = [square(9,center=true)];
+  assert(are_regions_equal(union(R1,R2), [square(10,center=true)]));
+}
+test_union();
+
+
+module test_intersection() {
+  R1 = [square(10,center=true), square(9,center=true)];
+  R6 = [square(9.5,center=true), square(9,center=true)];
+  assert(are_regions_equal(intersection(R6,R1), R6));
+  assert(are_regions_equal(intersection(R1,R6), R6));
+}
+test_intersection();
+
+
+module test_difference() {
+  R5 = [square(10,center=true), square(9,center=true),square(4,center=true)];
+  R4 = [square(9,center=true), square(3,center=true)];
+  assert(are_regions_equal(difference(R5,R4),
+                         [square(10,center=true), square(9, center=true), square(3,center=true)]));
+}
+test_difference();
+
+
+