Index: trunk/Ohana/src/opihi/cmd.data/cut.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/cut.c	(revision 41164)
+++ trunk/Ohana/src/opihi/cmd.data/cut.c	(revision 41340)
@@ -1,71 +1,10 @@
 # include "data.h"
 
-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;
-}
+enum {SUM, MEAN, MEDIAN};
 
 int cut (int argc, char **argv) {
   
   int i, j, N, Nx, Ny, Mode;
-  float *Vin, *Vbuf;
+  float *Vin, *Vbuf, value;
   int sx, sy, nx, ny;
   Vector *xvec, *yvec;
@@ -73,8 +12,4 @@
 
   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);
@@ -84,12 +19,4 @@
     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;
   }
 
@@ -109,7 +36,6 @@
   Ny = buf[0].matrix.Naxis[1];
 
-  // 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);
+  if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
+    gprint (GP_ERR, "region out of range\n");
     return (FALSE);
   }
@@ -131,31 +57,21 @@
 
     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;
+      /* accumulate values */
+      Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i; 
+      for (j = 0; j < ny; j++, Vin += Nx) {
+	Vbuf[j] = *Vin;
       }
-      if (i + sx >= Nx) {
-	yvec[0].elements.Flt[i] = NAN;
-	continue;
+      /* 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; }
       }
-
-      /* 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) {
-	if (!isfinite(*Vin)) continue;
-	Vbuf[Npix] = *Vin;
-	Npix ++;
-      }
-      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
+      yvec[0].elements.Flt[i] = value;
     }
     free (Vbuf);
@@ -173,31 +89,21 @@
 
     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;
+      /* accumulate values */
+      Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx; 
+      for (j = 0; j < nx; j++, Vin ++) {
+	Vbuf[j] = *Vin;
       }
-      if (i + sy >= Ny) {
-	yvec[0].elements.Flt[i] = NAN;
-	continue;
+      /* 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; }
       }
-
-      /* 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 ++) {
-	if (!isfinite(*Vin)) continue;
-	Vbuf[Npix] = *Vin;
-	Npix ++;
-      }
-      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
+      yvec[0].elements.Flt[i] = value;
     }
     free (Vbuf);