Changeset 35321

Timestamp:

Mar 19, 2013, 2:45:32 PM (13 years ago)

Author:

eugene

Message:

add mask and mark values to model guess code; add new angle guess code

Location:

branches/eam_branches/ipp-20130307/psModules/src/objects

Files:

• models/pmModel_DEV.c (modified) (1 diff)
• models/pmModel_DEV.h (modified) (1 diff)
• models/pmModel_EXP.c (modified) (1 diff)
• models/pmModel_EXP.h (modified) (1 diff)
• models/pmModel_GAUSS.c (modified) (1 diff)
• models/pmModel_GAUSS.h (modified) (1 diff)
• models/pmModel_PGAUSS.c (modified) (1 diff)
• models/pmModel_PGAUSS.h (modified) (1 diff)
• models/pmModel_PS1_V1.c (modified) (1 diff)
• models/pmModel_PS1_V1.h (modified) (1 diff)
• models/pmModel_QGAUSS.c (modified) (1 diff)
• models/pmModel_QGAUSS.h (modified) (1 diff)
• models/pmModel_RGAUSS.c (modified) (1 diff)
• models/pmModel_RGAUSS.h (modified) (1 diff)
• models/pmModel_SERSIC.c (modified) (1 diff)
• models/pmModel_SERSIC.h (modified) (1 diff)
• models/pmModel_TRAIL.c (modified) (4 diffs)
• models/pmModel_TRAIL.h (modified) (1 diff)
• pmModelFuncs.h (modified) (1 diff)
• pmPSFtryFitEXT.c (modified) (1 diff)
• pmSourceFitPCM.c (modified) (1 diff)
• pmSourceUtils.c (modified) (3 diffs)
• pmSourceUtils.h (modified) (1 diff)

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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;

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/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

branches/eam_branches/ipp-20130307/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);

branches/eam_branches/ipp-20130307/psModules/src/objects/pmModelFuncs.h

@@ -107 +107 @@
 //  This function provides the model guess parameters based on the details of
 //  the given source.
-typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 
 //  This function constructs the PSF model for the given source based on the

branches/eam_branches/ipp-20130307/psModules/src/objects/pmPSFtryFitEXT.c

@@ -83 +83 @@
         }
 
-        source->modelEXT = pmSourceModelGuess (source, options->type);
+        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
         if (source->modelEXT == NULL) {
             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);

branches/eam_branches/ipp-20130307/psModules/src/objects/pmSourceFitPCM.c

@@ -145 +145 @@
 bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    if (!pcm->modelConv->modelGuess(pcm->modelConv, source)) {
+  if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) {
         return false;
     }

branches/eam_branches/ipp-20130307/psModules/src/objects/pmSourceUtils.c

@@ -46 +46 @@
 
 /******************************************************************************
-    pmSourceModelGuess(source, model): This function allocates a new
-    pmModel structure based on the given modelType specified in the argument list.
-    The corresponding pmModelGuess function is returned, and used to
-    supply the values of the params array in the pmModel structure.
+    pmSourceModelGuess(source, model, maskVal, markVal): This function allocates a new
+    pmModel structure based on the given modelType specified in the argument list.  The
+    corresponding pmModelGuess function is returned, and used to supply the values of the
+    params array in the pmModel structure.
  
     XXX: Many parameters are based on the src->moments structure, which is in
@@ -55 +55 @@
     functions will be in image, not subImage coords.  Remember this.
 *****************************************************************************/
-pmModel *pmSourceModelGuess(pmSource *source,
-                            pmModelType modelType)
+pmModel *pmSourceModelGuess(pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal)
 {
     psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
@@ -65 +64 @@
     pmModel *model = pmModelAlloc(modelType);
 
-    if (!model->modelGuess(model, source)) {
+    if (!model->modelGuess(model, source, maskVal, markVal)) {
         psFree (model);
         return NULL;

branches/eam_branches/ipp-20130307/psModules/src/objects/pmSourceUtils.h

@@ -27 +27 @@
 pmModel *pmSourceModelGuess(
     pmSource *source,   ///< The input pmSource
-    pmModelType model   ///< The type of model to be created.
+    pmModelType model,   ///< The type of model to be created.
+    psImageMaskType maskVal, 
+    psImageMaskType markVal
 );