diff --git a/math.scad b/math.scad
index f8dc707..5ec3209 100644
--- a/math.scad
+++ b/math.scad
@@ -438,8 +438,11 @@ function lcm(a,b=[]) =
 //   sum([1,2,3]);  // returns 6.
 //   sum([[1,2,3], [3,4,5], [5,6,7]]);  // returns [9, 12, 15]
 function sum(v, dflt=0) =
+    is_vector(v) ? [for(i=v) 1]*v :
     assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
-    len(v) == 0 ? dflt : _sum(v,v[0]*0);
+    is_vector(v[0]) ? [for(i=v) 1]*v :
+    len(v) == 0 ? dflt :
+                  _sum(v,v[0]*0);
 
 function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
 
@@ -896,7 +899,13 @@ function count_true(l, nmax=undef, i=0, cnt=0) =
 //   data[len(data)-1].  This function uses a symetric derivative approximation
 //   for internal points, f'(t) = (f(t+h)-f(t-h))/2h.  For the endpoints (when closed=false) the algorithm
 //   uses a two point method if sufficient points are available: f'(t) = (3*(f(t+h)-f(t)) - (f(t+2*h)-f(t+h)))/2h.
+// 
+//   If `h` is a vector then it is assumed to be nonuniform, with h[i] giving the sampling distance
+//   between data[i+1] and data[i], and the data values will be linearly resampled at each corner
+//   to produce a uniform spacing for the derivative estimate.  At the endpoints a single point method
+//   is used: f'(t) = (f(t+h)-f(t))/h.  
 function deriv(data, h=1, closed=false) =
+    is_vector(h) ? _deriv_nonuniform(data, h, closed=closed) :
     let( L = len(data) )
     closed? [
         for(i=[0:1:L-1])
@@ -916,6 +925,28 @@ function deriv(data, h=1, closed=false) =
     ];
 
 
+function _dnu_calc(f1,fc,f2,h1,h2) =
+    let(
+        f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 , 
+        f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
+    )
+    (f2-f1) / 2 / min([h1,h2]);
+
+
+function _deriv_nonuniform(data, h, closed) =
+    assert(len(h) == len(data)-(closed?0:1),str("Vector valued h must be length ",len(data)-(closed?0:1)))
+    let(
+      L = len(data)
+    )
+    closed? [for(i=[0:1:L-1])
+    	        _dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
+    : [
+        (data[1]-data[0])/h[0],
+        for(i=[1:1:L-2]) _dnu_calc(data[i-1],data[i],data[i+1], h[i-1],h[i]),
+        (data[L-1]-data[L-2])/h[L-2]                            
+      ];
+
+
 // Function: deriv2()
 // Usage: deriv2(data, [h], [closed])
 // Description:
diff --git a/paths.scad b/paths.scad
index ad8b732..654a427 100644
--- a/paths.scad
+++ b/paths.scad
@@ -149,6 +149,21 @@ function path_length(path,closed=false) =
     sum([for (i = [0:1:len(path)-2]) norm(path[i+1]-path[i])])+(closed?norm(path[len(path)-1]-path[0]):0);
 
 
+// Function: path_segment_lengths()
+// Usage:
+//   path_segment_lengths(path,[closed])
+// Description:
+//   Returns list of the length of each segment in a path
+// Arguments:
+//   path = path to measure
+//   closed = true if the path is closed.  Default: false
+function path_segment_lengths(path, closed=false) =
+    [
+      for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
+      if (closed) norm(path[0]-path[len(path)-1])
+     ]; 
+
+
 // Function: path_pos_from_start()
 // Usage:
 //   pos = path_pos_from_start(path,length,[closed]);
@@ -280,13 +295,35 @@ function path_closest_point(path, pt) =
 
 
 // Function: path_tangents()
-// Usage: path_tangents(path, [closed])
+// Usage: path_tangents(path, [closed], [uniform])
 // Description:
 //   Compute the tangent vector to the input path.  The derivative approximation is described in deriv().
-//   The returns vectors will be normalized to length 1.
-function path_tangents(path, closed=false) =
+//   The returns vectors will be normalized to length 1.  If any derivatives are zero then
+//   the function fails with an error.  If you set `uniform` to false then the sampling is
+//   assumed to be non-uniform and the derivative is computed with adjustments to produce corrected
+//   values.
+// Arguments:
+//   path = path to find the tagent vectors for
+//   closed = set to true of the path is closed.  Default: false
+//   uniform = set to false to correct for non-uniform sampling.  Default: true
+// Example: A shape with non-uniform sampling gives distorted derivatives that may be undesirable
+//   rect = square([10,3]);
+//   tangents = path_tangents(rect,closed=true);
+//   stroke(rect,closed=true, width=0.1);
+//   color("purple")
+//       for(i=[0:len(tangents)-1])
+//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.1, endcap2="arrow2");
+// Example: A shape with non-uniform sampling gives distorted derivatives that may be undesirable
+//   rect = square([10,3]);
+//   tangents = path_tangents(rect,closed=true,uniform=false);
+//   stroke(rect,closed=true, width=0.1);
+//   color("purple")
+//       for(i=[0:len(tangents)-1])
+//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.1, endcap2="arrow2");
+function path_tangents(path, closed=false, uniform=true) =
     assert(is_path(path))
-    [for(t=deriv(path,closed=closed)) unit(t)];
+    !uniform ? [for(t=deriv(path,closed=closed, h=path_segment_lengths(path,closed))) unit(t)]
+             : [for(t=deriv(path,closed=closed)) unit(t)];
 
 
 // Function: path_normals()
diff --git a/vnf.scad b/vnf.scad
index 6021669..c9d56c0 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -357,17 +357,13 @@ module vnf_polyhedron(vnf, convexity=2) {
 //   Returns the volume enclosed by the given manifold VNF.   The VNF must describe a valid polyhedron with consistent face direction and
 //   no holes; otherwise the results are undefined.  Returns a positive volume if face direction is clockwise and a negative volume
 //   if face direction is counter-clockwise.
+
+// Divide the polyhedron into tetrahedra with the origin as one vertex and sum up the signed volume.
 function vnf_volume(vnf) =
     let(verts = vnf[0])
     sum([
          for(face=vnf[1], j=[1:1:len(face)-2])
-             let(
-                 v0 = verts[face[0]],
-                 v1 = verts[face[j]],
-                 v2 = verts[face[j+1]],
-                 n = cross(v2-v0,v1-v0)
-             )
-             v0 * n
+             cross(verts[face[j+1]], verts[face[j]]) * verts[face[0]]
     ])/6;
 
 
@@ -378,27 +374,22 @@ function vnf_volume(vnf) =
 //   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.
 
-// Algorithm from: https://wwwf.imperial.ac.uk/~rn/centroid.pdf
+// 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) =
     let(
         verts = vnf[0],
-        val=sum([
-            for(face=vnf[1], j=[1:1:len(face)-2])
-            let(
-                 v0 = verts[face[0]],
-                 v1 = verts[face[j]],
-                 v2 = verts[face[j+1]],
-                 n = cross(v2-v0,v1-v0)
-            ) [
-                v0 * n,
-                vmul(n,
-                    sqr(v0 + v1) +
-                    sqr(v0 + v2) +
-                    sqr(v1 + v2)
-                )
-            ]
-        ])
-    ) val[1]/val[0]/8;
+        val = sum([ for(face=vnf[1], j=[1:1:len(face)-2])
+                        let(
+                            v0  = verts[face[0]],
+                            v1  = verts[face[j]],
+                            v2  = verts[face[j+1]],
+                            vol = cross(v2,v1)*v0
+                        )
+                        [ vol, (v0+v1+v2)*vol ]
+                  ])    
+    )
+    val[1]/val[0]/4;
 
 
 function _triangulate_planar_convex_polygons(polys) =