Changeset 27840 for branches/simtest_nebulous_branches/psphot/src/psphotChoosePSF.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

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

branches/simtest_nebulous_branches/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/simtest_nebulous_branches/psphot/src/psphotChoosePSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-void psphotCountPSFStars (psArray *sources) {
-
-    int nPSF = 0;
-
-    // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
-            nPSF ++;
-        }
-    }
-
-    fprintf (stderr, "N PSF: %d\n", nPSF);
+// 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
@@ -73 +109 @@
     // get the fixed PSF fit radius
     // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    assert (status);
-
-    // 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);
+    // options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
+    // assert (status);
+
+    // 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");
+    options->fitRadius = (int)(fitScale*gaussSigma);
+    options->apRadius = (int)(apScale*gaussSigma);
+
+    // XXX use the same radii for standard analysis as for the PSF creation
+    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
@@ -99 +146 @@
 
     psArray *stars = psArrayAllocEmpty (sources->n);
-
-    // psphotCountPSFStars (sources);
 
     // select the candidate PSF stars (pointers to original sources)
@@ -115 +160 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     // check that the identified psf stars sufficiently cover the region; if not, extend the
     // limits somewhat
     psphotCheckStarDistribution (stars, sources, options);
 
-    // psphotCountPSFStars (sources);
-
     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
@@ -164 +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");
@@ -239 +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");
@@ -289 +337 @@
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) {
-        psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal);
+        psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal);
     }
 
     // save only the best model;
-    // XXX we are not saving the fitted sources
-    // XXX do we want to keep them so we may optionally write them out?
     pmPSF *psf = psMemIncrRefCounter(try->psf);
     psFree (models);
 
-    if (!psphotPSFstats (readout, recipe, psf)) {
+    if (!psphotPSFstats (readout, psf)) {
         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
         psFree(options);
@@ -305 +351 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     char *modelName = pmModelClassGetName (psf->type);
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
@@ -312 +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;
@@ -322 +377 @@
 
     PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
     PS_ASSERT_PTR_NON_NULL(psf, false);
 
@@ -341 +395 @@
             // create modelPSF from this model
             pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
-            if (!modelPSF) continue;
+            if (!modelPSF) {
+                fprintf (stderr, "?");
+                continue;
+            }
 
             // get the model full-width at half-max
@@ -354 +411 @@
 
             float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue;
+            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
+                fprintf (stderr, "!");
+                continue;
+            }
             psVectorAppend (fwhmMajor, FWHM_MAJOR);
             psVectorAppend (fwhmMinor, FWHM_MINOR);
@@ -366 +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)) {
@@ -375 +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);
@@ -391 +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
@@ -405 +475 @@
 
     PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
     PS_ASSERT_PTR_NON_NULL(sources, false);
 
@@ -433 +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;

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack