diff --git a/arrays.scad b/arrays.scad
index 74d2ea90..728cc47 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -53,6 +53,150 @@ function _same_type(a,b, depth) =
           && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); 
   
 
+
+
+// Function: list_shortest()
+// Usage:
+//   llen = list_shortest(array);
+// Topics: List Handling
+// See Also: list_longest()
+// Description:
+//   Returns the length of the shortest sublist in a list of lists.
+// Arguments:
+//   array = A list of lists.
+// Example:
+//   slen = list_shortest([[3,4,5],[6,7,8,9]]);  // Returns: 3
+function list_shortest(array) =
+    assert(is_list(array), "Invalid input." )
+    min([for (v = array) len(v)]);
+
+
+// Function: list_longest()
+// Usage:
+//   llen = list_longest(array);
+// Topics: List Handling
+// See Also: list_shortest()
+// Description:
+//   Returns the length of the longest sublist in a list of lists.
+// Arguments:
+//   array = A list of lists.
+// Example:
+//   llen = list_longest([[3,4,5],[6,7,8,9]]);  // Returns: 4
+function list_longest(array) =
+    assert(is_list(array), "Invalid input." )
+    max([for (v = array) len(v)]);
+
+
+
+// Function: in_list()
+// Usage:
+//   bool = in_list(val, list, [idx]);
+// Topics: List Handling
+// Description:
+//   Returns true if value `val` is in list `list`. When `val==NAN` the answer will be false for any list.
+// Arguments:
+//   val = The simple value to search for.
+//   list = The list to search.
+//   idx = If given, searches the given columns for matches for `val`.
+// Example:
+//   a = in_list("bar", ["foo", "bar", "baz"]);  // Returns true.
+//   b = in_list("bee", ["foo", "bar", "baz"]);  // Returns false.
+//   c = in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1);  // Returns true.
+function in_list(val,list,idx) = 
+    assert( is_list(list) && (is_undef(idx) || is_finite(idx)),
+            "Invalid input." )
+    let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
+    s==[] || s==[[]] ? false
+    : is_undef(idx) ? val==list[s] 
+    : val==list[s][idx];
+
+
+// Function: find_first_match()
+// Topics: List Handling
+// See Also: in_list()
+// Usage:
+//   idx = find_first_match(val, list, [start=], [eps=]);
+//   indices = find_first_match(val, list, all=true, [start=], [eps=]);
+// Description:
+//   Finds the first item in `list` that matches `val`, returning the index.
+// Arguments:
+//   val = The value to search for.  If given a function literal of signature `function (x)`, uses that function to check list items.  Returns true for a match.
+//   list = The list to search through.
+//   ---
+//   start = The index to start searching from.
+//   all = If true, returns a list of all matching item indices.
+//   eps = The maximum allowed floating point rounding error for numeric comparisons.
+function find_first_match(val, list, start=0, all=false, eps=EPSILON) =
+    all? [
+        for (i=[start:1:len(list)-1])
+        if (
+            (!is_func(val) && approx(val, list[i], eps=eps)) ||
+            (is_func(val) && val(list[i]))
+        ) i
+    ] :
+    __find_first_match(val, list, eps=eps, i=start);
+
+function __find_first_match(val, list, eps, i=0) =
+    i >= len(list)? undef :
+    (
+        (!is_func(val) && approx(val, list[i], eps=eps)) ||
+        (is_func(val) && val(list[i]))
+    )? i : __find_first_match(val, list, eps=eps, i=i+1);
+
+
+// Function: list_increasing()
+// Usage:
+//    bool = list_increasing(list);
+// Topics: List Handling
+// See Also: max_index(), min_index(), list_decreasing()
+// Description:
+//   Returns true if the list is (non-strictly) increasing
+// Example:
+//   a = list_increasing([1,2,3,4]);  // Returns: true
+//   b = list_increasing([1,3,2,4]);  // Returns: false
+//   c = list_increasing([4,3,2,1]);  // Returns: false
+function list_increasing(list) =
+    assert(is_list(list)||is_string(list))
+    len([for (p=pair(list)) if(p.x>p.y) true])==0;
+
+
+// Function: list_decreasing()
+// Usage:
+//   bool = list_decreasing(list);
+// Topics: List Handling
+// See Also: max_index(), min_index(), list_increasing()
+// Description:
+//   Returns true if the list is (non-strictly) decreasing
+// Example:
+//   a = list_decreasing([1,2,3,4]);  // Returns: false
+//   b = list_decreasing([4,2,3,1]);  // Returns: false
+//   c = list_decreasing([4,3,2,1]);  // Returns: true
+function list_decreasing(list) =
+    assert(is_list(list)||is_string(list))
+    len([for (p=pair(list)) if(p.x<p.y) true])==0;
+
+
+// Function: add_scalar()
+// Usage:  
+//   v = add_scalar(v, s);
+// Topics: List Handling
+// Description:
+//   Given a list and a scalar, returns the list with the scalar added to each item in it.
+//   If given a list of arrays, recursively adds the scalar to the each array.
+// Arguments:
+//   v = The initial array.
+//   s = A scalar value to add to every item in the array.
+// Example:
+//   a = add_scalar([1,2,3],3);            // Returns: [4,5,6]
+//   b = add_scalar([[1,2,3],[3,4,5]],3);  // Returns: [[4,5,6],[6,7,8]]
+function add_scalar(v,s) = 
+    is_finite(s) ? [for (x=v) is_list(x)? add_scalar(x,s) : is_finite(x) ? x+s: x] : v;
+
+
+
+
+// Section: List Indexing
+
 // Function: select()
 // Topics: List Handling
 // Description:
@@ -70,7 +214,7 @@ function _same_type(a,b, depth) =
 //   list = The list to get the portion of.
 //   start = Either the index of the first item or an index range or a list of indices.
 //   end = The index of the last item when `start` is a number. When `start` is a list or a range, `end` should not be given.
-// See Also: slice(), subindex(), last()
+// See Also: slice(), columns(), last()
 // Example:
 //   l = [3,4,5,6,7,8,9];
 //   a = select(l, 5, 6);   // Returns [8,9]
@@ -111,7 +255,7 @@ function select(list, start, end) =
 //   list = The list to get the slice of.
 //   s = The index of the first item to return.
 //   e = The index of the last item to return.
-// See Also: select(), subindex(), last()
+// See Also: select(), columns(), last()
 // Example:
 //   a = slice([3,4,5,6,7,8,9], 3, 5);   // Returns [6,7,8]
 //   b = slice([3,4,5,6,7,8,9], 2, -1);  // Returns [5,6,7,8,9]
@@ -136,7 +280,7 @@ function slice(list,s=0,e=-1) =
 // Usage:
 //   item = last(list);
 // Topics: List Handling
-// See Also: select(), slice(), subindex()
+// See Also: select(), slice(), columns()
 // Description:
 //   Returns the last element of a list, or undef if empty.
 // Arguments:
@@ -203,6 +347,32 @@ function list_tail(list, from=1) =
    list;
 
 
+
+// Function: bselect()
+// Usage:
+//   array = bselect(array, index);
+// Topics: List Handling
+// See Also: list_bset()
+// Description:
+//   Returns the items in `array` whose matching element in `index` is true.
+// Arguments:
+//   array = Initial list to extract items from.
+//   index = List of booleans.
+// Example:
+//   a = bselect([3,4,5,6,7], [false,true,true,false,true]);  // Returns: [4,5,7]
+function bselect(array,index) =
+    assert(is_list(array)||is_string(array), "Improper array." )
+    assert(is_list(index) && len(index)>=len(array) , "Improper index list." )
+    is_string(array)? str_join(bselect( [for (x=array) x], index)) :
+    [for(i=[0:len(array)-1]) if (index[i]) array[i]];
+
+
+
+
+
+// Section: List Construction
+
+
 // Function: list()
 // Topics: List Handling, Type Conversion
 // Usage:
@@ -246,155 +416,6 @@ function force_list(value, n=1, fill) =
     is_undef(fill)? [for (i=[1:1:n]) value] : [value, for (i=[2:1:n]) fill];
 
 
-// Function: add_scalar()
-// Usage:  
-//   v = add_scalar(v, s);
-// Topics: List Handling
-// Description:
-//   Given a list and a scalar, returns the list with the scalar added to each item in it.
-//   If given a list of arrays, recursively adds the scalar to the each array.
-// Arguments:
-//   v = The initial array.
-//   s = A scalar value to add to every item in the array.
-// Example:
-//   a = add_scalar([1,2,3],3);            // Returns: [4,5,6]
-//   b = add_scalar([[1,2,3],[3,4,5]],3);  // Returns: [[4,5,6],[6,7,8]]
-function add_scalar(v,s) = 
-    is_finite(s) ? [for (x=v) is_list(x)? add_scalar(x,s) : is_finite(x) ? x+s: x] : v;
-
-
-// Function: in_list()
-// Usage:
-//   bool = in_list(val, list, [idx]);
-// Topics: List Handling
-// Description:
-//   Returns true if value `val` is in list `list`. When `val==NAN` the answer will be false for any list.
-// Arguments:
-//   val = The simple value to search for.
-//   list = The list to search.
-//   idx = If given, searches the given subindex for matches for `val`.
-// Example:
-//   a = in_list("bar", ["foo", "bar", "baz"]);  // Returns true.
-//   b = in_list("bee", ["foo", "bar", "baz"]);  // Returns false.
-//   c = in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1);  // Returns true.
-function in_list(val,list,idx) = 
-    assert( is_list(list) && (is_undef(idx) || is_finite(idx)),
-            "Invalid input." )
-    let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
-    s==[] || s==[[]] ? false
-    : is_undef(idx) ? val==list[s] 
-    : val==list[s][idx];
-
-
-// Function: find_first_match()
-// Topics: List Handling
-// See Also: in_list()
-// Usage:
-//   idx = find_first_match(val, list, [start=], [eps=]);
-//   indices = find_first_match(val, list, all=true, [start=], [eps=]);
-// Description:
-//   Finds the first item in `list` that matches `val`, returning the index.
-// Arguments:
-//   val = The value to search for.  If given a function literal of signature `function (x)`, uses that function to check list items.  Returns true for a match.
-//   list = The list to search through.
-//   ---
-//   start = The index to start searching from.
-//   all = If true, returns a list of all matching item indices.
-//   eps = The maximum allowed floating point rounding error for numeric comparisons.
-function find_first_match(val, list, start=0, all=false, eps=EPSILON) =
-    all? [
-        for (i=[start:1:len(list)-1])
-        if (
-            (!is_func(val) && approx(val, list[i], eps=eps)) ||
-            (is_func(val) && val(list[i]))
-        ) i
-    ] :
-    __find_first_match(val, list, eps=eps, i=start);
-
-function __find_first_match(val, list, eps, i=0) =
-    i >= len(list)? undef :
-    (
-        (!is_func(val) && approx(val, list[i], eps=eps)) ||
-        (is_func(val) && val(list[i]))
-    )? i : __find_first_match(val, list, eps=eps, i=i+1);
-
-
-// Function: min_index()
-// Usage:
-//   idx = min_index(vals);
-//   idxlist = min_index(vals, all=true);
-// Topics: List Handling
-// See Also: max_index(), list_increasing(), list_decreasing()
-// Description:
-//   Returns the index of the first occurrence of the minimum value in the given list. 
-//   If `all` is true then returns a list of all indices where the minimum value occurs.
-// Arguments:
-//   vals = vector of values
-//   all = set to true to return indices of all occurences of the minimum.  Default: false
-// Example:
-//   a = min_index([5,3,9,6,2,7,8,2,1]); // Returns: 8
-//   b = min_index([5,3,9,6,2,7,8,2,7],all=true); // Returns: [4,7]
-function min_index(vals, all=false) =
-    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
-    all ? search(min(vals),vals,0) : search(min(vals), vals)[0];
-
-
-// Function: max_index()
-// Usage:
-//   idx = max_index(vals);
-//   idxlist = max_index(vals, all=true);
-// Topics: List Handling
-// See Also: min_index(), list_increasing(), list_decreasing()
-// Description:
-//   Returns the index of the first occurrence of the maximum value in the given list. 
-//   If `all` is true then returns a list of all indices where the maximum value occurs.
-// Arguments:
-//   vals = vector of values
-//   all = set to true to return indices of all occurences of the maximum.  Default: false
-// Example:
-//   max_index([5,3,9,6,2,7,8,9,1]); // Returns: 2
-//   max_index([5,3,9,6,2,7,8,9,1],all=true); // Returns: [2,7]
-function max_index(vals, all=false) =
-    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
-    all ? search(max(vals),vals,0) : search(max(vals), vals)[0];
-
-
-// Function: list_increasing()
-// Usage:
-//    bool = list_increasing(list);
-// Topics: List Handling
-// See Also: max_index(), min_index(), list_decreasing()
-// Description:
-//   Returns true if the list is (non-strictly) increasing
-// Example:
-//   a = list_increasing([1,2,3,4]);  // Returns: true
-//   b = list_increasing([1,3,2,4]);  // Returns: false
-//   c = list_increasing([4,3,2,1]);  // Returns: false
-function list_increasing(list) =
-    assert(is_list(list)||is_string(list))
-    len([for (p=pair(list)) if(p.x>p.y) true])==0;
-
-
-// Function: list_decreasing()
-// Usage:
-//   bool = list_decreasing(list);
-// Topics: List Handling
-// See Also: max_index(), min_index(), list_increasing()
-// Description:
-//   Returns true if the list is (non-strictly) decreasing
-// Example:
-//   a = list_decreasing([1,2,3,4]);  // Returns: false
-//   b = list_decreasing([4,2,3,1]);  // Returns: false
-//   c = list_decreasing([4,3,2,1]);  // Returns: true
-function list_decreasing(list) =
-    assert(is_list(list)||is_string(list))
-    len([for (p=pair(list)) if(p.x<p.y) true])==0;
-
-
-
-// Section: Basic List Generation
-
-
 // Function: repeat()
 // Usage:
 //   list = repeat(val, n);
@@ -442,7 +463,38 @@ function count(n,s=0,step=1,reverse=false) = let(n=is_list(n) ? len(n) : n)
 
 
 
-// Section: List Manipulation
+// Function: list_bset()
+// Usage:
+//   arr = list_bset(indexset, valuelist, [dflt]);
+// Topics: List Handling
+// See Also: bselect()
+// Description:
+//   Opposite of `bselect()`.  Returns a list the same length as `indexlist`, where each item will
+//   either be 0 if the corresponding item in `indexset` is false, or the next sequential value
+//   from `valuelist` if the item is true.  The number of `true` values in `indexset` must be equal 
+//   to the length of `valuelist`.
+// Arguments:
+//   indexset = A list of boolean values.
+//   valuelist = The list of values to set into the returned list.
+//   dflt = Default value to store when the indexset item is false.
+// Example:
+//   a = list_bset([false,true,false,true,false], [3,4]);  // Returns: [0,3,0,4,0]
+//   b = list_bset([false,true,false,true,false], [3,4], dflt=1);  // Returns: [1,3,1,4,1]
+function list_bset(indexset, valuelist, dflt=0) =
+    assert(is_list(indexset), "The index set is not a list." )
+    assert(is_list(valuelist), "The `valuelist` is not a list." )
+    let( trueind = search([true], indexset,0)[0] )
+    assert( !(len(trueind)>len(valuelist)), str("List `valuelist` too short; its length should be ",len(trueind)) )
+    assert( !(len(trueind)<len(valuelist)), str("List `valuelist` too long; its length should be ",len(trueind)) )
+    concat(
+        list_set([],trueind, valuelist, dflt=dflt),    // Fill in all of the values
+        repeat(dflt,len(indexset)-max(trueind)-1)  // Add trailing values so length matches indexset
+    );
+
+
+
+
+// Section: List Modification
 
 // Function: reverse()
 // Usage:
@@ -756,87 +808,9 @@ function list_remove_values(list,values=[],all=false) =
     ) list_remove(list,uidxs);
 
 
-// Function: bselect()
-// Usage:
-//   array = bselect(array, index);
-// Topics: List Handling
-// See Also: list_bset()
-// Description:
-//   Returns the items in `array` whose matching element in `index` is true.
-// Arguments:
-//   array = Initial list to extract items from.
-//   index = List of booleans.
-// Example:
-//   a = bselect([3,4,5,6,7], [false,true,true,false,true]);  // Returns: [4,5,7]
-function bselect(array,index) =
-    assert(is_list(array)||is_string(array), "Improper array." )
-    assert(is_list(index) && len(index)>=len(array) , "Improper index list." )
-    is_string(array)? str_join(bselect( [for (x=array) x], index)) :
-    [for(i=[0:len(array)-1]) if (index[i]) array[i]];
-
-
-// Function: list_bset()
-// Usage:
-//   arr = list_bset(indexset, valuelist, [dflt]);
-// Topics: List Handling
-// See Also: bselect()
-// Description:
-//   Opposite of `bselect()`.  Returns a list the same length as `indexlist`, where each item will
-//   either be 0 if the corresponding item in `indexset` is false, or the next sequential value
-//   from `valuelist` if the item is true.  The number of `true` values in `indexset` must be equal 
-//   to the length of `valuelist`.
-// Arguments:
-//   indexset = A list of boolean values.
-//   valuelist = The list of values to set into the returned list.
-//   dflt = Default value to store when the indexset item is false.
-// Example:
-//   a = list_bset([false,true,false,true,false], [3,4]);  // Returns: [0,3,0,4,0]
-//   b = list_bset([false,true,false,true,false], [3,4], dflt=1);  // Returns: [1,3,1,4,1]
-function list_bset(indexset, valuelist, dflt=0) =
-    assert(is_list(indexset), "The index set is not a list." )
-    assert(is_list(valuelist), "The `valuelist` is not a list." )
-    let( trueind = search([true], indexset,0)[0] )
-    assert( !(len(trueind)>len(valuelist)), str("List `valuelist` too short; its length should be ",len(trueind)) )
-    assert( !(len(trueind)<len(valuelist)), str("List `valuelist` too long; its length should be ",len(trueind)) )
-    concat(
-        list_set([],trueind, valuelist, dflt=dflt),    // Fill in all of the values
-        repeat(dflt,len(indexset)-max(trueind)-1)  // Add trailing values so length matches indexset
-    );
-
 
 // Section: List Length Manipulation
 
-// Function: list_shortest()
-// Usage:
-//   llen = list_shortest(array);
-// Topics: List Handling
-// See Also: list_longest()
-// Description:
-//   Returns the length of the shortest sublist in a list of lists.
-// Arguments:
-//   array = A list of lists.
-// Example:
-//   slen = list_shortest([[3,4,5],[6,7,8,9]]);  // Returns: 3
-function list_shortest(array) =
-    assert(is_list(array), "Invalid input." )
-    min([for (v = array) len(v)]);
-
-
-// Function: list_longest()
-// Usage:
-//   llen = list_longest(array);
-// Topics: List Handling
-// See Also: list_shortest()
-// Description:
-//   Returns the length of the longest sublist in a list of lists.
-// Arguments:
-//   array = A list of lists.
-// Example:
-//   llen = list_longest([[3,4,5],[6,7,8,9]]);  // Returns: 4
-function list_longest(array) =
-    assert(is_list(array), "Invalid input." )
-    max([for (v = array) len(v)]);
-
 
 // Function: list_pad()
 // Usage:
@@ -1253,8 +1227,49 @@ function unique_count(list) =
         )
         [ select(list,ind), deltas( concat(ind,[len(list)]) ) ];
 
-    
-// Section: List Iteration Helpers
+
+// Function: group_data()
+// Usage:
+//   groupings = group_data(groups, values);
+// Topics: Array Handling
+// See Also: zip(), zip_long(), array_group()
+// Description:
+//   Given a list of integer group numbers, and an equal-length list of values,
+//   returns a list of groups with the values sorted into the corresponding groups.
+//   Ie: if you have a groups index list of [2,3,2] and values of ["A","B","C"], then
+//   the values "A" and "C" will be put in group 2, and "B" will be in group 3.
+//   Groups that have no values grouped into them will be an empty list.  So the
+//   above would return [[], [], ["A","C"], ["B"]]
+// Arguments:
+//   groups = A list of integer group index numbers.
+//   values = A list of values to sort into groups.
+// Example:
+//   groups = group_data([1,2,0], ["A","B","C"]);  // Returns [["B"],["C"],["A"]]
+// Example:
+//   groups = group_data([1,3,1], ["A","B","C"]);  // Returns [[],["A","C"],[],["B"]]
+function group_data(groups, values) =
+    assert(all_integer(groups) && all_nonnegative(groups))
+    assert(is_list(values))
+    assert(len(groups)==len(values),
+           "The groups and values arguments should be lists of matching length.")
+    let( sorted = _group_sort_by_index(zip(groups,values),0) )
+    // retrieve values and insert []
+    [
+        for (i = idx(sorted))
+        let(
+            a  = i==0? 0 : sorted[i-1][0][0]+1,
+            g0 = sorted[i]
+        )
+        each [
+            for (j = [a:1:g0[0][0]-1]) [],
+            [for (g1 = g0) g1[1]]
+        ]
+    ];
+
+
+
+
+// Section: Iteration Helpers
 
 // Function: idx()
 // Usage:
@@ -1293,7 +1308,7 @@ function idx(list, s=0, e=-1, step=1) =
 //   Something like: `[[0,l[0]], [1,l[1]], [2,l[2]], ...]`
 // Arguments:
 //   l = List to enumerate.
-//   idx = If given, enumerates just the given subindex items of `l`.
+//   idx = If given, enumerates just the given columns items of `l`.
 // Example:
 //   enumerate(["a","b","c"]);  // Returns: [[0,"a"], [1,"b"], [2,"c"]]
 //   enumerate([[88,"a"],[76,"b"],[21,"c"]], idx=1);  // Returns: [[0,"a"], [1,"b"], [2,"c"]]
@@ -1562,9 +1577,9 @@ function set_intersection(a, b) =
 
 // Section: Array Manipulation
 
-// Function: subindex()
+// Function: columns()
 // Usage:
-//   list = subindex(M, idx);
+//   list = columns(M, idx);
 // Topics: Array Handling, List Handling
 // See Also: select(), slice()
 // Description:
@@ -1577,13 +1592,13 @@ function set_intersection(a, b) =
 //   idx = The index, list of indices, or range of indices to fetch.
 // Example:
 //   M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-//   a = subindex(M,2);      // Returns [3, 7, 11, 15]
-//   b = subindex(M,[2]);    // Returns [[3], [7], [11], [15]]
-//   c = subindex(M,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
-//   d = subindex(M,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
+//   a = columns(M,2);      // Returns [3, 7, 11, 15]
+//   b = columns(M,[2]);    // Returns [[3], [7], [11], [15]]
+//   c = columns(M,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
+//   d = columns(M,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
 //   N = [ [1,2], [3], [4,5], [6,7,8] ];
-//   e = subindex(N,[0,1]);  // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
-function subindex(M, idx) =
+//   e = columns(N,[0,1]);  // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
+function columns(M, idx) =
     assert( is_list(M), "The input is not a list." )
     assert( !is_undef(idx) && _valid_idx(idx,0,1/0), "Invalid index input." ) 
     is_finite(idx)
@@ -1595,7 +1610,7 @@ function subindex(M, idx) =
 // Usage:
 //   mat = submatrix(M, idx1, idx2);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), block_matrix(), submatrix_set()
+// See Also: columns(), block_matrix(), submatrix_set()
 // Description:
 //   The input must be a list of lists (a matrix or 2d array).  Returns a submatrix by selecting the rows listed in idx1 and columns listed in idx2.
 // Arguments:
@@ -1628,10 +1643,10 @@ function submatrix(M,idx1,idx2) =
 //   A = hstack(M1, M2, M3)
 //   A = hstack([M1, M2, M3, ...])
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix(), block_matrix()
+// See Also: columns(), submatrix(), block_matrix()
 // Description:
 //   Constructs a matrix by horizontally "stacking" together compatible matrices or vectors.  Vectors are treated as columsn in the stack.
-//   This command is the inverse of subindex.  Note: strings given in vectors are broken apart into lists of characters.  Strings given
+//   This command is the inverse of `columns`.  Note: strings given in vectors are broken apart into lists of characters.  Strings given
 //   in matrices are preserved as strings.  If you need to combine vectors of strings use array_group as shown below to convert the
 //   vector into a column matrix.  Also note that vertical stacking can be done directly with concat.  
 // Arguments:
@@ -1650,7 +1665,7 @@ function submatrix(M,idx1,idx2) =
 //   c = hstack([M,v1,M]);  // Returns [[1, 0, 0, 2, 1, 0, 0],
 //                          //          [0, 1, 0, 3, 0, 1, 0],
 //                          //          [0, 0, 1, 4, 0, 0, 1]]
-//   d = hstack(subindex(M,0), subindex(M,[1 2]));  // Returns M
+//   d = hstack(columns(M,0), columns(M,[1 2]));  // Returns M
 //   strvec = ["one","two"];
 //   strmat = [["three","four"], ["five","six"]];
 //   e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
@@ -1680,7 +1695,7 @@ function hstack(M1, M2, M3) =
 // Usage:
 //    bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix()
+// See Also: columns(), submatrix()
 // Description:
 //    Create a block matrix by supplying a matrix of matrices, which will
 //    be combined into one unified matrix.  Every matrix in one row
@@ -1724,7 +1739,7 @@ function block_matrix(M) =
 // Usage:
 //   mat = diagonal_matrix(diag, [offdiag]);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix()
+// See Also: columns(), submatrix()
 // Description:
 //   Creates a square matrix with the items in the list `diag` on
 //   its diagonal.  The off diagonal entries are set to offdiag,
@@ -1741,7 +1756,7 @@ function diagonal_matrix(diag, offdiag=0) =
 // Usage:
 //   mat = submatrix_set(M, A, [m], [n]);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix()
+// See Also: columns(), submatrix()
 // Description:
 //   Sets a submatrix of M equal to the matrix A.  By default the top left corner of M is set to A, but
 //   you can specify offset coordinates m and n.  If A (as adjusted by m and n) extends beyond the bounds
@@ -1772,7 +1787,7 @@ function submatrix_set(M,A,m=0,n=0) =
 //   Takes a flat array of values, and groups items in sets of `cnt` length.
 //   The opposite of this is `flatten()`.
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix(), hstack(), flatten(), full_flatten()
+// See Also: columns(), submatrix(), hstack(), flatten(), full_flatten()
 // Arguments:
 //   v = The list of items to group.
 //   cnt = The number of items to put in each grouping.  Default:2
@@ -1786,50 +1801,12 @@ function array_group(v, cnt=2, dflt=0) =
     [for (i = [0:cnt:len(v)-1]) [for (j = [0:1:cnt-1]) default(v[i+j], dflt)]];
 
 
-// Function: group_data()
-// Usage:
-//   groupings = group_data(groups, values);
-// Topics: Array Handling
-// See Also: zip(), zip_long(), array_group()
-// Description:
-//   Given a list of integer group numbers, and an equal-length list of values,
-//   returns a list of groups with the values sorted into the corresponding groups.
-//   Ie: if you have a groups index list of [2,3,2] and values of ["A","B","C"], then
-//   the values "A" and "C" will be put in group 2, and "B" will be in group 3.
-//   Groups that have no values grouped into them will be an empty list.  So the
-//   above would return [[], [], ["A","C"], ["B"]]
-// Arguments:
-//   groups = A list of integer group index numbers.
-//   values = A list of values to sort into groups.
-// Example:
-//   groups = group_data([1,2,0], ["A","B","C"]);  // Returns [["B"],["C"],["A"]]
-// Example:
-//   groups = group_data([1,3,1], ["A","B","C"]);  // Returns [[],["A","C"],[],["B"]]
-function group_data(groups, values) =
-    assert(all_integer(groups) && all_nonnegative(groups))
-    assert(is_list(values))
-    assert(len(groups)==len(values),
-           "The groups and values arguments should be lists of matching length.")
-    let( sorted = _group_sort_by_index(zip(groups,values),0) )
-    // retrieve values and insert []
-    [
-        for (i = idx(sorted))
-        let(
-            a  = i==0? 0 : sorted[i-1][0][0]+1,
-            g0 = sorted[i]
-        )
-        each [
-            for (j = [a:1:g0[0][0]-1]) [],
-            [for (g1 = g0) g1[1]]
-        ]
-    ];
-
 
 // Function: flatten()
 // Usage:
 //   list = flatten(l);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix(), hstack(), full_flatten()
+// See Also: columns(), submatrix(), hstack(), full_flatten()
 // Description:
 //   Takes a list of lists and flattens it by one level.
 // Arguments:
@@ -1845,7 +1822,7 @@ function flatten(l) =
 // Usage:
 //   list = full_flatten(l);
 // Topics: Matrices, Array Handling
-// See Also: subindex(), submatrix(), hstack(), flatten()
+// See Also: columns(), submatrix(), hstack(), flatten()
 // Description: 
 //   Collects in a list all elements recursively found in any level of the given list.
 //   The output list is ordered in depth first order.
@@ -1974,6 +1951,8 @@ function transpose(arr, reverse=false) =
            arr;
 
 
+// Section: Matrices
+
 // Function: is_matrix_symmetric()
 // Usage:
 //   b = is_matrix_symmetric(A, [eps])
@@ -1987,7 +1966,6 @@ function is_matrix_symmetric(A,eps=1e-12) =
     approx(A,transpose(A), eps);
 
 
-
 // Function&Module: echo_matrix()
 // Usage:
 //    echo_matrix(M, [description=], [sig=], [eps=]);
diff --git a/attachments.scad b/attachments.scad
index d75464e..7b1f510 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -1594,9 +1594,9 @@ function _find_anchor(anchor, geom) =
             hits = [
                 for (face = faces) let(
                     verts = select(rpts, face),
-                    xs = subindex(verts,0),
-                    ys = subindex(verts,1),
-                    zs = subindex(verts,2)
+                    xs = columns(verts,0),
+                    ys = columns(verts,1),
+                    zs = columns(verts,2)
                 ) if (
                     max(xs) >= -eps &&
                     max(ys) >= -eps &&
@@ -1615,7 +1615,7 @@ function _find_anchor(anchor, geom) =
         )
         assert(len(hits)>0, "Anchor vector does not intersect with the shape.  Attachment failed.")
         let(
-            furthest = max_index(subindex(hits,0)),
+            furthest = max_index(columns(hits,0)),
             dist = hits[furthest][0],
             pos = hits[furthest][2],
             hitnorms = [for (hit = hits) if (approx(hit[0],dist,eps=eps)) hit[1]],
@@ -1640,7 +1640,7 @@ function _find_anchor(anchor, geom) =
         ) vnf==EMPTY_VNF? [anchor, [0,0,0], unit(anchor), 0] :
         let(
             rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), vnf[0]),
-            maxx = max(subindex(rpts,0)),
+            maxx = max(columns(rpts,0)),
             idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
             mm = pointlist_bounds(select(rpts,idxs)),
             avgy = (mm[0].y+mm[1].y)/2,
@@ -1681,7 +1681,7 @@ function _find_anchor(anchor, geom) =
                 )
                 if(!is_undef(isect) && !approx(isect,t[0])) [norm(isect), isect, n2]
             ],
-            maxidx = max_index(subindex(isects,0)),
+            maxidx = max_index(columns(isects,0)),
             isect = isects[maxidx],
             pos = point2d(cp) + isect[1],
             vec = unit(isect[2],[0,1])
@@ -1691,7 +1691,7 @@ function _find_anchor(anchor, geom) =
             path = geom[1],
             anchor = point2d(anchor),
             rpath = rot(from=anchor, to=RIGHT, p=move(point2d(-cp), p=path)),
-            maxx = max(subindex(rpath,0)),
+            maxx = max(columns(rpath,0)),
             idxs = [for (i = idx(rpath)) if (approx(rpath[i].x, maxx)) i],
             miny = min([for (i=idxs) rpath[i].y]),
             maxy = max([for (i=idxs) rpath[i].y]),
@@ -1717,7 +1717,7 @@ function _find_anchor(anchor, geom) =
                 if(!is_undef(isect) && !approx(isect,t[0]))
                 [norm(isect), isect, n2]
             ],
-            maxidx = max_index(subindex(isects,0)),
+            maxidx = max_index(columns(isects,0)),
             isect = isects[maxidx],
             pos = point3d(cp) + point3d(isect[1]) + unit([0,0,anchor.z],CENTER)*l/2,
             xyvec = unit(isect[2],[0,1]),
@@ -1730,7 +1730,7 @@ function _find_anchor(anchor, geom) =
             anchor = point3d(anchor),
             xyanch = point2d(anchor),
             rpath = rot(from=xyanch, to=RIGHT, p=move(point2d(-cp), p=path)),
-            maxx = max(subindex(rpath,0)),
+            maxx = max(columns(rpath,0)),
             idxs = [for (i = idx(rpath)) if (approx(rpath[i].x, maxx)) i],
             ys = [for (i=idxs) rpath[i].y],
             avgy = (min(ys)+max(ys))/2,
diff --git a/beziers.scad b/beziers.scad
index d71d37b..4cc5ee8 100644
--- a/beziers.scad
+++ b/beziers.scad
@@ -1043,9 +1043,9 @@ function bezier_patch_points(patch, u, v) =
     assert(is_num(u) || !is_undef(u[0]))
     assert(is_num(v) || !is_undef(v[0]))
     let(
-        vbezes = [for (i = idx(patch[0])) bezier_points(subindex(patch,i), is_num(u)? [u] : u)]
+        vbezes = [for (i = idx(patch[0])) bezier_points(columns(patch,i), is_num(u)? [u] : u)]
     )
-    [for (i = idx(vbezes[0])) bezier_points(subindex(vbezes,i), is_num(v)? [v] : v)];
+    [for (i = idx(vbezes[0])) bezier_points(columns(vbezes,i), is_num(v)? [v] : v)];
 
 
 // Function: bezier_triangle_point()
@@ -1357,7 +1357,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
     all(row_degen) && all(col_degen) ?  // fully degenerate case
         [EMPTY_VNF, repeat([patch[0][0]],4)] :
     all(row_degen) ?                         // degenerate to a line (top to bottom)
-        let(pts = bezier_points(subindex(patch,0), samplepts))
+        let(pts = bezier_points(columns(patch,0), samplepts))
         [EMPTY_VNF, [pts,pts,[pts[0]],[last(pts)]]] :
     all(col_degen) ?                         // degenerate to a line (left to right)
         let(pts = bezier_points(patch[0], samplepts))
@@ -1366,7 +1366,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
        let(pts = bezier_patch_points(patch, samplepts, samplepts))
        [
         vnf_vertex_array(pts, reverse=!reverse),
-        [subindex(pts,0), subindex(pts,len(pts)-1), pts[0], last(pts)]
+        [columns(pts,0), columns(pts,len(pts)-1), pts[0], last(pts)]
        ] :
     top_degen && bot_degen ?
        let(
@@ -1375,17 +1375,17 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
                         if (splinesteps%2==0) splinesteps+1,
                         each reverse(list([3:2:splinesteps]))
                        ],
-            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)],
+            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
             pts = [
                   [bpatch[0][0]],
-                  for(j=[0:splinesteps-2]) bezier_points(subindex(bpatch,j+1), lerpn(0,1,rowcount[j])),
+                  for(j=[0:splinesteps-2]) bezier_points(columns(bpatch,j+1), lerpn(0,1,rowcount[j])),
                   [last(bpatch[0])]
                   ],
             vnf = vnf_tri_array(pts, reverse=!reverse)
          ) [
             vnf,
             [
-             subindex(pts,0),
+             columns(pts,0),
              [for(row=pts) last(row)],
              pts[0],
              last(pts),
@@ -1404,16 +1404,16 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
            full_degen = len(patch)>=4 && all(select(row_degen,1,ceil(len(patch)/2-1))),
            rowmax = full_degen ? count(splinesteps+1) :
                                  [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps],
-           bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)],
+           bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
            pts = [
                   [bpatch[0][0]],
-                  for(j=[1:splinesteps]) bezier_points(subindex(bpatch,j), lerpn(0,1,rowmax[j]+1))
+                  for(j=[1:splinesteps]) bezier_points(columns(bpatch,j), lerpn(0,1,rowmax[j]+1))
                  ],
            vnf = vnf_tri_array(pts, reverse=!reverse)
         ) [
             vnf,
             [
-             subindex(pts,0),
+             columns(pts,0),
              [for(row=pts) last(row)],
              pts[0],
              last(pts),
diff --git a/bottlecaps.scad b/bottlecaps.scad
index 7d88ac0..2451b7e 100644
--- a/bottlecaps.scad
+++ b/bottlecaps.scad
@@ -1121,7 +1121,7 @@ module sp_neck(diam,type,wall,id,style="L",bead=false, anchor, spin, orient)
 
     isect400 = [for(seg=pair(beadpts)) let(segisect = line_intersection([[T/2,0],[T/2,1]] , seg, LINE, SEGMENT)) if (is_def(segisect)) segisect.y];
 
-    extra_bot = type==400 && bead ? -min(subindex(beadpts,1))+max(isect400) : 0;
+    extra_bot = type==400 && bead ? -min(columns(beadpts,1))+max(isect400) : 0;
     bead_shift = type==400 ? H+max(isect400) : entry[5]+W/2;  // entry[5] is L
     
     attachable(anchor,spin,orient,r=bead ? beadmax : T/2, l=H+extra_bot){
diff --git a/geometry.scad b/geometry.scad
index c328c87..5a4a76c 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -1878,26 +1878,7 @@ function polygon_shift(poly, i) =
 //   move_copies(concat(circ,pent)) circle(r=.1,$fn=32);
 //   color("red") move_copies([pent[0],circ[0]]) circle(r=.1,$fn=32);
 //   color("blue") translate(reindexed[0])circle(r=.1,$fn=32);
-function old_reindex_polygon(reference, poly, return_error=false) =
-    assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
-           "Invalid polygon(s) or incompatible dimensions. " )
-    assert(len(reference)==len(poly), "The polygons must have the same length.")
-    let(
-        dim = len(reference[0]),
-        N = len(reference),
-        fixpoly = dim != 2? poly :
-                  is_polygon_clockwise(reference)
-                  ? clockwise_polygon(poly)
-                  : ccw_polygon(poly),
-        I   = [for(i=reference) 1],
-        val = [ for(k=[0:N-1])
-                    [for(i=[0:N-1])
-                      (reference[i]*poly[(i+k)%N]) ] ]*I,
-        optimal_poly = polygon_shift(fixpoly, max_index(val))
-    )
-    return_error? [optimal_poly, min(poly*(I*poly)-2*val)] :
-    optimal_poly;
-function reindex_polygon(reference, poly, return_error=false) =
+    function reindex_polygon(reference, poly, return_error=false) =
     assert(is_path(reference) && is_path(poly,dim=len(reference[0])),
            "Invalid polygon(s) or incompatible dimensions. " )
     assert(len(reference)==len(poly), "The polygons must have the same length.")
@@ -1919,68 +1900,70 @@ function reindex_polygon(reference, poly, return_error=false) =
     optimal_poly;
 
 
-
 // Function: align_polygon()
 // Usage:
-//   newpoly = align_polygon(reference, poly, angles, [cp]);
+//   newpoly = align_polygon(reference, poly, [angles], [cp], [tran], [return_ind]);
 // Topics: Geometry, Polygons
 // Description:
-//   Tries the list or range of angles to find a rotation of the specified 2D polygon that best aligns
-//   with the reference 2D polygon.  For each angle, the polygon is reindexed, which is a costly operation
-//   so if run time is a problem, use a smaller sampling of angles.  Returns the rotated and reindexed
-//   polygon.
+//   Find the best alignment of a specified 2D polygon with a reference 2D polygon over a set of
+//   transformations.  You can specify a list or range of angles and a centerpoint or you can
+//   give a list of arbitrary 2d transformation matrices.  For each transformation or angle, the polygon is
+//   reindexed, which is a costly operation so if run time is a problem, use a smaller sampling of angles or
+//   transformations.  By default returns the rotated and reindexed polygon.  You can also request that
+//   the best angle or the index into the transformation list be returned.  When you specify an angle
 // Arguments:
 //   reference = reference polygon
 //   poly = polygon to rotate into alignment with the reference
 //   angles = list or range of angles to test
 //   cp = centerpoint for rotations
-// Example(2D): The original hexagon in yellow is not well aligned with the pentagon.  Turning it so the faces line up gives an optimal alignment, shown in red.
-//   $fn=32;
-//   pentagon = subdivide_path(pentagon(side=2),60);
-//   hexagon = subdivide_path(hexagon(side=2.7),60);
-//   color("red") move_copies(scale(1.4,p=align_polygon(pentagon,hexagon,[0:10:359]))) circle(r=.1);
-//   move_copies(concat(pentagon,hexagon))circle(r=.1);
-function old_align_polygon(reference, poly, angles, cp) =
+//   ---
+//   tran = list of 2D transformation matrices to optimize over
+//   return_ind = if true, return the best angle (if you specified angles) or the index into tran otherwise of best alignment
+// Example(2D): Rotating the poorly aligned light gray triangle by 105 degrees produces the best alignment, shown in blue:
+//   ellipse = yscale(3,circle(r=10, $fn=32));
+//   tri = move([-50/3,-9],
+//              subdivide_path([[0,0], [50,0], [0,27]], 32));
+//   aligned = align_polygon(ellipse,tri, [0:5:180]);
+//   color("white")stroke(tri,width=.5,closed=true);
+//   stroke(ellipse, width=.5, closed=true);
+//   color("blue")stroke(aligned,width=.5,closed=true);
+// Example(2D,NoAxes): Translating a triangle (light gray) to the best alignment (blue)
+//   ellipse = yscale(2,circle(r=10, $fn=32));
+//   tri = subdivide_path([[0,0], [27,0], [-7,50]], 32);
+//   T = [for(x=[-10:0], y=[-30:-15]) move([x,y])];
+//   aligned = align_polygon(ellipse,tri, trans=T);
+//   color("white")stroke(tri,width=.5,closed=true);
+//   stroke(ellipse, width=.5, closed=true);
+//   color("blue")stroke(aligned,width=.5,closed=true);
+function align_polygon(reference, poly, angles, cp, trans, return_ind=false) =
+    assert(is_undef(trans) || (is_undef(angles) && is_undef(cp)), "Cannot give both angles/cp and trans as input")
+    let(
+        trans = is_def(trans) ? trans :
+            assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
+                "The `angle` parameter must be a range or a non void list of numbers.")
+            [for(angle=angles) zrot(angle,cp=cp)]
+    )
     assert(is_path(reference,dim=2) && is_path(poly,dim=2),
            "Invalid polygon(s). " )
     assert(len(reference)==len(poly), "The polygons must have the same length.")
-    assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
-            "The `angle` parameter must be a range or a non void list of numbers.")
     let(     // alignments is a vector of entries of the form: [polygon, error]
         alignments = [
-            for(angle=angles)
-            reindex_polygon(
-                reference,
-                zrot(angle,p=poly,cp=cp),
-                return_error=true
-            )
+            for(T=trans)
+              reindex_polygon(
+                  reference,
+                  apply(T,poly),
+                  return_error=true
+              )
         ],
-        best = min_index(subindex(alignments,1))
-    ) alignments[best][0];
-    
-    
-function align_polygon(reference, poly, angles, cp) =
-    assert(is_path(reference,dim=2) && is_path(poly,dim=2),
-           "Invalid polygon(s). " )
-    assert(len(reference)==len(poly), "The polygons must have the same length.")
-    assert( (is_vector(angles) && len(angles)>0) || valid_range(angles),
-            "The `angle` parameter must be a range or a non void list of numbers.")
-    let(     // alignments is a vector of entries of the form: [polygon, error]
-        alignments = [
-            for(angle=angles)
-            reindex_polygon(
-                reference,
-                zrot(angle,p=poly,cp=cp),
-                return_error=true
-            )
-        ],
-        scores = subindex(alignments,1),
+        scores = columns(alignments,1),
         minscore = min(scores),
         minind = [for(i=idx(scores)) if (scores[i]<minscore+EPSILON) i],
-        f=echo(best_angles = select(list(angles), minind)),
+        dummy = is_def(angles) ? echo(best_angles = select(list(angles), minind)):0,
         best = minind[0]
-    ) alignments[best][0];
-
+    )
+    return_ind ? (is_def(angles) ? list(angles)[best] : best)
+    : alignments[best][0];
+    
 
 // Function: are_polygons_equal()
 // Usage:
diff --git a/math.scad b/math.scad
index 360989d..6ac428d 100644
--- a/math.scad
+++ b/math.scad
@@ -994,7 +994,7 @@ function null_space(A,eps=1e-12) =
         zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
     )
     len(zrow)==0 ? [] :
-    transpose(subindex(Q_R[0],zrow));
+    transpose(columns(Q_R[0],zrow));
 
 
 // Function: qr_factor()
diff --git a/paths.scad b/paths.scad
index 74451c8..94ec945 100644
--- a/paths.scad
+++ b/paths.scad
@@ -427,7 +427,7 @@ function resample_path(path, N, spacing, closed=false) =
        distlist = lerpn(0,length,N,false), 
        cuts = _path_cut_points(path, distlist, closed=closed)
    )
-   [ each subindex(cuts,0),
+   [ each columns(cuts,0),
      if (!closed) last(path)     // Then add last point here
    ];
 
@@ -1174,7 +1174,7 @@ function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
    len(fragments)==0? _finished :
     let(
         minxidx = min_index([
-            for (frag=fragments) min(subindex(frag,0))
+            for (frag=fragments) min(columns(frag,0))
         ]),
         result_l = _assemble_a_path_from_fragments(
             fragments=fragments,
diff --git a/regions.scad b/regions.scad
index 6777abb..1a4e690 100644
--- a/regions.scad
+++ b/regions.scad
@@ -267,7 +267,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru
          cornerpts = [for(i=[0:1])
                          [for(k=vector_search(points[i],eps,points[i]))
                              each if (len(k)>1) select(ptind[i],k)]],
-         risect = [for(i=[0:1]) concat(subindex(intersections,i), cornerpts[i])],
+         risect = [for(i=[0:1]) concat(columns(intersections,i), cornerpts[i])],
          counts = [count(len(region1)), count(len(region2))],
          pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
        )
@@ -381,82 +381,6 @@ function region_parts(region) =
 
 // Section: Region Extrusion and VNFs
 
-function _path_path_closest_vertices(path1,path2) =
-    let(
-        dists = [for (i=idx(path1)) let(j=closest_point(path1[i],path2)) [j,norm(path2[j]-path1[i])]],
-        i1 = min_index(subindex(dists,1)),
-        i2 = dists[i1][0]
-    ) [dists[i1][1], i1, i2];
-
-
-function _join_paths_at_vertices(path1,path2,v1,v2) =
-    let(
-        repeat_start = !approx(path1[v1],path2[v2]),
-        path1 = clockwise_polygon(polygon_shift(path1,v1)),
-        path2 = ccw_polygon(polygon_shift(path2,v2))
-    )
-    [
-        each path1,
-        if (repeat_start) path1[0],
-        each path2,
-        if (repeat_start) path2[0],
-    ];                      
-
-        
-// Given a region that is connected and has its outer border in region[0],
-// produces a polygon with the same points that has overlapping connected paths
-// to join internal holes to the outer border.  Output is a single path.  
-function _cleave_connected_region(region) =
-    len(region)==0? [] :
-    len(region)<=1? clockwise_polygon(region[0]) :
-    let(
-        dists = [
-            for (i=[1:1:len(region)-1])
-            _path_path_closest_vertices(region[0],region[i])
-        ],
-        idxi = min_index(subindex(dists,0)),
-        newoline = _join_paths_at_vertices(
-            region[0], region[idxi+1],
-            dists[idxi][1], dists[idxi][2]
-        )
-    ) len(region)==2? clockwise_polygon(newoline) :
-    let(
-        orgn = [
-            newoline,
-            for (i=idx(region))
-                if (i>0 && i!=idxi+1)
-                    region[i]
-        ]
-    )
-    assert(len(orgn)<len(region))
-    _cleave_connected_region(orgn);
-
-
-
-// Function: region_faces()
-// Usage:
-//   vnf = region_faces(region, [transform], [reverse], [vnf]);
-// Description:
-//   Given a region, applies the given transformation matrix to it and makes a VNF of
-//   faces for that region, reversed if necessary.
-// Arguments:
-//   region = The region to make faces for.
-//   transform = If given, a transformation matrix to apply to the faces generated from the region.  Default: No transformation applied.
-//   reverse = If true, reverse the normals of the faces generated from the region.  An untransformed region will have face normals pointing `UP`.  Default: false
-//   vnf = If given, the faces are added to this VNF.  Default: `EMPTY_VNF`
-function region_faces(region, transform, reverse=false, vnf=EMPTY_VNF) =
-    let (
-        regions = region_parts(force_region(region)),
-        vnfs = [
-            if (vnf != EMPTY_VNF) vnf,
-            for (rgn = regions) let(
-                cleaved = path3d(_cleave_connected_region(rgn)),
-                face = is_undef(transform)? cleaved : apply(transform,cleaved),
-                faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i]
-            ) [face, [faceidxs]]
-        ],
-        outvnf = vnf_merge(vnfs)
-    ) outvnf;
 
 
 // Function&Module: linear_sweep()
@@ -564,8 +488,8 @@ function linear_sweep(region, height=1, center, twist=0, scale=1, slices,
                     ) scale([sc,sc,1], p=rot(ang, p=path3d(path,h)))
                 ]
             ) vnf_vertex_array(verts, caps=false, col_wrap=true, style=style),
-            for (rgn = regions) region_faces(rgn, move([0,0,-height/2]), reverse=true),
-            for (rgn = trgns) region_faces(rgn, move([0,0, height/2]), reverse=false)
+            for (rgn = regions) vnf_from_region(rgn, down(height/2), reverse=true),
+            for (rgn = trgns) vnf_from_region(rgn, up(height/2), reverse=false)
         ])
     ) reorient(anchor,spin,orient, cp=cp, vnf=vnf, extent=!anchor_isect, p=vnf);
 
diff --git a/rounding.scad b/rounding.scad
index 3fca7d5..19a6c28 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -1270,7 +1270,7 @@ module convex_offset_extrude(
         // The entry r[i] is [radius,z] for a given layer
         r = move([0,bottom_height],p=concat(
                           reverse(offsets_bot), [[0,0], [0,middle]], move([0,middle], p=offsets_top)));
-        delta = [for(val=deltas(subindex(r,0))) sign(val)];
+        delta = [for(val=deltas(columns(r,0))) sign(val)];
         below=[-thickness,0];
         above=[0,thickness];
            // layers is a list of pairs of the relative positions for each layer, e.g. [0,thickness]
@@ -1937,8 +1937,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
      verify_vert =
        [for(i=[0:N-1],j=[0:4])
          let(
-               vline = concat(select(subindex(top_patch[i],j),2,4),
-                              select(subindex(bot_patch[i],j),2,4))
+               vline = concat(select(columns(top_patch[i],j),2,4),
+                              select(columns(bot_patch[i],j),2,4))
              )
          if (!is_collinear(vline)) [i,j]],
      //verify horiz edges
@@ -1955,8 +1955,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
           "Roundovers interfere with each other on bottom face: either input is self intersecting or top joint length is too large")
     assert(debug || (verify_vert==[] && verify_horiz==[]), "Curvature continuity failed")
     let( 
-        vnf = vnf_merge([ each subindex(top_samples,0),
-                          each subindex(bot_samples,0),
+        vnf = vnf_merge([ each columns(top_samples,0),
+                          each columns(bot_samples,0),
                           for(pts=edge_points) vnf_vertex_array(pts),
                           debug ? vnf_from_polygons(faces) 
                                 : vnf_triangulate(vnf_from_polygons(faces))
@@ -2114,7 +2114,7 @@ function _circle_mask(r) =
 //           ]),
 //           radius = [0,0,each repeat(slotradius,4),0,0], closed=false
 //       )
-//   ) apply(left(max(subindex(slot,0))/2)*fwd(min(subindex(slot,1))), slot);
+//   ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
 //   stroke(slot(15,29,7));
 // Example: A cylindrical container with rounded edges and a rounded finger slot.
 //   function slot(slotwidth, slotheight, slotradius) = let(
@@ -2138,7 +2138,7 @@ function _circle_mask(r) =
 //           ]),
 //           radius = [0,0,each repeat(slotradius,4),0,0], closed=false
 //       )
-//   ) apply(left(max(subindex(slot,0))/2)*fwd(min(subindex(slot,1))), slot);
+//   ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
 //   diam = 80;
 //   wall = 4;
 //   height = 40;
@@ -2162,12 +2162,12 @@ module bent_cutout_mask(r, thickness, path, radius, convexity=10)
   path = clockwise_polygon(path);
   curvepoints = arc(d=thickness, angle = [-180,0]);
   profiles = [for(pt=curvepoints) _cyl_hole(r+pt.x,apply(xscale((r+pt.x)/r), offset(path,delta=thickness/2+pt.y,check_valid=false,closed=true)))];
-  pathx = subindex(path,0);
+  pathx = columns(path,0);
   minangle = (min(pathx)-thickness/2)*360/(2*PI*r);
   maxangle = (max(pathx)+thickness/2)*360/(2*PI*r);
   mindist = (r+thickness/2)/cos((maxangle-minangle)/2);
   assert(maxangle-minangle<180,"Cutout angle span is too large.  Must be smaller than 180.");
-  zmean = mean(subindex(path,1));
+  zmean = mean(columns(path,1));
   innerzero = repeat([0,0,zmean], len(path));
   outerpt = repeat( [1.5*mindist*cos((maxangle+minangle)/2),1.5*mindist*sin((maxangle+minangle)/2),zmean], len(path));
   vnf_polyhedron(vnf_vertex_array([innerzero, each profiles, outerpt],col_wrap=true),convexity=convexity);
diff --git a/shapes2d.scad b/shapes2d.scad
index bfeed33..b48d117 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -363,7 +363,7 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false
                     )
                     each arc(N=steps, cp=p, r=rounding, start=a+180/n, angle=-360/n)
                 ],
-                maxx_idx = max_index(subindex(path2,0)),
+                maxx_idx = max_index(columns(path2,0)),
                 path3 = polygon_shift(path2,maxx_idx)
             ) path3
         ),
@@ -961,7 +961,7 @@ function teardrop2d(r, ang=45, cap_h, d, anchor=CENTER, spin=0) =
                 [-cap_w,cap_h]
             ], closed=true
         ),
-        maxx_idx = max_index(subindex(path,0)),
+        maxx_idx = max_index(columns(path,0)),
         path2 = polygon_shift(path,maxx_idx)
     ) reorient(anchor,spin, two_d=true, path=path2, p=path2);
 
@@ -1016,7 +1016,7 @@ function glued_circles(r, spread=10, tangent=30, d, anchor=CENTER, spin=0) =
             [for (i=[0:1:lobesegs]) let(a=sa1+i*lobestep+180) r  * [cos(a),sin(a)] + cp1],
             tangent==0? [] : [for (i=[0:1:arcsegs])  let(a=ea2-i*arcstep)      r2 * [cos(a),sin(a)] + cp2]
         ),
-        maxx_idx = max_index(subindex(path,0)),
+        maxx_idx = max_index(columns(path,0)),
         path2 = reverse_polygon(polygon_shift(path,maxx_idx))
     ) reorient(anchor,spin, two_d=true, path=path2, extent=true, p=path2);
 
diff --git a/shapes3d.scad b/shapes3d.scad
index 6bc29f1..8d63d63 100644
--- a/shapes3d.scad
+++ b/shapes3d.scad
@@ -2211,7 +2211,7 @@ module path_text(path, text, font, size, thickness, lettersize, offset=0, revers
   usernorm = is_def(normal);
   usetop = is_def(top);
   
-  normpts = is_undef(normal) ? (reverse?1:-1)*subindex(pts,3) : _cut_interp(pts,path, normal);
+  normpts = is_undef(normal) ? (reverse?1:-1)*columns(pts,3) : _cut_interp(pts,path, normal);
   toppts = is_undef(top) ? undef : _cut_interp(pts,path,top);
   for(i=idx(text))
       let( tangent = pts[i][2] )
diff --git a/skin.scad b/skin.scad
index e2ad45a..fcc3f55 100644
--- a/skin.scad
+++ b/skin.scad
@@ -255,14 +255,16 @@
 //   skin([regular_ngon(n=4, r=4), regular_ngon(n=5,r=5)], z=[0,4], refine=10, slices=10);
 // Example(FlatSpin,VPD=35): Connecting square to shifted pentagon using "direct" method.
 //   skin([regular_ngon(n=4, r=4), right(4,p=regular_ngon(n=5,r=5))], z=[0,4], refine=10, slices=10);
-// Example(FlatSpin,VPD=35): To improve the look, you can actually rotate the polygons for a more symmetric pattern of lines.   You have to resample yourself before calling `align_polygon` and you should choose a length that is a multiple of both polygon lengths.
-//   sq = subdivide_path(regular_ngon(n=4, r=4),40);
-//   pent = subdivide_path(regular_ngon(n=5,r=5),40);
-//   skin([sq, align_polygon(sq,pent,[0:1:360/5])], z=[0,4], slices=10);
-// Example(FlatSpin,VPD=35): For the shifted pentagon we can also align, making sure to pass an appropriate centerpoint to `align_polygon`.
-//   sq = subdivide_path(regular_ngon(n=4, r=4),40);
-//   pent = right(4,p=subdivide_path(regular_ngon(n=5,r=5),40));
-//   skin([sq, align_polygon(sq,pent,[0:1:360/5],cp=[4,0])], z=[0,4], refine=10, slices=10);
+// Example(FlatSpin,VPD=185): In this example reindexing does not fix the orientation of the triangle because it happens in 3d within skin(), so we have to reverse the triangle manually
+//   ellipse = yscale(3,circle(r=10, $fn=32));
+//   tri = move([-50/3,-9],[[0,0], [50,0], [0,27]]);
+//   skin([ellipse, reverse(tri)], z=[0,20], slices=20, method="reindex");
+// Example(FlatSpin,VPD=185): You can get a nicer transition by rotating the polygons for better alignment.  You have to resample yourself before calling `align_polygon`. The orientation is fixed so we do not need to reverse.  
+//   ellipse = yscale(3,circle(r=10, $fn=32));
+//   tri = move([-50/3,-9],
+//              subdivide_path([[0,0], [50,0], [0,27]], 32));
+//   aligned = align_polygon(ellipse,tri, [0:5:180]);
+//   skin([ellipse, aligned], z=[0,20], slices=20);
 // Example(FlatSpin,VPD=35): The "distance" method is a completely different approach.
 //   skin([regular_ngon(n=4, r=4), regular_ngon(n=5,r=5)], z=[0,4], refine=10, slices=10, method="distance");
 // Example(FlatSpin,VPD=35,VPT=[0,0,4]): Connecting pentagon to heptagon inserts two triangular faces on each side
@@ -983,7 +985,7 @@ function path_sweep2d(shape, path, closed=false, caps, quality=1, style="min_edg
                      [for(pt = profile)
                         let( 
                             ofs = offset(path, delta=-flip*pt.x, return_faces=true,closed=closed, quality=quality),
-                            map = subindex(_ofs_vmap(ofs,closed=closed),1)
+                            map = columns(_ofs_vmap(ofs,closed=closed),1)
                         ) 
                         select(path3d(ofs[0],pt.y),map)
                       ]
@@ -1111,8 +1113,8 @@ function sweep(shape, transforms, closed=false, caps, style="min_edge") =
             for (rgn=regions) each [
                 for (path=rgn)
                     sweep(path, transforms, closed=closed, caps=false),
-                if (fullcaps[0]) region_faces(rgn, transform=transforms[0], reverse=true),
-                if (fullcaps[1]) region_faces(rgn, transform=last(transforms)),
+                if (fullcaps[0]) vnf_from_region(rgn, transform=transforms[0], reverse=true),
+                if (fullcaps[1]) vnf_from_region(rgn, transform=last(transforms)),
             ],
         ],
         vnf = vnf_merge(vnfs)
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 3aec4b7..4a62f0f 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -433,14 +433,14 @@ module test_add_scalar() {
 test_add_scalar();
 
 
-module test_subindex() {
+module test_columns() {
     v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-    assert(subindex(v,2) == [3, 7, 11, 15]);
-    assert(subindex(v,[2]) == [[3], [7], [11], [15]]);
-    assert(subindex(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
-    assert(subindex(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
+    assert(columns(v,2) == [3, 7, 11, 15]);
+    assert(columns(v,[2]) == [[3], [7], [11], [15]]);
+    assert(columns(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
+    assert(columns(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
 }
-test_subindex();
+test_columns();
 
 
 // Need decision about behavior for out of bounds ranges, empty ranges
@@ -532,7 +532,7 @@ module test_hstack() {
     a = hstack(v1,v2);   
     b = hstack(v1,v2,v3);
     c = hstack([M,v1,M]);
-    d = hstack(subindex(M,0), subindex(M,[1, 2]));
+    d = hstack(columns(M,0), columns(M,[1, 2]));
     assert_equal(a,[[2, 5], [3, 6], [4, 7]]);
     assert_equal(b,[[2, 5, 8], [3, 6, 9], [4, 7, 10]]);
     assert_equal(c,[[1, 0, 0, 2, 1, 0, 0], [0, 1, 0, 3, 0, 1, 0], [0, 0, 1, 4, 0, 0, 1]]);
diff --git a/tests/test_math.scad b/tests/test_math.scad
index a5a5843..0411b07 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -994,7 +994,7 @@ module test_linear_solve(){
                                                          -8.378819388897780e-01,
                                                           2.330507118860985e-01,
                                                           8.511278195488737e-01]);
-  assert_approx(linear_solve(subindex(M,[0:2]), [2,4,4,4]),
+  assert_approx(linear_solve(columns(M,[0:2]), [2,4,4,4]),
                  [-2.457142857142859e-01,
                    5.200000000000000e-01,
                    7.428571428571396e-02]);
diff --git a/threading.scad b/threading.scad
index 40fcc93..e9f210c 100644
--- a/threading.scad
+++ b/threading.scad
@@ -926,8 +926,8 @@ module generic_threaded_rod(
     assert(higang1 < twist/2);
     assert(higang2 < twist/2);
     prof3d = path3d(profile);
-    pdepth = -min(subindex(profile,1));
-    pmax = pitch * max(subindex(profile,1));
+    pdepth = -min(columns(profile,1));
+    pmax = pitch * max(columns(profile,1));
     rmax = max(_r1,_r2)+pmax;
     depth = pdepth * pitch;
     dummy1 = assert(_r1>depth && _r2>depth, "Screw profile deeper than rod radius");
@@ -1087,7 +1087,7 @@ module generic_threaded_nut(
     bevel1 = first_defined([bevel1,bevel,false]);
     bevel2 = first_defined([bevel2,bevel,false]);
     dummy1 = assert(is_num(pitch) && pitch>0);
-    depth = -pitch*min(subindex(profile,1));
+    depth = -pitch*min(columns(profile,1));
     attachable(anchor,spin,orient, size=[od/cos(30),od,h]) {
         difference() {
             cyl(d=od/cos(30), h=h, center=true, $fn=6,chamfer1=bevel1?depth:undef,chamfer2=bevel2?depth:undef);
diff --git a/vectors.scad b/vectors.scad
index e75216a..ac9c382 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -263,11 +263,50 @@ function vector_axis(v1,v2=undef,v3=undef) =
 
 
 
+// Function: min_index()
+// Usage:
+//   idx = min_index(vals);
+//   idxlist = min_index(vals, all=true);
+// Topics: List Handling
+// See Also: max_index(), list_increasing(), list_decreasing()
+// Description:
+//   Returns the index of the first occurrence of the minimum value in the given list. 
+//   If `all` is true then returns a list of all indices where the minimum value occurs.
+// Arguments:
+//   vals = vector of values
+//   all = set to true to return indices of all occurences of the minimum.  Default: false
+// Example:
+//   a = min_index([5,3,9,6,2,7,8,2,1]); // Returns: 8
+//   b = min_index([5,3,9,6,2,7,8,2,7],all=true); // Returns: [4,7]
+function min_index(vals, all=false) =
+    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
+    all ? search(min(vals),vals,0) : search(min(vals), vals)[0];
+
+
+// Function: max_index()
+// Usage:
+//   idx = max_index(vals);
+//   idxlist = max_index(vals, all=true);
+// Topics: List Handling
+// See Also: min_index(), list_increasing(), list_decreasing()
+// Description:
+//   Returns the index of the first occurrence of the maximum value in the given list. 
+//   If `all` is true then returns a list of all indices where the maximum value occurs.
+// Arguments:
+//   vals = vector of values
+//   all = set to true to return indices of all occurences of the maximum.  Default: false
+// Example:
+//   max_index([5,3,9,6,2,7,8,9,1]); // Returns: 2
+//   max_index([5,3,9,6,2,7,8,9,1],all=true); // Returns: [2,7]
+function max_index(vals, all=false) =
+    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
+    all ? search(max(vals),vals,0) : search(max(vals), vals)[0];
+
+
 
 // Section: Vector Searching
 
 
-
 // Function: closest_point()
 // Usage:
 //   index = closest_point(pt, points);
@@ -504,8 +543,8 @@ function vector_nearest(query, k, target) =
             "More results are requested than the number of points.")
     tgpts
     ?   let( tree = _bt_tree(target, count(len(target))) )
-        subindex(_bt_nearest( query, k, target,  tree),0)
-    :   subindex(_bt_nearest( query, k, target[0], target[1]),0);
+        columns(_bt_nearest( query, k, target,  tree),0)
+    :   columns(_bt_nearest( query, k, target[0], target[1]),0);
 
 
 //Ball tree nearest
diff --git a/vnf.scad b/vnf.scad
index 24fc232..9ba0e0d 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -287,9 +287,9 @@ function vnf_tri_array(points, row_wrap=false, reverse=false, vnf=EMPTY_VNF) =
 //   and eliminates any faces with fewer than 3 vertices.  
 //   (Unreferenced vertices of the input VNFs are not dropped.)
 // Arguments:
-//   vnfs - a list of the VNFs to merge in one VNF.
-//   cleanup - when true, consolidates the duplicate vertices of the merge. Default: false
-//   eps - the tolerance in finding duplicates when cleanup=true. Default: EPSILON
+//   vnfs = a list of the VNFs to merge in one VNF.
+//   cleanup = when true, consolidates the duplicate vertices of the merge. Default: false
+//   eps = the tolerance in finding duplicates when cleanup=true. Default: EPSILON
 function vnf_merge(vnfs, cleanup=false, eps=EPSILON) =
     is_vnf(vnfs) ? vnf_merge([vnfs], cleanup, eps) :
     assert( is_vnf_list(vnfs) , "Improper vnf or vnf list")  
@@ -350,6 +350,95 @@ function vnf_from_polygons(polygons) =
 
 
 
+
+function _path_path_closest_vertices(path1,path2) =
+    let(
+        dists = [for (i=idx(path1)) let(j=closest_point(path1[i],path2)) [j,norm(path2[j]-path1[i])]],
+        i1 = min_index(columns(dists,1)),
+        i2 = dists[i1][0]
+    ) [dists[i1][1], i1, i2];
+
+
+function _join_paths_at_vertices(path1,path2,v1,v2) =
+    let(
+        repeat_start = !approx(path1[v1],path2[v2]),
+        path1 = clockwise_polygon(polygon_shift(path1,v1)),
+        path2 = ccw_polygon(polygon_shift(path2,v2))
+    )
+    [
+        each path1,
+        if (repeat_start) path1[0],
+        each path2,
+        if (repeat_start) path2[0],
+    ];                      
+
+        
+// Given a region that is connected and has its outer border in region[0],
+// produces a polygon with the same points that has overlapping connected paths
+// to join internal holes to the outer border.  Output is a single path.  
+function _cleave_connected_region(region) =
+    len(region)==0? [] :
+    len(region)<=1? clockwise_polygon(region[0]) :
+    let(
+        dists = [
+            for (i=[1:1:len(region)-1])
+            _path_path_closest_vertices(region[0],region[i])
+        ],
+        idxi = min_index(columns(dists,0)),
+        newoline = _join_paths_at_vertices(
+            region[0], region[idxi+1],
+            dists[idxi][1], dists[idxi][2]
+        )
+    ) len(region)==2? clockwise_polygon(newoline) :
+    let(
+        orgn = [
+            newoline,
+            for (i=idx(region))
+                if (i>0 && i!=idxi+1)
+                    region[i]
+        ]
+    )
+    assert(len(orgn)<len(region))
+    _cleave_connected_region(orgn);
+
+
+
+// Function: vnf_from_region()
+// Usage:
+//   vnf = vnf_from_region(region, [transform], [reverse], [vnf]);
+// Description:
+//   Given a (two-dimensional) region, applies the given transformation matrix to it and makes a triangulated VNF of
+//   faces for that region, reversed if desired. 
+// Arguments:
+//   region = The region to conver to a vnf.
+//   transform = If given, a transformation matrix to apply to the faces generated from the region.  Default: No transformation applied.
+//   reverse = If true, reverse the normals of the faces generated from the region.  An untransformed region will have face normals pointing `UP`.  Default: false
+//   vnf = If given, the faces are added to this VNF.  Default: `EMPTY_VNF`
+// Example(3D):
+//   region = [square([20,10],center=true),
+//             right(5,square(4,center=true)),
+//             left(5,square(6,center=true))];
+//   vnf = vnf_from_region(region);
+//   color("gray")down(.125)
+//        linear_extrude(height=.125)region(region);
+//   vnf_wireframe(vnf,width=.25);
+function vnf_from_region(region, transform, reverse=false, vnf=EMPTY_VNF) =
+    let (
+        regions = region_parts(force_region(region)),
+        vnfs = [
+            if (vnf != EMPTY_VNF) vnf,
+            for (rgn = regions) let(
+                cleaved = path3d(_cleave_connected_region(rgn)),
+                face = is_undef(transform)? cleaved : apply(transform,cleaved),
+                faceidxs = reverse? [for (i=[len(face)-1:-1:0]) i] : [for (i=[0:1:len(face)-1]) i]
+            ) [face, [faceidxs]]
+        ],
+        outvnf = vnf_merge(vnfs)
+    )
+    vnf_triangulate(outvnf);
+
+
+
 // Section: VNF Testing and Access
 
 
@@ -390,7 +479,7 @@ function vnf_faces(vnf) = vnf[1];
 // Usage:
 //   rvnf = vnf_reverse_faces(vnf);
 // Description:
-//   Reverses the facing of all the faces in the given VNF.
+//   Reverses the orientation of all the faces in the given VNF.
 function vnf_reverse_faces(vnf) =
     [vnf[0], [for (face=vnf[1]) reverse(face)]];
 
@@ -407,14 +496,47 @@ function vnf_quantize(vnf,q=pow(2,-12)) =
     [[for (pt = vnf[0]) quant(pt,q)], vnf[1]];
 
 
+// Function: vnf_drop_unused_points()
+// Usage:
+//   clean_vnf=vnf_drop_unused_points(vnf);
+// Description:
+//   Remove all unreferenced vertices from a VNF.  Note that in most
+//   cases unreferenced vertices cause no harm, and this function may
+//   be slow on large VNFs.
+function vnf_drop_unused_points(vnf) =
+    let(
+        flat = flatten(vnf[1]),
+        ind  = _link_indicator(flat,0,len(vnf[0])-1),
+        verts = [for(i=idx(vnf[0])) if(ind[i]==1) vnf[0][i] ], 
+        map   = cumsum(ind) 
+    )
+    [ verts, [for(face=vnf[1]) [for(v=face) map[v]-1 ] ] ];
+
+function _link_indicator(l,imin,imax) =
+    len(l) == 0  ? repeat(imax-imin+1,0) :
+    imax-imin<100 || len(l)<400 ? [for(si=search(list([imin:1:imax]),l,1)) si!=[] ? 1: 0 ] : 
+    let( 
+        pivot   = floor((imax+imin)/2),
+        lesser  = [ for(li=l) if( li< pivot) li ],
+        greater = [ for(li=l) if( li> pivot) li ] 
+    )
+    concat( _link_indicator(lesser ,imin,pivot-1), 
+            search(pivot,l,1) ? 1 : 0 ,
+            _link_indicator(greater,pivot+1,imax) ) ;
+
 // Function: vnf_triangulate()
 // Usage:
 //   vnf2 = vnf_triangulate(vnf);
 // Description:
-//   Triangulates faces in the VNF that have more than 3 vertices.  
+//   Triangulates faces in the VNF that have more than 3 vertices.
+// Example:
+//   include <BOSL2/polyhedra.scad>
+//   vnf = zrot(33,regular_polyhedron_info("vnf", "dodecahedron", side=12));
+//   vnf_polyhedron(vnf);
+//   triangulated = vnf_triangulate(vnf);
+//   color("red")vnf_wireframe(triangulated,width=.3);
 function vnf_triangulate(vnf) =
     let(
-        vnf = is_vnf_list(vnf)? vnf_merge(vnf) : vnf,
         verts = vnf[0],
         faces = [for (face=vnf[1]) each len(face)==3 ? [face] : 
                                          polygon_triangulate(verts, face)]
@@ -446,7 +568,7 @@ function vnf_slice(vnf,dir,cuts) =
 
 function _split_polygon_at_x(poly, x) =
     let(
-        xs = subindex(poly,0)
+        xs = columns(poly,0)
     ) (min(xs) >= x || max(xs) <= x)? [poly] :
     let(
         poly2 = [
@@ -709,7 +831,7 @@ function vnf_halfspace(plane, vnf, closed=true) =
       let(
            M = project_plane(plane),
            faceregion = [for(path=newpaths) path2d(apply(M,select(newvert,path)))],
-           facevnf = region_faces(faceregion,transform=rot_inverse(M),reverse=true)
+           facevnf = vnf_from_region(faceregion,transform=rot_inverse(M),reverse=true)
       )
       vnf_merge([[newvert, faces_edges_vertices[0]], facevnf]);