Changeset 35322 for branches/eam_branches/ipp-20130307/psphot/src/psphotModelTestReadout.c

Timestamp:

Mar 19, 2013, 2:50:27 PM (13 years ago)

Author:

eugene

Message:

add mark & mask to model guess code; add -moments-radius and -psf to psphotModelTest; pass the fit tolerance info to the extended source fits; allow psphotModelTest to do trails and psf-convolved fits

File:

• branches/eam_branches/ipp-20130307/psphot/src/psphotModelTestReadout.c (modified) (10 diffs)

branches/eam_branches/ipp-20130307/psphot/src/psphotModelTestReadout.c

@@ -5 +5 @@
 
     bool status;
+
+    psTimerStart ("modelTest");
 
     pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
@@ -41 +43 @@
     }
 
+    // load the psf model, if suppled.  FWHM_MAJ,FWHM_MIN,etc are determined and saved on
+    // readout->analysis. NOTE: this function currently only loads from PSPHOT.PSF.LOAD
+    if (!psphotLoadPSF (config, view, filerule)) { // ??? need to supply 2 ?
+        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
     float MIN_KRON_RADIUS = 5.0;
 
@@ -53 +62 @@
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
 
     // find the currently selected readout
@@ -84 +97 @@
 
     // define the source of interest
-   float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
+    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
@@ -99 +112 @@
     source->type = PM_SOURCE_TYPE_EXTENDED;
 
+    bool TEST_FIT_CONVOLVED = psMetadataLookupBool (&status, recipe, "TEST_FIT_CONVOLVED");
+
     // find the local sky
     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     if (!status) psAbort("pmSourceLocalSky error");
 
-    { 
+    bool testMomentsRadius = psMetadataLookupBool (&status, config->arguments, "TEST_MOMENTS_RADIUS");
+    if (testMomentsRadius) { 
         // XXX I want to test an iterative aperture for brighter sources
         float radius = mRADIUS;
@@ -109 +125 @@
 
             // get the source moments
-            status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+            status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal);
             if (!status) psAbort("psSourceMoments error");
 
@@ -121 +137 @@
 
     // get the source moments
-    status = pmSourceMoments (source, mRADIUS, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+    status = pmSourceMoments (source, mRADIUS, 0.25*mRADIUS, 0.0, MIN_KRON_RADIUS, maskVal);
     if (!status) psAbort("psSourceMoments error");
 
@@ -135 +151 @@
     fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
 
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    if (psf) {
+      // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
+      source->peak->rawFlux = source->moments->Peak;
+      float Io = source->moments->Peak;
+      source->modelPSF = pmModelFromPSFforXY(psf, xObj, yObj, Io);
+    }
+
     // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
+    pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
     source->modelEXT = model;
 
@@ -159 +183 @@
     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
 
-    // list model input shape
-    psEllipseShape shape;
-    shape.sx  = 1.4 / model->params->data.F32[4];
-    shape.sy  = 1.4 / model->params->data.F32[5];
-    shape.sxy = model->params->data.F32[6];
-    axes = psEllipseShapeToAxes (shape, 20.0);
-
-    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        fprintf (stderr, "guess: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]);
+    } else {
+        // list model input shape
+        psEllipseShape shape;
+        shape.sx  = 1.4 / model->params->data.F32[4];
+        shape.sy  = 1.4 / model->params->data.F32[5];
+        shape.sxy = model->params->data.F32[6];
+        axes = psEllipseShapeToAxes (shape, 20.0);
+        fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+    }
 
     fprintf (stderr, "input parameters: \n");
@@ -190 +217 @@
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
 
-    status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        fitOptions->mode = PM_SOURCE_FIT_NO_INDEX;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        fitOptions->mode = PM_SOURCE_FIT_TRAIL;
+    }
+
+    if (TEST_FIT_CONVOLVED) {
+      pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+      if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        psphotSersicModelClassGuessPCM (pcm, source);
+      } else {
+        pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
+      }
+      pmPCMupdate(pcm, source, fitOptions, model);
+      pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
+      psFree (pcm);
+    } else {
+      status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    }
 
     // measure the source mags