From 9e4d46eecf02e6a30174962ccd1fd12bfeffd85e Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Thu, 21 Oct 2021 15:33:25 -0400
Subject: [PATCH] Add median function

---
 arrays.scad | 24 ++++++++++++++++++++++++
 math.scad   | 34 ++++++++++++++--------------------
 2 files changed, 38 insertions(+), 20 deletions(-)

diff --git a/arrays.scad b/arrays.scad
index 1e0d44e..0e3ab97 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -1244,6 +1244,30 @@ function group_sort(list, idx) =
         
 
 
+// Function: list_smallest()
+// Usage:
+//   small = list_smallest(list, k)
+// Description:
+//   Returns a set of the k smallest items in list in arbitrary order.  The items must be
+//   mutually comparable with native OpenSCAD comparison operations.
+// Arguments:
+//   list = list to process
+//   k = number of items to return
+function list_smallest(list, k) =
+    assert(is_list(list))
+    assert(is_finite(k) && k>=0, "k must be nonnegative")
+    let( 
+        v       = list[rand_int(0,len(list)-1,1)[0]],
+        smaller = [for(li=list) if(li<v) li ],
+        equal   = [for(li=list) if(li==v) li ]
+    )
+    len(smaller)   == k ? smaller :
+    len(smaller)<k && len(smaller)+len(equal) >= k ? [ each smaller, for(i=[1:k-len(smaller)]) v ] :
+    len(smaller)   >  k ? list_smallest(smaller, k) :
+    let( bigger  = [for(li=list) if(li>v) li ] )
+    concat(smaller, equal, list_smallest(bigger, k-len(smaller) -len(equal)));
+
+
 // Function: group_data()
 // Usage:
 //   groupings = group_data(groups, values);
diff --git a/math.scad b/math.scad
index 9d57201..b50e230 100644
--- a/math.scad
+++ b/math.scad
@@ -827,27 +827,21 @@ function mean(v) =
     sum(v)/len(v);
 
 
-// Function: ninther()
-// Usage:
-//    med = ninther(v)
-// Description:
-//    Finds a value in the input list of numbers `v` that is the median of  a 
-//    sample of 9 entries of `v`.
-//    It is a much faster approximation of the true median computation.
-// Arguments:
-//    v = an array of numbers
-function ninther(v) = 
-    let( l=len(v) )
-    l<=4 ? l<=2 ? v[0] : _med3(v[0], v[1], v[2]) : 
-    l==5 ? _med3(v[0], _med3(v[1], v[2], v[3]), v[4]) :
-    _med3(_med3(v[0],v[floor(l/6)],v[floor(l/3)]),
-          _med3(v[floor(l/3)],v[floor(l/2)],v[floor(2*l/3)]),
-          _med3(v[floor(2*l/3)],v[floor((5*l/3 -1)/2)],v[l-1]) );
 
-// the median of a triple
-function _med3(a,b,c) =
-    a < c ? a < b ? min(b,c) : min(a,c) :
-    b < c ? min(a,c) : min(a,b);
+// Function: median()
+// Usage:
+//   middle = median(v)
+// Description:
+//   Returns the median of the given vector.  
+function median(v) =
+    assert(is_vector(v), "Input to median must be a vector")
+    len(v)%2 ? max( list_smallest(v, ceil(len(v)/2)) ) :
+    let( lowest = list_smallest(v, len(v)/2 + 1),
+         max  = max(lowest),
+         imax = search(max,lowest,1),
+         max2 = max([for(i=idx(lowest)) if(i!=imax[0]) lowest[i] ])
+    )
+    (max+max2)/2;
 
 
 // Function: convolve()