Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/cut.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/cut.c	(revision 40811)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/cut.c	(revision 40812)
@@ -1,10 +1,68 @@
 # include "data.h"
 
-enum {SUM, MEAN, MEDIAN};
+enum {SUM, MEAN, MEDIAN, INNER};
+
+double cutstat (float *value, int Nvalue, int mode) {
+
+  double output = 0.0;
+
+  if (Nvalue == 0) return NAN;
+
+  if (mode == MEDIAN) {
+    fsort (value, Nvalue);
+    if (Nvalue % 2) {
+      int Ncenter = Nvalue / 2;
+      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue), "cutstat");
+      output = value[Ncenter];
+    } else {
+      int Ncenter = Nvalue / 2 - 1;
+      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue - 1), "cutstat");
+      output = 0.5*(value[Ncenter] + value[Ncenter + 1]);
+    }
+    return output;
+  } 
+  if (mode == INNER) {
+    fsort (value, Nvalue);
+
+    int Ns = 0, Ne = 0;
+    if (Nvalue % 2) {
+      // for an odd number of points take the same number below
+      // and above the center value
+      int Ncenter = Nvalue / 2;
+      int Nquarter = 0.25*Nvalue;
+      Ns = Ncenter - Nquarter;
+      Ne = Ncenter + Nquarter;
+    } else {
+      // for an even number, the middle lies between two points
+      // take the same number below as above
+      int Ncenter = (int)(Nvalue / 2) - 1;
+      int Nquarter = 0.25*Nvalue;
+      Ns = Ncenter     - Nquarter;
+      Ne = Ncenter + 1 + Nquarter;
+    }
+    int Nv = 0;
+    for (int i = Ns; i <= Ne; i++) {
+      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
+      output += value[i];
+      Nv ++;
+    }
+    output /= Nv;
+    return output;
+  }
+  if ((mode == MEAN) || (mode == SUM)) {
+    for (int i = 0; i < Nvalue; i++) {
+      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
+      output += value[i];
+    }
+    if (mode == MEAN) { output /= Nvalue; }
+    return output;
+  }
+  return NAN;
+}
 
 int cut (int argc, char **argv) {
   
   int i, j, N, Nx, Ny, Mode;
-  float *Vin, *Vbuf, value;
+  float *Vin, *Vbuf;
   int sx, sy, nx, ny;
   Vector *xvec, *yvec;
@@ -19,4 +77,8 @@
     remove_argument (N, &argc, argv);
     Mode = MEAN;
+  }
+  if ((N = get_argument (argc, argv, "-inner"))) {
+    remove_argument (N, &argc, argv);
+    Mode = INNER;
   }
 
@@ -65,21 +127,5 @@
 	Npix ++;
       }
-      if (Npix == 0) {
-	value = NAN;
-	goto set_yvalue;
-      }
-      /* apply stat of choice */
-      if (Mode == MEDIAN) {
-	fsort (Vbuf, Npix);
-	value = Vbuf[(int)(0.5*Npix)];
-      } else {
-	value = 0;
-	for (j = 0; j < Npix; j++) {
-	  value += Vbuf[j];
-	}
-	if (Mode == MEAN) { value /= Npix; }
-      }
-    set_yvalue:
-      yvec[0].elements.Flt[i] = value;
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     }
     free (Vbuf);
@@ -105,21 +151,5 @@
 	Npix ++;
       }
-      if (Npix == 0) {
-	value = NAN;
-	goto set_xvalue;
-      }
-      /* apply stat of choice */
-      if (Mode == MEDIAN) {
-	fsort (Vbuf, Npix);
-	value = Vbuf[(int)(0.5*Npix)];
-      } else {
-	value = 0;
-	for (j = 0; j < Npix; j++) {
-	  value += Vbuf[j];
-	}
-	if (Mode == MEAN) { value /= Npix; }
-      }
-    set_xvalue:
-      yvec[0].elements.Flt[i] = value;
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     }
     free (Vbuf);