Changeset 18865

Timestamp:

Aug 1, 2008, 2:50:07 PM (18 years ago)

Author:

Paul Price

Message:

Adding PSF_FIT_ITER, PSF_FIT_TOL, EXT_FIT_ITER and EXT_FIT_TOL recipe
parameters to control the non-linear source fitting. We want to
separate the PSF fits and the general extended source fits: we can
afford to spend more time on the PSF fits to make sure they're right
(and that they succeed), whereas we're willing to throw away other
sources if they don't fit.

Location:

trunk/psphot/src

Files:

• psphotChoosePSF.c (modified) (3 diffs)
• psphotReadout.c (modified) (7 diffs)

trunk/psphot/src/psphotChoosePSF.c

@@ -64 +64 @@
     assert (status);
 
-    // XXX ROBUST seems to be too agressive given the small numbers.  
+    // XXX ROBUST seems to be too agressive given the small numbers.
     // options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     options->stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
@@ -74 +74 @@
     options->psfFieldYo = readout->image->row0;
 
-    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
+    int fitIter = psMetadataLookupS32(&status, recipe, "PSF_FIT_ITER"); // Maximum number of fit iterations
+    if (!status || fitIter <= 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_ITER is not positive");
+        return false;
+    }
+    float fitTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitTol) || fitTol <= 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_TOL is not positive");
+        return false;
+    }
+    pmSourceFitModelInit(fitIter, fitTol, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
 
     psArray *stars = psArrayAllocEmpty (sources->n);
@@ -188 +198 @@
         bool mdok;                      // Status of MD lookup
         if (!psMetadataLookupBool(&mdok, recipe, PSPHOT_RECIPE_PSF_FAKE_ALLOW)) {
-            psFree (modelNames);
-            psFree(options);
+            psFree (modelNames);
+            psFree(options);
             return NULL;
         }
 
         // generate a psf model using the first selection
-        psLogMsg ("psphot.pspsf", PS_LOG_INFO, "Using guess PSF model");
+        psLogMsg ("psphot.pspsf", PS_LOG_INFO, "Using guess PSF model");
 
         // unset the PSFSTAR flags (none are used):

trunk/psphot/src/psphotReadout.c

@@ -62 +62 @@
         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
     }
- 
+
     // XXX test write out the footprint image
     if (0) {
-        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
-        psphotSaveImage (NULL, footprintImage, "footprints.1.fits");
-        psFree (footprintImage);
+        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
+        psphotSaveImage (NULL, footprintImage, "footprints.1.fits");
+        psFree (footprintImage);
     }
 
@@ -112 +112 @@
     }
 
+    // Define source fitting parameters for everything that follows PSF fits
+    // This allows different parameters to be used from the PSF fitting
+    {
+        bool status = true;             // Status of MD lookup
+        int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations
+        if (!status || fitIter <= 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "EXT_FIT_ITER is not positive");
+            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+        }
+        float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitTol) || fitTol <= 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "EXT_FIT_TOL is not positive");
+            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+        }
+        bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors?
+        if (!status) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "POISSON.ERRORS.PHOT.LMM is not defined");
+            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+        }
+        float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
+        if (!status || !isfinite(skySig) || skySig <= 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "SKY_SIG is not positive");
+            return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+        }
+        pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson);
+    }
+
     // include externally-supplied sources
     psphotLoadExtSources (config, view, sources);
@@ -124 +151 @@
     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
 
-    if (0) { 
-        FILE *out = fopen ("out.pass1.dat", "w");
-        
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            if (!source) continue;
-            pmPeak *peak = source->peak;
-            if (!peak) continue;
-            pmModel *model = source->modelPSF;
-            if (!model) continue;
-            if (!model->params) continue;
-            
-            fprintf (out, "%d %f %f  %f %f  %f %f\n", 
-                     i, peak->xf, peak->yf, peak->flux, peak->SN,
-                     model->params->data.F32[PM_PAR_I0],
-                     model->dparams->data.F32[PM_PAR_I0]);
-        }
-        fclose (out);
+    if (0) {
+        FILE *out = fopen ("out.pass1.dat", "w");
+
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            if (!source) continue;
+            pmPeak *peak = source->peak;
+            if (!peak) continue;
+            pmModel *model = source->modelPSF;
+            if (!model) continue;
+            if (!model->params) continue;
+
+            fprintf (out, "%d %f %f  %f %f  %f %f\n",
+                     i, peak->xf, peak->yf, peak->flux, peak->SN,
+                     model->params->data.F32[PM_PAR_I0],
+                     model->dparams->data.F32[PM_PAR_I0]);
+        }
+        fclose (out);
     }
 
@@ -178 +205 @@
     // XXX test write out the footprint image
     if (0) {
-        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
-        psphotSaveImage (NULL, footprintImage, "footprints.2.fits");
-        psFree (footprintImage);
+        psImage *footprintImage = pmSetFootprintArrayIDs (detections->footprints, true);
+        psphotSaveImage (NULL, footprintImage, "footprints.2.fits");
+        psFree (footprintImage);
     }
 
@@ -208 +235 @@
     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
 
-    if (0) { 
-        FILE *out = fopen ("out.pass2.dat", "w");
-        
-        for (int i = 0; i < sources->n; i++) {
-            pmSource *source = sources->data[i];
-            if (!source) continue;
-            pmPeak *peak = source->peak;
-            if (!peak) continue;
-            pmModel *model = source->modelPSF;
-            if (!model) continue;
-            if (!model->params) continue;
-            
-            fprintf (out, "%d %f %f  %f %f  %f %f\n", 
-                     i, peak->xf, peak->yf, peak->flux, peak->SN,
-                     model->params->data.F32[PM_PAR_I0],
-                     model->dparams->data.F32[PM_PAR_I0]);
-        }
-        fclose (out);
+    if (0) {
+        FILE *out = fopen ("out.pass2.dat", "w");
+
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            if (!source) continue;
+            pmPeak *peak = source->peak;
+            if (!peak) continue;
+            pmModel *model = source->modelPSF;
+            if (!model) continue;
+            if (!model->params) continue;
+
+            fprintf (out, "%d %f %f  %f %f  %f %f\n",
+                     i, peak->xf, peak->yf, peak->flux, peak->SN,
+                     model->params->data.F32[PM_PAR_I0],
+                     model->dparams->data.F32[PM_PAR_I0]);
+        }
+        fclose (out);
     }
 
@@ -234 +261 @@
 
     if (0) {
-        psphotSaveImage (NULL, readout->mask, "mask.fits");
-    }
-        
+        psphotSaveImage (NULL, readout->mask, "mask.fits");
+    }
+
     psphotExtendedSourceAnalysis (readout, sources, recipe);
 
@@ -243 +270 @@
 finish:
 
-    // plot positive sources 
-    // psphotSourcePlots (readout, sources, recipe); 
+    // plot positive sources
+    // psphotSourcePlots (readout, sources, recipe);
 
     // measure aperture photometry corrections