diff --git a/math.scad b/math.scad
index e125ba38..28b9b6cc 100644
--- a/math.scad
+++ b/math.scad
@@ -585,26 +585,35 @@ function _ceilall(data) =
 // Section: Constraints and Modulos
 
 // Function: constrain()
-// Synopsis: Returns a value constrained between `minval` and `maxval`, inclusive.
+// Summary: Limit (clamp) a number or array of numbers to a specified range of values.
 // Topics: Math
 // See Also: posmod(), modang()
 // Usage:
-//   val = constrain(v, minval, maxval);
+//   vals = constrain(v, minval, maxval);
 // Description:
-//   Constrains value to a range of values between minval and maxval, inclusive.
+//   Returns the value(s) in `v` limited to the range defined by `minval` and `maxval`.
+//   This operation is also known as "clamping" in other computer languages.
 // Arguments:
-//   v = value to constrain.
-//   minval = minimum value to return, if out of range.
-//   maxval = maximum value to return, if out of range.
+//   m = Value(s) to constrain. Can be a numerical value, a 1D vector, a 2D rectangular matrix, or a list of different-length vectors.
+//   minval = Minimum value to return. Set to `-INF` to unrestrict the minimum.
+//   maxval = Maximum value to return. Set to `INF` to unrestrict the maximum.
 // Example:
 //   a = constrain(-5, -1, 1);   // Returns: -1
 //   b = constrain(5, -1, 1);    // Returns: 1
 //   c = constrain(0.3, -1, 1);  // Returns: 0.3
 //   d = constrain(9.1, 0, 9);   // Returns: 9
-//   e = constrain(-0.1, 0, 9);  // Returns: 0
-function constrain(v, minval, maxval) = 
-    assert( is_finite(v+minval+maxval), "\nInput must be finite number(s).")
-    min(maxval, max(minval, v));
+//   e = constrain([1,2,3,4,5,6,7,8,9], 3, 7);          // Returns: [3,3,3,4,5,6,7,7,7]
+//   f = constrain([[1,2,3], [4,5,6], [7,8,9]], 3, 7);  // Returns: [[3,3,3], [4,5,6], [7,7,7]]
+//   g = constrain([[1,2,3,4], [5,6,7], [8,9]], 3, 7);  // Returns: [[3,3,3,4], [5,6,7], [7,7]]
+function constrain(v, minval, maxval) =
+    is_num(v) ? max(minval, min(v, maxval))
+    : is_vector(v) ? [for(f=v) max(minval, min(f, maxval))]
+    : is_matrix(v) ? let( // for a matrix, this should be more efficient than indexing
+        mflat = flatten(v),
+        clamped = [ for(f=mflat) max(minval, min(f, maxval)) ] 
+    ) list_to_matrix(clamped, len(v[0]), 0)
+    : is_list(v) ? [ for(vec=v) [ for(f=vec) max(minval, min(f, maxval)) ] ]
+    : assert(false, "\nIn constrain(), v must be a number, 1D vector, rectangular matrix, or list of vectors.");
 
 
 // Function: posmod()
diff --git a/vnf.scad b/vnf.scad
index 37d6f812..abc82f02 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -1575,17 +1575,23 @@ function vnf_bounds(vnf,fast=false) =
 // Topics: VNF Manipulation
 // See Also: vnf_halfspace()
 // Usage:
-//   region = projection(vnf, [cut]);
+//   region = projection(vnf, [cut], [z]);
 // Description:
-//   When `cut=false`, which is the default, projects the input VNF
-//   onto the XY plane, returning a region.  As currently implemented,  this operation
+//   Project a VNF object onto the xy plane at position `z`, returning a region.
+//   .
+//   The default action (`cut=false`) is to projects the input VNF
+//   onto the XY plane, returning a region.  As currently implemented, this operation
 //   involves the 2D union of all the projected faces and can be
 //   slow if the VNF has many faces.  Minimize the face count of the VNF for best performance. 
 //   .
 //   When `cut=true`, returns the intersection of the VNF with the
-//   XY plane, which is again a region.  If the VNF does not intersect
-//   the XY plane then returns the empty set.  This operation is
+//   XY plane at the position given by `z` (default `z=0`), which is again a region.
+//   If the VNF does not intersect the plane, then returns the empty set.  This operation is
 //   much faster than `cut=false`.
+// Arguments:
+//   vnf = The VNF object to project to a plane.
+//   cut = When true, returns a region containing intersection of the VNF with the plane. When false, projects the entire VNF onto the plane. Default: false
+//   z = Optional z position of the XY plane, useful when `cut=true` to get a specific slice position. Ignored if `cut=false`. Default: 0
 // Example(3D): Here's a VNF with two linked toruses and a small cube
 //   vnf = vnf_join([
 //            xrot(90,torus(id=15,od=24,$fn=5)),
@@ -1608,12 +1614,16 @@ function vnf_bounds(vnf,fast=false) =
 //   vnf = xrot(35,torus(id=4,od=12,$fn=32));
 //   reg = projection(vnf,cut=true);
 //   region(reg);
+// Example(2D): Projection of tilted torus using `cut=true` at a different z position for the XY plane.
+//   vnf = xrot(35,torus(id=4,od=12,$fn=32));
+//   reg = projection(vnf,cut=true,z=0.3);
+//   region(reg);
 
-function projection(vnf,cut=false,eps=EPSILON) =
+function projection(vnf,cut=false,z=0,eps=EPSILON) =
    assert(is_vnf(vnf))
    cut ?
          let(
-              vnf_bdy = vnf_halfspace([0,0,1,0],vnf, boundary=true),
+              vnf_bdy = vnf_halfspace([0,0,1,cut?z:0],vnf, boundary=true),
               ind = vnf_bdy[1],
               pts = path2d(vnf_bdy[0][0])
          )