Changeset 27838 for branches/tap_branches/psphot/src/psphotChoosePSF.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotChoosePSF.c (modified) (16 diffs)

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotChoosePSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// generate a PSF model for inputs without PSF models already loaded
+bool psphotChoosePSF (pmConfig *config, const pmFPAview *view)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
 // try PSF models and select best option
-pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
 
     bool status;
 
     psTimerStart ("psphot.choose.psf");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    // do not generate a PSF if we already were supplied one
+    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index);
+        return true;
+    }
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->newSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model");
+        return false;
+    }
 
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
+    psAssert (maskVal, "missing mask value?");
 
     // bit-mask to mark pixels not used in analysis
     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
+    psAssert (markVal, "missing mark value?");
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -61 +112 @@
     // assert (status);
 
-    // We have calculated a Gaussian window function, use that for both the PSF fit radius and
-    // the aperture radius (scaling SIGMA)
-    float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    // values on readout->analysis if we have calculated a Gaussian window function, use that
+    // for both the PSF fit radius and the aperture radius (scaling SIGMA), otherwise base the value on the recipe value for MOMENTS_GAUSS_SIGMA
+    float gaussSigma = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
+    if (!status) {
+        gaussSigma = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    }
     float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
     float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
@@ -70 +124 @@
 
     // XXX use the same radii for standard analysis as for the PSF creation
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
-
-    // 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);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
 
     // dimensions of the field for which the PSF is defined
@@ -116 +166 @@
     psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
 
+    if (stars->n < 50) {
+        // ROBUST is too agressive if we only have a small number of PSF stars
+        options->stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    } else {
+        options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    }
+
     // get the list pointers for the PSF_MODEL entries
     psArray *modelNames = NULL;
@@ -155 +212 @@
         bool status = true;
         status &= psphotMakeFluxScale (readout->image, recipe, psf);
-        status &= psphotPSFstats (readout, recipe, psf);
+        status &= psphotPSFstats (readout, psf);
         if (!status) {
             psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF");
@@ -230 +287 @@
         bool status = true;
         status &= psphotMakeFluxScale (readout->image, recipe, psf);
-        status &= psphotPSFstats (readout, recipe, psf);
+        status &= psphotPSFstats (readout, psf);
         if (!status) {
             psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF");
@@ -287 +344 @@
     psFree (models);
 
-    if (!psphotPSFstats (readout, recipe, psf)) {
+    if (!psphotPSFstats (readout, psf)) {
         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
         psFree(options);
@@ -299 +356 @@
 
     psFree (options);
-    return (psf);
+
+    // save PSF on readout->analysis
+    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
+        return false;
+    }
+    psFree (psf); // XXX double-check
+
+    // move into psphotChoosePSF
+    psphotVisualShowPSFModel (readout, psf);
+
+    return true;
 }
 
 // measure average parameters of the PSF model
-bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
+bool psphotPSFstats (pmReadout *readout, pmPSF *psf) {
 
     psEllipseShape shape;
@@ -309 +377 @@
 
     PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
     PS_ASSERT_PTR_NON_NULL(psf, false);
 
@@ -329 +396 @@
             pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
             if (!modelPSF) {
-                fprintf (stderr, "?");
-                continue;
-            }
+                fprintf (stderr, "?");
+                continue;
+            }
 
             // get the model full-width at half-max
@@ -345 +412 @@
             float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
             if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
-                fprintf (stderr, "!");
-                continue;
-            }
+                fprintf (stderr, "!");
+                continue;
+            }
             psVectorAppend (fwhmMajor, FWHM_MAJOR);
             psVectorAppend (fwhmMinor, FWHM_MINOR);
@@ -359 +426 @@
     }
 
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
+
+    float fwhmMaj = stats->sampleMean;  // FWHM on major axis
 
     if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
@@ -368 +437 @@
         return false;
     }
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
-    psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
+
+    float fwhmMin = stats->sampleMean;  // FWHM on minor axis
+    if (readout->parent) {
+        pmChip *chip = readout->parent->parent; // Parent chip
+        psAssert(chip, "Cell should be attached to a chip.");
+        psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip
+        item->data.F32 = 0.5 * (fwhmMaj + fwhmMin);
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
 
     psFree (fwhmMajor);
@@ -384 +461 @@
 
 // determine approximate PSF shape parameters based on the moments
-bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) {
+bool psphotMomentsStats (pmReadout *readout, psArray *sources) {
 
     // without the PSF model, we can only rely on the source->moments
@@ -398 +475 @@
 
     PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
     PS_ASSERT_PTR_NON_NULL(sources, false);
 
@@ -426 +502 @@
     FWHM_T /= FWHM_N;
 
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           FWHM_T);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0);
-    psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false);
+    if (readout->parent) {
+        pmChip *chip = readout->parent->parent; // Parent chip
+        psAssert(chip, "Cell should be attached to a chip.");
+        psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip
+        item->data.F32 = 0.5 * (FWHM_X + FWHM_Y);
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           FWHM_T);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0);
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false);
 
     return true;

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack