1//////////////////////////////////////////////////////////////////////
2// LibFile: coords.scad
3// Coordinate transformations and coordinate system conversions.
4// Includes:
5// include <BOSL2/std.scad>
6// FileGroup: Math
7// FileSummary: Conversions between coordinate systems.
8// FileFootnotes: STD=Included in std.scad
9//////////////////////////////////////////////////////////////////////
10
11
12// Section: Coordinate Manipulation
13
14// Function: point2d()
15// Usage:
16// pt = point2d(p, [fill]);
17// Topics: Coordinates, Points
18// See Also: path2d(), point3d(), path3d()
19// Description:
20// Returns a 2D vector/point from a 2D or 3D vector. If given a 3D point, removes the Z coordinate.
21// Arguments:
22// p = The coordinates to force into a 2D vector/point.
23// fill = Value to fill missing values in vector with. Default: 0
24function point2d(p, fill=0) = assert(is_list(p)) [for (i=[0:1]) (p[i]==undef)? fill : p[i]];
25
26
27// Function: path2d()
28// Usage:
29// pts = path2d(points);
30// Topics: Coordinates, Points, Paths
31// See Also: point2d(), point3d(), path3d()
32// Description:
33// Returns a list of 2D vectors/points from a list of 2D, 3D or higher dimensional vectors/points.
34// Removes the extra coordinates from higher dimensional points. The input must be a path, where
35// every vector has the same length.
36// Arguments:
37// points = A list of 2D or 3D points/vectors.
38function path2d(points) =
39 assert(is_path(points,dim=undef,fast=true),"Input to path2d is not a path")
40 let (result = points * concat(ident(2), repeat([0,0], len(points[0])-2)))
41 assert(is_def(result), "Invalid input to path2d")
42 result;
43
44
45// Function: point3d()
46// Usage:
47// pt = point3d(p, [fill]);
48// Topics: Coordinates, Points
49// See Also: path2d(), point2d(), path3d()
50// Description:
51// Returns a 3D vector/point from a 2D or 3D vector.
52// Arguments:
53// p = The coordinates to force into a 3D vector/point.
54// fill = Value to fill missing values in vector with. Default: 0
55function point3d(p, fill=0) =
56 assert(is_list(p))
57 [for (i=[0:2]) (p[i]==undef)? fill : p[i]];
58
59
60// Function: path3d()
61// Usage:
62// pts = path3d(points, [fill]);
63// Topics: Coordinates, Points, Paths
64// See Also: point2d(), path2d(), point3d()
65// Description:
66// Returns a list of 3D vectors/points from a list of 2D or higher dimensional vectors/points
67// by removing extra coordinates or adding the z coordinate.
68// Arguments:
69// points = A list of 2D, 3D or higher dimensional points/vectors.
70// fill = Value to fill missing values in vectors with (in the 2D case). Default: 0
71function path3d(points, fill=0) =
72 assert(is_num(fill))
73 assert(is_path(points, dim=undef, fast=true), "Input to path3d is not a path")
74 let (
75 change = len(points[0])-3,
76 M = change < 0? [[1,0,0],[0,1,0]] :
77 concat(ident(3), repeat([0,0,0],change)),
78 result = points*M
79 )
80 assert(is_def(result), "Input to path3d is invalid")
81 fill == 0 || change>=0 ? result : result + repeat([0,0,fill], len(result));
82
83
84// Function: point4d()
85// Usage:
86// pt = point4d(p, [fill]);
87// Topics: Coordinates, Points
88// See Also: point2d(), path2d(), point3d(), path3d(), path4d()
89// Description:
90// Returns a 4D vector/point from a 2D or 3D vector.
91// Arguments:
92// p = The coordinates to force into a 4D vector/point.
93// fill = Value to fill missing values in vector with. Default: 0
94function point4d(p, fill=0) = assert(is_list(p))
95 [for (i=[0:3]) (p[i]==undef)? fill : p[i]];
96
97
98// Function: path4d()
99// Usage:
100// pt = path4d(points, [fill]);
101// Topics: Coordinates, Points, Paths
102// See Also: point2d(), path2d(), point3d(), path3d(), point4d()
103// Description:
104// Returns a list of 4D vectors/points from a list of 2D or 3D vectors/points.
105// Arguments:
106// points = A list of 2D or 3D points/vectors.
107// fill = Value to fill missing values in vectors with. Default: 0
108function path4d(points, fill=0) =
109 assert(is_num(fill) || is_vector(fill))
110 assert(is_path(points, dim=undef, fast=true), "Input to path4d is not a path")
111 let (
112 change = len(points[0])-4,
113 M = change < 0 ? select(ident(4), 0, len(points[0])-1) :
114 concat(ident(4), repeat([0,0,0,0],change)),
115 result = points*M
116 )
117 assert(is_def(result), "Input to path4d is invalid")
118 fill == 0 || change >= 0 ? result :
119 let(
120 addition = is_list(fill) ? concat(0*points[0],fill) :
121 concat(0*points[0],repeat(fill,-change))
122 )
123 assert(len(addition) == 4, "Fill is the wrong length")
124 result + repeat(addition, len(result));
125
126
127
128// Section: Coordinate Systems
129
130// Function: polar_to_xy()
131// Usage:
132// pt = polar_to_xy(r, theta);
133// pt = polar_to_xy([r, theta]);
134// Topics: Coordinates, Points, Paths
135// See Also: xy_to_polar(), xyz_to_cylindrical(), cylindrical_to_xyz(), xyz_to_spherical(), spherical_to_xyz()
136// Description:
137// Convert polar coordinates to 2D cartesian coordinates.
138// Returns [X,Y] cartesian coordinates.
139// Arguments:
140// r = distance from the origin.
141// theta = angle in degrees, counter-clockwise of X+.
142// Example:
143// xy = polar_to_xy(20,45); // Returns: ~[14.1421365, 14.1421365]
144// xy = polar_to_xy(40,30); // Returns: ~[34.6410162, 15]
145// xy = polar_to_xy([40,30]); // Returns: ~[34.6410162, 15]
146// Example(2D):
147// r=40; ang=30; $fn=36;
148// pt = polar_to_xy(r,ang);
149// stroke(circle(r=r), closed=true, width=0.5);
150// color("black") stroke([[r,0], [0,0], pt], width=0.5);
151// color("black") stroke(arc(r=15, angle=ang), width=0.5);
152// color("red") move(pt) circle(d=3);
153function polar_to_xy(r,theta=undef) = let(
154 rad = theta==undef? r[0] : r,
155 t = theta==undef? r[1] : theta
156 ) rad*[cos(t), sin(t)];
157
158
159// Function: xy_to_polar()
160// Usage:
161// r_theta = xy_to_polar(x,y);
162// r_theta = xy_to_polar([X,Y]);
163// Topics: Coordinates, Points, Paths
164// See Also: polar_to_xy(), xyz_to_cylindrical(), cylindrical_to_xyz(), xyz_to_spherical(), spherical_to_xyz()
165// Description:
166// Convert 2D cartesian coordinates to polar coordinates.
167// Returns [radius, theta] where theta is the angle counter-clockwise of X+.
168// Arguments:
169// x = X coordinate.
170// y = Y coordinate.
171// Example:
172// plr = xy_to_polar(20,30);
173// plr = xy_to_polar([40,60]);
174// Example(2D):
175// pt = [-20,30]; $fn = 36;
176// rt = xy_to_polar(pt);
177// r = rt[0]; ang = rt[1];
178// stroke(circle(r=r), closed=true, width=0.5);
179// zrot(ang) stroke([[0,0],[r,0]],width=0.5);
180// color("red") move(pt) circle(d=3);
181function xy_to_polar(x,y=undef) = let(
182 xx = y==undef? x[0] : x,
183 yy = y==undef? x[1] : y
184 ) [norm([xx,yy]), atan2(yy,xx)];
185
186
187// Function: project_plane()
188// Usage:
189// xy = project_plane(plane, p);
190// Usage: To get a transform matrix
191// M = project_plane(plane)
192// Description:
193// Maps the provided 3d point(s) from 3D coordinates to a 2d coordinate system defined by `plane`. Points that are not
194// on the specified plane will be projected orthogonally onto the plane. This coordinate system is useful if you need
195// to perform 2d operations on a coplanar set of data. After those operations are done you can return the data
196// to 3d with `lift_plane()`. You could also use this to force approximately coplanar data to be exactly coplanar.
197// The parameter p can be a point, path, region, bezier patch or VNF.
198// The plane can be specified as
199// - 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.
200// - A list of coplanar points that define a plane (not-collinear)
201// - 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.
202// .
203// If you omit the point specification then `project_plane()` returns a rotation matrix that maps the specified plane to the XY plane.
204// Note that if you apply this transformation to data lying on the plane it will produce 3D points with the Z coordinate of zero.
205// Topics: Coordinates, Points, Paths
206// Arguments:
207// plane = plane specification or point list defining the plane
208// p = 3D point, path, region, VNF or bezier patch to project
209// Example:
210// pt = [5,-5,5];
211// a=[0,0,0]; b=[10,-10,0]; c=[10,0,10];
212// xy = project_plane([a,b,c],pt);
213// Example(3D): The yellow points in 3D project onto the red points in 2D
214// M = [[-1, 2, -1, -2], [-1, -3, 2, -1], [2, 3, 4, 53], [0, 0, 0, 1]];
215// data = apply(M,path3d(circle(r=10, $fn=20)));
216// move_copies(data) sphere(r=1);
217// color("red") move_copies(project_plane(data, data)) sphere(r=1);
218// Example:
219// xyzpath = move([10,20,30], p=yrot(25, p=path3d(circle(d=100))));
220// mat = project_plane(xyzpath);
221// xypath = path2d(apply(mat, xyzpath));
222// #stroke(xyzpath,closed=true);
223// stroke(xypath,closed=true);
224function project_plane(plane,p) =
225 is_matrix(plane,3,3) && is_undef(p) ? // no data, 3 points given
226 assert(!is_collinear(plane),"Points defining the plane must not be collinear")
227 let(
228 v = plane[2]-plane[0],
229 y = unit(plane[1]-plane[0]), // y axis goes to point b
230 x = unit(v-(v*y)*y) // x axis
231 )
232 frame_map(x,y) * move(-plane[0])
233 : is_vector(plane,4) && is_undef(p) ? // no data, plane given in "plane"
234 assert(_valid_plane(plane), "Plane is not valid")
235 let(
236 n = point3d(plane),
237 cp = n * plane[3] / (n*n)
238 )
239 rot(from=n, to=UP) * move(-cp)
240 : is_path(plane,3) && is_undef(p) ? // no data, generic point list plane
241 assert(len(plane)>=3, "Need three points to define a plane")
242 let(plane = plane_from_points(plane))
243 assert(is_def(plane), "Point list is not coplanar")
244 project_plane(plane)
245 : assert(is_def(p), str("Invalid plane specification: ",plane))
246 is_vnf(p) ? [project_plane(plane,p[0]), p[1]]
247 : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ? // bezier patch or region
248 [for(plist=p) project_plane(plane,plist)]
249 : assert(is_vector(p,3) || is_path(p,3),str("Data must be a 3d point, path, region, vnf or bezier patch",p))
250 is_matrix(plane,3,3) ?
251 assert(!is_collinear(plane),"Points defining the plane must not be collinear")
252 let(
253 v = plane[2]-plane[0],
254 y = unit(plane[1]-plane[0]), // y axis goes to point b
255 x = unit(v-(v*y)*y) // x axis
256 ) move(-plane[0],p) * transpose([x,y])
257 : is_vector(p) ? point2d(apply(project_plane(plane),p))
258 : path2d(apply(project_plane(plane),p));
259
260
261
262// Function: lift_plane()
263// Usage:
264// xyz = lift_plane(plane, p);
265// Usage: to get transform matrix
266// M = lift_plane(plane);
267// Topics: Coordinates, Points, Paths
268// See Also: project_plane()
269// Description:
270// Converts the given 2D point on the plane to 3D coordinates of the specified plane.
271// The parameter p can be a point, path, region, bezier patch or VNF.
272// The plane can be specified as
273// - 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.
274// - A list of coplanar points that define a plane (not-collinear)
275// - 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.
276// 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.
277// This matrix is a rotation, the inverse of the one produced by project_plane.
278// Arguments:
279// plane = Plane specification or list of points to define a plane
280// p = points, path, region, VNF, or bezier patch to transform.
281function lift_plane(plane, p) =
282 is_matrix(plane,3,3) && is_undef(p) ? // no data, 3 p given
283 let(
284 v = plane[2]-plane[0],
285 y = unit(plane[1]-plane[0]), // y axis goes to point b
286 x = unit(v-(v*y)*y) // x axis
287 )
288 move(plane[0]) * frame_map(x,y,reverse=true)
289 : is_vector(plane,4) && is_undef(p) ? // no data, plane given in "plane"
290 assert(_valid_plane(plane), "Plane is not valid")
291 let(
292 n = point3d(plane),
293 cp = n * plane[3] / (n*n)
294 )
295 move(cp) * rot(from=UP, to=n)
296 : is_path(plane,3) && is_undef(p) ? // no data, generic point list plane
297 assert(len(plane)>=3, "Need three p to define a plane")
298 let(plane = plane_from_points(plane))
299 assert(is_def(plane), "Point list is not coplanar")
300 lift_plane(plane)
301 : is_vnf(p) ? [lift_plane(plane,p[0]), p[1]]
302 : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ? // bezier patch or region
303 [for(plist=p) lift_plane(plane,plist)]
304 : assert(is_vector(p,2) || is_path(p,2),"Data must be a 2d point, path, region, vnf or bezier patch")
305 is_matrix(plane,3,3) ?
306 let(
307 v = plane[2]-plane[0],
308 y = unit(plane[1]-plane[0]), // y axis goes to point b
309 x = unit(v-(v*y)*y) // x axis
310 ) move(plane[0],p * [x,y])
311 : apply(lift_plane(plane),is_vector(p) ? point3d(p) : path3d(p));
312
313
314// Function: cylindrical_to_xyz()
315// Usage:
316// pt = cylindrical_to_xyz(r, theta, z);
317// pt = cylindrical_to_xyz([r, theta, z]);
318// Topics: Coordinates, Points, Paths
319// See Also: xyz_to_cylindrical(), xyz_to_spherical(), spherical_to_xyz()
320// Description:
321// Convert cylindrical coordinates to 3D cartesian coordinates. Returns [X,Y,Z] cartesian coordinates.
322// Arguments:
323// r = distance from the Z axis.
324// theta = angle in degrees, counter-clockwise of X+ on the XY plane.
325// z = Height above XY plane.
326// Example:
327// xyz = cylindrical_to_xyz(20,30,40);
328// xyz = cylindrical_to_xyz([40,60,50]);
329function cylindrical_to_xyz(r,theta=undef,z=undef) = let(
330 rad = theta==undef? r[0] : r,
331 t = theta==undef? r[1] : theta,
332 zed = theta==undef? r[2] : z
333 ) [rad*cos(t), rad*sin(t), zed];
334
335
336// Function: xyz_to_cylindrical()
337// Usage:
338// rtz = xyz_to_cylindrical(x,y,z);
339// rtz = xyz_to_cylindrical([X,Y,Z]);
340// Topics: Coordinates, Points, Paths
341// See Also: cylindrical_to_xyz(), xyz_to_spherical(), spherical_to_xyz()
342// Description:
343// Convert 3D cartesian coordinates to cylindrical coordinates. Returns [radius,theta,Z].
344// Theta is the angle counter-clockwise of X+ on the XY plane. Z is height above the XY plane.
345// Arguments:
346// x = X coordinate.
347// y = Y coordinate.
348// z = Z coordinate.
349// Example:
350// cyl = xyz_to_cylindrical(20,30,40);
351// cyl = xyz_to_cylindrical([40,50,70]);
352function xyz_to_cylindrical(x,y=undef,z=undef) = let(
353 p = is_num(x)? [x, default(y,0), default(z,0)] : point3d(x)
354 ) [norm([p.x,p.y]), atan2(p.y,p.x), p.z];
355
356
357// Function: spherical_to_xyz()
358// Usage:
359// pt = spherical_to_xyz(r, theta, phi);
360// pt = spherical_to_xyz([r, theta, phi]);
361// Description:
362// Convert spherical coordinates to 3D cartesian coordinates. Returns [X,Y,Z] cartesian coordinates.
363// Topics: Coordinates, Points, Paths
364// See Also: cylindrical_to_xyz(), xyz_to_spherical(), xyz_to_cylindrical()
365// Arguments:
366// r = distance from origin.
367// theta = angle in degrees, counter-clockwise of X+ on the XY plane.
368// phi = angle in degrees from the vertical Z+ axis.
369// Example:
370// xyz = spherical_to_xyz(20,30,40);
371// xyz = spherical_to_xyz([40,60,50]);
372function spherical_to_xyz(r,theta=undef,phi=undef) = let(
373 rad = theta==undef? r[0] : r,
374 t = theta==undef? r[1] : theta,
375 p = theta==undef? r[2] : phi
376 ) rad*[sin(p)*cos(t), sin(p)*sin(t), cos(p)];
377
378
379// Function: xyz_to_spherical()
380// Usage:
381// r_theta_phi = xyz_to_spherical(x,y,z)
382// r_theta_phi = xyz_to_spherical([X,Y,Z])
383// Topics: Coordinates, Points, Paths
384// See Also: cylindrical_to_xyz(), spherical_to_xyz(), xyz_to_cylindrical()
385// Description:
386// Convert 3D cartesian coordinates to spherical coordinates. Returns [r,theta,phi], where phi is
387// the angle from the Z+ pole, and theta is degrees counter-clockwise of X+ on the XY plane.
388// Arguments:
389// x = X coordinate.
390// y = Y coordinate.
391// z = Z coordinate.
392// Example:
393// sph = xyz_to_spherical(20,30,40);
394// sph = xyz_to_spherical([40,50,70]);
395function xyz_to_spherical(x,y=undef,z=undef) = let(
396 p = is_num(x)? [x, default(y,0), default(z,0)] : point3d(x)
397 ) [norm(p), atan2(p.y,p.x), atan2(norm([p.x,p.y]),p.z)];
398
399
400// Function: altaz_to_xyz()
401// Usage:
402// pt = altaz_to_xyz(alt, az, r);
403// pt = altaz_to_xyz([alt, az, r]);
404// Topics: Coordinates, Points, Paths
405// See Also: cylindrical_to_xyz(), xyz_to_spherical(), spherical_to_xyz(), xyz_to_cylindrical(), xyz_to_altaz()
406// Description:
407// Convert altitude/azimuth/range coordinates to 3D cartesian coordinates.
408// Returns [X,Y,Z] cartesian coordinates.
409// Arguments:
410// alt = altitude angle in degrees above the XY plane.
411// az = azimuth angle in degrees clockwise of Y+ on the XY plane.
412// r = distance from origin.
413// Example:
414// xyz = altaz_to_xyz(20,30,40);
415// xyz = altaz_to_xyz([40,60,50]);
416function altaz_to_xyz(alt,az=undef,r=undef) = let(
417 p = az==undef? alt[0] : alt,
418 t = 90 - (az==undef? alt[1] : az),
419 rad = az==undef? alt[2] : r
420 ) rad*[cos(p)*cos(t), cos(p)*sin(t), sin(p)];
421
422
423// Function: xyz_to_altaz()
424// Usage:
425// alt_az_r = xyz_to_altaz(x,y,z);
426// alt_az_r = xyz_to_altaz([X,Y,Z]);
427// Topics: Coordinates, Points, Paths
428// See Also: cylindrical_to_xyz(), xyz_to_spherical(), spherical_to_xyz(), xyz_to_cylindrical(), altaz_to_xyz()
429// Description:
430// Convert 3D cartesian coordinates to altitude/azimuth/range coordinates.
431// Returns [altitude,azimuth,range], where altitude is angle above the
432// XY plane, azimuth is degrees clockwise of Y+ on the XY plane, and
433// range is the distance from the origin.
434// Arguments:
435// x = X coordinate.
436// y = Y coordinate.
437// z = Z coordinate.
438// Example:
439// aa = xyz_to_altaz(20,30,40);
440// aa = xyz_to_altaz([40,50,70]);
441function xyz_to_altaz(x,y=undef,z=undef) = let(
442 p = is_num(x)? [x, default(y,0), default(z,0)] : point3d(x)
443 ) [atan2(p.z,norm([p.x,p.y])), atan2(p.x,p.y), norm(p)];
444
445
446
447// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap