diff --git a/geometry.scad b/geometry.scad
index 692e9b5..44a113f 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -1533,6 +1533,13 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
     assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
             "The point and polygon should be in 2D. The polygon should have more that 2 points." )
     assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
+    // Check bounding box
+    let(
+        box = pointlist_bounds(poly)
+    )
+    point.x<box[0].x-eps || point.x>box[1].x+eps
+        || point.y<box[0].y-eps || point.y>box[1].y+eps  ? -1
+    :
     // Does the point lie on any edges?  If so return 0.
     let(
         on_brd = [
diff --git a/regions.scad b/regions.scad
index e0b6d61..4d031b0 100644
--- a/regions.scad
+++ b/regions.scad
@@ -35,6 +35,13 @@
 function is_region(x) = is_list(x) && is_path(x.x);
 
 
+// Function: force_region()
+// Usage:
+//   region = force_region(path)
+// Description:
+//   If the input is a path then return it as a region.  Otherwise return it unaltered.
+function force_region(path) = is_path(path) ? [path] : path;
+
 
 // Function: check_and_fix_path()
 // Usage:
@@ -160,6 +167,8 @@ function is_region_simple(region, eps=EPSILON) =
    ] ==[];
 */
 
+function _clockwise_region(r) = [for(p=r) clockwise_polygon(p)];
+
 // Function: are_regions_equal()
 // Usage:
 //    b = are_regions_equal(region1, region2, [eps])
@@ -170,10 +179,16 @@ function is_region_simple(region, eps=EPSILON) =
 //    region1 = first region
 //    region2 = second region
 //    eps = tolerance for comparison
-function are_regions_equal(region1, region2) =
+function are_regions_equal(region1, region2, either_winding=false) =
+    let(
+        region1=force_region(region1),
+        region2=force_region(region2)
+    )
     assert(is_region(region1) && is_region(region2))
     len(region1) != len(region2)? false :
-    __are_regions_equal(region1, region2, 0);
+    __are_regions_equal(either_winding?_clockwise_region(region1):region1,
+                        either_winding?_clockwise_region(region2):region2,
+                        0);
 
 function __are_regions_equal(region1, region2, i) =
     i >= len(region1)? true :
@@ -272,6 +287,102 @@ function _path_region_intersections(path, region, closed=true, eps=EPSILON, extr
     ]);
 
 
+
+// Returns a list [reg1,reg2] such that reg1[i] is a list of intersection points for path i
+// in region1 having the form [seg, u].  
+function _region_region_intersections(region1, region2, closed1=true,closed2=true, eps=EPSILON) =
+   let(
+       intersections =   [
+           for(p1=idx(region1))
+              let(
+                  path = closed1?close_path(region1[p1]):region1[p1]
+              )
+              for(i = [0:1:len(path)-2])
+                  let(
+                      a1 = path[i],
+                      a2 = path[i+1],
+                      nrm = norm(a1-a2)
+                  )
+                  if( nrm>eps )  // ignore zero-length path edges
+                       let( 
+                           seg_normal = [-(a2-a1).y, (a2-a1).x]/nrm,
+                           ref = a1*seg_normal
+                       )
+                           // `signs[j]` is the sign of the signed distance from
+                           // poly vertex j to the line [a1,a2] where near zero
+                           // distances are snapped to zero;  poly edges 
+                           //  with equal signs at its vertices cannot intersect
+                           // the path edge [a1,a2] or they are collinear and 
+                           // further tests can be discarded.
+                       for(p2=idx(region2))
+                           let(
+                               poly  = closed2?close_path(region2[p2]):region2[p2],
+                               signs = [for(v=poly*seg_normal) v-ref> eps ? 1 : v-ref<-eps ? -1 : 0]
+                           ) 
+                           if(max(signs)>=0 && min(signs)<=0) // some edge edge intersects line [a1,a2]
+                               for(j=[0:1:len(poly)-2]) 
+                                   if(signs[j]!=signs[j+1])
+                                        let( // exclude non-crossing and collinear segments
+                                            b1 = poly[j],
+                                            b2 = poly[j+1],
+                                            isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
+                                        )
+                                        if (isect 
+                                            && isect[1]>= -eps 
+                                            && isect[1]<= 1+eps 
+                                            && isect[2]>= -eps
+                                            && isect[2]<= 1+eps)       
+                                         [[p1,i,isect[1]], [p2,j,isect[2]]]
+         ],
+         regions=[region1,region2],
+         // Create a flattened index list corresponding to the points in region1 and region2
+         // that gives each point as an intersection point
+         ptind = [for(i=[0:1])   
+                    [for(p=idx(regions[i]))
+                       for(j=idx(regions[i][p])) [p,j,0]]],
+         points = [for(i=[0:1]) flatten(regions[i])],
+         // Corner points are those points where the region touches itself, hence duplicate
+         // points in the region's point set
+         cornerpts = [for(i=[0:1])
+                         [for(k=vector_search(points[i],eps,points[i]))
+                             each if (len(k)>1) select(ptind[i],k)]],
+         risect = [for(i=[0:1]) concat(subindex(intersections,i), cornerpts[i])],
+         counts = [count(len(region1)), count(len(region2))],
+         pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
+       )
+       [for(i=[0:1]) [for(j=counts[i]) _sort_vectors(select(risect[i],pathind[i][j]))]];
+         
+
+
+function split_region_at_region_crossings(region1, region2, closed1=true, closed2=true, eps=EPSILON) = 
+    let(
+        xings = _region_region_intersections(region1, region2, closed1, closed2, eps),
+        regions = [region1,region2],
+        closed = [closed1,closed2]
+    )
+    [for(i=[0:1])
+      [for(p=idx(xings[i]))
+        let(
+            crossings = deduplicate([
+                                     [p,0,0],
+                                     each xings[i][p],
+                                     [p,len(regions[i][p])-(closed[i]?1:2), 1],
+                                    ],eps=eps),
+            subpaths = [
+                for (frag = pair(crossings)) 
+                    deduplicate(
+                        _path_select(regions[i][p], frag[0][1], frag[0][2], frag[1][1], frag[1][2], closed=closed[i]),
+                        eps=eps
+                    )
+            ]
+        )
+        [for(s=subpaths) if (len(s)>1) s]
+       ]
+    ];
+                
+                
+
+
 // Function: split_path_at_region_crossings()
 // Usage:
 //   paths = split_path_at_region_crossings(path, region, [eps]);
@@ -715,13 +826,14 @@ function _point_dist(path,pathseg_unit,pathseg_len,pt) =
 
 function _offset_region(region, r, delta, chamfer, check_valid, quality,closed,return_faces,firstface_index,flip_faces) =
     let(
-         reglist = [for(R=region_parts(region)) is_path(R) ? [R] : R],
+         reglist = [for(R=region_parts(region)) force_region(R)],
          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)]]
+             difference([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,
@@ -1020,10 +1132,41 @@ function _tag_subpaths(region1, region2, keep, eps=EPSILON) =
 
 
 
-function _tagged_region(region1,region2,keep1,keep2,eps=EPSILON) =
-    _assemble_path_fragments(concat(_tag_subpaths(region1, region2, keep1, eps=eps),
-                                    _tag_subpaths(region2, region1, keep2, eps=eps)),
-                             eps=eps);
+
+function _keep_some_region_parts(region1, region2, keep1, keep2, eps=EPSILON) = 
+    // 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(
+        keep = [keep1,keep2],
+        subpaths = split_region_at_region_crossings(region1,region2,eps=eps),
+        regions=[region1,region2]
+    )        
+    _assemble_path_fragments(
+    [for(i=[0:1])
+      let(
+          keepS = search("S",keep[i])!=[],
+          keepU = search("U",keep[i])!=[],        
+          keepoutside = search("O",keep[i]) !=[],
+          keepinside = search("I",keep[i]) !=[],
+          all_subpaths = flatten(subpaths[i])
+      )
+      for (subpath = all_subpaths)
+          let(
+                midpt = mean([subpath[0], subpath[1]]),
+                rel = point_in_region(midpt,regions[1-i],eps=eps),
+                keepthis = rel<0 ? keepoutside
+                         : rel>0 ? keepinside
+                         : !(keepS || keepU) ? false
+                         : let(
+                               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 ? keepS : keepU
+            )
+            if (keepthis) subpath
+    ]);
 
 
 // Function&Module: union()
@@ -1049,7 +1192,7 @@ function union(regions=[],b=undef,c=undef,eps=EPSILON) =
     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],"OS", "O", eps=eps),
+           _keep_some_region_parts(regions[0],regions[1],"OS", "O", eps=eps),           
             for (i=[2:1:len(regions)-1]) regions[i]
           ],
           eps=eps
@@ -1081,7 +1224,7 @@ function difference(regions=[],b=undef,c=undef,eps=EPSILON) =
     regions[0]==[] ? [] : 
     let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
     difference([
-                _tagged_region(regions[0],regions[1],"OU", "I", eps=eps),
+                _keep_some_region_parts(regions[0],regions[1],"OU", "I", eps=eps),                
                 for (i=[2:1:len(regions)-1]) regions[i]
                ],
                eps=eps
@@ -1112,7 +1255,7 @@ function intersection(regions=[],b=undef,c=undef,eps=EPSILON) =
    : regions[0]==[] || regions[1]==[] ? []   
    : let(regions=[for (r=regions) quant(is_path(r)? [r] : r, 1/65536)])
          intersection([
-                       _tagged_region(regions[0],regions[1],"IS","I",eps=eps),
+                       _keep_some_region_parts(regions[0],regions[1],"IS","I",eps=eps),                       
                        for (i=[2:1:len(regions)-1]) regions[i]
                       ],
                       eps=eps
@@ -1150,7 +1293,7 @@ function exclusive_or(regions=[],b=undef,c=undef,eps=EPSILON) =
     len(regions)==1? regions[0] :
     let(regions=[for (r=regions) is_path(r)? [r] : r])
     exclusive_or([
-                  _tagged_region(regions[0],regions[1],"IO","IO",eps=eps),
+                  _keep_some_region_parts(regions[0],regions[1],"IO","IO",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 4a2f396..90ada2e 100644
--- a/tests/test_regions.scad
+++ b/tests/test_regions.scad
@@ -17,7 +17,12 @@ 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)]));
-  assert(are_regions_equal(union(R2,R1), [square(10,center=true)]));  
+  assert(are_regions_equal(union(R2,R1), [square(10,center=true)]));
+  R8 = [right(8,square(10,center=true)), left(8,square(10,center=true))];
+  R9 = [back(8,square(10,center=true)), fwd(8,square(10,center=true))];
+  assert(are_regions_equal(union(R9,R8), [[[-5, -5], [-13, -5], [-13, 5], [-5, 5], [-5, 13], [5, 13], [5, 5], [13, 5], [13, -5], [5, -5], [5, -13], [-5, -13]], [[-3, 3], [-3, -3], [3, -3], [3, 3]]]));
+  assert(are_regions_equal(union(R8,R9), [[[-5, -5], [-13, -5], [-13, 5], [-5, 5], [-5, 13], [5, 13], [5, 5], [13, 5], [13, -5], [5, -5], [5, -13], [-5, -13]], [[-3, 3], [-3, -3], [3, -3], [3, 3]]]));
+
 }
 test_union();
 
@@ -27,6 +32,12 @@ module test_intersection() {
   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));
+  R8 = [right(8,square(10,center=true)), left(8,square(10,center=true))];
+  R9 = [back(8,square(10,center=true)), fwd(8,square(10,center=true))];
+  assert(are_regions_equal(intersection(R9,R8),[[[-3, -5], [-5, -5], [-5, -3], [-3, -3]], [[-5, 5], [-3, 5], [-3, 3], [-5, 3]], [[5, -5], [3, -5], [3, -3], [5, -3]], [[3, 3], [3, 5], [5, 5], [5, 3]]]));
+  assert(are_regions_equal(intersection(R8,R9),[[[-3, -5], [-5, -5], [-5, -3], [-3, -3]], [[-5, 5], [-3, 5], [-3, 3], [-5, 3]], [[5, -5], [3, -5], [3, -3], [5, -3]], [[3, 3], [3, 5], [5, 5], [5, 3]]]));
+
+
 }
 test_intersection();
 
@@ -36,23 +47,60 @@ module test_difference() {
   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)]));
-
   pathA = [
     [-9,12], [-6,2], [-3,12], [0,2], [3,10], [5,10], [19,-4], [-8,-4], [-12,0]
   ];
-
   pathB = [
     [-12,8], [7,8], [9,6], [7,5], [-3,5], [-5,-6], [-2,-6], [0,-4],
     [6,-4], [2,-8], [-7,-8], [-15,0]
   ];
 
-
   right=[[[-10, 8], [-9, 12], [-7.8, 8]], [[0, -4], [-4.63636363636, -4], [-3, 5], [-0.9, 5], [0, 2], [1.125, 5], [7, 5], [9, 6], [19, -4], [6, -4]], [[-4.2, 8], [-1.8, 8], [-3, 12]], [[2.25, 8], [3, 10], [5, 10], [7, 8]]];
   assert(are_regions_equal(difference(pathA,pathB),right));
 
+  R8 = [right(8,square(10,center=true)), left(8,square(10,center=true))];
+  R9 = [back(8,square(10,center=true)), fwd(8,square(10,center=true))];
 
+  assert(are_regions_equal(difference(R9,R8), [[[-5, 5], [-5, 13], [5, 13], [5, 5], [3, 5], [3, 3], [-3, 3], [-3, 5]], [[5, -13], [-5, -13], [-5, -5], [-3, -5], [-3, -3], [3, -3], [3, -5], [5, -5]]]));
+
+  assert(are_regions_equal(difference(R8,R9),[[[-5, -5], [-13, -5], [-13, 5], [-5, 5], [-5, 3], [-3, 3], [-3, -3], [-5, -3]], [[3, -3], [3, 3], [5, 3], [5, 5], [13, 5], [13, -5], [5, -5], [5, -3]]]));
+
+  
 }
 test_difference();
 
 
 
+module test_exclusive_or() {
+  R8 = [right(8,square(10,center=true)), left(8,square(10,center=true))];
+  R9 = [back(8,square(10,center=true)), fwd(8,square(10,center=true))];
+  assert(are_regions_equal(exclusive_or(R8,R9),[[[-5, -5], [-13, -5], [-13, 5], [-5, 5], [-5, 3], [-3, 3], [-3, -3], [-5, -3]], [[-3, -5], [-5, -5], [-5, -13], [5, -13], [5, -5], [3, -5], [3, -3], [-3, -3]], [[-5, 5], [-3, 5], [-3, 3], [3, 3], [3, 5], [5, 5], [5, 13], [-5, 13]], [[3, -3], [3, 3], [5, 3], [5, 5], [13, 5], [13, -5], [5, -5], [5, -3]]],either_winding=true));
+  assert(are_regions_equal(exclusive_or(R9,R8),[[[-5, -5], [-13, -5], [-13, 5], [-5, 5], [-5, 3], [-3, 3], [-3, -3], [-5, -3]], [[-3, -5], [-5, -5], [-5, -13], [5, -13], [5, -5], [3, -5], [3, -3], [-3, -3]], [[-5, 5], [-3, 5], [-3, 3], [3, 3], [3, 5], [5, 5], [5, 13], [-5, 13]], [[3, -3], [3, 3], [5, 3], [5, 5], [13, 5], [13, -5], [5, -5], [5, -3]]],either_winding=true));  
+
+  p = turtle(["move",100,"left",144], repeat=4);
+  p2 = move(-polygon_centroid(p),p);
+  p3 = polygon_parts(p2);
+  p4 = exclusive_or(p3,square(51,center=true));
+
+  star_square = [[[-50, -16.2459848116], [-25.5, -16.2459848116],
+  [-25.5, 1.55430712449]], [[-7.45841874701, 25.5], [-30.9016994375,
+  42.5325404176], [-25.3674915789, 25.5]], [[-19.0983005625,
+  6.20541401733], [-25.5, 1.55430712449], [-25.5, 25.5],
+  [-25.3674915789, 25.5]], [[-11.803398875, -16.2459848116],
+  [-19.0983005625, 6.20541401733], [-3.5527136788e-15, 20.0811415886],
+  [19.0983005625, 6.20541401733], [11.803398875, -16.2459848116]],
+  [[7.45841874701, 25.5], [0, 20.0811415886], [-7.45841874701, 25.5]],
+  [[25.3674915789, 25.5], [7.45841874701, 25.5], [30.9016994375,
+  42.5325404176]], [[25.5, 1.55430712449], [19.0983005625,
+  6.20541401733], [25.3674915789, 25.5], [25.5, 25.5]], [[25.5,
+  -16.2459848116], [25.5, 1.55430712449], [50, -16.2459848116]],
+  [[8.79658707105, -25.5], [11.803398875, -16.2459848116], [25.5,
+  -16.2459848116], [25.5, -25.5]], [[-8.79658707105, -25.5],
+  [8.79658707105, -25.5], [0, -52.5731112119]], [[-25.5,
+  -16.2459848116], [-11.803398875, -16.2459848116], [-8.79658707105,
+  -25.5], [-25.5, -25.5]]];
+  assert(are_regions_equal(exclusive_or(p3,square(51,center=true)),star_square,either_winding=true));
+  assert(are_regions_equal(exclusive_or(square(51,center=true),p3),star_square,either_winding=true));
+
+}
+test_exclusive_or();