From 493ef62826e0b7d926f8de2c72fdf12cc3764f07 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Wed, 21 Apr 2021 22:49:06 -0400
Subject: [PATCH 1/2] normalized project_plane and lift_plane to match other
 transform functions.

---
 coords.scad            | 207 ++++++++++++++++++++++-------------------
 geometry.scad          |  35 +------
 hull.scad              |   2 +-
 paths.scad             |   4 +-
 rounding.scad          |   2 +-
 shapes2d.scad          |   8 +-
 tests/test_coords.scad |  33 ++++++-
 vnf.scad               |   7 +-
 8 files changed, 158 insertions(+), 140 deletions(-)

diff --git a/coords.scad b/coords.scad
index 1843c7a..b32ad73 100644
--- a/coords.scad
+++ b/coords.scad
@@ -179,112 +179,131 @@ function xy_to_polar(x,y=undef) = let(
 
 
 // Function: project_plane()
-// Usage: With the plane defined by 3 Points
-//   pt = project_plane(point, a, b, c);
-// Usage: With the plane defined by Pointlist
-//   pt = project_plane(point, POINTLIST);
-// Usage: With the plane defined by Plane Definition [A,B,C,D] Where Ax+By+Cz=D
-//   pt = project_plane(point, PLANE);
-// Topics: Coordinates, Points, Paths
-// See Also: lift_plane()
+// Usage: 
+//   xy = project_plane(plane, p);
+// Usage: To get a transform matrix
+//   M = project_plane(plane)
 // Description:
-//   Converts the given 3D points from global coordinates to the 2D planar coordinates of the closest
-//   points on the plane.  This coordinate system can be useful in taking a set of nearly coplanar
-//   points, and converting them to a pure XY set of coordinates for manipulation, before converting
-//   them back to the original 3D plane. The parameter `point` may be a single point or a list of points
-//   The plane may be given in one of three ways:
-//   - by three points, `a`, `b`, and `c`, the planar coordinate system will have `[0,0]` at point `a`, and the Y+ axis will be towards point `b`.
-//   - by a list of points passed by `a`, finds three reasonably spaced non-collinear points in the list and uses them as points `a`, `b`, and `c` as above.
-//   - by a plane definition `[A,B,C,D]` passed by `a` where `Ax+By+Cz=D`, the closest point on that plane to the global origin at `[0,0,0]` will be the planar coordinate origin `[0,0]`.
+//   Maps the provided 3d point(s) from 3D coordinates to a 2d coordinate system defined by `plane`.  Points that are not
+//   on the specified plane will be projected orthogonally onto the plane.  This coordinate system is useful if you need
+//   to perform 2d operations on a coplanar set of data.  After those operations are done you can return the data
+//   to 3d with `lift_plane()`.  You could also use this to force approximately coplanar data to be exactly coplanar.
+//   The parameter p can be a point, path, region, bezier patch or VNF.
+//   The plane can be specified as
+//   - A list of three points.  The planar coordinate system will have [0,0] at plane[0], and plane[1] will lie on the Y+ axis.
+//   - A list of coplanar points that define a plane (not-collinear)
+//   - A plane definition `[A,B,C,D]` where `Ax+By+CZ=D`.  The closest point on that plane to the origin will map to the origin in the new coordinate system.
+//   .
+//   If you omit the point specification then `project_plane()` returns a rotation matrix that maps the specified plane to the XY plane.
+//   Note that if you apply this transformation to data lying on the plane it will produce 3D points with the Z coordinate of zero.
+// Topics: Coordinates, Points, Paths
+// See Also: project_plane(), projection_on_plane()
 // Arguments:
-//   point = The 3D point, or list of 3D points to project into the plane's 2D coordinate system.
-//   a = A 3D point that the plane passes through or a list of points or a plane definition vector.
-//   b = A 3D point that the plane passes through.  Used to define the plane.
-//   c = A 3D point that the plane passes through.  Used to define the plane.
+//   plane = plane specification or point list defining the plane
+//   p = 3D point, path, region, VNF or bezier patch to project
 // Example:
 //   pt = [5,-5,5];
 //   a=[0,0,0];  b=[10,-10,0];  c=[10,0,10];
-//   xy = project_plane(pt, a, b, c);
-//   xy2 = project_plane(pt, [a,b,c]);
-// Example(3D):
-//   points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100, $fn=36))));
-//   plane = plane_from_normal([1,0,1]);
-//   proj = project_plane(points,plane);
-//   n = plane_normal(plane);
-//   cp = centroid(proj);
-//   color("red") move_copies(points) sphere(d=2,$fn=12);
-//   color("blue") rot(from=UP,to=n,cp=cp) move_copies(proj) sphere(d=2,$fn=12);
-//   move(cp) {
-//       rot(from=UP,to=n) {
-//           anchor_arrow(30);
-//           %cube([120,150,0.1],center=true);
-//       }
-//   }
-function project_plane(point, a, b, c) =
-    is_undef(b) && is_undef(c) && is_list(a)? let(
-        mat = is_vector(a,4)? plane_transform(a) :
-            assert(is_path(a) && len(a)>=3)
-            plane_transform(plane_from_points(a)),
-        pts = is_vector(point)? point2d(apply(mat,point)) :
-            is_path(point)? path2d(apply(mat,point)) :
-            is_region(point)? [for (x=point) path2d(apply(mat,x))] :
-            assert(false, "point must be a 3D point, path, or region.")
-    ) pts :
-    assert(is_vector(a))
-    assert(is_vector(b))
-    assert(is_vector(c))
-    assert(is_vector(point)||is_path(point))
-    let(
-        u = unit(b-a),
-        v = unit(c-a),
-        n = unit(cross(u,v)),
-        w = unit(cross(n,u)),
-        relpoint = apply(move(-a),point)
-    ) relpoint * transpose([w,u]);
+//   xy = project_plane([a,b,c],pt);
+// Example(3D): The yellow points in 3D project onto the red points in 2D
+//   M = [[-1, 2, -1, -2], [-1, -3, 2, -1], [2, 3, 4, 53], [0, 0, 0, 1]];
+//   data = apply(M,path3d(circle(r=10, $fn=20)));
+//   move_copies(data) sphere(r=1);
+//   color("red") move_copies(project_plane(data, data)) sphere(r=1);
+// Example:
+//   xyzpath = move([10,20,30], p=yrot(25, p=path3d(circle(d=100))));
+//   mat = project_plane(xyzpath);
+//   xypath = path2d(apply(mat, xyzpath));
+//   #stroke(xyzpath,closed=true);
+//   stroke(xypath,closed=true);
+function project_plane(plane,p) =
+      is_matrix(plane,3,3) && is_undef(p) ? // no data, 3 points given
+          assert(!collinear(plane),"Points defining the plane must not be collinear")
+          let(
+              v = plane[2]-plane[0],
+              y = unit(plane[1]-plane[0]),        // y axis goes to point b
+              x = unit(v-(v*y)*y)   // x axis 
+          )            
+          affine3d_frame_map(x,y) * move(-plane[0])
+    : is_vector(plane,4) && is_undef(p) ?            // no data, plane given in "plane"
+          assert(_valid_plane(plane), "Plane is not valid")
+          let(
+               n = point3d(plane),
+               cp = n * plane[3] / (n*n)
+          )
+          rot(from=n, to=UP) * move(-cp)
+    : is_path(plane,3) && is_undef(p) ?               // no data, generic point list plane
+          assert(len(plane)>=3, "Need three points to define a plane")
+          let(plane = plane_from_points(plane))
+          assert(is_def(plane), "Point list is not coplanar")
+          project_plane(plane)
+    : assert(is_def(p), str("Invalid plane specification",plane))
+      is_vnf(p) ? [project_plane(plane,p[0]), p[1]] 
+    : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ?  // bezier patch or region
+           [for(plist=p) project_plane(plane,plist)]
+    : assert(is_vector(p,3) || is_path(p,3),str("Data must be a 3d point, path, region, vnf or bezier patch",p))
+      is_matrix(plane,3,3) ?
+          assert(!collinear(plane),"Points defining the plane must not be collinear")
+          let(
+              v = plane[2]-plane[0],
+              y = unit(plane[1]-plane[0]),        // y axis goes to point b
+              x = unit(v-(v*y)*y)  // x axis 
+          ) move(-plane[0],p) * transpose([x,y])
+    : is_vector(p) ? point2d(apply(project_plane(plane),p))
+    : path2d(apply(project_plane(plane),p));
+
 
 
 // Function: lift_plane()
-// Usage: With 3 Points
-//   xyz = lift_plane(point, a, b, c);
-// Usage: With Pointlist
-//   xyz = lift_plane(point, POINTLIST);
-// Usage: With Plane Definition [A,B,C,D] Where Ax+By+Cz=D
-//   xyz = lift_plane(point, PLANE);
+// Usage: 
+//   xyz = lift_plane(plane, p);
+// Usage: to get transform matrix
+//   M =  lift_plane(plane);
 // Topics: Coordinates, Points, Paths
 // See Also: project_plane()
 // Description:
-//   Converts the given 2D point from planar coordinates to the global 3D coordinates of the point on the plane.
-//   Can be called one of three ways:
-//   - Given three points, `a`, `b`, and `c`, the planar coordinate system will have `[0,0]` at point `a`, and the Y+ axis will be towards point `b`.
-//   - Given a list of points, finds three non-collinear points in the list and uses them as points `a`, `b`, and `c` as above.
-//   - Given a plane definition `[A,B,C,D]` where `Ax+By+Cz=D`, the closest point on that plane to the global origin at `[0,0,0]` will be the planar coordinate origin `[0,0]`.
+//   Converts the given 2D point on the plane to 3D coordinates of the specified plane.
+//   The parameter p can be a point, path, region, bezier patch or VNF.
+//   The plane can be specified as
+//   - A list of three points.  The planar coordinate system will have [0,0] at plane[0], and plane[1] will lie on the Y+ axis.
+//   - A list of coplanar points that define a plane (not-collinear)
+//   - A plane definition `[A,B,C,D]` where `Ax+By+CZ=D`.  The closest point on that plane to the origin will map to the origin in the new coordinate system.
+// If you do not supply `p` then you get a transformation matrix which operates in 3D, assuming that the Z coordinate of the points is zero.
+// This matrix is a rotation, the inverse of the one produced by project_plane.
 // Arguments:
-//   point = The 2D point, or list of 2D points in the plane's coordinate system to get the 3D position of.
-//   a = A 3D point that the plane passes through.  Used to define the plane.
-//   b = A 3D point that the plane passes through.  Used to define the plane.
-//   c = A 3D point that the plane passes through.  Used to define the plane.
-function lift_plane(point, a, b, c) =
-    is_undef(b) && is_undef(c) && is_list(a)? let(
-        mat = is_vector(a,4)? plane_transform(a) :
-            assert(is_path(a) && len(a)>=3)
-            plane_transform(plane_from_points(a)),
-        imat = matrix_inverse(mat),
-        pts = is_vector(point)? apply(imat,point3d(point)) :
-            is_path(point)? apply(imat,path3d(point)) :
-            is_region(point)? [for (x=point) apply(imat,path3d(x))] :
-            assert(false, "point must be a 2D point, path, or region.")
-    ) pts :
-    assert(is_vector(a))
-    assert(is_vector(b))
-    assert(is_vector(c))
-    assert(is_vector(point)||is_path(point))
-    let(
-        u = unit(b-a),
-        v = unit(c-a),
-        n = unit(cross(u,v)),
-        w = unit(cross(n,u)),
-        remapped = point*[w,u]
-    ) apply(move(a),remapped);
+//   plane = Plane specification or list of points to define a plane
+//   p = points, path, region, VNF, or bezier patch to transform. 
+function lift_plane(plane, p) =
+      is_matrix(plane,3,3) && is_undef(p) ? // no data, 3 p given
+          let(
+              v = plane[2]-plane[0],
+              y = unit(plane[1]-plane[0]),        // y axis goes to point b
+              x = unit(v-(v*y)*y)   // x axis 
+          )            
+          move(plane[0]) * affine3d_frame_map(x,y,reverse=true)
+    : is_vector(plane,4) && is_undef(p) ?            // no data, plane given in "plane"
+          assert(_valid_plane(plane), "Plane is not valid")
+          let(
+               n = point3d(plane),
+               cp = n * plane[3] / (n*n)
+          )
+          move(cp) * rot(from=UP, to=n)
+    : is_path(plane,3) && is_undef(p) ?               // no data, generic point list plane
+          assert(len(plane)>=3, "Need three p to define a plane")
+          let(plane = plane_from_points(plane))
+          assert(is_def(plane), "Point list is not coplanar")
+          lift_plane(plane)
+    : is_vnf(p) ? [lift_plane(plane,p[0]), p[1]] 
+    : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ?  // bezier patch or region
+           [for(plist=p) lift_plane(plane,plist)]
+    : assert(is_vector(p,2) || is_path(p,2),"Data must be a 2d point, path, region, vnf or bezier patch")
+      is_matrix(plane,3,3) ?
+          let(
+              v = plane[2]-plane[0],
+              y = unit(plane[1]-plane[0]),        // y axis goes to point b
+              x = unit(v-(v*y)*y)  // x axis 
+          ) move(plane[0],p * [x,y])
+    : apply(lift_plane(plane),is_vector(p) ? point3d(p) : path3d(p));
 
 
 // Function: cylindrical_to_xyz()
diff --git a/geometry.scad b/geometry.scad
index 408393d..24fbac6 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -1024,31 +1024,6 @@ function plane_offset(plane) =
     plane[3]/norm([plane.x, plane.y, plane.z]);
 
 
-// Function: plane_transform()
-// Usage:
-//   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 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.
-// Example(3D):
-//   xyzpath = move([10,20,30], p=yrot(25, p=path3d(circle(d=100))));
-//   plane = plane_from_points(xyzpath);
-//   mat = plane_transform(plane);
-//   xypath = path2d(apply(mat, xyzpath));
-//   #stroke(xyzpath,closed=true);
-//   stroke(xypath,closed=true);
-function plane_transform(plane) =
-    let(
-        plane = normalize_plane(plane),
-        n = point3d(plane),
-        cp = n * plane[3]
-        )
-    rot(from=n, to=UP) * move(-cp);
-
 
 // Function: projection_on_plane()
 // Usage:
@@ -1227,8 +1202,8 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
             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, plane)),
-            line2d = project_plane([lp1,lp2], plane),
+            poly2d = clockwise_polygon(project_plane(plane, poly)),
+            line2d = project_plane(plane, [lp1,lp2]),
             parts = split_path_at_region_crossings(line2d, [poly2d], closed=false),
             inside = [for (part = parts)
                           if (point_in_polygon(mean(part), poly2d)>0) part
@@ -1236,15 +1211,15 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
         )
         !inside? undef :
         let(
-            isegs = [for (seg = inside) lift_plane(seg, plane) ]
+            isegs = [for (seg = inside) lift_plane(plane, seg) ]
         )
         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)
+            proj = clockwise_polygon(project_plane([p1, p2, p3], poly)),
+            pt = project_plane([p1, p2, p3], res[0])
         )
         point_in_polygon(pt, proj) < 0 ? undef : res[0];
 
diff --git a/hull.scad b/hull.scad
index 6961e1a..f18e4ab 100644
--- a/hull.scad
+++ b/hull.scad
@@ -173,7 +173,7 @@ function hull3d_faces(points) =
     d == len(points)
   ? /* all coplanar*/
     let (
-        pts2d = [ for (p = points) project_plane(p, points[a], points[b], points[c]) ],
+        pts2d =  project_plane([points[a], points[b], points[c]],points),
         hull2d = hull2d_path(pts2d)
     ) hull2d
   : let(
diff --git a/paths.scad b/paths.scad
index 9fa8d73..645824e 100644
--- a/paths.scad
+++ b/paths.scad
@@ -1222,14 +1222,14 @@ module path_extrude(path, convexity=10, clipsize=100) {
 //     path_spread(wedge,n=5,spacing=3) fwd(.1) rect([1,4],anchor=FRONT);
 //   }
 // Example(Spin,VPD=115): 3d example, with children rotated into the plane of the path
-//   tilted_circle = lift_plane(regular_ngon(n=64, or=12), [0,0,0], [5,0,5], [0,2,3]);
+//   tilted_circle = lift_plane([[0,0,0], [5,0,5], [0,2,3]],regular_ngon(n=64, or=12));
 //   path_sweep(regular_ngon(n=16,or=.1),tilted_circle);
 //   path_spread(tilted_circle, n=15,closed=true) {
 //      color("blue") cyl(h=3,r=.2, anchor=BOTTOM);      // z-aligned cylinder
 //      color("red") xcyl(h=10,r=.2, anchor=FRONT+LEFT); // x-aligned cylinder
 //   }
 // Example(Spin,VPD=115): 3d example, with rotate_children set to false
-//   tilted_circle = lift_plane(regular_ngon(n=64, or=12), [0,0,0], [5,0,5], [0,2,3]);
+//   tilted_circle = lift_plane([[0,0,0], [5,0,5], [0,2,3]], regular_ngon(n=64, or=12));
 //   path_sweep(regular_ngon(n=16,or=.1),tilted_circle);
 //   path_spread(tilted_circle, n=25,rotate_children=false,closed=true) {
 //      color("blue") cyl(h=3,r=.2, anchor=BOTTOM);       // z-aligned cylinder
diff --git a/rounding.scad b/rounding.scad
index 0b980b3..eaf8e66 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -1827,7 +1827,7 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
    assert(len(bottom[0])==3 || is_num(height),"Must give height/length with 2d polygon input")
    let(
      // 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)),
+     bot_proj = len(bottom[0])==2 ? bottom :  project_plane(select(bottom,0,2),bottom),
      bottom_sign = polygon_is_clockwise(bot_proj) ? 1 : -1,
      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) :
diff --git a/shapes2d.scad b/shapes2d.scad
index e14467e..2fa1ea1 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -569,11 +569,11 @@ function arc(N, r, angle, d, cp, points, width, thickness, start, wedge=false, l
                 assert(!(cw || ccw), "(Counter)clockwise isn't meaningful in 3d, so `cw` and `ccw` must be false")
                 assert(is_undef(cp) || is_vector(cp,3),"points are 3d so cp must be 3d")
         let(
-            thirdpoint = is_def(cp) ? cp : points[2],
-            center2d = is_def(cp) ? project_plane(cp,thirdpoint,points[0],points[1]) : undef,
-            points2d = project_plane(points,thirdpoint,points[0],points[1])
+            plane = [is_def(cp) ? cp : points[2], points[0], points[1]],
+            center2d = is_def(cp) ? project_plane(plane,cp) : undef,
+            points2d = project_plane(plane, points)
         )
-        lift_plane(arc(N,cp=center2d,points=points2d,wedge=wedge,long=long),thirdpoint,points[0],points[1])
+        lift_plane(plane,arc(N,cp=center2d,points=points2d,wedge=wedge,long=long))
     ) : is_def(cp)? (
         // Arc defined by center plus two points, will have radius defined by center and points[0]
         // and extent defined by direction of point[1] from the center
diff --git a/tests/test_coords.scad b/tests/test_coords.scad
index 9fccd65..64343f8 100644
--- a/tests/test_coords.scad
+++ b/tests/test_coords.scad
@@ -83,15 +83,40 @@ test_xy_to_polar();
 
 
 module test_project_plane() {
-    assert(approx(project_plane([-5,0,-5], [-10,0,-10], [0,0,0], [0,-10,-10]),[0,10*sqrt(2)/2]));
-    assert(approx(project_plane([0,-5,-5], [-10,0,-10], [0,0,0], [0,-10,-10]),[6.12372, 10.6066],eps=1e-5));
+    assert(approx(project_plane([[-10,0,-10], [0,0,0], [0,-10,-10]],[-5,0,-5]),[0,10*sqrt(2)/2]));
+    assert(approx(project_plane([[-10,0,-10], [0,0,0], [0,-10,-10]],[0,-5,-5]),[6.12372, 10.6066],eps=1e-5));
+    assert_approx(project_plane([[3,4,5],[1,3,9],[4,7,13]], [[3,4,5],[1,3,9],[5,3,2]]),[[0,0],[0,4.58257569496],[-0.911684611677,-3.27326835354]]);
+    assert_approx(project_plane([[3,4,5],[1,3,9],[4,7,13]], [[3,4,5],[1,3,9],[4,7,13]]),[[0,0],[0,4.58257569496],[6.26783170528,5.89188303637]]);
+
+    assert_approx(project_plane([2,3,4,2], [4,2,3]),[2.33181857677,-0.502272134844]);
+    assert_approx(project_plane([2,3,4,2], [[1,1,1],[0,0,0]]),[[0.430748825729,0.146123238594],[0,0]]);
+    assert_approx(project_plane([2,3,4,2]),[[0.920855800833,-0.11871629875,-0.371390676354,0],[-0.11871629875,0.821925551875,-0.557086014531,-2.77555756156e-17],[0.371390676354,0.557086014531,0.742781352708,-0.371390676354],[0,0,0,1]]);
+    assert_approx(project_plane([[1,1,1],[3,1,3],[1,1,4]]),[[-1/sqrt(2),1/sqrt(2),0,0],[0,0,1,-1],[1/sqrt(2),1/sqrt(2),0,-sqrt(2)],[0,0,0,1]]);
 }
 test_project_plane();
 
 
 module test_lift_plane() {
-    assert(approx(lift_plane([0,10*sqrt(2)/2], [-10,0,-10], [0,0,0], [0,-10,-10]),[-5,0,-5]));
-    assert(approx(lift_plane([6.12372, 10.6066], [-10,0,-10], [0,0,0], [0,-10,-10]),[0,-5,-5],eps=1e-5));
+    assert(approx(lift_plane([[-10,0,-10], [0,0,0], [0,-10,-10]],[0,10*sqrt(2)/2]),[-5,0,-5]));
+    assert(approx(lift_plane([[-10,0,-10], [0,0,0], [0,-10,-10]],[6.12372, 10.6066]),[0,-5,-5],eps=1e-5));
+
+    assert_approx(lift_plane([[3,4,5],[1,3,9],[4,7,13]], [[0,0],[0,4.58257569496],[6.26783170528,5.89188303637]]),[[3,4,5],[1,3,9],[4,7,13]]);
+
+    assert_approx(project_plane([2,3,4,2]),[[0.920855800833,-0.11871629875,-0.371390676354,0],[-0.11871629875,0.821925551875,-0.557086014531,-2.77555756156e-17],[0.371390676354,0.557086014531,0.742781352708,-0.371390676354],[0,0,0,1]]);
+    assert_approx(project_plane([[1,1,1],[3,1,3],[1,1,4]]),[[-1/sqrt(2),1/sqrt(2),0,0],[0,0,1,-1],[1/sqrt(2),1/sqrt(2),0,-sqrt(2)],[0,0,0,1]]);
+
+    N=30;
+    data2 = array_group(rands(0,10,3*N,seed=77),3);
+    data3 = [for (d=data2) [d.x,d.y,d.x*3+d.y*5+2]];
+    planept = select(data3,0,N-4);
+    testpt = select(data3, N-3,-1);
+    newdata = project_plane(planept,testpt);
+    assert_approx( lift_plane(planept, newdata), testpt);
+    assert_approx( lift_plane(planept, project_plane(planept, last(testpt))), last(testpt));
+    assert_approx( lift_plane(planept) * project_plane(planept) , ident(4));
+    assert_approx( lift_plane([1,2,3,4]) * project_plane([1,2,3,4]) , ident(4));
+    assert_approx( lift_plane([[1,1,1],[3,1,3],[1,1,4]]) * project_plane([[1,1,1],[3,1,3],[1,1,4]]) , ident(4));        
+    
 }
 test_lift_plane();
 
diff --git a/vnf.scad b/vnf.scad
index c9ee06c..2613126 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -1084,11 +1084,10 @@ function vnf_halfspace(plane, vnf, closed=true) =
     len(newpaths)<=1 ? [newvert, concat(faces_edges_vertices[0], newpaths)] 
     :
       let(
-           faceregion = [for(p=newpaths) project_plane(select(newvert,p), plane)],
-           facevnf = region_faces(faceregion,reverse=true),
-           faceverts = lift_plane(facevnf[0], plane)
+           faceregion = project_plane(plane, newpaths),
+           facevnf = region_faces(faceregion,reverse=true)
       )
-      vnf_merge([[newvert, faces_edges_vertices[0]], [faceverts, facevnf[1]]]);
+      vnf_merge([[newvert, faces_edges_vertices[0]], lift_plane(plane, facevnf)]);
 
 
 function _assemble_paths(vertices, edges, paths=[],i=0) =

From d91a48ef2871d8a6b88c9032aacca5e0098a1d4f Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Thu, 22 Apr 2021 06:28:23 -0400
Subject: [PATCH 2/2] fix tests

---
 tests/test_coords.scad   |  5 +++++
 tests/test_geometry.scad | 11 -----------
 2 files changed, 5 insertions(+), 11 deletions(-)

diff --git a/tests/test_coords.scad b/tests/test_coords.scad
index 64343f8..0022efa 100644
--- a/tests/test_coords.scad
+++ b/tests/test_coords.scad
@@ -92,6 +92,11 @@ module test_project_plane() {
     assert_approx(project_plane([2,3,4,2], [[1,1,1],[0,0,0]]),[[0.430748825729,0.146123238594],[0,0]]);
     assert_approx(project_plane([2,3,4,2]),[[0.920855800833,-0.11871629875,-0.371390676354,0],[-0.11871629875,0.821925551875,-0.557086014531,-2.77555756156e-17],[0.371390676354,0.557086014531,0.742781352708,-0.371390676354],[0,0,0,1]]);
     assert_approx(project_plane([[1,1,1],[3,1,3],[1,1,4]]),[[-1/sqrt(2),1/sqrt(2),0,0],[0,0,1,-1],[1/sqrt(2),1/sqrt(2),0,-sqrt(2)],[0,0,0,1]]);
+
+    normal = rands(-1,1,3,seed=3)+[2,0,0];
+    offset = rands(-1,1,1,seed=4)[0];
+    assert_approx(project_plane([0,0,1,offset]),move([0,0,-offset]) );
+    assert_approx(project_plane([0,1,0,offset]),xrot(90)*move([0,-offset,0]) );
 }
 test_project_plane();
 
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index ccb4ac7..04f5126 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -42,7 +42,6 @@ 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();
@@ -175,16 +174,6 @@ module test_plane_point_nearest_origin(){
 test_plane_point_nearest_origin();
 
 
-module test_plane_transform(){
-    normal = rands(-1,1,3)+[2,0,0];
-    offset = rands(-1,1,1)[0];
-    info = info_str([["normal = ",normal],["offset = ",offset]]);
-    assert_approx(plane_transform([0,0,1,offset]),move([0,0,-offset]),info );
-    assert_approx(plane_transform([0,1,0,offset]),xrot(90)*move([0,-offset,0]),info );
-}
-*test_plane_transform();
-
-
 module test_plane_offset(){
     plane = rands(-1,1,4)+[2,0,0,0]; // a valid plane
     info = info_str([["plane = ",plane]]);