Changeset 20234

Timestamp:

Oct 17, 2008, 12:59:40 PM (18 years ago)

Author:

eugene

Message:

calculate FWHM for grid of points; add FWHM sigmas and ranges

File:

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

trunk/psphot/src/psphotChoosePSF.c

@@ -318 +318 @@
 bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
 
-    double FWHM_X, FWHM_Y;
-
     psEllipseShape shape;
     psEllipseAxes axes;
@@ -330 +328 @@
     PS_ASSERT_PTR_NON_NULL(image, false);
 
-    // XXXX A Serious hack: the psf might not be determined in the field center.
-    // for the moment, just try a few locations until it is defined!
-
-    pmModel *modelPSF = NULL;
-    for (int ix = -1; ix <= +1; ix ++) {
-        for (int iy = -1; iy <= +1; iy ++) {
+    // 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);
+
+    for (float ix = -0.4; ix <= +0.4; ix += 0.1) {
+        for (float iy = -0.4; iy <= +0.4; iy += 0.1) {
 
             // use the center of the center pixel of the image
-            float xc = (0.5 + 0.3*ix)*image->numCols + image->col0 + 0.5;
-            float yc = (0.5 + 0.3*ix)*image->numRows + image->row0 + 0.5;
+            float xc = ix*image->numCols + 0.5*image->numCols + image->col0 + 0.5;
+            float yc = iy*image->numRows + 0.5*image->numRows + image->row0 + 0.5;
 
             // create modelPSF from this model
-            modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
-            if (modelPSF) goto got_model;
+            pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
+            if (!modelPSF) continue;
+
+            // get the model full-width at half-max
+            float FWHM_MAJOR = 2*modelPSF->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);
+            psFree (modelPSF);
+
+            float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue;
+            psVectorAppend (fwhmMajor, FWHM_MAJOR);
+            psVectorAppend (fwhmMinor, FWHM_MINOR);
         }
     }
-    psAssert (modelPSF, "Failed to estimate PSF model at image center (psf is invalid everywhere?)");
-
-got_model:
-    // get the model full-width at half-max
-    FWHM_X = 2*modelPSF->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);
-
-    FWHM_Y = FWHM_X * (axes.minor / axes.major);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
+    psVectorStats (stats, fwhmMajor, NULL, NULL, 0);
+    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 (upper quartile)", stats->sampleLQ);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleUQ);
+
+    psVectorStats (stats, fwhmMinor, NULL, NULL, 0);
+    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 (upper quartile)", stats->sampleLQ);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower 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);
-    psFree (modelPSF);
-
+
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (stats);
     return true;
 }