Developmental superskin() wrapper and various supporting stuff.

Code in a rough state at this point.

skin.scad

@@ -293,4 +293,307 @@ function skin(profiles, closed=false, caps=true, method="uniform") =
 	) vnfout;
 
 
+
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+//
+// Developmental skin wrapper, called superskin() for now, but this
+// is not meant to be the final name.
+
+// Function&Module: superskin()
+// Usage: As module:
+//   skin(profiles, [slices], [samples|refine], [method], [smethod], [caps], [closed], [z]);
+// Usage: As function:
+//   vnf = skin(profiles, [slices], [samples|refine], [method], [smethod], [caps], [closed], [z]);
+// Description:
+//   Given a list of two ore more path `profiles` in 3d space, produces faces to skin a surface between
+//   the profiles.  Optionally the first and last profiles can have endcaps, or the first and last profiles
+//   can be connected together.  Each profile should be roughly planar, but some variation is allowed.
+//   Each profile must rotate in the same clockwise direction.  If called as a function, returns a
+//   [VNF structure](vnf.scad) like `[VERTICES, FACES]`.  If called as a module, creates a polyhedron
+//    of the skined profiles.
+//   
+//   The profiles can be specified either as a list of 3d curves or they can be specified as
+//   2d curves with heights given in the `z` parameter.  
+//   
+//   For this operation to be well-defined, we the profiles must all have the same vertex count and
+//   we must assume that profiles are aligned so that vertex `i` links to vertex `i` on all polygons.  
+//   Many interesting cases do not comply with this restriction.  To handle these cases, you can
+//   specify various matching methods (listed below).  You can also adjust non-matching profiles
+//   by either resampling them using `subdivide_path` or by duplicating vertices using
+//   `repeat_entries`.  It is OK to pass a profile that has the same vertex repeated, such as
+//   a square with 5 points (two of which are identical), so that it can match up to a pentagon.
+//   Such a combination would create a triangular face at the location of the duplicated vertex.
+//   
+//   In order for skinned surfaces to look good it is usually necessary to use a fine sampling of
+//   points on all of the profiles, and a large number of extra interpolated slices between the
+//   profiles that you specify.  The `slices` parameter specifies the number of slices to insert
+//   between each pair of profiles, either a scalar to insert the same number everywhere, or a vector
+//   to insert a different number between each pair.  To resample the profiles you can specify the
+//   number of samples at each profiles with the `samples` argument or you can use `refine`.  The
+//   `refine` parameter specifies a multiplication factor relative to the largest profile, so
+//   if refine is 10 and the largest profile has length 6 then you will get a total of 60 points,
+//   or 10 points per side of the longest profile.  The default is `samples` equal to the size
+//   of the largest profile, which will do nothing if all profiles are the same size.  
+//   
+//   Two methods are available for resampling, `"length"` and `"segment"`.  Specify them using
+//   the `smethod` argument.  The length resampling method resamples proportional to length.
+//   The segment method divides each segment of a profile into the same number of points.
+//   A uniform division may be impossible, in which case the code computes an approximation.
+//   See `subdivide_path` for more details.  
+//   
+//   You can choose from four methods for specifying alignment for incomensurate profiles.
+//   The available methods are `"distance"`, `"tangent"`, `"uniform"` and `"align"`.
+//   The "distance" method finds the global minimum distance method for connecting two
+//   profiles.  This algorithm generally produces a good result when both profiles have
+//   a small number of vertices.  It is computationally intensive (O(N^3)) and may be
+//   slow on large inputs.  The `"tangent"` method generally produces good results when
+//   connecting a discrete polygon to a convex, finely sampled curve.  It works by finding
+//   a plane that passed through each edge of the polygon that is tangent to
+//   the curve.  The `"uniform"` method simply connects the vertices, after resampling
+//   if it is required.  The `"align"` method resamples the vertices and then reindexes
+//   to find the shortest distance alignment.  This will result in the faces with the
+//   smallest amount of twist.  The align algorithm has quadratic run time and can be slow
+//   with large profiles.  
+//   
+// Arguments:
+//   profiles = list of 2d or 3d profiles to be skinned.  (If 2d must also give `z`.)
+//   slices = scalar or vector number of slices to insert between each pair of profiles.  Default: 8.
+//   samples = resample each profile to this many points.  If `method` is distance default is undef, otherwise default is the length of longest profile.
+//   refine = resample profiles to this number of points per side.  If `method` is "distance" default is 10, otherwise undef. 
+//   smethod = sampling method, either "length" or "segment".  If `method` is "distance" or tangent default is "segment", otherwise "length".
+//   caps = true to create endcap faces.  Default is true if closed is false.
+//   method = method for aligning and connecting profiles
+//   closed = set to true to connect first and last profile.  Default: false
+//   z = array of height values for each profile if the profiles are 2d
+module superskin(profiles, slices=8, samples, refine, method="uniform", smethod, caps, closed=false, z)
+{
+  	vnf_polyhedron(superskin(profiles, slices, samples, refine, method, smethod, caps, closed, z));
+}        
+
+function superskin(profiles, slices=8, samples, refine, method="uniform", smethod, caps, closed=false, z) = 
+  let(
+    legal_methods = ["uniform","align","distance","tangent"],
+    caps = is_def(caps) ? caps :
+           closed ? false : true,
+    default_refine = 10,  
+    maxsize = list_longest(profiles),
+    samples = echo(at_sample_method=method)is_def(samples) && is_def(refine) ? undef :
+              is_def(samples) ? samples :
+              is_def(refine)  ? maxsize*refine :
+              method=="distance" ? maxsize*default_refine :
+                                   maxsize,
+    
+    methodok = is_list(method) || in_list(method, legal_methods),
+    methodlistok = is_list(method) ? [for(i=[0:len(method)-1]) if (!in_list(method[i], legal_methods)) i] : [],
+    method = is_string(method) ? replist(method, len(profiles)+ (closed?1:0)) : method,
+    smethod = is_def(smethod)? smethod :
+              all([for(m=method) m=="distance" || m=="tangent"]) ? "segment" : "length"
+    )
+
+  assert(methodok,str("method must be one of ",legal_methods,". Got ",method))
+  assert(methodlistok==[], str("method list contains invalid method at ",methodlistok))
+  assert(!closed || !caps, "Cannot make closed shape with caps")
+  assert(is_def(samples),"Specify only one of `refine` and `samples`")
+  assert(samples>=maxsize,str("Requested number of samples ",samples," is smaller than size of largest profile, ",maxsize))
+  let(
+    profile_dim=array_dim(profiles,2),
+    profiles_ok = (profile_dim==2 && is_list(z) && len(z)==len(profiles)) || profile_dim==3
+  )
+  assert(profiles_ok,"Profiles must all be 3d or must all be 2d, with matching length z parameter.")
+  assert(is_undef(z) || profile_dim==2, "Do not specify z with 3d profiles")
+  assert(profile_dim==3 || len(z)==len(profiles),"Length of z does not match length of profiles.")
+  let(
+    profiles = profile_dim==3 ? profiles :
+               [for(i=[0:len(profiles)-1]) path3d(profiles[i], z[i])],
+    full_list =
+      [for(i=[0:len(profiles)-(closed?1:2)])
+        let(
+          pair = 
+            method[i]=="distance" ? minimum_distance_match(profiles[i],select(profiles,i+1)) :
+            method[i]=="tangent" ? tangent_align(profiles[i],select(profiles,i+1)) :
+            /*method[i]=="align" || method[i]=="uniform" ?*/
+               let( p1 = subdivide_path(profiles[i],samples, method=smethod),
+                    p2 = subdivide_path(select(profiles,i+1),samples, method=smethod)
+               ) (method[i]=="uniform" ? [p1,p2] : [p1, reindex_polygon(p1, p2)])
+          )
+          each interp_and_slice(pair,slices, samples, submethod=smethod)]
+  )
+  skin(full_list, method="uniform");
+
+
+
+
+// plist is list of polygons, N is list or value for number of slices to insert
+// numpoints can be "max", "lcm" or a number
+function interp_and_slice(plist, N, numpoints="max", align=false,submethod="length") =
+  let(
+    maxsize = list_longest(plist),
+    numpoints = numpoints == "max" ? maxsize :
+                numpoints == "lcm" ? lcmlist([for(p=plist) len(p)]) :
+                is_num(numpoints) ? round(numpoints) : undef
+  )
+  assert(is_def(numpoints), "Parameter numpoints must be \"max\", \"lcm\" or a positive number")
+  assert(numpoints>=maxsize, "Number of points requested is smaller than largest profile")
+  let(fixpoly = [for(poly=plist) subdivide_path(poly, numpoints,method=submethod)])
+  add_slices(fixpoly, N);
+  
+
+
+
+function add_slices(plist,N) =
+  assert(is_num(N) || is_list(N))
+  let(listok = !is_list(N) || len(N)==len(plist)-1)
+  assert(listok, "Input N to add_slices is a list with the wrong length")
+  let(
+    count = is_num(N) ? replist(N,len(plist)-1) : N,
+    slicelist = [for (i=[0:len(plist)-2])
+      each [for(j = [0:count[i]]) lerp(plist[i],plist[i+1],j/(count[i]+1))]
+    ]
+  )
+  concat(slicelist, [plist[len(plist)-1]]);
+
+
+
+// Function: unique_count()
+// Usage:
+//   unique_count(arr);
+// Description:
+//   Returns `[sorted,counts]` where `sorted` is a sorted list of the unique items in `arr` and `counts` is a list such 
+//   that `count[i]` gives the number of times that `sorted[i]` appears in `arr`.  
+// Arguments:
+//   arr = The list to analyze. 
+function unique_count(arr) =
+	assert(is_list(arr)||is_string(list))
+	len(arr)==0 ? [[],[]] :
+        len(arr)==1 ? [arr,[1]] :
+        _unique_count(sort(arr), ulist=[], counts=[], ind=1, curtot=1);
+
+function _unique_count(arr, ulist, counts, ind, curtot) = 
+     ind == len(arr)+1 ? [ulist, counts] :
+     ind==len(arr) || arr[ind] != arr[ind-1] ? _unique_count(arr,concat(ulist,[arr[ind-1]]), concat(counts,[curtot]),ind+1,1) :
+     _unique_count(arr,ulist,counts,ind+1,curtot+1);
+
+///////////////////////////////////////////////////////
+//
+// Given inputs of a small polygon (`small`) and a larger polygon (`big`), computes an onto mapping of
+// the the vertices of `big` onto `small` that minimizes the sum of the distances between every matched
+// pair of vertices.  The algorithm uses quadratic programming to calculate the optimal mapping under
+// the assumption that big[0]->small[0] and big[len(big)-1] does NOT map to small[0].  We then
+// rotate through all the possible indexings of `big`.  The theoretical run time is quadratic
+// in len(big) and linear in len(small).  
+//
+// The top level function, nbest_dmatch() cycles through all the of the indexings of `big`, computes
+// all of the optimal values, and chooses the overall best result.  It then interprets the result to
+// produce the index mapping.  The function _qp_extract_map() threads back through the quadratic programming
+// array to identify the actual mapping.
+// 
+// The function _qp_distance_array builds up the rows of the quadratic programming matrix with reference
+// to the previous rows, where `tdist` holds the total distance for a given mapping, and `map`
+// holds the information about which path was optimal for each position.
+//
+// The function _qp_distance_row constructs each row of the quadratic programming matrix.  Note that
+// in this problem we can delete entries from `big` but we cannot insert.  This means we can only
+// move to the right, or diagonally, and not down.  This in turn means that only a portion of the
+// quadratic programming matrix is reachable, so we fill in the unreachable lefthand triangular portion
+// with zeros and we just don't compute the righthand portion (meaning that each row of the output
+// has a different length).
+
+// This function builds up the quadratic programming distance array where each entry in the
+// array gives the optimal distance for aligning the corresponding subparts of the two inputs.
+// When the array is fully populated, the bottom right corner gives the minimum distance
+// for matching the full input lists.  The `map` array contains a 0 when the optimal value came from
+// the left (a "deletion") which means you match the next vertex in `big` with the previous, already
+// used vertex of `small`, or a 1 when the optimal value came from the diagonal, which means you
+// match the next vertex of `big` with the next vertex of `small`.
+//
+// Return value is [min_distance, map], where map is the array that is used to extract the actual
+// vertex map.  
+function _qp_distance_array(small, big, small_ind=0, tdist=[], map=[]) =
+   let(
+       N = len(small),
+       M = len(big)
+   )  
+   small_ind == N ? [tdist[N-1][M-1], map] :
+   let(
+     row_results = small_ind == 0 ? [cumsum([for(i=[0:M-N+1]) norm(big[i]-small[0])]), replist(0,M-N+1)] :
+                   _qp_distance_row(small, big, small_ind, small_ind, tdist, replist(0,small_ind), replist(0, small_ind))
+   )
+   _qp_distance_array(small, big, small_ind+1, concat(tdist, [row_results[0]]), concat(map, [row_results[1]]));
+
+
+function _qp_distance_row(small,big,small_ind, big_ind, tdist, newrow, maprow) =
+    big_ind == len(big)-len(small) + small_ind + 1 ? [newrow,maprow] :
+    _qp_distance_row(small,big, small_ind, big_ind+1, tdist,
+                concat(newrow, [norm(small[small_ind]-big[big_ind]) +
+                                (small_ind==big_ind ? tdist[small_ind-1][big_ind-1] : min(tdist[small_ind-1][big_ind-1],newrow[big_ind-1]))]), 
+                concat(maprow, [small_ind!=big_ind && newrow[big_ind-1] < tdist[small_ind-1][big_ind-1] ? 0 : 1]));
+
+
+function _qp_extract_map(map,i,j,result) =
+  is_undef(i) ? _qp_extract_map(map,len(map)-1,len(select(map,-1))-1,[]) :
+  i==0 && j==0 ? concat([0], result) :
+  _qp_extract_map(map,i-map[i][j],j-1,concat([i],result));
+
+
+function minimum_distance_match(poly1,poly2) =
+   let(
+      swap = len(poly1)>len(poly2),
+      big = swap ? poly1 : poly2,
+      small = swap ? poly2 : poly1,
+      matchres = [for(i=[0:len(big)-1]) _qp_distance_array(small,polygon_shift(big,i))],
+      best = min_index(subindex(matchres,0)),
+      newbig = polygon_shift(big,best),
+      newsmall = repeat_entries(small,unique_count(_qp_extract_map(matchres[best][1]))[1])
+      )
+      swap ? [newbig, newsmall] : [newsmall,newbig];
+
+
+function tangent_align(poly1, poly2) =
+    let(
+        swap = len(poly1)>len(poly2),
+        big = swap ? poly1 : poly2,
+        small = swap ? poly2 : poly1,
+        cutpts = [for(i=[0:len(small)-1]) find_one_tangent(big, select(small,i,i+1))],
+        d=echo(cutpts = cutpts),
+        shift = select(cutpts,-1)+1, 
+        newbig = polygon_shift(big, shift),
+        repeat_counts = [for(i=[0:len(small)-1]) posmod(cutpts[i]-select(cutpts,i-1),len(big))],
+        newsmall = repeat_entries(small,repeat_counts)
+      )
+      assert(len(newsmall)==len(newbig), "Tangent alignment failed, probably because of insufficient points or a concave curve")
+      swap ? [newbig, newsmall] : [newsmall, newbig];
+
+
+function find_one_tangent(curve, edge, closed=true) =
+  let(
+   angles = 
+   [for(i=[0:len(curve)-(closed?1:2)])
+     let( 
+       plane = plane3pt( edge[0], edge[1], curve[i]),
+       tangent = [curve[i], select(curve,i+1)]
+       )
+   plane_line_angle(plane,tangent)],
+   zcross = [for(i=[0:len(curve)-(closed?1:2)]) if (sign(angles[i]) != sign(select(angles,i+1))) i],
+   d = [for(i=zcross) distance_from_line(edge, curve[i])]
+    )
+   zcross[min_index(d)];//zcross;
+
+
+function plane_line_angle(plane, line) =
+   let(
+        vect = line[1]-line[0],
+        zplane = select(plane,0,2),
+        sin_angle = vect*zplane/norm(zplane)/norm(vect)
+        )
+   asin(constrain(sin_angle,-1,1));
+
+
+
 // vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+
+
+