diff --git a/arrays.scad b/arrays.scad
index 40b232f..0df54af 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -59,14 +59,17 @@ function _same_type(a,b, depth) =
 //   Returns a portion of a list, wrapping around past the beginning, if end<start. 
 //   The first item is index 0. Negative indexes are counted back from the end.
 //   The last item is -1.  If only the `start` index is given, returns just the value
-//   at that position.
+//   at that position when `start` is a number or the selected list of entries when `start` is
+//   a list of indices or a range.
 // Usage:
 //   item = select(list,start);
+//   item = select(list,[s:d:e]);
+//   item = select(list,[i0,i1...,ik]);
 //   list = select(list,start,end);
 // Arguments:
 //   list = The list to get the portion of.
-//   start = The index of the first item.
-//   end = The index of the last item.
+//   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()
 // Example:
 //   l = [3,4,5,6,7,8,9];
@@ -78,7 +81,7 @@ function _same_type(a,b, depth) =
 //   f = select(l, 4);      // Returns 7
 //   g = select(l, -2);     // Returns 8
 //   h = select(l, [1:3]);  // Returns [4,5,6]
-//   i = select(l, [1,3]);  // Returns [4,6]
+//   i = select(l, [3,1]);  // Returns [6,4]
 function select(list, start, end) =
     assert( is_list(list) || is_string(list), "Invalid list.")
     let(l=len(list))
@@ -89,7 +92,7 @@ function select(list, start, end) =
               ? list[ (start%l+l)%l ]
               : assert( is_list(start) || is_range(start), "Invalid start parameter")
                 [for (i=start) list[ (i%l+l)%l ] ]
-          : assert(is_finite(start), "Invalid start parameter.")
+          : assert(is_finite(start), "When `end` is given, `start` parameter should be a number.")
             assert(is_finite(end), "Invalid end parameter.")
             let( s = (start%l+l)%l, e = (end%l+l)%l )
             (s <= e)
@@ -946,6 +949,31 @@ function shuffle(list,seed) =
     ) 
     concat(shuffle(left), shuffle(right));
 
+// idx should be an index of the arrays l[i]
+function _group_sort_by_index(l,idx) =
+	len(l) == 0 ? [] :
+	len(l) == 1 ? [l] : 
+	let( pivot   = l[floor(len(l)/2)][idx],
+       equal   = [ for(li=l) if( li[idx]==pivot) li ] ,
+       lesser  = [ for(li=l) if( li[idx]< pivot) li ] ,
+       greater = [ for(li=l) if( li[idx]> pivot) li ]
+  )
+  concat( _group_sort_by_index(lesser,idx), 
+          [equal], 
+          _group_sort_by_index(greater,idx) ) ;  
+
+function _group_sort(l) =
+	len(l) == 0 ? [] : 
+	len(l) == 1 ? [l] : 
+	let( pivot   = l[floor(len(l)/2)],
+       equal   = [ for(li=l) if( li==pivot) li ] ,
+       lesser  = [ for(li=l) if( li< pivot) li ] ,
+       greater = [ for(li=l) if( li> pivot) li ]
+  )
+  concat( _group_sort(lesser), 
+          [equal], 
+          _group_sort(greater) ) ;    
+
 
 // Sort a vector of scalar values with the native comparison operator
 // all elements should have the same type.
@@ -1007,7 +1035,7 @@ function _sort_general(arr, idx=undef, indexed=false) =
       
 // lexical sort using compare_vals()
 function _lexical_sort(arr) = 
-    arr==[] ? [] : len(arr)==1? arr : 
+    len(arr)<=1? arr : 
     let( pivot = arr[floor(len(arr)/2)] )
     let(
         lesser  = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ],
@@ -1034,7 +1062,7 @@ function _indexed_sort(arrind) =
 // Usage:
 //   slist = sort(list, <idx>);
 // Topics: List Handling
-// See Also: shuffle(), sortidx(), unique(), unique_count()
+// See Also: shuffle(), sortidx(), unique(), unique_count(), group_sort()
 // Description:
 //   Sorts the given list in lexicographic order. If the input is a homogeneous simple list or a homogeneous 
 //   list of vectors (see function is_homogeneous), the sorting method uses the native comparison operator and is faster. 
@@ -1075,7 +1103,7 @@ function sort(list, idx=undef) =
 // Usage:
 //   idxlist = sortidx(list, <idx>);
 // Topics: List Handling
-// See Also: shuffle(), sort(), unique(), unique_count()
+// See Also: shuffle(), sort(), group_sort(), unique(), unique_count()
 // Description:
 //   Given a list, sort it as function `sort()`, and returns
 //   a list of indexes into the original list in that sorted order.
@@ -1123,26 +1151,69 @@ function sortidx(list, idx=undef) =
     : _sort_general(list,idx,indexed=true);
         
 
+// Function: group_sort()
+// Usage:
+//   ulist = group_sort(list);
+// Topics: List Handling
+// See Also: shuffle(), sort(), sortidx(), unique(), unique_count()
+// Description:
+//   Given a list of values, returns the sorted list with all repeated items grouped in a list.
+//   When the list entries are themselves lists, the sorting may be done based on the `idx` entry
+//   of those entries, that should be numbers. 
+//   The result is always a list of lists. 
+// Arguments:
+//   list = The list to sort.
+//   idx = If given, do the comparison based just on the specified index. Default: zero.
+// Example:
+//   sorted = group_sort([5,2,8,3,1,3,8,7,5]);  // Returns: [[1],[2],[3,3],[5,5],[7],[8,8]]
+//   sorted2 = group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0);  // Returns: [[[2,"b"],[2,"e"]], [[5,"a"],[5,"c"]], [[3,"d"]] ]
+function group_sort(list, idx) = 
+    assert(is_list(list), "Input should be a list." )
+    assert(is_undef(idx) || (is_finite(idx) && idx>=0) , "Invalid index." )
+    len(list)<=1 ? [list] :
+		is_vector(list)? _group_sort(list) :
+		let( idx = is_undef(idx) ? 0 : idx )
+		assert( [for(entry=list) if(!is_list(entry) || len(entry)<idx || !is_num(entry[idx]) ) 1]==[],
+            "Some entry of the list is a list shorter than `idx` or the indexed entry of it is not a number."	)
+		_group_sort_by_index(list,idx);
+        
+
 // Function: unique()
 // Usage:
 //   ulist = unique(list);
 // Topics: List Handling
 // See Also: shuffle(), sort(), sortidx(), unique_count()
 // Description:
-//   Returns a sorted list with all repeated items removed.
+//   Given a string or a list returns the sorted string or the sorted list with all repeated items removed.
+//   The sorting order of non homogeneous lists is the function `sort` order.
 // Arguments:
 //   list = The list to uniquify.
 // Example:
 //   sorted = unique([5,2,8,3,1,3,8,7,5]);  // Returns: [1,2,3,5,7,8]
+//   sorted = unique("axdbxxc");  // Returns: "abcdx"
+//   sorted = unique([true,2,"xba",[1,0],true,[0,0],3,"a",[0,0],2]); // Returns: [true,2,3,"a","xba",[0,0],[1,0]]
 function unique(list) =
     assert(is_list(list)||is_string(list), "Invalid input." )
     is_string(list)? str_join(unique([for (x = list) x])) :
     len(list)<=1? list : 
-    let( sorted = sort(list))
-    [   for (i=[0:1:len(sorted)-1])
-            if (i==0 || (sorted[i] != sorted[i-1]))
-                sorted[i]
-    ];
+    is_homogeneous(list,1) && ! is_list(list[0])
+		?   _unique_sort(list)
+    :   let( sorted = sort(list))
+				[   for (i=[0:1:len(sorted)-1])
+								if (i==0 || (sorted[i] != sorted[i-1]))
+										sorted[i]
+				];
+
+function _unique_sort(l) =
+		len(l) <= 1 ? l : 
+		let( pivot   = l[floor(len(l)/2)],
+				 equal   = [ for(li=l) if( li==pivot) li ] ,
+				 lesser  = [ for(li=l) if( li<pivot ) li ] ,
+				 greater = [ for(li=l) if( li>pivot) li ]
+		)
+		concat( _unique_sort(lesser), 
+						equal[0], 
+						_unique_sort(greater) ) ;    
 
 
 // Function: unique_count()
@@ -1160,11 +1231,14 @@ function unique(list) =
 function unique_count(list) =
     assert(is_list(list) || is_string(list), "Invalid input." )
     list == [] ? [[],[]] : 
-    let( list=sort(list) )
-    let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] )
-    [ select(list,ind), deltas( concat(ind,[len(list)]) ) ];
-
+    is_homogeneous(list,1) && ! is_list(list[0])
+    ?		let( sorted = _group_sort(list) )
+				[ [for(s=sorted) s[0] ], [for(s=sorted) len(s) ] ]
+		:		let( list=sort(list) )
+				let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] )
+				[ select(list,ind), deltas( concat(ind,[len(list)]) ) ];
 
+		
 // Section: List Iteration Helpers
 
 // Function: idx()
@@ -1692,7 +1766,7 @@ function submatrix_set(M,A,m=0,n=0) =
 // Arguments:
 //   v = The list of items to group.
 //   cnt = The number of items to put in each grouping.  Default:2
-//   dflt = The default value to fill in with is the list is not a multiple of `cnt` items long.  Default: 0
+//   dflt = The default value to fill in with if the list is not a multiple of `cnt` items long.  Default: 0
 // Example:
 //   v = [1,2,3,4,5,6];
 //   a = array_group(v,2) returns [[1,2], [3,4], [5,6]]
diff --git a/coords.scad b/coords.scad
index b32ad73..f0b4272 100644
--- a/coords.scad
+++ b/coords.scad
@@ -237,7 +237,7 @@ function project_plane(plane,p) =
           let(plane = plane_from_points(plane))
           assert(is_def(plane), "Point list is not coplanar")
           project_plane(plane)
-    : assert(is_def(p), str("Invalid plane specification",plane))
+    : assert(is_def(p), str("Invalid plane specification: ",plane))
       is_vnf(p) ? [project_plane(plane,p[0]), p[1]] 
     : is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ?  // bezier patch or region
            [for(plist=p) project_plane(plane,plist)]
diff --git a/distributors.scad b/distributors.scad
index 6010b12..ee60d6b 100644
--- a/distributors.scad
+++ b/distributors.scad
@@ -34,7 +34,7 @@ module move_copies(a=[[0,0,0]])
     assert(is_list(a));
     for ($idx = idx(a)) {
         $pos = a[$idx];
-        assert(is_vector($pos));
+        assert(is_vector($pos),"move_copies offsets should be a 2d or 3d vector.");
         translate($pos) children();
     }
 }
diff --git a/geometry.scad b/geometry.scad
index 27adf7f..0221a2b 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -19,15 +19,8 @@
 //   edge = Array of two points forming the line segment to test against.
 //   eps = Tolerance in geometric comparisons.  Default: `EPSILON` (1e-9)
 function point_on_segment2d(point, edge, eps=EPSILON) =
-    assert( is_vector(point,2), "Invalid point." )
     assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
-    assert( _valid_line(edge,2,eps=eps), "Invalid segment." )
-    let( dp = point-edge[0],
-         de = edge[1]-edge[0],
-         ne = norm(de) )
-    ( dp*de >= -eps*ne )
-    && ( (dp-de)*de <= eps*ne )                  // point projects on the segment
-    && _dist2line(point-edge[0],unit(de))<eps;   // point is on the line
+    point_segment_distance(point, edge)<eps;   
 
 
 //Internal - distance from point `d` to the line passing through the origin with unit direction n
@@ -44,7 +37,7 @@ function _point_above_below_segment(point, edge) =
 //Internal
 function _valid_line(line,dim,eps=EPSILON) =
     is_matrix(line,2,dim)
-    && ! approx(norm(line[1]-line[0]), 0, eps);
+    && norm(line[1]-line[0])>eps*max(norm(line[1]),norm(line[0]));
 
 //Internal
 function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
@@ -87,7 +80,7 @@ function collinear(a, b, c, eps=EPSILON) =
 
 // Function: point_line_distance()
 // Usage:
-//   point_line_distance(pt, line);
+//   point_line_distance(line, pt);
 // Description:
 //   Finds the perpendicular distance of a point `pt` from the line `line`.
 // Arguments:
@@ -113,6 +106,8 @@ function point_line_distance(pt, line) =
 //   dist = point_segment_distance([3,8], [[-10,0], [10,0]]);  // Returns: 8
 //   dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]);  // Returns: 5
 function point_segment_distance(pt, seg) =
+    assert( is_matrix(concat([pt],seg),3),
+            "Input should be a point and a valid segment with the dimension equal to the point." )
     norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) :
     norm(pt-segment_closest_point(seg,pt));
 
@@ -129,34 +124,9 @@ function point_segment_distance(pt, seg) =
 //   dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]);  // Returns: 5
 //   dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]);  // Returns: 0
 function segment_distance(seg1, seg2) =
-    let(
-        dseg1 = seg1[1]-seg1[0],
-        dseg2 = seg2[1]-seg2[0],
-        A  = [ [dseg1*dseg1, -dseg1*dseg2],
-              [-dseg2*dseg1, dseg2*dseg2] ],
-        b  = -[ dseg1, -dseg2 ]*(seg1[0]-seg2[0]),
-        uv = linear_solve(A,b)
-    )
-    !uv ?
-        norm(dseg1)<EPSILON
-          ? norm(dseg2)<EPSILON
-              ? norm(seg1[0]-seg2[0])
-              : point_segment_distance(seg1[0],seg2)
-          : point_segment_distance(seg2[0],seg1) :
-    uv[0]>=0 && uv[0]<=1 ?
-        let( p1 = seg1[0] + uv[0]*dseg1 )
-        uv[1]>=0 && uv[1]<=1
-          ? norm(p1 - (seg2[0] + uv[1]*dseg2) )
-          : min(norm(p1-seg2[0]), norm(p1-seg2[1])) :
-        uv[1]>=0 && uv[1]<=1
-          ? let( p2 = seg2[0] + uv[1]*dseg2 )
-            min(norm(p2-seg1[0]), norm(p2-seg1[1]))
-          : min(
-                norm(seg1[0]-seg2[0]),
-                norm(seg1[0]-seg2[1]),
-                norm(seg1[1]-seg2[0]),
-                norm(seg1[1]-seg2[1])
-            );
+    assert( is_matrix(concat(seg1,seg2),4),
+            "Inputs should be two valid segments." )
+    convex_distance(seg1,seg2);
 
 
 // Function: line_normal()
@@ -494,17 +464,9 @@ function ray_closest_point(ray,pt) =
 //   color("blue") translate(pt) sphere(r=1,$fn=12);
 //   color("red") translate(p2) sphere(r=1,$fn=12);
 function segment_closest_point(seg,pt) =
-    assert(_valid_line(seg), "Invalid segment." )
-    assert(len(pt)==len(seg[0]), "Incompatible dimensions." )
-    approx(seg[0],seg[1])? seg[0] :
-    let(
-        seglen = norm(seg[1]-seg[0]),
-        segvec = (seg[1]-seg[0])/seglen,
-        projection = (pt-seg[0]) * segvec
-    )
-    projection<=0 ? seg[0] :
-    projection>=seglen ? seg[1] :
-    seg[0] + projection*segvec;
+    assert( is_matrix(concat([pt],seg),3) ,
+            "Invalid point or segment or incompatible dimensions." )
+    pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0];
 
 
 // Function: line_from_points()
@@ -512,7 +474,7 @@ function segment_closest_point(seg,pt) =
 //   line_from_points(points, [fast], [eps]);
 // Description:
 //   Given a list of 2 or more collinear points, returns a line containing them.
-//   If `fast` is false and the points are coincident, then `undef` is returned.
+//   If `fast` is false and the points are coincident or non-collinear, then `undef` is returned.
 //   if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
 // Arguments:
 //   points = The list of points to find the line through.
@@ -522,7 +484,7 @@ function line_from_points(points, fast=false, eps=EPSILON) =
     assert( is_path(points,dim=undef), "Improper point list." )
     assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
     let( pb = furthest_point(points[0],points) )
-    approx(norm(points[pb]-points[0]),0) ? undef :
+    norm(points[pb]-points[0])<eps*max(norm(points[pb]),norm(points[0])) ? undef :
     fast || collinear(points) ? [points[pb], points[0]] : undef;
 
 
@@ -1171,7 +1133,7 @@ function _general_plane_line_intersection(plane, line, eps=EPSILON) =
 // Description:
 //   Returns a new representation [A,B,C,D] of `plane` where norm([A,B,C]) is equal to one.
 function normalize_plane(plane) =
-    assert( _valid_plane(plane), "Invalid plane." )
+    assert( _valid_plane(plane), str("Invalid plane. ",plane ) )
     plane/norm(point3d(plane));
 
 
@@ -1179,12 +1141,12 @@ function normalize_plane(plane) =
 // Usage:
 //   angle = plane_line_angle(plane,line);
 // Description:
-//   Compute the angle between a plane [A, B, C, D] and a line, specified as a pair of points [p1,p2].
+//   Compute the angle between a plane [A, B, C, D] and a 3d line, specified as a pair of 3d points [p1,p2].
 //   The resulting angle is signed, with the sign positive if the vector p2-p1 lies on
 //   the same side of the plane as the plane's normal vector.
 function plane_line_angle(plane, line) =
     assert( _valid_plane(plane), "Invalid plane." )
-    assert( _valid_line(line), "Invalid line." )
+    assert( _valid_line(line,dim=3), "Invalid 3d line." )
     let(
         linedir   = unit(line[1]-line[0]),
         normal    = plane_normal(plane),
@@ -1209,7 +1171,7 @@ function plane_line_angle(plane, line) =
 //   eps = Tolerance in geometric comparisons.  Default: `EPSILON` (1e-9)
 function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
     assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
-    assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.")
+    assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.")
     assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
@@ -1240,7 +1202,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
     assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
     assert(is_path(poly,dim=3), "Invalid polygon." )
     assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).")
-    assert(_valid_line(line,dim=3,eps=eps), "Invalid line." )
+    assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." )
     let(
         bounded = is_list(bounded)? bounded : [bounded, bounded],
         poly = deduplicate(poly),
@@ -1694,7 +1656,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
 //   eps = epsilon used for identifying the case with one solution.  Default: 1e-9
 function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) =
   let(r=get_radius(r=r,d=d,dflt=undef))
-  assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.")
+  assert(_valid_line(line,2), "Invalid 2d line.")
   assert(is_vector(c,2), "Circle center must be a 2-vector")
   assert(is_num(r) && r>0, "Radius must be positive")
   assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition")
@@ -1738,14 +1700,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
         pb = points[b],
         nrm = norm(pa-pb)
         )
-    approx(nrm, 0)
+    nrm <= eps*max(norm(pa),norm(pb))
     ? assert(!error, "Cannot find three noncollinear points in pointlist.")
         []
     :   let(
             n = (pb-pa)/nrm,
             distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]
            )
-        max(distlist)<eps*nrm
+        max(distlist) < eps*nrm
         ?  assert(!error, "Cannot find three noncollinear points in pointlist.")
            []
         :  [0,b,max_index(distlist)];
@@ -1761,12 +1723,13 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
 // Arguments:
 //   pts = List of points.
 function pointlist_bounds(pts) =
-    assert(is_matrix(pts) && len(pts)>0 && len(pts[0])>0 , "Invalid pointlist." )
-    let(ptsT = transpose(pts))
-    [
-      [for(row=ptsT) min(row)],
-      [for(row=ptsT) max(row)]
-    ];
+    assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." )
+    let( 
+        select = ident(len(pts[0])),
+        spread = [for(i=[0:len(pts[0])-1])
+                      let( spreadi = pts*select[i] ) 
+                      [min(spreadi), max(spreadi)] ] )
+    transpose(spread);
 
 
 // Function: closest_point()
@@ -1805,7 +1768,7 @@ function furthest_point(pt, points) =
 //   area = polygon_area(poly);
 // Description:
 //   Given a 2D or 3D planar polygon, returns the area of that polygon.
-//   If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is [].
+//   If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`.
 //   When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon.
 // Arguments:
 //   poly = Polygon to compute the area of.
@@ -1817,53 +1780,17 @@ function polygon_area(poly, signed=false) =
       ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 )
         signed ? total : abs(total)
       : let( plane = plane_from_polygon(poly) )
-        plane==[]? [] :
+        plane==[]? undef :
         let(
             n = plane_normal(plane),  
-            total = sum([
-                for(i=[1:1:len(poly)-2])
-                    let(
-                        v1 = poly[i] - poly[0],
-                        v2 = poly[i+1] - poly[0]
-                    )
-                    cross(v1,v2) 
-                ])* n/2
+            total = 
+								sum([ for(i=[1:1:len(poly)-2])
+												cross(poly[i]-poly[0], poly[i+1]-poly[0]) 
+										]) * n/2
         ) 
         signed ? total : abs(total);
 
 
-// Function: is_convex_polygon()
-// Usage:
-//   is_convex_polygon(poly);
-// Description:
-//   Returns true if the given 2D or 3D polygon is convex.  
-//   The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar.
-//   If the points are collinear an error is generated.
-// Arguments:
-//   poly = Polygon to check.
-//   eps = Tolerance for the collinearity test. Default: EPSILON.
-// Example:
-//   is_convex_polygon(circle(d=50));  // Returns: true
-//   is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
-// Example:
-//   spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
-//   is_convex_polygon(spiral);  // Returns: false
-function is_convex_polygon(poly,eps=EPSILON) =
-    assert(is_path(poly), "The input should be a 2D or 3D polygon." )
-    let( lp = len(poly),
-         p0 = poly[0] )
-    assert( lp>=3 , "A polygon must have at least 3 points" )
-    let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
-    len(p0)==2
-    ?   assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" )
-        min(crosses) >=0 || max(crosses)<=0
-    :   let( prod = crosses*sum(crosses),
-             minc = min(prod),
-             maxc = max(prod) )
-        assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" )
-        minc>=0 || maxc<=0;
-
-
 // Function: polygon_shift()
 // Usage:
 //   polygon_shift(poly, i);
@@ -2030,9 +1957,9 @@ function centroid(poly, eps=EPSILON) =
 //   Returns -1 if the point is outside the polygon.
 //   Returns 0 if the point is on the boundary.
 //   Returns 1 if the point lies in the interior.
-//   The polygon does not need to be simple: it can have self-intersections.
+//   The polygon does not need to be simple: it may have self-intersections.
 //   But the polygon cannot have holes (it must be simply connected).
-//   Rounding error may give mixed results for points on or near the boundary.
+//   Rounding errors may give mixed results for points on or near the boundary.
 // Arguments:
 //   point = The 2D point to check position of.
 //   poly = The list of 2D path points forming the perimeter of the polygon.
@@ -2125,7 +2052,7 @@ function ccw_polygon(poly) =
 //   poly = The list of the path points for the perimeter of the polygon.
 function reverse_polygon(poly) =
     assert(is_path(poly), "Input should be a polygon")
-    let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]];
+   [poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ]; 
 
 
 // Function: polygon_normal()
@@ -2133,7 +2060,7 @@ function reverse_polygon(poly) =
 //   n = polygon_normal(poly);
 // Description:
 //   Given a 3D planar polygon, returns a unit-length normal vector for the
-//   clockwise orientation of the polygon. If the polygon points are collinear, returns [].
+//   clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`.
 //   It doesn't check for coplanarity.
 // Arguments:
 //   poly = The list of 3D path points for the perimeter of the polygon.
@@ -2141,7 +2068,7 @@ function polygon_normal(poly) =
     assert(is_path(poly,dim=3), "Invalid 3D polygon." )
     len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) :
     let( triple = sort(noncollinear_triple(poly,error=false)) )
-    triple==[] ? [] :
+    triple==[] ? undef :
     point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ;
 
 
@@ -2295,4 +2222,255 @@ function split_polygons_at_each_z(polys, zs, _i=0) =
     );
 
 
+
+// Section: Convex Sets
+
+
+// Function: is_convex_polygon()
+// Usage:
+//   is_convex_polygon(poly);
+// Description:
+//   Returns true if the given 2D or 3D polygon is convex.  
+//   The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar.
+//   If the points are collinear an error is generated.
+// Arguments:
+//   poly = Polygon to check.
+//   eps = Tolerance for the collinearity test. Default: EPSILON.
+// Example:
+//   is_convex_polygon(circle(d=50));  // Returns: true
+//   is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
+// Example:
+//   spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
+//   is_convex_polygon(spiral);  // Returns: false
+function is_convex_polygon(poly,eps=EPSILON) =
+    assert(is_path(poly), "The input should be a 2D or 3D polygon." )
+    let( lp = len(poly),
+         p0 = poly[0] )
+    assert( lp>=3 , "A polygon must have at least 3 points" )
+    let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
+    len(p0)==2
+    ?   assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" )
+        min(crosses) >=0 || max(crosses)<=0
+    :   let( prod = crosses*sum(crosses),
+             minc = min(prod),
+             maxc = max(prod) )
+        assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" )
+        minc>=0 || maxc<=0;
+
+
+// Function: convex_distance()
+// Usage:
+//   convex_distance(pts1, pts2,<eps=>);
+// See also: 
+//   convex_collision
+// Description:
+//   Returns the smallest distance between a point in convex hull of `points1` 
+//   and a point in the convex hull of `points2`. All the points in the lists
+//   should have the same dimension, either 2D or 3D. 
+//   A zero result means the hulls intercept whithin a tolerance `eps`.
+// Arguments:
+//   points1 - first list of 2d or 3d points.
+//   points2 - second list of 2d or 3d points.
+//   eps - tolerance in distance evaluations. Default: EPSILON.
+// Example(2D):
+//    pts1 = move([-3,0], p=square(3,center=true));
+//    pts2 = rot(a=45, p=square(2,center=true));
+//    pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
+//    polygon(pts1);
+//    polygon(pts2);
+//    polygon(pts3);
+//    echo(convex_distance(pts1,pts2)); // Returns: 0.0857864
+//    echo(convex_distance(pts2,pts3)); // Returns: 0
+// Example(3D):
+//    sphr1 = sphere(2,$fn=10);
+//    sphr2 = move([4,0,0], p=sphr1);
+//    sphr3 = move([4.5,0,0], p=sphr1);    
+//    vnf_polyhedron(sphr1);
+//    vnf_polyhedron(sphr2);
+//    echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0
+//    echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5
+function convex_distance(points1, points2, eps=EPSILON) =
+    assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), 
+           "The input list should be a consistent non empty list of points of same dimension.")
+    assert(len(points1[0])==2 || len(points1[0])==3 ,
+           "The input points should be 2d or 3d points.")
+    let( d = points1[0]-points2[0] )
+    norm(d)<eps ? 0 :
+    let( v = _support_diff(points1,points2,-d) ) 
+    norm(_GJK_distance(points1, points2, eps, 0, v, [v]));
+
+
+// Finds the vector difference between the hulls of the two pointsets by the GJK algorithm
+// Based on:
+// http://www.dtecta.com/papers/jgt98convex.pdf
+function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) = 
+    let( nrd = norm(d) ) // distance upper bound
+    nrd<eps ? d :
+    let(
+        v     = _support_diff(points1,points2,-d),
+        lbd   = max(lbd, d*v/nrd), // distance lower bound
+        close = (nrd-lbd <= eps*nrd) 
+    )
+    // v already in the simplex is a degenerence due to numerical errors 
+    // and may produce a non-stopping loop
+    close || [for(nv=norm(v), s=simplex) if(norm(s-v)<=eps*nv) 1]!=[] ? d :  
+    let( newsplx = _closest_simplex(concat(simplex,[v]),eps) ) 
+    _GJK_distance(points1, points2, eps, lbd, newsplx[0], newsplx[1]);
+
+
+// Function: convex_collision()
+// Usage:
+//   convex_collision(pts1, pts2,<eps=>);
+// See also: 
+//   convex_distance
+// Description:
+//   Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2`
+//   otherwise, `false`.
+//   All the points in the lists should have the same dimension, either 2D or 3D.
+//   This function is tipically faster than `convex_distance` to find a non-collision.
+// Arguments:
+//   points1 - first list of 2d or 3d points.
+//   points2 - second list of 2d or 3d points.
+//   eps - tolerance for the intersection tests. Default: EPSILON.
+// Example(2D):
+//    pts1 = move([-3,0], p=square(3,center=true));
+//    pts2 = rot(a=45, p=square(2,center=true));
+//    pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
+//    polygon(pts1);
+//    polygon(pts2);
+//    polygon(pts3);
+//    echo(convex_collision(pts1,pts2)); // Returns: false
+//    echo(convex_collision(pts2,pts3)); // Returns: true
+// Example(3D):
+//    sphr1 = sphere(2,$fn=10);
+//    sphr2 = move([4,0,0], p=sphr1);
+//    sphr3 = move([4.5,0,0], p=sphr1);    
+//    vnf_polyhedron(sphr1);
+//    vnf_polyhedron(sphr2);
+//    echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true
+//    echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false
+//
+function convex_collision(points1, points2, eps=EPSILON) =
+    assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])), 
+           "The input list should be a consistent non empty list of points of same dimension.")
+    assert(len(points1[0])==2 || len(points1[0])==3 ,
+           "The input points should be 2d or 3d points.")
+    let( d = points1[0]-points2[0] )
+    norm(d)<eps ? true :
+    let( v = _support_diff(points1,points2,-d) ) 
+    _GJK_collide(points1, points2, v, [v], eps);
+
+
+// Based on the GJK collision algorithms found in:
+// http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
+// or
+// http://www.dtecta.com/papers/jgt98convex.pdf
+function _GJK_collide(points1, points2, d, simplex, eps=EPSILON) = 
+    norm(d) < eps ? true :          // does collide 
+    let( v = _support_diff(points1,points2,-d) )
+    v*d > eps ? false : // no collision
+    let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
+    _GJK_collide(points1, points2, newsplx[0], newsplx[1], eps);
+
+
+// given a simplex s, returns a pair:
+//  - the point of the s closest to the origin
+//  - the smallest sub-simplex of s that contains that point
+function _closest_simplex(s,eps=EPSILON) =
+    assert(len(s)>=2 && len(s)<=4, "Internal error.")
+    len(s)==2 ? _closest_s1(s,eps) : 
+    len(s)==3 ? _closest_s2(s,eps)
+              : _closest_s3(s,eps);
+
+
+// find the closest to a 1-simplex
+// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
+function _closest_s1(s,eps=EPSILON) =
+    norm(s[1]-s[0])<eps*(norm(s[0])+norm(s[1]))/2 ? [ s[0], [s[0]] ] :
+    let( 
+        c = s[1]-s[0],
+        t = -s[0]*c/(c*c)
+    ) 
+    t<0 ? [ s[0], [s[0]] ] :
+    t>1 ? [ s[1], [s[1]] ] :
+    [ s[0]+t*c, s ];
+
+
+// find the closest to a 2-simplex
+// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
+function _closest_s2(s,eps=EPSILON) =
+    let(
+        dim = len(s[0]),
+        a  = dim==3 ? s[0]: [ each s[0], 0] ,
+        b  = dim==3 ? s[1]: [ each s[1], 0] ,
+        c  = dim==3 ? s[2]: [ each s[2], 0] ,
+        ab = norm(a-b),
+        bc = norm(b-c),
+        ca = norm(c-a),
+        nr = cross(b-a,c-a)
+    )
+    norm(nr) <= eps*max(ab,bc,ca) // degenerate case
+    ?   let( i = max_index([ab, bc, ca]) )
+        _closest_s1([s[i],s[(i+1)%3]],eps)
+// considering that s[2] was the last inserted vertex in s,
+// the only possible outcomes are :
+//    s, [s[0],s[2]] and [s[1],s[2]]
+    :   let(
+            class =   (cross(nr,a-b)*a<0 ? 1 : 0 )
+                    + (cross(nr,c-a)*a<0 ? 2 : 0 )
+                    + (cross(nr,b-c)*b<0 ? 4 : 0 )
+        )
+        assert( class!=1, "Internal error" )  
+        class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri
+//        class==1 ? _closest_s1([s[0],s[1]]) :
+        class==2 ? _closest_s1([s[0],s[2]],eps) :  
+        class==4 ? _closest_s1([s[1],s[2]],eps) :
+//        class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) : 
+        class==3 ? _closest_s1([s[0],s[2]],eps) : 
+//        class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) :
+        class==5 ? _closest_s1([s[1],s[2]],eps) :
+        c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps);
+
+
+// find the closest to a 3-simplex
+// it seems that degenerate 3-simplices are correctly manage without extra code
+function _closest_s3(s,eps=EPSILON) =
+    assert( len(s[0])==3 && len(s)==4, "Internal error." )
+    let( nr = cross(s[1]-s[0],s[2]-s[0]),
+         sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] )
+    norm(nr)<eps*max(sz) 
+    ?   let( i = max_index(sz) )
+        _closest_s2([ s[i], s[(i+1)%3], s[3] ], eps) // degenerate case
+    // considering that s[3] was the last inserted vertex in s,
+    // the only possible outcomes will be:
+    //    s or some of the 3 triangles of s containing s[3]
+    :   let(
+            tris = [ [s[0], s[1], s[3]],
+                     [s[1], s[2], s[3]],
+                     [s[2], s[0], s[3]] ],
+            cntr = sum(s)/4,
+            // indicator of the tris facing the origin
+            facing = [for(i=[0:2])
+                        let( nrm = _tri_normal(tris[i]) )
+                        if( ((nrm*(s[i]-cntr))>0)==(nrm*s[i]<0) ) i ]
+        )
+        len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex
+        len(facing)==1 ? _closest_s2(tris[facing[0]], eps) :
+        let( // look for the origin-facing tri closest to the origin
+            closest = [for(i=facing) _closest_s2(tris[i], eps) ],
+            dist    = [for(cl=closest) norm(cl[0]) ],
+            nearest = min_index(dist)
+        )
+        closest[nearest];
+
+
+function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]);
+
+
+function _support_diff(p1,p2,d) =
+    let( p1d = p1*d, p2d = p2*d )
+    p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]];
+
+
+
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/hull.scad b/hull.scad
index 1dc147c..4adb7c3 100644
--- a/hull.scad
+++ b/hull.scad
@@ -202,17 +202,15 @@ function hull3d_faces(points) =
 
 
 // Adds the remaining points one by one to the convex hull
-function _hull3d_iterative(points, triangles, planes, remaining, _i=0) =
+function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = //let( EPSILON=1e-12 )
     _i >= len(remaining) ? triangles : 
     let (
         // pick a point
         i = remaining[_i],
         // evaluate the triangle plane equations at point i
-//        xx=[for(i=[0:len(planes)-1],p=[planes[i]]) if(len(p)!=4) echo(i=i,len(p))0 ],//echo([each points[i], -1]),
-//        planeq_val = [for(p=planes) p*[each points[i], -1]],
         planeq_val = planes*[each points[i], -1],
         // find the triangles that are in conflict with the point (point not inside)
-        conflicts = [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]>EPSILON) i ],
+        conflicts = [for (i = [0:1:len(planeq_val)-1]) if (planeq_val[i]>EPSILON) i ],
         // collect the halfedges of all triangles that are in conflict 
         halfedges = [ 
             for(c = conflicts, i = [0:2])
@@ -222,18 +220,15 @@ function _hull3d_iterative(points, triangles, planes, remaining, _i=0) =
         horizon = _remove_internal_edges(halfedges),
         // generate new triangles connecting point i to each horizon halfedge vertices
         tri2add = [ for (h = horizon) concat(h,i) ],
-//        w=[for(t=tri2add) if(collinear(points[t[0]],points[t[1]],points[t[2]])) echo(t)],
         // add tria2add and remove conflict triangles
         new_triangles = 
             concat( tri2add,
                     [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) triangles[i] ] 
                   ),
-//        y=[for(t=tri2add) if([]==plane3pt_indexed(points,t[0],t[1],t[2])) echo(tri2add=t,pts=[for(ti=t) points[ti]])],
         // add the plane equations of new added triangles and remove the plane equations of the conflict ones
         new_planes = 
-            concat( [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ],
-                    [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ] 
-                  )
+            [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ,
+              for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ] 
     ) _hull3d_iterative(
         points,
         new_triangles,
@@ -249,7 +244,6 @@ function _remove_internal_edges(halfedges) = [
             h
 ];
 
-
 function _find_first_noncoplanar(plane, points, i=0) = 
     (i >= len(points) || !points_on_plane([points[i]],plane))? i :
     _find_first_noncoplanar(plane, points, i+1);
diff --git a/math.scad b/math.scad
index 96bd475..30b585e 100644
--- a/math.scad
+++ b/math.scad
@@ -640,17 +640,19 @@ function sum_of_sines(a, sines) =
 // Usage:
 //   delts = deltas(v);
 // Description:
-//   Returns a list with the deltas of adjacent entries in the given list.
+//   Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front.
 //   The list should be a consistent list of numeric components (numbers, vectors, matrix, etc).
 //   Given [a,b,c,d], returns [b-a,c-b,d-c].
+//   
 // Arguments:
 //   v = The list to get the deltas of.
+//   wrap = If true, wrap back to the start from the end.  ie: return the difference between the last and first items as the last delta.  Default: false
 // Example:
 //   deltas([2,5,9,17]);  // returns [3,4,8].
 //   deltas([[1,2,3], [3,6,8], [4,8,11]]);  // returns [[2,4,5], [1,2,3]]
-function deltas(v) = 
+function deltas(v, wrap=false) = 
     assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
-    [for (p=pair(v)) p[1]-p[0]] ;
+    [for (p=pair(v,wrap)) p[1]-p[0]] ;
 
 
 // Function: product()
@@ -771,21 +773,31 @@ function _med3(a,b,c) =
 // Usage:
 //   x = convolve(p,q);
 // Description:
-//   Given two vectors, finds the convolution of them.
-//   The length of the returned vector is len(p)+len(q)-1 .
+//   Given two vectors, or one vector and a path or
+//   two paths of the same dimension, finds the convolution of them.
+//   If both parameter are vectors, returns the vector convolution.
+//   If one parameter is a vector and the other a path,
+//   convolves using products by scalars and returns a path. 
+//   If both parameters are paths, convolve using scalar products
+//   and returns a vector.
+//   The returned vector or path has length len(p)+len(q)-1.
 // Arguments:
-//   p = The first vector.
-//   q = The second vector.
+//   p = The first vector or path.
+//   q = The second vector or path.
 // Example:
 //   a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1]
 //   b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3]
+//   c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]]
+//   d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns:  [3,9,11,7]
 function convolve(p,q) =
     p==[] || q==[] ? [] :
-    assert( is_vector(p) && is_vector(q), "The inputs should be vectors.")
+    assert( (is_vector(p) || is_matrix(p))
+		        && ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) ,
+		        "The inputs should be vectors or paths all of the same dimension.")
     let( n = len(p),
          m = len(q))
     [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
-       [for(j=[k1:k2]) p[i-j] ] * [for(j=[k1:k2]) q[j] ] 
+       sum([for(j=[k1:k2]) p[i-j]*q[j] ]) 
     ];
 
 
@@ -1694,7 +1706,7 @@ function polynomial(p,z,k,total) =
 //   x = polymult(p,q)
 //   x = polymult([p1,p2,p3,...])
 // Description:
-//   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
+//   Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first, 
 //   computes the coefficient list of the product polynomial.  
 function poly_mult(p,q) = 
   is_undef(q) ?
@@ -1714,8 +1726,8 @@ function poly_mult(p,q) =
 // Description:
 //    Computes division of the numerator polynomial by the denominator polynomial and returns
 //    a list of two polynomials, [quotient, remainder].  If the division has no remainder then
-//    the zero polynomial [] is returned for the remainder.  Similarly if the quotient is zero
-//    the returned quotient will be [].  
+//    the zero polynomial [0] is returned for the remainder.  Similarly if the quotient is zero
+//    the returned quotient will be [0].  
 function poly_div(n,d) =
     assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
     let( d = _poly_trim(d), 
@@ -1740,7 +1752,7 @@ function _poly_div(n,d,q) =
 ///    _poly_trim(p,[eps])
 /// Description:
 ///    Removes leading zero terms of a polynomial.  By default zeros must be exact,
-///    or give epsilon for approximate zeros.  
+///    or give epsilon for approximate zeros. Returns [0] for a zero polynomial.
 function _poly_trim(p,eps=0) =
     let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
     len(nz)==0 ? [0] : list_tail(p,nz[0]);
diff --git a/paths.scad b/paths.scad
index 88d9ec7..3c79b97 100644
--- a/paths.scad
+++ b/paths.scad
@@ -997,7 +997,7 @@ module jittered_poly(path, dist=1/512) {
 
 // Module: extrude_from_to()
 // Description:
-//   Extrudes a 2D shape between the points pt1 and pt2.  Takes as children a set of 2D shapes to extrude.
+//   Extrudes a 2D shape between the 3d points pt1 and pt2.  Takes as children a set of 2D shapes to extrude.
 // Arguments:
 //   pt1 = starting point of extrusion.
 //   pt2 = ending point of extrusion.
@@ -1010,6 +1010,7 @@ module jittered_poly(path, dist=1/512) {
 //       xcopies(3) circle(3, $fn=32);
 //   }
 module extrude_from_to(pt1, pt2, convexity, twist, scale, slices) {
+		assert( is_path([pt1,pt2],3), "The points should be 3d points");
     rtp = xyz_to_spherical(pt2-pt1);
     translate(pt1) {
         rotate([0, rtp[2], rtp[1]]) {
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 8569d7f..a66aeed 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -358,11 +358,21 @@ module test_sortidx() {
 }
 test_sortidx();
 
+module test_group_sort() {
+    assert_equal(group_sort([]), [[]]);
+    assert_equal(group_sort([8]), [[8]]);
+    assert_equal(group_sort([7,3,9,4,3,1,8]), [[1], [3, 3], [4], [7], [8], [9]]);
+    assert_equal(group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0), [[[2, "b"], [2, "e"]], [[3, "d"]], [[5, "a"], [5, "c"]]]);
+    assert_equal(group_sort([["a",5],["b",6], ["c",1], ["d",2], ["e",6] ], idx=1), [[["c", 1]], [["d", 2]], [["a", 5]], [["b", 6], ["e", 6]]] );
+}
+test_group_sort();
+
 
 module test_unique() {
-    assert(unique([]) == []);
-    assert(unique([8]) == [8]);
-    assert(unique([7,3,9,4,3,1,8]) == [1,3,4,7,8,9]);
+    assert_equal(unique([]), []);
+    assert_equal(unique([8]), [8]);
+    assert_equal(unique([7,3,9,4,3,1,8]), [1,3,4,7,8,9]);
+    assert_equal(unique(["A","B","R","A","C","A","D","A","B","R","A"]), ["A", "B", "C", "D", "R"]);
 }
 test_unique();
 
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 9cf7c40..992a681 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -90,7 +90,8 @@ test_cleanup_path();
 test_simplify_path();
 test_simplify_path_indexed();
 test_is_region();
-
+test_convex_distance();
+test_convex_collision();
 
 // to be used when there are two alternative symmetrical outcomes 
 // from a function like a plane output; v must be a vector
@@ -301,7 +302,7 @@ module test_line_from_points() {
 }
 *test_line_from_points();
 
-module test_point_on_segment2d() {
+module test_point_on_segment2d() { 
     assert(point_on_segment2d([-15,0], [[-10,0], [10,0]]) == false);
     assert(point_on_segment2d([-10,0], [[-10,0], [10,0]]) == true);
     assert(point_on_segment2d([-5,0], [[-10,0], [10,0]]) == true);
@@ -1075,6 +1076,44 @@ module test_is_region() {
 }
 *test_is_region();
 
+module test_convex_distance() {
+// 2D
+    c1 = circle(10,$fn=24);
+    c2 = move([15,0], p=c1);
+    assert(convex_distance(c1, c2)==0);
+    c3 = move([22,0],c1);
+    assert_approx(convex_distance(c1, c3),2);
+// 3D
+    s1 = sphere(10,$fn=4);
+    s2 = move([15,0], p=s1);
+    assert_approx(convex_distance(s1[0], s2[0]), 0.857864376269);
+    s3 = move([25.3,0],s1);
+    assert_approx(convex_distance(s1[0], s3[0]), 11.1578643763);
+    s4 = move([30,25],s1);    
+    assert_approx(convex_distance(s1[0], s4[0]), 28.8908729653);
+    s5 = move([10*sqrt(2),0],s1);    
+    assert_approx(convex_distance(s1[0], s5[0]), 0);
+}
+*test_convex_distance();
 
+module test_convex_collision() {
+// 2D
+    c1 = circle(10,$fn=24);
+    c2 = move([15,0], p=c1);
+    assert(convex_collision(c1, c2));
+    c3 = move([22,0],c1);
+    assert(!convex_collision(c1, c3));
+// 3D
+    s1 = sphere(10,$fn=4);
+    s2 = move([15,0], p=s1);
+    assert(!convex_collision(s1[0], s2[0]));
+    s3 = move([25.3,0],s1);
+    assert(!convex_collision(s1[0], s3[0]));
+    s4 = move([5,0],s1);    
+    assert(convex_collision(s1[0], s4[0]));
+    s5 = move([10*sqrt(2),0],s1);    
+    assert(convex_collision(s1[0], s5[0]));
+}
+*test_convex_distance();
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/test_math.scad b/tests/test_math.scad
index adcfb6d..b5c25b8 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -546,7 +546,9 @@ test_sum_of_sines();
 
 module test_deltas() {
     assert_equal(deltas([2,5,9,17]), [3,4,8]);
+    assert_equal(deltas([2,5,9,17],wrap=true), [3,4,8,-15]);
     assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]]), [[2,4,5], [1,2,3]]);
+    assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]],wrap=true), [[2,4,5], [1,2,3], [-3,-6,-8]]);
 }
 test_deltas();
 
@@ -582,6 +584,10 @@ module test_convolve() {
     assert_equal(convolve([1,1],[]), []);
     assert_equal(convolve([1,1],[1,2,1]), [1,3,3,1]);
     assert_equal(convolve([1,2,3],[1,2,1]), [1,4,8,8,3]);
+    assert_equal(convolve([1,2,3],[[1],[2],[1]]),  [[1], [4], [8], [8], [3]]);
+    assert_equal(convolve([[1],[2],[3]],[[1],[2],[1]]), [1,4,8,8,3]);
+    assert_equal(convolve([[1,0],[2,1],[3,2]],[[1,0],[2,1],[1,2]]), [1,4,9,12,7]);
+    assert_equal(convolve([1,2,3],[[1,0],[2,1],[1,2]]), [[1,0],[4,1],[8,4],[8,7],[3,6]]);
 }
 test_convolve();