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