Changeset 31990 for branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Aug 3, 2011, 10:30:02 AM (15 years ago)

Author:

eugene

Message:

add new window function using first radial moment; use new sersic guess function

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (modified) (15 diffs)

branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c

@@ -45 +45 @@
     int Nfail = 0;
 
+    psphotSersicModelClassInit();
+
     psTimerStart ("psphot.extended");
 
@@ -103 +105 @@
 
       if (item->type != PS_DATA_METADATA) {
+        // XXX we could cull the bad entries or build a validated model folder
         psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
-        // XXX we could cull the bad entries or build a validated model folder
       }
 
@@ -247 +249 @@
     psFree(AnalysisRegion);
 
+    psphotSersicModelClassCleanup();
+
     psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
@@ -265 +269 @@
     int Nplain = 0;
     int NplainPass = 0;
-    bool savePics = false;
-    float radius;
+    float fitRadius, windowRadius;
     psScalar *scalar = NULL;
     pmMoments psfMoments;
@@ -279 +282 @@
     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
 
-    // pthread_t tid = pthread_self();     // Thread identifier
-
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
@@ -309 +310 @@
         if (source->peak->y > region->y1) continue;
 
-        // if model is NULL, we don't have a starting guess
-        // XXX use the parameters guessed from moments
-        // if (source->modelEXT == NULL) continue;
-
-        // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
-
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -321 +316 @@
         Next ++;
 
-        // set the radius based on the footprint (also sets the mask pixels)
-        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
+        // set the radius based on the first radial moment (also sets the mask pixels)
+        if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) {
             fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf);
             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error of some kind
+            // XXX raise an error of some kind or set a flag bit
             continue;
         }
 
-        // XXX note that this changes the source moments that are published...
-        // recalculate the source moments using the larger extended-source moments radius
-        // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-        // this uses the footprint to judge both radius and aperture?
-        // XXX save the psf-based moments for output
+        // Recalculate the source moments using the larger extended-source moments radius.  At
+        // this stage, skip Gaussian windowing, and do not clip pixels by S/N.  Save the
+        // psf-based moments for output
         psfMoments = *source->moments;
-        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
+        if (!pmSourceMoments (source, windowRadius, 0.0, 0.0, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
             // subtract the best fit from the object, leave local sky
             fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
@@ -344 +337 @@
         }
 
-        // save the modelFlux here in case we need to subtract it (for failure)
+        // save the PSF-based modelFlux here in case we need to subtract it (for failure).  If
+        // it does not exist, attempt to generate it.  Give up if we cannot even do that.
         psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
         if (!modelFluxStart) {
@@ -358 +352 @@
         }
         
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xp.fits");
-        }
-
         // array to store the pointers to the model flux images while the models are being fitted
         psArray *modelFluxes = psArrayAllocEmpty (models->list->n);
@@ -368 +358 @@
         if (source->modelFits == NULL) {
             source->modelFits = psArrayAllocEmpty (models->list->n);
+        }
+        // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
+        if (source->extFitPars == NULL) {
+            source->extFitPars = psArrayAllocEmpty (models->list->n);
         }
 
@@ -397 +391 @@
           bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE");
           assert (status);
-
-          // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-          // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
-
-          // XXX this does not make sense in a threaded context
-          psTimerStart ("psphot.extended.fit");
 
           // fit the model as convolved or not
@@ -439 +427 @@
           }
 
-          psLogMsg ("psphot", PS_LOG_INFO, "model fit: %7.5f sec, %5d npix, %2d iter, %s type\n", psTimerMark ("psphot.extended.fit"), modelFit->nPix, modelFit->nIter, pmModelClassGetName (modelFit->type));
-
-          // save each of the model flux images and store the best
+          // save each of the model flux images for now, and (below) store the best
           psArrayAdd (modelFluxes, 4, source->modelFlux);
 
+          pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
+          extFitPars->Mxx = source->moments->Mxx;
+          extFitPars->Mxy = source->moments->Mxy;
+          extFitPars->Myy = source->moments->Myy;
+          extFitPars->Mrf = source->moments->Mrf;
+          extFitPars->Mrh = source->moments->Mrh;
+          extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
+          
+          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
+          extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
+
+          psArrayAdd (source->extFitPars, 4, extFitPars);
+          psFree (extFitPars);
+
           // test for fit quality / result
-          modelFit->fitRadius = radius;
+          modelFit->fitRadius = fitRadius;
           psArrayAdd (source->modelFits, 4, modelFit);
 
@@ -470 +470 @@
         }
 
+        // clear the circular mask
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
         if (minModel == -1) {
-          // re-subtract the object, leave local sky
+          // no valid extended fit; re-subtract the object, leave local sky
           psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
 
@@ -483 +486 @@
           psFree (modelFluxes);
 
-          if (savePics) {
-              psphotSaveImage (NULL, readout->image, "image.xp.fits");
-              char key[10];
-              fprintf (stdout, "continue? ");
-              if (!fgets (key, 8, stdin)) {
-                  psWarning("Couldn't read anything.");
-              } else if (key[0] == 'n') {
-                  savePics = false;
-              }
-          }
           continue;
         }
@@ -516 +509 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-
-        if (savePics) {
-          psphotSaveImage (NULL, readout->image, "image.xm.fits");
-          char key[10];
-          fprintf (stdout, "continue? ");
-          if (!fgets (key, 8, stdin)) {
-              psWarning("Couldn't read anything.");
-          } else if (key[0] == 'n') {
-              savePics = false;
-          }
-        }
     }
     psFree (fitOptions);
-
-    // fprintf (stderr, "xfit : tried %ld objects\n", sources->n);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)