normalized project_plane and lift_plane to match other transform

functions.
This commit is contained in:
Adrian Mariano 2021-04-21 22:49:06 -04:00
parent 94abf65857
commit 493ef62826
8 changed files with 158 additions and 140 deletions

View File

@ -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()

View File

@ -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];

View File

@ -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(

View File

@ -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

View File

@ -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) :

View File

@ -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

View File

@ -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();

View File

@ -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) =