Changeset 33415 for branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotExtendedSourceFits.c (modified) (18 diffs)

branches/meh_branches/ppstack_test/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceFits.c

@@ -6 +6 @@
     bool status = true;
 
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Fits ---");
+
     // select the appropriate recipe information
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -12 +15 @@
     // perform full extended source non-linear fits?
     if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n");
         return true;
     }
@@ -45 +48 @@
     int Nfail = 0;
 
+    psphotSersicModelClassInit();
+
     psTimerStart ("psphot.extended");
 
@@ -103 +108 @@
 
       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
       }
 
@@ -213 +218 @@
 
         // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
+        if (!psThreadPoolWait (false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             psFree(AnalysisRegion);
@@ -247 +252 @@
     psFree(AnalysisRegion);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psphotSersicModelClassCleanup();
+
+    psphotVisualShowResidualImage (readout, false);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "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);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
@@ -265 +274 @@
     int Nplain = 0;
     int NplainPass = 0;
-    bool savePics = false;
-    float radius;
+    float fitRadius, windowRadius;
     psScalar *scalar = NULL;
-    pmMoments psfMoments;
 
     // arguments: readout, sources, models, region, psfSize, maskVal, markVal
@@ -275 +282 @@
     psMetadata *models      = job->args->data[2];
     psRegion *region        = job->args->data[3];
-    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[4],S32);
     psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     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
@@ -299 +304 @@
 
         // skip PSF-like and non-astronomical objects
-        if (source->type == PM_SOURCE_TYPE_STAR) continue;
+        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
@@ -309 +314 @@
         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 +320 @@
         Next ++;
 
-        // set the radius based on the footprint (also sets the mask pixels)
-        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) {
-            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
-            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
-        psfMoments = *source->moments;
-        if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
-            // subtract the best fit from the object, leave local sky
-            fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
-            *source->moments = psfMoments;
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-            // XXX raise an error flag of some kind
-            continue;
-        }
-
-        // save the modelFlux here in case we need to subtract it (for failure)
-        psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-        if (!modelFluxStart) {
-            pmSourceCacheModel (source, maskVal);
-            modelFluxStart = psMemIncrRefCounter (source->modelFlux);
-            if (!modelFluxStart) {
-                fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf);
-                *source->moments = psfMoments;
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                // XXX raise an error of some kind?
-                continue;
-            }
-        }
-        
-        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);
+        // set the fit radius based on the first radial moment (also sets the mask pixels)
+        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
+
+        // UPDATE : we have changed the moments calculation.  There is now an iteration within 
+        // psphotKronMasked to determine moments appropriate for a larger object.  The values
+        // Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius.  The other
+        // values are left at the psf-based values.
 
         // allocate the array to store the model fits
@@ -380 +343 @@
 
           // check the SNlim and skip model if source is too faint
-          float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
+          float FIT_SN_LIM = psMetadataLookupF32 (&status, model, "SNLIM_VALUE");
           assert (status);
 
           // limit selection to some SN limit
           // assert (source->peak); // how can a source not have a peak?
-          if (sqrt(source->peak->detValue) < SNlim) {
+          // limit selection to some SN limit
+          bool skipSource = false;
+          if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+              skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
+          } else {
+              skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
+          }
+          if (skipSource) {
               Nfaint ++;
               continue;
@@ -398 +368 @@
           assert (status);
 
-          // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-          // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
-
           // fit the model as convolved or not
           pmModel *modelFit = NULL;
           if (convolved) {
-              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize);
+              // NOTE : 4 more allocs to here
+              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize); // NOTE : 2313 allocs in here
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
@@ -410 +378 @@
                   continue;
               }
-              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
+                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nconvolve ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -425 +394 @@
                   continue;
               }
-              pmSourceCacheModel (source, maskVal);
-              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
+                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nplain ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -434 +403 @@
           }
 
-          // save each of the model flux images and store the best
-          psArrayAdd (modelFluxes, 4, source->modelFlux);
+          // XXX deprecate?
+          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
+          // XXX extFitPars->Mxx = source->moments->Mxx;
+          // XXX extFitPars->Mxy = source->moments->Mxy;
+          // XXX extFitPars->Myy = source->moments->Myy;
+          // XXX extFitPars->Mrf = source->moments->Mrf;
+          // XXX extFitPars->Mrh = source->moments->Mrh;
+          // XXX 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)
+          // XXX 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);
 
@@ -463 +444 @@
         }
 
+        // 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);
 
-          // replace original model, subtract it
-          psFree (source->modelFlux);
-          source->modelFlux = modelFluxStart;
-
-          *source->moments = psfMoments;
+          // ensure the modelEXT is cached
+          pmSourceCacheModel (source, maskVal);
+
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-          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;
         }
@@ -493 +463 @@
         source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
 
-        // save the modelFlux for the best model
-        psFree (source->modelFlux);
-        source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]);
-
-        // replace the original moments
-        *source->moments = psfMoments;
+        // adjust the window so the subtraction covers the faint wings
+        psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
+
+        // cache the model flux
+        if (source->modelEXT->isPCM) {
+            pmPCMCacheModel (source, maskVal, psfSize);
+        } else {
+            pmSourceCacheModel (source, maskVal);
+        }
 
         // subtract the best fit from the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-        // the initial model flux is no longer needed
-        psFree (modelFluxStart);
-        psFree (modelFluxes);
-
         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)

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest