diff --git a/affine.scad b/affine.scad
index 07c662c..c222ffe 100644
--- a/affine.scad
+++ b/affine.scad
@@ -249,7 +249,7 @@ function rot_decode(M,long=false) =
 // Usage:
 //   B = rot_inverse(A)
 // Description:
-//   Inverts a 2d or 3d rotation matrix.  The matrix can be a rotation around any center,
+//   Inverts a 2d (3x3) or 3d (4x4) rotation matrix.  The matrix can be a rotation around any center,
 //   so it may include a translation.
 function rot_inverse(T) =
     assert(is_matrix(T,square=true),"Matrix must be square")
diff --git a/debug.scad b/debug.scad
index 5d02c7e..8f25187 100644
--- a/debug.scad
+++ b/debug.scad
@@ -6,266 +6,6 @@
 //////////////////////////////////////////////////////////////////////
 
 
-// Section: Debugging Paths and Polygons
-
-// Module: trace_path()
-// Usage:
-//   trace_path(path, [closed=], [showpts=], [N=], [size=], [color=]);
-// Description:
-//   Renders lines between each point of a path.
-//   Can also optionally show the individual vertex points.
-// Arguments:
-//   path = The list of points in the path.
-//   ---
-//   closed = If true, draw the segment from the last vertex to the first.  Default: false
-//   showpts = If true, draw vertices and control points.
-//   N = Mark the first and every Nth vertex after in a different color and shape.
-//   size = Diameter of the lines drawn.
-//   color = Color to draw the lines (but not vertices) in.
-// Example(FlatSpin,VPD=44.4):
-//   path = [for (a=[0:30:210]) 10*[cos(a), sin(a), sin(a)]];
-//   trace_path(path, showpts=true, size=0.5, color="lightgreen");
-module trace_path(path, closed=false, showpts=false, N=1, size=1, color="yellow") {
-    assert(is_path(path),"Invalid path argument");
-    sides = segs(size/2);
-    path = closed? close_path(path) : path;
-    if (showpts) {
-        for (i = [0:1:len(path)-1]) {
-            translate(path[i]) {
-                if (i % N == 0) {
-                    color("blue") sphere(d=size*2.5, $fn=8);
-                } else {
-                    color("red") {
-                        cylinder(d=size/2, h=size*3, center=true, $fn=8);
-                        xrot(90) cylinder(d=size/2, h=size*3, center=true, $fn=8);
-                        yrot(90) cylinder(d=size/2, h=size*3, center=true, $fn=8);
-                    }
-                }
-            }
-        }
-    }
-    if (N!=3) {
-        color(color) stroke(path3d(path), width=size, $fn=8);
-    } else {
-        for (i = [0:1:len(path)-2]) {
-            if (N != 3 || (i % N) != 1) {
-                color(color) extrude_from_to(path[i], path[i+1]) circle(d=size, $fn=sides);
-            }
-        }
-    }
-}
-
-
-// Module: debug_polygon()
-// Usage:
-//   debug_polygon(points, paths, [convexity=], [size=]);
-// Description:
-//   A drop-in replacement for `polygon()` that renders and labels the path points.
-// Arguments:
-//   points = The array of 2D polygon vertices.
-//   paths = The path connections between the vertices.
-//   ---
-//   convexity = The max number of walls a ray can pass through the given polygon paths.
-//   size = The base size of the line and labels.
-// Example(Big2D):
-//   debug_polygon(
-//       points=concat(
-//           regular_ngon(or=10, n=8),
-//           regular_ngon(or=8, n=8)
-//       ),
-//       paths=[
-//           [for (i=[0:7]) i],
-//           [for (i=[15:-1:8]) i]
-//       ]
-//   );
-module debug_polygon(points, paths, convexity=2, size=1)
-{
-    paths = is_undef(paths)? [[for (i=[0:1:len(points)-1]) i]] :
-        is_num(paths[0])? [paths] :
-        paths;
-    echo(points=points);
-    echo(paths=paths);
-    linear_extrude(height=0.01, convexity=convexity, center=true) {
-        polygon(points=points, paths=paths, convexity=convexity);
-    }
-    for (i = [0:1:len(points)-1]) {
-        color("red") {
-            up(0.2) {
-                translate(points[i]) {
-                    linear_extrude(height=0.1, convexity=10, center=true) {
-                        text(text=str(i), size=size, halign="center", valign="center");
-                    }
-                }
-            }
-        }
-    }
-    for (j = [0:1:len(paths)-1]) {
-        path = paths[j];
-        translate(points[path[0]]) {
-            color("cyan") up(0.1) cylinder(d=size*1.5, h=0.01, center=false, $fn=12);
-        }
-        translate(points[path[len(path)-1]]) {
-            color("pink") up(0.11) cylinder(d=size*1.5, h=0.01, center=false, $fn=4);
-        }
-        for (i = [0:1:len(path)-1]) {
-            midpt = (points[path[i]] + points[path[(i+1)%len(path)]])/2;
-            color("blue") {
-                up(0.2) {
-                    translate(midpt) {
-                        linear_extrude(height=0.1, convexity=10, center=true) {
-                            text(text=str(chr(65+j),i), size=size/2, halign="center", valign="center");
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-
-// Section: Debugging Polyhedrons
-
-
-// Module: debug_vertices()
-// Usage:
-//   debug_vertices(vertices, [size], [disabled=]);
-// Description:
-//   Draws all the vertices in an array, at their 3D position, numbered by their
-//   position in the vertex array.  Also draws any children of this module with
-//   transparency.
-// Arguments:
-//   vertices = Array of point vertices.
-//   size = The size of the text used to label the vertices.  Default: 1
-//   ---
-//   disabled = If true, don't draw numbers, and draw children without transparency.  Default = false.
-// Example:
-//   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
-//   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
-//   debug_vertices(vertices=verts, size=2) {
-//       polyhedron(points=verts, faces=faces);
-//   }
-module debug_vertices(vertices, size=1, disabled=false) {
-    if (!disabled) {
-        color("blue") {
-            dups = vector_search(vertices, EPSILON, vertices);
-            for (ind = dups){
-                numstr = str_join([for(i=ind) str(i)],",");
-                v = vertices[ind[0]];
-                translate(v) {
-                    up(size/8) zrot($vpr[2]) xrot(90) {
-                        linear_extrude(height=size/10, center=true, convexity=10) {
-                            text(text=numstr, size=size, halign="center");
-                        }
-                    }
-                    sphere(size/10);
-                }
-            }
-        }
-    }
-    if ($children > 0) {
-        if (!disabled) {
-            color([0.2, 1.0, 0, 0.5]) children();
-        } else {
-            children();
-        }
-    }
-}
-
-
-
-// Module: debug_faces()
-// Usage:
-//   debug_faces(vertices, faces, [size=], [disabled=]);
-// Description:
-//   Draws all the vertices at their 3D position, numbered in blue by their
-//   position in the vertex array.  Each face will have their face number drawn
-//   in red, aligned with the center of face.  All children of this module are drawn
-//   with transparency.
-// Arguments:
-//   vertices = Array of point vertices.
-//   faces = Array of faces by vertex numbers.
-//   ---
-//   size = The size of the text used to label the faces and vertices.  Default: 1
-//   disabled = If true, don't draw numbers, and draw children without transparency.  Default: false.
-// Example(EdgesMed):
-//   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
-//   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
-//   debug_faces(vertices=verts, faces=faces, size=2) {
-//       polyhedron(points=verts, faces=faces);
-//   }
-module debug_faces(vertices, faces, size=1, disabled=false) {
-    if (!disabled) {
-        vlen = len(vertices);
-        color("red") {
-            for (i = [0:1:len(faces)-1]) {
-                face = faces[i];
-                if (face[0] < 0 || face[1] < 0 || face[2] < 0 || face[0] >= vlen || face[1] >= vlen || face[2] >= vlen) {
-                    echo("BAD FACE: ", vlen=vlen, face=face);
-                } else {
-                    verts = select(vertices,face);
-                    c = mean(verts);
-                    v0 = verts[0];
-                    v1 = verts[1];
-                    v2 = verts[2];
-                    dv0 = unit(v1 - v0);
-                    dv1 = unit(v2 - v0);
-                    nrm0 = cross(dv0, dv1);
-                    nrm1 = UP;
-                    axis = vector_axis(nrm0, nrm1);
-                    ang = vector_angle(nrm0, nrm1);
-                    theta = atan2(nrm0[1], nrm0[0]);
-                    translate(c) {
-                        rotate(a=180-ang, v=axis) {
-                            zrot(theta-90)
-                            linear_extrude(height=size/10, center=true, convexity=10) {
-                                union() {
-                                    text(text=str(i), size=size, halign="center");
-                                    text(text=str("_"), size=size, halign="center");
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-    debug_vertices(vertices, size=size, disabled=disabled) {
-        children();
-    }
-    if (!disabled) {
-        echo(faces=faces);
-    }
-}
-
-
-
-// Module: debug_vnf()
-// Usage:
-//   debug_vnf(vnfs, [convexity=], [txtsize=], [disabled=]);
-// Description:
-//   A drop-in module to replace `vnf_polyhedron()` and help debug vertices and faces.
-//   Draws all the vertices at their 3D position, numbered in blue by their
-//   position in the vertex array.  Each face will have its face number drawn
-//   in red, aligned with the center of face.  All given faces are drawn with
-//   transparency. All children of this module are drawn with transparency.
-//   Works best with Thrown-Together preview mode, to see reversed faces.
-// Arguments:
-//   vnf = vnf to display
-//   ---
-//   convexity = The max number of walls a ray can pass through the given polygon paths.
-//   txtsize = The size of the text used to label the faces and vertices.
-//   disabled = If true, act exactly like `polyhedron()`.  Default = false.
-// Example(EdgesMed):
-//   verts = [for (z=[-10,10], a=[0:120:359.9]) [10*cos(a),10*sin(a),z]];
-//   faces = [[0,1,2], [5,4,3], [0,3,4], [0,4,1], [1,4,5], [1,5,2], [2,5,3], [2,3,0]];
-//   debug_vnf([verts,faces], txtsize=2);
-module debug_vnf(vnf, convexity=6, txtsize=1, disabled=false) {
-    debug_faces(vertices=vnf[0], faces=vnf[1], size=txtsize, disabled=disabled) {
-        vnf_polyhedron(vnf, convexity=convexity);
-    }
-}
-
 
 // Function: standard_anchors()
 // Usage:
@@ -605,73 +345,6 @@ function random_polygon(n=3,size=1, seed) =
     [for(i=count(n)) rads[i]*[cos(angs[i]), sin(angs[i])] ];
 
 
-// Function: random_points()
-// Usage:
-//    points = random_points(n, dim, scale, [seed]);
-// See Also: random_polygon(), gaussian_random_points(), spherical_random_points()
-// Topics: Random, Points
-// Description:
-//    Generate `n` random points of dimension `dim` with coordinates absolute value less than `scale`.
-//    The `scale` may be a number or a vector with dimension `dim`.
-// Arguments:
-//    n = number of points to generate. 
-//    dim = dimension of the points. Default: 2
-//    scale = the scale of the point coordinates. Default: 1
-//    seed = an optional seed for the random generation.
-function random_points(n, dim=2, scale=1, seed) =
-    assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
-    assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
-    assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.")
-    let( 
-        rnds =   is_undef(seed) 
-                ? rands(-1,1,n*dim)
-                : rands(-1,1,n*dim, seed) )
-    is_num(scale) 
-    ? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]
-    : [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ];
-
-
-// Function: gaussian_random_points()
-// Usage:
-//    points = gaussian_random_points(n, dim, mean, stddev, [seed]);
-// See Also: random_polygon(), random_points(), spherical_random_points()
-// Topics: Random, Points
-// Description:
-//    Generate `n` random points of dimension `dim` with coordinates absolute value less than `scale`.
-//    The gaussian distribution of all the coordinates of the points will have a mean `mean` and 
-//    standard deviation `stddev` 
-// Arguments:
-//    n = number of points to generate. 
-//    dim = dimension of the points. Default: 2
-//    mean = the gaussian mean of the point coordinates. Default: 0
-//    stddev = the gaussian standard deviation of the point coordinates. Default: 0
-//    seed = an optional seed for the random generation.
-function gaussian_random_points(n, dim=2, mean=0, stddev=1, seed) =
-    assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
-    assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
-    let( rnds = gaussian_rands(mean, stddev, n*dim, seed=seed) )
-    [for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ];
-
-
-// Function: spherical_random_points()
-// Usage:
-//    points = spherical_random_points(n, radius, [seed]);
-// See Also: random_polygon(), random_points(), gaussian_random_points()
-// Topics: Random, Points
-// Description:
-//    Generate `n` 3D random points lying on a sphere centered at the origin with radius equal to `radius`.
-// Arguments:
-//    n = number of points to generate. 
-//    radius = the sphere radius. Default: 1
-//    seed = an optional seed for the random generation.
-function spherical_random_points(n, radius=1, seed) =
-    assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
-    assert( is_num(radius) && radius>0, "The radius should be a non-negative number.")
-    let( rnds = is_undef(seed) 
-                ? rands(-1,1,n*2)
-                : rands(-1,1,n*2, seed) )
-    [for(i=[0:1:n-1]) spherical_to_xyz(radius, theta=180*rnds[2*i], phi=180*rnds[2*i+1])  ];
-
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/drawing.scad b/drawing.scad
index 21daad8..1ff77d9 100644
--- a/drawing.scad
+++ b/drawing.scad
@@ -445,12 +445,14 @@ module stroke(
 // Topics: Paths, Drawing Tools
 // See Also: stroke(), path_cut()
 // Description:
-//   Given a path and a dash pattern, creates a dashed line that follows that
-//   path with the given dash pattern.
+//   Given a path (or region) and a dash pattern, creates a dashed line that follows that
+//   path or region boundary with the given dash pattern.
 //   - When called as a function, returns a list of dash sub-paths.
 //   - When called as a module, draws all those subpaths using `stroke()`.
+//   When called as a module the dash pattern is multiplied by the line width.  When called as
+//   a function the dash pattern applies as you specify it.  
 // Arguments:
-//   path = The path to subdivide into dashes.
+//   path = The path or region to subdivide into dashes.
 //   dashpat = A list of alternating dash lengths and space lengths for the dash pattern.  This will be scaled by the width of the line.
 //   ---
 //   width = The width of the dashed line to draw.  Module only.  Default: 1
@@ -466,6 +468,7 @@ module stroke(
 //   path = [for (a=[-180:5:180]) [a/3, 20*cos(3*a), 20*sin(3*a)]];
 //   dashed_stroke(path, [3,2], width=1);
 function dashed_stroke(path, dashpat=[3,3], closed=false) =
+    is_region(path) ? [for(p=path) each dashed_stroke(p,dashpat,closed=true)] : 
     let(
         path = closed? close_path(path) : path,
         dashpat = len(dashpat)%2==0? dashpat : concat(dashpat,[0]),
@@ -491,6 +494,55 @@ module dashed_stroke(path, dashpat=[3,3], width=1, closed=false) {
 }
 
 
+
+// Module: trace_path()
+// Usage:
+//   trace_path(path, [closed=], [showpts=], [N=], [size=], [color=]);
+// Description:
+//   Renders lines between each point of a path.
+//   Can also optionally show the individual vertex points.
+// Arguments:
+//   path = The list of points in the path.
+//   ---
+//   closed = If true, draw the segment from the last vertex to the first.  Default: false
+//   showpts = If true, draw vertices and control points.
+//   N = Mark the first and every Nth vertex after in a different color and shape.
+//   size = Diameter of the lines drawn.
+//   color = Color to draw the lines (but not vertices) in.
+// Example(FlatSpin,VPD=44.4):
+//   path = [for (a=[0:30:210]) 10*[cos(a), sin(a), sin(a)]];
+//   trace_path(path, showpts=true, size=0.5, color="lightgreen");
+module trace_path(path, closed=false, showpts=false, N=1, size=1, color="yellow") {
+    assert(is_path(path),"Invalid path argument");
+    sides = segs(size/2);
+    path = closed? close_path(path) : path;
+    if (showpts) {
+        for (i = [0:1:len(path)-1]) {
+            translate(path[i]) {
+                if (i % N == 0) {
+                    color("blue") sphere(d=size*2.5, $fn=8);
+                } else {
+                    color("red") {
+                        cylinder(d=size/2, h=size*3, center=true, $fn=8);
+                        xrot(90) cylinder(d=size/2, h=size*3, center=true, $fn=8);
+                        yrot(90) cylinder(d=size/2, h=size*3, center=true, $fn=8);
+                    }
+                }
+            }
+        }
+    }
+    if (N!=3) {
+        color(color) stroke(path3d(path), width=size, $fn=8);
+    } else {
+        for (i = [0:1:len(path)-2]) {
+            if (N != 3 || (i % N) != 1) {
+                color(color) extrude_from_to(path[i], path[i+1]) circle(d=size, $fn=sides);
+            }
+        }
+    }
+}
+
+
 // Section: Computing paths
 
 // Function&Module: arc()
diff --git a/geometry.scad b/geometry.scad
index 7dab6e8..7ef4bc0 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -1865,6 +1865,54 @@ function align_polygon(reference, poly, angles, cp) =
     ) alignments[best][0];
 
 
+// Function: are_polygons_equal()
+// Usage:
+//    b = are_polygons_equal(poly1, poly2, [eps])
+// Description:
+//    Returns true if poly1 and poly2 are the same polongs
+//    within given epsilon tolerance.
+// Arguments:
+//    poly1 = first polygon
+//    poly2 = second polygon
+//    eps = tolerance for comparison
+// Example(NORENDER):
+//    are_polygons_equal(pentagon(r=4),
+//                   rot(360/5, p=pentagon(r=4))); // returns true
+//    are_polygons_equal(pentagon(r=4),
+//                   rot(90, p=pentagon(r=4)));    // returns false
+function are_polygons_equal(poly1, poly2, eps=EPSILON) =
+    let(
+        poly1 = cleanup_path(poly1),
+        poly2 = cleanup_path(poly2),
+        l1 = len(poly1),
+        l2 = len(poly2)
+    ) l1 != l2 ? false :
+    let( maybes = find_first_match(poly1[0], poly2, eps=eps, all=true) )
+    maybes == []? false :
+    [for (i=maybes) if (_are_polygons_equal(poly1, poly2, eps, i)) 1] != [];
+
+function _are_polygons_equal(poly1, poly2, eps, st) =
+    max([for(d=poly1-select(poly2,st,st-1)) d*d])<eps*eps;
+
+
+// Function: is_polygon_in_list()
+// Topics: Polygons, Comparators
+// See Also: are_polygons_equal(), are_regions_equal()
+// Usage:
+//   bool = is_polygon_in_list(poly, polys);
+// Description:
+//   Returns true if one of the polygons in `polys` is equivalent to the polygon `poly`.
+// Arguments:
+//   poly = The polygon to search for.
+//   polys = The list of polygons to look for the polygon in.
+function is_polygon_in_list(poly, polys) =
+    __is_polygon_in_list(poly, polys, 0);
+
+function __is_polygon_in_list(poly, polys, i) =
+    i >= len(polys)? false :
+    are_polygons_equal(poly, polys[i])? true :
+    __is_polygon_in_list(poly, polys, i+1);
+
 
 
 // Section: Convex Sets
diff --git a/math.scad b/math.scad
index 0e15a4d..fe3beae 100644
--- a/math.scad
+++ b/math.scad
@@ -496,6 +496,33 @@ function rand_int(minval, maxval, N, seed=undef) =
     [for(entry = rvect) floor(entry)];
 
 
+// Function: random_points()
+// Usage:
+//    points = random_points(n, dim, scale, [seed]);
+// See Also: random_polygon(), gaussian_random_points(), spherical_random_points()
+// Topics: Random, Points
+// Description:
+//    Generate `n` uniform random points of dimension `dim` with data ranging from -scale to +scale.  
+//    The `scale` may be a number, in which case the random data lies in a cube,
+//    or a vector with dimension `dim`, in which case each dimension has its own scale.  
+// Arguments:
+//    n = number of points to generate. 
+//    dim = dimension of the points. Default: 2
+//    scale = the scale of the point coordinates. Default: 1
+//    seed = an optional seed for the random generation.
+function random_points(n, dim=2, scale=1, seed) =
+    assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
+    assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
+    assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.")
+    let( 
+        rnds =   is_undef(seed) 
+                ? rands(-1,1,n*dim)
+                : rands(-1,1,n*dim, seed) )
+    is_num(scale) 
+    ? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]
+    : [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ];
+
+
 // Function: gaussian_rands()
 // Usage:
 //   arr = gaussian_rands(mean, stddev, [N], [seed]);
@@ -512,6 +539,56 @@ function gaussian_rands(mean, stddev, N=1, seed=undef) =
     [for (i = count(N,0,2)) mean + stddev*sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])];
 
 
+// Function: gaussian_random_points()
+// Usage:
+//    points = gaussian_random_points(n, dim, mean, stddev, [seed]);
+// See Also: random_polygon(), random_points(), spherical_random_points()
+// Topics: Random, Points
+// Description:
+//    Generate `n` random points of dimension `dim` with coordinates absolute value less than `scale`.
+//    The gaussian distribution of all the coordinates of the points will have a mean `mean` and 
+//    standard deviation `stddev` 
+// Arguments:
+//    n = number of points to generate. 
+//    dim = dimension of the points. Default: 2
+//    mean = the gaussian mean of the point coordinates. Default: 0
+//    stddev = the gaussian standard deviation of the point coordinates. Default: 0
+//    seed = an optional seed for the random generation.
+function gaussian_random_points(n, dim=2, mean=0, stddev=1, seed) =
+    assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
+    assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
+    let( rnds = gaussian_rands(mean, stddev, n*dim, seed=seed) )
+    [for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ];
+
+
+// Function: spherical_random_points()
+// Usage:
+//    points = spherical_random_points(n, radius, [seed]);
+// See Also: random_polygon(), random_points(), gaussian_random_points()
+// Topics: Random, Points
+// Description:
+//    Generate `n` 3D uniformly distributed random points lying on a sphere centered at the origin with radius equal to `radius`.
+// Arguments:
+//    n = number of points to generate. 
+//    radius = the sphere radius. Default: 1
+//    seed = an optional seed for the random generation.
+
+// See https://mathworld.wolfram.com/SpherePointPicking.html
+function spherical_random_points(n, radius=1, seed) =
+    assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
+    assert( is_num(radius) && radius>0, "The radius should be a non-negative number.")
+    let( theta = is_undef(seed) 
+                ? rands(0,360,n)
+                : rands(0,360,n, seed),
+         cosphi = rands(-1,1,n))
+    [for(i=[0:1:n-1]) let(
+                          sin_phi=sqrt(1-cosphi[i]*cosphi[i])
+                      )
+                      radius*[sin_phi*cos(theta[i]),sin_phi*sin(theta[i]), cosphi[i]]];
+
+
+
+
 // Function: log_rands()
 // Usage:
 //   num = log_rands(minval, maxval, factor, [N], [seed]);
diff --git a/paths.scad b/paths.scad
index d346b96..21243e5 100644
--- a/paths.scad
+++ b/paths.scad
@@ -125,55 +125,6 @@ function path_merge_collinear(path, closed=false, eps=EPSILON) =
     ) [for (i=indices) path[i]];
 
 
-// Function: are_polygons_equal()
-// Usage:
-//    b = are_polygons_equal(poly1, poly2, [eps])
-// Description:
-//    Returns true if poly1 and poly2 are the same polongs
-//    within given epsilon tolerance.
-// Arguments:
-//    poly1 = first polygon
-//    poly2 = second polygon
-//    eps = tolerance for comparison
-// Example(NORENDER):
-//    are_polygons_equal(pentagon(r=4),
-//                   rot(360/5, p=pentagon(r=4))); // returns true
-//    are_polygons_equal(pentagon(r=4),
-//                   rot(90, p=pentagon(r=4)));    // returns false
-function are_polygons_equal(poly1, poly2, eps=EPSILON) =
-    let(
-        poly1 = cleanup_path(poly1),
-        poly2 = cleanup_path(poly2),
-        l1 = len(poly1),
-        l2 = len(poly2)
-    ) l1 != l2 ? false :
-    let( maybes = find_first_match(poly1[0], poly2, eps=eps, all=true) )
-    maybes == []? false :
-    [for (i=maybes) if (_are_polygons_equal(poly1, poly2, eps, i)) 1] != [];
-
-function _are_polygons_equal(poly1, poly2, eps, st) =
-    max([for(d=poly1-select(poly2,st,st-1)) d*d])<eps*eps;
-
-
-// Function: is_polygon_in_list()
-// Topics: Polygons, Comparators
-// See Also: are_polygons_equal(), are_regions_equal()
-// Usage:
-//   bool = is_polygon_in_list(poly, polys);
-// Description:
-//   Returns true if one of the polygons in `polys` is equivalent to the polygon `poly`.
-// Arguments:
-//   poly = The polygon to search for.
-//   polys = The list of polygons to look for the polygon in.
-function is_polygon_in_list(poly, polys) =
-    __is_polygon_in_list(poly, polys, 0);
-
-function __is_polygon_in_list(poly, polys, i) =
-    i >= len(polys)? false :
-    are_polygons_equal(poly, polys[i])? true :
-    __is_polygon_in_list(poly, polys, i+1);
-
-
 
 // Section: Path length calculation
 
diff --git a/regions.scad b/regions.scad
index 83ad84a..76dd863 100644
--- a/regions.scad
+++ b/regions.scad
@@ -35,26 +35,6 @@
 function is_region(x) = is_list(x) && is_path(x.x);
 
 
-// Function: close_region()
-// Usage:
-//   close_region(region);
-// Description:
-//   Closes all paths within a given region.
-function close_region(region, eps=EPSILON) = [for (path=region) close_path(path, eps=eps)];
-
-
-// Function: cleanup_region()
-// Usage:
-//   cleanup_region(region);
-// Description:
-//   For all paths in the given region, if the last point coincides with the first point, removes the last point.
-// Arguments:
-//   region = The region to clean up.  Given as a list of polygon paths.
-//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
-function cleanup_region(region, eps=EPSILON) =
-    [for (path=region) cleanup_path(path, eps=eps)];
-
-
 
 // Function: check_and_fix_path()
 // Usage:
@@ -329,7 +309,24 @@ function _join_paths_at_vertices(path1,path2,seg1,seg2) =
     ) cleanup_path(deduplicate([each path1, each path2]));
 
 
-function _cleave_simple_region(region) =
+function new_join_paths_at_vertices(path1,path2,v1,v2) =
+    let(
+        repeat_start = !approx(path1[v1],path2[v2]),
+        path1 = clockwise_polygon(polygon_shift(path1,v1)),
+        path2 = ccw_polygon(polygon_shift(path2,v2))
+    )
+    [
+        each path1,
+        if (repeat_start) path1[0],
+        each path2,
+        if (repeat_start) path2[0],
+    ];                      
+
+        
+// Given a region that is connected and has its outer border in region[0],
+// produces a polygon with the same points that has overlapping connected paths
+// to join internal holes to the outer border.  
+function _cleave_connected_region(region) =
     len(region)==0? [] :
     len(region)<=1? clockwise_polygon(region[0]) :
     let(
@@ -338,7 +335,36 @@ function _cleave_simple_region(region) =
             _path_path_closest_vertices(region[0],region[i])
         ],
         idxi = min_index(subindex(dists,0)),
-        newoline = _join_paths_at_vertices(
+        outline = _join_paths_at_vertices(
+            region[0], region[idxi+1],
+            dists[idxi][1], dists[idxi][2]
+        )
+    )
+    len(region)==2? clockwise_polygon(outline) :   // We joined 2 regions, so we're done
+    let(
+        newregion = [
+            outline,
+            for (i=idx(region))
+                if (i>0 && i!=idxi+1)
+                    region[i]
+        ]
+    )
+    assert(len(newregion)<len(region))
+    _cleave_connected_region(newregion);
+
+
+
+
+function new_cleave_connected_region(region) =
+    len(region)==0? [] :
+    len(region)<=1? clockwise_polygon(region[0]) :
+    let(
+        dists = [
+            for (i=[1:1:len(region)-1])
+            _path_path_closest_vertices(region[0],region[i])
+        ],
+        idxi = min_index(subindex(dists,0)),
+        newoline = new_join_paths_at_vertices(
             region[0], region[idxi+1],
             dists[idxi][1], dists[idxi][2]
         )
@@ -352,8 +378,8 @@ function _cleave_simple_region(region) =
         ]
     )
     assert(len(orgn)<len(region))
-    _cleave_simple_region(orgn);
-
+    new_cleave_connected_region(orgn);
+        
 
 // Function: region_faces()
 // Usage:
@@ -372,7 +398,7 @@ function region_faces(region, transform, reverse=false, vnf=EMPTY_VNF) =
         vnfs = [
             if (vnf != EMPTY_VNF) vnf,
             for (rgn = regions) let(
-                cleaved = path3d(_cleave_simple_region(rgn)),
+                cleaved = path3d(_cleave_connected_region(rgn)),
                 face = is_undef(transform)? cleaved : apply(transform,cleaved),
                 faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i]
             ) [face, [faceidxs]]
diff --git a/shapes2d.scad b/shapes2d.scad
index d99ed29..00e1445 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -1639,4 +1639,76 @@ function mask2d_ogee(pattern, excess=0.01, anchor=CENTER, spin=0) =
 
 
 
+// Section: Debugging polygons
+
+// Module: debug_polygon()
+// Usage:
+//   debug_polygon(points, paths, [convexity=], [size=]);
+// Description:
+//   A drop-in replacement for `polygon()` that renders and labels the path points.
+// Arguments:
+//   points = The array of 2D polygon vertices.
+//   paths = The path connections between the vertices.
+//   ---
+//   convexity = The max number of walls a ray can pass through the given polygon paths.
+//   size = The base size of the line and labels.
+// Example(Big2D):
+//   debug_polygon(
+//       points=concat(
+//           regular_ngon(or=10, n=8),
+//           regular_ngon(or=8, n=8)
+//       ),
+//       paths=[
+//           [for (i=[0:7]) i],
+//           [for (i=[15:-1:8]) i]
+//       ]
+//   );
+module debug_polygon(points, paths, convexity=2, size=1)
+{
+    paths = is_undef(paths)? [[for (i=[0:1:len(points)-1]) i]] :
+        is_num(paths[0])? [paths] :
+        paths;
+    echo(points=points);
+    echo(paths=paths);
+    linear_extrude(height=0.01, convexity=convexity, center=true) {
+        polygon(points=points, paths=paths, convexity=convexity);
+    }
+    for (i = [0:1:len(points)-1]) {
+        color("red") {
+            up(0.2) {
+                translate(points[i]) {
+                    linear_extrude(height=0.1, convexity=10, center=true) {
+                        text(text=str(i), size=size, halign="center", valign="center");
+                    }
+                }
+            }
+        }
+    }
+    for (j = [0:1:len(paths)-1]) {
+        path = paths[j];
+        translate(points[path[0]]) {
+            color("cyan") up(0.1) cylinder(d=size*1.5, h=0.01, center=false, $fn=12);
+        }
+        translate(points[path[len(path)-1]]) {
+            color("pink") up(0.11) cylinder(d=size*1.5, h=0.01, center=false, $fn=4);
+        }
+        for (i = [0:1:len(path)-1]) {
+            midpt = (points[path[i]] + points[path[(i+1)%len(path)]])/2;
+            color("blue") {
+                up(0.2) {
+                    translate(midpt) {
+                        linear_extrude(height=0.1, convexity=10, center=true) {
+                            text(text=str(chr(65+j),i), size=size/2, halign="center", valign="center");
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+
+
+
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/vnf.scad b/vnf.scad
index 4d9177e..db9ca92 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -928,6 +928,150 @@ function _split_polygons_at_each_y(polys, ys, _i=0) =
 
 // Section: Debugging VNFs
 
+// Section: Debugging Polyhedrons
+
+
+// Module: debug_vertices()
+// Usage:
+//   debug_vertices(vertices, [size], [disabled=]);
+// Description:
+//   Draws all the vertices in an array, at their 3D position, numbered by their
+//   position in the vertex array.  Also draws any children of this module with
+//   transparency.
+// Arguments:
+//   vertices = Array of point vertices.
+//   size = The size of the text used to label the vertices.  Default: 1
+//   ---
+//   disabled = If true, don't draw numbers, and draw children without transparency.  Default = false.
+// Example:
+//   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
+//   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
+//   debug_vertices(vertices=verts, size=2) {
+//       polyhedron(points=verts, faces=faces);
+//   }
+module debug_vertices(vertices, size=1, disabled=false) {
+    if (!disabled) {
+        color("blue") {
+            dups = vector_search(vertices, EPSILON, vertices);
+            for (ind = dups){
+                numstr = str_join([for(i=ind) str(i)],",");
+                v = vertices[ind[0]];
+                translate(v) {
+                    up(size/8) zrot($vpr[2]) xrot(90) {
+                        linear_extrude(height=size/10, center=true, convexity=10) {
+                            text(text=numstr, size=size, halign="center");
+                        }
+                    }
+                    sphere(size/10);
+                }
+            }
+        }
+    }
+    if ($children > 0) {
+        if (!disabled) {
+            color([0.2, 1.0, 0, 0.5]) children();
+        } else {
+            children();
+        }
+    }
+}
+
+
+
+// Module: debug_faces()
+// Usage:
+//   debug_faces(vertices, faces, [size=], [disabled=]);
+// Description:
+//   Draws all the vertices at their 3D position, numbered in blue by their
+//   position in the vertex array.  Each face will have their face number drawn
+//   in red, aligned with the center of face.  All children of this module are drawn
+//   with transparency.
+// Arguments:
+//   vertices = Array of point vertices.
+//   faces = Array of faces by vertex numbers.
+//   ---
+//   size = The size of the text used to label the faces and vertices.  Default: 1
+//   disabled = If true, don't draw numbers, and draw children without transparency.  Default: false.
+// Example(EdgesMed):
+//   verts = [for (z=[-10,10], y=[-10,10], x=[-10,10]) [x,y,z]];
+//   faces = [[0,1,2], [1,3,2], [0,4,5], [0,5,1], [1,5,7], [1,7,3], [3,7,6], [3,6,2], [2,6,4], [2,4,0], [4,6,7], [4,7,5]];
+//   debug_faces(vertices=verts, faces=faces, size=2) {
+//       polyhedron(points=verts, faces=faces);
+//   }
+module debug_faces(vertices, faces, size=1, disabled=false) {
+    if (!disabled) {
+        vlen = len(vertices);
+        color("red") {
+            for (i = [0:1:len(faces)-1]) {
+                face = faces[i];
+                if (face[0] < 0 || face[1] < 0 || face[2] < 0 || face[0] >= vlen || face[1] >= vlen || face[2] >= vlen) {
+                    echo("BAD FACE: ", vlen=vlen, face=face);
+                } else {
+                    verts = select(vertices,face);
+                    c = mean(verts);
+                    v0 = verts[0];
+                    v1 = verts[1];
+                    v2 = verts[2];
+                    dv0 = unit(v1 - v0);
+                    dv1 = unit(v2 - v0);
+                    nrm0 = cross(dv0, dv1);
+                    nrm1 = UP;
+                    axis = vector_axis(nrm0, nrm1);
+                    ang = vector_angle(nrm0, nrm1);
+                    theta = atan2(nrm0[1], nrm0[0]);
+                    translate(c) {
+                        rotate(a=180-ang, v=axis) {
+                            zrot(theta-90)
+                            linear_extrude(height=size/10, center=true, convexity=10) {
+                                union() {
+                                    text(text=str(i), size=size, halign="center");
+                                    text(text=str("_"), size=size, halign="center");
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    debug_vertices(vertices, size=size, disabled=disabled) {
+        children();
+    }
+    if (!disabled) {
+        echo(faces=faces);
+    }
+}
+
+
+
+// Module: debug_vnf()
+// Usage:
+//   debug_vnf(vnfs, [convexity=], [txtsize=], [disabled=]);
+// Description:
+//   A drop-in module to replace `vnf_polyhedron()` and help debug vertices and faces.
+//   Draws all the vertices at their 3D position, numbered in blue by their
+//   position in the vertex array.  Each face will have its face number drawn
+//   in red, aligned with the center of face.  All given faces are drawn with
+//   transparency. All children of this module are drawn with transparency.
+//   Works best with Thrown-Together preview mode, to see reversed faces.
+// Arguments:
+//   vnf = vnf to display
+//   ---
+//   convexity = The max number of walls a ray can pass through the given polygon paths.
+//   txtsize = The size of the text used to label the faces and vertices.
+//   disabled = If true, act exactly like `polyhedron()`.  Default = false.
+// Example(EdgesMed):
+//   verts = [for (z=[-10,10], a=[0:120:359.9]) [10*cos(a),10*sin(a),z]];
+//   faces = [[0,1,2], [5,4,3], [0,3,4], [0,4,1], [1,4,5], [1,5,2], [2,5,3], [2,3,0]];
+//   debug_vnf([verts,faces], txtsize=2);
+module debug_vnf(vnf, convexity=6, txtsize=1, disabled=false) {
+    debug_faces(vertices=vnf[0], faces=vnf[1], size=txtsize, disabled=disabled) {
+        vnf_polyhedron(vnf, convexity=convexity);
+    }
+}
+
+
+
 // Function&Module: vnf_validate()
 // Usage: As Function
 //   fails = vnf_validate(vnf);