Changeset 10202

Timestamp:

Nov 26, 2006, 12:37:44 PM (20 years ago)

Author:

eugene

Message:

added option to include weights in 2D parameter fits (PSF_PARAM_WEIGHTS); measure the stats of dSXX,dSYY to better judge psf-like fits

File:

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

trunk/psphot/src/psphotChoosePSF.c

@@ -26 +26 @@
     }
 
+    float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
+    if (!status) {
+        PSF_MIN_DS = 0.01;
+        psWarning("PSF_MIN_DS is not set --- defaulting to %f\n", PSF_MIN_DS);
+    }
+
     // supply the measured sky variance for optional constant errors (non-poissonian)
     float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV");
@@ -39 +45 @@
     }
     pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_STDEV), POISSON_ERRORS);
+
+    // use poissonian errors or local-sky errors
+    bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS");
+    if (!status) {
+        PSF_PARAM_WEIGHTS = false;
+        psWarning("PSF_PARAM_WEIGHTS is not set in the recipe --- defaulting to false.\n");
+    }
 
     // how to model the PSF variations across the field
@@ -98 +111 @@
         psMetadataItem *item = psListGetAndIncrement (iter);
         modelName = item->data.V;
-        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask);
+        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS);
     }
 
@@ -132 +145 @@
     try = models->data[bestN];
 
+    // measure and save the median value of dparams[PM_PAR_SXX],dparams[PM_PAR_SYY]
+    // these are used by psphotEvalPSF to identify extended sources
+    psVector *Sx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *Sy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    for (int i = 0; i < try->sources->n; i++) {
+        // masked for: bad model fit, outlier in parameters
+        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
+            continue;
+
+        pmSource *source = try->sources->data[i];
+        Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
+        Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
+        dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
+        dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
+        dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
+        Sx->n ++;
+        Sy->n ++;
+        dSN->n ++;
+        dSx->n ++;
+        dSy->n ++;
+    }
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    psVectorStats (stats, dSx, NULL, NULL, 0);
+    float dSxo = stats->sampleMean;
+    psLogMsg ("psphot.choosePSF", 3, "dSx: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
+
+    psVectorStats (stats, dSy, NULL, NULL, 0);
+    float dSyo = stats->sampleMean;
+    psLogMsg ("psphot.choosePSF", 3, "dSy: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
+
+    for (int i = 0; i < dSN->n; i++) {
+        dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
+        dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
+    }
+
+    psVectorStats (stats, dSx, NULL, NULL, 0);
+    float nSx = stats->sampleStdev;
+    psLogMsg ("psphot.choosePSF", 3, "dSx: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
+    psVectorStats (stats, dSy, NULL, NULL, 0);
+    float nSy = stats->sampleStdev;
+    psLogMsg ("psphot.choosePSF", 3, "dSy: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
+
+    psFree (Sx);
+    psFree (Sy);
+    psFree (dSx);
+    psFree (dSy);
+    psFree (dSN);
+    psFree (stats);
+
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_MEAN", PS_META_REPLACE, "mean offset of dSXX", dSxo);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_MEAN", PS_META_REPLACE, "mean offset of dSYY", dSyo);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_STDV", PS_META_REPLACE, "stdev of mean-corrected dSXX", nSx);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_STDV", PS_META_REPLACE, "stdev of mean-corrected dSYY", nSy);
+
+    // unset the PSFSTAR flag for stars not used for PSF model
+    for (int i = 0; i < try->sources->n; i++) {
+        pmSource *source = try->sources->data[i];
+        if (try->mask->data.U8[i]) {
+            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
+        }
+    }
+
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) { 
@@ -159 +238 @@
     }
 
-    // unset the PSFSTAR flag for stars not used for PSF model
-    for (int i = 0; i < try->sources->n; i++) {
-        pmSource *source = try->sources->data[i];
-        if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) {
-            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
-        }
-    }
-
     // save only the best model;
     // XXX we are not saving the fitted sources