Changeset 41516 for trunk/Ohana/src/opihi/cmd.data/densify.c

Timestamp:

Apr 2, 2021, 1:50:33 PM (5 years ago)

Author:

eugene

Message:

add -raw, -value, -binning options to densify (as in cdensify) and fix usage message; init Rmin for match1d; init chisq for vgauss

File:

• trunk/Ohana/src/opihi/cmd.data/densify.c (modified) (11 diffs)

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

@@ -13 +13 @@
   opihi_flt *x, *y;
 
-  // int Normalize = TRUE;
-  // if ((N = get_argument (argc, argv, "-raw"))) {
-  //   remove_argument (N, &argc, argv);
-  //   Normalize = FALSE;
-  // }
+  int Normalize = TRUE;
+  if ((N = get_argument (argc, argv, "-raw"))) {
+    remove_argument (N, &argc, argv);
+    Normalize = FALSE;
+  }
 
   UseGraph = FALSE;
@@ -25 +25 @@
   }
 
+  Vector *vv = NULL;
+  if ((N = get_argument (argc, argv, "-value"))) {
+    remove_argument (N, &argc, argv);
+    if ((vv = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); 
+    remove_argument (N, &argc, argv);
+  }
+
   float scale = 0.0;
   if ((N = get_argument (argc, argv, "-scale"))) {
     remove_argument (N, &argc, argv);
     scale = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int binning = 1;
+  if ((N = get_argument (argc, argv, "-binning"))) {
+    if (!UseGraph) {
+      gprint (GP_ERR, "-binning only valid for -graph option\n");
+      return FALSE;
+    }
+    remove_argument (N, &argc, argv);
+    binning = atoi(argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -47 +65 @@
     gprint (GP_ERR, "   OR: densify buffer x y -graph\n");
     gprint (GP_ERR, " option: -psf [dot] (circle) (square) (gauss)\n");
+    gprint (GP_ERR, " other options:\n");
+    gprint (GP_ERR, "   -raw : do not renormalize PSF\n");
+    gprint (GP_ERR, "   -scale : spatial scale factor for PSF (radius, half-width, or sigma)\n");
+    gprint (GP_ERR, "   -value (vector) : multiply sum by the given vector\n");
     return (FALSE);
   }
@@ -55 +77 @@
 
   if (vx[0].Nelements != vy[0].Nelements) return (FALSE);
+
+  if (vv) {
+    if (vv[0].Nelements != vx[0].Nelements) {
+      gprint (GP_ERR, "mis-match in vector lengths\n");
+      return FALSE;
+    }
+  }
 
   REQUIRE_VECTOR_FLT (vx, FALSE); 
@@ -68 +97 @@
     Ymax = graphmode.ymax;
     Ymin = graphmode.ymin;
-    dX = (Xmax - Xmin) / (Xpix - 1);
-    dY = (Ymax - Ymin) / (Ypix - 1);
+    dX = binning * (Xmax - Xmin) / (Xpix - 1);
+    dY = binning * (Ymax - Ymin) / (Ypix - 1);
   } else {
     Xmin = atof (argv[4]);
@@ -91 +120 @@
   float scaleY = (scale > 0.0) ? scale / dY : 3.0;
 
-  Nx = (Xmax - Xmin) / dX + 1;
-  Ny = (Ymax - Ymin) / dY + 1;
+  Nx = abs((Xmax - Xmin) / dX) + 1;
+  Ny = abs((Ymax - Ymin) / dY) + 1;
   
   gfits_free_matrix (&bf[0].matrix);
@@ -107 +136 @@
   x = vx[0].elements.Flt;
   y = vy[0].elements.Flt;
+
+  opihi_flt *Fs = vv ? vv[0].elements.Flt : NULL;
+  opihi_int *Is = vv ? vv[0].elements.Int : NULL;
+  int isFloatScale = (vv && vv[0].type == OPIHI_FLT);
+
   val = (float *)bf[0].matrix.buffer;
   for (i = 0; i < vx[0].Nelements; i++, x++, y++) {
     Xb = (*x - Xmin) / dX;
     Yb = (*y - Ymin) / dY;
+
+    float F = 1.0;
+    if (vv) { F = isFloatScale ? Fs[i] : Is[i]; }
+
     switch (PSFTYPE) {
       case IS_DOT:
@@ -117 +155 @@
         if (Xb < 0) continue;
         if (Yb < 0) continue;
-        val[Xb + Yb*Nx] ++;
+        if (vv) {
+          val[Xb + Yb*Nx] += F;
+        } else {
+          val[Xb + Yb*Nx] ++;
+        }
         break;
       case IS_SQUARE:
@@ -126 +168 @@
             if (ix < 0) continue;
             if (iy < 0) continue;
-            val[ix + iy*Nx] += fSquare;
+            if (vv) {
+              val[ix + iy*Nx] += Normalize ? fSquare*F : F;
+            } else {
+              val[ix + iy*Nx] += Normalize ? fSquare : 1.0;
+            }
           }
         }
@@ -141 +187 @@
             if (ix < 0) continue;
             if (iy < 0) continue;
-            val[ix + iy*Nx] += fCircle;
+            if (vv) {
+              val[Xb + Yb*Nx] += Normalize ? fCircle*F : F;
+            } else {
+              val[Xb + Yb*Nx] += Normalize ? fCircle : 1.0;
+            }
           }
         }
@@ -155 +205 @@
             if (ix < 0) continue;
             if (iy < 0) continue;
-            val[ix + iy*Nx] += fGauss*exp(-fSigma*r2);
+            if (vv) {
+              val[Xb + Yb*Nx] += F*fGauss*exp(-fSigma*r2);
+            } else {
+              val[Xb + Yb*Nx] += fGauss*exp(-fSigma*r2);
+            }
           }
         }