1//////////////////////////////////////////////////////////////////////
2// LibFile: geometry.scad
3// Perform calculations on lines, polygons, planes and circles, including
4// normals, intersections of objects, distance between objects, and tangent lines.
5// Throughout this library, lines can be treated as either unbounded lines, as rays with
6// a single endpoint or as segments, bounded by endpoints at both ends.
7// Includes:
8// include <BOSL2/std.scad>
9// FileGroup: Math
10// FileSummary: Geometrical calculations including intersections of lines, circles and planes, circle from 3 points
11// FileFootnotes: STD=Included in std.scad
12//////////////////////////////////////////////////////////////////////
13
14
15// Section: Lines, Rays, and Segments
16
17// Function: is_point_on_line()
18// Usage:
19// pt = is_point_on_line(point, line, [bounded], [eps]);
20// Topics: Geometry, Points, Segments
21// Description:
22// Determine if the point is on the line segment, ray or segment defined by the two between two points.
23// Returns true if yes, and false if not. If bounded is set to true it specifies a segment, with
24// both lines bounded at the ends. Set bounded to `[true,false]` to get a ray. You can use
25// the shorthands RAY and SEGMENT to set bounded.
26// Arguments:
27// point = The point to test.
28// line = Array of two points defining the line, ray, or segment to test against.
29// bounded = boolean or list of two booleans defining endpoint conditions for the line. If false treat the line as an unbounded line. If true treat it as a segment. If [true,false] treat as a ray, based at the first endpoint. Default: false
30// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
31function is_point_on_line(point, line, bounded=false, eps=EPSILON) =
32 assert(is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
33 assert(is_vector(point), "Point must be a vector")
34 assert(_valid_line(line, len(point),eps),"Given line is not valid")
35 _is_point_on_line(point, line, bounded,eps);
36
37function _is_point_on_line(point, line, bounded=false, eps=EPSILON) =
38 let(
39 v1 = (line[1]-line[0]),
40 v0 = (point-line[0]),
41 t = v0*v1/(v1*v1),
42 bounded = force_list(bounded,2)
43 )
44 abs(cross(v0,v1))<=eps*norm(v1)
45 && (!bounded[0] || t>=-eps)
46 && (!bounded[1] || t<1+eps) ;
47
48
49///Internal - distance from point `d` to the line passing through the origin with unit direction n
50///_dist2line works for any dimension
51function _dist2line(d,n) = norm(d-(d * n) * n);
52
53
54///Internal
55function _valid_line(line,dim,eps=EPSILON) =
56 is_matrix(line,2,dim)
57 && norm(line[1]-line[0])>eps*max(norm(line[1]),norm(line[0]));
58
59//Internal
60function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
61
62
63/// Internal Function: _is_at_left()
64/// Usage:
65/// pt = point_left_of_line2d(point, line);
66/// Topics: Geometry, Points, Lines
67/// Description:
68/// Return true iff a 2d point is on or at left of the line defined by `line`.
69/// Arguments:
70/// pt = The 2d point to check position of.
71/// line = Array of two 2d points forming the line segment to test against.
72/// eps = Tolerance in the geometrical tests.
73function _is_at_left(pt,line,eps=EPSILON) = _tri_class([pt,line[0],line[1]],eps) <= 0;
74
75
76/// Internal Function: _degenerate_tri()
77/// Usage:
78/// degen = _degenerate_tri(triangle);
79/// Topics: Geometry, Triangles
80/// Description:
81/// Return true for a specific kind of degeneracy: any two triangle vertices are equal
82/// Arguments:
83/// tri = A list of three 2d points
84/// eps = Tolerance in the geometrical tests.
85function _degenerate_tri(tri,eps) =
86 max(norm(tri[0]-tri[1]), norm(tri[1]-tri[2]), norm(tri[2]-tri[0])) < eps ;
87
88
89/// Internal Function: _tri_class()
90/// Usage:
91/// class = _tri_class(triangle);
92/// Topics: Geometry, Triangles
93/// Description:
94/// Return 1 if the triangle `tri` is CW.
95/// Return 0 if the triangle `tri` has colinear vertices.
96/// Return -1 if the triangle `tri` is CCW.
97/// Arguments:
98/// tri = A list of the three 2d vertices of a triangle.
99/// eps = Tolerance in the geometrical tests.
100function _tri_class(tri, eps=EPSILON) =
101 let( crx = cross(tri[1]-tri[2],tri[0]-tri[2]) )
102 abs( crx ) <= eps*norm(tri[1]-tri[2])*norm(tri[0]-tri[2]) ? 0 : sign( crx );
103
104
105/// Internal Function: _pt_in_tri()
106/// Usage:
107/// class = _pt_in_tri(point, tri);
108/// Topics: Geometry, Points, Triangles
109/// Description:
110// For CW triangles `tri` :
111/// return 1 if point is inside the triangle interior.
112/// return =0 if point is on the triangle border.
113/// return -1 if point is outside the triangle.
114/// Arguments:
115/// point = The point to check position of.
116/// tri = A list of the three 2d vertices of a triangle.
117/// eps = Tolerance in the geometrical tests.
118function _pt_in_tri(point, tri, eps=EPSILON) =
119 min( _tri_class([tri[0],tri[1],point],eps),
120 _tri_class([tri[1],tri[2],point],eps),
121 _tri_class([tri[2],tri[0],point],eps) );
122
123
124/// Internal Function: _point_left_of_line2d()
125/// Usage:
126/// pt = point_left_of_line2d(point, line);
127/// Topics: Geometry, Points, Lines
128/// Description:
129/// Return >0 if point is left of the line defined by `line`.
130/// Return =0 if point is on the line.
131/// Return <0 if point is right of the line.
132/// Arguments:
133/// point = The point to check position of.
134/// line = Array of two points forming the line segment to test against.
135function _point_left_of_line2d(point, line, eps=EPSILON) =
136 assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." )
137// cross(line[0]-point, line[1]-line[0]);
138 _tri_class([point,line[1],line[0]],eps);
139
140
141// Function: is_collinear()
142// Usage:
143// bool = is_collinear(a, [b, c], [eps]);
144// Topics: Geometry, Points, Collinearity
145// Description:
146// Returns true if the points `a`, `b` and `c` are co-linear or if the list of points `a` is collinear.
147// Arguments:
148// a = First point or list of points.
149// b = Second point or undef; it should be undef if `c` is undef
150// c = Third point or undef.
151// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
152function is_collinear(a, b, c, eps=EPSILON) =
153 assert( is_path([a,b,c],dim=undef)
154 || ( is_undef(b) && is_undef(c) && is_path(a,dim=undef) ),
155 "Input should be 3 points or a list of points with same dimension.")
156 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
157 let( points = is_def(c) ? [a,b,c]: a )
158 len(points)<3 ? true :
159 _noncollinear_triple(points,error=false,eps=eps) == [];
160
161
162// Function: point_line_distance()
163// Usage:
164// dist = point_line_distance(line, pt, [bounded]);
165// Topics: Geometry, Points, Lines, Distance
166// Description:
167// Finds the shortest distance from the point `pt` to the specified line, segment or ray.
168// The bounded parameter specifies the whether the endpoints give a ray or segment.
169// By default assumes an unbounded line.
170// Arguments:
171// line = A list of two points defining a line.
172// pt = A point to find the distance of from the line.
173// bounded = a boolean or list of two booleans specifiying whether each end is bounded. Default: false
174// Example:
175// dist1 = point_line_distance([3,8], [[-10,0], [10,0]]); // Returns: 8
176// dist2 = point_line_distance([3,8], [[-10,0], [10,0]],SEGMENT); // Returns: 8
177// dist3 = point_line_distance([14,3], [[-10,0], [10,0]],SEGMENT); // Returns: 5
178function point_line_distance(pt, line, bounded=false) =
179 assert(is_bool(bounded) || is_bool_list(bounded,2), "\"bounded\" is invalid")
180 assert( _valid_line(line) && is_vector(pt,len(line[0])),
181 "Invalid line, invalid point or incompatible dimensions." )
182 bounded == LINE ? _dist2line(pt-line[0],unit(line[1]-line[0]))
183 : norm(pt-line_closest_point(line,pt,bounded));
184
185
186// Function: segment_distance()
187// Usage:
188// dist = segment_distance(seg1, seg2, [eps]);
189// Topics: Geometry, Segments, Distance
190// See Also: convex_collision(), convex_distance()
191// Description:
192// Returns the closest distance of the two given line segments.
193// Arguments:
194// seg1 = The list of two points representing the first line segment to check the distance of.
195// seg2 = The list of two points representing the second line segment to check the distance of.
196// eps = tolerance for point comparisons
197// Example:
198// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5
199// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0
200function segment_distance(seg1, seg2,eps=EPSILON) =
201 assert( is_matrix(concat(seg1,seg2),4), "Inputs should be two valid segments." )
202 convex_distance(seg1,seg2,eps);
203
204
205// Function: line_normal()
206// Usage:
207// vec = line_normal([P1,P2])
208// vec = line_normal(p1,p2)
209// Topics: Geometry, Lines
210// Description:
211// Returns the 2D normal vector to the given 2D line. This is otherwise known as the perpendicular vector counter-clockwise to the given ray.
212// Arguments:
213// p1 = First point on 2D line.
214// p2 = Second point on 2D line.
215// Example(2D):
216// p1 = [10,10];
217// p2 = [50,30];
218// n = line_normal(p1,p2);
219// stroke([p1,p2], endcap2="arrow2");
220// color("green") stroke([p1,p1+10*n], endcap2="arrow2");
221// color("blue") move_copies([p1,p2]) circle(d=2, $fn=12);
222function line_normal(p1,p2) =
223 is_undef(p2)
224 ? assert( len(p1)==2 && !is_undef(p1[1]) , "Invalid input." )
225 line_normal(p1[0],p1[1])
226 : assert( _valid_line([p1,p2],dim=2), "Invalid line." )
227 unit([p1.y-p2.y,p2.x-p1.x]);
228
229
230// 2D Line intersection from two segments.
231// This function returns [p,t,u] where p is the intersection point of
232// the lines defined by the two segments, t is the proportional distance
233// of the intersection point along s1, and u is the proportional distance
234// of the intersection point along s2. The proportional values run over
235// the range of 0 to 1 for each segment, so if it is in this range, then
236// the intersection lies on the segment. Otherwise it lies somewhere on
237// the extension of the segment. If lines are parallel or coincident then
238// it returns undef.
239
240function _general_line_intersection(s1,s2,eps=EPSILON) =
241 let(
242 denominator = cross(s1[0]-s1[1],s2[0]-s2[1])
243 )
244 approx(denominator,0,eps=eps) ? undef :
245 let(
246 t = cross(s1[0]-s2[0],s2[0]-s2[1]) / denominator,
247 u = cross(s1[0]-s2[0],s1[0]-s1[1]) / denominator
248 )
249 [s1[0]+t*(s1[1]-s1[0]), t, u];
250
251
252// Function: line_intersection()
253// Usage:
254// pt = line_intersection(line1, line2, [bounded1], [bounded2], [bounded=], [eps=]);
255// Description:
256// Returns the intersection point of any two 2D lines, segments or rays. Returns undef
257// if they do not intersect. You specify a line by giving two distinct points on the
258// line. You specify rays or segments by giving a pair of points and indicating
259// bounded[0]=true to bound the line at the first point, creating rays based at l1[0] and l2[0],
260// or bounded[1]=true to bound the line at the second point, creating the reverse rays bounded
261// at l1[1] and l2[1]. If bounded=[true, true] then you have segments defined by their two
262// endpoints. By using bounded1 and bounded2 you can mix segments, rays, and lines as needed.
263// You can set the bounds parameters to true as a shorthand for [true,true] to sepcify segments.
264// Arguments:
265// line1 = List of two points in 2D defining the first line, segment or ray
266// line2 = List of two points in 2D defining the second line, segment or ray
267// bounded1 = boolean or list of two booleans defining which ends are bounded for line1. Default: [false,false]
268// bounded2 = boolean or list of two booleans defining which ends are bounded for line2. Default: [false,false]
269// ---
270// bounded = boolean or list of two booleans defining which ends are bounded for both lines. The bounded1 and bounded2 parameters override this if both are given.
271// eps = tolerance for geometric comparisons. Default: `EPSILON` (1e-9)
272// Example(2D): The segments do not intersect but the lines do in this example.
273// line1 = 10*[[9, 4], [5, 7]];
274// line2 = 10*[[2, 3], [6, 5]];
275// stroke(line1, endcaps="arrow2");
276// stroke(line2, endcaps="arrow2");
277// isect = line_intersection(line1, line2);
278// color("red") translate(isect) circle(r=1,$fn=12);
279// Example(2D): Specifying a ray and segment using the shorthand variables.
280// line1 = 10*[[0, 2], [4, 7]];
281// line2 = 10*[[10, 4], [3, 4]];
282// stroke(line1);
283// stroke(line2, endcap2="arrow2");
284// isect = line_intersection(line1, line2, SEGMENT, RAY);
285// color("red") translate(isect) circle(r=1,$fn=12);
286// Example(2D): Here we use the same example as above, but specify two segments using the bounded argument.
287// line1 = 10*[[0, 2], [4, 7]];
288// line2 = 10*[[10, 4], [3, 4]];
289// stroke(line1);
290// stroke(line2);
291// isect = line_intersection(line1, line2, bounded=true); // Returns undef
292function line_intersection(line1, line2, bounded1, bounded2, bounded, eps=EPSILON) =
293 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
294 assert( _valid_line(line1,dim=2,eps=eps), "First line invalid")
295 assert( _valid_line(line2,dim=2,eps=eps), "Second line invalid")
296 assert( is_undef(bounded) || is_bool(bounded) || is_bool_list(bounded,2), "Invalid value for \"bounded\"")
297 assert( is_undef(bounded1) || is_bool(bounded1) || is_bool_list(bounded1,2), "Invalid value for \"bounded1\"")
298 assert( is_undef(bounded2) || is_bool(bounded2) || is_bool_list(bounded2,2), "Invalid value for \"bounded2\"")
299 let(isect = _general_line_intersection(line1,line2,eps=eps))
300 is_undef(isect) ? undef :
301 let(
302 bounded1 = force_list(first_defined([bounded1,bounded,false]),2),
303 bounded2 = force_list(first_defined([bounded2,bounded,false]),2),
304 good = (!bounded1[0] || isect[1]>=0-eps)
305 && (!bounded1[1] || isect[1]<=1+eps)
306 && (!bounded2[0] || isect[2]>=0-eps)
307 && (!bounded2[1] || isect[2]<=1+eps)
308 )
309 good ? isect[0] : undef;
310
311
312// Function: line_closest_point()
313// Usage:
314// pt = line_closest_point(line, pt, [bounded]);
315// Topics: Geometry, Lines, Distance
316// Description:
317// Returns the point on the given line, segment or ray that is closest to the given point `pt`.
318// The inputs `line` and `pt` args should be of the same dimension. The parameter bounded indicates
319// whether the points of `line` should be treated as endpoints.
320// Arguments:
321// line = A list of two points that are on the unbounded line.
322// pt = The point to find the closest point on the line to.
323// bounded = boolean or list of two booleans indicating that the line is bounded at that end. Default: [false,false]
324// Example(2D):
325// line = [[-30,0],[30,30]];
326// pt = [-32,-10];
327// p2 = line_closest_point(line,pt);
328// stroke(line, endcaps="arrow2");
329// color("blue") translate(pt) circle(r=1,$fn=12);
330// color("red") translate(p2) circle(r=1,$fn=12);
331// Example(2D): If the line is bounded on the left you get the endpoint instead
332// line = [[-30,0],[30,30]];
333// pt = [-32,-10];
334// p2 = line_closest_point(line,pt,bounded=[true,false]);
335// stroke(line, endcap2="arrow2");
336// color("blue") translate(pt) circle(r=1,$fn=12);
337// color("red") translate(p2) circle(r=1,$fn=12);
338// Example(2D): In this case it doesn't matter how bounded is set. Using SEGMENT is the most restrictive option.
339// line = [[-30,0],[30,30]];
340// pt = [-5,0];
341// p2 = line_closest_point(line,pt,SEGMENT);
342// stroke(line);
343// color("blue") translate(pt) circle(r=1,$fn=12);
344// color("red") translate(p2) circle(r=1,$fn=12);
345// Example(2D): The result here is the same for a line or a ray.
346// line = [[-30,0],[30,30]];
347// pt = [40,25];
348// p2 = line_closest_point(line,pt,RAY);
349// stroke(line, endcap2="arrow2");
350// color("blue") translate(pt) circle(r=1,$fn=12);
351// color("red") translate(p2) circle(r=1,$fn=12);
352// Example(2D): But with a segment we get a different result
353// line = [[-30,0],[30,30]];
354// pt = [40,25];
355// p2 = line_closest_point(line,pt,SEGMENT);
356// stroke(line);
357// color("blue") translate(pt) circle(r=1,$fn=12);
358// color("red") translate(p2) circle(r=1,$fn=12);
359// Example(2D): The shorthand RAY uses the first point as the base of the ray. But you can specify a reversed ray directly, and in this case the result is the same as the result above for the segment.
360// line = [[-30,0],[30,30]];
361// pt = [40,25];
362// p2 = line_closest_point(line,pt,[false,true]);
363// stroke(line,endcap1="arrow2");
364// color("blue") translate(pt) circle(r=1,$fn=12);
365// color("red") translate(p2) circle(r=1,$fn=12);
366// Example(FlatSpin,VPD=200,VPT=[0,0,15]): A 3D example
367// line = [[-30,-15,0],[30,15,30]];
368// pt = [5,5,5];
369// p2 = line_closest_point(line,pt);
370// stroke(line, endcaps="arrow2");
371// color("blue") translate(pt) sphere(r=1,$fn=12);
372// color("red") translate(p2) sphere(r=1,$fn=12);
373function line_closest_point(line, pt, bounded=false) =
374 assert(_valid_line(line), "Invalid line")
375 assert(is_vector(pt, len(line[0])), "Invalid point or incompatible dimensions.")
376 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid value for \"bounded\"")
377 let(
378 bounded = force_list(bounded,2)
379 )
380 bounded==[false,false] ?
381 let( n = unit( line[0]- line[1]) )
382 line[1] + ((pt- line[1]) * n) * n
383 : bounded == [true,true] ?
384 pt + _closest_s1([line[0]-pt, line[1]-pt])[0]
385 :
386 let(
387 ray = bounded==[true,false] ? line : reverse(line),
388 seglen = norm(ray[1]-ray[0]),
389 segvec = (ray[1]-ray[0])/seglen,
390 projection = (pt-ray[0]) * segvec
391 )
392 projection<=0 ? ray[0] :
393 ray[0] + projection*segvec;
394
395
396// Function: line_from_points()
397// Usage:
398// line = line_from_points(points, [fast], [eps]);
399// Topics: Geometry, Lines, Points
400// Description:
401// Given a list of 2 or more collinear points, returns a line containing them.
402// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned.
403// if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
404// Arguments:
405// points = The list of points to find the line through.
406// fast = If true, don't verify that all points are collinear. Default: false
407// eps = How much variance is allowed in testing each point against the line. Default: `EPSILON` (1e-9)
408function line_from_points(points, fast=false, eps=EPSILON) =
409 assert( is_path(points), "Invalid point list." )
410 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
411 let( pb = furthest_point(points[0],points) )
412 norm(points[pb]-points[0])<eps*max(norm(points[pb]),norm(points[0])) ? undef :
413 fast || is_collinear(points)
414 ? [points[pb], points[0]]
415 : undef;
416
417
418
419// Section: Planes
420
421
422// Function: is_coplanar()
423// Usage:
424// bool = is_coplanar(points,[eps]);
425// Topics: Geometry, Coplanarity
426// Description:
427// Returns true if the given 3D points are non-collinear and are on a plane.
428// Arguments:
429// points = The points to test.
430// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
431function is_coplanar(points, eps=EPSILON) =
432 assert( is_path(points,dim=3) , "Input should be a list of 3D points." )
433 assert( is_finite(eps) && eps>=0, "The tolerance should be a non-negative value." )
434 len(points)<=2 ? false
435 : let( ip = _noncollinear_triple(points,error=false,eps=eps) )
436 ip == [] ? false :
437 let( plane = plane3pt(points[ip[0]],points[ip[1]],points[ip[2]]) )
438 _pointlist_greatest_distance(points,plane) < eps;
439
440
441
442// Function: plane3pt()
443// Usage:
444// plane = plane3pt(p1, p2, p3);
445// Topics: Geometry, Planes
446// Description:
447// Generates the normalized cartesian equation of a plane from three 3d points.
448// Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane.
449// Returns undef, if the points are collinear.
450// Arguments:
451// p1 = The first point on the plane.
452// p2 = The second point on the plane.
453// p3 = The third point on the plane.
454function plane3pt(p1, p2, p3) =
455 is_undef(p2) && is_undef(p3) && is_path(p1,dim=3) ? plane3pt(p1[0],p1[1],p1[2])
456 : assert( is_path([p1,p2,p3],dim=3) && len(p1)==3,
457 "Invalid points or incompatible dimensions." )
458 let(
459 crx = cross(p3-p1, p2-p1),
460 nrm = norm(crx)
461 ) approx(nrm,0) ? undef :
462 concat(crx, crx*p1)/nrm;
463
464
465// Function: plane3pt_indexed()
466// Usage:
467// plane = plane3pt_indexed(points, i1, i2, i3);
468// Topics: Geometry, Planes
469// Description:
470// Given a list of 3d points, and the indices of three of those points,
471// generates the normalized cartesian equation of a plane that those points all
472// lie on. If the points are not collinear, returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane.
473// If they are collinear, returns [].
474// Arguments:
475// points = A list of points.
476// i1 = The index into `points` of the first point on the plane.
477// i2 = The index into `points` of the second point on the plane.
478// i3 = The index into `points` of the third point on the plane.
479function plane3pt_indexed(points, i1, i2, i3) =
480 is_undef(i3) && is_undef(i2) && is_vector(i1) ? plane3pt_indexed(points, i1[0], i1[1], i1[2])
481 :
482 assert( is_vector([i1,i2,i3]) && min(i1,i2,i3)>=0 && is_list(points) && max(i1,i2,i3)<len(points),
483 "Invalid or out of range indices." )
484 assert( is_path([points[i1], points[i2], points[i3]],dim=3),
485 "Improper points or improper dimensions." )
486 let(
487 p1 = points[i1],
488 p2 = points[i2],
489 p3 = points[i3]
490 ) plane3pt(p1,p2,p3);
491
492
493// Function: plane_from_normal()
494// Usage:
495// plane = plane_from_normal(normal, [pt])
496// Topics: Geometry, Planes
497// Description:
498// Returns a plane defined by a normal vector and a point. If you omit `pt` you will get a plane
499// passing through the origin.
500// Arguments:
501// normal = Normal vector to the plane to find.
502// pt = Point 3D on the plane to find.
503// Example:
504// plane_from_normal([0,0,1], [2,2,2]); // Returns the xy plane passing through the point (2,2,2)
505function plane_from_normal(normal, pt=[0,0,0]) =
506 assert( is_matrix([normal,pt],2,3) && !approx(norm(normal),0),
507 "Inputs `normal` and `pt` should be 3d vectors/points and `normal` cannot be zero." )
508 concat(normal, normal*pt) / norm(normal);
509
510
511// Eigenvalues for a 3x3 symmetrical matrix in decreasing order
512// Based on: https://en.wikipedia.org/wiki/Eigenvalue_algorithm
513function _eigenvals_symm_3(M) =
514 let( p1 = pow(M[0][1],2) + pow(M[0][2],2) + pow(M[1][2],2) )
515 (p1<EPSILON)
516 ? -sort(-[ M[0][0], M[1][1], M[2][2] ]) // diagonal matrix: eigenvals in decreasing order
517 : let( q = (M[0][0]+M[1][1]+M[2][2])/3,
518 B = (M - q*ident(3)),
519 dB = [B[0][0], B[1][1], B[2][2]],
520 p2 = dB*dB + 2*p1,
521 p = sqrt(p2/6),
522 r = det3(B/p)/2,
523 ph = acos(constrain(r,-1,1))/3,
524 e1 = q + 2*p*cos(ph),
525 e3 = q + 2*p*cos(ph+120),
526 e2 = 3*q - e1 - e3 )
527 [ e1, e2, e3 ];
528
529
530// the i-th normalized eigenvector of a 3x3 symmetrical matrix M from its eigenvalues
531// using Cayley–Hamilton theorem according to:
532// https://en.wikipedia.org/wiki/Eigenvalue_algorithm
533function _eigenvec_symm_3(M,evals,i=0) =
534 let(
535 I = ident(3),
536 A = (M - evals[(i+1)%3]*I) * (M - evals[(i+2)%3]*I) ,
537 k = max_index( [for(i=[0:2]) norm(A[i]) ])
538 )
539 norm(A[k])<EPSILON ? I[k] : A[k]/norm(A[k]);
540
541
542// finds the eigenvector corresponding to the smallest eigenvalue of the covariance matrix of a pointlist
543// returns the mean of the points, the eigenvector and the greatest eigenvalue
544function _covariance_evec_eval(points) =
545 let( pm = sum(points)/len(points), // mean point
546 Y = [ for(i=[0:len(points)-1]) points[i] - pm ],
547 M = transpose(Y)*Y , // covariance matrix
548 evals = _eigenvals_symm_3(M), // eigenvalues in decreasing order
549 evec = _eigenvec_symm_3(M,evals,i=2) )
550 [pm, evec, evals[0] ];
551
552
553// Function: plane_from_points()
554// Usage:
555// plane = plane_from_points(points, [fast], [eps]);
556// Topics: Geometry, Planes, Points
557// See Also: plane_from_polygon()
558// Description:
559// Given a list of 3 or more coplanar 3D points, returns the coefficients of the normalized cartesian equation of a plane,
560// that is [A,B,C,D] where Ax+By+Cz=D is the equation of the plane and norm([A,B,C])=1.
561// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
562// If `fast` is true, the polygon coplanarity check is skipped and a best fitting plane is returned.
563// It differs from `plane_from_polygon` as the plane normal is independent of the point order. It is faster, though.
564// Arguments:
565// points = The list of points to find the plane of.
566// fast = If true, don't verify the point coplanarity. Default: false
567// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
568// Example(3D):
569// points = rot(45, v=[-0.3,1,0], p=path3d(random_points(25,2,scale=55,seed=47), 70));
570// plane = plane_from_points(points);
571// #move_copies(points)sphere(d=3);
572// cp = mean(points);
573// move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(50);
574function plane_from_points(points, fast=false, eps=EPSILON) =
575 assert( is_path(points,dim=3), "Improper 3d point list." )
576 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
577 len(points) == 3
578 ? plane3pt(points[0],points[1],points[2])
579 : let(
580 covmix = _covariance_evec_eval(points),
581 pm = covmix[0],
582 evec = covmix[1],
583 eval0 = covmix[2],
584 plane = [ each evec, pm*evec]
585 )
586 !fast && _pointlist_greatest_distance(points,plane)>eps*eval0 ? undef :
587 plane ;
588
589
590// Function: plane_from_polygon()
591// Usage:
592// plane = plane_from_polygon(points, [fast], [eps]);
593// Topics: Geometry, Planes, Polygons
594// See Also: plane_from_points()
595// Description:
596// Given a 3D planar polygon, returns the normalized cartesian equation of its plane.
597// Returns [A,B,C,D] where Ax+By+Cz=D is the equation of the plane where norm([A,B,C])=1.
598// If not all the points in the polygon are coplanar, then [] is returned.
599// If `fast` is false and the points in the list are collinear or not coplanar, then `undef` is returned.
600// if `fast` is true, then the coplanarity test is skipped and a plane passing through 3 non-collinear arbitrary points is returned.
601// Arguments:
602// poly = The planar 3D polygon to find the plane of.
603// fast = If true, doesn't verify that all points in the polygon are coplanar. Default: false
604// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
605// Example(3D):
606// xyzpath = rot(45, v=[0,1,0], p=path3d(star(n=5,step=2,d=100), 70));
607// plane = plane_from_polygon(xyzpath);
608// #stroke(xyzpath,closed=true,width=3);
609// cp = centroid(xyzpath);
610// move(cp) rot(from=UP,to=plane_normal(plane)) anchor_arrow(45);
611function plane_from_polygon(poly, fast=false, eps=EPSILON) =
612 assert( is_path(poly,dim=3), "Invalid polygon." )
613 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
614 let(
615 poly_normal = polygon_normal(poly)
616 )
617 is_undef(poly_normal) ? undef :
618 let(
619 plane = plane_from_normal(poly_normal, poly[0])
620 )
621 fast? plane: are_points_on_plane(poly, plane, eps=eps)? plane: undef;
622
623
624// Function: plane_normal()
625// Usage:
626// vec = plane_normal(plane);
627// Topics: Geometry, Planes
628// Description:
629// Returns the unit length normal vector for the given plane.
630// Arguments:
631// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
632function plane_normal(plane) =
633 assert( _valid_plane(plane), "Invalid input plane." )
634 unit([plane.x, plane.y, plane.z]);
635
636
637// Function: plane_offset()
638// Usage:
639// d = plane_offset(plane);
640// Topics: Geometry, Planes
641// Description:
642// Returns coeficient D of the normalized plane equation `Ax+By+Cz=D`, or the scalar offset of the plane from the origin.
643// This value may be negative.
644// The absolute value of this coefficient is the distance of the plane from the origin.
645// Arguments:
646// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
647function plane_offset(plane) =
648 assert( _valid_plane(plane), "Invalid input plane." )
649 plane[3]/norm([plane.x, plane.y, plane.z]);
650
651
652
653// Returns [POINT, U] if line intersects plane at one point, where U is zero at line[0] and 1 at line[1]
654// Returns [LINE, undef] if the line is on the plane.
655// Returns undef if line is parallel to, but not on the given plane.
656function _general_plane_line_intersection(plane, line, eps=EPSILON) =
657 let(
658 a = plane*[each line[0],-1], // evaluation of the plane expression at line[0]
659 b = plane*[each(line[1]-line[0]),0] // difference between the plane expression evaluation at line[1] and at line[0]
660 )
661 approx(b,0,eps) // is (line[1]-line[0]) "parallel" to the plane ?
662 ? approx(a,0,eps) // is line[0] on the plane ?
663 ? [line,undef] // line is on the plane
664 : undef // line is parallel but not on the plane
665 : [ line[0]-a/b*(line[1]-line[0]), -a/b ];
666
667
668/// Internal Function: normalize_plane()
669// Usage:
670// nplane = normalize_plane(plane);
671/// Topics: Geometry, Planes
672// Description:
673// Returns a new representation [A,B,C,D] of `plane` where norm([A,B,C]) is equal to one.
674function _normalize_plane(plane) =
675 assert( _valid_plane(plane), str("Invalid plane. ",plane ) )
676 plane/norm(point3d(plane));
677
678
679// Function: plane_line_intersection()
680// Usage:
681// pt = plane_line_intersection(plane, line, [bounded], [eps]);
682// Topics: Geometry, Planes, Lines, Intersection
683// Description:
684// Takes a line, and a plane [A,B,C,D] where the equation of that plane is `Ax+By+Cz=D`.
685// If `line` intersects `plane` at one point, then that intersection point is returned.
686// If `line` lies on `plane`, then the original given `line` is returned.
687// If `line` is parallel to, but not on `plane`, then undef is returned.
688// Arguments:
689// plane = The [A,B,C,D] values for the equation of the plane.
690// line = A list of two distinct 3D points that are on the line.
691// bounded = If false, the line is considered unbounded. If true, it is treated as a bounded line segment. If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray. Default: false (unbounded)
692// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
693function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
694 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
695 assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.")
696 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
697 let(
698 bounded = is_list(bounded)? bounded : [bounded, bounded],
699 res = _general_plane_line_intersection(plane, line, eps=eps)
700 ) is_undef(res) ? undef :
701 is_undef(res[1]) ? res[0] :
702 bounded[0] && res[1]<0 ? undef :
703 bounded[1] && res[1]>1 ? undef :
704 res[0];
705
706
707
708// Function: plane_intersection()
709// Usage:
710// line = plane_intersection(plane1, plane2)
711// pt = plane_intersection(plane1, plane2, plane3)
712// Topics: Geometry, Planes, Intersection
713// Description:
714// Compute the point which is the intersection of the three planes, or the line intersection of two planes.
715// If you give three planes the intersection is returned as a point. If you give two planes the intersection
716// is returned as a list of two points on the line of intersection. If any two input planes are parallel
717// or coincident then returns undef.
718// Arguments:
719// plane1 = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`.
720// plane2 = The [A,B,C,D] coefficients for the second plane equation `Ax+By+Cz=D`.
721// plane3 = The [A,B,C,D] coefficients for the third plane equation `Ax+By+Cz=D`.
722function plane_intersection(plane1,plane2,plane3) =
723 assert( _valid_plane(plane1) && _valid_plane(plane2) && (is_undef(plane3) ||_valid_plane(plane3)),
724 "The input must be 2 or 3 planes." )
725 is_def(plane3)
726 ? let(
727 matrix = [for(p=[plane1,plane2,plane3]) point3d(p)],
728 rhs = [for(p=[plane1,plane2,plane3]) p[3]]
729 )
730 linear_solve(matrix,rhs)
731 : let( normal = cross(plane_normal(plane1), plane_normal(plane2)) )
732 approx(norm(normal),0) ? undef :
733 let(
734 matrix = [for(p=[plane1,plane2]) point3d(p)],
735 rhs = [plane1[3], plane2[3]],
736 point = linear_solve(matrix,rhs)
737 )
738 point==[]? undef:
739 [point, point+normal];
740
741
742
743// Function: plane_line_angle()
744// Usage:
745// angle = plane_line_angle(plane,line);
746// Topics: Geometry, Planes, Lines, Angle
747// Description:
748// Compute the angle between a plane [A, B, C, D] and a 3d line, specified as a pair of 3d points [p1,p2].
749// The resulting angle is signed, with the sign positive if the vector p2-p1 lies above the plane, on
750// the same side of the plane as the plane's normal vector.
751function plane_line_angle(plane, line) =
752 assert( _valid_plane(plane), "Invalid plane." )
753 assert( _valid_line(line,dim=3), "Invalid 3d line." )
754 let(
755 linedir = unit(line[1]-line[0]),
756 normal = plane_normal(plane),
757 sin_angle = linedir*normal,
758 cos_angle = norm(cross(linedir,normal))
759 ) atan2(sin_angle,cos_angle);
760
761
762
763// Function: plane_closest_point()
764// Usage:
765// pts = plane_closest_point(plane, points);
766// Topics: Geometry, Planes, Projection
767// Description:
768// Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, and a list of 2d or
769// 3d points, return the closest 3D orthogonal projection of the points on the plane.
770// In other words, for every point given, returns the closest point to it on the plane.
771// If points is a single point then returns a single point result.
772// Arguments:
773// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
774// points = List of points to project
775// Example(FlatSpin,VPD=500,VPT=[2,20,10]):
776// points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100, $fn=36))));
777// plane = plane_from_normal([1,0,1]);
778// proj = plane_closest_point(plane,points);
779// color("red") move_copies(points) sphere(d=4,$fn=12);
780// color("blue") move_copies(proj) sphere(d=4,$fn=12);
781// move(centroid(proj)) {
782// rot(from=UP,to=plane_normal(plane)) {
783// anchor_arrow(50);
784// %cube([120,150,0.1],center=true);
785// }
786// }
787function plane_closest_point(plane, points) =
788 is_vector(points,3) ? plane_closest_point(plane,[points])[0] :
789 assert( _valid_plane(plane), "Invalid plane." )
790 assert( is_matrix(points,undef,3), "Must supply 3D points.")
791 let(
792 plane = _normalize_plane(plane),
793 n = point3d(plane)
794 )
795 [for(pi=points) pi - (pi*n - plane[3])*n];
796
797
798// Function: point_plane_distance()
799// Usage:
800// dist = point_plane_distance(plane, point)
801// Topics: Geometry, Planes, Distance
802// Description:
803// Given a plane as [A,B,C,D] where the cartesian equation for that plane
804// is Ax+By+Cz=D, determines how far from that plane the given point is.
805// The returned distance will be positive if the point is above the
806// plane, meaning on the side where the plane normal points.
807// If the point is below the plane, then the distance returned
808// will be negative. The normal of the plane is [A,B,C].
809// Arguments:
810// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
811// point = The distance evaluation point.
812function point_plane_distance(plane, point) =
813 assert( _valid_plane(plane), "Invalid input plane." )
814 assert( is_vector(point,3), "The point should be a 3D point." )
815 let( plane = _normalize_plane(plane) )
816 point3d(plane)* point - plane[3];
817
818
819
820// the maximum distance from points to the plane
821function _pointlist_greatest_distance(points,plane) =
822 let(
823 normal = [plane[0],plane[1],plane[2]],
824 pt_nrm = points*normal
825 )
826 max( max(pt_nrm) - plane[3], -min(pt_nrm) + plane[3]) / norm(normal);
827
828
829// Function: are_points_on_plane()
830// Usage:
831// bool = are_points_on_plane(points, plane, [eps]);
832// Topics: Geometry, Planes, Points
833// Description:
834// Returns true if the given 3D points are on the given plane.
835// Arguments:
836// plane = The plane to test the points on.
837// points = The list of 3D points to test.
838// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
839function are_points_on_plane(points, plane, eps=EPSILON) =
840 assert( _valid_plane(plane), "Invalid plane." )
841 assert( is_matrix(points,undef,3) && len(points)>0, "Invalid pointlist." ) // using is_matrix it accepts len(points)==1
842 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
843 _pointlist_greatest_distance(points,plane) < eps;
844
845
846/// Internal Function: is_point_above_plane()
847/// Usage:
848/// bool = _is_point_above_plane(plane, point);
849/// Topics: Geometry, Planes
850// Description:
851/// Given a plane as [A,B,C,D] where the cartesian equation for that plane
852/// is Ax+By+Cz=D, determines if the given 3D point is on the side of that
853/// plane that the normal points towards. The normal of the plane is the
854/// same as [A,B,C].
855/// Arguments:
856/// plane = The [A,B,C,D] coefficients for the first plane equation `Ax+By+Cz=D`.
857/// point = The 3D point to test.
858function _is_point_above_plane(plane, point) =
859 point_plane_distance(plane, point) > EPSILON;
860
861
862
863// Section: Circle Calculations
864
865// Function: circle_line_intersection()
866// Usage:
867// pts = circle_line_intersection(r|d=, cp, line, [bounded], [eps=]);
868// Topics: Geometry, Circles, Lines, Intersection
869// Description:
870// Find intersection points between a 2D circle and a line, ray or segment specified by two points.
871// By default the line is unbounded. Returns the list of zero or more intersection points.
872// Arguments:
873// r = Radius of circle
874// cp = Center of circle
875// line = Two points defining the line
876// bounded = False for unbounded line, true for a segment, or a vector [false,true] or [true,false] to specify a ray with the first or second end unbounded. Default: false
877// ---
878// d = Diameter of circle
879// eps = Epsilon used for identifying the case with one solution. Default: `1e-9`
880// Example(2D): Standard intersection returns two points.
881// line = [[-15,2], [15,7]];
882// cp = [1,2]; r = 10;
883// translate(cp) circle(r=r);
884// color("black") stroke(line, endcaps="arrow2", width=0.5);
885// isects = circle_line_intersection(r=r, cp=cp, line=line);
886// color("red") move_copies(isects) circle(d=1);
887// Example(2D): Tangent intersection returns one point.
888// line = [[-10,12], [10,12]];
889// cp = [1,2]; r = 10;
890// translate(cp) circle(r=r);
891// color("black") stroke(line, endcaps="arrow2", width=0.5);
892// isects = circle_line_intersection(r=r, cp=cp, line=line);
893// color("#f44") move_copies(isects) circle(d=1);
894// Example(2D): A bounded ray might only intersect in one direction.
895// line = [[-5,2], [5,7]];
896// extended = [line[0], line[0]+22*unit(line[1]-line[0])];
897// cp = [1,2]; r = 10;
898// translate(cp) circle(r=r);
899// color("gray") dashed_stroke(extended, width=0.2);
900// color("black") stroke(line, endcap2="arrow2", width=0.5);
901// isects = circle_line_intersection(r=r, cp=cp, line=line, bounded=[true,false]);
902// color("#f44") move_copies(isects) circle(d=1);
903// Example(2D): If they don't intersect at all, then an empty list is returned.
904// line = [[-12,12], [12,8]];
905// cp = [-5,-2]; r = 10;
906// translate(cp) circle(r=r);
907// color("black") stroke(line, endcaps="arrow2", width=0.5);
908// isects = circle_line_intersection(r=r, cp=cp, line=line);
909// color("#f44") move_copies(isects) circle(d=1);
910function circle_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
911 assert(_valid_line(line,2), "Invalid 2d line.")
912 assert(is_vector(cp,2), "Circle center must be a 2-vector")
913 _circle_or_sphere_line_intersection(r, cp, line, bounded, d, eps);
914
915
916
917function _circle_or_sphere_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
918 let(r=get_radius(r=r,d=d,dflt=undef))
919 assert(is_num(r) && r>0, "Radius must be positive")
920 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition")
921 let(
922 bounded = force_list(bounded,2),
923 closest = line_closest_point(line,cp),
924 d = norm(closest-cp)
925 )
926 d > r ? [] :
927 let(
928 isect = approx(d,r,eps) ? [closest] :
929 let( offset = sqrt(r*r-d*d),
930 uvec=unit(line[1]-line[0])
931 ) [closest-offset*uvec, closest+offset*uvec]
932 )
933 [for(p=isect)
934 if ((!bounded[0] || (p-line[0])*(line[1]-line[0])>=0)
935 && (!bounded[1] || (p-line[1])*(line[0]-line[1])>=0)) p];
936
937
938// Function: circle_circle_intersection()
939// Usage:
940// pts = circle_circle_intersection(r1|d1=, cp1, r2|d2=, cp2, [eps]);
941// Topics: Geometry, Circles
942// Description:
943// Compute the intersection points of two circles. Returns a list of the intersection points, which
944// will contain two points in the general case, one point for tangent circles, or will be empty
945// if the circles do not intersect.
946// Arguments:
947// r1 = Radius of the first circle.
948// cp1 = Centerpoint of the first circle.
949// r2 = Radius of the second circle.
950// cp2 = Centerpoint of the second circle.
951// eps = Tolerance for detecting tangent circles. Default: EPSILON
952// ---
953// d1 = Diameter of the first circle.
954// d2 = Diameter of the second circle.
955// Example(2D,NoAxes): Circles intersect in two points.
956// $fn=32;
957// cp1 = [4,4]; r1 = 3;
958// cp2 = [7,7]; r2 = 2;
959// pts = circle_circle_intersection(r1, cp1, r2, cp2);
960// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
961// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
962// color("red") move_copies(pts) circle(r=.3);
963// Example(2D,NoAxes): Circles are tangent, so one intersection point:
964// $fn=32;
965// cp1 = [4,4]; r1 = 4;
966// cp2 = [4,10]; r2 = 2;
967// pts = circle_circle_intersection(r1, cp1, r2, cp2);
968// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
969// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
970// color("red") move_copies(pts) circle(r=.3);
971// Example(2D,NoAxes): Another tangent example:
972// $fn=32;
973// cp1 = [4,4]; r1 = 4;
974// cp2 = [5,5]; r2 = 4-sqrt(2);
975// pts = circle_circle_intersection(r1, cp1, r2, cp2);
976// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
977// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
978// color("red") move_copies(pts) circle(r=.3);
979// Example(2D,NoAxes): Circles do not intersect. Returns empty list.
980// $fn=32;
981// cp1 = [3,4]; r1 = 2;
982// cp2 = [7,10]; r2 = 3;
983// pts = circle_circle_intersection(r1, cp1, r2, cp2);
984// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
985// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
986// color("red") move_copies(pts) circle(r=.3);
987function circle_circle_intersection(r1, cp1, r2, cp2, eps=EPSILON, d1, d2) =
988 assert( is_path([cp1,cp2],dim=2), "Invalid center point(s)." )
989 let(
990 r1 = get_radius(r1=r1,d1=d1),
991 r2 = get_radius(r1=r2,d1=d2),
992 d = norm(cp2-cp1),
993 a = (cp2-cp1)/d,
994 b = [-a.y,a.x],
995 L = (r1^2-r2^2+d^2)/2/d,
996 hsqr = r1^2-L^2
997 )
998 approx(hsqr,0,eps) ? [L*a+cp1]
999 : hsqr<0 ? []
1000 : let(h=sqrt(hsqr))
1001 [L*a+h*b+cp1, L*a-h*b+cp1];
1002
1003
1004// Function: circle_2tangents()
1005// Usage:
1006// circ = circle_2tangents(r|d=, pt1, pt2, pt3, [tangents=]);
1007// circ = circle_2tangents(r|d=, [PT1, PT2, PT3], [tangents=]);
1008// Topics: Geometry, Circles, Tangents
1009// Description:
1010// Given a pair of rays with a common origin, and a known circle radius/diameter, finds
1011// the centerpoint for the circle of that size that touches both rays tangentally.
1012// Both rays start at `pt2`, one passing through `pt1`, and the other through `pt3`.
1013// .
1014// When called with collinear rays, returns `undef`.
1015// Otherwise, when called with `tangents=false`, returns `[CP,NORMAL]`.
1016// Otherwise, when called with `tangents=true`, returns `[CP,NORMAL,TANPT1,TANPT2]`.
1017// - CP is the centerpoint of the circle.
1018// - NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
1019// - TANPT1 is the point where the circle is tangent to the ray `[pt2,pt1]`.
1020// - TANPT2 is the point where the circle is tangent to the ray `[pt2,pt3]`.
1021// Figure(3D,Med,NoAxes,VPD=130,VPT=[29,19,3],VPR=[55,0,25]):
1022// pts = [[45,10,-5], [10,5,10], [15,40,5]];
1023// rad = 15;
1024// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2], tangents=true);
1025// cp = circ[0]; n = circ[1]; tp1 = circ[2]; tp2 = circ[3];
1026// color("yellow") stroke(pts, endcaps="arrow2");
1027// color("purple") move_copies([cp,tp1,tp2]) sphere(d=2, $fn=12);
1028// color("lightgray") stroke([cp,tp2], width=0.5);
1029// stroke([cp,cp+n*20], endcap2="arrow2");
1030// labels = [
1031// ["pt1", "blue", 2.5, [ 4, 0, 1], pts[0]],
1032// ["pt2", "blue", 2.5, [-4, 0,-3], pts[1]],
1033// ["pt3", "blue", 2.5, [ 4, 0, 1], pts[2]],
1034// ["r", "blue", 2.5, [ 0,-2, 2], (cp+tp2)/2],
1035// ["CP", "brown", 2.5, [ 6,-4, 3], cp],
1036// ["Normal", "brown", 2.0, [ 5, 2, 1], cp+20*n],
1037// ["TanPt1", "brown", 2.0, [-5,-4, 0], tp1],
1038// ["TanPt2", "brown", 2.0, [-5, 0, 2], tp2],
1039// ];
1040// for(l=labels)
1041// color(l[1]) move(l[4]+l[3]) rot([55,0,25])
1042// linear_extrude(height=0.1)
1043// text(text=l[0], size=l[2], halign="center", valign="center");
1044// color("green",0.5) move(cp) cyl(h=0.1, r=rad, orient=n, $fn=36);
1045// Arguments:
1046// r = The radius of the circle to find.
1047// pt1 = A point that the first ray passes though.
1048// pt2 = The starting point of both rays.
1049// pt3 = A point that the second ray passes though.
1050// ---
1051// d = The diameter of the circle to find.
1052// tangents = If true, extended information about the tangent points is calculated and returned. Default: false
1053// Example(2D):
1054// pts = [[40,40], [10,10], [55,5]]; rad = 10;
1055// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2]);
1056// stroke(pts, endcaps="arrow2");
1057// color("red") move(circ[0]) circle(r=rad);
1058// Example(2D):
1059// pts = [[20,40], [10,10], [55,20]]; rad = 10;
1060// circ = circle_2tangents(r=rad, pt1=pts[0], pt2=pts[1], pt3=pts[2], tangents=true);
1061// stroke(pts, endcaps="arrow2");
1062// color("red") move(circ[0]) circle(r=rad);
1063// color("blue") move_copies(select(circ,2,3)) circle(d=2);
1064// Example(3D): Fit into 3D path corner.
1065// pts = [[45,5,10], [10,10,15], [30,40,30]]; rad = 10;
1066// circ = circle_2tangents(rad, [pts[0], pts[1], pts[2]]);
1067// stroke(pts, endcaps="arrow2");
1068// color("red") move(circ[0]) cyl(h=10, r=rad, orient=circ[1]);
1069// Example(3D):
1070// path = yrot(20, p=path3d(star(d=100, n=5, step=2)));
1071// stroke(path, closed=true);
1072// for (i = [0:1:5]) {
1073// crn = select(path, i*2-1, i*2+1);
1074// ci = circle_2tangents(5, crn[0], crn[1], crn[2]);
1075// move(ci[0]) cyl(h=10,r=5,orient=ci[1]);
1076// }
1077function circle_2tangents(r, pt1, pt2, pt3, tangents=false, d) =
1078 let(r = get_radius(r=r, d=d, dflt=undef))
1079 assert(r!=undef, "Must specify either r or d.")
1080 assert( ( is_path(pt1) && len(pt1)==3 && is_undef(pt2) && is_undef(pt3))
1081 || (is_matrix([pt1,pt2,pt3]) && (len(pt1)==2 || len(pt1)==3) ),
1082 "Invalid input points." )
1083 is_undef(pt2)
1084 ? circle_2tangents(r, pt1[0], pt1[1], pt1[2], tangents=tangents)
1085 : is_collinear(pt1, pt2, pt3)? undef :
1086 let(
1087 v1 = unit(pt1 - pt2),
1088 v2 = unit(pt3 - pt2),
1089 vmid = unit(mean([v1, v2])),
1090 n = vector_axis(v1, v2),
1091 a = vector_angle(v1, v2),
1092 hyp = r / sin(a/2),
1093 cp = pt2 + hyp * vmid
1094 )
1095 !tangents ? [cp, n] :
1096 let(
1097 x = hyp * cos(a/2),
1098 tp1 = pt2 + x * v1,
1099 tp2 = pt2 + x * v2
1100 )
1101 [cp, n, tp1, tp2];
1102
1103
1104// Function: circle_3points()
1105// Usage:
1106// circ = circle_3points(pt1, pt2, pt3);
1107// circ = circle_3points([PT1, PT2, PT3]);
1108// Topics: Geometry, Circles
1109// Description:
1110// Returns the [CENTERPOINT, RADIUS, NORMAL] of the circle that passes through three non-collinear
1111// points where NORMAL is the normal vector of the plane that the circle is on (UP or DOWN if the points are 2D).
1112// The centerpoint will be a 2D or 3D vector, depending on the points input. If all three
1113// points are 2D, then the resulting centerpoint will be 2D, and the normal will be UP ([0,0,1]).
1114// If any of the points are 3D, then the resulting centerpoint will be 3D. If the three points are
1115// collinear, then `[undef,undef,undef]` will be returned. The normal will be a normalized 3D
1116// vector with a non-negative Z axis. Instead of 3 arguments, it is acceptable to input the 3 points
1117// as a list given in `pt1`, leaving `pt2`and `pt3` as undef.
1118// Arguments:
1119// pt1 = The first point.
1120// pt2 = The second point.
1121// pt3 = The third point.
1122// Example(2D):
1123// pts = [[60,40], [10,10], [65,5]];
1124// circ = circle_3points(pts[0], pts[1], pts[2]);
1125// translate(circ[0]) color("green") stroke(circle(r=circ[1]),closed=true,$fn=72);
1126// translate(circ[0]) color("red") circle(d=3, $fn=12);
1127// move_copies(pts) color("blue") circle(d=3, $fn=12);
1128function circle_3points(pt1, pt2, pt3) =
1129 (is_undef(pt2) && is_undef(pt3) && is_list(pt1))
1130 ? circle_3points(pt1[0], pt1[1], pt1[2])
1131 : assert( is_vector(pt1) && is_vector(pt2) && is_vector(pt3)
1132 && max(len(pt1),len(pt2),len(pt3))<=3 && min(len(pt1),len(pt2),len(pt3))>=2,
1133 "Invalid point(s)." )
1134 is_collinear(pt1,pt2,pt3)? [undef,undef,undef] :
1135 let(
1136 v = [ point3d(pt1), point3d(pt2), point3d(pt3) ], // triangle vertices
1137 ed = [for(i=[0:2]) v[(i+1)%3]-v[i] ], // triangle edge vectors
1138 pm = [for(i=[0:2]) v[(i+1)%3]+v[i] ]/2, // edge mean points
1139 es = sortidx( [for(di=ed) norm(di) ] ),
1140 e1 = ed[es[1]], // take the 2 longest edges
1141 e2 = ed[es[2]],
1142 n0 = vector_axis(e1,e2), // normal standardization
1143 n = n0.z<0? -n0 : n0,
1144 sc = plane_intersection(
1145 [ each e1, e1*pm[es[1]] ], // planes orthogonal to 2 edges
1146 [ each e2, e2*pm[es[2]] ],
1147 [ each n, n*v[0] ]
1148 ), // triangle plane
1149 cp = len(pt1)+len(pt2)+len(pt3)>6 ? sc : [sc.x, sc.y],
1150 r = norm(sc-v[0])
1151 ) [ cp, r, n ];
1152
1153
1154
1155// Function: circle_point_tangents()
1156// Usage:
1157// tangents = circle_point_tangents(r|d=, cp, pt);
1158// Topics: Geometry, Circles, Tangents
1159// Description:
1160// Given a 2d circle and a 2d point outside that circle, finds the 2d tangent point(s) on the circle for a
1161// line passing through the point. Returns a list of zero or more 2D tangent points.
1162// Arguments:
1163// r = Radius of the circle.
1164// cp = The coordinates of the 2d circle centerpoint.
1165// pt = The coordinates of the 2d external point.
1166// ---
1167// d = Diameter of the circle.
1168// Example(2D):
1169// cp = [-10,-10]; r = 30; pt = [30,10];
1170// tanpts = circle_point_tangents(r=r, cp=cp, pt=pt);
1171// color("yellow") translate(cp) circle(r=r);
1172// color("cyan") for(tp=tanpts) {stroke([tp,pt]); stroke([tp,cp]);}
1173// color("red") move_copies(tanpts) circle(d=3,$fn=12);
1174// color("blue") move_copies([cp,pt]) circle(d=3,$fn=12);
1175function circle_point_tangents(r, cp, pt, d) =
1176 assert(is_finite(r) || is_finite(d), "Invalid radius or diameter." )
1177 assert(is_path([cp, pt],dim=2), "Invalid center point or external point.")
1178 let(
1179 r = get_radius(r=r, d=d, dflt=1),
1180 delta = pt - cp,
1181 dist = norm(delta),
1182 baseang = atan2(delta.y,delta.x)
1183 ) dist < r? [] :
1184 approx(dist,r)? [pt] :
1185 let(
1186 relang = acos(r/dist),
1187 angs = [baseang + relang, baseang - relang]
1188 ) [for (ang=angs) cp + r*[cos(ang),sin(ang)]];
1189
1190
1191// Function: circle_circle_tangents()
1192// Usage:
1193// segs = circle_circle_tangents(r1|d1=, cp1, r2|d2=, cp2);
1194// Topics: Geometry, Circles, Tangents
1195// Description:
1196// Computes 2d lines tangents to a pair of circles in 2d. Returns a list of line endpoints [p1,p2] where
1197// p2 is the tangent point on circle 1 and p2 is the tangent point on circle 2.
1198// If four tangents exist then the first one the left hand exterior tangent as regarded looking from
1199// circle 1 toward circle 2. The second value is the right hand exterior tangent. The third entry
1200// gives the interior tangent that starts on the left of circle 1 and crosses to the right side of
1201// circle 2. And the fourth entry is the last interior tangent that starts on the right side of
1202// circle 1. If the circles intersect then the interior tangents don't exist and the function
1203// returns only two entries. If one circle is inside the other one then no tangents exist
1204// so the function returns the empty set. When the circles are tangent a degenerate tangent line
1205// passes through the point of tangency of the two circles: this degenerate line is NOT returned.
1206// Arguments:
1207// r1 = Radius of the first circle.
1208// cp1 = Centerpoint of the first circle.
1209// r2 = Radius of the second circle.
1210// cp2 = Centerpoint of the second circle.
1211// ---
1212// d1 = Diameter of the first circle.
1213// d2 = Diameter of the second circle.
1214// Example(2D,NoAxes): Four tangents, first in green, second in black, third in blue, last in red.
1215// $fn=32;
1216// cp1 = [3,4]; r1 = 2;
1217// cp2 = [7,10]; r2 = 3;
1218// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1219// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1220// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1221// colors = ["green","black","blue","red"];
1222// for(i=[0:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1223// Example(2D,NoAxes): Circles overlap so only exterior tangents exist.
1224// $fn=32;
1225// cp1 = [4,4]; r1 = 3;
1226// cp2 = [7,7]; r2 = 2;
1227// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1228// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1229// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1230// colors = ["green","black","blue","red"];
1231// for(i=[0:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1232// Example(2D,NoAxes): Circles are tangent. Only exterior tangents are returned. The degenerate internal tangent is not returned.
1233// $fn=32;
1234// cp1 = [4,4]; r1 = 4;
1235// cp2 = [4,10]; r2 = 2;
1236// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1237// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1238// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1239// colors = ["green","black","blue","red"];
1240// for(i=[0:1:len(pts)-1]) color(colors[i]) stroke(pts[i],width=0.2);
1241// Example(2D,NoAxes): One circle is inside the other: no tangents exist. If the interior circle is tangent the single degenerate tangent will not be returned.
1242// $fn=32;
1243// cp1 = [4,4]; r1 = 4;
1244// cp2 = [5,5]; r2 = 2;
1245// pts = circle_circle_tangents(r1, cp1, r2, cp2);
1246// move(cp1) stroke(circle(r=r1), width=0.2, closed=true);
1247// move(cp2) stroke(circle(r=r2), width=0.2, closed=true);
1248// echo(pts); // Returns []
1249function circle_circle_tangents(r1, cp1, r2, cp2, d1, d2) =
1250 assert( is_path([cp1,cp2],dim=2), "Invalid center point(s)." )
1251 let(
1252 r1 = get_radius(r1=r1,d1=d1),
1253 r2 = get_radius(r1=r2,d1=d2),
1254 Rvals = [r2-r1, r2-r1, -r2-r1, -r2-r1]/norm(cp1-cp2),
1255 kvals = [-1,1,-1,1],
1256 ext = [1,1,-1,-1],
1257 N = 1-sqr(Rvals[2])>=0 ? 4 :
1258 1-sqr(Rvals[0])>=0 ? 2 : 0,
1259 coef= [
1260 for(i=[0:1:N-1]) [
1261 [Rvals[i], -kvals[i]*sqrt(1-sqr(Rvals[i]))],
1262 [kvals[i]*sqrt(1-sqr(Rvals[i])), Rvals[i]]
1263 ] * unit(cp2-cp1)
1264 ]
1265 ) [
1266 for(i=[0:1:N-1]) let(
1267 pt = [
1268 cp1-r1*coef[i],
1269 cp2-ext[i]*r2*coef[i]
1270 ]
1271 ) if (pt[0]!=pt[1]) pt
1272 ];
1273
1274
1275
1276/// Internal Function: _noncollinear_triple()
1277/// Usage:
1278/// bool = _noncollinear_triple(points);
1279/// Topics: Geometry, Noncollinearity
1280/// Description:
1281/// Finds the indices of three non-collinear points from the pointlist `points`.
1282/// It selects two well separated points to define a line and chooses the third point
1283/// to be the point farthest off the line. The points do not necessarily having the
1284/// same winding direction as the polygon so they cannot be used to determine the
1285/// winding direction or the direction of the normal.
1286/// If all points are collinear returns [] when `error=true` or an error otherwise .
1287/// Arguments:
1288/// points = List of input points.
1289/// error = Defines the behaviour for collinear input points. When `true`, produces an error, otherwise returns []. Default: `true`.
1290/// eps = Tolerance for collinearity test. Default: EPSILON.
1291function _noncollinear_triple(points,error=true,eps=EPSILON) =
1292 assert( is_path(points), "Invalid input points." )
1293 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
1294 len(points)<3 ? [] :
1295 let(
1296 pa = points[0],
1297 b = furthest_point(pa, points),
1298 pb = points[b],
1299 nrm = norm(pa-pb)
1300 )
1301 nrm <= eps ?
1302 assert(!error, "Cannot find three noncollinear points in pointlist.") [] :
1303 let(
1304 n = (pb-pa)/nrm,
1305 distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]
1306 )
1307 max(distlist) < eps*nrm ?
1308 assert(!error, "Cannot find three noncollinear points in pointlist.") [] :
1309 [0, b, max_index(distlist)];
1310
1311
1312
1313// Section: Sphere Calculations
1314
1315
1316// Function: sphere_line_intersection()
1317// Usage:
1318// isect = sphere_line_intersection(r|d=, cp, line, [bounded], [eps=]);
1319// Topics: Geometry, Spheres, Lines, Intersection
1320// Description:
1321// Find intersection points between a sphere and a line, ray or segment specified by two points.
1322// By default the line is unbounded.
1323// Arguments:
1324// r = Radius of sphere
1325// cp = Centerpoint of sphere
1326// line = Two points defining the line
1327// bounded = false for unbounded line, true for a segment, or a vector [false,true] or [true,false] to specify a ray with the first or second end unbounded. Default: false
1328// ---
1329// d = diameter of sphere
1330// eps = epsilon used for identifying the case with one solution. Default: 1e-9
1331// Example(3D):
1332// cp = [10,20,5]; r = 40;
1333// line = [[-50,-10,25], [70,0,40]];
1334// isects = sphere_line_intersection(r=r, cp=cp, line=line);
1335// color("cyan") stroke(line);
1336// move(cp) sphere(r=r, $fn=72);
1337// color("red") move_copies(isects) sphere(d=3, $fn=12);
1338function sphere_line_intersection(r, cp, line, bounded=false, d, eps=EPSILON) =
1339 assert(_valid_line(line,3), "Invalid 3d line.")
1340 assert(is_vector(cp,3), "Sphere center must be a 3-vector")
1341 _circle_or_sphere_line_intersection(r, cp, line, bounded, d, eps);
1342
1343
1344
1345
1346// Section: Polygons
1347
1348// Function: polygon_area()
1349// Usage:
1350// area = polygon_area(poly, [signed]);
1351// Topics: Geometry, Polygons, Area
1352// Description:
1353// Given a 2D or 3D simple planar polygon, returns the area of that polygon.
1354// If the polygon is non-planar the result is `undef.` If the polygon is self-intersecting
1355// then the return will be a meaningless number.
1356// When `signed` is true and the polygon is 2d, a signed area is returned: a positive area indicates a counter-clockwise polygon.
1357// The area of 3d polygons is always nonnegative.
1358// Arguments:
1359// poly = Polygon to compute the area of.
1360// signed = If true, a signed area is returned. Default: false.
1361function polygon_area(poly, signed=false) =
1362 assert(is_path(poly), "Invalid polygon." )
1363 len(poly)<3 ? 0 :
1364 len(poly)==3 ?
1365 let( total= len(poly[0])==2 ? 0.5*cross(poly[2]-poly[0],poly[2]-poly[1]) : 0.5*norm(cross(poly[2]-poly[0],poly[2]-poly[1])))
1366 signed ? total : abs(total) :
1367 len(poly[0])==2
1368 ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 )
1369 signed ? total : abs(total)
1370 : let( plane = plane_from_polygon(poly) )
1371 is_undef(plane) ? undef :
1372 let(
1373 n = plane_normal(plane),
1374 total =
1375 -sum([ for(i=[1:1:len(poly)-2])
1376 cross(poly[i]-poly[0], poly[i+1]-poly[0])
1377 ]) * n/2
1378 )
1379 signed ? total : abs(total);
1380
1381
1382// Function: centroid()
1383// Usage:
1384// c = centroid(object, [eps]);
1385// Topics: Geometry, Polygons, Centroid
1386// Description:
1387// Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
1388// Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
1389// If you provide a non-planar or collinear polygon you will get an error. For self-intersecting
1390// polygons you may get an error or you may get meaningless results.
1391// .
1392// Given a [[region|regions.scad]], returns the 2D coordinates of the region's centroid.
1393// .
1394// Given a manifold [[VNF|vnf.scad]] then returns the 3D centroid of the polyhedron. The VNF must
1395// describe a valid polyhedron with consistent face direction and no holes in the mesh; otherwise
1396// the results are undefined.
1397// Arguments:
1398// object = object to compute the centroid of
1399// eps = epsilon value for identifying degenerate cases
1400// Example(2D):
1401// path = [
1402// [-10,10], [-5,15], [15,15], [20,0],
1403// [15,-5], [25,-20], [25,-27], [15,-20],
1404// [0,-30], [-15,-25], [-5,-5]
1405// ];
1406// linear_extrude(height=0.01) polygon(path);
1407// cp = centroid(path);
1408// color("red") move(cp) sphere(d=2);
1409function centroid(object,eps=EPSILON) =
1410 assert(is_finite(eps) && (eps>=0), "The tolerance should a non-negative value." )
1411 is_vnf(object) ? _vnf_centroid(object,eps)
1412 : is_path(object,[2,3]) ? _polygon_centroid(object,eps)
1413 : is_region(object) ? (len(object)==1 ? _polygon_centroid(object[0],eps) : _region_centroid(object,eps))
1414 : assert(false, "Input must be a VNF, a region, or a 2D or 3D polygon");
1415
1416
1417/// Internal Function: _region_centroid()
1418/// Compute centroid of region
1419function _region_centroid(region,eps=EPSILON) =
1420 let(
1421 region=force_region(region),
1422 parts = region_parts(region),
1423 // Rely on region_parts returning all outside polygons clockwise
1424 // and inside (hole) polygons counterclockwise, so areas have reversed sign
1425 cent_area = [for(R=parts, p=R)
1426 let(A=polygon_area(p,signed=true))
1427 [A*_polygon_centroid(p),A]],
1428 total = sum(cent_area)
1429 )
1430 total[0]/total[1];
1431
1432
1433/// Internal Function: _polygon_centroid()
1434/// Usage:
1435/// cpt = _polygon_centroid(poly);
1436/// Topics: Geometry, Polygons, Centroid
1437/// Description:
1438/// Given a simple 2D polygon, returns the 2D coordinates of the polygon's centroid.
1439/// Given a simple 3D planar polygon, returns the 3D coordinates of the polygon's centroid.
1440/// Collinear points produce an error. The results are meaningless for self-intersecting
1441/// polygons or an error is produced.
1442/// Arguments:
1443/// poly = Points of the polygon from which the centroid is calculated.
1444/// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1445function _polygon_centroid(poly, eps=EPSILON) =
1446 assert( is_path(poly,dim=[2,3]), "The input must be a 2D or 3D polygon." )
1447 let(
1448 n = len(poly[0])==2 ? 1 :
1449 let( plane = plane_from_points(poly, fast=false))
1450 assert(!is_undef(plane), "The polygon must be planar." )
1451 plane_normal(plane),
1452 v0 = poly[0] ,
1453 val = sum([
1454 for(i=[1:len(poly)-2])
1455 let(
1456 v1 = poly[i],
1457 v2 = poly[i+1],
1458 area = cross(v2-v0,v1-v0)*n
1459 ) [ area, (v0+v1+v2)*area ]
1460 ])
1461 )
1462 assert(!approx(val[0],0, eps), "The polygon is self-intersecting or its points are collinear.")
1463 val[1]/val[0]/3;
1464
1465
1466
1467// Function: polygon_normal()
1468// Usage:
1469// vec = polygon_normal(poly);
1470// Topics: Geometry, Polygons
1471// Description:
1472// Given a 3D simple planar polygon, returns a unit normal vector for the polygon. The vector
1473// is oriented so that if the normal points towards the viewer, the polygon winds in the clockwise
1474// direction. If the polygon has zero area, returns `undef`. If the polygon is self-intersecting
1475// the the result is undefined. It doesn't check for coplanarity.
1476// Arguments:
1477// poly = The list of 3D path points for the perimeter of the polygon.
1478// Example(3D):
1479// path = rot([0,30,15], p=path3d(star(n=5, d=100, step=2)));
1480// stroke(path, closed=true);
1481// n = polygon_normal(path);
1482// rot(from=UP, to=n)
1483// color("red")
1484// stroke([[0,0,0], [0,0,20]], endcap2="arrow2");
1485function polygon_normal(poly) =
1486 assert(is_path(poly,dim=3), "Invalid 3D polygon." )
1487 let(
1488 area_vec = sum([for(i=[1:len(poly)-2])
1489 cross(poly[i]-poly[0],
1490 poly[i+1]-poly[i])])
1491 )
1492 unit(-area_vec, error=undef);
1493
1494
1495// Function: point_in_polygon()
1496// Usage:
1497// bool = point_in_polygon(point, poly, [nonzero], [eps])
1498// Topics: Geometry, Polygons
1499// Description:
1500// This function tests whether the given 2D point is inside, outside or on the boundary of
1501// the specified 2D polygon.
1502// The polygon is given as a list of 2D points, not including the repeated end point.
1503// Returns -1 if the point is outside the polygon.
1504// Returns 0 if the point is on the boundary.
1505// Returns 1 if the point lies in the interior.
1506// The polygon does not need to be simple: it may have self-intersections.
1507// But the polygon cannot have holes (it must be simply connected).
1508// Rounding errors may give mixed results for points on or near the boundary.
1509// .
1510// When polygons intersect themselves different definitions exist for determining which points
1511// are inside the polygon. The figure below shows the difference.
1512// OpenSCAD uses the Even-Odd rule when creating polygons, where membership in overlapping regions
1513// depends on how many times they overlap. The Nonzero rule considers point inside the polygon if
1514// the polygon overlaps them any number of times. For more information see
1515// https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Even–odd_rule.
1516// Figure(2D,Med,NoAxes):
1517// a=20;
1518// b=30;
1519// ofs = 17;
1520// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)-20]];
1521// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1522// xflip(offset(curve,r=ofs)),
1523// xflip(reverse(curve)),
1524// curve
1525// ));
1526// left(40){
1527// polygon(path);
1528// color("red")stroke(path, width=1, closed=true);
1529// color("red")back(28/(2/3))text("Even-Odd", size=5/(2/3), halign="center");
1530// }
1531// right(40){
1532// dp = polygon_parts(path,nonzero=true);
1533// region(dp);
1534// color("red"){stroke(path,width=1,closed=true);
1535// back(28/(2/3))text("Nonzero", size=5/(2/3), halign="center");
1536// }
1537// }
1538// Arguments:
1539// point = The 2D point to check
1540// poly = The list of 2D points forming the perimeter of the polygon.
1541// nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd". Default: false (Even-Odd)
1542// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1543// Example(2D): With nonzero set to false (the default), we get this result. Green dots are inside the polygon and red are outside:
1544// a=20*2/3;
1545// b=30*2/3;
1546// ofs = 17*2/3;
1547// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)]];
1548// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1549// xflip(offset(curve,r=ofs)),
1550// xflip(reverse(curve)),
1551// curve
1552// ));
1553// stroke(path,closed=true);
1554// pts = [[0,0],[10,0],[0,20]];
1555// for(p=pts){
1556// color(point_in_polygon(p,path)==1 ? "green" : "red")
1557// move(p)circle(r=1.5, $fn=12);
1558// }
1559// Example(2D): With nonzero set to true, one dot changes color:
1560// a=20*2/3;
1561// b=30*2/3;
1562// ofs = 17*2/3;
1563// curve = [for(theta=[0:10:140]) [a * theta/360*2*PI - b*sin(theta), a-b*cos(theta)]];
1564// path = deduplicate(concat( reverse(offset(curve,r=ofs)),
1565// xflip(offset(curve,r=ofs)),
1566// xflip(reverse(curve)),
1567// curve
1568// ));
1569// stroke(path,closed=true);
1570// pts = [[0,0],[10,0],[0,20]];
1571// for(p=pts){
1572// color(point_in_polygon(p,path,nonzero=true)==1 ? "green" : "red")
1573// move(p)circle(r=1.5, $fn=12);
1574// }
1575
1576// Internal function for point_in_polygon
1577
1578function _point_above_below_segment(point, edge) =
1579 let( edge = edge - [point, point] )
1580 edge[0].y <= 0
1581 ? (edge[1].y > 0 && cross(edge[0], edge[1]-edge[0]) > 0) ? 1 : 0
1582 : (edge[1].y <= 0 && cross(edge[0], edge[1]-edge[0]) < 0) ? -1 : 0;
1583
1584
1585function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) =
1586 // Original algorithms from http://geomalgorithms.com/a03-_inclusion.html
1587 assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
1588 "The point and polygon should be in 2D. The polygon should have more that 2 points." )
1589 assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
1590 // Check bounding box
1591 let(
1592 box = pointlist_bounds(poly)
1593 )
1594 point.x<box[0].x-eps || point.x>box[1].x+eps
1595 || point.y<box[0].y-eps || point.y>box[1].y+eps ? -1
1596 :
1597 // Does the point lie on any edges? If so return 0.
1598 let(
1599 segs = pair(poly,true),
1600 on_border = [for (seg=segs)
1601 if (norm(seg[0]-seg[1])>eps && _is_point_on_line(point, seg, SEGMENT, eps=eps)) 1]
1602 )
1603 on_border != [] ? 0 :
1604 nonzero // Compute winding number and return 1 for interior, -1 for exterior
1605 ? let(
1606 winding = [
1607 for(seg=segs)
1608 let(
1609 p0=seg[0]-point,
1610 p1=seg[1]-point
1611 )
1612 if (norm(p0-p1)>eps)
1613 p0.y <=0
1614 ? p1.y > 0 && cross(p0,p1-p0)>0 ? 1 : 0
1615 : p1.y <=0 && cross(p0,p1-p0)<0 ? -1: 0
1616 ]
1617 )
1618 sum(winding) != 0 ? 1 : -1
1619 : // or compute the crossings with the ray [point, point+[1,0]]
1620 let(
1621 cross = [
1622 for(seg=segs)
1623 let(
1624 p0 = seg[0]-point,
1625 p1 = seg[1]-point
1626 )
1627 if (
1628 ( (p1.y>eps && p0.y<=eps) || (p1.y<=eps && p0.y>eps) )
1629 && -eps < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y)
1630 )
1631 1
1632 ]
1633 )
1634 2*(len(cross)%2)-1;
1635
1636
1637
1638// Function: polygon_line_intersection()
1639// Usage:
1640// pt = polygon_line_intersection(poly, line, [bounded], [nonzero], [eps]);
1641// Topics: Geometry, Polygons, Lines, Intersection
1642// Description:
1643// Takes a possibly bounded line, and a 2D or 3D planar polygon, and finds their intersection. Note the polygon is
1644// treated as its boundary and interior, so the intersection may include both points and line segments.
1645// If the line does not intersect the polygon returns `undef`.
1646// In 3D if the line is not on the plane of the polygon but intersects it then you get a single intersection point.
1647// Otherwise the polygon and line are in the same plane, or when your input is 2D, you will get a list of segments and
1648// single point lists. Use `is_vector` to distinguish these two cases.
1649// .
1650// In the 2D case, a common result is a list containing a single segment, which lists the two intersection points
1651// with the boundary of the polygon.
1652// When single points are in the intersection (the line just touches a polygon corner) they appear on the segment
1653// list as lists of a single point
1654// (like single point segments) so a single point intersection in 2D has the form `[[[x,y,z]]]` as compared
1655// to a single point intersection in 3D which has the form `[x,y,z]`. You can identify whether an entry in the
1656// segment list is a true segment by checking its length, which will be 2 for a segment and 1 for a point.
1657// Arguments:
1658// poly = The 3D planar polygon to find the intersection with.
1659// line = A list of two distinct 3D points on the line.
1660// bounded = If false, the line is considered unbounded. If true, it is treated as a bounded line segment. If given as `[true, false]` or `[false, true]`, the boundedness of the points are specified individually, allowing the line to be treated as a half-bounded ray. Default: false (unbounded)
1661// nonzero = set to true to use the nonzero rule for determining it points are in a polygon. See point_in_polygon. Default: false.
1662// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
1663// Example(3D): The line intersects the 3d hexagon in a single point.
1664// hex = zrot(140,p=rot([-45,40,20],p=path3d(hexagon(r=15))));
1665// line = [[5,0,-13],[-3,-5,13]];
1666// isect = polygon_line_intersection(hex,line);
1667// stroke(hex,closed=true);
1668// stroke(line);
1669// color("red")move(isect)sphere(r=1,$fn=12);
1670// Example(2D): In 2D things are more complicated. The output is a list of intersection parts, in the simplest case a single segment.
1671// hex = hexagon(r=15);
1672// line = [[-20,10],[25,-7]];
1673// isect = polygon_line_intersection(hex,line);
1674// stroke(hex,closed=true);
1675// stroke(line,endcaps="arrow2");
1676// color("red")
1677// for(part=isect)
1678// if(len(part)==1)
1679// move(part[0]) sphere(r=1);
1680// else
1681// stroke(part);
1682// Example(2D): Here the line is treated as a ray.
1683// hex = hexagon(r=15);
1684// line = [[0,0],[25,-7]];
1685// isect = polygon_line_intersection(hex,line,RAY);
1686// stroke(hex,closed=true);
1687// stroke(line,endcap2="arrow2");
1688// color("red")
1689// for(part=isect)
1690// if(len(part)==1)
1691// move(part[0]) circle(r=1,$fn=12);
1692// else
1693// stroke(part);
1694// Example(2D): Here the intersection is a single point, which is returned as a single point "path" on the path list.
1695// hex = hexagon(r=15);
1696// line = [[15,-10],[15,13]];
1697// isect = polygon_line_intersection(hex,line,RAY);
1698// stroke(hex,closed=true);
1699// stroke(line,endcap2="arrow2");
1700// color("red")
1701// for(part=isect)
1702// if(len(part)==1)
1703// move(part[0]) circle(r=1,$fn=12);
1704// else
1705// stroke(part);
1706// Example(2D): Another way to get a single segment
1707// hex = hexagon(r=15);
1708// line = rot(30,p=[[15,-10],[15,25]],cp=[15,0]);
1709// isect = polygon_line_intersection(hex,line,RAY);
1710// stroke(hex,closed=true);
1711// stroke(line,endcap2="arrow2");
1712// color("red")
1713// for(part=isect)
1714// if(len(part)==1)
1715// move(part[0]) circle(r=1,$fn=12);
1716// else
1717// stroke(part);
1718// Example(2D): Single segment again
1719// star = star(r=15,n=8,step=2);
1720// line = [[20,-5],[-5,20]];
1721// isect = polygon_line_intersection(star,line,RAY);
1722// stroke(star,closed=true);
1723// stroke(line,endcap2="arrow2");
1724// color("red")
1725// for(part=isect)
1726// if(len(part)==1)
1727// move(part[0]) circle(r=1,$fn=12);
1728// else
1729// stroke(part);
1730// Example(2D): Solution is two points
1731// star = star(r=15,n=8,step=3);
1732// line = rot(22.5,p=[[15,-10],[15,20]],cp=[15,0]);
1733// isect = polygon_line_intersection(star,line,SEGMENT);
1734// stroke(star,closed=true);
1735// stroke(line);
1736// color("red")
1737// for(part=isect)
1738// if(len(part)==1)
1739// move(part[0]) circle(r=1,$fn=12);
1740// else
1741// stroke(part);
1742// Example(2D): Solution is list of three segments
1743// star = star(r=25,ir=9,n=8);
1744// line = [[-25,12],[25,12]];
1745// isect = polygon_line_intersection(star,line);
1746// stroke(star,closed=true);
1747// stroke(line,endcaps="arrow2");
1748// color("red")
1749// for(part=isect)
1750// if(len(part)==1)
1751// move(part[0]) circle(r=1,$fn=12);
1752// else
1753// stroke(part);
1754// Example(2D): Solution is a mixture of segments and points
1755// star = star(r=25,ir=9,n=7);
1756// line = [left(10,p=star[8]), right(50,p=star[8])];
1757// isect = polygon_line_intersection(star,line);
1758// stroke(star,closed=true);
1759// stroke(line,endcaps="arrow2");
1760// color("red")
1761// for(part=isect)
1762// if(len(part)==1)
1763// move(part[0]) circle(r=1,$fn=12);
1764// else
1765// stroke(part);
1766function polygon_line_intersection(poly, line, bounded=false, nonzero=false, eps=EPSILON) =
1767 assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
1768 assert(is_path(poly,dim=[2,3]), "Invalid polygon." )
1769 assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
1770 assert(_valid_line(line,dim=len(poly[0]),eps=eps), "Line invalid or does not match polygon dimension." )
1771 let(
1772 bounded = force_list(bounded,2),
1773 poly = deduplicate(poly)
1774 )
1775 len(poly[0])==2 ? // planar case
1776 let(
1777 linevec = unit(line[1] - line[0]),
1778 bound = 100*max(v_abs(flatten(pointlist_bounds(poly)))),
1779 boundedline = [line[0] + (bounded[0]? 0 : -bound) * linevec,
1780 line[1] + (bounded[1]? 0 : bound) * linevec],
1781 parts = split_region_at_region_crossings(boundedline, [poly], closed1=false)[0][0],
1782 inside = [
1783 if(point_in_polygon(parts[0][0], poly, nonzero=nonzero, eps=eps) == 0)
1784 [parts[0][0]], // Add starting point if it is on the polygon
1785 for(part = parts)
1786 if (point_in_polygon(mean(part), poly, nonzero=nonzero, eps=eps) >=0 )
1787 part
1788 else if(len(part)==2 && point_in_polygon(part[1], poly, nonzero=nonzero, eps=eps) == 0)
1789 [part[1]] // Add segment end if it is on the polygon
1790 ]
1791 )
1792 (len(inside)==0 ? undef : _merge_segments(inside, [inside[0]], eps))
1793 : // 3d case
1794 let(indices = _noncollinear_triple(poly))
1795 indices==[] ? undef : // Polygon is collinear
1796 let(
1797 plane = plane3pt(poly[indices[0]], poly[indices[1]], poly[indices[2]]),
1798 plane_isect = plane_line_intersection(plane, line, bounded, eps)
1799 )
1800 is_undef(plane_isect) ? undef :
1801 is_vector(plane_isect,3) ?
1802 let(
1803 poly2d = project_plane(plane,poly),
1804 pt2d = project_plane(plane, plane_isect)
1805 )
1806 (point_in_polygon(pt2d, poly2d, nonzero=nonzero, eps=eps) < 0 ? undef : plane_isect)
1807 : // Case where line is on the polygon plane
1808 let(
1809 poly2d = project_plane(plane, poly),
1810 line2d = project_plane(plane, line),
1811 segments = polygon_line_intersection(poly2d, line2d, bounded=bounded, nonzero=nonzero, eps=eps)
1812 )
1813 segments==undef ? undef
1814 : [for(seg=segments) len(seg)==2 ? lift_plane(plane,seg) : [lift_plane(plane,seg[0])]];
1815
1816function _merge_segments(insegs,outsegs, eps, i=1) =
1817 i==len(insegs) ? outsegs :
1818 approx(last(last(outsegs)), insegs[i][0], eps)
1819 ? _merge_segments(insegs, [each list_head(outsegs),[last(outsegs)[0],last(insegs[i])]], eps, i+1)
1820 : _merge_segments(insegs, [each outsegs, insegs[i]], eps, i+1);
1821
1822
1823
1824// Function: polygon_triangulate()
1825// Usage:
1826// triangles = polygon_triangulate(poly, [ind], [error], [eps])
1827// Description:
1828// Given a simple polygon in 2D or 3D, triangulates it and returns a list
1829// of triples indexing into the polygon vertices. When the optional argument `ind` is
1830// given, it is used as an index list into `poly` to define the polygon vertices. In that case,
1831// `poly` may have a length greater than `ind`. When `ind` is undefined, all points in `poly`
1832// are considered as vertices of the polygon.
1833// .
1834// For 2d polygons, the output triangles will have the same winding (CW or CCW) of
1835// the input polygon. For 3d polygons, the triangle windings will induce a normal
1836// vector with the same direction of the polygon normal.
1837// .
1838// The function produce correct triangulations for some non-twisted non-simple polygons.
1839// A polygon is non-twisted iff it is simple or it has a partition in
1840// simple polygons with the same winding such that the intersection of any two partitions is
1841// made of full edges and/or vertices of both partitions. These polygons may have "touching" vertices
1842// (two vertices having the same coordinates, but distinct adjacencies) and "contact" edges
1843// (edges whose vertex pairs have the same pairwise coordinates but are in reversed order) but has
1844// no self-crossing. See examples bellow. If all polygon edges are contact edges (polygons with
1845// zero area), it returns an empty list for 2d polygons and reports an error for 3d polygons.
1846// Triangulation errors are reported either by an assert error (when `error=true`) or by returning
1847// `undef` (when `error=false`). Invalid arguments always produce an assert error.
1848// .
1849// Twisted polygons have no consistent winding and when input to this function usually reports
1850// an error but when an error is not reported the outputs are not correct triangulations. The function
1851// can work for 3d non-planar polygons if they are close enough to planar but may otherwise
1852// report an error for this case.
1853// Arguments:
1854// poly = Array of the polygon vertices.
1855// ind = If given, a list of indices indexing the vertices of the polygon in `poly`. Default: use all the points of poly
1856// error = If false, returns `undef` when the polygon cannot be triangulated; otherwise, issues an assert error. Default: true.
1857// eps = A maximum tolerance in geometrical tests. Default: EPSILON
1858// Example(2D,NoAxes): a simple polygon; see from above
1859// poly = star(id=10, od=15,n=11);
1860// tris = polygon_triangulate(poly);
1861// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1862// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1863// color("magenta") up(2) stroke(poly,.25,closed=true);
1864// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1865// Example(2D,NoAxes): a polygon with a hole and one "contact" edge; see from above
1866// poly = [ [-10,0], [10,0], [0,10], [-10,0], [-4,4], [4,4], [0,2], [-4,4] ];
1867// tris = polygon_triangulate(poly);
1868// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1869// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1870// color("magenta") up(2) stroke(poly,.25,closed=true);
1871// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1872// Example(2D,NoAxes): a polygon with "touching" vertices and no holes; see from above
1873// poly = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ];
1874// tris = polygon_triangulate(poly);
1875// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1876// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1877// color("magenta") up(2) stroke(poly,.25,closed=true);
1878// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1879// Example(2D,NoAxes): a polygon with "contact" edges and no holes; see from above
1880// poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3],
1881// [7,7], [7,3], [3,3] ];
1882// tris = polygon_triangulate(poly);
1883// color("lightblue") for(tri=tris) polygon(select(poly,tri));
1884// color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); }
1885// color("magenta") up(2) stroke(poly,.25,closed=true);
1886// color("black") up(3) debug_vnf([path3d(poly),[]],faces=false,size=1);
1887// Example(3D):
1888// include <BOSL2/polyhedra.scad>
1889// vnf = regular_polyhedron_info(name="dodecahedron",side=5,info="vnf");
1890// vnf_polyhedron(vnf);
1891// vnf_tri = [vnf[0], [for(face=vnf[1]) each polygon_triangulate(vnf[0], face) ] ];
1892// color("blue")
1893// vnf_wireframe(vnf_tri, width=.15);
1894function polygon_triangulate(poly, ind, error=true, eps=EPSILON) =
1895 assert(is_path(poly) && len(poly)>=3, "Polygon `poly` should be a list of at least three 2d or 3d points")
1896 assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ),
1897 "Improper or out of bounds list of indices")
1898 let( ind = is_undef(ind) ? count(len(poly)) : ind )
1899 len(ind) <=2 ? [] :
1900 len(ind) == 3
1901 ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] :
1902 // non zero area
1903 let( degen = norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) < 2*eps )
1904 assert( ! error || ! degen, "The polygon vertices are collinear.")
1905 degen ? undef : [ind]
1906 : len(poly[ind[0]]) == 3
1907 ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane
1908 let(
1909 ind = deduplicate_indexed(poly, ind, eps)
1910 )
1911 len(ind)<3 ? [] :
1912 let(
1913 pts = select(poly,ind),
1914 nrm = -polygon_normal(pts)
1915 )
1916 assert( ! error || (nrm != undef),
1917 "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.")
1918 nrm == undef ? undef :
1919 let(
1920 imax = max_index([for(p=pts) norm(p-pts[0]) ]),
1921 v1 = unit( pts[imax] - pts[0] ),
1922 v2 = cross(v1,nrm),
1923 prpts = pts*transpose([v1,v2]) // the 2d projection of pts on the polygon plane
1924 )
1925 let( tris = _triangulate(prpts, count(len(ind)), error, eps) )
1926 tris == undef ? undef :
1927 [for(tri=tris) select(ind,tri) ]
1928 : is_polygon_clockwise(select(poly, ind))
1929 ? _triangulate( poly, ind, error, eps )
1930 : let( tris = _triangulate( poly, reverse(ind), error, eps ) )
1931 tris == undef ? undef :
1932 [for(tri=tris) reverse(tri) ];
1933
1934
1935// poly is supposed to be a 2d cw polygon
1936// implements a modified version of ear cut method for non-twisted polygons
1937// the polygons accepted by this function are those decomposable in simple
1938// CW polygons.
1939function _triangulate(poly, ind, error, eps=EPSILON, tris=[]) =
1940 len(ind)==3
1941 ? _degenerate_tri(select(poly,ind),eps)
1942 ? tris // if last 3 pts perform a degenerate triangle, ignore it
1943 : concat(tris,[ind]) // otherwise, include it
1944 : let( ear = _get_ear(poly,ind,eps) )
1945 assert( ! error || (ear != undef),
1946 "The polygon has twists or all its vertices are collinear or non coplanar.")
1947 ear == undef ? undef :
1948 is_list(ear) // is it a degenerate ear ?
1949 ? len(ind) <= 4 ? tris :
1950 _triangulate(poly, select(ind,ear[0]+3, ear[0]), error, eps, tris) // discard it
1951 : let(
1952 ear_tri = select(ind,ear,ear+2),
1953 indr = select(ind,ear+2, ear) // indices of the remaining path
1954 )
1955 _triangulate(poly, indr, error, eps, concat(tris,[ear_tri]));
1956
1957
1958// a returned ear will be:
1959// 1. a CW non-reflex triangle, made of subsequent poly vertices, without any other
1960// poly points inside except possibly at its own vertices
1961// 2. or a degenerate triangle where two vertices are coincident
1962// the returned ear is specified by the index of `ind` of its first vertex
1963function _get_ear(poly, ind, eps, _i=0) =
1964 let( lind = len(ind) )
1965 lind==3 ? 0 :
1966 let( // the _i-th ear candidate
1967 p0 = poly[ind[_i]],
1968 p1 = poly[ind[(_i+1)%lind]],
1969 p2 = poly[ind[(_i+2)%lind]]
1970 )
1971 // if vertex p1 is a convex candidate to be an ear,
1972 // check if the triangle [p0,p1,p2] contains any other point
1973 // except possibly p0 and p2
1974 // exclude the ear candidate central vertex p1 from the verts to check
1975 _tri_class([p0,p1,p2],eps) > 0
1976 && _none_inside(select(ind,_i+2, _i),poly,p0,p1,p2,eps) ? _i : // found an ear
1977 // otherwise check the next ear candidate
1978 _i<lind-1 ? _get_ear(poly, ind, eps, _i=_i+1) :
1979 // poly has no ears, look for wiskers
1980 let( wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] )
1981 wiskers==[] ? undef : [wiskers[0]];
1982
1983
1984
1985// returns false ASA it finds some reflex vertex of poly[idxs[.]]
1986// inside the triangle different from p0 and p2
1987// note: to simplify the expressions it is assumed that the input polygon has no twists
1988function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) =
1989 i>=len(idxs) ? true :
1990 let(
1991 vert = poly[idxs[i]],
1992 prev_vert = poly[select(idxs,i-1)],
1993 next_vert = poly[select(idxs,i+1)]
1994 )
1995 // check if vert prevent [p0,p1,p2] to be an ear
1996 // this conditions might have a simpler expression
1997 _tri_class([prev_vert, vert, next_vert],eps) <= 0 // reflex condition
1998 && ( // vert is a cw reflex poly vertex inside the triangle [p0,p1,p2]
1999 ( _tri_class([p0,p1,vert],eps)>0 &&
2000 _tri_class([p1,p2,vert],eps)>0 &&
2001 _tri_class([p2,p0,vert],eps)>=0 )
2002 // or it is equal to p1 and some of its adjacent edges cross the open segment (p0,p2)
2003 || ( norm(vert-p1) < eps
2004 && _is_at_left(p0,[prev_vert,p1],eps) && _is_at_left(p2,[p1,prev_vert],eps)
2005 && _is_at_left(p2,[p1,next_vert],eps) && _is_at_left(p0,[next_vert,p1],eps)
2006 )
2007 )
2008 ? false
2009 : _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1);
2010
2011
2012// Function: is_polygon_clockwise()
2013// Usage:
2014// bool = is_polygon_clockwise(poly);
2015// Topics: Geometry, Polygons, Clockwise
2016// See Also: clockwise_polygon(), ccw_polygon(), reverse_polygon()
2017// Description:
2018// Return true if the given 2D simple polygon is in clockwise order, false otherwise.
2019// Results for complex (self-intersecting) polygon are indeterminate.
2020// Arguments:
2021// poly = The list of 2D path points for the perimeter of the polygon.
2022
2023// For algorithm see 2.07 here: http://www.faqs.org/faqs/graphics/algorithms-faq/
2024function is_polygon_clockwise(poly) =
2025 assert(is_path(poly,dim=2), "Input should be a 2d path")
2026 let(
2027 minx = min(poly*[1,0]),
2028 lowind = search(minx, poly, 0, 0),
2029 lowpts = select(poly,lowind),
2030 miny = min(lowpts*[0,1]),
2031 extreme_sub = search(miny, lowpts, 1, 1)[0],
2032 extreme = lowind[extreme_sub]
2033 )
2034 cross(select(poly,extreme+1)-poly[extreme],
2035 select(poly,extreme-1)-poly[extreme])<0;
2036
2037
2038// Function: clockwise_polygon()
2039// Usage:
2040// newpoly = clockwise_polygon(poly);
2041// Topics: Geometry, Polygons, Clockwise
2042// See Also: is_polygon_clockwise(), ccw_polygon(), reverse_polygon()
2043// Description:
2044// Given a 2D polygon path, returns the clockwise winding version of that path.
2045// Arguments:
2046// poly = The list of 2D path points for the perimeter of the polygon.
2047function clockwise_polygon(poly) =
2048 assert(is_path(poly,dim=2), "Input should be a 2d polygon")
2049 is_polygon_clockwise(poly) ? poly : reverse_polygon(poly);
2050
2051
2052// Function: ccw_polygon()
2053// Usage:
2054// newpoly = ccw_polygon(poly);
2055// See Also: is_polygon_clockwise(), clockwise_polygon(), reverse_polygon()
2056// Topics: Geometry, Polygons, Clockwise
2057// Description:
2058// Given a 2D polygon poly, returns the counter-clockwise winding version of that poly.
2059// Arguments:
2060// poly = The list of 2D path points for the perimeter of the polygon.
2061function ccw_polygon(poly) =
2062 assert(is_path(poly,dim=2), "Input should be a 2d polygon")
2063 is_polygon_clockwise(poly) ? reverse_polygon(poly) : poly;
2064
2065
2066// Function: reverse_polygon()
2067// Usage:
2068// newpoly = reverse_polygon(poly)
2069// Topics: Geometry, Polygons, Clockwise
2070// See Also: is_polygon_clockwise(), ccw_polygon(), clockwise_polygon()
2071// Description:
2072// Reverses a polygon's winding direction, while still using the same start point.
2073// Arguments:
2074// poly = The list of the path points for the perimeter of the polygon.
2075function reverse_polygon(poly) =
2076 let(poly=force_path(poly,"poly"))
2077 assert(is_path(poly), "Input should be a polygon")
2078 [ poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ];
2079
2080
2081// Function: reindex_polygon()
2082// Usage:
2083// newpoly = reindex_polygon(reference, poly);
2084// Topics: Geometry, Polygons
2085// Description:
2086// Rotates and possibly reverses the point order of a 2d or 3d polygon path to optimize its pairwise point
2087// association with a reference polygon. The two polygons must have the same number of vertices and be the same dimension.
2088// The optimization is done by computing the distance, norm(reference[i]-poly[i]), between
2089// corresponding pairs of vertices of the two polygons and choosing the polygon point index rotation that
2090// makes the total sum over all pairs as small as possible. Returns the reindexed polygon. Note
2091// that the geometry of the polygon is not changed by this operation, just the labeling of its
2092// vertices. If the input polygon is 2d and is oriented opposite the reference then its point order is
2093// reversed.
2094// Arguments:
2095// reference = reference polygon path
2096// poly = input polygon to reindex
2097// Example(2D): The red dots show the 0th entry in the two input path lists. Note that the red dots are not near each other. The blue dot shows the 0th entry in the output polygon
2098// pent = subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],30);
2099// circ = circle($fn=30,r=2.2);
2100// reindexed = reindex_polygon(circ,pent);
2101// move_copies(concat(circ,pent)) circle(r=.1,$fn=32);
2102// color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
2103// color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
2104// Example(2D): The indexing that minimizes the total distance will not necessarily associate the nearest point of `poly` with the reference, as in this example where again the blue dot indicates the 0th entry in the reindexed result.
2105// pent = move([3.5,-1],p=subdivide_path([for(i=[0:4])[sin(72*i),cos(72*i)]],30));
2106// circ = circle($fn=30,r=2.2);
2107// reindexed = reindex_polygon(circ,pent);
2108// move_copies(concat(circ,pent)) circle(r=.1,$fn=32);
2109// color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
2110// color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
2111function reindex_polygon(reference, poly, return_error=false) =
2112 let(reference=force_path(reference,"reference"),
2113 poly=force_path(poly,"poly"))
2114 assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
2115 "Invalid polygon(s) or incompatible dimensions. " )
2116 assert(len(reference)==len(poly), "The polygons must have the same length.")
2117 let(
2118 dim = len(reference[0]),
2119 N = len(reference),
2120 fixpoly = dim != 2? poly :
2121 is_polygon_clockwise(reference)
2122 ? clockwise_polygon(poly)
2123 : ccw_polygon(poly),
2124 I = [for(i=reference) 1],
2125 val = [ for(k=[0:N-1])
2126 [for(i=[0:N-1])
2127 norm(reference[i]-fixpoly[(i+k)%N]) ] ]*I,
2128 min_ind = min_index(val),
2129 optimal_poly = list_rotate(fixpoly, min_ind)
2130 )
2131 return_error? [optimal_poly, val[min_ind]] :
2132 optimal_poly;
2133
2134
2135// Function: align_polygon()
2136// Usage:
2137// newpoly = align_polygon(reference, poly, [angles], [cp], [tran], [return_ind]);
2138// Topics: Geometry, Polygons
2139// Description:
2140// Find the best alignment of a specified 2D polygon with a reference 2D polygon over a set of
2141// transformations. You can specify a list or range of angles and a centerpoint or you can
2142// give a list of arbitrary 2d transformation matrices. For each transformation or angle, the polygon is
2143// reindexed, which is a costly operation so if run time is a problem, use a smaller sampling of angles or
2144// transformations. By default returns the rotated and reindexed polygon. You can also request that
2145// the best angle or the index into the transformation list be returned. When you specify an angle
2146// Arguments:
2147// reference = reference polygon
2148// poly = polygon to rotate into alignment with the reference
2149// angles = list or range of angles to test
2150// cp = centerpoint for rotations
2151// ---
2152// tran = list of 2D transformation matrices to optimize over
2153// return_ind = if true, return the best angle (if you specified angles) or the index into tran otherwise of best alignment
2154// Example(2D): Rotating the poorly aligned light gray triangle by 105 degrees produces the best alignment, shown in blue:
2155// ellipse = yscale(3,circle(r=10, $fn=32));
2156// tri = move([-50/3,-9],
2157// subdivide_path([[0,0], [50,0], [0,27]], 32));
2158// aligned = align_polygon(ellipse,tri, [0:5:180]);
2159// color("white")stroke(tri,width=.5,closed=true);
2160// stroke(ellipse, width=.5, closed=true);
2161// color("blue")stroke(aligned,width=.5,closed=true);
2162// Example(2D,NoAxes): Translating a triangle (light gray) to the best alignment (blue)
2163// ellipse = yscale(2,circle(r=10, $fn=32));
2164// tri = subdivide_path([[0,0], [27,0], [-7,50]], 32);
2165// T = [for(x=[-10:0], y=[-30:-15]) move([x,y])];
2166// aligned = align_polygon(ellipse,tri, trans=T);
2167// color("white")stroke(tri,width=.5,closed=true);
2168// stroke(ellipse, width=.5, closed=true);
2169// color("blue")stroke(aligned,width=.5,closed=true);
2170function align_polygon(reference, poly, angles, cp, trans, return_ind=false) =
2171 let(reference=force_path(reference,"reference"),
2172 poly=force_path(poly,"poly"))
2173 assert(is_undef(trans) || (is_undef(angles) && is_undef(cp)), "Cannot give both angles/cp and trans as input")
2174 let(
2175 trans = is_def(trans) ? trans :
2176 assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
2177 "The `angle` parameter must be a range or a non void list of numbers.")
2178 [for(angle=angles) zrot(angle,cp=cp)]
2179 )
2180 assert(is_path(reference,dim=2), "reference must be a 2D polygon")
2181 assert(is_path(poly,dim=2), "poly must be a 2D polygon")
2182 assert(len(reference)==len(poly), "The polygons must have the same length.")
2183 let( // alignments is a vector of entries of the form: [polygon, error]
2184 alignments = [
2185 for(T=trans)
2186 reindex_polygon(
2187 reference,
2188 apply(T,poly),
2189 return_error=true
2190 )
2191 ],
2192 scores = column(alignments,1),
2193 minscore = min(scores),
2194 minind = [for(i=idx(scores)) if (scores[i]<minscore+EPSILON) i],
2195 dummy = is_def(angles) ? echo(best_angles = select(list(angles), minind)):0,
2196 best = minind[0]
2197 )
2198 return_ind ? (is_def(angles) ? list(angles)[best] : best)
2199 : alignments[best][0];
2200
2201
2202// Function: are_polygons_equal()
2203// Usage:
2204// bool = are_polygons_equal(poly1, poly2, [eps])
2205// Description:
2206// Returns true if poly1 and poly2 are the same polongs
2207// within given epsilon tolerance.
2208// Arguments:
2209// poly1 = first polygon
2210// poly2 = second polygon
2211// eps = tolerance for comparison
2212// Example(NORENDER):
2213// are_polygons_equal(pentagon(r=4),
2214// rot(360/5, p=pentagon(r=4))); // returns true
2215// are_polygons_equal(pentagon(r=4),
2216// rot(90, p=pentagon(r=4))); // returns false
2217function are_polygons_equal(poly1, poly2, eps=EPSILON) =
2218 let(
2219 poly1 = cleanup_path(poly1),
2220 poly2 = cleanup_path(poly2),
2221 l1 = len(poly1),
2222 l2 = len(poly2)
2223 ) l1 != l2 ? false :
2224 let( maybes = find_approx(poly1[0], poly2, eps=eps, all=true) )
2225 maybes == []? false :
2226 [for (i=maybes) if (_are_polygons_equal(poly1, poly2, eps, i)) 1] != [];
2227
2228function _are_polygons_equal(poly1, poly2, eps, st) =
2229 max([for(d=poly1-select(poly2,st,st-1)) d*d])<eps*eps;
2230
2231
2232/// Function: _is_polygon_in_list()
2233/// Topics: Polygons, Comparators
2234/// See Also: are_polygons_equal(), are_regions_equal()
2235/// Usage:
2236/// bool = _is_polygon_in_list(poly, polys);
2237/// Description:
2238/// Returns true if one of the polygons in `polys` is equivalent to the polygon `poly`.
2239/// Arguments:
2240/// poly = The polygon to search for.
2241/// polys = The list of polygons to look for the polygon in.
2242function _is_polygon_in_list(poly, polys) =
2243 ___is_polygon_in_list(poly, polys, 0);
2244
2245function ___is_polygon_in_list(poly, polys, i) =
2246 i >= len(polys)? false :
2247 are_polygons_equal(poly, polys[i])? true :
2248 ___is_polygon_in_list(poly, polys, i+1);
2249
2250
2251// Section: Convex Hull
2252
2253// This section originally based on Oskar Linde's Hull:
2254// - https://github.com/openscad/scad-utils
2255
2256
2257// Function: hull()
2258// Usage:
2259// face_list_or_index_list = hull(points);
2260// Description:
2261// Takes a list of 2D or 3D points (but not both in the same list) and returns either the list of
2262// indexes into `points` that forms the 2D convex hull perimeter path, or the list of faces that
2263// form the 3d convex hull surface. Each face is a list of indexes into `points`. If the input
2264// points are co-linear, the result will be the indexes of the two extrema points. If the input
2265// points are co-planar, the results will be a simple list of vertex indices that will form a planar
2266// perimeter. Otherwise a list of faces will be returned, where each face is a simple list of
2267// vertex indices for the perimeter of the face.
2268// Arguments:
2269// points = The set of 2D or 3D points to find the hull of.
2270function hull(points) =
2271 assert(is_path(points),"Invalid input to hull")
2272 len(points[0]) == 2
2273 ? hull2d_path(points)
2274 : hull3d_faces(points);
2275
2276
2277// Module: hull_points()
2278// Usage:
2279// hull_points(points, [fast]);
2280// Description:
2281// If given a list of 2D points, creates a 2D convex hull polygon that encloses all those points.
2282// If given a list of 3D points, creates a 3D polyhedron that encloses all the points. This should
2283// handle about 4000 points in slow mode. If `fast` is set to true, this should be able to handle
2284// far more. When fast mode is off, 3d hulls that lie in a plane will produce a single face of a polyhedron, which can be viewed in preview but will not render.
2285// Arguments:
2286// points = The list of points to form a hull around.
2287// fast = If true for 3d case, uses a faster cheat that may handle more points, but also may emit warnings that can stop your script if you have "Halt on first warning" enabled. Ignored for the 2d case. Default: false
2288// Example(2D):
2289// pts = [[-10,-10], [0,10], [10,10], [12,-10]];
2290// hull_points(pts);
2291// Example(3D):
2292// pts = [for (phi = [30:60:150], theta = [0:60:359]) spherical_to_xyz(10, theta, phi)];
2293// hull_points(pts);
2294module hull_points(points, fast=false) {
2295 no_children($children);
2296 check = assert(is_path(points))
2297 assert(len(points)>=3, "Point list must contain 3 points");
2298 if (len(points[0])==2)
2299 hull() polygon(points=points);
2300 else {
2301 if (fast) {
2302 extra = len(points)%3;
2303 faces = [
2304 [for(i=[0:1:extra+2])i], // If vertex count not divisible by 3, combine extras with first 3
2305 for(i=[extra+3:3:len(points)-3])[i,i+1,i+2]
2306 ];
2307 hull() polyhedron(points=points, faces=faces);
2308 } else {
2309 faces = hull(points);
2310 if (is_num(faces[0])){
2311 if (len(faces)<=2) echo("Hull contains only two points");
2312 else polyhedron(points=points, faces=[faces]);
2313 }
2314 else polyhedron(points=points, faces=faces);
2315 }
2316 }
2317}
2318
2319
2320
2321function _backtracking(i,points,h,t,m,all) =
2322 m<t || _is_cw(points[i], points[h[m-1]], points[h[m-2]],all) ? m :
2323 _backtracking(i,points,h,t,m-1,all) ;
2324
2325// clockwise check (2d)
2326function _is_cw(a,b,c,all) =
2327 all ? cross(a-c,b-c)<=EPSILON*norm(a-c)*norm(b-c) :
2328 cross(a-c,b-c)<-EPSILON*norm(a-c)*norm(b-c);
2329
2330
2331// Function: hull2d_path()
2332// Usage:
2333// index_list = hull2d_path(points,all)
2334// Description:
2335// Takes a list of arbitrary 2D points, and finds the convex hull polygon to enclose them.
2336// Returns a path as a list of indices into `points`.
2337// When all==true, returns extra points that are on edges of the hull.
2338// Arguments:
2339// points = list of 2d points to get the hull of.
2340// all = when true, includes all points on the edges of the convex hull. Default: false.
2341// Example(2D):
2342// pts = [[-10,-10], [0,10], [10,10], [12,-10]];
2343// path = hull2d_path(pts);
2344// move_copies(pts) color("red") circle(1,$fn=12);
2345// polygon(points=pts, paths=[path]);
2346//
2347// Code based on this method:
2348// https://www.hackerearth.com/practice/math/geometry/line-sweep-technique/tutorial/
2349//
2350function hull2d_path(points, all=false) =
2351 assert(is_path(points,2),"Invalid input to hull2d_path")
2352 len(points) < 2 ? [] :
2353 let( n = len(points),
2354 ip = sortidx(points) )
2355 // lower hull points
2356 let( lh =
2357 [ for( i = 2,
2358 k = 2,
2359 h = [ip[0],ip[1]]; // current list of hull point indices
2360 i <= n;
2361 k = i<n ? _backtracking(ip[i],points,h,2,k,all)+1 : k,
2362 h = i<n ? [for(j=[0:1:k-2]) h[j], ip[i]] : [],
2363 i = i+1
2364 ) if( i==n ) h ][0] )
2365 // concat lower hull points with upper hull ones
2366 [ for( i = n-2,
2367 k = len(lh),
2368 t = k+1,
2369 h = lh; // current list of hull point indices
2370 i >= -1;
2371 k = i>=0 ? _backtracking(ip[i],points,h,t,k,all)+1 : k,
2372 h = [for(j=[0:1:k-2]) h[j], if(i>0) ip[i]],
2373 i = i-1
2374 ) if( i==-1 ) h ][0] ;
2375
2376
2377function _hull_collinear(points) =
2378 let(
2379 a = points[0],
2380 i = max_index([for(pt=points) norm(pt-a)]),
2381 n = points[i] - a
2382 )
2383 norm(n)==0 ? [0]
2384 :
2385 let(
2386 points1d = [ for(p = points) (p-a)*n ],
2387 min_i = min_index(points1d),
2388 max_i = max_index(points1d)
2389 ) [min_i, max_i];
2390
2391
2392
2393// Function: hull3d_faces()
2394// Usage:
2395// faces = hull3d_faces(points)
2396// Description:
2397// Takes a list of arbitrary 3D points, and finds the convex hull polyhedron to enclose
2398// them. Returns a list of triangular faces, where each face is a list of indexes into the given `points`
2399// list. The output will be valid for use with the polyhedron command, but may include vertices that are in the interior of a face of the hull, so it is not
2400// necessarily the minimal representation of the hull.
2401// If all points passed to it are coplanar, then the return is the list of indices of points
2402// forming the convex hull polygon.
2403// Example(3D):
2404// pts = [[-20,-20,0], [20,-20,0], [0,20,5], [0,0,20]];
2405// faces = hull3d_faces(pts);
2406// move_copies(pts) color("red") sphere(1);
2407// %polyhedron(points=pts, faces=faces);
2408function hull3d_faces(points) =
2409 assert(is_path(points,3),"Invalid input to hull3d_faces")
2410 len(points) < 3 ? count(len(points))
2411 : let ( // start with a single non-collinear triangle
2412 tri = _noncollinear_triple(points, error=false)
2413 )
2414 tri==[] ? _hull_collinear(points)
2415 : let(
2416 a = tri[0],
2417 b = tri[1],
2418 c = tri[2],
2419 plane = plane3pt_indexed(points, a, b, c),
2420 d = _find_first_noncoplanar(plane, points)
2421 )
2422 d == len(points)
2423 ? /* all coplanar*/
2424 let (
2425 pts2d = project_plane([points[a], points[b], points[c]],points),
2426 hull2d = hull2d_path(pts2d)
2427 ) hull2d
2428 : let(
2429 remaining = [for (i = [0:1:len(points)-1]) if (i!=a && i!=b && i!=c && i!=d) i],
2430 // Build an initial tetrahedron.
2431 // Swap b, c if d is in front of triangle t.
2432 ifop = _is_point_above_plane(plane, points[d]),
2433 bc = ifop? [c,b] : [b,c],
2434 b = bc[0],
2435 c = bc[1],
2436 triangles = [
2437 [a,b,c],
2438 [d,b,a],
2439 [c,d,a],
2440 [b,d,c]
2441 ],
2442 // calculate the plane equations
2443 planes = [ for (t = triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
2444 ) _hull3d_iterative(points, triangles, planes, remaining);
2445
2446
2447// Adds the remaining points one by one to the convex hull
2448function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = //let( EPSILON=1e-12 )
2449 _i >= len(remaining) ? triangles :
2450 let (
2451 // pick a point
2452 i = remaining[_i],
2453 // evaluate the triangle plane equations at point i
2454 planeq_val = planes*[each points[i], -1],
2455 // find the triangles that are in conflict with the point (point not inside)
2456 conflicts = [for (i = [0:1:len(planeq_val)-1]) if (planeq_val[i]>EPSILON) i ],
2457 // collect the halfedges of all triangles that are in conflict
2458 halfedges = [
2459 for(c = conflicts, i = [0:2])
2460 [triangles[c][i], triangles[c][(i+1)%3]]
2461 ],
2462 // find the outer perimeter of the set of conflicting triangles
2463 horizon = _remove_internal_edges(halfedges),
2464 // generate new triangles connecting point i to each horizon halfedge vertices
2465 tri2add = [ for (h = horizon) concat(h,i) ],
2466 // add tria2add and remove conflict triangles
2467 new_triangles =
2468 concat( tri2add,
2469 [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) triangles[i] ]
2470 ),
2471 // add the plane equations of new added triangles and remove the plane equations of the conflict ones
2472 new_planes =
2473 [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ,
2474 for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ]
2475 ) _hull3d_iterative(
2476 points,
2477 new_triangles,
2478 new_planes,
2479 remaining,
2480 _i+1
2481 );
2482
2483
2484function _remove_internal_edges(halfedges) = [
2485 for (h = halfedges)
2486 if (!in_list(reverse(h), halfedges))
2487 h
2488];
2489
2490function _find_first_noncoplanar(plane, points, i=0) =
2491 (i >= len(points) || !are_points_on_plane([points[i]],plane))? i :
2492 _find_first_noncoplanar(plane, points, i+1);
2493
2494
2495// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
2496
2497
2498
2499
2500
2501// Section: Convex Sets
2502
2503
2504// Function: is_polygon_convex()
2505// Usage:
2506// bool = is_polygon_convex(poly, [eps]);
2507// Topics: Geometry, Convexity, Test
2508// Description:
2509// Returns true if the given 2D or 3D polygon is convex.
2510// The result is meaningless if the polygon is not simple (self-crossing) or non coplanar.
2511// If the points are collinear or not coplanar an error may be generated.
2512// Arguments:
2513// poly = Polygon to check.
2514// eps = Tolerance for the collinearity and coplanarity tests. Default: EPSILON.
2515// Example:
2516// test1 = is_polygon_convex(circle(d=50)); // Returns: true
2517// test2 = is_polygon_convex(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
2518// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
2519// test = is_polygon_convex(spiral); // Returns: false
2520function is_polygon_convex(poly,eps=EPSILON) =
2521 assert(is_path(poly), "The input should be a 2D or 3D polygon." )
2522 let(
2523 lp = len(poly),
2524 p0 = poly[0]
2525 )
2526 assert( lp>=3 , "A polygon must have at least 3 points" )
2527 let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
2528 len(p0)==2
2529 ? let( size = max([for(p=poly) norm(p-p0)]), tol=pow(size,2)*eps )
2530 assert( size>eps, "The polygon is self-crossing or its points are collinear" )
2531 min(crosses) >=-tol || max(crosses)<=tol
2532 : let( ip = _noncollinear_triple(poly,error=false,eps=eps) )
2533 assert( ip!=[], "The points are collinear")
2534 let(
2535 crx = cross(poly[ip[1]]-poly[ip[0]],poly[ip[2]]-poly[ip[1]]),
2536 nrm = crx/norm(crx),
2537 plane = concat(nrm, nrm*poly[0]),
2538 prod = crosses*nrm,
2539 size = norm(poly[ip[1]]-poly[ip[0]]),
2540 tol = pow(size,2)*eps
2541 )
2542 assert(_pointlist_greatest_distance(poly,plane) < size*eps, "The polygon points are not coplanar")
2543 let(
2544 minc = min(prod),
2545 maxc = max(prod) )
2546 minc>=-tol || maxc<=tol;
2547
2548
2549// Function: convex_distance()
2550// Usage:
2551// dist = convex_distance(points1, points2,eps);
2552// Topics: Geometry, Convexity, Distance
2553// See also:
2554// convex_collision(), hull()
2555// Description:
2556// Returns the smallest distance between a point in convex hull of `points1`
2557// and a point in the convex hull of `points2`. All the points in the lists
2558// should have the same dimension, either 2D or 3D.
2559// A zero result means the hulls intercept whithin a tolerance `eps`.
2560// Arguments:
2561// points1 = first list of 2d or 3d points.
2562// points2 = second list of 2d or 3d points.
2563// eps = tolerance in distance evaluations. Default: EPSILON.
2564// Example(2D):
2565// pts1 = move([-3,0], p=square(3,center=true));
2566// pts2 = rot(a=45, p=square(2,center=true));
2567// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
2568// polygon(pts1);
2569// polygon(pts2);
2570// polygon(pts3);
2571// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864
2572// echo(convex_distance(pts2,pts3)); // Returns: 0
2573// Example(3D):
2574// sphr1 = sphere(2,$fn=10);
2575// sphr2 = move([4,0,0], p=sphr1);
2576// sphr3 = move([4.5,0,0], p=sphr1);
2577// vnf_polyhedron(sphr1);
2578// vnf_polyhedron(sphr2);
2579// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0
2580// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5
2581function convex_distance(points1, points2, eps=EPSILON) =
2582 assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
2583 "The input lists should be compatible consistent non empty lists of points.")
2584 assert(len(points1[0])==2 || len(points1[0])==3 ,
2585 "The input points should be 2d or 3d points.")
2586 let( d = points1[0]-points2[0] )
2587 norm(d)<eps ? 0 :
2588 let( v = _support_diff(points1,points2,-d) )
2589 norm(_GJK_distance(points1, points2, eps, 0, v, [v]));
2590
2591
2592// Finds the vector difference between the hulls of the two pointsets by the GJK algorithm
2593// Based on:
2594// http://www.dtecta.com/papers/jgt98convex.pdf
2595function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) =
2596 let( nrd = norm(d) ) // distance upper bound
2597 nrd<eps ? d :
2598 let(
2599 v = _support_diff(points1,points2,-d),
2600 lbd = max(lbd, d*v/nrd), // distance lower bound
2601 close = (nrd-lbd <= eps*nrd)
2602 )
2603 close ? d :
2604 let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
2605 _GJK_distance(points1, points2, eps, lbd, newsplx[0], newsplx[1]);
2606
2607
2608// Function: convex_collision()
2609// Usage:
2610// bool = convex_collision(points1, points2, [eps]);
2611// Topics: Geometry, Convexity, Collision, Intersection
2612// See also:
2613// convex_distance(), hull()
2614// Description:
2615// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2`
2616// otherwise, `false`.
2617// All the points in the lists should have the same dimension, either 2D or 3D.
2618// This function is tipically faster than `convex_distance` to find a non-collision.
2619// Arguments:
2620// points1 = first list of 2d or 3d points.
2621// points2 = second list of 2d or 3d points.
2622// eps - tolerance for the intersection tests. Default: EPSILON.
2623// Example(2D):
2624// pts1 = move([-3,0], p=square(3,center=true));
2625// pts2 = rot(a=45, p=square(2,center=true));
2626// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
2627// polygon(pts1);
2628// polygon(pts2);
2629// polygon(pts3);
2630// echo(convex_collision(pts1,pts2)); // Returns: false
2631// echo(convex_collision(pts2,pts3)); // Returns: true
2632// Example(3D):
2633// sphr1 = sphere(2,$fn=10);
2634// sphr2 = move([4,0,0], p=sphr1);
2635// sphr3 = move([4.5,0,0], p=sphr1);
2636// vnf_polyhedron(sphr1);
2637// vnf_polyhedron(sphr2);
2638// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true
2639// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false
2640//
2641function convex_collision(points1, points2, eps=EPSILON) =
2642 assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
2643 "The input lists should be compatible consistent non empty lists of points.")
2644 assert(len(points1[0])==2 || len(points1[0])==3 ,
2645 "The input points should be 2d or 3d points.")
2646 let( d = points1[0]-points2[0] )
2647 norm(d)<eps ? true :
2648 let( v = _support_diff(points1,points2,-d) )
2649 _GJK_collide(points1, points2, v, [v], eps);
2650
2651
2652// Based on the GJK collision algorithms found in:
2653// http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
2654// or
2655// http://www.dtecta.com/papers/jgt98convex.pdf
2656function _GJK_collide(points1, points2, d, simplex, eps=EPSILON) =
2657 norm(d) < eps ? true : // does collide
2658 let( v = _support_diff(points1,points2,-d) )
2659 v*d > eps*eps ? false : // no collision
2660 let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
2661 norm(v-newsplx[0])<eps ? norm(v)<eps :
2662 _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps);
2663
2664
2665// given a simplex s, returns a pair:
2666// - the point of the s closest to the origin
2667// - the smallest sub-simplex of s that contains that point
2668function _closest_simplex(s,eps=EPSILON) =
2669 len(s)==2 ? _closest_s1(s,eps) :
2670 len(s)==3 ? _closest_s2(s,eps) :
2671 len(s)==4 ? _closest_s3(s,eps) :
2672 assert(false, "Internal error.");
2673
2674
2675// find the point of a 1-simplex closest to the origin
2676function _closest_s1(s,eps=EPSILON) =
2677 norm(s[1]-s[0])<=eps*(norm(s[0])+norm(s[1]))/2 ? [ s[0], [s[0]] ] :
2678 let(
2679 c = s[1]-s[0],
2680 t = -s[0]*c/(c*c)
2681 )
2682 t<0 ? [ s[0], [s[0]] ] :
2683 t>1 ? [ s[1], [s[1]] ] :
2684 [ s[0]+t*c, s ];
2685
2686
2687// find the point of a 2-simplex closest to the origin
2688function _closest_s2(s, eps=EPSILON) =
2689 // considering that s[2] was the last inserted vertex in s by GJK,
2690 // the plane orthogonal to the triangle [ origin, s[0], s[1] ] that
2691 // contains [s[0],s[1]] have the origin and s[2] on the same side;
2692 // that reduces the cases to test and the only possible simplex
2693 // outcomes are s, [s[0],s[2]] and [s[1],s[2]]
2694 let(
2695 area = cross(s[2]-s[0], s[1]-s[0]),
2696 area2 = area*area // tri area squared
2697 )
2698 area2<=eps*max([for(si=s) pow(si*si,2)]) // degenerate tri
2699 ? norm(s[2]-s[0]) < norm(s[2]-s[1])
2700 ? _closest_s1([s[1],s[2]])
2701 : _closest_s1([s[0],s[2]])
2702 : let(
2703 crx1 = cross(s[0], s[2])*area,
2704 crx2 = cross(s[1], s[0])*area,
2705 crx0 = cross(s[2], s[1])*area
2706 )
2707 // all have the same signal -> origin projects inside the tri
2708 max(crx1, crx0, crx2) < 0 || min(crx1, crx0, crx2) > 0
2709 ? // baricentric coords of projection
2710 [ [abs(crx0),abs(crx1),abs(crx2)]*s/area2, s ]
2711 : let(
2712 cl12 = _closest_s1([s[1],s[2]]),
2713 cl02 = _closest_s1([s[0],s[2]])
2714 )
2715 norm(cl12[0])<norm(cl02[0]) ? cl12 : cl02;
2716
2717
2718// find the point of a 3-simplex closest to the origin
2719function _closest_s3(s,eps=EPSILON) =
2720 let( nr = cross(s[1]-s[0],s[2]-s[0]),
2721 sz = [ norm(s[0]-s[1]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] )
2722 norm(nr)<=eps*pow(max(sz),2)
2723 ? let( i = max_index(sz) )
2724 _closest_s2([ s[i], s[(i+1)%3], s[3] ], eps) // degenerate case
2725 : // considering that s[3] was the last inserted vertex in s by GJK,
2726 // the only possible outcomes will be:
2727 // s or some of the 3 faces of s containing s[3]
2728 let(
2729 tris = [ [s[0], s[1], s[3]],
2730 [s[1], s[2], s[3]],
2731 [s[2], s[0], s[3]] ],
2732 cntr = sum(s)/4,
2733 // indicator of the tris facing the origin
2734 facing = [for(i=[0:2])
2735 let( nrm = _tri_normal(tris[i]) )
2736 if( ((nrm*(s[i]-cntr))>0)==(nrm*s[i]<0) ) i ]
2737 )
2738 len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex
2739 len(facing)==1 ? _closest_s2(tris[facing[0]], eps) :
2740 let( // look for the origin-facing tri closest to the origin
2741 closest = [for(i=facing) _closest_s2(tris[i], eps) ],
2742 dist = [for(cl=closest) norm(cl[0]) ],
2743 nearest = min_index(dist)
2744 )
2745 closest[nearest];
2746
2747
2748function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]);
2749
2750
2751function _support_diff(p1,p2,d) =
2752 let( p1d = p1*d, p2d = p2*d )
2753 p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]];
2754
2755
2756// Section: Rotation Decoding
2757
2758// Function: rot_decode()
2759// Usage:
2760// info = rot_decode(rotation,[long]); // Returns: [angle,axis,cp,translation]
2761// Topics: Affine, Matrices, Transforms
2762// Description:
2763// Given an input 3D rigid transformation operator (one composed of just rotations and translations) represented
2764// as a 4x4 matrix, compute the rotation and translation parameters of the operator. Returns a list of the
2765// four parameters, the angle, in the interval [0,180], the rotation axis as a unit vector, a centerpoint for
2766// the rotation, and a translation. If you set `parms = rot_decode(rotation)` then the transformation can be
2767// reconstructed from parms as `move(parms[3]) * rot(a=parms[0],v=parms[1],cp=parms[2])`. This decomposition
2768// makes it possible to perform interpolation. If you construct a transformation using `rot` the decoding
2769// may flip the axis (if you gave an angle outside of [0,180]). The returned axis will be a unit vector, and
2770// the centerpoint lies on the plane through the origin that is perpendicular to the axis. It may be different
2771// than the centerpoint you used to construct the transformation.
2772// .
2773// If you set `long` to true then return the reversed rotation, with the angle in [180,360].
2774// Arguments:
2775// rotation = rigid transformation to decode
2776// long = if true return the "long way" around, with the angle in [180,360]. Default: false
2777// Example:
2778// info = rot_decode(rot(45));
2779// // Returns: [45, [0,0,1], [0,0,0], [0,0,0]]
2780// info = rot_decode(rot(a=37, v=[1,2,3], cp=[4,3,-7])));
2781// // Returns: [37, [0.26, 0.53, 0.80], [4.8, 4.6, -4.6], [0,0,0]]
2782// info = rot_decode(left(12)*xrot(-33));
2783// // Returns: [33, [-1,0,0], [0,0,0], [-12,0,0]]
2784// info = rot_decode(translate([3,4,5]));
2785// // Returns: [0, [0,0,1], [0,0,0], [3,4,5]]
2786function rot_decode(M,long=false) =
2787 assert(is_matrix(M,4,4) && approx(M[3],[0,0,0,1]), "Input matrix must be a 4x4 matrix representing a 3d transformation")
2788 let(R = submatrix(M,[0:2],[0:2]))
2789 assert(approx(det3(R),1) && approx(norm_fro(R * transpose(R)-ident(3)),0),"Input matrix is not a rotation")
2790 let(
2791 translation = [for(row=[0:2]) M[row][3]], // translation vector
2792 largest = max_index([R[0][0], R[1][1], R[2][2]]),
2793 axis_matrix = R + transpose(R) - (matrix_trace(R)-1)*ident(3), // Each row is on the rotational axis
2794 // Construct quaternion q = c * [x sin(theta/2), y sin(theta/2), z sin(theta/2), cos(theta/2)]
2795 q_im = axis_matrix[largest],
2796 q_re = R[(largest+2)%3][(largest+1)%3] - R[(largest+1)%3][(largest+2)%3],
2797 c_sin = norm(q_im), // c * sin(theta/2) for some c
2798 c_cos = abs(q_re) // c * cos(theta/2)
2799 )
2800 approx(c_sin,0) ? [0,[0,0,1],[0,0,0],translation] :
2801 let(
2802 angle = 2*atan2(c_sin, c_cos), // This is supposed to be more accurate than acos or asin
2803 axis = (q_re>=0 ? 1:-1)*q_im/c_sin,
2804 tproj = translation - (translation*axis)*axis, // Translation perpendicular to axis determines centerpoint
2805 cp = (tproj + cross(axis,tproj)*c_cos/c_sin)/2
2806 )
2807 [long ? 360-angle:angle,
2808 long? -axis : axis,
2809 cp,
2810 (translation*axis)*axis];
2811
2812
2813
2814
2815// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap