1//////////////////////////////////////////////////////////////////////
2// LibFile: math.scad
3// Math helper functions.
4// To use, add the following lines to the beginning of your file:
5// ```
6// use <BOSL/math.scad>
7// ```
8//////////////////////////////////////////////////////////////////////
9
10/*
11BSD 2-Clause License
12
13Copyright (c) 2017, Revar Desmera
14All rights reserved.
15
16Redistribution and use in source and binary forms, with or without
17modification, are permitted provided that the following conditions are met:
18
19* Redistributions of source code must retain the above copyright notice, this
20 list of conditions and the following disclaimer.
21
22* Redistributions in binary form must reproduce the above copyright notice,
23 this list of conditions and the following disclaimer in the documentation
24 and/or other materials provided with the distribution.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
27AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
29DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
30FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
31DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
32SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
33CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
34OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36*/
37
38
39include <constants.scad>
40include <compat.scad>
41
42
43// Section: Math Constants
44
45PHI = (1+sqrt(5))/2; // The golden ratio phi.
46
47EPSILON = 1e-9; // A really small value useful in comparing FP numbers. ie: abs(a-b)<EPSILON
48
49
50
51// Function: Cpi()
52// Status: DEPRECATED, use `PI` instead.
53// Description:
54// Returns the value of pi.
55function Cpi() = PI; // Deprecated! Use the variable PI instead.
56
57
58// Section: Simple Calculations
59
60// Function: quant()
61// Description:
62// Quantize a value `x` to an integer multiple of `y`, rounding to the nearest multiple.
63// Arguments:
64// x = The value to quantize.
65// y = The multiple to quantize to.
66function quant(x,y) = floor(x/y+0.5)*y;
67
68
69// Function: quantdn()
70// Description:
71// Quantize a value `x` to an integer multiple of `y`, rounding down to the previous multiple.
72// Arguments:
73// x = The value to quantize.
74// y = The multiple to quantize to.
75function quantdn(x,y) = floor(x/y)*y;
76
77
78// Function: quantup()
79// Description:
80// Quantize a value `x` to an integer multiple of `y`, rounding up to the next multiple.
81// Arguments:
82// x = The value to quantize.
83// y = The multiple to quantize to.
84function quantup(x,y) = ceil(x/y)*y;
85
86
87// Function: constrain()
88// Usage:
89// constrain(v, minval, maxval);
90// Description:
91// Constrains value to a range of values between minval and maxval, inclusive.
92// Arguments:
93// v = value to constrain.
94// minval = minimum value to return, if out of range.
95// maxval = maximum value to return, if out of range.
96function constrain(v, minval, maxval) = min(maxval, max(minval, v));
97
98
99// Function: min_index()
100// Usage:
101// min_index(vals);
102// Description:
103// Returns the index of the minimal value in the given list.
104function min_index(vals, _minval, _minidx, _i=0) =
105 _i>=len(vals)? _minidx :
106 min_index(
107 vals,
108 ((_minval == undef || vals[_i] < _minval)? vals[_i] : _minval),
109 ((_minval == undef || vals[_i] < _minval)? _i : _minidx),
110 _i+1
111 );
112
113
114// Function: max_index()
115// Usage:
116// max_index(vals);
117// Description:
118// Returns the index of the maximum value in the given list.
119function max_index(vals, _maxval, _maxidx, _i=0) =
120 _i>=len(vals)? _maxidx :
121 max_index(
122 vals,
123 ((_maxval == undef || vals[_i] > _maxval)? vals[_i] : _maxval),
124 ((_maxval == undef || vals[_i] > _maxval)? _i : _maxidx),
125 _i+1
126 );
127
128
129// Function: posmod()
130// Usage:
131// posmod(x,m)
132// Description:
133// Returns the positive modulo `m` of `x`. Value returned will be in the range 0 ... `m`-1.
134// This if useful for normalizing angles to 0 ... 360.
135// Arguments:
136// x = The value to constrain.
137// m = Modulo value.
138function posmod(x,m) = (x%m+m)%m;
139
140
141// Function: modrange()
142// Usage:
143// modrange(x, y, m, [step])
144// Description:
145// Returns a normalized list of values from `x` to `y`, by `step`, modulo `m`. Wraps if `x` > `y`.
146// Arguments:
147// x = The start value to constrain.
148// y = The end value to constrain.
149// m = Modulo value.
150// step = Step by this amount.
151// Examples:
152// modrange(90,270,360, step=45); // Outputs [90,135,180,225,270]
153// modrange(270,90,360, step=45); // Outputs [270,315,0,45,90]
154// modrange(90,270,360, step=-45); // Outputs [90,45,0,315,270]
155// modrange(270,90,360, step=-45); // Outputs [270,225,180,135,90]
156function modrange(x, y, m, step=1) =
157 let(
158 a = posmod(x, m),
159 b = posmod(y, m),
160 c = step>0? (a>b? b+m : b) : (a<b? b-m : b)
161 ) [for (i=[a:step:c]) (i%m+m)%m];
162
163
164// Function: gaussian_rand()
165// Usage:
166// gaussian_rand(mean, stddev)
167// Description:
168// Returns a random number with a gaussian/normal distribution.
169// Arguments:
170// mean = The average random number returned.
171// stddev = The standard deviation of the numbers to be returned.
172function gaussian_rand(mean, stddev) = let(s=rands(0,1,2)) mean + stddev*sqrt(-2*ln(s.x))*cos(360*s.y);
173
174
175// Function: log_rand()
176// Usage:
177// log_rand(minval, maxval, factor);
178// Description:
179// Returns a single random number, with a logarithmic distribution.
180// Arguments:
181// minval = Minimum value to return.
182// maxval = Maximum value to return. `minval` <= X < `maxval`.
183// factor = Log factor to use. Values of X are returned `factor` times more often than X+1.
184function log_rand(minval, maxval, factor) = -ln(1-rands(1-1/pow(factor,minval), 1-1/pow(factor,maxval), 1)[0])/ln(factor);
185
186
187// Function: segs()
188// Description:
189// Calculate the standard number of sides OpenSCAD would give a circle based on `$fn`, `$fa`, and `$fs`.
190// Arguments:
191// r = Radius of circle to get the number of segments for.
192function segs(r) =
193 $fn>0? ($fn>3? $fn : 3) :
194 ceil(max(5, min(360/$fa, abs(r)*2*PI/$fs)));
195
196
197// Function: lerp()
198// Description: Interpolate between two values or vectors.
199// Arguments:
200// a = First value.
201// b = Second value.
202// u = The proportion from `a` to `b` to calculate. Valid range is 0.0 to 1.0, inclusive.
203function lerp(a,b,u) = (1-u)*a + u*b;
204
205
206// Function: hypot()
207// Description: Calculate hypotenuse length of a 2D or 3D triangle.
208// Arguments:
209// x = Length on the X axis.
210// y = Length on the Y axis.
211// z = Length on the Z axis.
212function hypot(x,y,z=0) = norm([x,y,z]);
213
214
215// Function: hypot3()
216// Status: DEPRECATED, use `norm([x,y,z])` instead.
217// Description: Calculate hypotenuse length of 3D triangle.
218// Arguments:
219// x = Length on the X axis.
220// y = Length on the Y axis.
221// z = Length on the Z axis.
222function hypot3(x,y,z) = norm([x,y,z]);
223
224
225// Function: distance()
226// Status: DEPRECATED, use `norm(p2-p1)` instead. It's shorter.
227// Description: Returns the distance between a pair of 2D or 3D points.
228function distance(p1, p2) = norm(point3d(p2)-point3d(p1));
229
230
231// Function: sinh()
232// Description: Takes a value `x`, and returns the hyperbolic sine of it.
233function sinh(x) = (exp(x)-exp(-x))/2;
234
235
236// Function: cosh()
237// Description: Takes a value `x`, and returns the hyperbolic cosine of it.
238function cosh(x) = (exp(x)+exp(-x))/2;
239
240
241// Function: tanh()
242// Description: Takes a value `x`, and returns the hyperbolic tangent of it.
243function tanh(x) = sinh(x)/cosh(x);
244
245
246// Function: asinh()
247// Description: Takes a value `x`, and returns the inverse hyperbolic sine of it.
248function asinh(x) = ln(x+sqrt(x*x+1));
249
250
251// Function: acosh()
252// Description: Takes a value `x`, and returns the inverse hyperbolic cosine of it.
253function acosh(x) = ln(x+sqrt(x*x-1));
254
255
256// Function: atanh()
257// Description: Takes a value `x`, and returns the inverse hyperbolic tangent of it.
258function atanh(x) = ln((1+x)/(1-x))/2;
259
260
261// Function: sum()
262// Description:
263// Returns the sum of all entries in the given array.
264// If passed an array of vectors, returns a vector of sums of each part.
265// Arguments:
266// v = The vector to get the sum of.
267// Example:
268// sum([1,2,3]); // returns 6.
269// sum([[1,2,3], [3,4,5], [5,6,7]]); // returns [9, 12, 15]
270function sum(v, i=0, tot=undef) = i>=len(v)? tot : sum(v, i+1, ((tot==undef)? v[i] : tot+v[i]));
271
272
273// Function: sum_of_squares()
274// Description:
275// Returns the sum of the square of each element of a vector.
276// Arguments:
277// v = The vector to get the sum of.
278// Example:
279// sum_of_squares([1,2,3]); // returns 14.
280function sum_of_squares(v, i=0, tot=0) = sum(vmul(v,v));
281
282
283// Function: sum_of_sines()
284// Usage:
285// sum_of_sines(a,sines)
286// Description:
287// Gives the sum of a series of sines, at a given angle.
288// Arguments:
289// a = Angle to get the value for.
290// sines = List of [amplitude, frequency, offset] items, where the frequency is the number of times the cycle repeats around the circle.
291function sum_of_sines(a, sines) =
292 sum([
293 for (s = sines) let(
294 ss=point3d(s),
295 v=ss.x*sin(a*ss.y+ss.z)
296 ) v
297 ]);
298
299
300// Function: mean()
301// Description:
302// Returns the mean of all entries in the given array.
303// If passed an array of vectors, returns a vector of mean of each part.
304// Arguments:
305// v = The list of values to get the mean of.
306// Example:
307// mean([2,3,4]); // returns 3.
308// mean([[1,2,3], [3,4,5], [5,6,7]]); // returns [3, 4, 5]
309function mean(v) = sum(v)/len(v);
310
311
312// Section: Comparisons and Logic
313
314
315// Function: compare_vals()
316// Usage:
317// compare_vals(a, b);
318// Description:
319// Compares two values. Lists are compared recursively.
320// Results are undefined if the two values are not of similar types.
321// Arguments:
322// a = First value to compare.
323// b = Second value to compare.
324function compare_vals(a, b) =
325 (a==b)? 0 :
326 (a==undef)? -1 :
327 (b==undef)? 1 :
328 ((a==[] || a=="" || a[0]!=undef) && (b==[] || b=="" || b[0]!=undef))? (
329 compare_lists(a, b)
330 ) : (a<b)? -1 :
331 (a>b)? 1 : 0;
332
333
334// Function: compare_lists()
335// Usage:
336// compare_lists(a, b)
337// Description:
338// Compare contents of two lists.
339// Returns <0 if `a`<`b`.
340// Returns 0 if `a`==`b`.
341// Returns >0 if `a`>`b`.
342// Results are undefined if elements are not of similar types.
343// Arguments:
344// a = First list to compare.
345// b = Second list to compare.
346function compare_lists(a, b, n=0) =
347 let(
348 // This curious construction enables tail recursion optimization.
349 cmp = (a==b)? 0 :
350 (len(a)<=n)? -1 :
351 (len(b)<=n)? 1 :
352 (a==a[n] || b==b[n])? (
353 a<b? -1 : a>b? 1 : 0
354 ) : compare_vals(a[n], b[n])
355 )
356 (cmp != 0 || a==b)? cmp :
357 compare_lists(a, b, n+1);
358
359
360// Function: any()
361// Description:
362// Returns true if any item in list `l` evaluates as true.
363// If `l` is a lists of lists, `any()` is applied recursively to each sublist.
364// Arguments:
365// l = The list to test for true items.
366// Example:
367// any([0,false,undef]); // Returns false.
368// any([1,false,undef]); // Returns true.
369// any([1,5,true]); // Returns true.
370// any([[0,0], [0,0]]); // Returns false.
371// any([[0,0], [1,0]]); // Returns true.
372function any(l, i=0, succ=false) =
373 (i>=len(l) || succ)? succ :
374 any(
375 l, i=i+1, succ=(
376 is_array(l[i])? any(l[i]) :
377 !(!l[i])
378 )
379 );
380
381
382// Function: all()
383// Description:
384// Returns true if all items in list `l` evaluate as true.
385// If `l` is a lists of lists, `all()` is applied recursively to each sublist.
386// Arguments:
387// l = The list to test for true items.
388// Example:
389// all([0,false,undef]); // Returns false.
390// all([1,false,undef]); // Returns false.
391// all([1,5,true]); // Returns true.
392// all([[0,0], [0,0]]); // Returns false.
393// all([[0,0], [1,0]]); // Returns false.
394// all([[1,1], [1,1]]); // Returns true.
395function all(l, i=0, fail=false) =
396 (i>=len(l) || fail)? (!fail) :
397 all(
398 l, i=i+1, fail=(
399 is_array(l[i])? !all(l[i]) :
400 !l[i]
401 )
402 );
403
404
405// Function: count_true()
406// Usage:
407// count_true(l)
408// Description:
409// Returns the number of items in `l` that evaluate as true.
410// If `l` is a lists of lists, this is applied recursively to each
411// sublist. Returns the total count of items that evaluate as true
412// in all recursive sublists.
413// Arguments:
414// l = The list to test for true items.
415// nmax = If given, stop counting if `nmax` items evaluate as true.
416// Example:
417// count_true([0,false,undef]); // Returns 0.
418// count_true([1,false,undef]); // Returns 1.
419// count_true([1,5,false]); // Returns 2.
420// count_true([1,5,true]); // Returns 3.
421// count_true([[0,0], [0,0]]); // Returns 0.
422// count_true([[0,0], [1,0]]); // Returns 1.
423// count_true([[1,1], [1,1]]); // Returns 4.
424// count_true([[1,1], [1,1]], nmax=3); // Returns 3.
425function count_true(l, nmax=undef, i=0, cnt=0) =
426 (i>=len(l) || (nmax!=undef && cnt>=nmax))? cnt :
427 count_true(
428 l=l, nmax=nmax, i=i+1, cnt=cnt+(
429 is_array(l[i])? count_true(l[i], nmax=nmax-cnt) :
430 (l[i]? 1 : 0)
431 )
432 );
433
434
435
436// Section: List/Array Operations
437
438
439// Function: cdr()
440// Status: DEPRECATED, use `slice(list,1,-1)` instead.
441// Description: Returns all but the first item of a given array.
442// Arguments:
443// list = The list to get the tail of.
444function cdr(list) = len(list)<=1? [] : [for (i=[1:len(list)-1]) list[i]];
445
446
447// Function: replist()
448// Usage:
449// replist(val, n)
450// Description:
451// Generates a list or array of `n` copies of the given `list`.
452// If the count `n` is given as a list of counts, then this creates a
453// multi-dimensional array, filled with `val`.
454// Arguments:
455// val = The value to repeat to make the list or array.
456// n = The number of copies to make of `val`.
457// Example:
458// replist(1, 4); // Returns [1,1,1,1]
459// replist(8, [2,3]); // Returns [[8,8,8], [8,8,8]]
460// replist(0, [2,2,3]); // Returns [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]
461// replist([1,2,3],3); // Returns [[1,2,3], [1,2,3], [1,2,3]]
462function replist(val, n, i=0) =
463 is_scalar(n)? [for(j=[1:n]) val] :
464 (i>=len(n))? val :
465 [for (j=[1:n[i]]) replist(val, n, i+1)];
466
467
468// Function: in_list()
469// Description: Returns true if value `x` is in list `l`.
470// Arguments:
471// x = The value to search for.
472// l = The list to search.
473// idx = If given, searches the given subindexes for matches for `x`.
474// Example:
475// in_list("bar", ["foo", "bar", "baz"]); // Returns true.
476// in_list("bee", ["foo", "bar", "baz"]); // Returns false.
477// in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1); // Returns true.
478function in_list(x,l,idx=undef) = search([x], l, num_returns_per_match=1, index_col_num=idx) != [[]];
479
480
481// Function: slice()
482// Description:
483// Returns a slice of a list. The first item is index 0.
484// Negative indexes are counted back from the end. The last item is -1.
485// Arguments:
486// arr = The array/list to get the slice of.
487// st = The index of the first item to return.
488// end = The index after the last item to return, unless negative, in which case the last item to return.
489// Example:
490// slice([3,4,5,6,7,8,9], 3, 5); // Returns [6,7]
491// slice([3,4,5,6,7,8,9], 2, -1); // Returns [5,6,7,8,9]
492// slice([3,4,5,6,7,8,9], 1, 1); // Returns []
493// slice([3,4,5,6,7,8,9], 6, -1); // Returns [9]
494// slice([3,4,5,6,7,8,9], 2, -2); // Returns [5,6,7,8]
495function slice(arr,st,end) = let(
496 s=st<0?(len(arr)+st):st,
497 e=end<0?(len(arr)+end+1):end
498 ) (s==e)? [] : [for (i=[s:e-1]) if (e>s) arr[i]];
499
500
501// Function: wrap_range()
502// Status: DEPRECATED, use `select()` instead.
503// Description:
504// Returns a portion of a list, wrapping around past the beginning, if end<start.
505// The first item is index 0. Negative indexes are counted back from the end.
506// The last item is -1. If only the `start` index is given, returns just the value
507// at that position.
508// Usage:
509// wrap_range(list,start)
510// wrap_range(list,start,end)
511// Arguments:
512// list = The list to get the portion of.
513// start = The index of the first item.
514// end = The index of the last item.
515function wrap_range(list, start, end=undef) = select(list,start,end);
516
517
518// Function: select()
519// Description:
520// Returns a portion of a list, wrapping around past the beginning, if end<start.
521// The first item is index 0. Negative indexes are counted back from the end.
522// The last item is -1. If only the `start` index is given, returns just the value
523// at that position.
524// Usage:
525// select(list,start)
526// select(list,start,end)
527// Arguments:
528// list = The list to get the portion of.
529// start = The index of the first item.
530// end = The index of the last item.
531// Example:
532// l = [3,4,5,6,7,8,9];
533// select(l, 5, 6); // Returns [8,9]
534// select(l, 5, 8); // Returns [8,9,3,4]
535// select(l, 5, 2); // Returns [8,9,3,4,5]
536// select(l, -3, -1); // Returns [7,8,9]
537// select(l, 3, 3); // Returns [6]
538// select(l, 4); // Returns 7
539// select(l, -2); // Returns 8
540// select(l, [1:3]); // Returns [4,5,6]
541// select(l, [1,3]); // Returns [4,6]
542function select(list, start, end=undef) =
543 let(l=len(list))
544 (list==[])? [] :
545 end==undef? (
546 is_scalar(start)?
547 let(s=(start%l+l)%l) list[s] :
548 [for (i=start) list[(i%l+l)%l]]
549 ) : (
550 let(s=(start%l+l)%l, e=(end%l+l)%l)
551 (s<=e)?
552 [for (i = [s:e]) list[i]] :
553 concat([for (i = [s:l-1]) list[i]], [for (i = [0:e]) list[i]])
554 );
555
556
557// Function: reverse()
558// Description: Reverses a list/array.
559// Arguments:
560// list = The list to reverse.
561// Example:
562// reverse([3,4,5,6]); // Returns [6,5,4,3]
563function reverse(list) = [ for (i = [len(list)-1 : -1 : 0]) list[i] ];
564
565
566// Function: array_subindex()
567// Description:
568// For each array item, return the indexed subitem.
569// Returns a list of the values of each vector at the specfied
570// index list or range. If the index list or range has
571// only one entry the output list is flattened.
572// Arguments:
573// v = The given list of lists.
574// idx = The index, list of indices, or range of indices to fetch.
575// Example:
576// v = [[[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
577// array_subindex(v,2); // Returns [3, 7, 11, 15]
578// array_subindex(v,[2,1]); // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
579// array_subindex(v,[1:3]); // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
580function array_subindex(v, idx) = [
581 for(val=v) let(value=[for(i=idx) val[i]])
582 len(value)==1 ? value[0] : value
583];
584
585
586// Function: list_range()
587// Usage:
588// list_range(n, [s], [e], [step])
589// list_range(e, [step])
590// list_range(s, e, [step])
591// Description:
592// Returns a list, counting up from starting value `s`, by `step` increments,
593// until either `n` values are in the list, or it reaches the end value `e`.
594// Arguments:
595// n = Desired number of values in returned list, if given.
596// s = Starting value. Default: 0
597// e = Ending value to stop at, if given.
598// step = Amount to increment each value. Default: 1
599// Example:
600// list_range(4); // Returns [0,1,2,3]
601// list_range(n=4, step=2); // Returns [0,2,4,6]
602// list_range(n=4, s=3, step=3); // Returns [3,6,9,12]
603// list_range(n=5, s=0, e=10); // Returns [0, 2.5, 5, 7.5, 10]
604// list_range(e=3); // Returns [0,1,2,3]
605// list_range(e=6, step=2); // Returns [0,2,4,6]
606// list_range(s=3, e=5); // Returns [3,4,5]
607// list_range(s=3, e=8, step=2); // Returns [3,5,7]
608// list_range(s=4, e=8, step=2); // Returns [4,6,8]
609// list_range(n=4, s=[3,4], step=[2,3]); // Returns [[3,4], [5,7], [7,10], [9,13]]
610function list_range(n=undef, s=0, e=undef, step=1) =
611 (n!=undef && e!=undef)? [for (i=[0:1:n-1]) s+(e-s)*i/(n-1)] :
612 (n!=undef)? [for (i=[0:n-1]) let(v=s+step*i) v] :
613 (e!=undef)? [for (v=[s:step:e]) v] :
614 assertion(false, "Must supply one of `n` or `e`.");
615
616
617// Function: array_shortest()
618// Description:
619// Returns the length of the shortest sublist in a list of lists.
620// Arguments:
621// vecs = A list of lists.
622function array_shortest(vecs) = min([for (v = vecs) len(v)]);
623
624
625// Function: array_longest()
626// Description:
627// Returns the length of the longest sublist in a list of lists.
628// Arguments:
629// vecs = A list of lists.
630function array_longest(vecs) = max([for (v = vecs) len(v)]);
631
632
633// Function: array_pad()
634// Description:
635// If the list `v` is shorter than `minlen` length, pad it to length with the value given in `fill`.
636// Arguments:
637// v = A list.
638// minlen = The minimum length to pad the list to.
639// fill = The value to pad the list with.
640function array_pad(v, minlen, fill=undef) = let(l=len(v)) [for (i=[0:max(l,minlen)-1]) i<l? v[i] : fill];
641
642
643// Function: array_trim()
644// Description:
645// If the list `v` is longer than `maxlen` length, truncates it to be `maxlen` items long.
646// Arguments:
647// v = A list.
648// minlen = The minimum length to pad the list to.
649function array_trim(v, maxlen) = maxlen<1? [] : [for (i=[0:min(len(v),maxlen)-1]) v[i]];
650
651
652// Function: array_fit()
653// Description:
654// If the list `v` is longer than `length` items long, truncates it to be exactly `length` items long.
655// If the list `v` is shorter than `length` items long, pad it to length with the value given in `fill`.
656// Arguments:
657// v = A list.
658// minlen = The minimum length to pad the list to.
659// fill = The value to pad the list with.
660function array_fit(v, length, fill) = let(l=len(v)) (l==length)? v : (l>length)? array_trim(v,length) : array_pad(v,length,fill);
661
662
663// Function: enumerate()
664// Description:
665// Returns a list, with each item of the given list `l` numbered in a sublist.
666// Something like: `[[0,l[0]], [1,l[1]], [2,l[2]], ...]`
667// Arguments:
668// l = List to enumerate.
669// idx = If given, enumerates just the given subindex items of `l`.
670// Example:
671// enumerate(["a","b","c"]); // Returns: [[0,"a"], [1,"b"], [2,"c"]]
672// enumerate([[88,"a"],[76,"b"],[21,"c"]], idx=1); // Returns: [[0,"a"], [1,"b"], [2,"c"]]
673// enumerate([["cat","a",12],["dog","b",10],["log","c",14]], idx=[1:2]); // Returns: [[0,"a",12], [1,"b",10], [2,"c",14]]
674function enumerate(l,idx=undef) =
675 (l==[])? [] :
676 (idx==undef)?
677 [for (i=[0:len(l)-1]) [i,l[i]]] :
678 [for (i=[0:len(l)-1]) concat([i], [for (j=idx) l[i][j]])];
679
680
681// Function: array_zip()
682// Usage:
683// array_zip(v1, v2, v3, [fit], [fill]);
684// array_zip(vecs, [fit], [fill]);
685// Description:
686// Zips together corresponding items from two or more lists.
687// Returns a list of lists, where each sublist contains corresponding
688// items from each of the input lists. `[[A1, B1, C1], [A2, B2, C2], ...]`
689// Arguments:
690// vecs = A list of two or more lists to zipper together.
691// fit = If `fit=="short"`, the zips together up to the length of the shortest list in vecs. If `fit=="long"`, then pads all lists to the length of the longest, using the value in `fill`. If `fit==false`, then requires all lists to be the same length. Default: false.
692// fill = The default value to fill in with if one or more lists if short. Default: undef
693// Example:
694// v1 = [1,2,3,4];
695// v2 = [5,6,7];
696// v3 = [8,9,10,11];
697// array_zip(v1,v3); // returns [[1,8], [2,9], [3,10], [4,11]]
698// array_zip([v1,v3]); // returns [[1,8], [2,9], [3,10], [4,11]]
699// array_zip([v1,v2], fit="short"); // returns [[1,5], [2,6], [3,7]]
700// array_zip([v1,v2], fit="long"); // returns [[1,5], [2,6], [3,7], [4,undef]]
701// array_zip([v1,v2], fit="long, fill=0); // returns [[1,5], [2,6], [3,7], [4,0]]
702// array_zip([v1,v2,v3], fit="long"); // returns [[1,5,8], [2,6,9], [3,7,10], [4,undef,11]]
703// Example:
704// v1 = [[1,2,3], [4,5,6], [7,8,9]];
705// v2 = [[20,19,18], [17,16,15], [14,13,12]];
706// array_zip(v1,v2); // Returns [[1,2,3,20,19,18], [4,5,6,17,16,15], [7,8,9,14,13,12]]
707function array_zip(vecs, v2, v3, fit=false, fill=undef) =
708 (v3!=undef)? array_zip([vecs,v2,v3], fit=fit, fill=fill) :
709 (v2!=undef)? array_zip([vecs,v2], fit=fit, fill=fill) :
710 let(
711 dummy1 = assert_in_list("fit", fit, [false, "short", "long"]),
712 minlen = array_shortest(vecs),
713 maxlen = array_longest(vecs),
714 dummy2 = (fit==false)? assertion(minlen==maxlen, "Input vectors must have the same length") : 0
715 ) (fit == "long")?
716 [for(i=[0:maxlen-1]) [for(v=vecs) for(x=(i<len(v)? v[i] : (fill==undef)? [fill] : fill)) x] ] :
717 [for(i=[0:minlen-1]) [for(v=vecs) for(x=v[i]) x] ];
718
719
720
721// Function: array_group()
722// Description:
723// Takes a flat array of values, and groups items in sets of `cnt` length.
724// The opposite of this is `flatten()`.
725// Arguments:
726// v = The list of items to group.
727// cnt = The number of items to put in each grouping.
728// dflt = The default value to fill in with is the list is not a multiple of `cnt` items long.
729// Example:
730// v = [1,2,3,4,5,6];
731// array_group(v,2) returns [[1,2], [3,4], [5,6]]
732// array_group(v,3) returns [[1,2,3], [4,5,6]]
733// array_group(v,4,0) returns [[1,2,3,4], [5,6,0,0]]
734function array_group(v, cnt=2, dflt=0) = [for (i = [0:cnt:len(v)-1]) [for (j = [0:cnt-1]) default(v[i+j], dflt)]];
735
736
737// Function: flatten()
738// Description: Takes a list of lists and flattens it by one level.
739// Arguments:
740// l = List to flatten.
741// Example:
742// flatten([[1,2,3], [4,5,[6,7,8]]]) returns [1,2,3,4,5,[6,7,8]]
743function flatten(l) = [for (a = l) for (b = a) b];
744
745
746// Function: sort()
747// Usage:
748// sort(arr, [idx])
749// Description:
750// Sorts the given list using `compare_vals()`. Results are undefined if list elements are not of similar type.
751// Arguments:
752// arr = The list to sort.
753// idx = If given, the index, range, or list of indices of sublist items to compare.
754// Example:
755// l = [45,2,16,37,8,3,9,23,89,12,34];
756// sorted = sort(l); // Returns [2,3,8,9,12,16,23,34,37,45,89]
757function sort(arr, idx=undef) =
758 (len(arr)<=1) ? arr :
759 let(
760 pivot = arr[floor(len(arr)/2)],
761 pivotval = idx==undef? pivot : [for (i=idx) pivot[i]],
762 compare = [
763 for (entry = arr) let(
764 val = idx==undef? entry : [for (i=idx) entry[i]],
765 cmp = compare_vals(val, pivotval)
766 ) cmp
767 ],
768 lesser = [ for (i = [0:len(arr)-1]) if (compare[i] < 0) arr[i] ],
769 equal = [ for (i = [0:len(arr)-1]) if (compare[i] ==0) arr[i] ],
770 greater = [ for (i = [0:len(arr)-1]) if (compare[i] > 0) arr[i] ]
771 )
772 concat(sort(lesser,idx), equal, sort(greater,idx));
773
774
775// Function: sortidx()
776// Description:
777// Given a list, calculates the sort order of the list, and returns
778// a list of indexes into the original list in that sorted order.
779// If you iterate the returned list in order, and use the list items
780// to index into the original list, you will be iterating the original
781// values in sorted order.
782// Example:
783// lst = ["d","b","e","c"];
784// idxs = sortidx(lst); // Returns: [1,3,0,2]
785// ordered = [for (i=idxs) lst[i]]; // Returns: ["b", "c", "d", "e"]
786// Example:
787// lst = [
788// ["foo", 88, [0,0,1], false],
789// ["bar", 90, [0,1,0], true],
790// ["baz", 89, [1,0,0], false],
791// ["qux", 23, [1,1,1], true]
792// ];
793// idxs1 = sortidx(lst, idx=1); // Returns: [3,0,2,1]
794// idxs2 = sortidx(lst, idx=0); // Returns: [1,2,0,3]
795// idxs3 = sortidx(lst, idx=[1,3]); // Returns: [3,0,2,1]
796function sortidx(l, idx=undef) =
797 (l==[])? [] :
798 let(
799 ll=enumerate(l,idx=idx),
800 sidx = [1:len(ll[0])-1]
801 )
802 array_subindex(sort(ll, idx=sidx), 0);
803
804
805// Function: unique()
806// Usage:
807// unique(arr);
808// Description:
809// Returns a sorted list with all repeated items removed.
810// Arguments:
811// arr = The list to uniquify.
812function unique(arr) =
813 len(arr)<=1? arr : let(
814 sorted = sort(arr)
815 ) [
816 for (i=[0:len(sorted)-1])
817 if (i==0 || (sorted[i] != sorted[i-1]))
818 sorted[i]
819 ];
820
821
822
823// Function: list_remove()
824// Usage:
825// list_remove(list, elements)
826// Description:
827// Remove all items from `list` whose indexes are in `elements`.
828// Arguments:
829// list = The list to remove items from.
830// elements = The list of indexes of items to remove.
831function list_remove(list, elements) = [
832 for (i = [0:len(list)-1]) if (!search(i, elements)) list[i]
833];
834
835
836
837// Internal. Not exposed.
838function _array_dim_recurse(v) =
839 !is_list(v[0])? (
840 sum( [for(entry=v) is_list(entry) ? 1 : 0]) == 0 ? [] : [undef]
841 ) : let(
842 firstlen = len(v[0]),
843 first = sum( [for(entry = v) len(entry) == firstlen ? 0 : 1] ) == 0 ? firstlen : undef,
844 leveldown = flatten(v)
845 ) is_list(leveldown[0])? (
846 concat([first],_array_dim_recurse(leveldown))
847 ) : [first];
848
849
850// Function: array_dim()
851// Usage:
852// array_dim(v, [depth])
853// Description:
854// Returns the size of a multi-dimensional array. Returns a list of
855// dimension lengths. The length of `v` is the dimension `0`. The
856// length of the items in `v` is dimension `1`. The length of the
857// items in the items in `v` is dimension `2`, etc. For each dimension,
858// if the length of items at that depth is inconsistent, `undef` will
859// be returned. If no items of that dimension depth exist, `0` is
860// returned. Otherwise, the consistent length of items in that
861// dimensional depth is returned.
862// Arguments:
863// v = Array to get dimensions of.
864// depth = Dimension to get size of. If not given, returns a list of dimension lengths.
865// Examples:
866// array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]]); // Returns [2,2,3]
867// array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0); // Returns 2
868// array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2); // Returns 3
869// array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]); // Returns [2,undef,3]
870function array_dim(v, depth=undef) =
871 (depth == undef)? (
872 concat([len(v)], _array_dim_recurse(v))
873 ) : (depth == 0)? (
874 len(v)
875 ) : (
876 let(dimlist = _array_dim_recurse(v))
877 (depth > len(dimlist))? 0 : dimlist[depth-1]
878 );
879
880
881
882// Section: Vector Manipulation
883
884// Function: vmul()
885// Description:
886// Element-wise vector multiplication. Multiplies each element of vector `v1` by
887// the corresponding element of vector `v2`. Returns a vector of the products.
888// Arguments:
889// v1 = The first vector.
890// v2 = The second vector.
891// Example:
892// vmul([3,4,5], [8,7,6]); // Returns [24, 28, 30]
893function vmul(v1, v2) = [for (i = [0:len(v1)-1]) v1[i]*v2[i]];
894
895
896// Function: vdiv()
897// Description:
898// Element-wise vector division. Divides each element of vector `v1` by
899// the corresponding element of vector `v2`. Returns a vector of the quotients.
900// Arguments:
901// v1 = The first vector.
902// v2 = The second vector.
903// Example:
904// vdiv([24,28,30], [8,7,6]); // Returns [3, 4, 5]
905function vdiv(v1, v2) = [for (i = [0:len(v1)-1]) v1[i]/v2[i]];
906
907
908// Function: vabs()
909// Description: Returns a vector of the absolute value of each element of vector `v`.
910// Arguments:
911// v = The vector to get the absolute values of.
912function vabs(v) = [for (x=v) abs(x)];
913
914
915// Function: normalize()
916// Description:
917// Returns unit length normalized version of vector v.
918// Arguments:
919// v = The vector to normalize.
920function normalize(v) = v/norm(v);
921
922
923// Function: vector2d_angle()
924// Status: DEPRECATED, use `vector_angle()` instead.
925// Usage:
926// vector2d_angle(v1,v2);
927// Description:
928// Returns angle in degrees between two 2D vectors.
929// Arguments:
930// v1 = First 2D vector.
931// v2 = Second 2D vector.
932function vector2d_angle(v1,v2) = vector_angle(v1,v2);
933
934
935// Function: vector3d_angle()
936// Status: DEPRECATED, use `vector_angle()` instead.
937// Usage:
938// vector3d_angle(v1,v2);
939// Description:
940// Returns angle in degrees between two 3D vectors.
941// Arguments:
942// v1 = First 3D vector.
943// v2 = Second 3D vector.
944function vector3d_angle(v1,v2) = vector_angle(v1,v2);
945
946
947// Function: vector_angle()
948// Usage:
949// vector_angle(v1,v2);
950// Description:
951// Returns angle in degrees between two vectors of similar dimensions.
952// Arguments:
953// v1 = First vector.
954// v2 = Second vector.
955// NOTE: constrain() corrects crazy FP rounding errors that exceed acos()'s domain.
956function vector_angle(v1,v2) = acos(constrain((v1*v2)/(norm(v1)*norm(v2)), -1, 1));
957
958
959// Function: vector_axis()
960// Usage:
961// vector_xis(v1,v2);
962// Description:
963// Returns the vector perpendicular to both of the given vectors.
964// Arguments:
965// v1 = First vector.
966// v2 = Second vector.
967function vector_axis(v1,v2) =
968 let(
969 eps = 1e-6,
970 v1 = point3d(v1/norm(v1)),
971 v2 = point3d(v2/norm(v2)),
972 v3 = (norm(v1-v2) > eps && norm(v1+v2) > eps)? v2 :
973 (norm(vabs(v2)-V_UP) > eps)? V_UP :
974 V_RIGHT
975 ) normalize(cross(v1,v3));
976
977
978// Section: Coordinates Manipulation
979
980// Function: point2d()
981// Description:
982// Returns a 2D vector/point from a 2D or 3D vector.
983// If given a 3D point, removes the Z coordinate.
984// Arguments:
985// p = The coordinates to force into a 2D vector/point.
986function point2d(p) = [for (i=[0:1]) (p[i]==undef)? 0 : p[i]];
987
988
989// Function: path2d()
990// Description:
991// Returns a list of 2D vectors/points from a list of 2D or 3D vectors/points.
992// If given a 3D point list, removes the Z coordinates from each point.
993// Arguments:
994// points = A list of 2D or 3D points/vectors.
995function path2d(points) = [for (point = points) point2d(point)];
996
997
998// Function: point3d()
999// Description:
1000// Returns a 3D vector/point from a 2D or 3D vector.
1001// Arguments:
1002// p = The coordinates to force into a 3D vector/point.
1003function point3d(p) = [for (i=[0:2]) (p[i]==undef)? 0 : p[i]];
1004
1005
1006// Function: path3d()
1007// Description:
1008// Returns a list of 3D vectors/points from a list of 2D or 3D vectors/points.
1009// Arguments:
1010// points = A list of 2D or 3D points/vectors.
1011function path3d(points) = [for (point = points) point3d(point)];
1012
1013
1014// Function: translate_points()
1015// Usage:
1016// translate_points(pts, v);
1017// Description:
1018// Moves each point in an array by a given amount.
1019// Arguments:
1020// pts = List of points to translate.
1021// v = Amount to translate points by.
1022function translate_points(pts, v=[0,0,0]) = [for (pt = pts) pt+v];
1023
1024
1025// Function: scale_points()
1026// Usage:
1027// scale_points(pts, v, [cp]);
1028// Description:
1029// Scales each point in an array by a given amount, around a given centerpoint.
1030// Arguments:
1031// pts = List of points to scale.
1032// v = A vector with a scaling factor for each axis.
1033// cp = Centerpoint to scale around.
1034function scale_points(pts, v=[0,0,0], cp=[0,0,0]) = [for (pt = pts) [for (i = [0:len(pt)-1]) (pt[i]-cp[i])*v[i]+cp[i]]];
1035
1036
1037// Function: rotate_points2d()
1038// Usage:
1039// rotate_points2d(pts, ang, [cp]);
1040// Description:
1041// Rotates each 2D point in an array by a given amount, around an optional centerpoint.
1042// Arguments:
1043// pts = List of 3D points to rotate.
1044// ang = Angle to rotate by.
1045// cp = 2D Centerpoint to rotate around. Default: `[0,0]`
1046function rotate_points2d(pts, ang, cp=[0,0]) = let(
1047 m = matrix3_zrot(ang)
1048 ) [for (pt = pts) m*point3d(pt-cp)+cp];
1049
1050
1051// Function: rotate_points3d()
1052// Usage:
1053// rotate_points3d(pts, v, [cp], [reverse]);
1054// rotate_points3d(pts, v, axis, [cp], [reverse]);
1055// rotate_points3d(pts, from, to, v, [cp], [reverse]);
1056// Description:
1057// Rotates each 3D point in an array by a given amount, around a given centerpoint.
1058// Arguments:
1059// pts = List of points to rotate.
1060// v = Rotation angle(s) in degrees.
1061// axis = If given, axis vector to rotate around.
1062// cp = Centerpoint to rotate around.
1063// from = If given, the vector to rotate something from. Used with `to`.
1064// to = If given, the vector to rotate something to. Used with `from`.
1065// reverse = If true, performs an exactly reversed rotation.
1066function rotate_points3d(pts, v=0, cp=[0,0,0], axis=undef, from=undef, to=undef, reverse=false) =
1067 let(
1068 dummy = assertion(is_def(from)==is_def(to), "`from` and `to` must be given together."),
1069 mrot = reverse? (
1070 is_def(from)? (
1071 let (
1072 from = from / norm(from),
1073 to = to / norm(from),
1074 ang = vector_angle(from, to),
1075 axis = vector_axis(from, to)
1076 )
1077 matrix4_rot_by_axis(from, -v) * matrix4_rot_by_axis(axis, -ang)
1078 ) : is_def(axis)? (
1079 matrix4_rot_by_axis(axis, -v)
1080 ) : is_scalar(v)? (
1081 matrix4_zrot(-v)
1082 ) : (
1083 matrix4_xrot(-v.x) * matrix4_yrot(-v.y) * matrix4_zrot(-v.z)
1084 )
1085 ) : (
1086 is_def(from)? (
1087 let (
1088 from = from / norm(from),
1089 to = to / norm(from),
1090 ang = vector_angle(from, to),
1091 axis = vector_axis(from, to)
1092 )
1093 matrix4_rot_by_axis(axis, ang) * matrix4_rot_by_axis(from, v)
1094 ) : is_def(axis)? (
1095 matrix4_rot_by_axis(axis, v)
1096 ) : is_scalar(v)? (
1097 matrix4_zrot(v)
1098 ) : (
1099 matrix4_zrot(v.z) * matrix4_yrot(v.y) * matrix4_xrot(v.x)
1100 )
1101 ),
1102 m = matrix4_translate(cp) * mrot * matrix4_translate(-cp)
1103 ) [for (pt = pts) point3d(m*concat(point3d(pt),[1]))];
1104
1105
1106
1107// Function: rotate_points3d_around_axis()
1108// Status: DEPRECATED, use `rotate_points3d(pts, v=ang, axis=u, cp=cp)` instead.
1109// Usage:
1110// rotate_points3d_around_axis(pts, ang, u, [cp])
1111// Description:
1112// Rotates each 3D point in an array by a given amount, around a given centerpoint and axis.
1113// Arguments:
1114// pts = List of 3D points to rotate.
1115// ang = Angle to rotate by.
1116// u = Vector of the axis to rotate around.
1117// cp = 3D Centerpoint to rotate around.
1118function rotate_points3d_around_axis(pts, ang, u=[0,0,0], cp=[0,0,0]) = let(
1119 m = matrix4_rot_by_axis(u, ang)
1120 ) [for (pt = pts) m*concat(point3d(pt)-cp, 0)+cp];
1121
1122
1123// Section: Coordinate Systems
1124
1125// Function: polar_to_xy()
1126// Usage:
1127// polar_to_xy(r, theta);
1128// polar_to_xy([r, theta]);
1129// Description:
1130// Convert polar coordinates to 2D cartesian coordinates.
1131// Returns [X,Y] cartesian coordinates.
1132// Arguments:
1133// r = distance from the origin.
1134// theta = angle in degrees, counter-clockwise of X+.
1135// Examples:
1136// xy = polar_to_xy(20,30);
1137// xy = polar_to_xy([40,60]);
1138function polar_to_xy(r,theta=undef) = let(
1139 rad = theta==undef? r[0] : r,
1140 t = theta==undef? r[1] : theta
1141 ) rad*[cos(t), sin(t)];
1142
1143
1144// Function: xy_to_polar()
1145// Usage:
1146// xy_to_polar(x,y);
1147// xy_to_polar([X,Y]);
1148// Description:
1149// Convert 2D cartesian coordinates to polar coordinates.
1150// Returns [radius, theta] where theta is the angle counter-clockwise of X+.
1151// Arguments:
1152// x = X coordinate.
1153// y = Y coordinate.
1154// Examples:
1155// plr = xy_to_polar(20,30);
1156// plr = xy_to_polar([40,60]);
1157function xy_to_polar(x,y=undef) = let(
1158 xx = y==undef? x[0] : x,
1159 yy = y==undef? x[1] : y
1160 ) [norm([xx,yy]), atan2(yy,xx)];
1161
1162
1163// Function: xyz_to_planar()
1164// Usage:
1165// xyz_to_planar(point, a, b, c);
1166// Description:
1167// Given three points defining a plane, returns the projected planar
1168// [X,Y] coordinates of the closest point to a 3D `point`. The origin
1169// of the planar coordinate system [0,0] will be at point `a`, and the
1170// Y+ axis direction will be towards point `b`. This coordinate system
1171// can be useful in taking a set of nearly coplanar points, and converting
1172// them to a pure XY set of coordinates for manipulation, before convering
1173// them back to the original 3D plane.
1174function xyz_to_planar(point, a, b, c) = let(
1175 u = normalize(b-a),
1176 v = normalize(c-a),
1177 n = normalize(cross(u,v)),
1178 w = normalize(cross(n,u)),
1179 relpoint = point-a
1180) [relpoint * w, relpoint * u];
1181
1182
1183// Function: planar_to_xyz()
1184// Usage:
1185// planar_to_xyz(point, a, b, c);
1186// Description:
1187// Given three points defining a plane, converts a planar [X,Y]
1188// coordinate to the actual corresponding 3D point on the plane.
1189// The origin of the planar coordinate system [0,0] will be at point
1190// `a`, and the Y+ axis direction will be towards point `b`.
1191function planar_to_xyz(point, a, b, c) = let(
1192 u = normalize(b-a),
1193 v = normalize(c-a),
1194 n = normalize(cross(u,v)),
1195 w = normalize(cross(n,u))
1196) a + point.x * w + point.y * u;
1197
1198
1199// Function: cylindrical_to_xyz()
1200// Usage:
1201// cylindrical_to_xyz(r, theta, z)
1202// cylindrical_to_xyz([r, theta, z])
1203// Description:
1204// Convert cylindrical coordinates to 3D cartesian coordinates. Returns [X,Y,Z] cartesian coordinates.
1205// Arguments:
1206// r = distance from the Z axis.
1207// theta = angle in degrees, counter-clockwise of X+ on the XY plane.
1208// z = Height above XY plane.
1209// Examples:
1210// xyz = cylindrical_to_xyz(20,30,40);
1211// xyz = cylindrical_to_xyz([40,60,50]);
1212function cylindrical_to_xyz(r,theta=undef,z=undef) = let(
1213 rad = theta==undef? r[0] : r,
1214 t = theta==undef? r[1] : theta,
1215 zed = theta==undef? r[2] : z
1216 ) [rad*cos(t), rad*sin(t), zed];
1217
1218
1219// Function: xyz_to_cylindrical()
1220// Usage:
1221// xyz_to_cylindrical(x,y,z)
1222// xyz_to_cylindrical([X,Y,Z])
1223// Description:
1224// Convert 3D cartesian coordinates to cylindrical coordinates.
1225// Returns [radius,theta,Z]. Theta is the angle counter-clockwise
1226// of X+ on the XY plane. Z is height above the XY plane.
1227// Arguments:
1228// x = X coordinate.
1229// y = Y coordinate.
1230// z = Z coordinate.
1231// Examples:
1232// cyl = xyz_to_cylindrical(20,30,40);
1233// cyl = xyz_to_cylindrical([40,50,70]);
1234function xyz_to_cylindrical(x,y=undef,z=undef) = let(
1235 p = is_scalar(x)? [x, default(y,0), default(z,0)] : point3d(x)
1236 ) [norm([p.x,p.y]), atan2(p.y,p.x), p.z];
1237
1238
1239// Function: spherical_to_xyz()
1240// Usage:
1241// spherical_to_xyz(r, theta, phi);
1242// spherical_to_xyz([r, theta, phi]);
1243// Description:
1244// Convert spherical coordinates to 3D cartesian coordinates.
1245// Returns [X,Y,Z] cartesian coordinates.
1246// Arguments:
1247// r = distance from origin.
1248// theta = angle in degrees, counter-clockwise of X+ on the XY plane.
1249// phi = angle in degrees from the vertical Z+ axis.
1250// Examples:
1251// xyz = spherical_to_xyz(20,30,40);
1252// xyz = spherical_to_xyz([40,60,50]);
1253function spherical_to_xyz(r,theta=undef,phi=undef) = let(
1254 rad = theta==undef? r[0] : r,
1255 t = theta==undef? r[1] : theta,
1256 p = theta==undef? r[2] : phi
1257 ) rad*[sin(p)*cos(t), sin(p)*sin(t), cos(p)];
1258
1259
1260// Function: xyz_to_spherical()
1261// Usage:
1262// xyz_to_spherical(x,y,z)
1263// xyz_to_spherical([X,Y,Z])
1264// Description:
1265// Convert 3D cartesian coordinates to spherical coordinates.
1266// Returns [r,theta,phi], where phi is the angle from the Z+ pole,
1267// and theta is degrees counter-clockwise of X+ on the XY plane.
1268// Arguments:
1269// x = X coordinate.
1270// y = Y coordinate.
1271// z = Z coordinate.
1272// Examples:
1273// sph = xyz_to_spherical(20,30,40);
1274// sph = xyz_to_spherical([40,50,70]);
1275function xyz_to_spherical(x,y=undef,z=undef) = let(
1276 p = is_scalar(x)? [x, default(y,0), default(z,0)] : point3d(x)
1277 ) [norm(p), atan2(p.y,p.x), atan2(norm([p.x,p.y]),p.z)];
1278
1279
1280// Function: altaz_to_xyz()
1281// Usage:
1282// altaz_to_xyz(alt, az, r);
1283// altaz_to_xyz([alt, az, r]);
1284// Description:
1285// Convert altitude/azimuth/range coordinates to 3D cartesian coordinates.
1286// Returns [X,Y,Z] cartesian coordinates.
1287// Arguments:
1288// alt = altitude angle in degrees above the XY plane.
1289// az = azimuth angle in degrees clockwise of Y+ on the XY plane.
1290// r = distance from origin.
1291// Examples:
1292// xyz = altaz_to_xyz(20,30,40);
1293// xyz = altaz_to_xyz([40,60,50]);
1294function altaz_to_xyz(alt,az=undef,r=undef) = let(
1295 p = az==undef? alt[0] : alt,
1296 t = 90 - (az==undef? alt[1] : az),
1297 rad = az==undef? alt[2] : r
1298 ) rad*[cos(p)*cos(t), cos(p)*sin(t), sin(p)];
1299
1300
1301// Function: xyz_to_altaz()
1302// Usage:
1303// xyz_to_altaz(x,y,z);
1304// xyz_to_altaz([X,Y,Z]);
1305// Description:
1306// Convert 3D cartesian coordinates to altitude/azimuth/range coordinates.
1307// Returns [altitude,azimuth,range], where altitude is angle above the
1308// XY plane, azimuth is degrees clockwise of Y+ on the XY plane, and
1309// range is the distance from the origin.
1310// Arguments:
1311// x = X coordinate.
1312// y = Y coordinate.
1313// z = Z coordinate.
1314// Examples:
1315// aa = xyz_to_altaz(20,30,40);
1316// aa = xyz_to_altaz([40,50,70]);
1317function xyz_to_altaz(x,y=undef,z=undef) = let(
1318 p = is_scalar(x)? [x, default(y,0), default(z,0)] : point3d(x)
1319 ) [atan2(p.z,norm([p.x,p.y])), atan2(p.x,p.y), norm(p)];
1320
1321
1322// Section: Matrix Manipulation
1323
1324// Function: ident()
1325// Description: Create an `n` by `n` identity matrix.
1326// Arguments:
1327// n = The size of the identity matrix square, `n` by `n`.
1328function ident(n) = [for (i = [0:n-1]) [for (j = [0:n-1]) (i==j)?1:0]];
1329
1330
1331
1332// Function: matrix_transpose()
1333// Description: Returns the transposition of the given matrix.
1334// Example:
1335// m = [
1336// [11,12,13,14],
1337// [21,22,23,24],
1338// [31,32,33,34],
1339// [41,42,43,44]
1340// ];
1341// tm = matrix_transpose(m);
1342// // Returns:
1343// // [
1344// // [11,21,31,41],
1345// // [12,22,32,42],
1346// // [13,23,33,43],
1347// // [14,24,34,44]
1348// // ]
1349function matrix_transpose(m) = [for (i=[0:len(m[0])-1]) [for (j=[0:len(m)-1]) m[j][i]]];
1350
1351
1352
1353// Function: mat3_to_mat4()
1354// Description: Takes a 3x3 matrix and returns its 4x4 equivalent.
1355function mat3_to_mat4(m) = concat(
1356 [for (r = [0:2])
1357 concat(
1358 [for (c = [0:2]) m[r][c]],
1359 [0]
1360 )
1361 ],
1362 [[0, 0, 0, 1]]
1363);
1364
1365
1366
1367// Function: matrix3_translate()
1368// Description:
1369// Returns the 3x3 matrix to perform a 2D translation.
1370// Arguments:
1371// v = 2D Offset to translate by. [X,Y]
1372function matrix3_translate(v) = [
1373 [1, 0, v.x],
1374 [0, 1, v.y],
1375 [0 ,0, 1]
1376];
1377
1378
1379// Function: matrix4_translate()
1380// Description:
1381// Returns the 4x4 matrix to perform a 3D translation.
1382// Arguments:
1383// v = 3D offset to translate by. [X,Y,Z]
1384function matrix4_translate(v) = [
1385 [1, 0, 0, v.x],
1386 [0, 1, 0, v.y],
1387 [0, 0, 1, v.z],
1388 [0 ,0, 0, 1]
1389];
1390
1391
1392// Function: matrix3_scale()
1393// Description:
1394// Returns the 3x3 matrix to perform a 2D scaling transformation.
1395// Arguments:
1396// v = 2D vector of scaling factors. [X,Y]
1397function matrix3_scale(v) = [
1398 [v.x, 0, 0],
1399 [ 0, v.y, 0],
1400 [ 0, 0, 1]
1401];
1402
1403
1404// Function: matrix4_scale()
1405// Description:
1406// Returns the 4x4 matrix to perform a 3D scaling transformation.
1407// Arguments:
1408// v = 3D vector of scaling factors. [X,Y,Z]
1409function matrix4_scale(v) = [
1410 [v.x, 0, 0, 0],
1411 [ 0, v.y, 0, 0],
1412 [ 0, 0, v.z, 0],
1413 [ 0, 0, 0, 1]
1414];
1415
1416
1417// Function: matrix3_zrot()
1418// Description:
1419// Returns the 3x3 matrix to perform a rotation of a 2D vector around the Z axis.
1420// Arguments:
1421// ang = Number of degrees to rotate.
1422function matrix3_zrot(ang) = [
1423 [cos(ang), -sin(ang), 0],
1424 [sin(ang), cos(ang), 0],
1425 [ 0, 0, 1]
1426];
1427
1428
1429// Function: matrix4_xrot()
1430// Description:
1431// Returns the 4x4 matrix to perform a rotation of a 3D vector around the X axis.
1432// Arguments:
1433// ang = number of degrees to rotate.
1434function matrix4_xrot(ang) = [
1435 [1, 0, 0, 0],
1436 [0, cos(ang), -sin(ang), 0],
1437 [0, sin(ang), cos(ang), 0],
1438 [0, 0, 0, 1]
1439];
1440
1441
1442// Function: matrix4_yrot()
1443// Description:
1444// Returns the 4x4 matrix to perform a rotation of a 3D vector around the Y axis.
1445// Arguments:
1446// ang = Number of degrees to rotate.
1447function matrix4_yrot(ang) = [
1448 [ cos(ang), 0, sin(ang), 0],
1449 [ 0, 1, 0, 0],
1450 [-sin(ang), 0, cos(ang), 0],
1451 [ 0, 0, 0, 1]
1452];
1453
1454
1455// Function: matrix4_zrot()
1456// Usage:
1457// matrix4_zrot(ang)
1458// Description:
1459// Returns the 4x4 matrix to perform a rotation of a 3D vector around the Z axis.
1460// Arguments:
1461// ang = number of degrees to rotate.
1462function matrix4_zrot(ang) = [
1463 [cos(ang), -sin(ang), 0, 0],
1464 [sin(ang), cos(ang), 0, 0],
1465 [ 0, 0, 1, 0],
1466 [ 0, 0, 0, 1]
1467];
1468
1469
1470// Function: matrix4_rot_by_axis()
1471// Usage:
1472// matrix4_rot_by_axis(u, ang);
1473// Description:
1474// Returns the 4x4 matrix to perform a rotation of a 3D vector around an axis.
1475// Arguments:
1476// u = 3D axis vector to rotate around.
1477// ang = number of degrees to rotate.
1478function matrix4_rot_by_axis(u, ang) = let(
1479 u = normalize(u),
1480 c = cos(ang),
1481 c2 = 1-c,
1482 s = sin(ang)
1483) [
1484 [u[0]*u[0]*c2+c , u[0]*u[1]*c2-u[2]*s, u[0]*u[2]*c2+u[1]*s, 0],
1485 [u[1]*u[0]*c2+u[2]*s, u[1]*u[1]*c2+c , u[1]*u[2]*c2-u[0]*s, 0],
1486 [u[2]*u[0]*c2-u[1]*s, u[2]*u[1]*c2+u[0]*s, u[2]*u[2]*c2+c , 0],
1487 [ 0, 0, 0, 1]
1488];
1489
1490
1491// Function: matrix3_skew()
1492// Usage:
1493// matrix3_skew(xa, ya)
1494// Description:
1495// Returns the 3x3 matrix to skew a 2D vector along the XY plane.
1496// Arguments:
1497// xa = Skew angle, in degrees, in the direction of the X axis.
1498// ya = Skew angle, in degrees, in the direction of the Y axis.
1499function matrix3_skew(xa, ya) = [
1500 [1, tan(xa), 0],
1501 [tan(ya), 1, 0],
1502 [0, 0, 1]
1503];
1504
1505
1506
1507// Function: matrix4_skew_xy()
1508// Usage:
1509// matrix4_skew_xy(xa, ya)
1510// Description:
1511// Returns the 4x4 matrix to perform a skew transformation along the XY plane..
1512// Arguments:
1513// xa = Skew angle, in degrees, in the direction of the X axis.
1514// ya = Skew angle, in degrees, in the direction of the Y axis.
1515function matrix4_skew_xy(xa, ya) = [
1516 [1, 0, tan(xa), 0],
1517 [0, 1, tan(ya), 0],
1518 [0, 0, 1, 0],
1519 [0, 0, 0, 1]
1520];
1521
1522
1523
1524// Function: matrix4_skew_xz()
1525// Usage:
1526// matrix4_skew_xz(xa, za)
1527// Description:
1528// Returns the 4x4 matrix to perform a skew transformation along the XZ plane.
1529// Arguments:
1530// xa = Skew angle, in degrees, in the direction of the X axis.
1531// za = Skew angle, in degrees, in the direction of the Z axis.
1532function matrix4_skew_xz(xa, za) = [
1533 [1, tan(xa), 0, 0],
1534 [0, 1, 0, 0],
1535 [0, tan(za), 1, 0],
1536 [0, 0, 0, 1]
1537];
1538
1539
1540// Function: matrix4_skew_yz()
1541// Usage:
1542// matrix4_skew_yz(ya, za)
1543// Description:
1544// Returns the 4x4 matrix to perform a skew transformation along the YZ plane.
1545// Arguments:
1546// ya = Skew angle, in degrees, in the direction of the Y axis.
1547// za = Skew angle, in degrees, in the direction of the Z axis.
1548function matrix4_skew_yz(ya, za) = [
1549 [ 1, 0, 0, 0],
1550 [tan(ya), 1, 0, 0],
1551 [tan(za), 0, 1, 0],
1552 [ 0, 0, 0, 1]
1553];
1554
1555
1556// Function: matrix3_mult()
1557// Usage:
1558// matrix3_mult(matrices)
1559// Description:
1560// Returns a 3x3 transformation matrix which results from applying each matrix in `matrices` in order.
1561// Arguments:
1562// matrices = A list of 3x3 matrices.
1563// m = Optional starting matrix to apply everything to.
1564function matrix3_mult(matrices, m=ident(3), i=0) =
1565 (i>=len(matrices))? m :
1566 let (newmat = is_def(m)? matrices[i] * m : matrices[i])
1567 matrix3_mult(matrices, m=newmat, i=i+1);
1568
1569
1570// Function: matrix4_mult()
1571// Usage:
1572// matrix4_mult(matrices)
1573// Description:
1574// Returns a 4x4 transformation matrix which results from applying each matrix in `matrices` in order.
1575// Arguments:
1576// matrices = A list of 4x4 matrices.
1577// m = Optional starting matrix to apply everything to.
1578function matrix4_mult(matrices, m=ident(4), i=0) =
1579 (i>=len(matrices))? m :
1580 let (newmat = is_def(m)? matrices[i] * m : matrices[i])
1581 matrix4_mult(matrices, m=newmat, i=i+1);
1582
1583
1584// Function: matrix3_apply()
1585// Usage:
1586// matrix3_apply(pts, matrices)
1587// Description:
1588// Given a list of transformation matrices, applies them in order to the points in the point list.
1589// Arguments:
1590// pts = A list of 2D points to transform.
1591// matrices = A list of 3x3 matrices to apply, in order.
1592// Example:
1593// npts = matrix3_apply(
1594// pts = [for (x=[0:3]) [5*x,0]],
1595// matrices =[
1596// matrix3_scale([3,1]),
1597// matrix3_rot(90),
1598// matrix3_translate([5,5])
1599// ]
1600// ); // Returns [[5,5], [5,20], [5,35], [5,50]]
1601function matrix3_apply(pts, matrices) = let(m = matrix3_mult(matrices)) [for (p = pts) point2d(m * concat(point2d(p),[1]))];
1602
1603
1604// Function: matrix4_apply()
1605// Usage:
1606// matrix4_apply(pts, matrices)
1607// Description:
1608// Given a list of transformation matrices, applies them in order to the points in the point list.
1609// Arguments:
1610// pts = A list of 3D points to transform.
1611// matrices = A list of 4x4 matrices to apply, in order.
1612// Example:
1613// npts = matrix4_apply(
1614// pts = [for (x=[0:3]) [5*x,0,0]],
1615// matrices =[
1616// matrix4_scale([2,1,1]),
1617// matrix4_zrot(90),
1618// matrix4_translate([5,5,10])
1619// ]
1620// ); // Returns [[5,5,10], [5,15,10], [5,25,10], [5,35,10]]
1621
1622function matrix4_apply(pts, matrices) = let(m = matrix4_mult(matrices)) [for (p = pts) point3d(m * concat(point3d(p),[1]))];
1623
1624
1625// Section: Geometry
1626
1627// Function: point_on_segment()
1628// Usage:
1629// point_on_segment(point, edge);
1630// Description:
1631// Determine if the point is on the line segment between two points.
1632// Returns true if yes, and false if not.
1633// Arguments:
1634// point = The point to check colinearity of.
1635// edge = Array of two points forming the line segment to test against.
1636function point_on_segment(point, edge) =
1637 point==edge[0] || point==edge[1] || // The point is an endpoint
1638 sign(edge[0].x-point.x)==sign(point.x-edge[1].x) // point is in between the
1639 && sign(edge[0].y-point.y)==sign(point.y-edge[1].y) // edge endpoints
1640 && point_left_of_segment(point, edge)==0; // and on the line defined by edge
1641
1642
1643// Function: point_left_of_segment()
1644// Usage:
1645// point_left_of_segment(point, edge);
1646// Description:
1647// Return >0 if point is left of the line defined by edge.
1648// Return =0 if point is on the line.
1649// Return <0 if point is right of the line.
1650// Arguments:
1651// point = The point to check position of.
1652// edge = Array of two points forming the line segment to test against.
1653function point_left_of_segment(point, edge) =
1654 (edge[1].x-edge[0].x) * (point.y-edge[0].y) - (point.x-edge[0].x) * (edge[1].y-edge[0].y);
1655
1656
1657// Internal non-exposed function.
1658function _point_above_below_segment(point, edge) =
1659 edge[0].y <= point.y? (
1660 (edge[1].y > point.y && point_left_of_segment(point, edge) > 0)? 1 : 0
1661 ) : (
1662 (edge[1].y <= point.y && point_left_of_segment(point, edge) < 0)? -1 : 0
1663 );
1664
1665
1666// Function: point_in_polygon()
1667// Usage:
1668// point_in_polygon(point, path)
1669// Description:
1670// This function tests whether the given point is inside, outside or on the boundary of
1671// the specified polygon using the Winding Number method. (http://geomalgorithms.com/a03-_inclusion.html)
1672// The polygon is given as a list of points, not including the repeated end point.
1673// Returns -1 if the point is outside the polyon.
1674// Returns 0 if the point is on the boundary.
1675// Returns 1 if the point lies in the interior.
1676// The polygon does not need to be simple: it can have self-intersections.
1677// But the polygon cannot have holes (it must be simply connected).
1678// Rounding error may give mixed results for points on or near the boundary.
1679// Arguments:
1680// point = The point to check position of.
1681// path = The list of 2D path points forming the perimeter of the polygon.
1682function point_in_polygon(point, path) =
1683 // Does the point lie on any edges? If so return 0.
1684 sum([for(i=[0:len(path)-1]) point_on_segment(point, select(path, i, i+1))?1:0])>0 ? 0 :
1685 // Otherwise compute winding number and return 1 for interior, -1 for exterior
1686 sum([for(i=[0:len(path)-1]) _point_above_below_segment(point, select(path, i, i+1))]) != 0 ? 1 : -1;
1687
1688
1689// Function: pointlist_bounds()
1690// Usage:
1691// pointlist_bounds(pts);
1692// Description:
1693// Finds the bounds containing all the points in pts.
1694// Returns [[minx, miny, minz], [maxx, maxy, maxz]]
1695// Arguments:
1696// pts = List of points.
1697function pointlist_bounds(pts) = [
1698 [for (a=[0:2]) min([ for (x=pts) point3d(x)[a] ]) ],
1699 [for (a=[0:2]) max([ for (x=pts) point3d(x)[a] ]) ]
1700];
1701
1702
1703// Function: triangle_area2d()
1704// Usage:
1705// triangle_area2d(a,b,c);
1706// Description:
1707// Returns the area of a triangle formed between three vertices.
1708// Result will be negative if the points are in clockwise order.
1709// Examples:
1710// triangle_area2d([0,0], [5,10], [10,0]); // Returns -50
1711// triangle_area2d([10,0], [5,10], [0,0]); // Returns 50
1712function triangle_area2d(a,b,c) =
1713 (
1714 a.x * (b.y - c.y) +
1715 b.x * (c.y - a.y) +
1716 c.x * (a.y - b.y)
1717 ) / 2;
1718
1719
1720// Function: right_of_line2d()
1721// Usage:
1722// right_of_line2d(line, pt)
1723// Description:
1724// Returns true if the given point is to the left of the given line.
1725// Arguments:
1726// line = A list of two points.
1727// pt = The point to test.
1728function right_of_line2d(line, pt) =
1729 triangle_area2d(line[0], line[1], pt) < 0;
1730
1731
1732// Function: collinear()
1733// Usage:
1734// collinear(a, b, c, [eps]);
1735// Description:
1736// Returns true if three points are co-linear.
1737// Arguments:
1738// a = First point.
1739// b = Second point.
1740// c = Third point.
1741// eps = Acceptable max angle variance. Default: EPSILON (1e-9) degrees.
1742function collinear(a, b, c, eps=EPSILON) =
1743 abs(vector_angle(b-a,c-a)) < eps;
1744
1745
1746// Function: collinear_indexed()
1747// Usage:
1748// collinear_indexed(points, a, b, c, [eps]);
1749// Description:
1750// Returns true if three points are co-linear.
1751// Arguments:
1752// points = A list of points.
1753// a = Index in `points` of first point.
1754// b = Index in `points` of second point.
1755// c = Index in `points` of third point.
1756// eps = Acceptable max angle variance. Default: EPSILON (1e-9) degrees.
1757function collinear_indexed(points, a, b, c, eps=EPSILON) =
1758 let(
1759 p1=points[a],
1760 p2=points[b],
1761 p3=points[c]
1762 ) abs(vector_angle(p2-p1,p3-p1)) < eps;
1763
1764
1765// Function: plane3pt()
1766// Usage:
1767// plane3pt(p1, p2, p3);
1768// Description:
1769// Generates the cartesian equation of a plane from three non-colinear points on the plane.
1770// Returns [A,B,C,D] where Ax+By+Cz+D=0 is the equation of a plane.
1771// Arguments:
1772// p1 = The first point on the plane.
1773// p2 = The second point on the plane.
1774// p3 = The third point on the plane.
1775function plane3pt(p1, p2, p3) =
1776 let(normal = normalize(cross(p3-p1, p2-p1))) concat(normal, [normal*p1]);
1777
1778
1779// Function: plane3pt_indexed()
1780// Usage:
1781// plane3pt_indexed(points, i1, i2, i3);
1782// Description:
1783// Given a list of points, and the indexes of three of those points,
1784// generates the cartesian equation of a plane that those points all
1785// lie on. Requires that the three indexed points be non-collinear.
1786// Returns [A,B,C,D] where Ax+By+Cz+D=0 is the equation of a plane.
1787// Arguments:
1788// points = A list of points.
1789// i1 = The index into `points` of the first point on the plane.
1790// i2 = The index into `points` of the second point on the plane.
1791// i3 = The index into `points` of the third point on the plane.
1792function plane3pt_indexed(points, i1, i2, i3) =
1793 let(
1794 p1 = points[i1],
1795 p2 = points[i2],
1796 p3 = points[i3],
1797 normal = normalize(cross(p3-p1, p2-p1))
1798 ) concat(normal, [normal*p1]);
1799
1800
1801// Function: distance_from_plane()
1802// Usage:
1803// distance_from_plane(plane, point)
1804// Description:
1805// Given a plane as [A,B,C,D] where the cartesian equation for that plane
1806// is Ax+By+Cz+D=0, determines how far from that plane the given point is.
1807// The returned distance will be positive if the point is in front of the
1808// plane; on the same side of the plane as the normal of that plane points
1809// towards. If the point is behind the plane, then the distance returned
1810// will be negative. The normal of the plane is the same as [A,B,C].
1811// Arguments:
1812// plane = The [A,B,C,D] values for the equation of the plane.
1813// point = The point to test.
1814function distance_from_plane(plane, point) =
1815 [plane.x, plane.y, plane.z] * point - plane[3];
1816
1817
1818// Function: coplanar()
1819// Usage:
1820// coplanar(plane, point);
1821// Description:
1822// Given a plane as [A,B,C,D] where the cartesian equation for that plane
1823// is Ax+By+Cz+D=0, determines if the given point is on that plane.
1824// Returns true if the point is on that plane.
1825// Arguments:
1826// plane = The [A,B,C,D] values for the equation of the plane.
1827// point = The point to test.
1828function coplanar(plane, point) =
1829 abs(distance_from_plane(plane, point)) <= EPSILON;
1830
1831
1832// Function: in_front_of_plane()
1833// Usage:
1834// in_front_of_plane(plane, point);
1835// Description:
1836// Given a plane as [A,B,C,D] where the cartesian equation for that plane
1837// is Ax+By+Cz+D=0, determines if the given point is on the side of that
1838// plane that the normal points towards. The normal of the plane is the
1839// same as [A,B,C].
1840// Arguments:
1841// plane = The [A,B,C,D] values for the equation of the plane.
1842// point = The point to test.
1843function in_front_of_plane(plane, point) =
1844 distance_from_plane(plane, point) > EPSILON;
1845
1846
1847// Function: simplify_path()
1848// Description:
1849// Takes a path and removes unnecessary collinear points.
1850// Usage:
1851// simplify_path(path, [eps])
1852// Arguments:
1853// path = A list of 2D path points.
1854// eps = Largest angle variance allowed. Default: EPSILON (1-e9) degrees.
1855function simplify_path(path, eps=EPSILON, _a=0, _b=2, _acc=[]) =
1856 (_b >= len(path))? concat([path[0]], _acc, [path[len(path)-1]]) :
1857 simplify_path(
1858 path, eps,
1859 (collinear_indexed(path, _a, _b-1, _b, eps=eps)? _a : _b-1),
1860 _b+1,
1861 (collinear_indexed(path, _a, _b-1, _b, eps=eps)? _acc : concat(_acc, [path[_b-1]]))
1862 );
1863
1864
1865// Function: simplify_path_indexed()
1866// Description:
1867// Takes a list of points, and a path as a list of indexes into `points`,
1868// and removes all path points that are unecessarily collinear.
1869// Usage:
1870// simplify_path_indexed(path, eps)
1871// Arguments:
1872// points = A list of points.
1873// path = A list of indexes into `points` that forms a path.
1874// eps = Largest angle variance allowed. Default: EPSILON (1-e9) degrees.
1875function simplify_path_indexed(points, path, eps=EPSILON, _a=0, _b=2, _acc=[]) =
1876 (_b >= len(path))? concat([path[0]], _acc, [path[len(path)-1]]) :
1877 simplify_path_indexed(
1878 points, path, eps,
1879 (collinear_indexed(points, path[_a], path[_b-1], path[_b], eps=eps)? _a : _b-1),
1880 _b+1,
1881 (collinear_indexed(points, path[_a], path[_b-1], path[_b], eps=eps)? _acc : concat(_acc, [path[_b-1]]))
1882 );
1883
1884
1885// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap