Changeset 41340 for trunk/Ohana/src/opihi/cmd.data/cut.c

Timestamp:

Apr 16, 2020, 1:54:47 PM (6 years ago)

Author:

tdeboer

Message:

revert to working Ohana build

File:

• trunk/Ohana/src/opihi/cmd.data/cut.c (modified) (6 diffs)

trunk/Ohana/src/opihi/cmd.data/cut.c

@@ -1 +1 @@
 # 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 +12 @@
 
   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 +19 @@
     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 +36 @@
   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 +57 @@
 
     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 +89 @@
 
     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);