Changeset 32348 for trunk/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Sep 6, 2011, 1:32:31 PM (15 years ago)

Author:

eugene

Message:

• clarified the output to show the stages and their timing more clearly
• psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
• psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
• for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
• extended source analysis (petrosians) is threaded
• iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
• the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
• dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
• fake sources inserted for the efficiency analysis are removed at the end of the analysis
• radial apertures is now threaded
• edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
• new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
• define the fit region and window region based on the 1st radial moments
• use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
• remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotExtendedSourceFits.c (modified) (17 diffs)

trunk/psphot/src

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

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

trunk/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
       }
 
@@ -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);
+
+        // 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;
@@ -397 +367 @@
           bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE");
           assert (status);
-
-          // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-          // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6);
 
           // fit the model as convolved or not
@@ -410 +377 @@
                   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 +393 @@
                   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 +402 @@
           }
 
-          // 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 +443 @@
         }
 
+        // 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 +462 @@
         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