diff --git a/attachments.scad b/attachments.scad
index c2c969e..48ed522 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -1235,17 +1235,20 @@ module show(tags="")
 //   }
 module diff(neg, pos=undef, keep=undef)
 {
-    difference() {
-        if (pos != undef) {
-            show(pos) children();
-        } else {
-            if (keep == undef) {
-                hide(neg) children();
+    // Don't perform the operation if the current tags are hidden
+    if (attachment_is_shown($tags)) {
+        difference() {
+            if (pos != undef) {
+                show(pos) children();
             } else {
-                hide(str(neg," ",keep)) children();
+                if (keep == undef) {
+                    hide(neg) children();
+                } else {
+                    hide(str(neg," ",keep)) children();
+                }
             }
+            show(neg) children();
         }
-        show(neg) children();
     }
     if (keep!=undef) {
         show(keep) children();
@@ -1280,17 +1283,20 @@ module diff(neg, pos=undef, keep=undef)
 //   }
 module intersect(a, b=undef, keep=undef)
 {
-    intersection() {
-        if (b != undef) {
-            show(b) children();
-        } else {
-            if (keep == undef) {
-                hide(a) children();
+    // Don't perform the operation if the current tags are hidden
+    if (attachment_is_shown($tags)) {
+        intersection() {
+            if (b != undef) {
+                show(b) children();
             } else {
-                hide(str(a," ",keep)) children();
+                if (keep == undef) {
+                    hide(a) children();
+                } else {
+                    hide(str(a," ",keep)) children();
+                }
             }
+            show(a) children();
         }
-        show(a) children();
     }
     if (keep!=undef) {
         show(keep) children();
diff --git a/common.scad b/common.scad
index eebd141..db4f3bb 100644
--- a/common.scad
+++ b/common.scad
@@ -129,11 +129,6 @@ function is_list_of(list,pattern) =
     is_list(list) &&
     []==[for(entry=0*list) if (entry != pattern) entry];
 
-function _list_pattern(list) =
-  is_list(list) ? [for(entry=list) is_list(entry) ? _list_pattern(entry) : 0]
-                : 0;
-
-
 
 // Function: is_consistent()
 // Usage:
@@ -198,11 +193,11 @@ function first_defined(v,recursive=false,_i=0) =
             is_undef(first_defined(v[_i],recursive=recursive))
         )
     )? first_defined(v,recursive=recursive,_i=_i+1) : v[_i];
-
+    
 
 // Function: one_defined()
 // Usage:
-//   one_defined(vars, names, [required])
+//   one_defined(vars, names, <required>)
 // Description:
 //   Examines the input list `vars` and returns the entry which is not `undef`.  If more
 //   than one entry is `undef` then issues an assertion specifying "Must define exactly one of" followed
@@ -221,8 +216,7 @@ function one_defined(vars, names, required=true) =
 
 // Function: num_defined()
 // Description: Counts how many items in list `v` are not `undef`.
-function num_defined(v,_i=0,_cnt=0) = _i>=len(v)? _cnt : num_defined(v,_i+1,_cnt+(is_undef(v[_i])? 0 : 1));
-
+function num_defined(v) = len([for(vi=v) if(!is_undef(vi)) 1]);
 
 // Function: any_defined()
 // Description:
@@ -239,8 +233,8 @@ function any_defined(v,recursive=false) = first_defined(v,recursive=recursive) !
 // Arguments:
 //   v = The list whose items are being checked.
 //   recursive = If true, any sublists are evaluated recursively.
-function all_defined(v,recursive=false) = max([for (x=v) is_undef(x)||(recursive&&is_list(x)&&!all_defined(x))? 1 : 0])==0;
-
+function all_defined(v,recursive=false) = 
+    []==[for (x=v) if(is_undef(x)||(recursive && is_list(x) && !all_defined(x,recursive))) 0 ];
 
 
 
@@ -249,7 +243,7 @@ function all_defined(v,recursive=false) = max([for (x=v) is_undef(x)||(recursive
 
 // Function: get_anchor()
 // Usage:
-//   get_anchor(anchor,center,[uncentered],[dflt]);
+//   get_anchor(anchor,center,<uncentered>,<dflt>);
 // Description:
 //   Calculated the correct anchor from `anchor` and `center`.  In order:
 //   - If `center` is not `undef` and `center` evaluates as true, then `CENTER` (`[0,0,0]`) is returned.
@@ -270,7 +264,7 @@ function get_anchor(anchor,center,uncentered=BOT,dflt=CENTER) =
 
 // Function: get_radius()
 // Usage:
-//   get_radius([r1], [r2], [r], [d1], [d2], [d], [dflt]);
+//   get_radius(<r1>, <r2>, <r>, <d1>, <d2>, <d>, <dflt>);
 // Description:
 //   Given various radii and diameters, returns the most specific radius.
 //   If a diameter is most specific, returns half its value, giving the radius.
@@ -288,19 +282,23 @@ function get_anchor(anchor,center,uncentered=BOT,dflt=CENTER) =
 //   r = Most general radius.
 //   d = Most general diameter.
 //   dflt = Value to return if all other values given are `undef`.
-function get_radius(r1=undef, r2=undef, r=undef, d1=undef, d2=undef, d=undef, dflt=undef) = (
-    !is_undef(r1)? assert(is_undef(r2)&&is_undef(d1)&&is_undef(d2), "Conflicting or redundant radius/diameter arguments given.") r1 :
-    !is_undef(r2)? assert(is_undef(d1)&&is_undef(d2), "Conflicting or redundant radius/diameter arguments given.") r2 :
-    !is_undef(d1)? d1/2 :
-    !is_undef(d2)? d2/2 :
-    !is_undef(r)? assert(is_undef(d), "Conflicting or redundant radius/diameter arguments given.") r :
-    !is_undef(d)? d/2 :
-    dflt
-);
+function get_radius(r1=undef, r2=undef, r=undef, d1=undef, d2=undef, d=undef, dflt=undef) = 
+    assert(num_defined([r1,d1,r2,d2])<2, "Conflicting or redundant radius/diameter arguments given.")
+    !is_undef(r1) ?   assert(is_finite(r1), "Invalid radius r1." ) r1 
+    : !is_undef(r2) ? assert(is_finite(r2), "Invalid radius r2." ) r2
+    : !is_undef(d1) ? assert(is_finite(d1), "Invalid diameter d1." ) d1/2
+    : !is_undef(d2) ? assert(is_finite(d2), "Invalid diameter d2." ) d2/2
+    : !is_undef(r)
+      ? assert(is_undef(d), "Conflicting or redundant radius/diameter arguments given.")
+        assert(is_finite(r) || is_vector(r,1) || is_vector(r,2), "Invalid radius r." )
+        r 
+    : !is_undef(d) ? assert(is_finite(d) || is_vector(d,1) || is_vector(d,2), "Invalid diameter d." ) d/2
+    : dflt;
+
 
 // Function: get_height()
 // Usage:
-//   get_height([h],[l],[height],[dflt])
+//   get_height(<h>,<l>,<height>,<dflt>)
 // Description:
 //   Given several different parameters for height check that height is not multiply defined
 //   and return a single value.  If the three values `l`, `h`, and `height` are all undefined
@@ -317,7 +315,7 @@ function get_height(h=undef,l=undef,height=undef,dflt=undef) =
 
 // Function: scalar_vec3()
 // Usage:
-//   scalar_vec3(v, [dflt]);
+//   scalar_vec3(v, <dflt>);
 // Description:
 //   If `v` is a scalar, and `dflt==undef`, returns `[v, v, v]`.
 //   If `v` is a scalar, and `dflt!=undef`, returns `[v, dflt, dflt]`.
@@ -369,7 +367,7 @@ function _valstr(x) =
 
 // Module: assert_approx()
 // Usage:
-//   assert_approx(got, expected, [info]);
+//   assert_approx(got, expected, <info>);
 // Description:
 //   Tests if the value gotten is what was expected.  If not, then
 //   the expected and received values are printed to the console and
@@ -396,7 +394,7 @@ module assert_approx(got, expected, info) {
 
 // Module: assert_equal()
 // Usage:
-//   assert_equal(got, expected, [info]);
+//   assert_equal(got, expected, <info>);
 // Description:
 //   Tests if the value gotten is what was expected.  If not, then
 //   the expected and received values are printed to the console and
@@ -423,7 +421,7 @@ module assert_equal(got, expected, info) {
 
 // Module: shape_compare()
 // Usage:
-//   shape_compare([eps]) {test_shape(); expected_shape();}
+//   shape_compare(<eps>) {test_shape(); expected_shape();}
 // Description:
 //   Compares two child shapes, returning empty geometry if they are very nearly the same shape and size.
 //   Returns the differential geometry if they are not nearly the same shape and size.
diff --git a/geometry.scad b/geometry.scad
index 443c026..dc86187 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -21,85 +21,71 @@
 //   edge = Array of two points forming the line segment to test against.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function point_on_segment2d(point, edge, eps=EPSILON) =
-    approx(point,edge[0],eps=eps) || approx(point,edge[1],eps=eps) ||  // The point is an endpoint
-    sign(edge[0].x-point.x)==sign(point.x-edge[1].x)  // point is in between the
-        && sign(edge[0].y-point.y)==sign(point.y-edge[1].y)  // edge endpoints
-        && approx(point_left_of_segment2d(point, edge),0,eps=eps);  // and on the line defined by edge
+    assert( is_vector(point,2), "Invalid point." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(edge,2,eps=eps), "Invalid segment." )
+    let( dp = point-edge[0],
+         de = edge[1]-edge[0],
+         ne = norm(de) ) 
+    ( dp*de >= -eps*ne ) 
+    && ( (dp-de)*de <= eps*ne )                  // point projects on the segment
+    && _dist2line(point-edge[0],unit(de))<eps;   // point is on the line 
+    
+    
+//Internal - distance from point `d` to the line passing through the origin with unit direction n
+//_dist2line works for any dimension
+function _dist2line(d,n) = norm(d-(d * n) * n);
+
+// Internal non-exposed function.
+function _point_above_below_segment(point, edge) =
+    let( edge = edge - [point, point] )
+    edge[0].y <= 0 
+    ?   (edge[1].y >  0 && cross(edge[0], edge[1]-edge[0]) > 0) ?  1 : 0
+    :   (edge[1].y <= 0 && cross(edge[0], edge[1]-edge[0]) < 0) ? -1 : 0 ;
+
+//Internal
+function _valid_line(line,dim,eps=EPSILON) = 
+    is_matrix(line,2,dim) 
+    && ! approx(norm(line[1]-line[0]), 0, eps); 
+    
+//Internal
+function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
 
 
-// Function: point_left_of_segment2d()
+// Function: point_left_of_line2d()
 // Usage:
-//   point_left_of_segment2d(point, edge);
+//   point_left_of_line2d(point, line);
 // Description:
-//   Return >0 if point is left of the line defined by edge.
+//   Return >0 if point is left of the line defined by `line`.
 //   Return =0 if point is on the line.
 //   Return <0 if point is right of the line.
 // Arguments:
 //   point = The point to check position of.
-//   edge = Array of two points forming the line segment to test against.
-function point_left_of_segment2d(point, edge) =
-    (edge[1].x-edge[0].x) * (point.y-edge[0].y) - (point.x-edge[0].x) * (edge[1].y-edge[0].y);
-
-
-// Internal non-exposed function.
-function _point_above_below_segment(point, edge) =
-    edge[0].y <= point.y? (
-        (edge[1].y > point.y && point_left_of_segment2d(point, edge) > 0)? 1 : 0
-    ) : (
-        (edge[1].y <= point.y && point_left_of_segment2d(point, edge) < 0)? -1 : 0
-    );
+//   line  = Array of two points forming the line segment to test against.
+function point_left_of_line2d(point, line) =
+    assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." )
+    cross(line[0]-point, line[1]-line[0]);
 
 
 // Function: collinear()
 // Usage:
-//   collinear(a, b, c, [eps]);
+//   collinear(a, [b, c], [eps]);
 // Description:
-//   Returns true if three points are co-linear.
+//   Returns true if the points `a`, `b` and `c` are co-linear or if the list of points `a` is collinear.
 // Arguments:
-//   a = First point.
-//   b = Second point.
-//   c = Third point.
+//   a = First point or list of points.
+//   b = Second point or undef; it should be undef if `c` is undef
+//   c = Third point or undef.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function collinear(a, b, c, eps=EPSILON) =
-    approx(a,b,eps=eps)? true :
-    distance_from_line([a,b], c) < eps;
-
-
-// Function: collinear_indexed()
-// Usage:
-//   collinear_indexed(points, a, b, c, [eps]);
-// Description:
-//   Returns true if three points are co-linear.
-// Arguments:
-//   points = A list of points.
-//   a = Index in `points` of first point.
-//   b = Index in `points` of second point.
-//   c = Index in `points` of third point.
-//   eps = Acceptable max angle variance.  Default: EPSILON (1e-9) degrees.
-function collinear_indexed(points, a, b, c, eps=EPSILON) =
-    let(
-        p1=points[a],
-        p2=points[b],
-        p3=points[c]
-    ) collinear(p1, p2, p3, eps);
-
-
-// Function: points_are_collinear()
-// Usage:
-//   points_are_collinear(points);
-// Description:
-//   Given a list of points, returns true if all points in the list are collinear.
-// Arguments:
-//   points = The list of points to test.
-//   eps = How much variance is allowed in testing that each point is on the same line.  Default: `EPSILON` (1e-9)
-function points_are_collinear(points, eps=EPSILON) =
-    let(
-        a = furthest_point(points[0], points),
-        b = furthest_point(points[a], points),
-        pa = points[a],
-        pb = points[b]
-    ) all([for (pt = points) collinear(pa, pb, pt, eps=eps)]);
-
+    assert( is_path([a,b,c],dim=undef)
+            || ( is_undef(b) && is_undef(c) && is_path(a,dim=undef) ), 
+            "Input should be 3 points or a list of points with same dimension.")
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( points = is_def(c) ? [a,b,c]: a )
+    len(points)<3 ? true
+    : noncollinear_triple(points,error=false,eps=eps)==[];
+ 
 
 // Function: distance_from_line()
 // Usage:
@@ -112,10 +98,11 @@ function points_are_collinear(points, eps=EPSILON) =
 // Example:
 //   distance_from_line([[-10,0], [10,0]], [3,8]);  // Returns: 8
 function distance_from_line(line, pt) =
-    let(a=line[0], n=unit(line[1]-a), d=a-pt)
-    norm(d - ((d * n) * n));
-
-
+    assert( _valid_line(line) && is_vector(pt,len(line[0])), 
+            "Invalid line, invalid point or incompatible dimensions." )
+    _dist2line(pt-line[0],unit(line[1]-line[0]));
+    
+    
 // Function: line_normal()
 // Usage:
 //   line_normal([P1,P2])
@@ -133,9 +120,11 @@ function distance_from_line(line, pt) =
 //   color("green") stroke([p1,p1+10*n], endcap2="arrow2");
 //   color("blue") move_copies([p1,p2]) circle(d=2, $fn=12);
 function line_normal(p1,p2) =
-    is_undef(p2)?
-        assert(is_path(p1,2)) line_normal(p1[0],p1[1]) :
-        assert(is_vector(p1,2)&&is_vector(p2,2)) unit([p1.y-p2.y,p2.x-p1.x]);
+    is_undef(p2)
+    ?   assert( len(p1)==2 && !is_undef(p1[1]) , "Invalid input." ) 
+        line_normal(p1[0],p1[1]) 
+    :   assert( _valid_line([p1,p2],dim=2), "Invalid line." ) 
+        unit([p1.y-p2.y,p2.x-p1.x]);
 
 
 // 2D Line intersection from two segments.
@@ -166,7 +155,10 @@ function _general_line_intersection(s1,s2,eps=EPSILON) =
 //   l2 = Second 2D line, given as a list of two 2D points on the line.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_intersection(l1,l2,eps=EPSILON) =
-    let(isect = _general_line_intersection(l1,l2,eps=eps)) isect[0];
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(l1,dim=2,eps=eps) &&_valid_line(l2,dim=2,eps=eps), "Invalid line(s)." )
+    let(isect = _general_line_intersection(l1,l2,eps=eps)) 
+    isect[0];
 
 
 // Function: line_ray_intersection()
@@ -180,9 +172,12 @@ function line_intersection(l1,l2,eps=EPSILON) =
 //   ray = The 2D ray, given as a list `[START,POINT]` of the 2D start-point START, and a 2D point POINT on the ray.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_ray_intersection(line,ray,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(line,dim=2,eps=eps) && _valid_line(ray,dim=2,eps=eps), "Invalid line or ray." )
     let(
         isect = _general_line_intersection(line,ray,eps=eps)
-    ) isect[2]<0-eps? undef : isect[0];
+    ) 
+    (isect[2]<0-eps) ? undef : isect[0];
 
 
 // Function: line_segment_intersection()
@@ -196,6 +191,8 @@ function line_ray_intersection(line,ray,eps=EPSILON) =
 //   segment = The bounded 2D line segment, given as a list of the two 2D endpoints of the segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function line_segment_intersection(line,segment,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(line,  dim=2,eps=eps) &&_valid_line(segment,dim=2,eps=eps), "Invalid line or segment." )
     let(
         isect = _general_line_intersection(line,segment,eps=eps)
     ) isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
@@ -212,9 +209,12 @@ function line_segment_intersection(line,segment,eps=EPSILON) =
 //   r2 = Second 2D ray, given as a list `[START,POINT]` of the 2D start-point START, and a 2D point POINT on the ray.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function ray_intersection(r1,r2,eps=EPSILON) =
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert( _valid_line(r1,dim=2,eps=eps) && _valid_line(r2,dim=2,eps=eps), "Invalid ray(s)." )
     let(
         isect = _general_line_intersection(r1,r2,eps=eps)
-    ) isect[1]<0-eps || isect[2]<0-eps? undef : isect[0];
+    ) 
+    isect[1]<0-eps || isect[2]<0-eps ? undef : isect[0];
 
 
 // Function: ray_segment_intersection()
@@ -228,9 +228,16 @@ function ray_intersection(r1,r2,eps=EPSILON) =
 //   segment = The bounded 2D line segment, given as a list of the two 2D endpoints of the segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function ray_segment_intersection(ray,segment,eps=EPSILON) =
+    assert( _valid_line(ray,dim=2,eps=eps) && _valid_line(segment,dim=2,eps=eps), "Invalid ray or segment." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         isect = _general_line_intersection(ray,segment,eps=eps)
-    ) isect[1]<0-eps || isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
+    ) 
+    isect[1]<0-eps 
+    || isect[2]<0-eps 
+    || isect[2]>1+eps
+    ? undef 
+    : isect[0];
 
 
 // Function: segment_intersection()
@@ -244,9 +251,17 @@ function ray_segment_intersection(ray,segment,eps=EPSILON) =
 //   s2 = Second 2D segment, given as a list of the two 2D endpoints of the line segment.
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 function segment_intersection(s1,s2,eps=EPSILON) =
+    assert( _valid_line(s1,dim=2,eps=eps) && _valid_line(s2,dim=2,eps=eps), "Invalid segment(s)." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         isect = _general_line_intersection(s1,s2,eps=eps)
-    ) isect[1]<0-eps || isect[1]>1+eps || isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
+    ) 
+    isect[1]<0-eps 
+    || isect[1]>1+eps 
+    || isect[2]<0-eps 
+    || isect[2]>1+eps 
+    ? undef 
+    : isect[0];
 
 
 // Function: line_closest_point()
@@ -301,15 +316,10 @@ function segment_intersection(s1,s2,eps=EPSILON) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function line_closest_point(line,pt) =
-    assert(is_path(line)&&len(line)==2)
-    assert(same_shape(pt,line[0]))
-    assert(!approx(line[0],line[1]))
-    let(
-        seglen = norm(line[1]-line[0]),
-        segvec = (line[1]-line[0])/seglen,
-        projection = (pt-line[0]) * segvec
-    )
-    line[0] + projection*segvec;
+    assert(_valid_line(line), "Invalid line." )
+    assert( is_vector(pt,len(line[0])), "Invalid point or incompatible dimensions." )
+    let( n = unit( line[0]- line[1]) )
+    line[1]+((pt- line[1]) * n) * n;
 
 
 // Function: ray_closest_point()
@@ -364,9 +374,8 @@ function line_closest_point(line,pt) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function ray_closest_point(ray,pt) =
-    assert(is_path(ray)&&len(ray)==2)
-    assert(same_shape(pt,ray[0]))
-    assert(!approx(ray[0],ray[1]))
+    assert( _valid_line(ray), "Invalid ray." )
+    assert(is_vector(pt,len(ray[0])), "Invalid point or incompatible dimensions." )
     let(
         seglen = norm(ray[1]-ray[0]),
         segvec = (ray[1]-ray[0])/seglen,
@@ -428,8 +437,8 @@ function ray_closest_point(ray,pt) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function segment_closest_point(seg,pt) =
-    assert(is_path(seg)&&len(seg)==2)
-    assert(same_shape(pt,seg[0]))
+    assert(_valid_line(seg), "Invalid segment." )
+    assert(len(pt)==len(seg[0]), "Incompatible dimensions." )
     approx(seg[0],seg[1])? seg[0] :
     let(
         seglen = norm(seg[1]-seg[0]),
@@ -440,6 +449,26 @@ function segment_closest_point(seg,pt) =
     projection>=seglen ? seg[1] :
     seg[0] + projection*segvec;
 
+    
+// Function: line_from_points()
+// Usage:
+//   line_from_points(points, [fast], [eps]);
+// Description:
+//   Given a list of 2 or more colinear points, returns a line containing them.
+//   If `fast` is false and the points are coincident, then `undef` is returned.
+//   if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
+// Arguments:
+//   points = The list of points to find the line through.
+//   fast = If true, don't verify that all points are collinear.  Default: false
+//   eps = How much variance is allowed in testing each point against the line.  Default: `EPSILON` (1e-9)
+function line_from_points(points, fast=false, eps=EPSILON) =
+    assert( is_path(points,dim=undef), "Improper point list." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( pb = furthest_point(points[0],points) )
+    approx(norm(points[pb]-points[0]),0) ? undef :
+    fast || collinear(points) ? [points[pb], points[0]] : undef;
+
+
 
 // Section: 2D Triangles
 
@@ -486,21 +515,30 @@ function segment_closest_point(seg,pt) =
 //   ang = tri_calc(adj=20,hyp=30)[3];
 //   ang2 = tri_calc(adj=20,hyp=40)[4];
 function tri_calc(ang,ang2,adj,opp,hyp) =
-    assert(ang==undef || ang2==undef,"You cannot specify both ang and ang2.")
-    assert(num_defined([ang,ang2,adj,opp,hyp])==2, "You must specify exactly two arguments.")
+    assert(ang==undef || ang2==undef,"At most one angle is allowed.")
+    assert(num_defined([ang,ang2,adj,opp,hyp])==2, "Exactly two arguments must be given.")
     let(
-        ang = ang!=undef? assert(ang>0&&ang<90) ang :
-            ang2!=undef? (90-ang2) :
-            adj==undef? asin(constrain(opp/hyp,-1,1)) :
-            opp==undef? acos(constrain(adj/hyp,-1,1)) :
-            atan2(opp,adj),
-        ang2 = ang2!=undef? assert(ang2>0&&ang2<90) ang2 : (90-ang),
-        adj = adj!=undef? assert(adj>0) adj :
-            (opp!=undef? (opp/tan(ang)) : (hyp*cos(ang))),
-        opp = opp!=undef? assert(opp>0) opp :
-            (adj!=undef? (adj*tan(ang)) : (hyp*sin(ang))),
-        hyp = hyp!=undef? assert(hyp>0) assert(adj<hyp) assert(opp<hyp) hyp :
-            (adj!=undef? (adj/cos(ang)) : (opp/sin(ang)))
+        ang   = ang!=undef
+                ? assert(ang>0&&ang<90, "The input angles should be acute angles." ) ang 
+                : ang2!=undef ? (90-ang2) 
+                : adj==undef ? asin(constrain(opp/hyp,-1,1)) 
+                : opp==undef ? acos(constrain(adj/hyp,-1,1)) 
+                : atan2(opp,adj),
+        ang2 =  ang2!=undef
+                ? assert(ang2>0&&ang2<90, "The input angles should be acute angles." ) ang2 
+                : (90-ang),
+        adj  =  adj!=undef
+                ? assert(adj>0, "Triangle side lengths should be positive." ) adj 
+                : (opp!=undef? (opp/tan(ang)) : (hyp*cos(ang))),
+        opp  =  opp!=undef
+                ? assert(opp>0, "Triangle side lengths should be positive." )  opp 
+                : (adj!=undef? (adj*tan(ang)) : (hyp*sin(ang))),
+        hyp  =  hyp!=undef
+                ? assert(hyp>0, "Triangle side lengths should be positive." ) 
+                  assert(adj<hyp && opp<hyp, "Hyphotenuse length should be greater than the other sides." )
+                  hyp 
+                : (adj!=undef? (adj/cos(ang)) 
+                : (opp/sin(ang)))
     )
     [adj, opp, hyp, ang, ang2];
 
@@ -517,8 +555,8 @@ function tri_calc(ang,ang2,adj,opp,hyp) =
 // Example:
 //   hyp = hyp_opp_to_adj(5,3);  // Returns: 4
 function hyp_opp_to_adj(hyp,opp) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(hyp+opp) && hyp>=0 && opp>=0, 
+           "Triangle side lengths should be a positive numbers." )
     sqrt(hyp*hyp-opp*opp);
 
 
@@ -534,8 +572,8 @@ function hyp_opp_to_adj(hyp,opp) =
 // Example:
 //   adj = hyp_ang_to_adj(8,60);  // Returns: 4
 function hyp_ang_to_adj(hyp,ang) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(hyp) && hyp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     hyp*cos(ang);
 
 
@@ -551,8 +589,8 @@ function hyp_ang_to_adj(hyp,ang) =
 // Example:
 //   adj = opp_ang_to_adj(8,30);  // Returns: 4
 function opp_ang_to_adj(opp,ang) =
-    assert(is_num(opp)&&opp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(opp) && opp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     opp/tan(ang);
 
 
@@ -567,8 +605,8 @@ function opp_ang_to_adj(opp,ang) =
 // Example:
 //   opp = hyp_adj_to_opp(5,4);  // Returns: 3
 function hyp_adj_to_opp(hyp,adj) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(adj)&&adj>=0)
+    assert(is_finite(hyp) && hyp>=0 && is_finite(adj) && adj>=0, 
+           "Triangle side lengths should be a positive numbers." )
     sqrt(hyp*hyp-adj*adj);
 
 
@@ -583,8 +621,8 @@ function hyp_adj_to_opp(hyp,adj) =
 // Example:
 //   opp = hyp_ang_to_opp(8,30);  // Returns: 4
 function hyp_ang_to_opp(hyp,ang) =
-    assert(is_num(hyp)&&hyp>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(hyp)&&hyp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     hyp*sin(ang);
 
 
@@ -599,8 +637,8 @@ function hyp_ang_to_opp(hyp,ang) =
 // Example:
 //   opp = adj_ang_to_opp(8,45);  // Returns: 8
 function adj_ang_to_opp(adj,ang) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(ang)&&ang>0&&ang<90)
+    assert(is_finite(adj)&&adj>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     adj*tan(ang);
 
 
@@ -615,8 +653,8 @@ function adj_ang_to_opp(adj,ang) =
 // Example:
 //   hyp = adj_opp_to_hyp(3,4);  // Returns: 5
 function adj_opp_to_hyp(adj,opp) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(opp) && opp>=0 && is_finite(adj) && adj>=0, 
+           "Triangle side lengths should be a positive numbers." )
     norm([opp,adj]);
 
 
@@ -631,8 +669,8 @@ function adj_opp_to_hyp(adj,opp) =
 // Example:
 //   hyp = adj_ang_to_hyp(4,60);  // Returns: 8
 function adj_ang_to_hyp(adj,ang) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(ang)&&ang>=0&&ang<90)
+    assert(is_finite(adj) && adj>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     adj/cos(ang);
 
 
@@ -647,8 +685,8 @@ function adj_ang_to_hyp(adj,ang) =
 // Example:
 //   hyp = opp_ang_to_hyp(4,30);  // Returns: 8
 function opp_ang_to_hyp(opp,ang) =
-    assert(is_num(opp)&&opp>=0)
-    assert(is_num(ang)&&ang>0&&ang<=90)
+    assert(is_finite(opp) && opp>=0, "Triangle side length should be a positive number." )
+    assert(is_finite(ang) && ang>0 && ang<90, "The angle should be an acute angle." )
     opp/sin(ang);
 
 
@@ -663,8 +701,8 @@ function opp_ang_to_hyp(opp,ang) =
 // Example:
 //   ang = hyp_adj_to_ang(8,4);  // Returns: 60 degrees
 function hyp_adj_to_ang(hyp,adj) =
-    assert(is_num(hyp)&&hyp>0)
-    assert(is_num(adj)&&adj>=0)
+    assert(is_finite(hyp) && hyp>0 && is_finite(adj) && adj>=0, 
+            "Triangle side lengths should be positive numbers." )
     acos(adj/hyp);
 
 
@@ -679,8 +717,8 @@ function hyp_adj_to_ang(hyp,adj) =
 // Example:
 //   ang = hyp_opp_to_ang(8,4);  // Returns: 30 degrees
 function hyp_opp_to_ang(hyp,opp) =
-    assert(is_num(hyp)&&hyp>0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(hyp+opp) && hyp>0 && opp>=0, 
+            "Triangle side lengths should be positive numbers." )
     asin(opp/hyp);
 
 
@@ -695,8 +733,8 @@ function hyp_opp_to_ang(hyp,opp) =
 // Example:
 //   ang = adj_opp_to_ang(sqrt(3)/2,0.5);  // Returns: 30 degrees
 function adj_opp_to_ang(adj,opp) =
-    assert(is_num(adj)&&adj>=0)
-    assert(is_num(opp)&&opp>=0)
+    assert(is_finite(adj+opp) && adj>0 && opp>=0, 
+            "Triangle side lengths should be positive numbers." )
     atan2(opp,adj);
 
 
@@ -709,55 +747,63 @@ function adj_opp_to_ang(adj,opp) =
 // Examples:
 //   triangle_area([0,0], [5,10], [10,0]);  // Returns -50
 //   triangle_area([10,0], [5,10], [0,0]);  // Returns 50
-function triangle_area(a,b,c) =
-    len(a)==3? 0.5*norm(cross(c-a,c-b)) : (
-        a.x * (b.y - c.y) +
-        b.x * (c.y - a.y) +
-        c.x * (a.y - b.y)
-    ) / 2;
+function triangle_area(a,b,c) = 
+    assert( is_path([a,b,c]), "Invalid points or incompatible dimensions." )    
+    len(a)==3 
+    ? 0.5*norm(cross(c-a,c-b)) 
+    : 0.5*cross(c-a,c-b);
 
 
 
 // Section: Planes
 
+
 // Function: plane3pt()
 // Usage:
 //   plane3pt(p1, p2, p3);
 // Description:
-//   Generates the cartesian equation of a plane from three non-collinear points on the plane.
+//   Generates the cartesian equation of a plane from three 3d points.
 //   Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane.
+//   Returns [], if the points are collinear.
 // Arguments:
 //   p1 = The first point on the plane.
 //   p2 = The second point on the plane.
 //   p3 = The third point on the plane.
 function plane3pt(p1, p2, p3) =
+    assert( is_path([p1,p2,p3],dim=3) && len(p1)==3,
+            "Invalid points or incompatible dimensions." )    
     let(
-        p1=point3d(p1),
-        p2=point3d(p2),
-        p3=point3d(p3),
-        normal = unit(cross(p3-p1, p2-p1))
-    ) concat(normal, [normal*p1]);
+        crx = cross(p3-p1, p2-p1),
+        nrm = norm(crx)
+    ) 
+    approx(nrm,0) ? [] :
+    concat(crx/nrm, [crx*p1]/nrm);
 
 
 // Function: plane3pt_indexed()
 // Usage:
 //   plane3pt_indexed(points, i1, i2, i3);
 // Description:
-//   Given a list of points, and the indices of three of those points,
+//   Given a list of 3d points, and the indices of three of those points,
 //   generates the cartesian equation of a plane that those points all
-//   lie on.  Requires that the three indexed points be non-collinear.
-//   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
+//   lie on. If the points are not collinear, returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
+//   If they are collinear, returns [].
 // Arguments:
 //   points = A list of points.
 //   i1 = The index into `points` of the first point on the plane.
 //   i2 = The index into `points` of the second point on the plane.
 //   i3 = The index into `points` of the third point on the plane.
 function plane3pt_indexed(points, i1, i2, i3) =
+    assert( is_vector([i1,i2,i3]) && min(i1,i2,i3)>=0 && is_list(points) && max(i1,i2,i3)<len(points),
+            "Invalid or out of range indices." )
+    assert( is_path([points[i1], points[i2], points[i3]],dim=3),
+            "Improper points or improper dimensions." )
     let(
         p1 = points[i1],
         p2 = points[i2],
         p3 = points[i3]
-    ) plane3pt(p1,p2,p3);
+    ) 
+    plane3pt(p1,p2,p3);
 
 
 // Function: plane_from_normal()
@@ -768,18 +814,19 @@ function plane3pt_indexed(points, i1, i2, i3) =
 // Example:
 //   plane_from_normal([0,0,1], [2,2,2]);  // Returns the xy plane passing through the point (2,2,2)
 function plane_from_normal(normal, pt=[0,0,0]) =
+  assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0), 
+          "Inputs `normal` and `pt` should 3d vectors/points and `normal` cannot be zero." )
   concat(normal, [normal*pt]);
 
 
 // Function: plane_from_points()
 // Usage:
-//   plane_from_points(points, [fast], [eps]);
+//   plane_from_points(points, <fast>, <eps>);
 // Description:
-//   Given a list of 3 or more coplanar points, returns the cartesian equation of a plane.
-//   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
-//   If not all the points in the points list are coplanar, then `undef` is returned.
-//   If `fast` is true, then a list where not all points are coplanar will result
-//   in an invalid plane value, as all coplanar checks are skipped.
+//   Given a list of 3 or more coplanar 3D points, returns the coefficients of the cartesian equation of a plane,
+//   that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
+//   If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
+//   if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned.
 // Arguments:
 //   points = The list of points to find the plane of.
 //   fast = If true, don't verify that all points in the list are coplanar.  Default: false
@@ -791,31 +838,33 @@ function plane_from_normal(normal, pt=[0,0,0]) =
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
 function plane_from_points(points, fast=false, eps=EPSILON) =
+    assert( is_path(points,dim=3), "Improper 3d point list." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         points = deduplicate(points),
-        indices = sort(find_noncollinear_points(points)),
+        indices = noncollinear_triple(points,error=false)
+    )
+    indices==[] ? undef :
+    let(
         p1 = points[indices[0]],
         p2 = points[indices[1]],
         p3 = points[indices[2]],
-        plane = plane3pt(p1,p2,p3),
-        all_coplanar = fast || all([
-            for (pt = points) coplanar(plane,pt,eps=eps)
-        ])
-    ) all_coplanar? plane : undef;
+        plane = plane3pt(p1,p2,p3)
+    )
+    fast || points_on_plane(points,plane,eps=eps) ? plane : undef;
 
 
 // Function: plane_from_polygon()
 // Usage:
 //   plane_from_polygon(points, [fast], [eps]);
 // Description:
-//   Given a 3D planar polygon, returns the cartesian equation of a plane.
+//   Given a 3D planar polygon, returns the cartesian equation of its plane.
 //   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane.
-//   If not all the points in the polygon are coplanar, then `undef` is returned.
-//   If `fast` is true, then a polygon where not all points are coplanar will
-//   result in an invalid plane value, as all coplanar checks are skipped.
+//   If not all the points in the polygon are coplanar, then [] is returned.
+//   If `fast` is true, the polygon coplanarity check is skipped and the plane may not contain all polygon points.
 // Arguments:
 //   poly = The planar 3D polygon to find the plane of.
-//   fast = If true, don't verify that all points in the polygon are coplanar.  Default: false
+//   fast = If true, doesn't verify that all points in the polygon are coplanar.  Default: false
 //   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
 // Example(3D):
 //   xyzpath = rot(45, v=[0,1,0], p=path3d(star(n=5,step=2,d=100), 70));
@@ -824,17 +873,14 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
 function plane_from_polygon(poly, fast=false, eps=EPSILON) =
+    assert( is_path(poly,dim=3), "Invalid polygon." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     let(
         poly = deduplicate(poly),
         n = polygon_normal(poly),
         plane = [n.x, n.y, n.z, n*poly[0]]
-    ) fast? plane : let(
-        all_coplanar = [
-            for (pt = poly)
-            if (!coplanar(plane,pt,eps=eps)) 1
-        ] == []
-    ) all_coplanar? plane :
-    undef;
+    ) 
+    fast? plane: coplanar(poly,eps=eps)? plane: [];
 
 
 // Function: plane_normal()
@@ -842,7 +888,9 @@ function plane_from_polygon(poly, fast=false, eps=EPSILON) =
 //   plane_normal(plane);
 // Description:
 //   Returns the unit length normal vector for the given plane.
-function plane_normal(plane) = unit([for (i=[0:2]) plane[i]]);
+function plane_normal(plane) = 
+    assert( _valid_plane(plane), "Invalid input plane." )
+    unit([plane.x, plane.y, plane.z]);
 
 
 // Function: plane_offset()
@@ -851,7 +899,9 @@ function plane_normal(plane) = unit([for (i=[0:2]) plane[i]]);
 // Description:
 //   Returns D, or the scalar offset of the plane from the origin. This can be a negative value.
 //   The absolute value of this is the distance of the plane from the origin at its closest approach.
-function plane_offset(plane) = plane[3];
+function plane_offset(plane) =  
+    assert( _valid_plane(plane), "Invalid input plane." )
+    plane[3]/norm([plane.x, plane.y, plane.z]);
 
 
 // Function: plane_transform()
@@ -859,8 +909,8 @@ function plane_offset(plane) = plane[3];
 //   mat = plane_transform(plane);
 // Description:
 //   Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, returns a 3D affine
-//   transformation matrix that will rotate and translate from points on that plane
-//   to points on the XY plane.  You can generally then use `path2d()` to drop the
+//   transformation matrix that will linear transform points on that plane
+//   into points on the XY plane.  You can generally then use `path2d()` to drop the
 //   Z coordinates, so you can work with the points in 2D.
 // Arguments:
 //   plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
@@ -875,7 +925,34 @@ function plane_transform(plane) =
     let(
         n = plane_normal(plane),
         cp = n * plane[3]
-    ) rot(from=n, to=UP) * move(-cp);
+        ) 
+    rot(from=n, to=UP) * move(-cp);
+
+
+// Function: projection_on_plane()
+// Usage:
+//   projection_on_plane(points);
+// Description:
+//   Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, and a list of 2d or 3d points, return the 3D orthogonal 
+//   projection of the points on the plane.
+// Arguments:
+//   plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
+//   points = List of points to project
+// Example(3D):
+//   points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100))));
+//   plane = plane3pt([1,0,0],[0,1,0],[0,0,1]);
+//   proj = projection_on_plane(plane,points);
+function projection_on_plane(plane, points) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( is_path(points), "Invalid list of points or dimension." )
+    let( 
+        p  = len(points[0])==2
+             ? [for(pi=points) point3d(pi) ]
+             : points, 
+        plane = plane/norm([plane.x,plane.y,plane.z]),
+        n = [plane.x,plane.y,plane.z]
+        ) 
+    [for(pi=p) pi - (pi*n - plane[3])*n];
 
 
 // Function: plane_point_nearest_origin()
@@ -899,9 +976,12 @@ function plane_point_nearest_origin(plane) =
 //   will be negative.  The normal of the plane is the same as [A,B,C].
 // Arguments:
 //   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
+//   point = The distance evaluation point.
 function distance_from_plane(plane, point) =
-    [plane.x, plane.y, plane.z] * point3d(point) - plane[3];
+    assert( _valid_plane(plane), "Invalid input plane." )
+    assert( is_vector(point,3), "The point should be a 3D point." )
+    let( nrml = [plane.x, plane.y, plane.z] )
+    ( nrml* point - plane[3])/norm(nrml);
 
 
 // Function: closest_point_on_plane()
@@ -911,34 +991,27 @@ function distance_from_plane(plane, point) =
 //   Takes a point, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
 //   Returns the coordinates of the closest point on that plane to the given `point`.
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
+//   plane = The [A,B,C,D] coefficients for the equation of the plane.
 //   point = The 3D point to find the closest point to.
 function closest_point_on_plane(plane, point) =
+    assert( _valid_plane(plane), "Invalid input plane." )
+    assert( is_vector(point,3), "Invalid point." )
     let(
-        n = unit(plane_normal(plane)),
+        n = unit([plane.x, plane.y, plane.z]),
         d = distance_from_plane(plane, point)
-    ) point - n*d;
+        ) 
+    point - n*d;
 
 
 // Returns [POINT, U] if line intersects plane at one point.
 // Returns [LINE, undef] if the line is on the plane.
 // Returns undef if line is parallel to, but not on the given plane.
 function _general_plane_line_intersection(plane, line, eps=EPSILON) =
-    let(
-        p0 = line[0],
-        p1 = line[1],
-        n = plane_normal(plane),
-        u = p1 - p0,
-        d = n * u
-    ) abs(d)<eps? (
-        coplanar(plane, p0)? [line,undef] :  // Line on plane
-        undef  // Line parallel to plane
-    ) : let(
-        v0 = closest_point_on_plane(plane, [0,0,0]),
-        w = p0 - v0,
-        s1 = (-n * w) / d,
-        pt = s1 * u + p0
-    ) [pt, s1];
+    let( a = plane*[each line[0],-1],
+         b = plane*[each(line[1]-line[0]),-1] )
+    approx(b,0,eps) 
+    ? points_on_plane(line[0],plane,eps)? [line,undef]: undef
+    : [ line[0]+a/b*(line[1]-line[0]), a/b ];
 
 
 // Function: plane_line_angle()
@@ -948,38 +1021,41 @@ function _general_plane_line_intersection(plane, line, eps=EPSILON) =
 //   The resulting angle is signed, with the sign positive if the vector p2-p1 lies on 
 //   the same side of the plane as the plane's normal vector.  
 function plane_line_angle(plane, line) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( _valid_line(line), "Invalid line." )
     let(
         vect = line[1]-line[0],
         zplane = plane_normal(plane),
         sin_angle = vect*zplane/norm(zplane)/norm(vect)
-    ) asin(constrain(sin_angle,-1,1));
+        ) 
+    asin(constrain(sin_angle,-1,1));
 
 
 // Function: plane_line_intersection()
 // Usage:
-//   pt = plane_line_intersection(plane, line, [eps]);
+//   pt = plane_line_intersection(plane, line, [bounded], [eps]);
 // Description:
 //   Takes a line, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
 //   If `line` intersects `plane` at one point, then that intersection point is returned.
 //   If `line` lies on `plane`, then the original given `line` is returned.
-//   If `line` is parallel to, but not on `plane`, then `undef` is returned.
+//   If `line` is parallel to, but not on `plane`, then undef is returned.
 // Arguments:
 //   plane = The [A,B,C,D] values for the equation of the plane.
-//   line = A list of two 3D points that are on the line.
+//   line = A list of two distinct 3D points that are on the line.
 //   bounded = If false, the line is considered unbounded.  If true, it is treated as a bounded line segment.  If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray.  Default: false (unbounded)
-//   eps = The epsilon error value to determine whether the line is too close to parallel to the plane.  Default: `EPSILON` (1e-9)
+//   eps = The tolerance value in determining whether the line is parallel to the plane.  Default: `EPSILON` (1e-9)
 function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
-    assert(is_vector(plane)&&len(plane)==4, "Invalid plane value.")
-    assert(is_path(line)&&len(line)==2, "Invalid line value.")
-    assert(!approx(line[0],line[1]), "The two points defining the line must not be the same point.")
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.")
+    assert(is_bool(bounded) || (is_list(bounded) && len(bounded)==2), "Invalid bound condition(s).")
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
         res = _general_plane_line_intersection(plane, line, eps=eps)
     )
-    is_undef(res)? undef :
-    is_undef(res[1])? res[0] :
-    bounded[0]&&res[1]<0? undef :
-    bounded[1]&&res[1]>1? undef :
+    is_undef(res) ? undef :
+    is_undef(res[1]) ? res[0] :
+    bounded[0] && res[1]<0 ? undef :
+    bounded[1] && res[1]>1 ? undef :
     res[0];
 
 
@@ -988,22 +1064,28 @@ function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
 //   pt = polygon_line_intersection(poly, line, [bounded], [eps]);
 // Description:
 //   Takes a possibly bounded line, and a 3D planar polygon, and finds their intersection point.
-//   If the line is on the plane as the polygon, and intersects, then a list of 3D line
-//   segments is returned, one for each section of the line that is inside the polygon.
-//   If the line is not on the plane of the polygon, but intersects, then the 3D intersection
-//   point is returned.  If the line does not intersect the polygon, then `undef` is returned.
+//   If the line and the polygon are on the same plane then returns a list, possibly empty, of 3D line
+//   segments, one for each section of the line that is inside the polygon.
+//   If the line is not on the plane of the polygon, but intersects it, then returns the 3D intersection
+//   point.  If the line does not intersect the polygon, then `undef` is returned.
 // Arguments:
 //   poly = The 3D planar polygon to find the intersection with.
-//   line = A list of two 3D points that are on the line.
+//   line = A list of two distinct 3D points on the line.
 //   bounded = If false, the line is considered unbounded.  If true, it is treated as a bounded line segment.  If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray.  Default: false (unbounded)
-//   eps = The epsilon error value to determine whether the line is too close to parallel to the plane.  Default: `EPSILON` (1e-9)
+//   eps = The tolerance value in determining whether the line is parallel to the plane.  Default: `EPSILON` (1e-9)
 function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
-    assert(is_path(poly))
-    assert(is_path(line)&&len(line)==2)
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    assert(is_path(poly,dim=3), "Invalid polygon." )
+    assert(is_bool(bounded) || (is_list(bounded) && len(bounded)==2), "Invalid bound condition(s).")
+    assert(_valid_line(line,dim=3,eps=eps), "Invalid line." )
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
         poly = deduplicate(poly),
-        indices = sort(find_noncollinear_points(poly)),
+        indices = noncollinear_triple(poly)
+    )
+    indices==[] ? undef :
+    let(
+        indices = sort(indices),
         p1 = poly[indices[0]],
         p2 = poly[indices[1]],
         p3 = poly[indices[2]],
@@ -1011,34 +1093,31 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
         res = _general_plane_line_intersection(plane, line, eps=eps)
     )
     is_undef(res)? undef :
-    is_undef(res[1])? (
-        let(
-            // Line is on polygon plane.
+    is_undef(res[1])
+    ? ( let(// Line is on polygon plane.
             linevec = unit(line[1] - line[0]),
             lp1 = line[0] + (bounded[0]? 0 : -1000000) * linevec,
             lp2 = line[1] + (bounded[1]? 0 :  1000000) * linevec,
             poly2d = clockwise_polygon(project_plane(poly, p1, p2, p3)),
             line2d = project_plane([lp1,lp2], p1, p2, p3),
             parts = split_path_at_region_crossings(line2d, [poly2d], closed=false),
-            inside = [
-                for (part = parts)
-                if (point_in_polygon(mean(part), poly2d)>0) part
-            ]
-        ) !inside? undef :
+            inside = [for (part = parts)
+                          if (point_in_polygon(mean(part), poly2d)>0) part
+                     ]
+        ) 
+        !inside? undef :
         let(
-            isegs = [
-                for (seg = inside)
-                lift_plane(seg, p1, p2, p3)
-            ]
-        ) isegs
-    ) :
-    bounded[0]&&res[1]<0? undef :
-    bounded[1]&&res[1]>1? undef :
-    let(
-        proj = clockwise_polygon(project_plane(poly, p1, p2, p3)),
-        pt = project_plane(res[0], p1, p2, p3)
-    ) point_in_polygon(pt, proj) < 0? undef :
-    res[0];
+            isegs = [for (seg = inside) lift_plane(seg, p1, p2, p3) ]
+        ) 
+        isegs
+    ) 
+    :   bounded[0]&&res[1]<0? [] :
+        bounded[1]&&res[1]>1? [] :
+        let(
+            proj = clockwise_polygon(project_plane(poly, p1, p2, p3)),
+            pt = project_plane(res[0], p1, p2, p3)
+        ) 
+        point_in_polygon(pt, proj) < 0 ? undef : res[0];
 
 
 // Function: plane_intersection()
@@ -1047,53 +1126,62 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
 // Description:
 //   Compute the point which is the intersection of the three planes, or the line intersection of two planes.
 //   If you give three planes the intersection is returned as a point.  If you give two planes the intersection
-//   is returned as a list of two points on the line of intersection.  If any of the input planes are parallel
-//   then returns undef.  
+//   is returned as a list of two points on the line of intersection.  If any two input planes are parallel
+//   or coincident then returns undef.  
 function plane_intersection(plane1,plane2,plane3) =
-    is_def(plane3)? let(
-        matrix = [for(p=[plane1,plane2,plane3]) select(p,0,2)],
-        rhs = [for(p=[plane1,plane2,plane3]) p[3]]
-    ) linear_solve(matrix,rhs) :
-    let(
-        normal = cross(plane_normal(plane1), plane_normal(plane2))
-    ) approx(norm(normal),0) ? undef :
-    let(
-        matrix = [for(p=[plane1,plane2]) select(p,0,2)],
-        rhs = [for(p=[plane1,plane2]) p[3]],
-        point = linear_solve(matrix,rhs)
-    ) is_undef(point)? undef :
-    [point, point+normal];
+    assert( _valid_plane(plane1) && _valid_plane(plane2) && (is_undef(plane3) ||_valid_plane(plane3)),
+                "The input must be 2 or 3 planes." )
+    is_def(plane3)
+    ?   let(
+          matrix = [for(p=[plane1,plane2,plane3]) select(p,0,2)],
+          rhs = [for(p=[plane1,plane2,plane3]) p[3]]
+        ) 
+        linear_solve(matrix,rhs)
+    :   let( normal = cross(plane_normal(plane1), plane_normal(plane2)) ) 
+        approx(norm(normal),0) ? undef :
+        let(
+            matrix = [for(p=[plane1,plane2]) select(p,0,2)],
+            rhs = [for(p=[plane1,plane2]) p[3]],
+            point = linear_solve(matrix,rhs)
+        ) 
+        point==[]? undef: [point, point+normal];
 
 
 // Function: coplanar()
 // Usage:
-//   coplanar(plane, point);
+//   coplanar(points,<eps>);
 // Description:
-//   Given a plane as [A,B,C,D] where the cartesian equation for that plane
-//   is Ax+By+Cz=D, determines if the given point is on that plane.
-//   Returns true if the point is on that plane.
+//   Returns true if the given 3D points are non-collinear and are on a plane.
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
-//   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
-function coplanar(plane, point, eps=EPSILON) =
-    abs(distance_from_plane(plane, point)) <= eps;
+//   points = The points to test.
+//   eps = How much variance is allowed in the planarity test.  Default: `EPSILON` (1e-9)
+function coplanar(points, eps=EPSILON) =
+    assert( is_path(points,dim=3) , "Input should be a list of 3D points." )
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a non-negative number." )
+    len(points)<=2 ? false
+    :   let( ip = noncollinear_triple(points,error=false,eps=eps) ) 
+        ip == [] ? false : 
+        let( plane  = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]), 
+             normal = point3d(plane) )
+        max( points*normal ) - plane[3]< eps*norm(normal);
 
-
-// Function: points_are_coplanar()
+    
+// Function: points_on_plane()
 // Usage:
-//   points_are_coplanar(points, [eps]);
+//   points_on_plane(points, plane, <eps>);
 // Description:
-//   Given a list of points, returns true if all points in the list are coplanar.
+//   Returns true if the given 3D points are on the given plane.
 // Arguments:
-//   points = The list of points to test.
-//   eps = How much variance is allowed in testing that each point is on the same plane.  Default: `EPSILON` (1e-9)
-function points_are_coplanar(points, eps=EPSILON) =
-    points_are_collinear(points, eps=eps)? true :
-    let(
-        plane = plane_from_points(points, fast=true, eps=eps)
-    ) all([for (pt = points) coplanar(plane, pt, eps=eps)]);
-
+//   plane = The plane to test the points on.
+//   points = The list of 3D points to test.
+//   eps = How much variance is allowed in the planarity testing.  Default: `EPSILON` (1e-9)
+function points_on_plane(points, plane, eps=EPSILON) =
+    assert( _valid_plane(plane), "Invalid plane." )
+    assert( is_matrix(points,undef,3) && len(points)>0, "Invalid pointlist." ) // using is_matrix it accepts len(points)==1
+    assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
+    let( normal = point3d(plane),
+         pt_nrm = points*normal )
+    abs(max( max(pt_nrm) - plane[3], -min(pt_nrm)+plane[3]))< eps*norm(normal);
 
 
 // Function: in_front_of_plane()
@@ -1101,12 +1189,12 @@ function points_are_coplanar(points, eps=EPSILON) =
 //   in_front_of_plane(plane, point);
 // Description:
 //   Given a plane as [A,B,C,D] where the cartesian equation for that plane
-//   is Ax+By+Cz=D, determines if the given point is on the side of that
+//   is Ax+By+Cz=D, determines if the given 3D point is on the side of that
 //   plane that the normal points towards.  The normal of the plane is the
 //   same as [A,B,C].
 // Arguments:
-//   plane = The [A,B,C,D] values for the equation of the plane.
-//   point = The point to test.
+//   plane = The [A,B,C,D] coefficients for the equation of the plane.
+//   point = The 3D point to test.
 function in_front_of_plane(plane, point) =
     distance_from_plane(plane, point) > EPSILON;
 
@@ -1116,7 +1204,7 @@ function in_front_of_plane(plane, point) =
 
 // Function: find_circle_2tangents()
 // Usage:
-//   find_circle_2tangents(pt1, pt2, pt3, r|d, [tangents]);
+//   find_circle_2tangents(pt1, pt2, pt3, r|d, <tangents>);
 // Description:
 //   Given a pair of rays with a common origin, and a known circle radius/diameter, finds
 //   the centerpoint for the circle of that size that touches both rays tangentally.
@@ -1156,37 +1244,46 @@ function in_front_of_plane(plane, point) =
 function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) =
     let(r = get_radius(r=r, d=d, dflt=undef))
     assert(r!=undef, "Must specify either r or d.")
-    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))? find_circle_2tangents(pt1[0], pt1[1], pt1[2], r=r) :
-    collinear(pt1, pt2, pt3)? undef :
-    let(
-        v1 = unit(pt1 - pt2),
-        v2 = unit(pt3 - pt2),
-        vmid = unit(mean([v1, v2])),
-        n = vector_axis(v1, v2),
-        a = vector_angle(v1, v2),
-        hyp = r / sin(a/2),
-        cp = pt2 + hyp * vmid
-    ) !tangents? [cp, n] :
-    let(
-        x = hyp * cos(a/2),
-        tp1 = pt2 + x * v1,
-        tp2 = pt2 + x * v2,
-        fff=echo(tp1=tp1,cp=cp,pt2=pt2),
-        dang1 = vector_angle(tp1-cp,pt2-cp),
-        dang2 = vector_angle(tp2-cp,pt2-cp)
-    ) [cp, n, tp1, tp2, dang1, dang2];
+    assert( ( is_path(pt1) && len(pt1)==3 && is_undef(pt2) && is_undef(pt3)) 
+            || (is_matrix([pt1,pt2,pt3]) && (len(pt1)==2 || len(pt1)==3) ),
+            "Invalid input points." )
+    is_undef(pt2) 
+    ? find_circle_2tangents(pt1[0], pt1[1], pt1[2], r=r, tangents=tangents) 
+    : collinear(pt1, pt2, pt3)? undef :
+        let(
+            v1 = unit(pt1 - pt2),
+            v2 = unit(pt3 - pt2),
+            vmid = unit(mean([v1, v2])),
+            n = vector_axis(v1, v2),
+            a = vector_angle(v1, v2),
+            hyp = r / sin(a/2),
+            cp = pt2 + hyp * vmid
+            ) 
+        !tangents ? [cp, n] :
+        let(
+            x = hyp * cos(a/2),
+            tp1 = pt2 + x * v1,
+            tp2 = pt2 + x * v2,
+//            fff=echo(tp1=tp1,cp=cp,pt2=pt2),
+            dang1 = vector_angle(tp1-cp,pt2-cp),
+            dang2 = vector_angle(tp2-cp,pt2-cp)
+            ) 
+        [cp, n, tp1, tp2, dang1, dang2];
 
 
 // Function: find_circle_3points()
 // Usage:
 //   find_circle_3points(pt1, pt2, pt3);
+//   find_circle_3points([pt1, pt2, pt3]);
 // Description:
 //   Returns the [CENTERPOINT, RADIUS, NORMAL] of the circle that passes through three non-collinear
-//   points.  The centerpoint will be a 2D or 3D vector, depending on the points input.  If all three
+//   points where NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
+//   The centerpoint will be a 2D or 3D vector, depending on the points input.  If all three
 //   points are 2D, then the resulting centerpoint will be 2D, and the normal will be UP ([0,0,1]).
 //   If any of the points are 3D, then the resulting centerpoint will be 3D.  If the three points are
 //   collinear, then `[undef,undef,undef]` will be returned.  The normal will be a normalized 3D
 //   vector with a non-negative Z axis.
+//   Instead of 3 arguments, it is acceptable to input the 3 points in a list `pt1`, leaving `pt2`and `pt3` as undef.
 // Arguments:
 //   pt1 = The first point.
 //   pt2 = The second point.
@@ -1198,49 +1295,43 @@ function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) =
 //   translate(circ[0]) color("red") circle(d=3, $fn=12);
 //   move_copies(pts) color("blue") circle(d=3, $fn=12);
 function find_circle_3points(pt1, pt2, pt3) =
-    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))? find_circle_3points(pt1[0], pt1[1], pt1[2]) :
-    collinear(pt1,pt2,pt3)? [undef,undef,undef] :
-    let(
-        v1 = pt1-pt2,
-        v2 = pt3-pt2,
-        n = vector_axis(v1,v2),
-        n2 = n.z<0? -n : n
-    ) len(pt1)+len(pt2)+len(pt3)>6? (
+    (is_undef(pt2) && is_undef(pt3) && is_list(pt1))
+    ? find_circle_3points(pt1[0], pt1[1], pt1[2]) 
+    :   assert( is_vector(pt1) && is_vector(pt2) && is_vector(pt3) 
+                && max(len(pt1),len(pt2),len(pt3))<=3 && min(len(pt1),len(pt2),len(pt3))>=2,
+                "Invalid point(s)." )
+        collinear(pt1,pt2,pt3)? [undef,undef,undef] :
         let(
-            a = project_plane(pt1, pt1, pt2, pt3),
-            b = project_plane(pt2, pt1, pt2, pt3),
-            c = project_plane(pt3, pt1, pt2, pt3),
-            res = find_circle_3points(a, b, c)
-        ) res[0]==undef? [undef,undef,undef] : let(
-            cp = lift_plane(res[0], pt1, pt2, pt3),
-            r = norm(pt2-cp)
-        ) [cp, r, n2]
-    ) : let(
-        mp1 = pt2 + v1/2,
-        mp2 = pt2 + v2/2,
-        mpv1 = rot(90, v=n, p=v1),
-        mpv2 = rot(90, v=n, p=v2),
-        l1 = [mp1, mp1+mpv1],
-        l2 = [mp2, mp2+mpv2],
-        isect = line_intersection(l1,l2)
-    ) is_undef(isect)? [undef,undef,undef] : let(
-        r = norm(pt2-isect)
-    ) [isect, r, n2];
-
-
+            v  = [ point3d(pt1), point3d(pt2), point3d(pt3) ], // triangle vertices
+            ed = [for(i=[0:2]) v[(i+1)%3]-v[i] ],    // triangle edge vectors
+            pm = [for(i=[0:2]) v[(i+1)%3]+v[i] ]/2,  // edge mean points
+            es = sortidx( [for(di=ed) norm(di) ] ),   
+            e1 = ed[es[1]],                          // take the 2 longest edges
+            e2 = ed[es[2]],
+            n0 = vector_axis(e1,e2),                 // normal standardization 
+            n  = n0.z<0? -n0 : n0,
+            sc = plane_intersection(                 
+                    [ each e1, e1*pm[es[1]] ],       // planes orthogonal to 2 edges
+                    [ each e2, e2*pm[es[2]] ],
+                    [ each n,  n*v[0] ] ) ,          // triangle plane
+            cp = len(pt1)+len(pt2)+len(pt3)>6 ? sc: [sc.x, sc.y], 
+            r  = norm(sc-v[0])
+            )
+        [ cp, r, n ];
+     
 
 // Function: circle_point_tangents()
 // Usage:
 //   tangents = circle_point_tangents(r|d, cp, pt);
 // Description:
-//   Given a circle and a point outside that circle, finds the tangent point(s) on the circle for a
+//   Given a 2d circle and a 2d point outside that circle, finds the 2d tangent point(s) on the circle for a
 //   line passing through the point.  Returns list of zero or more sublists of [ANG, TANGPT]
 // Arguments:
 //   r = Radius of the circle.
 //   d = Diameter of the circle.
-//   cp = The coordinates of the circle centerpoint.
-//   pt = The coordinates of the external point.
-// Example(2D):
+//   cp = The coordinates of the 2d circle centerpoint.
+//   pt = The coordinates of the 2d external point.
+// Example:
 //   cp = [-10,-10];  r = 30;  pt = [30,10];
 //   tanpts = subindex(circle_point_tangents(r=r, cp=cp, pt=pt),1);
 //   color("yellow") translate(cp) circle(r=r);
@@ -1248,9 +1339,8 @@ function find_circle_3points(pt1, pt2, pt3) =
 //   color("red") move_copies(tanpts) circle(d=3,$fn=12);
 //   color("blue") move_copies([cp,pt]) circle(d=3,$fn=12);
 function circle_point_tangents(r, d, cp, pt) =
-    assert(is_num(r) || is_num(d))
-    assert(is_vector(cp))
-    assert(is_vector(pt))
+    assert(is_finite(r) || is_finite(d), "Invalid radius or diameter." )
+    assert(is_path([cp, pt],dim=2), "Invalid center point or external point.")
     let(
         r = get_radius(r=r, d=d, dflt=1),
         delta = pt - cp,
@@ -1264,11 +1354,10 @@ function circle_point_tangents(r, d, cp, pt) =
     ) [for (ang=angs) [ang, cp + r*[cos(ang),sin(ang)]]];
 
 
-
 // Function: circle_circle_tangents()
 // Usage: circle_circle_tangents(c1, r1|d1, c2, r2|d2)
 // Description:
-//   Computes lines tangents to a pair of circles.  Returns a list of line endpoints [p1,p2] where
+//   Computes 2d lines tangents to a pair of circles in 2d.  Returns a list of line endpoints [p1,p2] where
 //   p2 is the tangent point on circle 1 and p2 is the tangent point on circle 2.
 //   If four tangents exist then the first one the left hand exterior tangent as regarded looking from
 //   circle 1 toward circle 2.  The second value is the right hand exterior tangent.  The third entry
@@ -1314,6 +1403,7 @@ function circle_point_tangents(r, d, cp, pt) =
 //   move(c2) stroke(circle(r=r2), width=.1, closed=true);
 //   echo(pts);   // Returns []
 function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
+    assert( is_path([c1,c2],dim=2), "Invalid center point(s)." )
     let(
         r1 = get_radius(r1=r1,d1=d1),
         r2 = get_radius(r1=r2,d1=d2),
@@ -1342,24 +1432,32 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
 // Section: Pointlists
 
 
-// Function: find_noncollinear_points()
+// Function: noncollinear_triple()
 // Usage:
-//   find_noncollinear_points(points);
+//   noncollinear_triple(points);
 // Description:
 //   Finds the indices of three good non-collinear points from the points list `points`.
-function find_noncollinear_points(points,error=true,eps=EPSILON) =
+//   If all points are collinear, returns [].
+function noncollinear_triple(points,error=true,eps=EPSILON) =
+    assert( is_path(points), "Invalid input points." )
+    assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative number." )
     let(
         pa = points[0],
-        b = furthest_point(pa, points),
-        n = unit(points[b]-pa),
-        relpoints = [for(pt=points) pt-pa],
-        proj = relpoints * n,
-        distlist = [for(i=[0:len(points)-1]) norm(relpoints[i]-proj[i]*n)]
-    )
-    max(distlist)<eps
+        b  = furthest_point(pa, points),
+        pb = points[b],
+        nrm = norm(pa-pb)
+        )
+    approx(nrm, 0)
     ? assert(!error, "Cannot find three noncollinear points in pointlist.")
-      []
-    : [0,b,max_index(distlist)];
+        []
+    :   let(
+            n = (pb-pa)/nrm,
+            distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]   
+           )
+        max(distlist)<eps
+        ?  assert(!error, "Cannot find three noncollinear points in pointlist.")
+           []
+        :  [0,b,max_index(distlist)];    
 
 
 // Function: pointlist_bounds()
@@ -1372,13 +1470,14 @@ function find_noncollinear_points(points,error=true,eps=EPSILON) =
 // Arguments:
 //   pts = List of points.
 function pointlist_bounds(pts) =
-    assert(is_matrix(pts))
+    assert(is_matrix(pts) && len(pts)>0 && len(pts[0])>0 , "Invalid pointlist." ) 
     let(ptsT = transpose(pts))
     [
       [for(row=ptsT) min(row)],
       [for(row=ptsT) max(row)]
     ];
 
+
 // Function: closest_point()
 // Usage:
 //   closest_point(pt, points);
@@ -1388,6 +1487,8 @@ function pointlist_bounds(pts) =
 //   pt = The point to find the closest point to.
 //   points = The list of points to search.
 function closest_point(pt, points) =
+    assert( is_vector(pt), "Invalid point." )
+    assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
     min_index([for (p=points) norm(p-pt)]);
 
 
@@ -1400,6 +1501,8 @@ function closest_point(pt, points) =
 //   pt = The point to find the farthest point from.
 //   points = The list of points to search.
 function furthest_point(pt, points) =
+    assert( is_vector(pt), "Invalid point." )
+    assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
     max_index([for (p=points) norm(p-pt)]);
 
 
@@ -1410,37 +1513,48 @@ function furthest_point(pt, points) =
 // Usage:
 //   area = polygon_area(poly);
 // Description:
-//   Given a 2D or 3D planar polygon, returns the area of that polygon.  If the polygon is self-crossing, the results are undefined.
-function polygon_area(poly) =
-    len(poly)<3? 0 :
-    len(poly[0])==2? 0.5*sum([for(i=[0:1:len(poly)-1]) det2(select(poly,i,i+1))]) :
-    let(
-        plane = plane_from_points(poly)
-    ) plane==undef? undef :
-    let(
-        n = unit(plane_normal(plane)),
-        total = sum([for (i=[0:1:len(poly)-1]) cross(poly[i], select(poly,i+1))]),
-        res = abs(total * n) / 2
-    ) res;
+//   Given a 2D or 3D planar polygon, returns the area of that polygon.  
+//   If the polygon is self-crossing, the results are undefined. For non-planar points the result is undef.
+//   When `signed` is true, a signed area is returned; a positive area indicates a counterclockwise polygon.
+// Arguments:
+//   poly = polygon to compute the area of.
+//   signed = if true, a signed area is returned (default: false)
+function polygon_area(poly, signed=false) =
+    assert(is_path(poly), "Invalid polygon." )
+    len(poly)<3 ? 0 :
+    len(poly[0])==2
+    ? sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2
+    : let( plane = plane_from_points(poly) )
+      plane==undef? undef :
+      let( n = unit(plane_normal(plane)), 
+           total = sum([for(i=[1:1:len(poly)-1]) cross(poly[i]-poly[0],poly[i+1]-poly[0])*n ])/2
+          )
+      signed ? total : abs(total);
 
 
-// Function: polygon_is_convex()
+// Function: is_convex_polygon()
 // Usage:
-//   polygon_is_convex(poly);
+//   is_convex_polygon(poly);
 // Description:
-//   Returns true if the given polygon is convex.  Result is undefined if the polygon is self-intersecting.
+//   Returns true if the given 2D polygon is convex.  The result is meaningless if the polygon is not simple (self-intersecting).
+//   If the points are collinear the result is true. 
 // Example:
-//   polygon_is_convex(circle(d=50));  // Returns: true
+//   is_convex_polygon(circle(d=50));  // Returns: true
 // Example:
 //   spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
-//   polygon_is_convex(spiral);  // Returns: false
-function polygon_is_convex(poly) =
-    let(
-        l = len(poly),
-        c = [for (i=idx(poly)) cross(poly[(i+1)%l]-poly[i],poly[(i+2)%l]-poly[(i+1)%l])]
-    )
-    len([for (x=c) if(x>0) 1])==0 ||
-    len([for (x=c) if(x<0) 1])==0;
+//   is_convex_polygon(spiral);  // Returns: false
+function is_convex_polygon(poly) =
+    assert(is_path(poly,dim=2), "The input should be a 2D polygon." )
+    let( l = len(poly) )
+    len([for( i = l-1, 
+              c = cross(poly[(i+1)%l]-poly[i], poly[(i+2)%l]-poly[(i+1)%l]), 
+              s = sign(c);
+            i>=0 && sign(c)==s;
+              i = i-1, 
+              c = i<0? 0: cross(poly[(i+1)%l]-poly[i],poly[(i+2)%l]-poly[(i+1)%l]), 
+              s = s==0 ? sign(c) : s
+            ) i 
+        ])== l;
 
 
 // Function: polygon_shift()
@@ -1454,6 +1568,7 @@ function polygon_is_convex(poly) =
 // Example:
 //   polygon_shift([[3,4], [8,2], [0,2], [-4,0]], 2);   // Returns [[0,2], [-4,0], [3,4], [8,2]]
 function polygon_shift(poly, i) =
+    assert(is_path(poly), "Invalid polygon." )
     list_rotate(cleanup_path(poly), i);
 
 
@@ -1461,13 +1576,15 @@ function polygon_shift(poly, i) =
 // Usage:
 //   polygon_shift_to_closest_point(path, pt);
 // Description:
-//   Given a polygon `path`, rotates the point ordering so that the first point in the path is the one closest to the given point `pt`.
-function polygon_shift_to_closest_point(path, pt) =
+//   Given a polygon `poly`, rotates the point ordering so that the first point in the path is the one closest to the given point `pt`.
+function polygon_shift_to_closest_point(poly, pt) =
+    assert(is_vector(pt), "Invalid point." )
+    assert(is_path(poly,dim=len(pt)), "Invalid polygon or incompatible dimension with the point." )
     let(
-        path = cleanup_path(path),
-        dists = [for (p=path) norm(p-pt)],
+        poly = cleanup_path(poly),
+        dists = [for (p=poly) norm(p-pt)],
         closest = min_index(dists)
-    ) select(path,closest,closest+len(path)-1);
+    ) select(poly,closest,closest+len(poly)-1);
 
 
 // Function: reindex_polygon()
@@ -1500,38 +1617,32 @@ function polygon_shift_to_closest_point(path, pt) =
 //   color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
 //   color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
 function reindex_polygon(reference, poly, return_error=false) = 
-    assert(is_path(reference) && is_path(poly))
-    assert(len(reference)==len(poly), "Polygons must be the same length in reindex_polygon")
+    assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
+           "Invalid polygon(s) or incompatible dimensions. " )
+    assert(len(reference)==len(poly), "The polygons must have the same length.")
     let(
         dim = len(reference[0]),
         N = len(reference),
         fixpoly = dim != 2? poly :
-            polygon_is_clockwise(reference)? clockwise_polygon(poly) :
-            ccw_polygon(poly),
-        dist = [
-            // Matrix of all pairwise distances
-            for (p1=reference) [
-                for (p2=fixpoly) norm(p1-p2)
-            ]
-        ],
-        // Compute the sum of all distance pairs for a each shift
-        sums = [
-            for(shift=[0:1:N-1]) sum([
-                for(i=[0:1:N-1]) dist[i][(i+shift)%N]
-            ])
-        ],
-        optimal_poly = polygon_shift(fixpoly,min_index(sums))
-    )
-    return_error? [optimal_poly, min(sums)] :
+                  polygon_is_clockwise(reference)
+                  ? clockwise_polygon(poly)
+                  : ccw_polygon(poly),
+        I   = [for(i=[0:N-1]) 1],
+        val = [ for(k=[0:N-1]) 
+                  [for(i=[0:N-1]) 
+                     (reference[i]*poly[(i+k)%N]) ] ]*I,
+        optimal_poly = polygon_shift(fixpoly, max_index(val))
+      )
+    return_error? [optimal_poly, min(poly*(I*poly)-2*val)] :
     optimal_poly;
-
+    
 
 // Function: align_polygon()
 // Usage:
-//   newpoly = align_polygon(reference, poly, angles, [cp]);
+//   newpoly = align_polygon(reference, poly, angles, <cp>);
 // Description:
-//   Tries the list or range of angles to find a rotation of the specified polygon that best aligns
-//   with the reference polygon.  For each angle, the polygon is reindexed, which is a costly operation
+//   Tries the list or range of angles to find a rotation of the specified 2D polygon that best aligns
+//   with the reference 2D polygon.  For each angle, the polygon is reindexed, which is a costly operation
 //   so if run time is a problem, use a smaller sampling of angles.  Returns the rotated and reindexed
 //   polygon.
 // Arguments:
@@ -1546,9 +1657,11 @@ function reindex_polygon(reference, poly, return_error=false) =
 //   color("red") move_copies(scale(1.4,p=align_polygon(pentagon,hexagon,[0:10:359]))) circle(r=.1);
 //   move_copies(concat(pentagon,hexagon))circle(r=.1);
 function align_polygon(reference, poly, angles, cp) =
-    assert(is_path(reference) && is_path(poly))
-    assert(len(reference)==len(poly), "Polygons must be the same length to be aligned in align_polygon")
-    assert(is_num(angles[0]), "The `angle` parameter to align_polygon must be a range or vector")
+    assert(is_path(reference,dim=2) && is_path(poly,dim=2),
+           "Invalid polygon(s). " )
+    assert(len(reference)==len(poly), "The polygons must have the same length.")
+    assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
+            "The `angle` parameter must be a range or a non void list of numbers.")
     let(     // alignments is a vector of entries of the form: [polygon, error]
         alignments = [
             for(angle=angles) reindex_polygon(
@@ -1569,30 +1682,37 @@ function align_polygon(reference, poly, angles, cp) =
 //   Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
 //   If the polygon is self-intersecting, the results are undefined.
 function centroid(poly) =
-    len(poly[0])==2? (
-        sum([
+    assert( is_path(poly,dim=[2,3]), "The input must be a 2D or 3D polygon." )
+    len(poly[0])==2
+    ? sum([
             for(i=[0:len(poly)-1])
             let(segment=select(poly,i,i+1))
             det2(segment)*sum(segment)
-        ]) / 6 / polygon_area(poly)
-    ) : (
-        let(
-            n = plane_normal(plane_from_points(poly)),
-            p1 = vector_angle(n,UP)>15? vector_axis(n,UP) : vector_axis(n,RIGHT),
-            p2 = vector_axis(n,p1),
-            cp = mean(poly),
-            proj = project_plane(poly,cp,cp+p1,cp+p2),
-            cxy = centroid(proj)
-        ) lift_plane(cxy,cp,cp+p1,cp+p2)
-    );
-
+          ]) / 6 / polygon_area(poly)
+    : let( plane = plane_from_points(poly, fast=true) )
+      assert( !is_undef(plane), "The polygon must be planar." )
+      let(
+        n = plane_normal(plane),
+        val = sum([for(i=[1:len(poly)-2])
+                      let(
+                         v0 = poly[0],
+                         v1 = poly[i],
+                         v2 = poly[i+1],
+                         area = cross(v2-v0,v1-v0)*n
+                         )
+                      [ area, (v0+v1+v2)*area ]
+                  ] )
+          )
+      val[1]/val[0]/3;
+            
 
 // Function: point_in_polygon()
 // Usage:
-//   point_in_polygon(point, path, [eps])
+//   point_in_polygon(point, poly, <eps>)
 // Description:
 //   This function tests whether the given 2D point is inside, outside or on the boundary of
-//   the specified 2D polygon using the Winding Number method.
+//   the specified 2D polygon using either the Nonzero Winding rule or the Even-Odd rule.
+//   See https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Even–odd_rule.
 //   The polygon is given as a list of 2D points, not including the repeated end point.
 //   Returns -1 if the point is outside the polyon.
 //   Returns 0 if the point is on the boundary.
@@ -1602,52 +1722,81 @@ function centroid(poly) =
 //   Rounding error may give mixed results for points on or near the boundary.
 // Arguments:
 //   point = The 2D point to check position of.
-//   path = The list of 2D path points forming the perimeter of the polygon.
+//   poly = The list of 2D path points forming the perimeter of the polygon.
+//   nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd" (Default: true )
 //   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
-function point_in_polygon(point, path, eps=EPSILON) =
-    // Original algorithm from http://geomalgorithms.com/a03-_inclusion.html
+function point_in_polygon(point, poly, eps=EPSILON, nonzero=true) =
+    // Original algorithms from http://geomalgorithms.com/a03-_inclusion.html
+    assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
+            "The point and polygon should be in 2D. The polygon should have more that 2 points." )
+    assert( is_finite(eps) && eps>=0, "Invalid tolerance." )
     // Does the point lie on any edges?  If so return 0.
-    sum([for(i=[0:1:len(path)-1]) let(seg=select(path,i,i+1)) if(!approx(seg[0],seg[1],eps=eps)) point_on_segment2d(point, seg, eps=eps)?1:0]) > 0? 0 :
-    // Otherwise compute winding number and return 1 for interior, -1 for exterior
-    sum([for(i=[0:1:len(path)-1]) let(seg=select(path,i,i+1)) if(!approx(seg[0],seg[1],eps=eps)) _point_above_below_segment(point, seg)]) != 0? 1 : -1;
+    let(
+        on_brd = [for(i=[0:1:len(poly)-1]) 
+                    let( seg = select(poly,i,i+1) ) 
+                    if( !approx(seg[0],seg[1],eps=EPSILON) ) 
+                        point_on_segment2d(point, seg, eps=eps)? 1:0 ]
+        )
+    sum(on_brd) > 0
+    ? 0 
+    :   nonzero
+        ?    // Compute winding number and return 1 for interior, -1 for exterior
+            let(
+                windchk = [for(i=[0:1:len(poly)-1]) 
+                            let(seg=select(poly,i,i+1)) 
+                            if(!approx(seg[0],seg[1],eps=eps)) 
+                                _point_above_below_segment(point, seg)
+                          ]
+                )
+            sum(windchk) != 0 ? 1 : -1
+        :   // or compute the crossings with the ray [point, point+[1,0]]
+            let( 
+              n  = len(poly),
+              cross = 
+                [for(i=[0:n-1])
+                    let( 
+                      p0 = poly[i]-point, 
+                      p1 = poly[(i+1)%n]-point
+                      )
+                    if( ( (p1.y>eps && p0.y<=0) || (p1.y<=0 && p0.y>eps) )
+                        &&  0 < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y) )
+                    1
+                ]
+            )
+            2*(len(cross)%2)-1;;
 
 
 // Function: polygon_is_clockwise()
 // Usage:
-//   polygon_is_clockwise(path);
+//   polygon_is_clockwise(poly);
 // Description:
 //   Return true if the given 2D simple polygon is in clockwise order, false otherwise.
 //   Results for complex (self-intersecting) polygon are indeterminate.
 // Arguments:
-//   path = The list of 2D path points for the perimeter of the polygon.
-function polygon_is_clockwise(path) =
-    assert(is_path(path) && len(path[0])==2, "Input must be a 2d path")
-    let(
-        minx = min(subindex(path,0)),
-        lowind = search(minx, path, 0, 0),
-        lowpts = select(path, lowind),
-        miny = min(subindex(lowpts, 1)),
-        extreme_sub = search(miny, lowpts, 1, 1)[0],
-        extreme = select(lowind,extreme_sub)
-    ) det2([select(path,extreme+1)-path[extreme], select(path, extreme-1)-path[extreme]])<0;
+//   poly = The list of 2D path points for the perimeter of the polygon.
+function polygon_is_clockwise(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d path")
+    polygon_area(poly, signed=true)<0;
 
 
 // Function: clockwise_polygon()
 // Usage:
-//   clockwise_polygon(path);
+//   clockwise_polygon(poly);
 // Description:
 //   Given a 2D polygon path, returns the clockwise winding version of that path.
-function clockwise_polygon(path) =
-    polygon_is_clockwise(path)? path : reverse_polygon(path);
+function clockwise_polygon(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d polygon")
+    polygon_area(poly, signed=true)<0 ? poly : reverse_polygon(poly);
 
 
 // Function: ccw_polygon()
 // Usage:
-//   ccw_polygon(path);
+//   ccw_polygon(poly);
 // Description:
-//   Given a 2D polygon path, returns the counter-clockwise winding version of that path.
-function ccw_polygon(path) =
-    polygon_is_clockwise(path)? reverse_polygon(path) : path;
+//   Given a 2D polygon poly, returns the counter-clockwise winding version of that poly.
+function ccw_polygon(poly) =
+    assert(is_path(poly,dim=2), "Input should be a 2d polygon")
+    polygon_area(poly, signed=true)<0 ? reverse_polygon(poly) : poly;
 
 
 // Function: reverse_polygon()
@@ -1656,6 +1805,7 @@ function ccw_polygon(path) =
 // Description:
 //   Reverses a polygon's winding direction, while still using the same start point.
 function reverse_polygon(poly) =
+    assert(is_path(poly), "Input should be a polygon")
     let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]];
 
 
@@ -1666,8 +1816,9 @@ function reverse_polygon(poly) =
 //   Given a 3D planar polygon, returns a unit-length normal vector for the
 //   clockwise orientation of the polygon. 
 function polygon_normal(poly) =
+    assert(is_path(poly,dim=3), "Invalid 3D polygon." )
     let(
-        poly = path3d(cleanup_path(poly)),
+        poly = cleanup_path(poly),
         p0 = poly[0],
         n = sum([
             for (i=[1:1:len(poly)-2])
@@ -1772,6 +1923,9 @@ function _split_polygon_at_z(poly, z) =
 //   polys = A list of 3D polygons to split.
 //   xs = A list of scalar X values to split at.
 function split_polygons_at_each_x(polys, xs, _i=0) =
+    assert( is_consistent(polys) && is_path(poly[0],dim=3) ,
+            "The input list should contains only 3D polygons." )
+    assert( is_finite(xs), "The split value list should contain only numbers." )  
     _i>=len(xs)? polys :
     split_polygons_at_each_x(
         [
@@ -1779,7 +1933,7 @@ function split_polygons_at_each_x(polys, xs, _i=0) =
             each _split_polygon_at_x(poly, xs[_i])
         ], xs, _i=_i+1
     );
-
+    
 
 // Function: split_polygons_at_each_y()
 // Usage:
@@ -1790,6 +1944,9 @@ function split_polygons_at_each_x(polys, xs, _i=0) =
 //   polys = A list of 3D polygons to split.
 //   ys = A list of scalar Y values to split at.
 function split_polygons_at_each_y(polys, ys, _i=0) =
+//    assert( is_consistent(polys) && is_path(polys[0],dim=3) , // not all polygons should have the same length!!!
+  //          "The input list should contains only 3D polygons." )
+    assert( is_finite(ys) || is_vector(ys), "The split value list should contain only numbers." )  //***
     _i>=len(ys)? polys :
     split_polygons_at_each_y(
         [
@@ -1808,6 +1965,9 @@ function split_polygons_at_each_y(polys, ys, _i=0) =
 //   polys = A list of 3D polygons to split.
 //   zs = A list of scalar Z values to split at.
 function split_polygons_at_each_z(polys, zs, _i=0) =
+    assert( is_consistent(polys) && is_path(poly[0],dim=3) ,
+            "The input list should contains only 3D polygons." )
+    assert( is_finite(zs), "The split value list should contain only numbers." )  
     _i>=len(zs)? polys :
     split_polygons_at_each_z(
         [
@@ -1817,5 +1977,4 @@ function split_polygons_at_each_z(polys, zs, _i=0) =
     );
 
 
-
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/hull.scad b/hull.scad
index da7d2c9..874e2c8 100644
--- a/hull.scad
+++ b/hull.scad
@@ -92,7 +92,7 @@ function hull2d_path(points) =
     assert(is_path(points,2),"Invalid input to hull2d_path")
     len(points) < 2 ? []
   : len(points) == 2 ? [0,1]
-  : let(tri=find_noncollinear_points(points, error=false))
+  : let(tri=noncollinear_triple(points, error=false))
     tri == [] ? _hull_collinear(points)
   : let(
         remaining = [ for (i = [0:1:len(points)-1]) if (i != tri[0] && i!=tri[1] && i!=tri[2]) i ],
@@ -170,7 +170,7 @@ function hull3d_faces(points) =
     assert(is_path(points,3),"Invalid input to hull3d_faces")
     len(points) < 3 ? list_range(len(points))
   : let ( // start with a single non-collinear triangle
-          tri = find_noncollinear_points(points, error=false)
+          tri = noncollinear_triple(points, error=false)
         )
     tri==[] ? _hull_collinear(points)
   : let(
@@ -250,7 +250,7 @@ function _find_conflicts(point, planes) = [
 
 
 function _find_first_noncoplanar(plane, points, i) = 
-    (i >= len(points) || !coplanar(plane, points[i]))? i :
+    (i >= len(points) || !points_on_plane([points[i]],plane))? i :
     _find_first_noncoplanar(plane, points, i+1);
 
 
diff --git a/math.scad b/math.scad
index 158b1f2..888d048 100644
--- a/math.scad
+++ b/math.scad
@@ -675,7 +675,7 @@ function convolve(p,q) =
 // Usage: linear_solve(A,b)
 // Description:
 //   Solves the linear system Ax=b.  If A is square and non-singular the unique solution is returned.  If A is overdetermined
-//   the least squares solution is returned.  If A is underdetermined, the minimal norm solution is returned.
+//   the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
 //   If A is rank deficient or singular then linear_solve returns [].  If b is a matrix that is compatible with A
 //   then the problem is solved for the matrix valued right hand side and a matrix is returned.  Note that if you 
 //   want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then
@@ -686,7 +686,7 @@ function linear_solve(A,b) =
         m = len(A),
         n = len(A[0])
     )
-    assert(is_vector(b,m) || is_matrix(b,m),"Incompatible matrix and right hand side")
+    assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
     let (
         qr = m<n? qr_factor(transpose(A)) : qr_factor(A),
         maxdim = max(n,m),
@@ -727,7 +727,7 @@ function qr_factor(A) =
       n = len(A[0])
     )
     let(
-        qr =_qr_factor(A, Q=ident(m), column=0, m = m, n=n),
+        qr = _qr_factor(A, Q=ident(m), column=0, m = m, n=n),
         Rzero = 
           let( R = qr[1] )
           [ for(i=[0:m-1]) [
@@ -745,7 +745,13 @@ function _qr_factor(A,Q, column, m, n) =
         u = x - concat([alpha],repeat(0,m-1)),
         v = alpha==0 ? u : u / norm(u),
         Qc = ident(len(x)) - 2*outer_product(v,v),
-        Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<column || j<column ? (i==j ? 1 : 0) : Qc[i-column][j-column]]]
+        Qf = [for(i=[0:m-1]) 
+                [for(j=[0:m-1]) 
+                    i<column || j<column 
+                    ? (i==j ? 1 : 0) 
+                    : Qc[i-column][j-column]
+                ]
+             ]
     )
     _qr_factor(Qf*A, Q*Qf, column+1, m, n);
 
diff --git a/rounding.scad b/rounding.scad
index b3bfa45..ae2aed0 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -1563,7 +1563,7 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
      // Determine which points are concave by making bottom 2d if necessary
      bot_proj = len(bottom[0])==2 ? bottom :  project_plane(bottom, select(bottom,0,2)),
      bottom_sign = polygon_is_clockwise(bot_proj) ? 1 : -1,
-     concave = [for(i=[0:N-1]) bottom_sign*sign(point_left_of_segment2d(select(bot_proj,i+1), select(bot_proj, i-1,i)))>0],
+     concave = [for(i=[0:N-1]) bottom_sign*sign(point_left_of_line2d(select(bot_proj,i+1), select(bot_proj, i-1,i)))>0],
      top = is_undef(top) ? path3d(bottom,height/2) :
            len(top[0])==2 ? path3d(top,height/2) :
            top,
@@ -1578,16 +1578,16 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
    assert(jsvecok || jssingleok,
           str("Argument joint_sides is invalid.  All entries must be nonnegative, and it must be a number, 2-vector, or a length ",N," list those."))
    assert(is_num(k_sides) || is_vector(k_sides,N), str("Curvature parameter k_sides must be a number or length ",N," vector"))
-   assert(points_are_coplanar(bottom))
-   assert(points_are_coplanar(top))
+   assert(coplanar(bottom))
+   assert(coplanar(top))
    assert(!is_num(k_sides) || (k_sides>=0 && k_sides<=1), "Curvature parameter k_sides must be in interval [0,1]")
    let(
-     non_coplanar=[for(i=[0:N-1]) if (!points_are_coplanar(concat(select(top,i,i+1), select(bottom,i,i+1)))) [i,(i+1)%N]],
+     non_coplanar=[for(i=[0:N-1]) if (!coplanar(concat(select(top,i,i+1), select(bottom,i,i+1)))) [i,(i+1)%N]],
      k_sides_vec = is_num(k_sides) ? repeat(k_sides, N) : k_sides,
      kbad = [for(i=[0:N-1]) if (k_sides_vec[i]<0 || k_sides_vec[i]>1) i],
      joint_sides_vec = jssingleok ? repeat(joint_sides,N) : joint_sides,
-     top_collinear = [for(i=[0:N-1]) if (points_are_collinear(select(top,i-1,i+1))) i],
-     bot_collinear = [for(i=[0:N-1]) if (points_are_collinear(select(bottom,i-1,i+1))) i]
+     top_collinear = [for(i=[0:N-1]) if (collinear(select(top,i-1,i+1))) i],
+     bot_collinear = [for(i=[0:N-1]) if (collinear(select(bottom,i-1,i+1))) i]
    )
    assert(non_coplanar==[], str("Side faces are non-coplanar at edges: ",non_coplanar))
    assert(top_collinear==[], str("Top has collinear or duplicated points at indices: ",top_collinear))
@@ -1653,14 +1653,14 @@ function rounded_prism(bottom, top, joint_bot, joint_top, joint_sides, k_bot, k_
                vline = concat(select(subindex(top_patch[i],j),2,4),
                               select(subindex(bot_patch[i],j),2,4))
              )
-         if (!points_are_collinear(vline)) [i,j]],
+         if (!collinear(vline)) [i,j]],
      //verify horiz edges
      verify_horiz=[for(i=[0:N-1], j=[0:4])
          let(
              hline_top = concat(select(top_patch[i][j],2,4), select(select(top_patch, i+1)[j],0,2)),
              hline_bot = concat(select(bot_patch[i][j],2,4), select(select(bot_patch, i+1)[j],0,2))
          )
-         if (!points_are_collinear(hline_top) || !points_are_collinear(hline_bot)) [i,j]]
+         if (!collinear(hline_top) || !collinear(hline_bot)) [i,j]]
     )
     assert(debug || top_intersections==[],
           "Roundovers interfere with each other on top face: either input is self intersecting or top joint length is too large")
@@ -1994,4 +1994,4 @@ module bent_cutout_mask(r, thickness, path, convexity=10)
 }
 
 
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
\ No newline at end of file
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 8f780f4..0950e1c 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -1,6 +1,146 @@
 include <../std.scad>
 
 
+//the commented lines are for tests to be written
+//the tests are ordered as they appear in geometry.scad
+
+test_point_on_segment2d();
+test_point_left_of_line2d();
+test_collinear();
+test_distance_from_line();
+test_line_normal();
+test_line_intersection();
+//test_line_ray_intersection();
+test_line_segment_intersection();
+//test_ray_intersection();
+//test_ray_segment_intersection();
+test_segment_intersection();
+test_line_closest_point();
+//test_ray_closest_point();
+test_segment_closest_point();
+test_line_from_points();
+test_tri_calc();
+//test_hyp_opp_to_adj();
+//test_hyp_ang_to_adj();
+//test_opp_ang_to_adj();
+//test_hyp_adj_to_opp();
+//test_hyp_ang_to_opp();
+//test_adj_ang_to_opp();
+//test_adj_opp_to_hyp();
+//test_adj_ang_to_hyp();
+//test_opp_ang_to_hyp();
+//test_hyp_adj_to_ang();
+//test_hyp_opp_to_ang();
+//test_adj_opp_to_ang();
+test_triangle_area();
+test_plane3pt();
+test_plane3pt_indexed();
+//test_plane_from_normal();
+test_plane_from_points();
+//test_plane_from_polygon();
+test_plane_normal();
+//test_plane_offset();
+//test_plane_transform();
+test_projection_on_plane();
+//test_plane_point_nearest_origin();
+test_distance_from_plane();
+
+test_find_circle_2tangents();
+test_find_circle_3points();
+test_circle_point_tangents();
+test_tri_functions();
+//test_closest_point_on_plane();
+//test__general_plane_line_intersection();
+//test_plane_line_angle();
+//test_plane_line_intersection();
+//test_polygon_line_intersection();
+//test_plane_intersection();
+test_coplanar();
+test_points_on_plane();
+test_in_front_of_plane();
+//test_find_circle_2tangents();
+//test_find_circle_3points();
+//test_circle_point_tangents();
+//test_circle_circle_tangents();
+test_noncollinear_triple();
+test_pointlist_bounds();
+test_closest_point();
+test_furthest_point();
+test_polygon_area();
+test_is_convex_polygon();
+test_polygon_shift();
+test_polygon_shift_to_closest_point();
+test_reindex_polygon();
+test_align_polygon();
+test_centroid();
+test_point_in_polygon();
+test_polygon_is_clockwise();
+test_clockwise_polygon();
+test_ccw_polygon();
+test_reverse_polygon();
+//test_polygon_normal();
+//test_split_polygons_at_each_x();
+//test_split_polygons_at_each_y();
+//test_split_polygons_at_each_z();
+
+//tests to migrate to other files
+test_is_path();
+test_is_closed_path();
+test_close_path();
+test_cleanup_path();
+test_simplify_path();
+test_simplify_path_indexed();
+test_is_region();
+
+// to be used when there are two alternative symmetrical outcomes 
+// from a function like a plane output.
+function standardize(v) = 
+    v==[]? [] :
+    sign([for(vi=v) if( ! approx(vi,0)) vi,0 ][0])*v;
+
+module assert_std(vc,ve) { assert(standardize(vc)==standardize(ve)); }
+
+module test_points_on_plane() {
+    pts     = [for(i=[0:40]) rands(-1,1,3) ];
+    dir     = rands(-10,10,3);
+    normal0 = unit([1,2,3]);
+    ang     = rands(0,360,1)[0];
+    normal  = rot(a=ang,p=normal0);
+    plane   = [each normal, normal*dir];
+    prj_pts = projection_on_plane(plane,pts);
+    assert(points_on_plane(prj_pts,plane));
+    assert(!points_on_plane(concat(pts,[normal-dir]),plane));
+}
+*test_points_on_plane();
+
+module test_projection_on_plane(){
+    ang     = rands(0,360,1)[0];
+    dir     = rands(-10,10,3);
+    normal0 = unit([1,2,3]);
+    normal  = rot(a=ang,p=normal0);
+    plane0  = [each normal0, 0];
+    plane   = [each normal,  0];
+    planem  = [each normal, normal*dir];
+    pts     = [for(i=[1:10]) rands(-1,1,3)];
+    assert_approx( projection_on_plane(plane,pts),
+                   projection_on_plane(plane,projection_on_plane(plane,pts)));
+    assert_approx( projection_on_plane(plane,pts),
+                   rot(a=ang,p=projection_on_plane(plane0,rot(a=-ang,p=pts))));    
+    assert_approx( move((-normal*dir)*normal,p=projection_on_plane(planem,pts)),
+                   projection_on_plane(plane,pts));
+    assert_approx( move((normal*dir)*normal,p=projection_on_plane(plane,pts)),
+                   projection_on_plane(planem,pts));
+}
+*test_projection_on_plane();
+
+module test_line_from_points() {
+    assert_approx(line_from_points([[1,0],[0,0],[-1,0]]),[[-1,0],[1,0]]);
+    assert_approx(line_from_points([[1,1],[0,1],[-1,1]]),[[-1,1],[1,1]]);
+    assert(line_from_points([[1,1],[0,1],[-1,0]])==undef);
+    assert(line_from_points([[1,1],[0,1],[-1,0]],fast=true)== [[-1,0],[1,1]]);
+}
+*test_line_from_points();
+
 module test_point_on_segment2d() {
     assert(point_on_segment2d([-15,0], [[-10,0], [10,0]]) == false);
     assert(point_on_segment2d([-10,0], [[-10,0], [10,0]]) == true);
@@ -29,42 +169,28 @@ module test_point_on_segment2d() {
     assert(point_on_segment2d([ 10, 10], [[-10,-10], [10,10]]) == true);
     assert(point_on_segment2d([ 15, 15], [[-10,-10], [10,10]]) == false);
 }
-test_point_on_segment2d();
+*test_point_on_segment2d();
 
 
-module test_point_left_of_segment() {
-    assert(point_left_of_segment2d([ -3,  0], [[-10,-10], [10,10]]) > 0);
-    assert(point_left_of_segment2d([  0,  0], [[-10,-10], [10,10]]) == 0);
-    assert(point_left_of_segment2d([  3,  0], [[-10,-10], [10,10]]) < 0);
+module test_point_left_of_line2d() {
+    assert(point_left_of_line2d([ -3,  0], [[-10,-10], [10,10]]) > 0);
+    assert(point_left_of_line2d([  0,  0], [[-10,-10], [10,10]]) == 0);
+    assert(point_left_of_line2d([  3,  0], [[-10,-10], [10,10]]) < 0);
 }
-test_point_left_of_segment();
-
+*test_point_left_of_line2d();
 
 module test_collinear() {
     assert(collinear([-10,-10], [-15, -16], [10,10]) == false);
+    assert(collinear([[-10,-10], [-15, -16], [10,10]]) == false);
     assert(collinear([-10,-10], [-15, -15], [10,10]) == true);
+    assert(collinear([[-10,-10], [-15, -15], [10,10]]) == true);
     assert(collinear([-10,-10], [ -3,   0], [10,10]) == false);
     assert(collinear([-10,-10], [  0,   0], [10,10]) == true);
     assert(collinear([-10,-10], [  3,   0], [10,10]) == false);
     assert(collinear([-10,-10], [ 15,  15], [10,10]) == true);
     assert(collinear([-10,-10], [ 15,  16], [10,10]) == false);
 }
-test_collinear();
-
-
-module test_collinear_indexed() {
-    pts = [
-        [-20,-20], [-10,-20], [0,-10], [10,0], [20,10], [20,20], [15,30]
-    ];
-    assert(collinear_indexed(pts, 0,1,2) == false);
-    assert(collinear_indexed(pts, 1,2,3) == true);
-    assert(collinear_indexed(pts, 2,3,4) == true);
-    assert(collinear_indexed(pts, 3,4,5) == false);
-    assert(collinear_indexed(pts, 4,5,6) == false);
-    assert(collinear_indexed(pts, 4,3,2) == true);
-    assert(collinear_indexed(pts, 0,5,6) == false);
-}
-test_collinear_indexed();
+*test_collinear();
 
 
 module test_distance_from_line() {
@@ -73,7 +199,7 @@ module test_distance_from_line() {
     assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [1,-1,0]) - sqrt(2)) < EPSILON);
     assert(abs(distance_from_line([[-10,-10,-10], [10,10,10]], [8,-8,0]) - 8*sqrt(2)) < EPSILON);
 }
-test_distance_from_line();
+*test_distance_from_line();
 
 
 module test_line_normal() {
@@ -97,7 +223,7 @@ module test_line_normal() {
         assert(approx(n2, n1));
     }
 }
-test_line_normal();
+*test_line_normal();
 
 
 module test_line_intersection() {
@@ -110,7 +236,7 @@ module test_line_intersection() {
     assert(line_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [-10, 10]]) == [0,0]);
     assert(line_intersection([[ -8,  0], [ 12,  4]], [[ 12,  0], [ -8,  4]]) == [2,2]);
 }
-test_line_intersection();
+*test_line_intersection();
 
 
 module test_segment_intersection() {
@@ -126,7 +252,7 @@ module test_segment_intersection() {
     assert(segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [-10, 10]]) == [0,0]);
     assert(segment_intersection([[ -8,  0], [ 12,  4]], [[ 12,  0], [ -8,  4]]) == [2,2]);
 }
-test_segment_intersection();
+*test_segment_intersection();
 
 
 module test_line_segment_intersection() {
@@ -141,7 +267,7 @@ module test_line_segment_intersection() {
     assert(line_segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [  1, -1]]) == undef);
     assert(line_segment_intersection([[-10,-10], [ 10, 10]], [[ 10,-10], [ -1,  1]]) == [0,0]);
 }
-test_line_segment_intersection();
+*test_line_segment_intersection();
 
 
 module test_line_closest_point() {
@@ -151,7 +277,7 @@ module test_line_closest_point() {
     assert(approx(line_closest_point([[-10,-20], [10,20]], [1,2]+[2,-1]), [1,2]));
     assert(approx(line_closest_point([[-10,-20], [10,20]], [13,31]), [15,30]));
 }
-test_line_closest_point();
+*test_line_closest_point();
 
 
 module test_segment_closest_point() {
@@ -162,10 +288,10 @@ module test_segment_closest_point() {
     assert(approx(segment_closest_point([[-10,-20], [10,20]], [13,31]), [10,20]));
     assert(approx(segment_closest_point([[-10,-20], [10,20]], [15,25]), [10,20]));
 }
-test_segment_closest_point();
-
+*test_segment_closest_point();
 
 module test_find_circle_2tangents() {
+//** missing tests with arg tangent=true
     assert(approx(find_circle_2tangents([10,10],[0,0],[10,-10],r=10/sqrt(2))[0],[10,0]));
     assert(approx(find_circle_2tangents([-10,10],[0,0],[-10,-10],r=10/sqrt(2))[0],[-10,0]));
     assert(approx(find_circle_2tangents([-10,10],[0,0],[10,10],r=10/sqrt(2))[0],[0,10]));
@@ -174,9 +300,9 @@ module test_find_circle_2tangents() {
     assert(approx(find_circle_2tangents([10,0],[0,0],[0,-10],r=10)[0],[10,-10]));
     assert(approx(find_circle_2tangents([0,-10],[0,0],[-10,0],r=10)[0],[-10,-10]));
     assert(approx(find_circle_2tangents([-10,0],[0,0],[0,10],r=10)[0],[-10,10]));
-    assert(approx(find_circle_2tangents(polar_to_xy(10,60),[0,0],[10,0],r=10)[0],polar_to_xy(20,30)));
+    assert_approx(find_circle_2tangents(polar_to_xy(10,60),[0,0],[10,0],r=10)[0],polar_to_xy(20,30));
 }
-test_find_circle_2tangents();
+*test_find_circle_2tangents();
 
 
 module test_find_circle_3points() {
@@ -291,7 +417,7 @@ module test_find_circle_3points() {
         }
     }
 }
-test_find_circle_3points();
+*test_find_circle_3points();
 
 
 module test_circle_point_tangents() {
@@ -304,7 +430,7 @@ module test_circle_point_tangents() {
         assert(approx(flatten(got), flatten(expected)));
     }
 }
-test_circle_point_tangents();
+*test_circle_point_tangents();
 
 
 module test_tri_calc() {
@@ -327,23 +453,9 @@ module test_tri_calc() {
         assert(approx(tri_calc(hyp=hyp, ang2=ang2), expected));
     }
 }
-test_tri_calc();
+*test_tri_calc();
 
 
-// Dummy modules to show up in coverage check script.
-module test_hyp_opp_to_adj();
-module test_hyp_ang_to_adj();
-module test_opp_ang_to_adj();
-module test_hyp_adj_to_opp();
-module test_hyp_ang_to_opp();
-module test_adj_ang_to_opp();
-module test_adj_opp_to_hyp();
-module test_adj_ang_to_hyp();
-module test_opp_ang_to_hyp();
-module test_hyp_adj_to_ang();
-module test_hyp_opp_to_ang();
-module test_adj_opp_to_ang();
-
 module test_tri_functions() {
     sides = rands(1,100,100,seed_value=8181);
     for (p = pair_wrap(sides)) {
@@ -365,7 +477,7 @@ module test_tri_functions() {
         assert_approx(adj_opp_to_ang(adj,opp), ang);
     }
 }
-test_tri_functions();
+*test_tri_functions();
 
 
 module test_triangle_area() {
@@ -373,55 +485,53 @@ module test_triangle_area() {
     assert(abs(triangle_area([0,0], [0,10], [0,15])) < EPSILON);
     assert(abs(triangle_area([0,0], [10,0], [0,10]) - 50) < EPSILON);
 }
-test_triangle_area();
+*test_triangle_area();
 
 
 module test_plane3pt() {
-    assert(plane3pt([0,0,20], [0,10,10], [0,0,0]) == [1,0,0,0]);
-    assert(plane3pt([2,0,20], [2,10,10], [2,0,0]) == [1,0,0,2]);
-    assert(plane3pt([0,0,0], [10,0,10], [0,0,20]) == [0,1,0,0]);
-    assert(plane3pt([0,2,0], [10,2,10], [0,2,20]) == [0,1,0,2]);
-    assert(plane3pt([0,0,0], [10,10,0], [20,0,0]) == [0,0,1,0]);
-    assert(plane3pt([0,0,2], [10,10,2], [20,0,2]) == [0,0,1,2]);
+    assert_std(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]);
+    assert_std(plane3pt([2,0,20], [2,10,10], [2,0,0]), [1,0,0,2]);
+    assert_std(plane3pt([0,0,0], [10,0,10], [0,0,20]), [0,1,0,0]);
+    assert_std(plane3pt([0,2,0], [10,2,10], [0,2,20]), [0,1,0,2]);
+    assert_std(plane3pt([0,0,0], [10,10,0], [20,0,0]), [0,0,1,0]);
+    assert_std(plane3pt([0,0,2], [10,10,2], [20,0,2]), [0,0,1,2]);
 }
-test_plane3pt();
-
+*test_plane3pt();
 
 module test_plane3pt_indexed() {
     pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ];
     s13 = sqrt(1/3);
-    assert(plane3pt_indexed(pts, 0,3,2) == [1,0,0,0]);
-    assert(plane3pt_indexed(pts, 0,2,3) == [-1,0,0,0]);
-    assert(plane3pt_indexed(pts, 0,1,3) == [0,1,0,0]);
-    assert(plane3pt_indexed(pts, 0,3,1) == [0,-1,0,0]);
-    assert(plane3pt_indexed(pts, 0,2,1) == [0,0,1,0]);
-    assert(plane3pt_indexed(pts, 0,1,2) == [0,0,-1,0]);
-    assert(plane3pt_indexed(pts, 3,2,1) == [s13,s13,s13,10*s13]);
-    assert(plane3pt_indexed(pts, 1,2,3) == [-s13,-s13,-s13,-10*s13]);
+    assert_std(plane3pt_indexed(pts, 0,3,2), [1,0,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]);
+    assert_std(plane3pt_indexed(pts, 0,2,1), [0,0,1,0]);
+    assert_approx(plane3pt_indexed(pts, 0,1,2), [0,0,-1,0]);
+    assert_approx(plane3pt_indexed(pts, 3,2,1), [s13,s13,s13,10*s13]);
+    assert_approx(plane3pt_indexed(pts, 1,2,3), [-s13,-s13,-s13,-10*s13]);
 }
-test_plane3pt_indexed();
-
+*test_plane3pt_indexed();
 
 module test_plane_from_points() {
-    assert(plane_from_points([[0,0,20], [0,10,10], [0,0,0], [0,5,3]]) == [1,0,0,0]);
-    assert(plane_from_points([[2,0,20], [2,10,10], [2,0,0], [2,3,4]]) == [1,0,0,2]);
-    assert(plane_from_points([[0,0,0], [10,0,10], [0,0,20], [5,0,7]]) == [0,1,0,0]);
-    assert(plane_from_points([[0,2,0], [10,2,10], [0,2,20], [4,2,3]]) == [0,1,0,2]);
-    assert(plane_from_points([[0,0,0], [10,10,0], [20,0,0], [8,3,0]]) == [0,0,1,0]);
-    assert(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]) == [0,0,1,2]);
+    assert_std(plane_from_points([[0,0,20], [0,10,10], [0,0,0], [0,5,3]]), [1,0,0,0]);
+    assert_std(plane_from_points([[2,0,20], [2,10,10], [2,0,0], [2,3,4]]), [1,0,0,2]);
+    assert_std(plane_from_points([[0,0,0], [10,0,10], [0,0,20], [5,0,7]]), [0,1,0,0]);
+    assert_std(plane_from_points([[0,2,0], [10,2,10], [0,2,20], [4,2,3]]), [0,1,0,2]);
+    assert_std(plane_from_points([[0,0,0], [10,10,0], [20,0,0], [8,3,0]]), [0,0,1,0]);
+    assert_std(plane_from_points([[0,0,2], [10,10,2], [20,0,2], [3,4,2]]), [0,0,1,2]);
 }
-test_plane_from_points();
+*test_plane_from_points();
 
 
 module test_plane_normal() {
-    assert(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])) == [1,0,0]);
-    assert(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])) == [1,0,0]);
-    assert(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])) == [0,1,0]);
-    assert(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])) == [0,1,0]);
-    assert(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])) == [0,0,1]);
-    assert(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])) == [0,0,1]);
+    assert_std(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]);
+    assert_std(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])), [1,0,0]);
+    assert_std(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])), [0,1,0]);
+    assert_std(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])), [0,1,0]);
+    assert_std(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])), [0,0,1]);
+    assert_std(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])), [0,0,1]);
 }
-test_plane_normal();
+*test_plane_normal();
 
 
 module test_distance_from_plane() {
@@ -429,20 +539,16 @@ module test_distance_from_plane() {
     assert(distance_from_plane(plane1, [0,0,5]) == 5);
     assert(distance_from_plane(plane1, [5,5,8]) == 8);
 }
-test_distance_from_plane();
+*test_distance_from_plane();
 
 
 module test_coplanar() {
-    plane = plane3pt([0,0,0], [0,10,10], [10,0,10]);
-    assert(coplanar(plane, [5,5,10]) == true);
-    assert(coplanar(plane, [10/3,10/3,20/3]) == true);
-    assert(coplanar(plane, [0,0,0]) == true);
-    assert(coplanar(plane, [1,1,0]) == false);
-    assert(coplanar(plane, [-1,1,0]) == true);
-    assert(coplanar(plane, [1,-1,0]) == true);
-    assert(coplanar(plane, [5,5,5]) == false);
+    assert(coplanar([ [5,5,1],[0,0,1],[-1,-1,1] ]) == false);
+    assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true);
+    assert(coplanar([ [0,0,0],[1,0,1],[1,1,1], [0,1,2] ]) == false);
+    assert(coplanar([ [0,0,0],[1,0,1],[1,1,2], [0,1,1] ]) == true);
 }
-test_coplanar();
+*test_coplanar();
 
 
 module test_in_front_of_plane() {
@@ -455,7 +561,7 @@ module test_in_front_of_plane() {
     assert(in_front_of_plane(plane, [0,0,5]) == true);
     assert(in_front_of_plane(plane, [0,0,-5]) == false);
 }
-test_in_front_of_plane();
+*test_in_front_of_plane();
 
 
 module test_is_path() {
@@ -470,35 +576,43 @@ module test_is_path() {
     assert(is_path([[1,2,3],[4,5,6]]));
     assert(is_path([[1,2,3],[4,5,6],[7,8,9]]));
 }
-test_is_path();
+*test_is_path();
 
 
 module test_is_closed_path() {
     assert(!is_closed_path([[1,2,3],[4,5,6],[1,8,9]]));
     assert(is_closed_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]));
 }
-test_is_closed_path();
+*test_is_closed_path();
 
 
 module test_close_path() {
     assert(close_path([[1,2,3],[4,5,6],[1,8,9]]) == [[1,2,3],[4,5,6],[1,8,9],[1,2,3]]);
     assert(close_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]) == [[1,2,3],[4,5,6],[1,8,9],[1,2,3]]);
 }
-test_close_path();
+*test_close_path();
 
 
 module test_cleanup_path() {
     assert(cleanup_path([[1,2,3],[4,5,6],[1,8,9]]) == [[1,2,3],[4,5,6],[1,8,9]]);
     assert(cleanup_path([[1,2,3],[4,5,6],[1,8,9],[1,2,3]]) == [[1,2,3],[4,5,6],[1,8,9]]);
 }
-test_cleanup_path();
+*test_cleanup_path();
 
 
 module test_polygon_area() {
     assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4));
     assert(approx(polygon_area(circle(r=50,$fn=1000)), -PI*50*50, eps=0.1));
 }
-test_polygon_area();
+*test_polygon_area();
+
+
+module test_is_convex_polygon() {
+    assert(is_convex_polygon([[1,1],[-1,1],[-1,-1],[1,-1]]));
+    assert(is_convex_polygon(circle(r=50,$fn=1000)));
+    assert(!is_convex_polygon([[1,1],[0,0],[-1,1],[-1,-1],[1,-1]]));
+}
+*test_is_convex_polygon();
 
 
 module test_polygon_shift() {
@@ -506,7 +620,7 @@ module test_polygon_shift() {
     assert(polygon_shift(path,1) == [[-1,1],[-1,-1],[1,-1],[1,1]]);
     assert(polygon_shift(path,2) == [[-1,-1],[1,-1],[1,1],[-1,1]]);
 }
-test_polygon_shift();
+*test_polygon_shift();
 
 
 module test_polygon_shift_to_closest_point() {
@@ -516,56 +630,45 @@ module test_polygon_shift_to_closest_point() {
     assert(polygon_shift_to_closest_point(path,[-1.1,-1.1]) == [[-1,-1],[1,-1],[1,1],[-1,1]]);
     assert(polygon_shift_to_closest_point(path,[1.1,-1.1]) == [[1,-1],[1,1],[-1,1],[-1,-1]]);
 }
-test_polygon_shift_to_closest_point();
+*test_polygon_shift_to_closest_point();
 
 
-/*
-module test_first_noncollinear(){
-    pts = [
-        [1,1], [2,2], [3,3], [4,4], [4,5], [5,6]
-    ];
-    assert(first_noncollinear(0,1,pts) == 4);
-    assert(first_noncollinear(1,0,pts) == 4);
-    assert(first_noncollinear(0,2,pts) == 4);
-    assert(first_noncollinear(2,0,pts) == 4);
-    assert(first_noncollinear(1,2,pts) == 4);
-    assert(first_noncollinear(2,1,pts) == 4);
-    assert(first_noncollinear(0,3,pts) == 4);
-    assert(first_noncollinear(3,0,pts) == 4);
-    assert(first_noncollinear(1,3,pts) == 4);
-    assert(first_noncollinear(3,1,pts) == 4);
-    assert(first_noncollinear(2,3,pts) == 4);
-    assert(first_noncollinear(3,2,pts) == 4);
-    assert(first_noncollinear(0,4,pts) == 1);
-    assert(first_noncollinear(4,0,pts) == 1);
-    assert(first_noncollinear(1,4,pts) == 0);
-    assert(first_noncollinear(4,1,pts) == 0);
-    assert(first_noncollinear(2,4,pts) == 0);
-    assert(first_noncollinear(4,2,pts) == 0);
-    assert(first_noncollinear(3,4,pts) == 0);
-    assert(first_noncollinear(4,3,pts) == 0);
-    assert(first_noncollinear(0,5,pts) == 1);
-    assert(first_noncollinear(5,0,pts) == 1);
-    assert(first_noncollinear(1,5,pts) == 0);
-    assert(first_noncollinear(5,1,pts) == 0);
-    assert(first_noncollinear(2,5,pts) == 0);
-    assert(first_noncollinear(5,2,pts) == 0);
-    assert(first_noncollinear(3,5,pts) == 0);
-    assert(first_noncollinear(5,3,pts) == 0);
-    assert(first_noncollinear(4,5,pts) == 0);
-    assert(first_noncollinear(5,4,pts) == 0);
+module test_reindex_polygon() {
+   pent = subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],5);
+   circ = circle($fn=5,r=2.2);
+   assert_approx(reindex_polygon(circ,pent), [[0.951056516295,0.309016994375],[0.587785252292,-0.809016994375],[-0.587785252292,-0.809016994375],[-0.951056516295,0.309016994375],[0,1]]);
+   poly = [[-1,1],[-1,-1],[1,-1],[1,1],[0,0]];
+   ref  = [for(i=[0:4])[sin(72*i),cos(72*i)]];
+   assert_approx(reindex_polygon(ref,poly),[[0,0],[1,1],[1,-1],[-1,-1],[-1,1]]);
 }
-test_first_noncollinear();
-*/
+*test_reindex_polygon();
 
 
-module test_find_noncollinear_points() {
-    assert(find_noncollinear_points([[1,1],[2,2],[3,3],[4,4],[4,5],[5,6]]) == [0,5,3]);
-    assert(find_noncollinear_points([[1,1],[2,2],[8,3],[4,4],[4,5],[5,6]]) == [0,2,5]);
+module test_align_polygon() {
+   pentagon = subdivide_path(pentagon(side=2),10);
+   hexagon  = subdivide_path(hexagon(side=2.7),10);
+   aligned =  [[2.7,0],[2.025,-1.16913429511],[1.35,-2.33826859022],
+               [-1.35,-2.33826859022],[-2.025,-1.16913429511],[-2.7,0],
+               [-2.025,1.16913429511],[-1.35,2.33826859022],[1.35,2.33826859022],
+               [2.025,1.16913429511]];
+   assert_approx(align_polygon(pentagon,hexagon,[0:10:359]), aligned);
+   aligned2 = [[1.37638192047,0],[1.37638192047,-1],[0.425325404176,-1.30901699437],
+               [-0.525731112119,-1.61803398875],[-1.11351636441,-0.809016994375],
+               [-1.7013016167,0],[-1.11351636441,0.809016994375],
+               [-0.525731112119,1.61803398875],[0.425325404176,1.30901699437],
+               [1.37638192047,1]];
+   assert_approx(align_polygon(hexagon,pentagon,[0:10:359]), aligned2);
+}
+*test_align_polygon();
+
+
+module test_noncollinear_triple() {
+    assert(noncollinear_triple([[1,1],[2,2],[3,3],[4,4],[4,5],[5,6]]) == [0,5,3]);
+    assert(noncollinear_triple([[1,1],[2,2],[8,3],[4,4],[4,5],[5,6]]) == [0,2,5]);
     u = unit([5,3]);
-    assert_equal(find_noncollinear_points([for(i = [2,3,4,5,7,12,15]) i * u], error=false),[]);
+    assert_equal(noncollinear_triple([for(i = [2,3,4,5,7,12,15]) i * u], error=false),[]);
 }
-test_find_noncollinear_points();
+*test_noncollinear_triple();
 
 
 module test_centroid() {
@@ -573,15 +676,18 @@ module test_centroid() {
     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 = [for(a=[0:90:360]) 
+              move([1,2.5,3.1], rot(p=[cos(a),sin(a),0],from=[0,0,1],to=[1,1,1])) ];
+    assert_approx(centroid(poly), [1,2.5,3.1]);
 }
-test_centroid();
+*test_centroid();
 
 
 module test_simplify_path() {
     path = [[-20,-20], [-10,-20], [0,-10], [10,0], [20,10], [20,20], [15,30]];
     assert(simplify_path(path) == [[-20,-20], [-10,-20], [20,10], [20,20], [15,30]]);
 }
-test_simplify_path();
+*test_simplify_path();
 
 
 module test_simplify_path_indexed() {
@@ -589,23 +695,29 @@ module test_simplify_path_indexed() {
     path = [4,6,1,0,3,2,5];
     assert(simplify_path_indexed(pts, path) == [4,6,3,2,5]);
 }
-test_simplify_path_indexed();
+*test_simplify_path_indexed();
 
 
 module test_point_in_polygon() {
     poly = [for (a=[0:30:359]) 10*[cos(a),sin(a)]];
+    poly2 = [ [-3,-3],[2,-3],[2,1],[-1,1],[-1,-1],[1,-1],[1,2],[-3,2] ];
     assert(point_in_polygon([0,0], poly) == 1);
     assert(point_in_polygon([20,0], poly) == -1);
+    assert(point_in_polygon([20,0], poly,EPSILON,nonzero=false) == -1);
     assert(point_in_polygon([5,5], poly) == 1);
     assert(point_in_polygon([-5,5], poly) == 1);
     assert(point_in_polygon([-5,-5], poly) == 1);
     assert(point_in_polygon([5,-5], poly) == 1);
+    assert(point_in_polygon([5,-5], poly,EPSILON,nonzero=false) == 1);
     assert(point_in_polygon([-10,-10], poly) == -1);
     assert(point_in_polygon([10,0], poly) == 0);
     assert(point_in_polygon([0,10], poly) == 0);
     assert(point_in_polygon([0,-10], poly) == 0);
+    assert(point_in_polygon([0,-10], poly,EPSILON,nonzero=false) == 0);
+    assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=true) == 1);
+    assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=false) == -1);
 }
-test_point_in_polygon();
+*test_point_in_polygon();
 
 
 module test_pointlist_bounds() {
@@ -626,16 +738,16 @@ module test_pointlist_bounds() {
     ];
     assert(pointlist_bounds(pts2d) == [[-63,-42],[84,42]]);
     pts5d = [
-        [-53,27,12,-53,12],
-        [-63,97,36,-63,36],
-        [84,-32,-5,84,-5], 
-        [63,-24,42,63,42], 
-        [23,57,-42,23,-42]
+        [-53, 27, 12,-53, 12],
+        [-63, 97, 36,-63, 36],
+        [ 84,-32, -5, 84, -5], 
+        [ 63,-24, 42, 63, 42], 
+        [ 23, 57,-42, 23,-42]
     ];
     assert(pointlist_bounds(pts5d) == [[-63,-32,-42,-63,-42],[84,97,42,84,42]]);
     assert(pointlist_bounds([[3,4,5,6]]), [[3,4,5,6],[3,4,5,6]]);
 }
-test_pointlist_bounds();
+*test_pointlist_bounds();
 
 
 module test_closest_point() {
@@ -648,7 +760,7 @@ module test_closest_point() {
         assert(mindist == dists[pidx]);
     }
 }
-test_closest_point();
+*test_closest_point();
 
 
 module test_furthest_point() {
@@ -661,7 +773,7 @@ module test_furthest_point() {
         assert(mindist == dists[pidx]);
     }
 }
-test_furthest_point();
+*test_furthest_point();
 
 
 module test_polygon_is_clockwise() {
@@ -670,7 +782,7 @@ module test_polygon_is_clockwise() {
     assert(polygon_is_clockwise(circle(d=100)));
     assert(polygon_is_clockwise(square(100)));
 }
-test_polygon_is_clockwise();
+*test_polygon_is_clockwise();
 
 
 module test_clockwise_polygon() {
@@ -679,7 +791,7 @@ module test_clockwise_polygon() {
     assert(clockwise_polygon(path) == path);
     assert(clockwise_polygon(rpath) == path);
 }
-test_clockwise_polygon();
+*test_clockwise_polygon();
 
 
 module test_ccw_polygon() {
@@ -688,7 +800,7 @@ module test_ccw_polygon() {
     assert(ccw_polygon(path) == rpath);
     assert(ccw_polygon(rpath) == rpath);
 }
-test_ccw_polygon();
+*test_ccw_polygon();
 
 
 module test_reverse_polygon() {
@@ -697,7 +809,7 @@ module test_reverse_polygon() {
     assert(reverse_polygon(path) == rpath);
     assert(reverse_polygon(rpath) == path);
 }
-test_reverse_polygon();
+*test_reverse_polygon();
 
 
 module test_is_region() {
@@ -709,7 +821,7 @@ module test_is_region() {
     assert(!is_region(true));
     assert(!is_region("foo"));
 }
-test_is_region();
+*test_is_region();
 
 
 
diff --git a/vnf.scad b/vnf.scad
index b72902e..05fd82f 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -403,6 +403,16 @@ function _triangulate_planar_convex_polygons(polys) =
         outtris = concat(tris, newtris, newtris2)
     ) outtris;
 
+//**
+// this function may produce degenerate triangles:
+//    _triangulate_planar_convex_polygons([ [for(i=[0:1]) [i,i],
+//                                           [1,-1], [-1,-1],
+//                                           for(i=[-1:0]) [i,i] ] ] )
+//    == [[[-1, -1], [ 0,  0], [0,  0]]
+//        [[-1, -1], [-1, -1], [0,  0]]
+//        [[ 1, -1], [-1, -1], [0,  0]]
+//        [[ 0,  0], [ 1,  1], [1, -1]] ]
+//
 
 // Function: vnf_bend()
 // Usage:
@@ -647,7 +657,7 @@ function vnf_validate(vnf, show_warns=true, check_isects=false) =
         nonplanars = unique([
             for (face = faces) let(
                 faceverts = [for (k=face) varr[k]]
-            ) if (!points_are_coplanar(faceverts)) [
+            ) if (!coplanar(faceverts)) [
                 "ERROR",
                 "NONPLANAR",
                 "Face vertices are not coplanar",