Index: trunk/Ohana/src/opihi/cmd.data/cut.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/cut.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.data/cut.c	(revision 41341)
@@ -1,10 +1,71 @@
 # include "data.h"
 
-enum {SUM, MEAN, MEDIAN};
+enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
+
+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 == NGOOD) {
+    return Nvalue;
+  }
+  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;
@@ -12,4 +73,8 @@
 
   Mode = SUM;
+  if ((N = get_argument (argc, argv, "-sum"))) {
+    remove_argument (N, &argc, argv);
+    Mode = SUM;
+  }
   if ((N = get_argument (argc, argv, "-median"))) {
     remove_argument (N, &argc, argv);
@@ -19,4 +84,12 @@
     remove_argument (N, &argc, argv);
     Mode = MEAN;
+  }
+  if ((N = get_argument (argc, argv, "-inner"))) {
+    remove_argument (N, &argc, argv);
+    Mode = INNER;
+  }
+  if ((N = get_argument (argc, argv, "-ngood"))) {
+    remove_argument (N, &argc, argv);
+    Mode = NGOOD;
   }
 
@@ -36,6 +109,7 @@
   Ny = buf[0].matrix.Naxis[1];
 
-  if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
-    gprint (GP_ERR, "region out of range\n");
+  // ny & nx do not need to be constrained by the buffer, but they do need to be sensible
+  if ((nx < 0) || (ny < 0) || (nx > 1e8) || (ny > 1e8)) {
+    gprint (GP_ERR, "warning : extraction size is crazy: %d,%d\n", nx, ny);
     return (FALSE);
   }
@@ -57,21 +131,31 @@
 
     for (i = 0; i < nx; i++) {
+      // for out-of-range areas, set the output pixels to NAN
+      if (i + sx < 0) {
+	yvec[0].elements.Flt[i] = NAN;
+	continue;
+      }
+      if (i + sx >= Nx) {
+	yvec[0].elements.Flt[i] = NAN;
+	continue;
+      }
+
       /* accumulate values */
+      // skip out-of-range areas in y
+      if (sy < 0) {
+	ny += sy;
+	sy = 0;
+      }
+      if (sy > Ny) sy = Ny;
+      if (sy + ny >= Ny) ny = Ny - sy;
+
       Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i; 
+      int Npix = 0;
       for (j = 0; j < ny; j++, Vin += Nx) {
-	Vbuf[j] = *Vin;
-      }
-      /* apply stat of choice */
-      if (Mode == MEDIAN) {
-	fsort (Vbuf, ny);
-	value = Vbuf[(int)(0.5*ny)];
-      } else {
-	value = 0;
-	for (j = 0; j < ny; j++) {
-	  value += Vbuf[j];
-	}
-	if (Mode == MEAN) { value /= ny; }
-      }
-      yvec[0].elements.Flt[i] = value;
+	if (!isfinite(*Vin)) continue;
+	Vbuf[Npix] = *Vin;
+	Npix ++;
+      }
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     }
     free (Vbuf);
@@ -89,21 +173,31 @@
 
     for (i = 0; i < ny; i++) {
+      // for out-of-range areas, set the output pixels to NAN
+      if (i + sy < 0) {
+	yvec[0].elements.Flt[i] = NAN;
+	continue;
+      }
+      if (i + sy >= Ny) {
+	yvec[0].elements.Flt[i] = NAN;
+	continue;
+      }
+
       /* accumulate values */
+      // skip out-of-range areas in x
+      if (sx < 0) {
+	nx += sx;
+	sx = 0;
+      }
+      if (sx > Nx) sx = Nx;
+      if (sx + nx >= Nx) nx = Nx - sx;
+
       Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx; 
+      int Npix = 0;
       for (j = 0; j < nx; j++, Vin ++) {
-	Vbuf[j] = *Vin;
-      }
-      /* apply stat of choice */
-      if (Mode == MEDIAN) {
-	fsort (Vbuf, nx);
-	value = Vbuf[(int)(0.5*nx)];
-      } else {
-	value = 0;
-	for (j = 0; j < nx; j++) {
-	  value += Vbuf[j];
-	}
-	if (Mode == MEAN) { value /= nx; }
-      }
-      yvec[0].elements.Flt[i] = value;
+	if (!isfinite(*Vin)) continue;
+	Vbuf[Npix] = *Vin;
+	Npix ++;
+      }
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
     }
     free (Vbuf);