Add vantage point tree for searching vector lists. Add "concave"

method to vnf_vertex_array.

math.scad

@@ -744,6 +744,29 @@ function mean(v) =
     sum(v)/len(v);
 
 
+// Function: ninther()
+// Usage:
+//    med = ninther(v)
+// Description:
+//    Finds a value in the input list of numbers `v` that is the median of  a 
+//    sample of 9 entries of `v`.
+//    It is a much faster approximation of the true median computation.
+// Arguments:
+//    v = an array of numbers
+function ninther(v) = 
+    let( l=len(v) )
+    l<=4 ? l<=2 ? v[0] : _med3(v[0], v[1], v[2]) : 
+    l==5 ? _med3(v[0], _med3(v[1], v[2], v[3]), v[4]) :
+    _med3(_med3(v[0],v[floor(l/6)],v[floor(l/3)]),
+          _med3(v[floor(l/3)],v[floor(l/2)],v[floor(2*l/3)]),
+          _med3(v[floor(2*l/3)],v[floor((5*l/3 -1)/2)],v[l-1]) );
+
+// the median of a triple
+function _med3(a,b,c) =
+    a < c ? a < b ? min(b,c) : min(a,c) :
+    b < c ? min(a,c) : min(a,b);
+
+
 // Function: convolve()
 // Usage:
 //   x = convolve(p,q);

vectors.scad

@@ -218,5 +218,141 @@ function vector_axis(v1,v2=undef,v3=undef) =
             ) unit(cross(w1,w3));
 
 
+// Section: Vector Searching
+
+// Function: vp_tree()
+// Usage:
+//    tree = vp_tree(points, <leafsize>)
+// Description:
+//    Organizes n-dimensional data into a Vantage Point Tree, which can be
+//    efficiently searched for for nearest matches.  The Vantage Point Tree
+//    is an effort to generalize binary search to n dimensions.  Constructing the
+//    tree should be O(n log n) and searches should be O(log n), though real life
+//    performance depends on how the data is distributed, and it will deteriorate
+//    for high data dimensions.  This data structure is useful when you will be
+//    performing many searches of the same data, so that the cost of constructing
+//    the tree is justified.  
+//    .
+//    The vantage point tree at a given level chooses vp, the
+//    "vantage point", and a radius, R, and divides the data based
+//    on distance to vp.  Points closer than R go in on branch
+//    of the tree and points farther than R go in the other branch.
+//    .
+//    The tree has the form [vp, R, inside, outside], where vp is
+//    the vantage point index, R is the radius, inside is a
+//    recursively computed tree for the inside points (distance less than
+//    or equal to R from the vantage point), and outside
+//    is a tree for the outside points (distance greater than R from the
+//    vantage point).
+//    .
+//    If the number of points is less than or equal to leafsize then
+//    vp_tree instead returns the list [ind] where ind is a list of
+//    the indices of the points.  This means the list has the form
+//    [[i0, i1, i2,...]], so tree[0] is a list of indices.  You can
+//    tell that a node is a leaf node by checking if tree[0] is a list.
+//    The leafsize parameter determines how many points can be
+//    store in the leaf nodes.  The default value of 25 was found
+//    emperically to be a reasonable option for 3d data searched with vp_search().
+//    .
+//    Vantage point tree is described here: http://web.cs.iastate.edu/~honavar/nndatastructures.pdf
+// Arguments:
+//    points = list of points to store in the tree
+//    leafsize = maximum number of points to store in the tree's leaf nodes.  Default: 25
+function vp_tree(points, leafsize=25) =
+  assert(is_matrix(points),"points must be a consistent list of data points")
+  _vp_tree(points, count(len(points)), leafsize);
+
+function _vp_tree(ptlist, ind, leafsize) =
+  len(ind)<=leafsize ? [ind] :
+  let(
+       center = mean(select(ptlist,ind)),
+       cdistances = [for(i=ind) norm(ptlist[i]-center)],
+       vpind = ind[max_index(cdistances)],
+       vp = ptlist[vpind],
+       vp_dist = [for(i=ind) norm(vp-ptlist[i])],
+       r = ninther(vp_dist),
+       inside = [for(i=idx(ind)) if (vp_dist[i]<=r && ind[i]!=vpind) ind[i]],
+       outside = [for(i=idx(ind)) if (vp_dist[i]>r) ind[i]]
+   )
+   [vpind, r, _vp_tree(ptlist,inside,leafsize),_vp_tree(ptlist,outside,leafsize)];
+
+
+// Function: vp_search()
+// Usage:
+//   indices = vp_search(points, tree, p, r);
+// Description:
+//   Search a vantage point tree for all points whose distance from p
+//   is less than or equal to r.  Returns a list of indices of the points it finds
+//   in arbitrary order.  The input points is a list of points to search and tree is the
+//   vantage point tree computed from that point list.  The search should be
+//   around O(log n).  
+// Arguments:
+//   points = points indexed by the vantage point tree
+//   tree = vantage point tree from vp_tree
+//   p = point to search for
+//   r = search radius
+function _vp_search(points, tree, p, r) =
+    is_list(tree[0]) ? [for(i=tree[0]) if (norm(points[i]-p)<=r) i]
+    :
+    let(
+        d = norm(p-points[tree[0]])  // dist to vantage point
+    )
+    [
+      if (d <= r) tree[0],
+      if (d-r <= tree[1]) each _vp_search(points, tree[2], p, r),
+      if (d+r > tree[1]) each _vp_search(points, tree[3], p, r)
+    ];
+
+function vp_search(points, tree, p, r) =
+    assert(is_list(tree) && (len(tree)==4 || (len(tree)==1 && is_list(tree[0]))), "Vantage point tree not valid")
+    assert(is_matrix(points), "Parameter points is not a consistent point list")
+    assert(is_vector(p,len(points[0])), "Query must be a vector whose length matches the point list")
+    assert(all_positive(r),"Radius r must be a positive number")
+    _vp_search(points, tree, p, r);
+  
+
+// Function: vp_nearest()
+// Usage:
+//    indices = vp_nearest(points, tree, p, k)
+// Description:
+//    Search the vantage point tree for the k points closest to point p.
+//    The input points is the list of points to search and tree is
+//    the vantage point tree computed from that point list.  The list is
+//    returned in sorted order, closest point first.  
+// Arguments:
+//    points = points indexed by the vantage point tree
+//    tree = vantage point tree from vp_tree
+//    p = point to search for
+//    k = number of neighbors to return
+function _insert_sorted(list, k, new) =
+      len(list)==k && new[1]>= last(list)[1] ? list
+    : [
+        for(entry=list) if (entry[1]<=new[1]) entry,
+        new,
+        for(i=[0:1:min(k-1,len(list))-1]) if (list[i][1]>new[1]) list[i]
+      ];
+
+function _insert_many(list, k, newlist,i=0) =
+  i==len(newlist) ? list :
+                    _insert_many(_insert_sorted(list,k,newlist[i]),k,newlist,i+1);
+
+function _vp_nearest(points, tree, p, k, answers=[]) =
+  is_list(tree[0]) ? _insert_many(answers, k, [for(entry=tree[0]) [entry, norm(points[entry]-p)]]) :
+  let(
+      d = norm(p-points[tree[0]]),
+      answers1 = _insert_sorted(answers, k, [tree[0],d]),
+      answers2 = d-last(answers1)[1] <= tree[1] ? _vp_nearest(points, tree[2], p, k, answers1) : answers1,
+      answers3 = d+last(answers2)[1] > tree[1] ? _vp_nearest(points, tree[3], p, k, answers2) : answers2
+   )
+   answers3;
+
+function vp_nearest(points, tree, p, k) =
+      assert(is_int(k) && k>0)
+      assert(k<=len(points), "You requested more results that contained in the set")
+      assert(is_matrix(points), "Parameter points is not a consistent point list")
+      assert(is_vector(p,len(points[0])), "Query must be a vector whose length matches the point list")
+      assert(is_list(tree) && (len(tree)==4 || (len(tree)==1 && is_list(tree[0]))), "Vantage point tree not valid")
+      subindex(_vp_nearest(points, tree, p, k),0);
+
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

vnf.scad

@@ -221,8 +221,8 @@ function vnf_triangulate(vnf) =
 //   triangles.  The default style is an arbitrary, systematic subdivision in the same direction.  The "alt" style
 //   is the uniform subdivision in the other (alternate) direction.  The "min_edge" style picks the shorter edge to
 //   subdivide for each quadrilateral, so the division may not be uniform across the shape.  The "quincunx" style
-//   adds a vertex in the center of each quadrilateral and creates four triangles, and the "convex" style
+//   adds a vertex in the center of each quadrilateral and creates four triangles, and the "convex" and "concave" styles
-//   chooses the locally convex subdivision.  
+//   chooses the locally convex/concave subdivision.  
 // Arguments:
 //   points = A list of vertices to divide into columns and rows.
 //   caps = If true, add endcap faces to the first AND last rows.
@@ -231,7 +231,7 @@ function vnf_triangulate(vnf) =
 //   col_wrap = If true, add faces to connect the last column to the first.
 //   row_wrap = If true, add faces to connect the last row to the first.
 //   reverse = If true, reverse all face normals.
-//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", "min_edge", "quincunx", and "convex".
+//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", "min_edge", "quincunx","convex" and "concave".
 //   vnf = If given, add all the vertices and faces to this existing VNF structure.
 // Example(3D):
 //   vnf = vnf_vertex_array(
@@ -297,7 +297,7 @@ function vnf_vertex_array(
 ) = 
     assert(!(any([caps,cap1,cap2]) && !col_wrap), "col_wrap must be true if caps are requested")
     assert(!(any([caps,cap1,cap2]) && row_wrap), "Cannot combine caps with row_wrap")
-    assert(in_list(style,["default","alt","quincunx", "convex","min_edge"]))
+    assert(in_list(style,["default","alt","quincunx", "convex","concave", "min_edge"]))
     assert(is_consistent(points), "Non-rectangular or invalid point array")
     let(
         pts = flatten(points),
@@ -358,6 +358,14 @@ function vnf_vertex_array(
                                                                       : [[i1,i3,i2],[i1,i4,i3]]
                               )
                               convexfaces
+                          : style=="concave"?  
+                              let(   // Find normal for 3 of the points.  Is the other point above or below?
+                                  n = (reverse?-1:1)*cross(pts[i2]-pts[i1],pts[i3]-pts[i1]),
+                                  concavefaces = n==0 ? [[i1,i4,i3]]
+                                              : n*pts[i4] <= n*pts[i1] ? [[i1,i4,i2],[i2,i4,i3]]
+                                                                      : [[i1,i3,i2],[i1,i4,i3]]
+                              )
+                              concavefaces
                           : [[i1,i3,i2],[i1,i4,i3]],
                        // remove degenerate faces 
                        culled_faces= [for(face=faces)