Changeset 35560 for trunk/psModules/src/objects/models/pmModel_TRAIL.c

Timestamp:

May 9, 2013, 12:19:21 PM (13 years ago)

Author:

eugene

Message:

do not crash psphot if variance is not supplied; include mask and mark in modelGuess functions (needed for trail angle guess); fix Reff sx,sy,sxx relationships; define new function to measure guess at the trail angle; tell pmSourceIO if MATCHED_REFS have been read (& skip); do not read sources for WCS type

File:

• trunk/psModules/src/objects/models/pmModel_TRAIL.c (modified) (4 diffs)

trunk/psModules/src/objects/models/pmModel_TRAIL.c

@@ -168 +168 @@
         dPAR[PM_PAR_THETA]  = PAR[PM_PAR_I0] * dPdT;
         dPAR[PM_PAR_SIGMA]  = 0;        // we don't actually allow this to vary, so we do not need to calculate it
+
+        for (int i = 0; i < 7; i++) {
+          if (isnan(dPAR[i])) {
+            fprintf (stderr, "*");
+          }
+        }
+    }
+    if (isnan(f)) {
+      fprintf (stderr, "!");
     }
     return(f);
@@ -228 +237 @@
 }
 
+# define NA 21
+# define NR 21
+static float flux[NA][NR];
+static float npix[NA][NR];
+
+bool pmTrailGetAngle (float *To, float *Io, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal, float sigma) {
+
+  float Xo, Yo;
+  if (!pmModelSetPosition(&Xo, &Yo, source)) return false;
+
+  psImage *image = source->pixels;
+  psImage *mask = source->maskObj;
+  psF32 **imData = image->data.F32;
+  psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
+
+  // do a loop over the pixels, generating (dX,dY) dot (cos(theta),sin(theta))
+  int dP;
+  int dX = Xo - mask->col0;
+  dP = mask->numCols - dX;
+  int DX = PS_MAX(dP, dX);
+  int NX = mask->numCols;
+
+  int dY = Yo - mask->row0;
+  dP = mask->numRows - dY;
+  int DY = PS_MAX(dP, dY);
+  int NY = mask->numRows;
+
+  // just hard wire this for now...
+  float radius = 10.0;
+  float radius2 = PS_SQR(radius);
+
+  // we have an array of Angles x Radii
+  float dT = M_PI / NA;
+  for (int na = 0; na < NA; na++) {
+    memset (flux[na], 0, NR*sizeof(float));
+    memset (npix[na], 0, NR*sizeof(float));
+  }
+
+  // we skip any pixels [real or virtual] outside of the specified radius (nominally the aperture radius)
+  // ix and iy track pixels relative to the centroid
+  for (int ix = -DX; ix < DX + 1; ix++) {
+    if (ix > radius) continue;
+    int mx = ix + dX;
+    for (int iy = -DY; iy < DY + 1; iy++) {
+      if (iy > radius) continue;
+      if (ix*ix + iy*iy > radius2) continue;
+      int my = iy + dY;
+      
+      // include count only the unmasked pixels within the image area
+      if (mx < 0) continue;
+      if (my < 0) continue;
+      if (mx >= NX) continue;
+      if (my >= NY) continue;
+      
+      // count pixels which are masked only with bad pixels
+      if (mkData[my][mx] & maskVal)continue;
+
+      // we have defined NA to be 21
+      int na = 0;
+      for (float angle = 0.0; na < NA; angle += dT, na ++) {
+
+        // XXX optimization : pre-compute the angle sines and cosines
+        float Rad = (ix * cos(angle)) + (iy * sin(angle));
+        int nr = PS_MAX (PS_MIN (NR, Rad + 0.5*NR), 0);
+
+        flux[na][nr] += imData[my][mx];
+        npix[na][nr] ++;
+      }
+    }
+  }
+
+  // generate a psf model (with integral of 1.0)
+  float Ao = 1.0 / sqrt(2*M_PI) / sigma;
+  float psf[NR];
+  for (int nr = 0; nr < NR; nr++) {
+    psf[nr] = Ao*exp(-0.5*PS_SQR((nr - 0.5*NR)/sigma));
+  }
+
+  float fangle[NA];
+  for (int na = 0; na < NA; na++) {
+    float fpsf = 0.0;
+    for (int nr = 0; nr < NR; nr++) {
+      fpsf += psf[nr]*flux[na][nr];
+    }
+    fangle[na] = fpsf;
+    // fprintf (stderr, "fpsf: %f, theta = %f\n", fpsf, PS_DEG_RAD*dT*na);
+  }
+
+  float peak = fangle[0];
+  int pbin = 0;
+  for (int na = 0; na < NA; na++) {
+    if (fangle[na] > peak) {
+      peak = fangle[na];
+      pbin = na;
+    }
+  }
+
+  // result is peak, pbin -> angle
+  *To = dT * pbin;
+  *Io = peak * Ao;
+  
+  return true;
+}
+
 // make an initial guess for parameters
 // 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
-bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
+bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
     psF32 *PAR  = model->params->data.F32;
@@ -260 +373 @@
     //else { size = psfAxes.major; }
 
+    float theta, peak;
+    pmTrailGetAngle (&theta, &peak, source, maskVal, markVal, psfAxes.major);
+
     // axes.major is a sigma in the major direction; scale to 
     PAR[PM_PAR_LENGTH] = 1.5*2.35*size; // a tophat of length L has L = 1.5 * 2.35 * sigma
@@ -266 +382 @@
 
     // set the model normalization
-    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
-      return false;
-    }
+    // if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+    //   return false;
+    // }
+    PAR[PM_PAR_I0] = peak;
 
     // set the model position