diff --git a/examples/boolean_geometry.scad b/examples/boolean_geometry.scad
new file mode 100644
index 0000000..151e361
--- /dev/null
+++ b/examples/boolean_geometry.scad
@@ -0,0 +1,48 @@
+include <BOSL2/std.scad>
+include <BOSL2/debug.scad>
+
+$fn = 36;
+
+rgn1 = [
+	square(100),
+	move([50,50], p=circle(d=60)),
+	[[40,40],[40,60],[60,60]]
+];
+
+rgn2 = [
+	[[0,0], [100,100], [100,0]],
+	[[27,10], [90,73], [90,10]],
+	move([70,30], p=circle(d=20))
+];
+
+
+module showit(label, rgn, poly=undef, outline=undef, width=0.75) {
+	move([-50,-50]) {
+		if(outline) color(outline) linear_extrude(height=max(0.1,1-width)) for(path=rgn) stroke(path, width=width, close=true);
+		if(poly) color(poly) linear_extrude(height=0.1) region(rgn);
+		color("black") right(50) fwd(7) linear_extrude(height=0.1) text(text=label, size=8, halign="center", valign="center");
+	}
+}
+
+
+ydistribute(-125) {
+	xdistribute(120) {
+		showit("Region A", rgn1, poly="green", outline="black");
+		showit("Region B", rgn2, poly="green", outline="black");
+		union() {
+			showit("A and B Overlaid", rgn1, outline="red",width=1);
+			showit("", rgn2, outline="blue",width=0.5);
+		}
+	}
+	xdistribute(120) {
+		showit("Union A+B", union(rgn1, rgn2), poly="green", outline="black");
+		showit("Difference A-B", difference(rgn1, rgn2), poly="green", outline="black");
+		showit("Difference B-A", difference(rgn2, rgn1), poly="green", outline="black");
+	}
+	xdistribute(120) {
+		showit("Intersection A&B", intersection(rgn1, rgn2), poly="green", outline="black");
+		showit("Exclusive OR A^B", exclusive_or(rgn1, rgn2), poly="green", outline="black");
+	}
+}
+
+// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/geometry.scad b/geometry.scad
index c0b976f..55f6fed 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -566,13 +566,18 @@ function close_region(region, eps=EPSILON) = [for (path=region) close_path(path,
 //   region = Region to test for crossings of.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function region_path_crossings(path, region, eps=EPSILON) = sort([
-	for (s1=enumerate(pair_wrap(path)), path=region, s2=pair_wrap(path)) let(
+	for (
+		s1=enumerate(pair(close_path(path))),
+		p=close_region(region),
+		s2=pair(p)
+	) let(
 		isect = _general_line_intersection(s1[1],s2,eps=eps)
 	) if (
 		!is_undef(isect) &&
-		isect[1] >= 0-eps && isect[1] < 1-eps &&
-		isect[2] >= 0-eps && isect[2] < 1-eps
-	) [s1[0], isect[1]]
+		isect[1] >= 0-eps && isect[1] < 1+eps &&
+		isect[2] >= 0-eps && isect[2] < 1+eps
+	)
+	[s1[0], isect[1]]
 ]);
 
 
@@ -580,15 +585,23 @@ function _split_path_at_region_crossings(path, region, eps=EPSILON) =
 	let(
 		path = deduplicate(path, eps=eps),
 		region = [for (path=region) deduplicate(path, eps=eps)],
+		xings = region_path_crossings(path, region, eps=eps),
 		crossings = deduplicate(
 			concat(
 				[[0,0]],
-				region_path_crossings(path, region),
+				xings,
 				[[len(path)-2,1]]
 			),
 			eps=eps
-		)
-	) [for (p = pair(crossings)) path_subselect(path, p[0][0], p[0][1], p[1][0], p[1][1])];
+		),
+		subpaths = [
+			for (p = pair(crossings))
+				deduplicate(eps=eps,
+					path_subselect(path, p[0][0], p[0][1], p[1][0], p[1][1])
+				)
+		]
+	)
+	subpaths;
 
 
 function _tag_subpaths(path, region, eps=EPSILON) =