Changeset 37767

Timestamp:

Dec 22, 2014, 8:48:02 AM (12 years ago)

Author:

eugene

Message:

create a new function psphotPSFstatsSources to use real sources to measure the psf stats

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotChoosePSF.c (modified) (4 diffs)
• psphotReadoutCleanup.c (modified) (1 diff)

trunk/psphot/src/psphot.h

@@ -172 +172 @@
 // in psphotChoosePSF.c:
 bool            psphotPSFstats (pmReadout *readout, pmPSF *psf);
+bool            psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psArray *sources);
 

trunk/psphot/src/psphotChoosePSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf);
 
 // generate a PSF model for inputs without PSF models already loaded
@@ -388 +389 @@
     psFree (models);
 
-    if (!psphotPSFstats (readout, psf)) {
+    // if (!psphotPSFstats (readout, psf)) {
+
+    if (!psphotPSFstatsSources (readout, sources, psf)) {
         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
         psFree(options);
@@ -466 +469 @@
 
             if (modelPSF->params->n > 7) {
-              psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
+                psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
             }
             if (modelPSF->params->n > 8) {
-              psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
+                psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
             }
             psFree (modelPSF);
 
-            // float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
-            // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
+            float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
+            fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
         }
     }
@@ -522 +525 @@
 
     if (psfExtra1->n) {
-      if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+        if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+            goto escape;
+        }
+        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    if (psfExtra2->n) {
+        if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
+            goto escape;
+        }
+        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    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);
+
+    // psLogMsg ("psphot", PS_LOG_DETAIL, "PSF angle: %f, nstars: %d", axes.theta, psf->nPSFstars);
+
+    char *psfModelName = pmModelClassGetName(psf->type);
+    psMetadataAddStr(readout->analysis,  PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName);
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK",   PS_META_REPLACE, "Valid PSF Model?", true);
+
+    int nParams = pmModelClassParameterCount(psf->type);
+    psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR",   PS_META_REPLACE, "Number of PSF parameters", nParams);
+
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
+    psFree (stats);
+    return true;
+
+escape:
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
+    psFree (stats);
+    return false;
+}
+
+// measure average parameters of the PSF model
+bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf) {
+
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(psf, false);
+
+    psImage *image = readout->image;
+    PS_ASSERT_PTR_NON_NULL(image, false);
+
+    // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc:
+    psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        float xc = source->peak->xf;
+        float yc = source->peak->yf;
+
+        // create modelPSF from this model
+        pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
+        if (!modelPSF) {
+            fprintf (stderr, "?");
+            continue;
+        }
+
+        // get the model full-width at half-max
+        float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5);
+
+        // XXX make sure this is consistent with the re-definition of PM_PAR_SXX
+        shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+        shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+        shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+        axes = psEllipseShapeToAxes (shape, 20.0);
+
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
+            fprintf (stderr, "!");
+            psFree (modelPSF);
+            continue;
+        }
+        psVectorAppend (fwhmMajor, FWHM_MAJOR);
+        psVectorAppend (fwhmMinor, FWHM_MINOR);
+
+        if (modelPSF->params->n > 7) {
+            psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
+        }
+        if (modelPSF->params->n > 8) {
+            psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
+        }
+        psFree (modelPSF);
+
+        float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
+        fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
+    }
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
+    if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
         goto escape;
-      }
-      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
-      psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    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)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
+        goto escape;
+    }
+    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) {
+
+        // we now have 2 possible measurements of the seeing : the PSF model based version and
+        // the Moments based value we need to define a definitive "CHIP.SEEING" value.  This
+        // value is used by the PSF-matching programs (ppSub and ppStack).  The moments-based
+        // value is probably more representative of the value needed by those tools (it is also
+        // more stable?)
+        bool status = false;
+        float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+        float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+
+        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 * (fwhmMajorMoments + fwhmMinorMoments);
+
+        psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments);
+    }
+
+    if (psfExtra1->n) {
+        if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+            goto escape;
+        }
+        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     }
 
     if (psfExtra2->n) {
-      if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
-        goto escape;
-      }
-      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
-      psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+        if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
+            goto escape;
+        }
+        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     }
 

trunk/psphot/src/psphotReadoutCleanup.c

@@ -57 +57 @@
     // use the psf-model to measure FWHM stats
     if (psf) {
-        if (!psphotPSFstats (readout, psf)) {
+      if (!psphotPSFstatsSources (readout, sources, psf)) {
             psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
             return false;