1//////////////////////////////////////////////////////////////////////
2// LibFile: vectors.scad
3// This file provides some mathematical operations that apply to each
4// entry in a vector. It provides normalization and angle computation, and
5// it provides functions for searching lists of vectors for matches to
6// a given vector.
7// Includes:
8// include <BOSL2/std.scad>
9// FileGroup: Math
10// FileSummary: Vector arithmetic, angle, and searching.
11// FileFootnotes: STD=Included in std.scad
12//////////////////////////////////////////////////////////////////////
13
14
15// Section: Vector Testing
16
17
18// Function: is_vector()
19// Synopsis: Returns true if the given value is a vector.
20// Topics: Vectors, Math
21// See Also: is_matrix(), is_path(), is_region()
22// Usage:
23// bool = is_vector(v, [length], [zero=], [all_nonzero=], [eps=]);
24// Description:
25// Returns true if v is a list of finite numbers.
26// Arguments:
27// v = The value to test to see if it is a vector.
28// length = If given, make sure the vector is `length` items long.
29// ---
30// zero = If false, require that the `norm()` of the vector is not approximately zero. If true, require the `norm()` of the vector to be approximately zero. Default: `undef` (don't check vector `norm()`.)
31// all_nonzero = If true, requires all elements of the vector to be more than `eps` different from zero. Default: `false`
32// eps = The minimum vector length that is considered non-zero. Default: `EPSILON` (`1e-9`)
33// Example:
34// is_vector(4); // Returns false
35// is_vector([4,true,false]); // Returns false
36// is_vector([3,4,INF,5]); // Returns false
37// is_vector([3,4,5,6]); // Returns true
38// is_vector([3,4,undef,5]); // Returns false
39// is_vector([3,4,5],3); // Returns true
40// is_vector([3,4,5],4); // Returns true
41// is_vector([]); // Returns false
42// is_vector([0,4,0],3,zero=false); // Returns true
43// is_vector([0,0,0],zero=false); // Returns false
44// is_vector([0,0,1e-12],zero=false); // Returns false
45// is_vector([0,1,0],all_nonzero=false); // Returns false
46// is_vector([1,1,1],all_nonzero=false); // Returns true
47// is_vector([],zero=false); // Returns false
48function is_vector(v, length, zero, all_nonzero=false, eps=EPSILON) =
49 is_list(v) && len(v)>0 && []==[for(vi=v) if(!is_finite(vi)) 0]
50 && (is_undef(length) || len(v)==length)
51 && (is_undef(zero) || ((norm(v) >= eps) == !zero))
52 && (!all_nonzero || all_nonzero(v)) ;
53
54
55
56// Section: Scalar operations on vectors
57
58// Function: add_scalar()
59// Synopsis: Adds a scalar value to every item in a vector.
60// Topics: Vectors, Math
61// See Also: add_scalar(), v_mul(), v_div()
62// Usage:
63// v_new = add_scalar(v, s);
64// Description:
65// Given a vector and a scalar, returns the vector with the scalar added to each item in it.
66// Arguments:
67// v = The initial array.
68// s = A scalar value to add to every item in the array.
69// Example:
70// a = add_scalar([1,2,3],3); // Returns: [4,5,6]
71function add_scalar(v,s) =
72 assert(is_vector(v), "Input v must be a vector")
73 assert(is_finite(s), "Input s must be a finite scalar")
74 [for(entry=v) entry+s];
75
76
77// Function: v_mul()
78// Synopsis: Returns the element-wise multiplication of two equal-length vectors.
79// Topics: Vectors, Math
80// See Also: add_scalar(), v_mul(), v_div()
81// Usage:
82// v3 = v_mul(v1, v2);
83// Description:
84// Element-wise multiplication. Multiplies each element of `v1` by the corresponding element of `v2`.
85// Both `v1` and `v2` must be the same length. Returns a vector of the products. Note that
86// the items in `v1` and `v2` can be anything that OpenSCAD will multiply.
87// Arguments:
88// v1 = The first vector.
89// v2 = The second vector.
90// Example:
91// v_mul([3,4,5], [8,7,6]); // Returns [24, 28, 30]
92function v_mul(v1, v2) =
93 assert( is_list(v1) && is_list(v2) && len(v1)==len(v2), "Incompatible input")
94 [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
95
96
97// Function: v_div()
98// Synopsis: Returns the element-wise division of two equal-length vectors.
99// Topics: Vectors, Math
100// See Also: add_scalar(), v_mul(), v_div()
101// Usage:
102// v3 = v_div(v1, v2);
103// Description:
104// Element-wise vector division. Divides each element of vector `v1` by
105// the corresponding element of vector `v2`. Returns a vector of the quotients.
106// Arguments:
107// v1 = The first vector.
108// v2 = The second vector.
109// Example:
110// v_div([24,28,30], [8,7,6]); // Returns [3, 4, 5]
111function v_div(v1, v2) =
112 assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
113 [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
114
115
116// Function: v_abs()
117// Synopsis: Returns the absolute values of the given vector.
118// Topics: Vectors, Math
119// See Also: v_abs(), v_floor(), v_ceil()
120// Usage:
121// v2 = v_abs(v);
122// Description: Returns a vector of the absolute value of each element of vector `v`.
123// Arguments:
124// v = The vector to get the absolute values of.
125// Example:
126// v_abs([-1,3,-9]); // Returns: [1,3,9]
127function v_abs(v) =
128 assert( is_vector(v), "Invalid vector" )
129 [for (x=v) abs(x)];
130
131
132// Function: v_floor()
133// Synopsis: Returns the values of the given vector, rounded down.
134// Topics: Vectors, Math
135// See Also: v_abs(), v_floor(), v_ceil()
136// Usage:
137// v2 = v_floor(v);
138// Description:
139// Returns the given vector after performing a `floor()` on all items.
140function v_floor(v) =
141 assert( is_vector(v), "Invalid vector" )
142 [for (x=v) floor(x)];
143
144
145// Function: v_ceil()
146// Synopsis: Returns the values of the given vector, rounded up.
147// Topics: Vectors, Math
148// See Also: v_abs(), v_floor(), v_ceil()
149// Usage:
150// v2 = v_ceil(v);
151// Description:
152// Returns the given vector after performing a `ceil()` on all items.
153function v_ceil(v) =
154 assert( is_vector(v), "Invalid vector" )
155 [for (x=v) ceil(x)];
156
157
158// Function: v_lookup()
159// Synopsis: Like `lookup()`, but it can interpolate between vector results.
160// Topics: Vectors, Math
161// See Also: v_abs(), v_floor(), v_ceil()
162// Usage:
163// v2 = v_lookup(x, v);
164// Description:
165// Works just like the built-in function [`lookup()`](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Mathematical_Functions#lookup), except that it can also interpolate between vector result values of the same length.
166// Arguments:
167// x = The scalar value to look up.
168// v = A list of [KEY,VAL] pairs. KEYs are scalars. VALs should either all be scalar, or all be vectors of the same length.
169// Example:
170// x = v_lookup(4.5, [[4, [3,4,5]], [5, [5,6,7]]]); // Returns: [4,5,6]
171function v_lookup(x, v) =
172 is_num(v[0][1])? lookup(x,v) :
173 let(
174 i = lookup(x, [for (i=idx(v)) [v[i].x,i]]),
175 vlo = v[floor(i)],
176 vhi = v[ceil(i)],
177 lo = vlo[1],
178 hi = vhi[1]
179 )
180 assert(is_vector(lo) && is_vector(hi),
181 "Result values must all be numbers, or all be vectors.")
182 assert(len(lo) == len(hi), "Vector result values must be the same length")
183 vlo.x == vhi.x? vlo[1] :
184 let( u = (x - vlo.x) / (vhi.x - vlo.x) )
185 lerp(lo,hi,u);
186
187
188// Section: Vector Properties
189
190
191// Function: unit()
192// Synopsis: Returns the unit length of a given vector.
193// Topics: Vectors, Math
194// See Also: v_abs(), v_floor(), v_ceil()
195// Usage:
196// v = unit(v, [error]);
197// Description:
198// Returns the unit length normalized version of vector v. If passed a zero-length vector,
199// asserts an error unless `error` is given, in which case the value of `error` is returned.
200// Arguments:
201// v = The vector to normalize.
202// error = If given, and input is a zero-length vector, this value is returned. Default: Assert error on zero-length vector.
203// Example:
204// v1 = unit([10,0,0]); // Returns: [1,0,0]
205// v2 = unit([0,10,0]); // Returns: [0,1,0]
206// v3 = unit([0,0,10]); // Returns: [0,0,1]
207// v4 = unit([0,-10,0]); // Returns: [0,-1,0]
208// v5 = unit([0,0,0],[1,2,3]); // Returns: [1,2,3]
209// v6 = unit([0,0,0]); // Asserts an error.
210function unit(v, error=[[["ASSERT"]]]) =
211 assert(is_vector(v), "Invalid vector")
212 norm(v)<EPSILON? (error==[[["ASSERT"]]]? assert(norm(v)>=EPSILON,"Cannot normalize a zero vector") : error) :
213 v/norm(v);
214
215
216// Function: v_theta()
217// Synopsis: Returns the angle counter-clockwise from X+ on the XY plane.
218// Topics: Vectors, Math
219// See Also: unit()
220// Usage:
221// theta = v_theta([X,Y]);
222// Description:
223// Given a vector, returns the angle in degrees counter-clockwise from X+ on the XY plane.
224function v_theta(v) =
225 assert( is_vector(v,2) || is_vector(v,3) , "Invalid vector")
226 atan2(v.y,v.x);
227
228
229
230// Function: vector_angle()
231// Synopsis: Returns the minor angle between two vectors.
232// Topics: Vectors, Math
233// See Also: unit(), v_theta()
234// Usage:
235// ang = vector_angle(v1,v2);
236// ang = vector_angle([v1,v2]);
237// ang = vector_angle(PT1,PT2,PT3);
238// ang = vector_angle([PT1,PT2,PT3]);
239// Description:
240// If given a single list of two vectors, like `vector_angle([V1,V2])`, returns the angle between the two vectors V1 and V2.
241// If given a single list of three points, like `vector_angle([A,B,C])`, returns the angle between the line segments AB and BC.
242// If given two vectors, like `vector_angle(V1,V2)`, returns the angle between the two vectors V1 and V2.
243// If given three points, like `vector_angle(A,B,C)`, returns the angle between the line segments AB and BC.
244// Arguments:
245// v1 = First vector or point.
246// v2 = Second vector or point.
247// v3 = Third point in three point mode.
248// Example:
249// ang1 = vector_angle(UP,LEFT); // Returns: 90
250// ang2 = vector_angle(RIGHT,LEFT); // Returns: 180
251// ang3 = vector_angle(UP+RIGHT,RIGHT); // Returns: 45
252// ang4 = vector_angle([10,10], [0,0], [10,-10]); // Returns: 90
253// ang5 = vector_angle([10,0,10], [0,0,0], [-10,10,0]); // Returns: 120
254// ang6 = vector_angle([[10,0,10], [0,0,0], [-10,10,0]]); // Returns: 120
255function vector_angle(v1,v2,v3) =
256 assert( ( is_undef(v3) && ( is_undef(v2) || same_shape(v1,v2) ) )
257 || is_consistent([v1,v2,v3]) ,
258 "Bad arguments.")
259 assert( is_vector(v1) || is_consistent(v1), "Bad arguments.")
260 let( vecs = ! is_undef(v3) ? [v1-v2,v3-v2] :
261 ! is_undef(v2) ? [v1,v2] :
262 len(v1) == 3 ? [v1[0]-v1[1], v1[2]-v1[1]]
263 : v1
264 )
265 assert(is_vector(vecs[0],2) || is_vector(vecs[0],3), "Bad arguments.")
266 let(
267 norm0 = norm(vecs[0]),
268 norm1 = norm(vecs[1])
269 )
270 assert(norm0>0 && norm1>0, "Zero length vector.")
271 // NOTE: constrain() corrects crazy FP rounding errors that exceed acos()'s domain.
272 acos(constrain((vecs[0]*vecs[1])/(norm0*norm1), -1, 1));
273
274
275// Function: vector_axis()
276// Synopsis: Returns the perpendicular axis between two vectors.
277// Topics: Vectors, Math
278// See Also: unit(), v_theta(), vector_angle()
279// Usage:
280// axis = vector_axis(v1,v2);
281// axis = vector_axis([v1,v2]);
282// axis = vector_axis(PT1,PT2,PT3);
283// axis = vector_axis([PT1,PT2,PT3]);
284// Description:
285// If given a single list of two vectors, like `vector_axis([V1,V2])`, returns the vector perpendicular the two vectors V1 and V2.
286// If given a single list of three points, like `vector_axis([A,B,C])`, returns the vector perpendicular to the plane through a, B and C.
287// If given two vectors, like `vector_axis(V1,V2)`, returns the vector perpendicular to the two vectors V1 and V2.
288// If given three points, like `vector_axis(A,B,C)`, returns the vector perpendicular to the plane through a, B and C.
289// Arguments:
290// v1 = First vector or point.
291// v2 = Second vector or point.
292// v3 = Third point in three point mode.
293// Example:
294// axis1 = vector_axis(UP,LEFT); // Returns: [0,-1,0] (FWD)
295// axis2 = vector_axis(RIGHT,LEFT); // Returns: [0,-1,0] (FWD)
296// axis3 = vector_axis(UP+RIGHT,RIGHT); // Returns: [0,1,0] (BACK)
297// axis4 = vector_axis([10,10], [0,0], [10,-10]); // Returns: [0,0,-1] (DOWN)
298// axis5 = vector_axis([10,0,10], [0,0,0], [-10,10,0]); // Returns: [-0.57735, -0.57735, 0.57735]
299// axis6 = vector_axis([[10,0,10], [0,0,0], [-10,10,0]]); // Returns: [-0.57735, -0.57735, 0.57735]
300function vector_axis(v1,v2=undef,v3=undef) =
301 is_vector(v3)
302 ? assert(is_consistent([v3,v2,v1]), "Bad arguments.")
303 vector_axis(v1-v2, v3-v2)
304 : assert( is_undef(v3), "Bad arguments.")
305 is_undef(v2)
306 ? assert( is_list(v1), "Bad arguments.")
307 len(v1) == 2
308 ? vector_axis(v1[0],v1[1])
309 : vector_axis(v1[0],v1[1],v1[2])
310 : assert( is_vector(v1,zero=false) && is_vector(v2,zero=false) && is_consistent([v1,v2])
311 , "Bad arguments.")
312 let(
313 eps = 1e-6,
314 w1 = point3d(v1/norm(v1)),
315 w2 = point3d(v2/norm(v2)),
316 w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2
317 : (norm(v_abs(w2)-UP) > eps)? UP
318 : RIGHT
319 ) unit(cross(w1,w3));
320
321
322// Function: vector_bisect()
323// Synopsis: Returns the vector that bisects two vectors.
324// Topics: Vectors, Math
325// See Also: unit(), v_theta(), vector_angle(), vector_axis()
326// Usage:
327// newv = vector_bisect(v1,v2);
328// Description:
329// Returns a unit vector that exactly bisects the minor angle between two given vectors.
330// If given two vectors that are directly opposed, returns `undef`.
331function vector_bisect(v1,v2) =
332 assert(is_vector(v1))
333 assert(is_vector(v2))
334 assert(!approx(norm(v1),0), "Zero length vector.")
335 assert(!approx(norm(v2),0), "Zero length vector.")
336 assert(len(v1)==len(v2), "Vectors are of different sizes.")
337 let( v1 = unit(v1), v2 = unit(v2) )
338 approx(v1,-v2)? undef :
339 let(
340 axis = vector_axis(v1,v2),
341 ang = vector_angle(v1,v2),
342 v3 = unit(rot(ang/2, v=axis, p=v1))
343 ) v3;
344
345
346
347// Section: Vector Searching
348
349
350// Function: pointlist_bounds()
351// Synopsis: Returns the min and max bounding coordinates for the given list of points.
352// Topics: Geometry, Bounding Boxes, Bounds
353// See Also: closest_point()
354// Usage:
355// pt_pair = pointlist_bounds(pts);
356// Description:
357// Finds the bounds containing all the points in `pts` which can be a list of points in any dimension.
358// Returns a list of two items: a list of the minimums and a list of the maximums. For example, with
359// 3d points `[[MINX, MINY, MINZ], [MAXX, MAXY, MAXZ]]`
360// Arguments:
361// pts = List of points.
362function pointlist_bounds(pts) =
363 assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." )
364 let(
365 select = ident(len(pts[0])),
366 spread = [
367 for(i=[0:len(pts[0])-1])
368 let( spreadi = pts*select[i] )
369 [ min(spreadi), max(spreadi) ]
370 ]
371 ) transpose(spread);
372
373
374
375// Function: closest_point()
376// Synopsis: Finds the closest point in a list of points.
377// Topics: Geometry, Points, Distance
378// See Also: pointlist_bounds(), furthest_point(), closest_point()
379// Usage:
380// index = closest_point(pt, points);
381// Description:
382// Given a list of `points`, finds the index of the closest point to `pt`.
383// Arguments:
384// pt = The point to find the closest point to.
385// points = The list of points to search.
386function closest_point(pt, points) =
387 assert( is_vector(pt), "Invalid point." )
388 assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
389 min_index([for (p=points) norm(p-pt)]);
390
391
392// Function: furthest_point()
393// Synopsis: Finds the furthest point in a list of points.
394// Topics: Geometry, Points, Distance
395// See Also: pointlist_bounds(), furthest_point(), closest_point()
396// Usage:
397// index = furthest_point(pt, points);
398// Description:
399// Given a list of `points`, finds the index of the furthest point from `pt`.
400// Arguments:
401// pt = The point to find the farthest point from.
402// points = The list of points to search.
403function furthest_point(pt, points) =
404 assert( is_vector(pt), "Invalid point." )
405 assert(is_path(points,dim=len(pt)), "Invalid pointlist or incompatible dimensions." )
406 max_index([for (p=points) norm(p-pt)]);
407
408
409// Function: vector_search()
410// Synopsis: Finds points in a list that are close to a given point.
411// Topics: Search, Points, Closest
412// See Also: vector_search_tree(), vector_nearest()
413// Usage:
414// indices = vector_search(query, r, target);
415// Description:
416// Given a list of query points `query` and a `target` to search,
417// finds the points in `target` that match each query point. A match holds when the
418// distance between a point in `target` and a query point is less than or equal to `r`.
419// The returned list will have a list for each query point containing, in arbitrary
420// order, the indices of all points that match that query point.
421// The `target` may be a simple list of points or a search tree.
422// When `target` is a large list of points, a search tree is constructed to
423// speed up the search with an order around O(log n) per query point.
424// For small point lists, a direct search is done dispensing a tree construction.
425// Alternatively, `target` may be a search tree built with `vector_search_tree()`.
426// In that case, that tree is parsed looking for matches.
427// An empty list of query points will return a empty output list.
428// An empty list of target points will return a output list with an empty list for each query point.
429// Arguments:
430// query = list of points to find matches for.
431// r = the search radius.
432// target = list of the points to search for matches or a search tree.
433// Example: A set of four queries to find points within 1 unit of the query. The circles show the search region and all have radius 1.
434// $fn=32;
435// k = 2000;
436// points = list_to_matrix(rands(0,10,k*2,seed=13333),2);
437// queries = [for(i=[3,7],j=[3,7]) [i,j]];
438// search_ind = vector_search(queries, points, 1);
439// move_copies(points) circle(r=.08);
440// for(i=idx(queries)){
441// color("blue")stroke(move(queries[i],circle(r=1)), closed=true, width=.08);
442// color("red") move_copies(select(points, search_ind[i])) circle(r=.08);
443// }
444// Example: when a series of searches with different radius are needed, its is faster to pre-compute the tree
445// $fn=32;
446// k = 2000;
447// points = list_to_matrix(rands(0,10,k*2),2,seed=13333);
448// queries1 = [for(i=[3,7]) [i,i]];
449// queries2 = [for(i=[3,7]) [10-i,i]];
450// r1 = 1;
451// r2 = .7;
452// search_tree = vector_search_tree(points);
453// search_1 = vector_search(queries1, r1, search_tree);
454// search_2 = vector_search(queries2, r2, search_tree);
455// move_copies(points) circle(r=.08);
456// for(i=idx(queries1)){
457// color("blue")stroke(move(queries1[i],circle(r=r1)), closed=true, width=.08);
458// color("red") move_copies(select(points, search_1[i])) circle(r=.08);
459// }
460// for(i=idx(queries2)){
461// color("green")stroke(move(queries2[i],circle(r=r2)), closed=true, width=.08);
462// color("red") move_copies(select(points, search_2[i])) circle(r=.08);
463// }
464function vector_search(query, r, target) =
465 query==[] ? [] :
466 is_list(query) && target==[] ? is_vector(query) ? [] : [for(q=query) [] ] :
467 assert( is_finite(r) && r>=0,
468 "The query radius should be a positive number." )
469 let(
470 tgpts = is_matrix(target), // target is a point list
471 tgtree = is_list(target) // target is a tree
472 && (len(target)==2)
473 && is_matrix(target[0])
474 && is_list(target[1])
475 && (len(target[1])==4 || (len(target[1])==1 && is_list(target[1][0])) )
476 )
477 assert( tgpts || tgtree,
478 "The target should be a list of points or a search tree compatible with the query." )
479 let(
480 dim = tgpts ? len(target[0]) : len(target[0][0]),
481 simple = is_vector(query, dim)
482 )
483 assert( simple || is_matrix(query,undef,dim),
484 "The query points should be a list of points compatible with the target point list.")
485 tgpts
486 ? len(target)<=400
487 ? simple ? [for(i=idx(target)) if(norm(target[i]-query)<=r) i ] :
488 [for(q=query) [for(i=idx(target)) if(norm(target[i]-q)<=r) i ] ]
489 : let( tree = _bt_tree(target, count(len(target)), leafsize=25) )
490 simple ? _bt_search(query, r, target, tree) :
491 [for(q=query) _bt_search(q, r, target, tree)]
492 : simple ? _bt_search(query, r, target[0], target[1]) :
493 [for(q=query) _bt_search(q, r, target[0], target[1])];
494
495
496//Ball tree search
497function _bt_search(query, r, points, tree) =
498 assert( is_list(tree)
499 && ( ( len(tree)==1 && is_list(tree[0]) )
500 || ( len(tree)==4 && is_num(tree[0]) && is_num(tree[1]) ) ),
501 "The tree is invalid.")
502 len(tree)==1
503 ? assert( tree[0]==[] || is_vector(tree[0]), "The tree is invalid." )
504 [for(i=tree[0]) if(norm(points[i]-query)<=r) i ]
505 : norm(query-points[tree[0]]) > r+tree[1] ? [] :
506 concat(
507 [ if(norm(query-points[tree[0]])<=r) tree[0] ],
508 _bt_search(query, r, points, tree[2]),
509 _bt_search(query, r, points, tree[3]) ) ;
510
511
512// Function: vector_search_tree()
513// Synopsis: Makes a distance search tree for a list of points.
514// Topics: Search, Points, Closest
515// See Also: vector_nearest(), vector_search()
516// Usage:
517// tree = vector_search_tree(points,leafsize);
518// Description:
519// Construct a search tree for the given list of points to be used as input
520// to the function `vector_search()`. The use of a tree speeds up the
521// search process. The tree construction stops branching when
522// a tree node represents a number of points less or equal to `leafsize`.
523// Search trees are ball trees. Constructing the
524// tree should be O(n log n) and searches should be O(log n), though real life
525// performance depends on how the data is distributed, and it will deteriorate
526// for high data dimensions. This data structure is useful when you will be
527// performing many searches of the same data, so that the cost of constructing
528// the tree is justified. (See https://en.wikipedia.org/wiki/Ball_tree)
529// For a small lists of points, the search with a tree may be more expensive
530// than direct comparisons. The argument `treemin` sets the minimum length of
531// point set for which a tree search will be done by `vector_search`.
532// For an empty list of points it returns an empty list.
533// Arguments:
534// points = list of points to store in the search tree.
535// leafsize = the size of the tree leaves. Default: 25
536// treemin = the minimum size of the point list for which a tree search is done. Default: 400
537// Example: A set of four queries to find points within 1 unit of the query. The circles show the search region and all have radius 1.
538// $fn=32;
539// k = 2000;
540// points = random_points(k, scale=10, dim=2,seed=13333);
541// queries = [for(i=[3,7],j=[3,7]) [i,j]];
542// search_tree = vector_search_tree(points);
543// search_ind = vector_search(queries,1,search_tree);
544// move_copies(points) circle(r=.08);
545// for(i=idx(queries)){
546// color("blue") stroke(move(queries[i],circle(r=1)), closed=true, width=.08);
547// color("red") move_copies(select(points, search_ind[i])) circle(r=.08);
548// }
549function vector_search_tree(points, leafsize=25, treemin=400) =
550 points==[] ? [] :
551 assert( is_matrix(points), "The input list entries should be points." )
552 assert( is_int(leafsize) && leafsize>=1,
553 "The tree leaf size should be an integer greater than zero.")
554 len(points)<treemin ? points :
555 [ points, _bt_tree(points, count(len(points)), leafsize) ];
556
557
558//Ball tree construction
559function _bt_tree(points, ind, leafsize=25) =
560 len(ind)<=leafsize ? [ind] :
561 let(
562 bounds = pointlist_bounds(select(points,ind)),
563 coord = max_index(bounds[1]-bounds[0]),
564 projc = [for(i=ind) points[i][coord] ],
565 meanpr = mean(projc),
566 pivot = min_index([for(p=projc) abs(p-meanpr)]),
567 radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]),
568 Lind = [for(i=idx(ind)) if(projc[i]<=meanpr && i!=pivot) ind[i] ],
569 Rind = [for(i=idx(ind)) if(projc[i] >meanpr && i!=pivot) ind[i] ]
570 )
571 [ ind[pivot], radius, _bt_tree(points, Lind, leafsize), _bt_tree(points, Rind, leafsize) ];
572
573
574// Function: vector_nearest()
575// Synopsis: Finds the `k` nearest points in a list to a given point.
576// Topics: Search, Points, Closest
577// See Also: vector_search(), vector_search_tree()
578// Usage:
579// indices = vector_nearest(query, k, target);
580// Description:
581// Search `target` for the `k` points closest to point `query`.
582// The input `target` is either a list of points to search or a search tree
583// pre-computed by `vector_search_tree(). A list is returned containing the indices
584// of the points found in sorted order, closest point first.
585// Arguments:
586// query = point to search for
587// k = number of neighbors to return
588// target = a list of points or a search tree to search in
589// Example: Four queries to find the 15 nearest points. The circles show the radius defined by the most distant query result. Note they are different for each query.
590// $fn=32;
591// k = 1000;
592// points = list_to_matrix(rands(0,10,k*2,seed=13333),2);
593// tree = vector_search_tree(points);
594// queries = [for(i=[3,7],j=[3,7]) [i,j]];
595// search_ind = [for(q=queries) vector_nearest(q, 15, tree)];
596// move_copies(points) circle(r=.08);
597// for(i=idx(queries)){
598// circle = circle(r=norm(points[last(search_ind[i])]-queries[i]));
599// color("red") move_copies(select(points, search_ind[i])) circle(r=.08);
600// color("blue") stroke(move(queries[i], circle), closed=true, width=.08);
601// }
602function vector_nearest(query, k, target) =
603 assert(is_int(k) && k>0)
604 assert(is_vector(query), "Query must be a vector.")
605 let(
606 tgpts = is_matrix(target,undef,len(query)), // target is a point list
607 tgtree = is_list(target) // target is a tree
608 && (len(target)==2)
609 && is_matrix(target[0],undef,len(query))
610 && (len(target[1])==4 || (len(target[1])==1 && is_list(target[1][0])) )
611 )
612 assert( tgpts || tgtree,
613 "The target should be a list of points or a search tree compatible with the query." )
614 assert((tgpts && (k<=len(target))) || (tgtree && (k<=len(target[0]))),
615 "More results are requested than the number of points.")
616 tgpts
617 ? let( tree = _bt_tree(target, count(len(target))) )
618 column(_bt_nearest( query, k, target, tree),0)
619 : column(_bt_nearest( query, k, target[0], target[1]),0);
620
621
622//Ball tree nearest
623function _bt_nearest(p, k, points, tree, answers=[]) =
624 assert( is_list(tree)
625 && ( ( len(tree)==1 && is_list(tree[0]) )
626 || ( len(tree)==4 && is_num(tree[0]) && is_num(tree[1]) ) ),
627 "The tree is invalid.")
628 len(tree)==1
629 ? _insert_many(answers, k, [for(entry=tree[0]) [entry, norm(points[entry]-p)]])
630 : let( d = norm(p-points[tree[0]]) )
631 len(answers)==k && ( d > last(answers)[1]+tree[1] ) ? answers :
632 let(
633 answers1 = _insert_sorted(answers, k, [tree[0],d]),
634 answers2 = _bt_nearest(p, k, points, tree[2], answers1),
635 answers3 = _bt_nearest(p, k, points, tree[3], answers2)
636 )
637 answers3;
638
639
640function _insert_sorted(list, k, new) =
641 (len(list)==k && new[1]>= last(list)[1]) ? list
642 : [
643 for(entry=list) if (entry[1]<=new[1]) entry,
644 new,
645 for(i=[0:1:min(k-1,len(list))-1]) if (list[i][1]>new[1]) list[i]
646 ];
647
648
649function _insert_many(list, k, newlist,i=0) =
650 i==len(newlist)
651 ? list
652 : assert(is_vector(newlist[i],2), "The tree is invalid.")
653 _insert_many(_insert_sorted(list,k,newlist[i]),k,newlist,i+1);
654
655
656
657// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap