diff --git a/arrays.scad b/arrays.scad
index 1e0d44e..5c0e7e7 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -1244,6 +1244,31 @@ function group_sort(list, idx) =
         
 
 
+// Function: list_smallest()
+// Usage:
+//   small = list_smallest(list, k)
+// Description:
+//   Returns a set of the k smallest items in list in arbitrary order.  The items must be
+//   mutually comparable with native OpenSCAD comparison operations.  You will get "undefined operation"
+//   errors if you provide invalid input. 
+// Arguments:
+//   list = list to process
+//   k = number of items to return
+function list_smallest(list, k) =
+    assert(is_list(list))
+    assert(is_finite(k) && k>=0, "k must be nonnegative")
+    let( 
+        v       = list[rand_int(0,len(list)-1,1)[0]],
+        smaller = [for(li=list) if(li<v) li ],
+        equal   = [for(li=list) if(li==v) li ]
+    )
+    len(smaller)   == k ? smaller :
+    len(smaller)<k && len(smaller)+len(equal) >= k ? [ each smaller, for(i=[1:k-len(smaller)]) v ] :
+    len(smaller)   >  k ? list_smallest(smaller, k) :
+    let( bigger  = [for(li=list) if(li>v) li ] )
+    concat(smaller, equal, list_smallest(bigger, k-len(smaller) -len(equal)));
+
+
 // Function: group_data()
 // Usage:
 //   groupings = group_data(groups, values);
diff --git a/attachments.scad b/attachments.scad
index 7b1f510..efa11ad 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -1117,7 +1117,7 @@ function reorient(
     two_d=false,
     axis=UP,
     p=undef
-) =
+) = 
     assert(is_undef(anchor) || is_vector(anchor) || is_string(anchor), str("Got: ",anchor))
     assert(is_undef(spin)   || is_vector(spin,3) || is_num(spin), str("Got: ",spin))
     assert(is_undef(orient) || is_vector(orient,3), str("Got: ",orient))
@@ -1269,7 +1269,7 @@ function _attach_geom(
     axis=UP
 ) =
     assert(is_bool(extent))
-    assert(is_vector(cp))
+    assert(is_vector(cp) || is_string(cp))
     assert(is_vector(offset))
     assert(is_list(anchors))
     assert(is_bool(two_d))
@@ -1496,6 +1496,21 @@ function _attach_transform(anchor, spin, orient, geom, p) =
     apply(m, p);
 
 
+function _get_cp(geom) =
+    let(cp=select(geom,-3))
+    is_vector(cp) ? cp
+  : let(
+        type = in_list(geom[0],["vnf_extent","vnf_isect"]) ? "vnf"
+             : in_list(geom[0],["path_extent","path_isect"]) ? "path"
+             : "other"
+    )
+    assert(type!="other", "Invalid cp value")
+    cp=="centroid" ? centroid(geom[1])
+  : let(points = type=="vnf"?geom[1][0]:geom[1])
+    cp=="mean" ? mean(points)
+  : cp=="box" ? mean(pointlist_bounds(points))
+  : assert(false,"Invalid cp specification");
+
 
 /// Internal Function: _find_anchor()
 // Usage:
@@ -1511,19 +1526,19 @@ function _attach_transform(anchor, spin, orient, geom, p) =
 //   anchor = Vector or named anchor string.
 //   geom = The geometry description of the shape.
 function _find_anchor(anchor, geom) =
-    let(
-        cp = select(geom,-3),
+    let( 
+        cp = _get_cp(geom),
         offset_raw = select(geom,-2),
         offset = [for (i=[0:2]) anchor[i]==0? 0 : offset_raw[i]],  // prevents bad centering.
         anchors = last(geom),
         type = geom[0]
     )
-    is_string(anchor)? (
-        anchor=="origin"? [anchor, CENTER, UP, 0] :
-        let(found = search([anchor], anchors, num_returns_per_match=1)[0])
-        assert(found!=[], str("Unknown anchor: ",anchor))
-        anchors[found]
-    ) :
+    is_string(anchor)? (  
+          anchor=="origin"? [anchor, CENTER, UP, 0]
+        : let(found = search([anchor], anchors, num_returns_per_match=1)[0])
+          assert(found!=[], str("Unknown anchor: ",anchor))
+          anchors[found]
+    ) : 
     assert(is_vector(anchor),str("anchor=",anchor))
     let(anchor = point3d(anchor))
     anchor==CENTER? [anchor, cp, UP, 0] :
@@ -1590,27 +1605,25 @@ function _find_anchor(anchor, geom) =
             eps = 1/2048,
             points = vnf[0],
             faces = vnf[1],
-            rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), points),
+            rpts = apply(rot(from=anchor, to=RIGHT) * move(-cp), points),
             hits = [
-                for (face = faces) let(
-                    verts = select(rpts, face),
-                    xs = columns(verts,0),
-                    ys = columns(verts,1),
-                    zs = columns(verts,2)
-                ) if (
-                    max(xs) >= -eps &&
-                    max(ys) >= -eps &&
-                    max(zs) >= -eps &&
-                    min(ys) <=  eps &&
-                    min(zs) <=  eps
-                ) let(
-                    poly = select(points, face),
-                    pt = polygon_line_intersection(poly, [cp,cp+anchor], bounded=[true,false], eps=eps)
-                ) if (!is_undef(pt)) let(
-                    plane = plane_from_polygon(poly),
-                    n = unit(plane_normal(plane))
-                )
-                [norm(pt-cp), n, pt]
+                for (face = faces)
+                    let(
+                        verts = select(rpts, face),
+                        ys = columns(verts,1),
+                        zs = columns(verts,2)
+                    )
+                    if (max(ys) >= -eps && max(zs) >= -eps &&
+                        min(ys) <=  eps &&  min(zs) <=  eps)
+                        let(
+                            poly = select(points, face),
+                            isect = polygon_line_intersection(poly, [cp,cp+anchor], eps=eps),
+                            ptlist = is_undef(isect) ? [] :
+                                     is_vector(isect) ? [isect]
+                                                      : flatten(isect),   // parallel to a face
+                            n = len(ptlist)>0 ? polygon_normal(poly) : undef
+                        )
+                        for(pt=ptlist) [anchor * (pt-cp), n, pt]
             ]
         )
         assert(len(hits)>0, "Anchor vector does not intersect with the shape.  Attachment failed.")
@@ -1619,17 +1632,17 @@ function _find_anchor(anchor, geom) =
             dist = hits[furthest][0],
             pos = hits[furthest][2],
             hitnorms = [for (hit = hits) if (approx(hit[0],dist,eps=eps)) hit[1]],
-            unorms = len(hitnorms) > 7
-              ? unique([for (nn = hitnorms) quant(nn,1e-9)])
-              : [
-                    for (i = idx(hitnorms)) let(
-                        nn = hitnorms[i],
-                        isdup = [
-                            for (j = [i+1:1:len(hitnorms)-1])
-                            if (approx(nn, hitnorms[j])) 1
-                        ] != []
-                    ) if (!isdup) nn
-                ],
+            unorms = [
+                      for (i = idx(hitnorms))
+                          let(
+                              thisnorm = hitnorms[i],
+                              isdup = [
+                                       for (j = [i+1:1:len(hitnorms)-1])
+                                           if (approx(thisnorm, hitnorms[j])) 1
+                                      ] != []
+                          )
+                          if (!isdup) thisnorm
+                     ],
             n = unit(sum(unorms)),
             oang = approx(point2d(n), [0,0])? 0 : atan2(n.y, n.x) + 90
         )
diff --git a/geometry.scad b/geometry.scad
index 112d8d3..21ee616 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -546,7 +546,7 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
 //   xyzpath = rot(45, v=[0,1,0], p=path3d(star(n=5,step=2,d=100), 70));
 //   plane = plane_from_polygon(xyzpath);
 //   #stroke(xyzpath,closed=true,width=3);
-//   cp = polygon_centroid(xyzpath);
+//   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(45);
 function plane_from_polygon(poly, fast=false, eps=EPSILON) =
     assert( is_path(poly,dim=3), "Invalid polygon." )
@@ -897,7 +897,7 @@ function plane_line_angle(plane, line) =
 //   proj = plane_closest_point(plane,points);
 //   color("red") move_copies(points) sphere(d=4,$fn=12);
 //   color("blue") move_copies(proj) sphere(d=4,$fn=12);
-//   move(polygon_centroid(proj)) {
+//   move(centroid(proj)) {
 //       rot(from=UP,to=plane_normal(plane)) {
 //           anchor_arrow(50);
 //           %cube([120,150,0.1],center=true);
@@ -1403,21 +1403,60 @@ function polygon_area(poly, signed=false) =
         signed ? total : abs(total);
 
 
-// Function: polygon_centroid()
+// Function: centroid()
 // Usage:
-//   cpt = polygon_centroid(poly);
+//   c = centroid(object, [eps]);
 // Topics: Geometry, Polygons, Centroid
 // Description:
 //   Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
 //   Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
-//   Collinear points produce an error.  The results are meaningless for self-intersecting
-//   polygons or an error is produced.
+//   If you provide a non-planar or collinear polygon you will get an error.  For self-intersecting
+//   polygons you may get an error or you may get meaningless results.
+//   .
+//   If object is a manifold VNF then returns the 3d centroid of the polyhedron.  The VNF must
+//   describe a valid polyhedron with consistent face direction and no holes in the mesh; otherwise
+//   the results are undefined.
 // Arguments:
-//   poly = Points of the polygon from which the centroid is calculated.
-//   eps = Tolerance in geometric comparisons.  Default: `EPSILON` (1e-9)
-function polygon_centroid(poly, eps=EPSILON) =
+//   object = object to compute the centroid of
+//   eps = epsilon value for identifying degenerate cases
+function centroid(object,eps=EPSILON) =
+    assert(is_finite(eps) && (eps>=0), "The tolerance should a non-negative value." )
+    is_vnf(object) ? _vnf_centroid(object,eps)
+  : is_path(object,[2,3]) ? _polygon_centroid(object,eps)
+  : is_region(object) ? (len(object)==1 ? _polygon_centroid(object[0],eps) : _region_centroid(object,eps))
+  : assert(false, "Input must be a VNF, a region, or a 2D or 3D polygon");
+
+
+/// Internal Function: _region_centroid()
+/// Compute centroid of region
+function _region_centroid(region,eps=EPSILON) =
+   let(
+       region=force_region(region),
+       parts = region_parts(region),
+       // Rely on region_parts returning all outside polygons clockwise
+       // and inside (hole) polygons counterclockwise, so areas have reversed sign
+       cent_area = [for(R=parts, p=R)
+                       let(A=polygon_area(p,signed=true))
+                       [A*_polygon_centroid(p),A]],
+       total = sum(cent_area)
+   )
+   total[0]/total[1];
+
+
+/// Function: _polygon_centroid()
+/// Usage:
+///   cpt = _polygon_centroid(poly);
+/// Topics: Geometry, Polygons, Centroid
+/// Description:
+///   Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
+///   Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
+///   Collinear points produce an error.  The results are meaningless for self-intersecting
+///   polygons or an error is produced.
+/// Arguments:
+///   poly = Points of the polygon from which the centroid is calculated.
+///   eps = Tolerance in geometric comparisons.  Default: `EPSILON` (1e-9)
+function _polygon_centroid(poly, eps=EPSILON) =
     assert( is_path(poly,dim=[2,3]), "The input must be a 2D or 3D polygon." )
-    assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
     let(
         n = len(poly[0])==2 ? 1 :
             let( plane = plane_from_points(poly, fast=false))
@@ -1633,21 +1672,21 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
 //   color("lightblue") for(tri=tris) polygon(select(poly,tri));
 //   color("blue")    up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
 //   color("magenta") up(2) stroke(poly,.25,closed=true); 
-//   color("black")   up(3) vnf_debug([poly,[]],faces=false,size=1);
+//   color("black")   up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
 // Example(2D,NoAxes): a polygon with a hole and one "contact" edge
 //   poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ];
 //   tris =  polygon_triangulate(poly);
 //   color("lightblue") for(tri=tris) polygon(select(poly,tri));
 //   color("blue")    up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
 //   color("magenta") up(2) stroke(poly,.25,closed=true); 
-//   color("black")   up(3) vnf_debug([poly,[]],faces=false,size=1);
+//   color("black")   up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
 // Example(2D,NoAxes): a polygon with "touching" vertices and no holes
 //   poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ];
 //   tris =  polygon_triangulate(poly);
 //   color("lightblue") for(tri=tris) polygon(select(poly,tri));
 //   color("blue")    up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
 //   color("magenta") up(2) stroke(poly,.25,closed=true); 
-//   color("black")   up(3) vnf_debug([poly,[]],faces=false,size=1);
+//   color("black")   up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
 // Example(2D,NoAxes): a polygon with "contact" edges and no holes
 //   poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3], 
 //            [7,7], [7,3], [3,3] ];
@@ -1655,7 +1694,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
 //   color("lightblue") for(tri=tris) polygon(select(poly,tri));
 //   color("blue")    up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
 //   color("magenta") up(2) stroke(poly,.25,closed=true); 
-//   color("black")   up(3) vnf_debug([poly,[]],faces=false,size=1);
+//   color("black")   up(3) vnf_debug([path3d(poly),[]],faces=false,size=1);
 // Example(3D): 
 //   include <BOSL2/polyhedra.scad>
 //   vnf = regular_polyhedron_info(name="dodecahedron",side=5,info="vnf");
diff --git a/hull.scad b/hull.scad
index 34a57d4..a2efb7d 100644
--- a/hull.scad
+++ b/hull.scad
@@ -93,8 +93,8 @@ function _is_cw(a,b,c,all) =
 //   Returns a path as a list of indices into `points`. 
 //   When all==true, returns extra points that are on edges of the hull.
 // Arguments:
-//   points - list of 2d points to get the hull of.
-//   all - when true, includes all points on the edges of the convex hull. Default: false.
+//   points = list of 2d points to get the hull of.
+//   all = when true, includes all points on the edges of the convex hull. Default: false.
 // Example(2D):
 //   pts = [[-10,-10], [0,10], [10,10], [12,-10]];
 //   path = hull2d_path(pts);
diff --git a/math.scad b/math.scad
index 6ac428d..a80aa82 100644
--- a/math.scad
+++ b/math.scad
@@ -827,27 +827,21 @@ function mean(v) =
     sum(v)/len(v);
 
 
-// Function: ninther()
-// Usage:
-//    med = ninther(v)
-// Description:
-//    Finds a value in the input list of numbers `v` that is the median of  a 
-//    sample of 9 entries of `v`.
-//    It is a much faster approximation of the true median computation.
-// Arguments:
-//    v = an array of numbers
-function ninther(v) = 
-    let( l=len(v) )
-    l<=4 ? l<=2 ? v[0] : _med3(v[0], v[1], v[2]) : 
-    l==5 ? _med3(v[0], _med3(v[1], v[2], v[3]), v[4]) :
-    _med3(_med3(v[0],v[floor(l/6)],v[floor(l/3)]),
-          _med3(v[floor(l/3)],v[floor(l/2)],v[floor(2*l/3)]),
-          _med3(v[floor(2*l/3)],v[floor((5*l/3 -1)/2)],v[l-1]) );
 
-// the median of a triple
-function _med3(a,b,c) =
-    a < c ? a < b ? min(b,c) : min(a,c) :
-    b < c ? min(a,c) : min(a,b);
+// Function: median()
+// Usage:
+//   middle = median(v)
+// Description:
+//   Returns the median of the given vector.  
+function median(v) =
+    assert(is_vector(v), "Input to median must be a vector")
+    len(v)%2 ? max( list_smallest(v, ceil(len(v)/2)) ) :
+    let( lowest = list_smallest(v, len(v)/2 + 1),
+         max  = max(lowest),
+         imax = search(max,lowest,1),
+         max2 = max([for(i=idx(lowest)) if(i!=imax[0]) lowest[i] ])
+    )
+    (max+max2)/2;
 
 
 // Function: convolve()
@@ -953,7 +947,7 @@ function linear_solve(A,b,pivot=true) =
 // Description:
 //    Compute the matrix inverse of the square matrix `A`.  If `A` is singular, returns `undef`.
 //    Note that if you just want to solve a linear system of equations you should NOT use this function.
-//    Instead use [[`linear_solve()`|linear_solve]], or use [[`qr_factor()`|qr_factor]].  The computation
+//    Instead use {{linear_solve()}}, or use {{qr_factor()}}.  The computation
 //    will be faster and more accurate.  
 function matrix_inverse(A) =
     assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
@@ -1002,7 +996,7 @@ function null_space(A,eps=1e-12) =
 //   qr = qr_factor(A,[pivot]);
 // Description:
 //   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P].  This factorization can be
-//   used to solve linear systems of equations.  The factorization is A = Q*R*transpose(P).  If pivot is false (the default)
+//   used to solve linear systems of equations.  The factorization is `A = Q*R*transpose(P)`.  If pivot is false (the default)
 //   then P is the identity matrix and A = Q*R.  If pivot is true then column pivoting results in an R matrix where the diagonal
 //   is non-decreasing.  The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
 //   for rank estimation or computation of the null space, but it may be slower.  
@@ -1083,7 +1077,7 @@ function _back_substitute(R, b, x=[]) =
 //   L = cholesky(A);
 // Description:
 //   Compute the cholesky factor, L, of the symmetric positive definite matrix A.
-//   The matrix L is lower triangular and L * transpose(L) = A.  If the A is
+//   The matrix L is lower triangular and `L * transpose(L) = A`.  If the A is
 //   not symmetric then an error is displayed.  If the matrix is symmetric but
 //   not positive definite then undef is returned.  
 function cholesky(A) =
diff --git a/regions.scad b/regions.scad
index 1a4e690..2084efe 100644
--- a/regions.scad
+++ b/regions.scad
@@ -390,20 +390,25 @@ function region_parts(region) =
 //   If called as a module, creates a polyhedron that is the linear extrusion of the given 2D region or path.
 //   If called as a function, returns a VNF that can be used to generate a polyhedron of the linear extrusion
 //   of the given 2D region or path.  The benefit of using this, over using `linear_extrude region(rgn)` is
-//   that you can use `anchor`, `spin`, `orient` and attachments with it.  Also, you can make more refined
+//   that it supports `anchor`, `spin`, `orient` and attachments.  You can also make more refined
 //   twisted extrusions by using `maxseg` to subsample flat faces.
+//   Note that the center option centers vertically using the named anchor "zcenter" whereas
+//   `anchor=CENTER` centers the entire shape relative to
+//   the shape's centroid, or other centerpoint you specify.  The centerpoint can be "centroid", "mean", "box" or
+//   a custom point location.  
 // Arguments:
 //   region = The 2D [Region](regions.scad) or path that is to be extruded.
 //   height = The height to extrude the region.  Default: 1
-//   center = If true, the created polyhedron will be vertically centered.  If false, it will be extruded upwards from the origin.  Default: `false`
+//   center = If true, the created polyhedron will be vertically centered.  If false, it will be extruded upwards from the XY plane.  Default: `false`
 //   slices = The number of slices to divide the shape into along the Z axis, to allow refinement of detail, especially when working with a twist.  Default: `twist/5`
 //   maxseg = If given, then any long segments of the region will be subdivided to be shorter than this length.  This can refine twisting flat faces a lot.  Default: `undef` (no subsampling)
 //   twist = The number of degrees to rotate the shape clockwise around the Z axis, as it rises from bottom to top.  Default: 0
 //   scale = The amount to scale the shape, from bottom to top.  Default: 1
 //   style = The style to use when triangulating the surface of the object.  Valid values are `"default"`, `"alt"`, or `"quincunx"`.
 //   convexity = Max number of surfaces any single ray could pass through.  Module use only.
+//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `"origin"`
 //   anchor_isect = If true, anchoring it performed by finding where the anchor vector intersects the swept shape.  Default: false
-//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
+//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 3D point.  Default: "centroid"
 //   spin = Rotate this many degrees around the Z axis after anchor.  See [spin](attachments.scad#spin).  Default: `0`
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
 // Example: Extruding a Compound Region.
@@ -427,18 +432,18 @@ function region_parts(region) =
 //   mrgn = union(rgn1,rgn2);
 //   orgn = difference(mrgn,rgn3);
 //   linear_sweep(orgn,height=20,convexity=16) show_anchors();
-module linear_sweep(region, height=1, center, twist=0, scale=1, slices, maxseg, style="default", convexity, anchor_isect=false, anchor, spin=0, orient=UP) {
+module linear_sweep(region, height=1, center, twist=0, scale=1, slices, maxseg, style="default", convexity, anchor_isect=false, anchor, spin=0, orient=UP, cp="centroid", anchor="origin") {
     region = force_region(region);
     dummy=assert(is_region(region),"Input is not a region");
-    cp = mean(pointlist_bounds(flatten(region)));
-    anchor = get_anchor(anchor, center, "origin", "origin");
+    anchor = center ? "zcenter" : anchor;
+    anchors = [named_anchor("zcenter", [0,0,height/2], UP)];
     vnf = linear_sweep(
         region, height=height,
         twist=twist, scale=scale,
         slices=slices, maxseg=maxseg,
         style=style
     );
-    attachable(anchor,spin,orient, cp=cp, vnf=vnf, extent=!anchor_isect) {
+    attachable(anchor,spin,orient, cp=cp, vnf=vnf, extent=!anchor_isect, anchors=anchors) {
         vnf_polyhedron(vnf, convexity=convexity);
         children();
     }
@@ -446,15 +451,15 @@ module linear_sweep(region, height=1, center, twist=0, scale=1, slices, maxseg,
 
 
 function linear_sweep(region, height=1, center, twist=0, scale=1, slices,
-                      maxseg, style="default", anchor_isect=false, anchor, spin=0, orient=UP) =
+                      maxseg, style="default", cp="centroid", anchor_isect=false, anchor, spin=0, orient=UP) =
     let(
         region = force_region(region)
     )
     assert(is_region(region), "Input is not a region")
     let(
-        anchor = get_anchor(anchor,center,BOT,BOT),
+        anchor = center ? "zcenter" : anchor,
+        anchors = [named_anchor("zcenter", [0,0,height/2], UP)],
         regions = region_parts(region),
-        cp = mean(pointlist_bounds(flatten(region))),
         slices = default(slices, floor(twist/5+1)),
         step = twist/slices,
         hstep = height/slices,
@@ -484,14 +489,14 @@ function linear_sweep(region, height=1, center, twist=0, scale=1, slices,
                     for (i=[0:1:slices]) let(
                         sc = lerp(1, scale, i/slices),
                         ang = i * step,
-                        h = i * hstep - height/2
+                        h = i * hstep //- height/2
                     ) scale([sc,sc,1], p=rot(ang, p=path3d(path,h)))
                 ]
             ) vnf_vertex_array(verts, caps=false, col_wrap=true, style=style),
-            for (rgn = regions) vnf_from_region(rgn, down(height/2), reverse=true),
-            for (rgn = trgns) vnf_from_region(rgn, up(height/2), reverse=false)
+            for (rgn = regions) vnf_from_region(rgn, ident(4), reverse=true),
+            for (rgn = trgns) vnf_from_region(rgn, up(height), reverse=false)
         ])
-    ) reorient(anchor,spin,orient, cp=cp, vnf=vnf, extent=!anchor_isect, p=vnf);
+    ) reorient(anchor,spin,orient, cp=cp, vnf=vnf, extent=!anchor_isect, p=vnf, anchors=anchors);
 
 
 
diff --git a/rounding.scad b/rounding.scad
index 19a6c28..69e2bf1 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -1041,7 +1041,7 @@ module offset_sweep(path, height,
                        quality=quality, check_valid=true, extra=extra, cut=cut, chamfer_width=chamfer_width,
                        chamfer_height=chamfer_height, joint=joint, k=k, angle=angle);
   
-    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
     {
         vnf_polyhedron(vnf,convexity=convexity);
         children();
@@ -1818,7 +1818,7 @@ module rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_bot
   result = rounded_prism(bottom=bottom, top=top, joint_bot=joint_bot, joint_top=joint_top, joint_sides=joint_sides,
                          k_bot=k_bot, k_top=k_top, k_sides=k_sides, k=k, splinesteps=splinesteps, h=h, length=length, height=height, l=l,debug=debug);
   vnf = debug ? result[1] : result;
-  attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+  attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
   {
     if (debug){
         vnf_polyhedron(vnf, convexity=convexity);
diff --git a/skin.scad b/skin.scad
index b7080a6..e019846 100644
--- a/skin.scad
+++ b/skin.scad
@@ -278,12 +278,12 @@
 // Example(FlatSpin,VPD=32,VPT=[1.2,4.3,2]): Another "distance" example:
 //   off = [0,2]; 
 //   shape = turtle(["right",45,"move", "left",45,"move", "left",45, "move", "jump", [.5+sqrt(2)/2,8]]);
-//   rshape = rot(180,cp=polygon_centroid(shape)+off, p=shape);
+//   rshape = rot(180,cp=centroid(shape)+off, p=shape);
 //   skin([shape,rshape],z=[0,4], method="distance",slices=10,refine=15);
 // Example(FlatSpin,VPD=32,VPT=[1.2,4.3,2]): Slightly shifting the profile changes the optimal linkage
 //   off = [0,1]; 
 //   shape = turtle(["right",45,"move", "left",45,"move", "left",45, "move", "jump", [.5+sqrt(2)/2,8]]);
-//   rshape = rot(180,cp=polygon_centroid(shape)+off, p=shape);
+//   rshape = rot(180,cp=centroid(shape)+off, p=shape);
 //   skin([shape,rshape],z=[0,4], method="distance",slices=10,refine=15);
 // Example(FlatSpin,VPD=444,VPT=[0,0,50]): This optimal solution doesn't look terrible:
 //   prof1 = path3d([[-50,-50], [-50,50], [50,50], [25,25], [50,0], [25,-25], [50,-50]]);
@@ -386,7 +386,7 @@ module skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=
             anchor="origin",cp,spin=0, orient=UP, extent=false)
 {   
     vnf = skin(profiles, slices, refine, method, sampling, caps, closed, z, style=style);
-    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
     {      
         vnf_polyhedron(vnf,convexity=convexity);
         children();
@@ -816,7 +816,7 @@ module path_sweep(shape, path, method="incremental", normal, closed=false, twist
 {
     vnf = path_sweep(shape, path, method, normal, closed, twist, twist_by_length,
                     symmetry, last_normal, tangent, relaxed, caps, style);
-    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
     {      
         vnf_polyhedron(vnf,convexity=convexity);
         children();
@@ -1001,7 +1001,7 @@ module path_sweep2d(profile, path, closed=false, caps, quality=1, style="min_edg
                     anchor="origin", cp, spin=0, orient=UP, extent=false)
 {
    vnf = path_sweep2d(profile, path, closed, caps, quality, style);
-   attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+   attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
     {      
         vnf_polyhedron(vnf,convexity=convexity);
         children();
@@ -1128,7 +1128,7 @@ module sweep(shape, transforms, closed=false, caps, style="min_edge", convexity=
              anchor="origin",cp,spin=0, orient=UP, extent=false)
 {
     vnf = sweep(shape, transforms, closed, caps, style);
-    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf))
+    attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : centroid(vnf))
     {      
         vnf_polyhedron(vnf,convexity=convexity);
         children();
@@ -1591,7 +1591,7 @@ function _skin_tangent_match(poly1, poly2) =
         swap = len(poly1)>len(poly2),
         big = swap ? poly1 : poly2,
         small = swap ? poly2 : poly1,
-        curve_offset = polygon_centroid(small)-polygon_centroid(big),
+        curve_offset = centroid(small)-centroid(big),
         cutpts = [for(i=[0:len(small)-1]) _find_one_tangent(big, select(small,i,i+1),curve_offset=curve_offset)],
         shift = last(cutpts)+1,
         newbig = polygon_shift(big, shift),
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index d7405df..7eff2bc 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -46,7 +46,7 @@ test_is_polygon_convex();
 test_polygon_shift();
 test_reindex_polygon();
 test_align_polygon();
-test_polygon_centroid();
+test_centroid();
 test_point_in_polygon();
 test_polygon_triangulate();
 test_is_polygon_clockwise();
@@ -835,15 +835,31 @@ module test_noncollinear_triple() {
 *test_noncollinear_triple();
 
 
-module test_polygon_centroid() {
+module test_centroid() {
+    // polygons
     $fn = 24;
-    assert_approx(polygon_centroid(circle(d=100)), [0,0]);
-    assert_approx(polygon_centroid(rect([40,60],rounding=10,anchor=LEFT)), [20,0]);
-    assert_approx(polygon_centroid(rect([40,60],rounding=10,anchor=FWD)), [0,30]);
+    assert_approx(centroid(circle(d=100)), [0,0]);
+    assert_approx(centroid(rect([40,60],rounding=10,anchor=LEFT)), [20,0]);
+    assert_approx(centroid(rect([40,60],rounding=10,anchor=FWD)), [0,30]);
     poly = move([1,2.5,3.1],p=rot([12,49,24], p=path3d(circle(10,$fn=33))));
-    assert_approx(polygon_centroid(poly), [1,2.5,3.1]);
+    assert_approx(centroid(poly), [1,2.5,3.1]);
+
+    // regions
+    R = [square(10), move([5,4],circle(r=3,$fn=32)),  right(15,square(7)), move([18,3],circle(r=2,$fn=5))];
+    assert_approx(centroid(R), [9.82836532809, 4.76313546433]);
+
+    // VNFs
+    assert_approx(centroid(cube(100, center=false)), [50,50,50]);
+    assert_approx(centroid(cube(100, center=true)), [0,0,0]);
+    assert_approx(centroid(cube(100, anchor=ALLPOS)), [-50,-50,-50]);
+    assert_approx(centroid(cube(100, anchor=BOT)), [0,0,50]);
+    assert_approx(centroid(cube(100, anchor=TOP)), [0,0,-50]);
+    assert_approx(centroid(sphere(d=100, anchor=CENTER, $fn=36)), [0,0,0]);
+    assert_approx(centroid(sphere(d=100, anchor=BOT, $fn=36)), [0,0,50]);
+    ellipse = xscale(2, p=circle($fn=24, r=3));
+    assert_approx(centroid(path_sweep(pentagon(r=1), path3d(ellipse), closed=true)),[0,0,0]);
 }
-*test_polygon_centroid();
+*test_centroid();
 
 
 
diff --git a/tests/test_regions.scad b/tests/test_regions.scad
index 90ada2e..9a2126b 100644
--- a/tests/test_regions.scad
+++ b/tests/test_regions.scad
@@ -78,7 +78,7 @@ module test_exclusive_or() {
   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);
+  p2 = move(-centroid(p),p);
   p3 = polygon_parts(p2);
   p4 = exclusive_or(p3,square(51,center=true));
 
diff --git a/tests/test_vnf.scad b/tests/test_vnf.scad
index 3c65c26..14675fb 100644
--- a/tests/test_vnf.scad
+++ b/tests/test_vnf.scad
@@ -43,18 +43,6 @@ module test_vnf_from_polygons() {
 test_vnf_from_polygons();
 
 
-module test_vnf_centroid() {
-    assert_approx(vnf_centroid(cube(100, center=false)), [50,50,50]);
-    assert_approx(vnf_centroid(cube(100, center=true)), [0,0,0]);
-    assert_approx(vnf_centroid(cube(100, anchor=ALLPOS)), [-50,-50,-50]);
-    assert_approx(vnf_centroid(cube(100, anchor=BOT)), [0,0,50]);
-    assert_approx(vnf_centroid(cube(100, anchor=TOP)), [0,0,-50]);
-    assert_approx(vnf_centroid(sphere(d=100, anchor=CENTER, $fn=36)), [0,0,0]);
-    assert_approx(vnf_centroid(sphere(d=100, anchor=BOT, $fn=36)), [0,0,50]);
-    ellipse = xscale(2, p=circle($fn=24, r=3));
-    assert_approx(vnf_centroid(path_sweep(pentagon(r=1), path3d(ellipse), closed=true)),[0,0,0]);}
-test_vnf_centroid();
-
 
 module test_vnf_volume() {
     assert_approx(vnf_volume(cube(100, center=false)), 1000000);
diff --git a/vectors.scad b/vectors.scad
index 4abdfe1..c57ebdd 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -494,7 +494,7 @@ function _bt_tree(points, ind, leafsize=25) =
         pmc    = mean(projc), 
         pivot  = min_index([for(p=projc) abs(p-pmc)]),
         radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]),
-        median = ninther(projc),
+        median = median(projc),
         Lind   = [for(i=idx(ind)) if(projc[i]<=median && i!=pivot) ind[i] ],
         Rind   = [for(i=idx(ind)) if(projc[i] >median && i!=pivot) ind[i] ]
       )
diff --git a/vnf.scad b/vnf.scad
index f6a27c0..3d03ae1 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -462,7 +462,7 @@ function is_vnf(x) =
     len(x)==2 &&
     is_list(x[0]) &&
     is_list(x[1]) &&
-    (x[0]==[] || (len(x[0])>=3 && is_vector(x[0][0]))) &&
+    (x[0]==[] || (len(x[0])>=3 && is_vector(x[0][0],3))) &&
     (x[1]==[] || is_vector(x[1][0]));
 
 
@@ -684,7 +684,7 @@ function _slice_3dpolygons(polys, dir, cuts) =
 //   orient = Vector to rotate top towards, after spin.  See [orient](attachments.scad#orient).  Default: `UP`
 module vnf_polyhedron(vnf, convexity=2, extent=true, cp=[0,0,0], anchor="origin", spin=0, orient=UP) {
     vnf = is_vnf_list(vnf)? vnf_merge(vnf) : vnf;
-    cp = is_def(cp) ? cp : vnf_centroid(vnf);
+    cp = is_def(cp) ? cp : centroid(vnf);
     attachable(anchor,spin,orient, vnf=vnf, extent=extent, cp=cp) {
         polyhedron(vnf[0], vnf[1], convexity=convexity);
         children();
@@ -760,17 +760,17 @@ function vnf_area(vnf) =
     sum([for(face=vnf[1]) polygon_area(select(verts,face))]);
 
 
-// Function: vnf_centroid()
-// Usage:
-//   vol = vnf_centroid(vnf);
-// Description:
-//   Returns the centroid of the given manifold VNF.  The VNF must describe a valid polyhedron with consistent face direction and
-//   no holes; otherwise the results are undefined.
+/// Function: _vnf_centroid()
+/// Usage:
+///   vol = _vnf_centroid(vnf);
+/// Description:
+///   Returns the centroid of the given manifold VNF.  The VNF must describe a valid polyhedron with consistent face direction and
+///   no holes; otherwise the results are undefined.
 
-// Divide the solid up into tetrahedra with the origin as one vertex.  
-// The centroid of a tetrahedron is the average of its vertices.
-// The centroid of the total is the volume weighted average.
-function vnf_centroid(vnf) =
+/// Divide the solid up into tetrahedra with the origin as one vertex.  
+/// The centroid of a tetrahedron is the average of its vertices.
+/// The centroid of the total is the volume weighted average.
+function _vnf_centroid(vnf,eps=EPSILON) =
     assert(is_vnf(vnf) && len(vnf[0])!=0 ) 
     let(
         verts = vnf[0],
@@ -784,7 +784,7 @@ function vnf_centroid(vnf) =
             [ vol, (v0+v1+v2)*vol ]
         ])
     )
-    assert(!approx(pos[0],0, EPSILON), "The vnf has self-intersections.")
+    assert(!approx(pos[0],0, eps), "The vnf has self-intersections.")
     pos[1]/pos[0]/4;