diff --git a/math.scad b/math.scad index 408afe2..73659ca 100644 --- a/math.scad +++ b/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); diff --git a/vectors.scad b/vectors.scad index f8bc2f5..b110fd7 100644 --- a/vectors.scad +++ b/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, ) +// 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 diff --git a/vnf.scad b/vnf.scad index 50950fd..dc5ffa0 100644 --- a/vnf.scad +++ b/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 -// chooses the locally convex subdivision. +// adds a vertex in the center of each quadrilateral and creates four triangles, and the "convex" and "concave" styles +// 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)