Index: trunk/Ohana/src/opihi/cmd.data/extract.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/extract.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.data/extract.c	(revision 41341)
@@ -1,11 +1,25 @@
 # 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 i, j;
-  float *Vin, *Vout;
-  int sx, sy, nx, ny, NX, NY;
-  int Sx, Sy, Nx, Ny;
+  int N;
   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) {
@@ -15,22 +29,24 @@
 
   if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
-  NX = in[0].matrix.Naxis[0];
-  NY = in[0].matrix.Naxis[1];
+  int NX = in[0].matrix.Naxis[0];
+  int NY = in[0].matrix.Naxis[1];
 
-  sx = atof (argv[3]);
-  sy = atof (argv[4]);
-  nx = atof (argv[5]);
-  ny = atof (argv[6]);
+  int sx = atof (argv[3]);
+  int sy = atof (argv[4]);
+  int nx = atof (argv[5]);
+  int ny = atof (argv[6]);
 
-  Sx = atof (argv[7]);
-  Sy = atof (argv[8]);
-  Nx = atof (argv[9]);
-  Ny = atof (argv[10]);
+  int Sx = atof (argv[7]);
+  int Sy = atof (argv[8]);
+  int Nx = atof (argv[9]);
+  int Ny = atof (argv[10]);
 
   if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
-    gprint (GP_ERR, "mismatch between source and dest regions\n");
+    gprint (GP_ERR, "source pixels extend beyond target pixels\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 */
@@ -39,14 +55,10 @@
       (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);
   }
 
@@ -76,13 +88,61 @@
   }
 
-  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++) {
+  // (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++) {
       if (i + sx < 0) continue;
       if (i + sx >= NX) continue;
-      *Vout = *Vin;
+      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];
     }
   }
@@ -91,3 +151,2 @@
 
 }
-