Changeset 35560 for trunk/psModules/src/objects/models

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

Location:

trunk/psModules/src/objects/models

Files:

• pmModel_DEV.c (modified) (1 diff)
• pmModel_DEV.h (modified) (1 diff)
• pmModel_EXP.c (modified) (1 diff)
• pmModel_EXP.h (modified) (1 diff)
• pmModel_GAUSS.c (modified) (1 diff)
• pmModel_GAUSS.h (modified) (1 diff)
• pmModel_PGAUSS.c (modified) (1 diff)
• pmModel_PGAUSS.h (modified) (1 diff)
• pmModel_PS1_V1.c (modified) (1 diff)
• pmModel_PS1_V1.h (modified) (1 diff)
• pmModel_QGAUSS.c (modified) (1 diff)
• pmModel_QGAUSS.h (modified) (1 diff)
• pmModel_RGAUSS.c (modified) (1 diff)
• pmModel_RGAUSS.h (modified) (1 diff)
• pmModel_SERSIC.c (modified) (1 diff)
• pmModel_SERSIC.h (modified) (1 diff)
• pmModel_TRAIL.c (modified) (4 diffs)
• pmModel_TRAIL.h (modified) (1 diff)

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

@@ -295 +295 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_DEV.h

@@ -5 +5 @@
 psF32 pmModelFunc_DEV(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_DEV(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_DEV(pmModel *model, pmSource *source);
+bool pmModelGuess_DEV(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_DEV(const psVector *params);
 psF64 pmModelRadius_DEV(const psVector *params, psF64 flux);

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

@@ -287 +287 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_EXP.h

@@ -5 +5 @@
 psF32 pmModelFunc_EXP(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_EXP(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_EXP(pmModel *model, pmSource *source);
+bool pmModelGuess_EXP(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_EXP(const psVector *params);
 psF64 pmModelRadius_EXP(const psVector *params, psF64 flux);

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

@@ -192 +192 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_GAUSS.h

@@ -5 +5 @@
 psF32 pmModelFunc_GAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_GAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_GAUSS(pmModel *model, pmSource *source);
+bool pmModelGuess_GAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_GAUSS(const psVector *params);
 psF64 pmModelRadius_GAUSS(const psVector *params, psF64 flux);

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

@@ -193 +193 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_PGAUSS.h

@@ -5 +5 @@
 psF32 pmModelFunc_PGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_PGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_PGAUSS(pmModel *model, pmSource *source);
+bool pmModelGuess_PGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_PGAUSS(const psVector *params);
 psF64 pmModelRadius_PGAUSS(const psVector *params, psF64 flux);

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

@@ -212 +212 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_PS1_V1.h

@@ -5 +5 @@
 psF32 pmModelFunc_PS1_V1(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_PS1_V1(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_PS1_V1(pmModel *model, pmSource *source);
+bool pmModelGuess_PS1_V1(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_PS1_V1(const psVector *params);
 psF64 pmModelRadius_PS1_V1(const psVector *params, psF64 flux);

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

@@ -213 +213 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_QGAUSS.h

@@ -5 +5 @@
 psF32 pmModelFunc_QGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_QGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_QGAUSS(pmModel *model, pmSource *source);
+bool pmModelGuess_QGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_QGAUSS(const psVector *params);
 psF64 pmModelRadius_QGAUSS(const psVector *params, psF64 flux);

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

@@ -202 +202 @@
 // 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;

trunk/psModules/src/objects/models/pmModel_RGAUSS.h

@@ -5 +5 @@
 psF32 pmModelFunc_RGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_RGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_RGAUSS(pmModel *model, pmSource *source);
+bool pmModelGuess_RGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_RGAUSS(const psVector *params);
 psF64 pmModelRadius_RGAUSS(const psVector *params, psF64 flux);

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

@@ -307 +307 @@
 // 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)
 {
     pmMoments *moments = source->moments;

trunk/psModules/src/objects/models/pmModel_SERSIC.h

@@ -5 +5 @@
 psF32 pmModelFunc_SERSIC(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_SERSIC(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_SERSIC(pmModel *model, pmSource *source);
+bool pmModelGuess_SERSIC(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_SERSIC(const psVector *params);
 psF64 pmModelRadius_SERSIC(const psVector *params, psF64 flux);

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

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

@@ -5 +5 @@
 psF32 pmModelFunc_TRAIL(psVector *deriv, const psVector *params, const psVector *pixcoord);
 bool pmModelLimits_TRAIL(psMinConstraintMode mode, int nParam, float *params, float *beta);
-bool pmModelGuess_TRAIL(pmModel *model, pmSource *source);
+bool pmModelGuess_TRAIL(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 psF64 pmModelFlux_TRAIL(const psVector *params);
 psF64 pmModelRadius_TRAIL(const psVector *params, psF64 flux);