Index: trunk/Ohana/src/opihi/cmd.data/extract.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/extract.c	(revision 41164)
+++ trunk/Ohana/src/opihi/cmd.data/extract.c	(revision 41340)
@@ -1,25 +1,11 @@
 # include "data.h"
-
-/* <from> : source image, must exist
-   <to>   : target image -- if it does not exist, it is created of size (Nx,Ny)
-   sx, sy : source starting coordinate -- need not be on a valid image pixel
-   nx, ny : number of source pixels -- currently must not go out of bounds -> allow out-of-bounds
-   Sx, Sy : target starting coordinate -- need not be on a valid image pixel
-   Nx, Ny : redundant information UNLESS <to> does exist (in which case it is required information) -> make Nx,Ny optional if <to> exists?
- */
 
 int extract (int argc, char **argv) {
   
-  int N;
+  int i, j;
+  float *Vin, *Vout;
+  int sx, sy, nx, ny, NX, NY;
+  int Sx, Sy, Nx, Ny;
   Buffer *in, *out;
-
-  float initValue = 0.0;
-  int initOutput = FALSE;
-  if ((N = get_argument (argc, argv, "-init"))) {
-    remove_argument (N, &argc, argv);
-    initValue = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    initOutput = TRUE;
-  }
 
   if (argc != 11) {
@@ -29,24 +15,22 @@
 
   if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
-  int NX = in[0].matrix.Naxis[0];
-  int NY = in[0].matrix.Naxis[1];
+  NX = in[0].matrix.Naxis[0];
+  NY = in[0].matrix.Naxis[1];
 
-  int sx = atof (argv[3]);
-  int sy = atof (argv[4]);
-  int nx = atof (argv[5]);
-  int ny = atof (argv[6]);
+  sx = atof (argv[3]);
+  sy = atof (argv[4]);
+  nx = atof (argv[5]);
+  ny = atof (argv[6]);
 
-  int Sx = atof (argv[7]);
-  int Sy = atof (argv[8]);
-  int Nx = atof (argv[9]);
-  int Ny = atof (argv[10]);
+  Sx = atof (argv[7]);
+  Sy = atof (argv[8]);
+  Nx = atof (argv[9]);
+  Ny = atof (argv[10]);
 
   if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
-    gprint (GP_ERR, "source pixels extend beyond target pixels\n");
+    gprint (GP_ERR, "mismatch between source and dest regions\n");
     gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
-    // return (FALSE);
+    return (FALSE);
   }
-
-  // XXX : allow source region to fall outside source image
 
   /* region is not on first image */
@@ -55,10 +39,14 @@
       (sy > in[0].matrix.Naxis[1])) {
     gprint (GP_ERR, "region outside of source image\n");
-    // return (FALSE);
+    return (FALSE);
   }
 
+  if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
+    gprint (GP_ERR, "source region larger than dest region\n");
+    return (FALSE);
+  }
   if ((Sx < 0) || (Sy < 0)) {
     gprint (GP_ERR, "dest region out of range\n");
-    // return (FALSE);
+    return (FALSE);
   }
 
@@ -88,61 +76,13 @@
   }
 
-  // (NX,NY) : source image dimensions
-  // (Nx,Ny) : target image dimensions
-
-  // allow (sx, sy), (Sx, Sy) to be off image
-  // allow (sx+nx, sy+ny) to be off image
-
-  // if (sx < 0) {
-  //   nx += sx;
-  //   sx = 0;
-  // }
-  // if (sx > NX) sx = NX;
-  // 
-  // if (Sx < 0) {
-  //   Nx += sx;
-  //   sx = 0;
-  // }
-  // if (sx > NX) sx = NX;
-  
-  float *Vin = (float *)(in[0].matrix.buffer);
-  float *Vout = (float *)(out[0].matrix.buffer);
-
-  if (initOutput) {
-    for (int j = 0; j < Ny; j++) {
-      for (int i = 0; i < Nx; i++) {
-	Vout[i + Nx*j] = initValue;
-      }
-    }
-  }
-
-  int Nps = NX*NY;
-  int Npt = Nx*Ny;
-
-  for (int j = 0; j < ny; j++) {
-    if (j + sy < 0) continue; // not yet on source image
-    if (j + Sy < 0) continue; // not yet on target image
-    if (j + sy >= NY) continue; // past edge of source image
-    if (j + Sy >= Ny) continue; // past edge of target image
-
-    // Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx;  
-    // Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
-
-    for (int i = 0; i < nx; i++) {
+  for (j = 0; j < ny; j++) {
+    if (j + sy < 0) continue;
+    if (j + sy >= NY) continue;
+    Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx;  
+    Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
+    for (i = 0; i < nx; i++, Vin++, Vout++) {
       if (i + sx < 0) continue;
       if (i + sx >= NX) continue;
-      if (i + Sx < 0) continue;
-      if (i + Sx >= Nx) continue;
-
-      int ps = i + sx + (j + sy)*NX;
-      int pt = i + Sx + (j + Sy)*Nx;
-
-      myAssert (pt >= 0, "oops 1");
-      myAssert (pt < Npt, "oops 2");
-
-      myAssert (ps >= 0, "oops 3");
-      myAssert (ps < Nps, "oops 4");
-
-      Vout[pt] = Vin[ps];
+      *Vout = *Vin;
     }
   }
@@ -151,2 +91,3 @@
 
 }
+