Revert "Review of geometry.scad for speed"

This reverts commit 49a3a166eb.

arrays.scad

@@ -41,7 +41,6 @@ function is_homogeneous(l, depth=10) =
     [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) )  0 ];
 
 function is_homogenous(l, depth=10) = is_homogeneous(l, depth);
-                 
 
 function _same_type(a,b, depth) = 
     (depth==0) ||
@@ -51,7 +50,7 @@ function _same_type(a,b, depth) =
     (is_string(a) && is_string(b)) ||
     (is_list(a) && is_list(b) && len(a)==len(b) 
           && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); 
-  
+
 
 // Function: select()
 // Topics: List Handling
@@ -98,6 +97,7 @@ function select(list, start, end) =
 
 
 // Function: slice()
+// Topics: List Handling
 // Usage:
 //   list = slice(list,s,e);
 // Description:
@@ -476,7 +476,7 @@ function reverse(x) =
 //   l9 = list_rotate([1,2,3,4,5],6);  // Returns: [2,3,4,5,1]
 function list_rotate(list,n=1) =
     assert(is_list(list)||is_string(list), "Invalid list or string.")
-    assert(is_int(n), "The rotation number should be integer")
+    assert(is_finite(n), "Invalid number")
     let (
         ll = len(list),
         n = ((n % ll) + ll) % ll,
@@ -990,7 +990,7 @@ function _sort_vectors(arr, idxlist, _i=0) =
         _sort_vectors(equal,   idxlist, _i+1), 
         _sort_vectors(greater, idxlist, _i  ) );
         
- 
+
 // sorting using compare_vals(); returns indexed list when `indexed==true`
 function _sort_general(arr, idx=undef, indexed=false) =
     (len(arr)<=1) ? arr :
@@ -1332,8 +1332,6 @@ function permutations(l,n=2) =
 //   pairs = zip(a,b);
 //   triples = zip(a,b,c);
 //   quads = zip([LIST1,LIST2,LIST3,LIST4]);
-// Topics: List Handling, Iteration
-// See Also: zip_long()
 // Description:
 //   Zips together two or more lists into a single list.  For example, if you have two
 //   lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]].
@@ -1359,8 +1357,6 @@ function zip(a,b,c) =
 //   pairs = zip_long(a,b);
 //   triples = zip_long(a,b,c);
 //   quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
-// Topics: List Handling, Iteration
-// See Also: zip()
 // Description:
 //   Zips together two or more lists into a single list.  For example, if you have two
 //   lists [3,4,5], and [8,7,6], and zip them together, you get [[3,8],[4,7],[5,6]].
@@ -1530,6 +1526,7 @@ function subindex(M, idx) =
 //        [[4,2],   91, false],
 //        [6,    [3,4], undef]];
 //   submatrix(A,[0,2],[1,2]);   // Returns [[17, "test"], [[3, 4], undef]]
+
 function submatrix(M,idx1,idx2) =
     [for(i=idx1) [for(j=idx2) M[i][j] ] ];
 
@@ -1632,6 +1629,7 @@ function block_matrix(M) =
     assert(badrows==[], "Inconsistent or invalid input")
     bigM;
 
+
 // Function: diagonal_matrix()
 // Usage:
 //   mat = diagonal_matrix(diag, <offdiag>);
@@ -1857,7 +1855,7 @@ function transpose(arr, reverse=false) =
 //   A = matrix to test
 //   eps = epsilon for comparing equality.  Default: 1e-12
 function is_matrix_symmetric(A,eps=1e-12) =
-    approx(A,transpose(A), eps);
+    approx(A,transpose(A));
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

common.scad

@@ -205,8 +205,7 @@ function is_func(x) = version_num()>20210000 && is_function(x);
 // Description:
 //   Tests whether input is a list of entries which all have the same list structure
 //   and are filled with finite numerical data.  You can optionally specify a required 
-//   list structure with the pattern argument.  
+//   list structure with the pattern argument.  It returns `true` for the empty list.
-//   It returns `true` for the empty list regardless the value of the `pattern`.
 // Arguments:
 //   list = list to check
 //   pattern = optional pattern required to match
@@ -294,7 +293,7 @@ function default(v,dflt=undef) = is_undef(v)? dflt : v;
 //   v = The list whose items are being checked.
 //   recursive = If true, sublists are checked recursively for defined values.  The first sublist that has a defined item is returned.
 // Examples:
-//   val = first_defined([undef,7,undef,true]);  // Returns: 7
+//   val = first_defined([undef,7,undef,true]);  // Returns: 1
 function first_defined(v,recursive=false,_i=0) =
     _i<len(v) && (
         is_undef(v[_i]) || (
@@ -606,15 +605,15 @@ function segs(r) =
 
 
 // Module: no_children()
+// Topics: Error Checking
 // Usage:
 //   no_children($children);
-// Topics: Error Checking
-// See Also: no_function(), no_module()
 // Description:
 //   Assert that the calling module does not support children.  Prints an error message to this effect and fails if children are present,
 //   as indicated by its argument.
 // Arguments:
 //   $children = number of children the module has.  
+// See Also: no_function(), no_module()
 // Example:
 //   module foo() {
 //       no_children($children);
@@ -677,7 +676,7 @@ function _valstr(x) =
 //   expected = The value that was expected.
 //   info = Extra info to print out to make the error clearer.
 // Example:
-//   assert_approx(1/3, 0.333333333333333, str("number=",1,", demon=",3));
+//   assert_approx(1/3, 0.333333333333333, str("numer=",1,", demon=",3));
 module assert_approx(got, expected, info) {
     no_children($children);
     if (!approx(got, expected)) {

geometry.scad

@@ -888,58 +888,22 @@ 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),
+  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." )
+          "Inputs `normal` and `pt` should 3d vectors/points and `normal` cannot be zero." )
-    concat(normal, normal*pt) / norm(normal);
+  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>);
 // Description:
 //   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 the points in the list are collinear or not coplanar, then `undef` is returned.
+//   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
 //   eps = Tolerance in geometric comparisons.  Default: `EPSILON` (1e-9)
 // Example(3D):
 //   xyzpath = rot(45, v=[-0.3,1,0], p=path3d(star(n=6,id=70,d=100), 70));
@@ -947,21 +911,20 @@ function _covariance_evals(points) =
 //   #stroke(xyzpath,closed=true);
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
-function plane_from_points(points,fast=false, eps=EPSILON) =  
+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." )
-    len(points) == 3 
+    let(
-    ?   let( plane = plane3pt(points[0],points[1],points[2]) )
+        indices = noncollinear_triple(points,error=false)
-        plane==[] ? undef : plane    
+    )
-    :   let(  
+    indices==[] ? undef :
-            cov_evals = _covariance_evals(points),
+    let(
-            pm        = cov_evals[0],
+        p1 = points[indices[0]],
-            evals     = cov_evals[1],
+        p2 = points[indices[1]],
-            M         = cov_evals[2],
+        p3 = points[indices[2]],
-            evec      = _eigenvec_symm_3(M,evals,i=2) )
+        plane = plane3pt(p1,p2,p3)
-// echo(error_points_plane= abs(max(points*evec)-pm*evec), limit=eps)
+    )
-        !fast && abs(max(points*evec)-pm*evec)>eps*evals[0] ? undef :
+    fast || points_on_plane(points,plane,eps=eps) ? plane : undef;
-        [ each evec, pm*evec] ;
 
 
 // Function: plane_from_polygon()
@@ -982,16 +945,17 @@ function plane_from_points(points,fast=false, eps=EPSILON) =
 //   #stroke(xyzpath,closed=true);
 //   cp = centroid(xyzpath);
 //   move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow();
-function plane_from_polygon(poly, 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." )
-		len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) :
+    let(
-    let( triple = sort(noncollinear_triple(poly,error=false)) )
+        poly = deduplicate(poly),
-    triple==[] ? [] :
+        n = polygon_normal(poly),
-    let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]]))
+        plane = [n.x, n.y, n.z, n*poly[0]]
-    fast? plane: points_on_plane(poly, plane, eps=eps)? plane: [];
+    )
- 
+    fast? plane: coplanar(poly,eps=eps)? plane: [];
- 
+
+
 // Function: plane_normal()
 // Usage:
 //   plane_normal(plane);
@@ -1288,11 +1252,9 @@ function coplanar(points, eps=EPSILON) =
     len(points)<=2 ? false
     :   let( ip = noncollinear_triple(points,error=false,eps=eps) )
         ip == [] ? false :
-        let( 
+        let( plane  = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]),
-						plane  = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]),
+             normal = point3d(plane) )
-            normal = point3d(plane),
+        max( points*normal ) - plane[3]< eps*norm(normal);
-						pt_nrm = points*normal )
-				abs(max(max(pt_nrm)-plane[3], -min(pt_nrm)+plane[3])) < eps;
 
 
 // Function: points_on_plane()
@@ -1703,11 +1665,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*nrm
+        max(distlist)<eps
         ?  assert(!error, "Cannot find three noncollinear points in pointlist.")
            []
         :  [0,b,max_index(distlist)];
-        
+
 
 // Function: pointlist_bounds()
 // Usage:
@@ -1784,9 +1746,9 @@ function polygon_area(poly, signed=false) =
                         v1 = poly[i] - poly[0],
                         v2 = poly[i+1] - poly[0]
                     )
-                    cross(v1,v2) 
+                    cross(v1,v2) * n
-                ])* n/2
+                ])/2
-        ) 
+        )
         signed ? total : abs(total);
 
 
@@ -1958,7 +1920,7 @@ function centroid(poly, eps=EPSILON) =
     let(
         n = len(poly[0])==2 ? 1 :
             let( 
-                plane = plane_from_points(poly) )
+                plane = plane_from_points(poly, fast=true) )
             assert( !is_undef(plane), "The polygon must be planar." )
             plane_normal(plane),
         v0 = poly[0] ,
@@ -2046,7 +2008,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)<-EPSILON;
+    polygon_area(poly, signed=true)<0;
 
 
 // Function: clockwise_polygon()

tests/test_geometry.scad

@@ -705,11 +705,11 @@ module test_plane3pt_indexed() {
 
 module test_plane_from_points() {
     assert_std(plane_from_points([[0,0,20], [0,10,10], [0,0,0], [0,5,3]]), [1,0,0,0]);
-    assert_approx(plane_from_points([[2,0,20], [2,10,10], [2,0,0], [2,3,4]]), [1,0,0,2]);
+    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]);  
+    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();
 
@@ -836,9 +836,7 @@ 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],
+    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));
-                               p=path3d(circle(r=50,$fn=1000),fill=23)),
-                               signed=true), -PI*50*50, eps=0.1));
 }
 *test_polygon_area();