Changeset 4574 for trunk/psphot/src/fit_galaxies.c

Timestamp:

Jul 16, 2005, 6:02:29 AM (21 years ago)

Author:

eugene

Message:

improved galaxy model guesses, better handling of errors

File:

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

trunk/psphot/src/fit_galaxies.c

@@ -3 +3 @@
 // fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
 
-bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources) 
+bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
 { 
     bool  status;
+    float x;
+    float y;
+    float sky;
     int   Nfit = 0;
+    int   Nfail = 0;
     int   Niter = 0;
 
+    float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
     float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
     float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
-    psModelType modelType = psModelSetType ("PS_MODEL_PGAUSS");
+    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
+    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+
+    float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
+
+    psModelType   modelType   = psModelSetType ("PS_MODEL_RGAUSS");
+    psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
+
+    psTraceSetLevel (".psModules.pmSourceMoments", 5);
 
     psTimerStart ("psphot");
@@ -19 +33 @@
         if (source->moments->SN < snFaint) continue;
 
+        // recalculate the source moments using the galaxy radius (larger)
+        status = pmSourceMoments (source, MOMENT_R);
+
+        // use the source moments, etc to guess basic model parameters
+        psModel  *model  = pmSourceModelGuess (source, modelType); 
+
+        x = model->params->data.F32[2];
+        y = model->params->data.F32[3];
+
         // need a better model guess and a better radius choice
         // when radius is not fixed, we will need to check if new radius fits on image
-# if (0)
+
+        // set the fit radius based on the object flux limit and the model
+        // FLUX_LIMIT should be set based on local sky model (not global median)
+        // model->radius = 25.0;
+        model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
+        if (isnan(model->radius)) {
+          fprintf (stderr, "error in radius\n");
+          continue;
+        }
         if (model->radius > OUTER) {
           // allocate image, noise, mask arrays for each peak (square of radius OUTER)
-          pmSourceDefinePixels (source, imdata, x, y, OUTER);
+          pmSourceDefinePixels (source, imdata, x, y, model->radius);
         }
-# endif
-        psModel  *model  = pmSourceModelGuess (source, modelType); 
 
         // fit as FLT, not PSF (skip poor fits)
-        if (!pmSourceFitModel (source, model, false)) continue;
+        psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
+        status = pmSourceFitModel (source, model, false);
+        psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
+        if (!status || (model->params->data.F32[1] < 0)) {
+          // if the fit fails, we need to change the classification
+          psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
+          source->type = PS_SOURCE_OTHER;  // better choice?
+          Nfail ++;
+          continue;
+        }
+
         source->modelFLT = model;
-
         Niter += model[0].nIter;
         Nfit++;
+
+        // subtract object, leave local sky
+        sky = model->params->data.F32[0];
+        model->params->data.F32[0] = 0;
+        pmSourceSubModel (source->pixels, source->mask, model, false);
+        model->params->data.F32[0] = sky;
+
     }
-    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
     return (true);
 }