Changeset 32348 for trunk/psphot/src/psphotStackPSF.c

Timestamp:

Sep 6, 2011, 1:32:31 PM (15 years ago)

Author:

eugene

Message:

• clarified the output to show the stages and their timing more clearly
• psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
• psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
• for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
• extended source analysis (petrosians) is threaded
• iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
• the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
• dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
• fake sources inserted for the efficiency analysis are removed at the end of the analysis
• radial apertures is now threaded
• edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
• new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
• define the fit region and window region based on the 1st radial moments
• use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
• remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotStackPSF.c (modified) (3 diffs)

trunk/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

trunk/psphot/src/psphotStackPSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask)
-{
+// determine the 1st target PSF (either AUTO or defined by PSPHOT.STACK.TARGET.PSF.FWHM)
+bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options)  {
+
     bool mdok = false;
-    pmPSF *psf = NULL;
+
+    int numCols = options->numCols;
+    int numRows = options->numRows;
+    psArray *psfs = options->psfs;
+    psVector *inputMask = options->inputMask;
 
     // Get the recipe values
@@ -17 +22 @@
         
     char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
+
+    options->targetSeeing = psVectorAllocEmpty(4, PS_TYPE_F32);
 
     if (autoPSF) {
@@ -38 +45 @@
 
         // Solve for the target PSF
-        psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
-        if (!psf) {
+        options->psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
+        if (!options->psf) {
             psError(PSPHOT_ERR_PSF, false, "Unable to determine output PSF.");
-            return NULL;
+            return false;
         }
 
-    } else {
-
-        // externally-defined PSF
-        // XXX need to test for compatibility of target with inputs
-
-        float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
-        if (!mdok) {
-            psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
-            psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
-            targetFWHM = fwhmValues->data.F32[0];
-        }
-
-        // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
-        // GAUSS  : 1.000
-        // PGAUSS : 1.006
-        // QGAUSS : 1.151
-        // RGAUSS : 0.883
-        // PS1_V1 : 1.134
-        
-        float scaleFactor = NAN;
-        if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
-            scaleFactor = 1.000;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
-            scaleFactor = 1.0006;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
-            scaleFactor = 1.151;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
-            scaleFactor = 0.883;
-        }
-        if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
-            scaleFactor = 1.134;
-        }
-        psAssert (isfinite(scaleFactor), "invalid model for PSF"); 
-
-        float Sxx = sqrt(2.0)*targetFWHM / 2.35 / scaleFactor;
-
-        // XXX probably should make the model type (and par 7) optional from recipe
-        // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
-        psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
-        if (!psf) {
-            psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
-            return NULL;
-        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+        float targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
+        psVectorAppend(options->targetSeeing, targetSeeing);
+        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (auto-scaled): %f\n", targetSeeing);
+        return true;
     }
 
-    return psf;
+    // externally-defined PSF
+    // XXX need to test for compatibility of target with inputs
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+    if (!mdok) {
+        psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+        psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+        for (int i = 0; i < fwhmValues->n; i++) {
+            psVectorAppend(options->targetSeeing, fwhmValues->data.F32[i]);
+        }           
+        targetSeeing = fwhmValues->data.F32[0];
+    } else {
+        psVectorAppend(options->targetSeeing, targetSeeing);
+    }
+
+    // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0))
+    // GAUSS  : 1.000
+    // PGAUSS : 1.006
+    // QGAUSS : 1.151
+    // RGAUSS : 0.883
+    // PS1_V1 : 1.134
+        
+    float scaleFactor = NAN;
+    if (!strcmp(psfModel, "PS_MODEL_GAUSS")) {
+        scaleFactor = 1.000;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) {
+        scaleFactor = 1.0006;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) {
+        scaleFactor = 1.151;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) {
+        scaleFactor = 0.883;
+    }
+    if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) {
+        scaleFactor = 1.134;
+    }
+    psAssert (isfinite(scaleFactor), "invalid model for PSF"); 
+
+    float Sxx = sqrt(2.0)*targetSeeing / 2.35 / scaleFactor;
+
+    // XXX probably should make the model type (and par 7) optional from recipe
+    // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
+    options->psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2);
+    if (!options->psf) {
+        psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
+        return false;
+    }
+
+    psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+    psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM (1 of %ld): %f\n", options->targetSeeing->n, targetSeeing);
+    return true;
 }

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest