Changeset 40830

Timestamp:

Jul 9, 2019, 4:40:10 PM (7 years ago)

Author:

eugene

Message:

allow the extraction and the cut to extend beyond the valid image pixels (and fill with nan values)

Location:

branches/eam_branches/ohana.20190329/src/opihi

Files:

• cmd.data/cut.c (modified) (7 diffs)
• cmd.data/extract.c (modified) (5 diffs)
• mana/deimos_mkslit.c (modified) (2 diffs)

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/cut.c

@@ -1 +1 @@
 # include "data.h"
 
-enum {SUM, MEAN, MEDIAN, INNER};
+enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
 
 double cutstat (float *value, int Nvalue, int mode) {
@@ -50 +50 @@
     return output;
   }
+  if (mode == NGOOD) {
+    return Nvalue;
+  }
   if ((mode == MEAN) || (mode == SUM)) {
     for (int i = 0; i < Nvalue; i++) {
@@ -70 +73 @@
 
   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);
@@ -81 +88 @@
     remove_argument (N, &argc, argv);
     Mode = INNER;
+  }
+  if ((N = get_argument (argc, argv, "-ngood"))) {
+    remove_argument (N, &argc, argv);
+    Mode = NGOOD;
   }
 
@@ -98 +109 @@
   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);
   }
@@ -119 +131 @@
 
     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) {
         if (!isfinite(*Vin)) continue;
-        Vbuf[j] = *Vin;
+        Vbuf[Npix] = *Vin;
         Npix ++;
       }
@@ -143 +173 @@
 
     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 ++) {
         if (!isfinite(*Vin)) continue;
-        Vbuf[j] = *Vin;
+        Vbuf[Npix] = *Vin;
         Npix ++;
       }

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/extract.c

@@ -1 +1 @@
 # 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 +29 @@
 
   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 +55 @@
       (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 +88 @@
   }
 
-  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 +151 @@
 
 }
-

branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkslit.c

@@ -47 +47 @@
   }
 
+  // either flux or fluxbins must be specified to define output size
+  int profilebins = 0;
+  if ((N = get_argument (argc, argv, "-profilebins"))) {
+    remove_argument (N, &argc, argv);
+    profilebins = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   // for a red vs blu spline, we need to specify the split point
   // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
@@ -74 +82 @@
 
   // define the output window
-  int Nx = 2*profile[0].Nelements + 1;
+  // if no profile half-width is specified, use supplied profile 
+  int Nx = profilebins ? 2*profilebins + 1 : 2*profile[0].Nelements + 1;
   int Ny = flux ? flux[0].Nelements : fluxbins;