Changeset 5772 for trunk/psphot/src/psphotFitGalaxies.c

Timestamp:

Dec 13, 2005, 6:43:44 AM (21 years ago)

Author:

eugene

Message:

substantial work on ensemble PSF fitting, better functions for evaluation, overall code cleanups

File:

• trunk/psphot/src/psphotFitGalaxies.c (modified) (2 diffs)

trunk/psphot/src/psphotFitGalaxies.c

@@ -3 +3 @@
 bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
 { 
-    bool  status, goodfit;
+    bool  status;
     float x;
     float y;
-    int   Nfit = 0;
-    int   Nfail = 0;
+    int   Nfit  = 0;
+    int   Nsub  = 0;
     int   Niter = 0;
 
-    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    psTimerStart ("psphot");
 
     float GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
-    float GAL_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA");
-    float GAL_FIT_PADDING  = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING");
     float GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+    char       *modelName  = psMetadataLookupPtr (&status, config, "GAL_MODEL");
+    pmModelType modelType  = pmModelSetType (modelName);
 
-    float FLUX_LIMIT       = GAL_FIT_NSIGMA * skyStats->sampleStdev;
-
-    char         *modelName   = psMetadataLookupPtr (&status, config, "GAL_MODEL");
-    pmModelType   modelType   = pmModelSetType (modelName);
-    pmModelRadius modelRadius = pmModelRadius_GetFunction (modelType);
-
-    psTimerStart ("psphot");
+    psphotInitRadiusFLT (config, skyStats, modelType);
 
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
 
-        // sources which should not be fitted
-        // skip all valid stars
-        if (source->type == PM_SOURCE_PSFSTAR) continue;
-        if (source->type == PM_SOURCE_SATSTAR) continue;
-        if (source->type == PM_SOURCE_GOODSTAR) continue;
+        // only fit suspected extended sources
+        if (source->type != PM_SOURCE_GALAXY) continue;
+        if (source->mode  & PM_SOURCE_BLEND) continue; 
 
-        // skip all likely defects
-        if (source->type == PM_SOURCE_DEFECT) continue;
-        if (source->type == PM_SOURCE_SATURATED) continue;
-
-        // skip poorly fitted stars
-        if (source->type == PM_SOURCE_FAINTSTAR) continue;
-        if (source->type == PM_SOURCE_POOR_FIT_PSF) continue;
-
-        // XXX EAM when do we pick these up again?
-        if (source->moments->SN < GAL_MIN_SN) {
-          source->type = PM_SOURCE_FAINT_GALAXY;  // better choice?
-          continue;
-        }
+        // skip possible extended sources below fit threshold
+        if (source->moments->SN < GAL_MIN_SN) continue;
 
         // recalculate the source moments using the larger galaxy moments radius
         status = pmSourceMoments (source, GAL_MOMENTS_RAD);
         if (!status) {
-          source->type = PM_SOURCE_DROP_GALAXY;  // better choice?
+          source->mode |= PM_SOURCE_FAIL;  // better choice?
           continue;
         }
@@ -57 +38 @@
         // use the source moments, etc to guess basic model parameters
         pmModel  *model  = pmSourceModelGuess (source, modelType); 
+        source->modelFLT = model;
         x = model->params->data.F32[2];
         y = model->params->data.F32[3];
 
-        // set the fit radius based on the object flux limit and the model
-        // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
-        model->radius = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING;
-        if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
-
-        // check if we need to redefine the pixels
-        // XXX EAM : a better test would examine the source pixels
-        if (model->radius > OUTER_RADIUS) {
-          // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
-          psphotDefinePixels (source, imdata, x, y, model->radius);
-        }
+        // sets the model radius (via source->model) and adjusts pixels as needed
+        psphotCheckRadiusFLT (imdata, source);
 
         // fit FLT (not PSF) model (set/unset the pixel mask)
         psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
-        status = pmSourceFitModel (source, model, false);
+        pmSourceFitModel (source, model, false);
         psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
-        if (!status) {
-          // if the fit fails, we need to change the classification
-          psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
-          source->type = PM_SOURCE_FAIL_FIT_GAL;  // better choice?
-          source->modelFLT = model;
-          Nfail ++;
-          continue;
-        }
 
-        // check if model fit is acceptable
-        goodfit = pmModelFitStatus (model);
-        if (!goodfit) {
-          // if the fit fails, we need to change the classification
-          psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
-          source->type = PM_SOURCE_POOR_FIT_GAL;  // better choice?
-          source->modelFLT = model;
-          Nfail ++;
-          continue;
-        }
-
-        // accept the model
-        source->type = PM_SOURCE_GALAXY;
-        source->modelFLT = model;
         Niter += model[0].nIter;
         Nfit++;
 
-        // subtract object, leave local sky
-        pmSourceSubModel (source->pixels, source->mask, model, false, false);
+        // check if model fit is acceptable
+        if (pmModelFitStatus (model)) {
+            pmSourceSubModel (source->pixels, source->mask, model, false, false);
+            source->mode |= PM_SOURCE_SUBTRACTED;
+            Nsub ++;
+        }
     }
-    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
+    psLogMsg ("psphot", 3, "fit GAL models: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 4, "subtracted %d GAL models\n", Nsub);
     return (true);
 }