diff --git a/geometry.scad b/geometry.scad
index 21b4f29..80f6756 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -888,11 +888,49 @@ 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) / norm(normal);
+    assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0),
+            "Inputs `normal` and `pt` should be 3d vectors/points and `normal` cannot be zero." )
+    concat(normal, normal*pt) / norm(normal);
 
 
+// Eigenvalues for a 3x3 symmetrical matrix in decreasing order
+// Based on: https://en.wikipedia.org/wiki/Eigenvalue_algorithm
+function _eigenvals_symm_3(M) =
+  let( p1 = pow(M[0][1],2) + pow(M[0][2],2) + pow(M[1][2],2) )
+  (p1<EPSILON)  
+  ? -sort(-[ M[0][0], M[1][1], M[2][2] ]) //  diagonal matrix: eigenvals in decreasing order
+  : let(  q  = (M[0][0]+M[1][1]+M[2][2])/3,
+          B  = (M - q*ident(3)),
+          dB = [B[0][0], B[1][1], B[2][2]],
+          p2 = dB*dB + 2*p1,
+          p  = sqrt(p2/6),
+          r  = det3(B/p)/2,
+          ph = acos(constrain(r,-1,1))/3,
+          e1 = q + 2*p*cos(ph),
+          e3 = q + 2*p*cos(ph+120),
+          e2 = 3*q - e1 - e3 ) 
+    [ e1, e2, e3 ];
+
+
+// i-th normalized eigenvector of 3x3 symmetrical matrix M from its eigenvalues
+// using Cayley–Hamilton theorem according to:
+// https://en.wikipedia.org/wiki/Eigenvalue_algorithm 
+function _eigenvec_symm_3(M,evals,i=0) = 
+  let(  A  = (M - evals[(i+1)%3]*ident(3)) * (M - evals[(i+2)%3]*ident(3)) ,
+        k = max_index( [for(i=[0:2]) norm(A[i]) ]) 
+     )
+  norm(A[k])<EPSILON ? ident(3)[k] : A[k]/norm(A[k]);
+
+
+// eigenvalues of the covariance matrix of points
+function _covariance_evals(points) =
+    let(  pm   = sum(points)/len(points), // mean point
+          Y    = [ for(i=[0:len(points)-1]) points[i] - pm ],
+          M    = transpose(Y)*Y ,     // covariance matrix
+          evals = _eigenvals_symm_3(M) )
+    [pm, evals, M ];
+          
+
 // Function: plane_from_points()
 // Usage:
 //   plane_from_points(points, <fast>, <eps>);
@@ -900,7 +938,7 @@ function plane_from_normal(normal, pt=[0,0,0]) =
 //   Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane,
 //   that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1.
 //   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.
+//   If `fast` is true, the polygon coplanarity check is skipped and a best fitted plane 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
@@ -914,17 +952,18 @@ function plane_from_normal(normal, pt=[0,0,0]) =
 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 non-negative value." )
-    let(
-        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)
-    )
-    fast || points_on_plane(points,plane,eps=eps) ? plane : undef;
+    len(points) == 3 
+    ?   let( plane = plane3pt(points[0],points[1],points[2]) )
+        plane==[] ? undef : plane    
+    :   let(  
+            cov_evals = _covariance_evals(points),
+            pm        = cov_evals[0],
+            evals     = cov_evals[1],
+            M         = cov_evals[2],
+            evec      = _eigenvec_symm_3(M,evals,i=2) )
+// echo(error_points_plane= abs(max(points*evec)-pm*evec), limit=eps)
+        !fast && abs(max(points*evec)-pm*evec)>eps*evals[0] ? undef :
+        [ each evec, pm*evec] ;
 
 
 // Function: plane_from_polygon()
@@ -934,7 +973,8 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
 //   Given a 3D planar polygon, returns the normalized cartesian equation of its plane.
 //   Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1.
 //   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.
+//   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:
 //   poly = The planar 3D polygon to find the plane of.
 //   fast = If true, doesn't verify that all points in the polygon are coplanar.  Default: false
@@ -948,14 +988,13 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
 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 non-negative value." )
-    let(
-        poly = deduplicate(poly),
-        n = polygon_normal(poly),
-        plane = [n.x, n.y, n.z, n*poly[0]]
-    )
-    fast? plane: coplanar(poly,eps=eps)? plane: [];
-
-
+    len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) :
+    let( triple = sort(noncollinear_triple(poly,error=false)) )
+    triple==[] ? [] :
+    let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]]))
+    fast? plane: points_on_plane(poly, plane, eps=eps)? plane: [];
+ 
+ 
 // Function: plane_normal()
 // Usage:
 //   plane_normal(plane);
@@ -1252,9 +1291,11 @@ function coplanar(points, eps=EPSILON) =
     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);
+        let( 
+            plane  = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]),
+            normal = point3d(plane),
+            pt_nrm = points*normal )
+        abs(max(max(pt_nrm)-plane[3], -min(pt_nrm)+plane[3])) < eps;
 
 
 // Function: points_on_plane()
@@ -1665,11 +1706,11 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
             n = (pb-pa)/nrm,
             distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]
            )
-        max(distlist)<eps
+        max(distlist)<eps*nrm
         ?  assert(!error, "Cannot find three noncollinear points in pointlist.")
            []
         :  [0,b,max_index(distlist)];
-
+        
 
 // Function: pointlist_bounds()
 // Usage:
@@ -1746,9 +1787,9 @@ function polygon_area(poly, signed=false) =
                         v1 = poly[i] - poly[0],
                         v2 = poly[i+1] - poly[0]
                     )
-                    cross(v1,v2) * n
-                ])/2
-        )
+                    cross(v1,v2) 
+                ])* n/2
+        ) 
         signed ? total : abs(total);
 
 
@@ -2008,7 +2049,7 @@ function point_in_polygon(point, poly, nonzero=true, eps=EPSILON) =
 //   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;
+    polygon_area(poly, signed=true)<-EPSILON;
 
 
 // Function: clockwise_polygon()
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 6fe00c1..69d4705 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -197,8 +197,8 @@ module test_plane_from_polygon(){
     poly1 = [ rands(-1,1,3), rands(-1,1,3)+[2,0,0], rands(-1,1,3)+[0,2,2] ];
     poly2 = concat(poly1, [sum(poly1)/3] );
     info = info_str([["poly1 = ",poly1],["poly2 = ",poly2]]);
-    assert_std(plane_from_polygon(poly1),plane3pt(poly1[0],poly1[1],poly1[2]),info);
-    assert_std(plane_from_polygon(poly2),plane3pt(poly1[0],poly1[1],poly1[2]),info);
+    assert_approx(plane_from_polygon(poly1),plane3pt(poly1[0],poly1[1],poly1[2]),info);
+    assert_approx(plane_from_polygon(poly2),plane3pt(poly1[0],poly1[1],poly1[2]),info);
 }
 *test_plane_from_polygon();
 
@@ -208,8 +208,7 @@ module test_plane_from_normal(){
     displ = normal*point;
     info = info_str([["normal = ",normal],["point = ",point],["displ = ",displ]]);
     assert_approx(plane_from_normal(normal,point)*[each point,-1],0,info);
-    assert_std(plane_from_normal(normal,point),normalize_plane([each normal,displ]),info);
-    assert_std(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]);
+    assert_approx(plane_from_normal([1,1,1],[1,2,3]),[0.57735026919,0.57735026919,0.57735026919,3.46410161514]);
 }
 *test_plane_from_normal();
 
@@ -680,23 +679,23 @@ module test_triangle_area() {
 
 
 module test_plane3pt() {
-    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]);
+    assert_approx(plane3pt([0,0,20], [0,10,10], [0,0,0]), [1,0,0,0]);
+    assert_approx(plane3pt([2,0,20], [2,10,10], [2,0,0]), [1,0,0,2]);
+    assert_approx(plane3pt([0,0,0], [10,0,10], [0,0,20]), [0,1,0,0]);
+    assert_approx(plane3pt([0,2,0], [10,2,10], [0,2,20]), [0,1,0,2]);
+    assert_approx(plane3pt([0,0,0], [10,10,0], [20,0,0]), [0,0,1,0]);
+    assert_approx(plane3pt([0,0,2], [10,10,2], [20,0,2]), [0,0,1,2]);
 }
 *test_plane3pt();
 
 module test_plane3pt_indexed() {
     pts = [ [0,0,0], [10,0,0], [0,10,0], [0,0,10] ];
     s13 = sqrt(1/3);
-    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,3,2), [1,0,0,0]);
+    assert_approx(plane3pt_indexed(pts, 0,2,3), [-1,0,0,0]);
+    assert_approx(plane3pt_indexed(pts, 0,1,3), [0,1,0,0]);
+    assert_approx(plane3pt_indexed(pts, 0,3,1), [0,-1,0,0]);
+    assert_approx(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]);
@@ -709,18 +708,18 @@ module test_plane_from_points() {
     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]);
+    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();
 
 
 module test_plane_normal() {
-    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]);
+    assert_approx(plane_normal(plane3pt([0,0,20], [0,10,10], [0,0,0])), [1,0,0]);
+    assert_approx(plane_normal(plane3pt([2,0,20], [2,10,10], [2,0,0])), [1,0,0]);
+    assert_approx(plane_normal(plane3pt([0,0,0], [10,0,10], [0,0,20])), [0,1,0]);
+    assert_approx(plane_normal(plane3pt([0,2,0], [10,2,10], [0,2,20])), [0,1,0]);
+    assert_approx(plane_normal(plane3pt([0,0,0], [10,10,0], [20,0,0])), [0,0,1]);
+    assert_approx(plane_normal(plane3pt([0,0,2], [10,10,2], [20,0,2])), [0,0,1]);
 }
 *test_plane_normal();
 
@@ -780,7 +779,7 @@ module test_coplanar() {
     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();
 
 
@@ -836,7 +835,9 @@ module 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),signed=true), -PI*50*50, eps=0.1));
-    assert(approx(polygon_area(rot([13,27,75],p=path3d(circle(r=50,$fn=1000),fill=23)),signed=true), PI*50*50, eps=0.1));
+    assert(approx(polygon_area(rot([13,27,75],
+                               p=path3d(circle(r=50,$fn=1000),fill=23)),
+                               signed=true), -PI*50*50, eps=0.1));
 }
 *test_polygon_area();