Index: /trunk/Ohana/src/opihi/cmd.data/stats-new.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/stats-new.c	(revision 15274)
+++ /trunk/Ohana/src/opihi/cmd.data/stats-new.c	(revision 15274)
@@ -0,0 +1,259 @@
+# include "data.h"
+
+double interpolateValue (int *hist, int Nbins, float min, float binsize, int bin, float value);
+double sampleMedian (Matrix *matrix, int Nsubset);
+Stats *robustMedian (Matrix *matrix, Region *region, float minValue, float maxValue, float sigma);
+
+typedef struct {
+  int sx;
+  int sy;
+  int nx;
+  int ny;
+} Region;
+
+typedef struct {
+  float median;
+  float mean;
+  float mode;
+  float sigma;
+} Stats;
+
+int stats (int argc, char **argv) {
+  
+  int i, j, Nmode, Imode;
+  double Npix, N1, N2, max, min, range, median, mode, IgnoreValue;
+  float *V;
+  int sx, sy, nx, ny, *hist, Nhist, bin;
+  int Ignore, Quiet, N;
+  Buffer *buf;
+
+  IgnoreValue = 0;
+  Ignore = FALSE;
+  if ((N = get_argument (argc, argv, "-ignore"))) {
+    Ignore = TRUE;
+    remove_argument (N, &argc, argv);
+    IgnoreValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  Quiet = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((argc != 2) && (argc != 6)) {
+    gprint (GP_ERR, "USAGE: stats <buffer> sx sy nx ny\n");
+    gprint (GP_ERR, "OR:    stats <buffer>\n");
+    return (FALSE);
+  }
+
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+
+  if (argc == 6) {
+    sx = strcmp (argv[2], "-") ? atof (argv[2]) : 0;
+    sy = strcmp (argv[3], "-") ? atof (argv[3]) : 0;
+    nx = strcmp (argv[4], "-") ? atof (argv[4]) : buf[0].matrix.Naxis[0];
+    ny = strcmp (argv[5], "-") ? atof (argv[5]) : buf[0].matrix.Naxis[1];
+  } else {
+    sx = 0;
+    sy = 0;
+    nx = buf[0].matrix.Naxis[0];
+    ny = buf[0].matrix.Naxis[1];
+  }
+
+  Npix = N1 = N2 = 0;
+  if ((sx < 0) || (sy < 0) || 
+      (sx+nx > buf[0].matrix.Naxis[0]) || 
+      (sy+ny > buf[0].matrix.Naxis[1])) {
+    gprint (GP_ERR, "region out of range\n");
+    return (FALSE);
+  }
+
+  max = min = *((float *)(buf[0].matrix.buffer) + sy*buf[0].matrix.Naxis[0] + sx);
+  for (j = sy; j < sy + ny; j++) {
+    V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + sx; 
+    for (i = 0; i < nx; i++, V++) {
+      if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue;
+      N1 += *V;
+      N2 += (*V)*(*V);
+      Npix += 1.0;
+      max = MAX (max, *V);
+      min = MIN (min, *V);
+    }
+  }
+  N1 = N1 / Npix;
+
+  sigma = sqrt (N2/Npix - N1*N1);
+
+  if (ROBUST) {
+    stats = robustMedian (matrix, min, max, sigma);
+    if (stats->sigma 
+  } 
+  
+  if (!Quiet) {
+    gprint (GP_LOG, "  mean    stdev    min     max   median   Npix   Total\n");
+    gprint (GP_LOG, "%7.4g %7.4g %7.4g %7.4g %7.4g %7.0f %7.4g\n", N1, sqrt (N2/Npix - N1*N1), 
+	    min, max, median, Npix, Npix*N1);
+  }
+
+  set_variable ("MIN",    min);
+  set_variable ("MAX",    max);
+  set_variable ("MEDIAN", median);
+  set_variable ("MEAN",   N1);
+  set_variable ("MODE",   mode);
+  set_variable ("TOTAL",  N1*Npix);
+  set_variable ("NPIX",   Npix);
+  set_variable ("SIGMA",  sqrt (N2/Npix - N1*N1));
+
+  return (TRUE);
+}
+
+double sampleMedian (Matrix *matrix, int Nsubset) {
+
+  int i, Nsample, Npix;
+  long A, B;
+  float *values, *buffer, median;
+
+  /* Generate a vector of Nsample elements */
+  Npix = matrix[0].Naxis[0]*matrix[0].Naxis[1]
+    Nsample = Nsubset;
+  if (Nsubset == 0) {
+    Nsample = Npix;
+  }
+
+  ALLOCATE (values, float, Nsample);
+    
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
+ 
+  *buffer = (float *) matrix[0].buffer;
+
+  for (i = 0; i < Nsample; i++) {
+    if (Nsubset) {
+      j = MIN(Npix - 1, MAX(0, Npix*drand48()));
+    } else {
+      j = i;
+    }
+    values[i] = buffer[j];
+  }
+    
+  fsort (values, Nsample);
+
+  if (Nsample == 0) return 0.0;
+
+  // these are all covered by the logic below
+  // if (Nsample == 1) return values[0];
+  // if (Nsample == 2) return (0.5*(values[0] + values[1]));
+  // if (Nsample == 3) return (values[1]);
+
+  if (Nsample % 2) {
+    median = 0.5*(values[(int)(0.5*Nsample)] + values[(int)(0.5*Nsample) - 1]);
+  } else {
+    median = values[(int)(0.5*Nsample)];
+  }
+  return (median);
+}
+
+Stats *robustMedian (Matrix *matrix, Region *region, float minValue, float maxValue, float sigma) {
+
+  Stats *stats;
+
+  ALLOCATE (stats, Stats, 1);
+
+  // no data range:
+  if (maxValue - minValue < FLT_EPSILON) {
+    stats->median = minValue;
+    stats->mean   = minValue;
+    stats->mode   = minValue;
+    stats->sigma  = 0.0;
+    return stats;
+  }
+
+  // generate a histogram ranging from min to max with step size of sigma / 4
+  // no more than 0x1000 bins, no less than 0x10
+
+  Nbins = 4.0 * (maxValue - minValue) / sigma;
+  Nbins = MIN (0x1000, MAX (0x10, Nbins));
+
+  ALLOCATE (hist, int, Nbins);
+  memset (hist, 0, Nbins*sizeof(int));
+  
+  delta = Nbins / (max - min);
+  binsize = 1.0 / delta;
+
+  Ntotal = 0;
+  for (j = region->sy; j < region->sy + region->ny; j++) {
+    V = (float *)(buf[0].matrix.buffer) + j*buf[0].matrix.Naxis[0] + region->sx; 
+    for (i = 0; i < region->nx; i++, V++) {
+      if (Ignore && (fabs (*V - IgnoreValue) < 1e-8)) continue;
+      bin = MIN (MAX (0, (*V - min) * delta), Nbins - 1);
+      hist[bin] ++;
+      Ntotal ++;
+    }
+  }
+
+  // generate the cumulative histogram & find mode
+  ALLOCATE (cumu, int, Nbins);
+  memset (cumu, 0, Nbins*sizeof(int));
+
+  Nhist = 0;
+  Imode = -1;
+  Vmode = 0;
+  for (i = 0; i < Nbins; i++) {
+    Nhist += hist[i];
+    cumu[i] = Nhist;
+    if (hist[i] > Vmode) {
+      Imode = i;
+      Vmode = hist[i];
+    }
+  }
+
+  Ioo = ibracket (cumu, Nbins, 0.500000*Nhist, 1);
+  Im1 = ibracket (cumu, Nbins, 0.308538*Nhist, 1);
+  Ip1 = ibracket (cumu, Nbins, 0.691462*Nhist, 1);
+  Im2 = ibracket (cumu, Nbins, 0.022481*Nhist, 1);
+  Ip2 = ibracket (cumu, Nbins, 0.977519*Nhist, 1);
+
+  Voo = interpolateValue (cumu, Nbins, min, delta, Ioo, 0.500000*Nhist);
+  Vm1 = interpolateValue (cumu, Nbins, min, delta, Im1, 0.308538*Nhist);
+  Vp1 = interpolateValue (cumu, Nbins, min, delta, Ip1, 0.691462*Nhist);
+  Vm2 = interpolateValue (cumu, Nbins, min, delta, Im2, 0.022481*Nhist);
+  Vp2 = interpolateValue (cumu, Nbins, min, delta, Ip2, 0.977519*Nhist);
+
+  sigma1 = (Vp1 - Vm1);
+  sigma2 = (Vp2 - Vm2) / 4.0;
+
+  Vs = Vn = 0;
+  for (i = Im1; i < Ip1 + 1; i++) {
+    Vn += hist[i];
+    Vs += hist[i] * (i * binsize + min);
+  }
+
+  stats->median = Voo;
+  stats->mean = Vs / Vn;
+  stats->mode = Imode * delta + min;
+  stats->sigma = MIN (sigma1, sigma2);
+
+  free (hist);
+  free (cumu);
+
+  return (stats);
+}
+
+// what is the fractional bin position for the given value
+double interpolateValue (int *hist, int Nbins, float min, float binsize, int bin, float value) {
+
+  V0 = bin * binsize + min;
+  V1 = (bin + 1) * binsize + min;
+
+  // fbin = (value - V0) / (V1 - V0) + bin;
+  fbin = (value - V0) / binsize + bin;
+
+  return fbin;
+}
Index: /trunk/Ohana/src/opihi/dvo/subpix.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/subpix.c	(revision 15273)
+++ /trunk/Ohana/src/opihi/dvo/subpix.c	(revision 15274)
@@ -1,5 +1,3 @@
 # include "dvoshell.h"
-
-int bracket (double *list, double value, int Nlist, int mode);
 
 int subpix (int argc, char **argv) {
@@ -65,6 +63,6 @@
 
   /* bracket the DEC range of interest */
-  Nlo = bracket (DEC, Dec - Radius, Nstars, FALSE);
-  Nhi = bracket (DEC, Dec + Radius, Nstars, TRUE);
+  Nlo = bracket (DEC, Nstars, FALSE, Dec - Radius);
+  Nhi = bracket (DEC, Nstars, TRUE,  Dec + Radius);
   ALLOCATE (entry, int, MAX (Nhi - Nlo, 1));
   Nentry = 0;
@@ -155,34 +153,2 @@
   return (TRUE);
 }
-
-/* fast operation to find an entry just below (0) or above (1) value */
-int bracket (double *list, double value, int Nlist, int mode) {
-
-  int Nlo, Nhi, N;
-
-  if (mode == 0) {
-    Nlo = 0; Nhi = Nlist;
-    while (Nhi - Nlo > 10) {
-      N = 0.5*(Nlo + Nhi);
-      if (list[N] < value) {
-	Nlo = N;
-      } else {
-	Nhi = N + 1;
-      }
-    }
-    return (Nlo);
-  }
-  if (mode == 1) {
-    Nlo = 0; Nhi = Nlist;
-    while (Nhi - Nlo > 10) {
-      N = 0.5*(Nlo + Nhi);
-      if (list[N] > value) {
-	Nhi = N;
-      } else {
-	Nlo = N - 1;
-      }
-    }
-    return (Nhi);
-  }
-  return (0);
-}
Index: /trunk/Ohana/src/opihi/include/data.h
===================================================================
--- /trunk/Ohana/src/opihi/include/data.h	(revision 15273)
+++ /trunk/Ohana/src/opihi/include/data.h	(revision 15274)
@@ -168,4 +168,7 @@
 void *db_getConnection ();
 
+int bracket (double *list, int Nlist, int mode, double value);
+int ibracket (int *list, int Nlist, int mode, double value);
+
 # endif
 
Index: /trunk/Ohana/src/opihi/lib.data/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/Makefile	(revision 15273)
+++ /trunk/Ohana/src/opihi/lib.data/Makefile	(revision 15274)
@@ -22,4 +22,5 @@
 $(SDIR)/svdcmp.$(ARCH).o		\
 $(SDIR)/convert.$(ARCH).o		\
+$(SDIR)/bracket.$(ARCH).o		\
 $(SDIR)/spline.$(ARCH).o		\
 $(SDIR)/mrqmin.$(ARCH).o		\
Index: /trunk/Ohana/src/opihi/lib.data/bracket.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/bracket.c	(revision 15274)
+++ /trunk/Ohana/src/opihi/lib.data/bracket.c	(revision 15274)
@@ -0,0 +1,65 @@
+# include "data.h"
+
+/* fast operation to find an entry just below (0) or above (1) value */
+/* input list must be sorted */
+int bracket (double *list, int Nlist, int mode, double value) {
+
+  int Nlo, Nhi, N;
+
+  if (mode == 0) {
+    Nlo = 0; Nhi = Nlist;
+    while (Nhi - Nlo > 10) {
+      N = 0.5*(Nlo + Nhi);
+      if (list[N] < value) {
+	Nlo = N;
+      } else {
+	Nhi = N + 1;
+      }
+    }
+    return (Nlo);
+  }
+  if (mode == 1) {
+    Nlo = 0; Nhi = Nlist;
+    while (Nhi - Nlo > 10) {
+      N = 0.5*(Nlo + Nhi);
+      if (list[N] > value) {
+	Nhi = N;
+      } else {
+	Nlo = N - 1;
+      }
+    }
+    return (Nhi);
+  }
+  return (0);
+}
+
+int ibracket (int *list, int Nlist, int mode, double value) {
+
+  int Nlo, Nhi, N;
+
+  if (mode == 0) {
+    Nlo = 0; Nhi = Nlist;
+    while (Nhi - Nlo > 10) {
+      N = 0.5*(Nlo + Nhi);
+      if (list[N] < value) {
+	Nlo = N;
+      } else {
+	Nhi = N + 1;
+      }
+    }
+    return (Nlo);
+  }
+  if (mode == 1) {
+    Nlo = 0; Nhi = Nlist;
+    while (Nhi - Nlo > 10) {
+      N = 0.5*(Nlo + Nhi);
+      if (list[N] > value) {
+	Nhi = N;
+      } else {
+	Nlo = N - 1;
+      }
+    }
+    return (Nhi);
+  }
+  return (0);
+}