diff --git a/tests/test_vnf.scad b/tests/test_vnf.scad
index bb171d4..b83775e 100644
--- a/tests/test_vnf.scad
+++ b/tests/test_vnf.scad
@@ -74,7 +74,8 @@ module test_vnf_centroid() {
     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();
 
 
diff --git a/vnf.scad b/vnf.scad
index 9c05e0f..7d1c021 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -450,18 +450,13 @@ 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.
 
-// Divide the solid up into tetrahedra with the origin as one vertex.  The centroid of a tetrahedron is the average of its vertices.
+// 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) =
+    assert(is_vnf(vnf) && len(vnf[0])!=0 ) 
     let(
         verts = vnf[0],
-        vol = 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]]
-            ) cross(v2,v1)*v0
-        ]),
         pos = sum([
             for(face=vnf[1], j=[1:1:len(face)-2]) let(
                 v0  = verts[face[0]],
@@ -469,10 +464,11 @@ function vnf_centroid(vnf) =
                 v2  = verts[face[j+1]],
                 vol = cross(v2,v1)*v0
             )
-            (v0+v1+v2)*vol
+            [ vol, (v0+v1+v2)*vol ]
         ])
     )
-    pos/vol/4;
+    assert(!approx(pos[0],0, EPSILON), "The vnf has self-intersections.")
+    pos[1]/pos[0]/4;
 
 
 function _triangulate_planar_convex_polygons(polys) =