Changeset 5772 for trunk/psphot/src/psphotApplyPSF.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/psphotApplyPSF.c (modified) (3 diffs)

trunk/psphot/src/psphotApplyPSF.c

@@ -2 +2 @@
 
 // fit psf model to all objects 
-// PSFSTAR objects will be refitted
-// run this function to a specific flux limit?
-
 bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
 { 
-    bool  status;
     float x;
     float y;
-    int   Nfit = 0;
-    int   Nsub = 0;
+    int   Nfit  = 0;
+    int   Nsub  = 0;
     int   Niter = 0;
 
     psTimerStart ("psphot");
 
-    // we may set this differently here from the value used to mark likely saturated stars
-    float SATURATION   = psMetadataLookupF32 (&status, config, "SATURATION");
-    float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
-
-    float PSF_MIN_SN       = psMetadataLookupF32 (&status, config, "PSF_MIN_SN");
-    float PSF_MAX_CHI      = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI");
-    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
-    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
-    float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
-
-    // set the object surface-brightness limit for fitted pixels
-    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
-    psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
-
-    // this function specifies the radius at this the model hits the given flux
-    pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
+    psphotInitLimitsPSF (config);
+    psphotInitRadiusPSF (config, sky, psf->type);
 
     for (int i = 0; i < sources->n; i++) {
@@ -38 +20 @@
 
         // skip non-astronomical objects (very likely defects)
-        // XXX EAM : should we try these anyway?
+        if (source->mode  & PM_SOURCE_BLEND) continue; 
         if (source->type == PM_SOURCE_DEFECT) continue; 
         if (source->type == PM_SOURCE_SATURATED) continue;
@@ -47 +29 @@
         // set PSF parameters for this model
         pmModel *model = pmModelFromPSF (modelFLT, psf);
+        psFree (modelFLT);
+
+        source->modelPSF = model;
+
+        // sets the model radius (via source->model) and adjusts pixels as needed
+        psphotCheckRadiusPSF (imdata, source);
+
         x = model->params->data.F32[2];
         y = model->params->data.F32[3];
-        psFree (modelFLT);
 
-        // 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) + PSF_FIT_PADDING;
-        if (isnan(model->radius)) {
-          psAbort ("apply_psf_model", "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);
-        }
-
-        // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
-            
         // fit PSF model (set/unset the pixel mask)
         psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
-        status = pmSourceFitModel (source, model, true);
+        pmSourceFitModel (source, model, true);
         psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
 
-        if (!status || (model->params->data.F32[1] < 0)) {
-          psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
-          source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
-          psFree (model);
-          continue;
-        }
-        // XXX EAM : this was an attempt to correct for fit-sky biases
-        // pmModelSkyOffset (model, x, y, model->radius);
-
-        // model succeeded : tentatively keep it
-        source->modelPSF = model;
         Niter += model[0].nIter;
         Nfit ++;
 
-        // is it a good model?
-        psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION);
-        if (psphotSubtractPSF (source)) {
-          Nsub ++;
+        // check if model fit is acceptable
+        if (psphotEvalPSF (source)) {
+            pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
+            source->mode |= PM_SOURCE_SUBTRACTED;
+            Nsub ++;
         }
     }
 
     psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
-    psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub);
+    psLogMsg ("psphot", 4, "subtracted %d PSF models\n", Nsub);
     return (true);
 }