Changeset 27876 for branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c

Timestamp:

May 6, 2010, 8:02:38 PM (16 years ago)

Author:

eugene

Message:

working on psphotStack

File:

• branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c (modified) (4 diffs)

branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass)
+bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view)
 {
     bool status = true;
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); 
+    psAssert(recipe, "We've thrown an error on this before.");
 
     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
 
-    // loop over the available readouts
+    // 'options' carries info needed to perform the stack matching
+    psphotStackOptions *options = psphotStackOptionsAlloc(num);
+
+    options->convolve = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.MATCH.PSF");
+    psAssert (status, "PSPHOT.STACK.MATCH.PSF not in recipe");
+
+    if (options->convolve) {
+        char *convolveSource = psMetadataLookupStr (&status, recipe, "PSPHOT.STACK.MATCH.PSF.SOURCE");
+        options->convolveSource = psphotStackConvolveSourceFromString (convolveSource);
+        if (options->convolveSource == PSPHOT_CNV_SRC_NONE) {
+            psError (PSPHOT_ERR_CONFIG, true, "stack convolution source not defined in recipe");
+            return false;
+        }
+    }
+
+    // loop over the available readouts (ignore chisq image)
     for (int i = 0; i < num; i++) {
-        if (!psphotStackMatchPSFsReadout (config, view, i)) {
+        if (!psphotStackMatchPSFsPrepare (config, view, options, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i);
+            return false;
+        }
+    }
+
+    // Generate target PSF
+    if (options->convolve) {
+        options->psf = ppStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask);
+        if (!options->psf) {
+            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
+            return false;
+        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf);
+        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target
+        psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
+
+        // XXX is this needed to supply the psf to psphot??
+        // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
+        // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf);
+        // outChip->data_exists = true;
+    }
+
+    // loop over the available readouts (ignore chisq image)
+    for (int i = 0; i < num; i++) { 
+        if (!psphotStackMatchPSFsReadout (config, view, options, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
             return false;
@@ -19 +62 @@
 
 // convolve the image to match desired PSF
-bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) {
-
-    bool status;
-    int pass;
-    float NSIGMA_PEAK = 25.0;
-    int NMAX = 0;
+bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) {
 
     // find the currently selected readout
-    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", index); // File of interest
     psAssert (fileRaw, "missing file?");
 
-    pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest
+    pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", index); // File of interest
     psAssert (fileCnv, "missing file?");
 
@@ -39 +77 @@
     psAssert (readoutCnv, "missing readout?");
 
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    psAssert (maskVal, "missing mask value?");
-
-    /***** set up recipe options *****/
-
-    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
-    psAssert(stackRecipe, "We've thrown an error on this before.");
-
-    // Look up appropriate values from the ppSub recipe
-    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
-    psAssert(subRecipe, "recipe missing");
-
-    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
-
-    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
-
-    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
-    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
-    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
-    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
-    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
-    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-
-    bool mdok;                          // Status of MD lookup
-    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
-    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
-
-    if (!pmReadoutMaskNonfinite(readout, maskVal)) {
+    // set NAN pixels to 'SAT'
+    psImageMaskType maskVal = pmConfigMaskGet("SAT", config);
+    if (!pmReadoutMaskNonfinite(readoutRaw, maskVal)) {
         psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
         return false;
@@ -74 +86 @@
     // Image Matching (PSFs or just flux)
     if (options->convolve) {
-      // Full match of PSFs
-        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
-
-        // Read previously produced kernel
-        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
-          loadKernel();
-        } else {
-          matchKernel();
-        } // !DEBUG.STACK
-
-        saveMatchData();
-
-        saveChiSquare();
-
-        renormKernel();
-
-        // Reject image completely if the maximum deconvolution fraction exceeds the limit
-        float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
-        if (deconv > deconvLimit) {
-            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
-            psFree(conv);
-            return NULL;
-        }
-
-        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
-
-        psFree(conv);
+        matchKernel(config, readoutCnv, readoutRaw, options, index);
+        saveMatchData(readoutCnv, options, index);
+        // renormKernel(readoutCnv, options, index);
     } else {
-        // only match the flux
-        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
-        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
-        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
+        // only match the flux (NO! not for multi-filter, at least!)
+        // XXX do not generate readoutCnv in this case?
+        // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
+        // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+        // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     }
 
-    // Ensure the background value is zero
-    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-        psWarning("Can't measure background for image.");
-        psErrorClear();
-    } else {
-        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-            (void)psBinaryOp(readout->image, readout->image, "-",
-                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-        }
-    }
+    rescaleData(readoutCnv, config, options, index);
 
-    if (!stackRenormaliseReadout(config, readout)) {
-        psFree(rng);
-        psFree(bg);
-        return false;
-    }
-
-    // Measure the variance level for the weighting
-    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
-        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
-            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
-            psFree(rng);
-            psFree(bg);
-            return false;
-        }
-        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
-    } else {
-        options->weightings->data.F32[index] = 1.0;
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
-             index, options->weightings->data.F32[index]);
-
-    psFree(rng);
-    psFree(bg);
-
-    dumpImage3();
+    dumpImage(readoutCnv, readoutRaw, index, "convolved");
 
     return true;
 }
+
+
+# if (0)
+// Read previously produced kernel
+if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
+    loadKernel(config, readoutCnv, options, index);
+} else {
+    matchKernel(config, readoutCnv, readoutRaw, options, index);
+}
+# endif