Changeset 41916 for branches/eam_branches/ipp-20211108/Ohana

Timestamp:

Nov 16, 2021, 8:28:36 AM (5 years ago)

Author:

eugene

Message:

clean up imfit, fix imfit-trail, add imfit tests

Location:

branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro

Files:

• imfit-fgauss.c (modified) (1 diff)
• imfit-qrgauss.c (modified) (1 diff)
• imfit-trail.c (modified) (9 diffs)
• imfit.c (modified) (9 diffs)
• test/imfit.sh (added)

branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-fgauss.c

@@ -1 +1 @@
 # include "imfit.h"
+
+/** fgauss : a real 2D Gaussian **/
 
 opihi_flt fgaussTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *);

branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-qrgauss.c

@@ -59 +59 @@
 
   r = 1.0 / (1 + fpar[0]*z + pow(z,par[7]));
-  f = par[5]*r + par[6];
+  f = par[5]*r + par[6]; // Io * f(r) + Sky
   q = par[5]*SQ(r)*(fpar[0] + par[7]*pow(z,(par[7]-1)));
 

branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit-trail.c

@@ -4 +4 @@
 void  trailCL ();
 
+// fitted parameters:
 # define PAR_X      0
 # define PAR_Y      1
 # define PAR_THETA  2
-# define PAR_SIGMA  3
-# define PAR_LENGTH 4
-# define PAR_I0     5
-# define PAR_SKY    6
+# define PAR_LENGTH 3
+# define PAR_I0     4
+# define PAR_SKY    5
+
+// fixed parameters:
+# define PAR_SIGMA  0
 
 void trail_setup (char *name) {
@@ -18 +21 @@
   fitfunc = trailTD;
   imfit_cleanup = trailCL;
-  Npar = 7;
-  Nfpar = 0;
+  Npar  = 6;
+  Nfpar = 1;
 
   /* allocate free and fixed parameters */
@@ -30 +33 @@
   par[PAR_Y      ] = get_variable_default ("Yg",      0.0);
   par[PAR_THETA  ] = get_variable_default ("Tg",      0.0);
-  par[PAR_SIGMA  ] = get_variable_default ("Wg",      2.0);
   par[PAR_LENGTH ] = get_variable_default ("Lg",     10.0);
   par[PAR_I0     ] = get_variable_default ("Zpk", 10000.0);
   par[PAR_SKY    ] = get_variable_default ("Sg",      0.0);
   sky = &par[PAR_SKY];
+
+  fpar[PAR_SIGMA ] = get_variable_default ("Wg",      2.0);
 }
 
@@ -40 +44 @@
   set_variable ("Xg",  par[PAR_X     ]);
   set_variable ("Yg",  par[PAR_Y     ]);
-  set_variable ("Wg",  par[PAR_SIGMA ]);
   set_variable ("Tg",  par[PAR_THETA ]);
   set_variable ("Lg",  par[PAR_LENGTH]);
   set_variable ("Zpk", par[PAR_I0    ]);
   set_variable ("Sg",  par[PAR_SKY   ]);
+
+  set_variable ("Wg",  fpar[PAR_SIGMA]);
 }
 
@@ -54 +59 @@
   opihi_flt Y = y - par[PAR_Y];
   
-  opihi_flt S2 = 2.0 * SQ(par[PAR_SIGMA]);
+  opihi_flt S2 = 2.0 * SQ(fpar[PAR_SIGMA]);
 
   opihi_flt ST = sin(RAD_DEG*par[PAR_THETA]);
@@ -78 +83 @@
 
     // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
+    // since X = X - par[PAR_X], dFoo/dXo = -dFoo/dX
     float dRdX = +ST;
     float dRdY = -CT;
-    float dRdT = -Y*ST - X*CT;
+    float dRdT = (-Y*ST - X*CT)*RAD_DEG;
+    // note PAR_THETA is in degrees
 
     float dGdX = dGdR * dRdX;
     float dGdY = dGdR * dRdY;
     float dGdT = dGdR * dRdT;
+    // dGdL is 0.0 because dRdL is 0.0 (R is not a function of L)
 
     // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
@@ -96 +104 @@
     float dZmdL = -0.5 / sqrt(S2);
 
-    float dZpdT = (-X*ST + Y*CT) / sqrt(S2);
+    // note PAR_THETA is in degrees
+    float dZpdT = (-X*ST + Y*CT) * RAD_DEG / sqrt(S2);
     float dZmdT = dZpdT; // dZpdT = dZmdT
 
@@ -119 +128 @@
     // dGdL is 0.0 because dRdL is 0.0
     float dPdL = Gxy * (dEpdL - dEmdL);
-
     float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT);
 
@@ -127 +135 @@
     dpar[PAR_LENGTH] = par[PAR_I0] * dPdL;
     dpar[PAR_THETA]  = par[PAR_I0] * dPdT;
-    dpar[PAR_SIGMA]  = 0;       // we don't actually allow this to vary, so we do not need to calculate it
   }
 

branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit.c

@@ -3 +3 @@
 int imfit (int argc, char **argv) {
 
-  int i, j, N, Npts, Save, VERBOSE;
-  int sx, sy, nx, ny, Nx, Ny;
-  float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold;
-  opihi_flt *x, *y, *z, *dz;
-  float *V;
+  int N;
   Buffer *buf;
 
-  Save = FALSE;
+  char *Save = NULL;
   if ((N = get_argument (argc, argv, "-save"))) {
     remove_argument (N, &argc, argv);
-    Save = TRUE;
-  }
-
-  // int ShapeVariation = FALSE;
-  // if ((N = get_argument (argc, argv, "-shapes"))) {
-  //   remove_argument (N, &argc, argv);
-  //   ShapeVariation = TRUE;
-  // }
-
-  SatThreshold = 0xffff;
+    Save = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int Insert = FALSE;
+  if ((N = get_argument (argc, argv, "-insert"))) {
+    remove_argument (N, &argc, argv);
+    Insert = TRUE;
+    if (Save) { gprint (GP_ERR, "-save and -insert are mutually exclusive\n"); free (Save); return (FALSE); }
+  }
+
+  int SatThreshold = 0xffff;
   if ((N = get_argument (argc, argv, "-sat"))) {
     remove_argument (N, &argc, argv);
@@ -30 +28 @@
 
   /* Gain in e/DN */
-  Gain = 1.0;
+  float Gain = 1.0;
   if ((N = get_argument (argc, argv, "-gain"))) {
     remove_argument (N, &argc, argv);
@@ -38 +36 @@
 
   /* RD noise in DN */
-  RDnoise = 0.0;
+  float RDnoise = 0.0;
   if ((N = get_argument (argc, argv, "-rdnoise"))) {
     remove_argument (N, &argc, argv);
@@ -45 +43 @@
   }
 
-  VERBOSE = FALSE;
+  int VERBOSE = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
     remove_argument (N, &argc, argv);
@@ -51 +49 @@
   }
 
-  /* set fitting function */
+  /* set fitting function : defines par, Npar, fitfunc, etc globals (imfit.h) */
   fgauss_setup ("fgauss");
   if ((N = get_argument (argc, argv, "-func"))) {
     fitfunc = NULL;
     remove_argument (N, &argc, argv);
-    fgauss_setup (argv[N]);
-    pgauss_setup (argv[N]);
-    pgauss_psf_setup (argv[N]);
+    fgauss_setup (argv[N]); // OK
+    pgauss_setup (argv[N]); // OK 
     sgauss_setup (argv[N]);
-    sgauss_psf_setup (argv[N]);
-    qgauss_setup (argv[N]);
-    qgauss_psf_setup (argv[N]);
-    qfgauss_setup (argv[N]);
+    qgauss_setup (argv[N]); // OK
+    qfgauss_setup (argv[N]); // OK
     qrgauss_setup (argv[N]);
     trail_setup (argv[N]);
+    pgauss_psf_setup (argv[N]);
+    sgauss_psf_setup (argv[N]);
+    qgauss_psf_setup (argv[N]);
     if (fitfunc == NULL) {
       gprint (GP_ERR, "unknown function %s\n", argv[N]);
+      FREE (Save);
       return (FALSE);
     }
@@ -74 +73 @@
 
   if (argc != 6) {
-    gprint (GP_ERR, "USAGE: imfit <buffer> sx sy nx ny\n");
+    gprint (GP_ERR, "USAGE: imfit <buffer> Xo Yo dX dY\n");
+    gprint (GP_ERR, "options: [-save buffer] [-insert] [-sat value] [-gain value] [-rdnoise value] [-v] [-func option]\n");
+    gprint (GP_ERR, "   (Xo,Yo) : center\n");
+    gprint (GP_ERR, "   (dX,dY) : window size\n");
+    FREE (Save);
     return (FALSE);
   }
 
   /* non-optional arguments */
-  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
-  sx = atof (argv[2]);
-  sy = atof (argv[3]);
-  nx = atof (argv[4]);
-  ny = atof (argv[5]);
-  Nx = buf[0].matrix.Naxis[0];
-  Ny = buf[0].matrix.Naxis[1];
-
-  /* check if region is valid */
-  if (sx + 0.5*nx < 0) goto range;
-  if (sy + 0.5*ny < 0) goto range;
-  if (sx + 0.5*nx >= Nx) goto range;
-  if (sy + 0.5*ny >= Ny) goto range;
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) { FREE (Save); return (FALSE); }
+  int Xo = atof (argv[2]);
+  int Yo = atof (argv[3]);
+  int dX = atof (argv[4]);
+  int dY = atof (argv[5]);
+  int Nx = buf[0].matrix.Naxis[0];
+  int Ny = buf[0].matrix.Naxis[1];
+
+  int sx = Xo - dX/2;
+  int sy = Yo - dY/2;
+
+  /* check if region is valid (center must be in range of image pixels) */
+  if (Xo < 0) goto range;
+  if (Yo < 0) goto range;
+  if (Xo >= Nx) goto range;
+  if (Yo >= Ny) goto range;
+
+  // image value in DN
+  // rdnoise in DN
+  // sigma_DN^2 = sigma_e^2 / gain^2
+  // sigma_e^2  = Ne
+  // sigma_e^2  = DN * gain
+  // sigma_DN^2 = DN * gain / gain^2 = DN / gain
+
+  if (Insert) {
+    float *Vi = (float *)buf[0].matrix.buffer;
+    for (int j = 0; j < dY; j++) {
+      for (int i = 0; i < dX; i++) {
+        float vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
+        Vi[(i+sx)+(j+sy)*Nx] += vf;
+      }
+    }
+    return TRUE;
+  }
 
   /* convert array z[x,y] to x[i], y[i], z[i] */
   N = 0;
-  Npts = nx*ny;
-  ALLOCATE (x,  opihi_flt, 2*Npts);
-  ALLOCATE (y,  opihi_flt, 2*Npts);
-  ALLOCATE (z,  opihi_flt, 2*Npts);
-  ALLOCATE (dz, opihi_flt, 2*Npts);
-  for (j = 0; j < ny; j++) {
+  int Npts = dX*dY;
+  ALLOCATE_PTR (x,  opihi_flt, 2*Npts);
+  ALLOCATE_PTR (y,  opihi_flt, 2*Npts);
+  ALLOCATE_PTR (z,  opihi_flt, 2*Npts);
+  ALLOCATE_PTR (dz, opihi_flt, 2*Npts);
+  for (int j = 0; j < dY; j++) {
     if (j + sy < 0) continue;
     if (j + sy >= Ny) continue;
-    V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; 
-    for (i = 0; i < nx; i++) {
+    float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; 
+    for (int i = 0; i < dX; i++) {
       if (i + sx < 0) continue;
       if (i + sx >= Nx) continue;
       if (*V > SatThreshold) goto next;
-      dz[N] = (SQ(RDnoise) + *V/Gain);
+      dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
       if (dz[N] <= 0) goto next;
       dz[N] = 1.0 / dz[N];
@@ -122 +146 @@
 
   /* run fit routine */
-  ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
-  dchisq = ochisq;
-  chisq  = ochisq;
-
-//for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
-
-  for (i = 0; (i < 25); i++) {
+  float ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
+  float dchisq = ochisq;
+  float chisq  = ochisq;
+
+  int Niter = 0;
+  // for (int i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
+
+  // keep iterating if chisq is increasing or 
+  for (Niter = 0; (Niter < 25) && ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar))); Niter++) {
     chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     dchisq = ochisq - chisq;
+    // fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
     ochisq = chisq;
-    fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
   }  
-  set_int_variable ("Niter",  i);
+  set_int_variable ("Niter",  Niter);
 
   /** create output image (keep in sky) **/
@@ -141 +167 @@
     float *Vi, *Vo, vr, vf;
 
-    if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
+    if ((out = SelectBuffer (Save, ANYBUFFER, TRUE)) == NULL) { free (Save); return (FALSE); }
     free (out[0].header.buffer);
     free (out[0].matrix.buffer);
 
     strcpy (out[0].file, "(empty)");
-    if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE;
-
-    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
+    if (!CreateBuffer (out, 2*dX, 2*dY, -32, 0.0, 1.0)) { free (Save); return FALSE; }
+
+    /* four panels: 1) raw image. 2) fit  3) raw - fit  4) absolute deviation */
     Vi = (float *)buf[0].matrix.buffer;
     Vo = (float *)out[0].matrix.buffer;
-    for (j = 0; j < ny; j++) {
-      for (i = 0; i < nx; i++) {
+    for (int j = 0; j < dY; j++) {
+      for (int i = 0; i < dX; i++) {
         vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
         vr = Vi[(i+sx)+(j+sy)*Nx];
-        Vo[(i   )+(j   )*2*nx] = vr;
-        Vo[(i+nx)+(j   )*2*nx] = vf;
-        Vo[(i   )+(j+ny)*2*nx] = vr - vf + *sky;
-        Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky;
+        Vo[(i   )+(j   )*2*dX] = vr;
+        Vo[(i+dX)+(j   )*2*dX] = vf;
+        Vo[(i   )+(j+dY)*2*dX] = vr - vf + *sky;
+        Vo[(i+dX)+(j+dY)*2*dX] = fabs(vr-vf) + *sky;
       }
     }
+    free (Save);
   }
 
@@ -169 +196 @@
 
   if (VERBOSE) {
-    for (i = 0; i < Npar; i++) {
+    for (int i = 0; i < Npar; i++) {
       gprint (GP_ERR, "%g ", par[i]);
     }