Changeset 30624 for trunk/psphot/src/psphotStackReadout.c

Timestamp:

Feb 13, 2011, 12:33:05 PM (15 years ago)

Author:

eugene

Message:

adding some of the metadata needed by PSPS to output headers; skip models with few valid pixels (eg, only the edge is showing); only do the linear fit on pixels within the fit radius; modify psf-match auto-scaling process; enable multiple target psfs for matched-psf aperture photometry; speed up analysis of the radial apertures

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src/psphotStackReadout.c (modified) (7 diffs)

trunk/psphot/src/psphotStackReadout.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// we have 3 possible real filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
+
+// we have 3 files on which we operate:
+// DET (detection image)       : nominally RAW (optionally CNV?)
+// SRC (source analysis image) : nominally CNV (optionally RAW)
+// OUT (psf-matched images)    : always OUT
+
+bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        psphotVisualShowResidualImage (readout, true);
+        psphotVisualShowObjectRegions (readout, recipe, sources);
+    }
+    return true;
+}
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
@@ -20 +58 @@
     PS_ASSERT_PTR_NON_NULL (breakPt, false);
 
+    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
+    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
+    char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
+
     // we have 3 relevant files: RAW, CNV, OUT 
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_OUT)) {
+    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -30 +73 @@
     // Generate the mask and weight images
     // XXX this should be done before we perform the convolutions
-    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -40 +83 @@
     // XXX I think this is not defined correctly for an array of images.
     // XXX probably need to subtract the model (same model?) for both RAW and OUT
-    if (!psphotModelBackground (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotModelBackground (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
-    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
-    if (!psphotLoadPSF (config, view, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "CHISQ")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
+    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
         // this only happens if we had an error in psphotFindDetections
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // copy the detections from RAW to OUT
-    if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // copy the detections from DET to SRC 
+    if (strcmp(STACK_SRC, STACK_DET)) {
+        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
+            psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+            return psphotReadoutCleanup (config, view, STACK_SRC);
+        }
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    // only run this on detections from the input images, not chisq image
-    if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    if (!strcasecmp (breakPt, "TEST1")) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // generate the objects (object unify the sources from the different images)
-    psArray *objects = psphotMatchSources (config, view, STACK_OUT);
+    // XXX this could just match the detections for the chisq image, and not bother measuring the
+    // source stats in that case...
+    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
+
+    if (!strcasecmp (breakPt, "TEST2")) {
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
 
     // construct sources for the newly-generated sources (from other images)
-    if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
+        psFree(objects);
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -98 +147 @@
     // if (!psphotDeblendSatstars (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
@@ -104 +153 @@
     // if (!psphotBasicDeblend (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_OUT)) {
+    if (!psphotRoughClass (config, view, STACK_SRC)) {
+        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     // this step is skipped
-    if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+        psFree(objects);
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        psFree(objects);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, view, STACK_OUT);
+    psphotGuessModels (config, view, STACK_SRC);
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_OUT);
+    psphotMergeSources (config, view, STACK_SRC);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinearStack (config, objects, FALSE);
+    psphotStackVisualFilerule(config, view, STACK_SRC);
 
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_OUT, TRUE);
+    psphotSourceSize (config, view, STACK_SRC, TRUE);
+
+    // XXX do we want to do a preliminary (unconvolved) model fit here, and then 
+    // do a second detection pass? (like standard psphot)
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_OUT)) {
+    if (!psphotApResid (config, view, STACK_SRC)) {
         psFree (objects);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     psphotStackObjectsUnifyPosition (objects);
 
-    // measure circular, radial apertures (objects sorted by S/N)
-    psphotRadialAperturesByObject (config, objects, view, STACK_OUT); 
-
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_OUT);
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    // create source children for the OUT filerule (for radial aperture photometry) 
+    psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    if (!objectsRadial) {
+        psFree(objects);
+        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    bool smoothAgain = true;
+    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+
+        // re-measure the PSF for the smoothed image (using entries in 'allSources')
+        psphotChoosePSF (config, view, STACK_OUT, false);
+
+        // this is necessary to update the models based on the new PSF
+        psphotResetModels (config, view, STACK_OUT);
+
+        // this is necessary to get the right normalization for the new models
+        psphotFitSourcesLinear (config, view, STACK_OUT, false);
+
+        // measure circular, radial apertures (objects sorted by S/N)
+        psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); 
+
+        // replace the flux in the image so it is returned to its original state
+        psphotReplaceAllSources (config, view, STACK_OUT);
+
+        // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
+        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+    }
 
     if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
@@ -176 +261 @@
 
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_RAW);
+    psphotSkyReplace (config, view, STACK_DET);
 
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (config, view, STACK_OUT);
-
-    // remove chisq image from config->file:PSPHOT.INPUT (why?)
-    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
+    // psphotSourceFreePixels (config, view, STACK_OUT);
+    psphotSourceFreePixels (config, view, STACK_SRC);
+
+    // remove chisq image from config->file:PSPHOT.INPUT
+    psphotStackRemoveChisqFromInputs(config, STACK_DET);
+    if (strcmp(STACK_SRC, STACK_DET)) {
+        psphotStackRemoveChisqFromInputs(config, STACK_SRC);
+    }
 
     psFree (objects);
+    psFree (objectsRadial);
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_OUT);
+    return psphotReadoutCleanup (config, view, STACK_SRC);
 }
 
+/* here is the process:
+
+ * we have three(*) images:
+ * RAW : unconvolved image stack 
+ * CNV : input convolved image
+
+ * OUT : psf-matched output image (there may be more than one of
+ * these.  we will generate the first matched image by selecting the
+ * target PSF and doing a full psf-maching process (as used by ppStack
+ * and ppSub).  But, additional target output files should use a
+ * simple gaussian convolution kernel determind from therms of the
+ * current and the target).
+
+ * the output should be / could be one of the matched images, but not
+ * all.  should we ensure the first gets written out, and ot save the
+ * others (or only optionally).
+
+ * by default, we probably only sve the cmf ffile outputs.
+
+ * load the RAW image (unconvolved stacks)
+ * add photcode to the output headers / readout->analysis
+ * generate mask and variance image (this is probably never needed in
+   practice: we always load an input mask & var.
+ * generate & subtract a model background for ?? (RAW? CNV? OUT? all?)
+ * load a PSF (probably not yet working)
+
+ * generate the CHISQ image from the RAW input images (why save on OUT?)
+
+ * find detections on RAW
+
+ * copy detections to OUT
+
+ * generate source stats (moments) for OUT
+
+ * match sources across inputs (on OUT?)
+
+ * generate source stats for the new constructions
+
+ * rough class (star, galaxy, cosmic, etc)
+
+ * Image quality
+
+ * generate PSF
+
+ * guess models
+
+ * merge sources (new -> old)
+
+ * linear fit to the psf
+
+ * find ApResid
+
+ * assign common positions
+
+ * radial apertures (** this should be on the PSF-matched images
+
+ * extended analysis (elliptical profile & petrosian)
+
+ * extended fits (sersic, etc)
+
+ * psphot magnitudes
+
+
+ ******
+
+ the above is all wrong:  first, we should be doing the full
+ morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or
+ RAW image (as desired optionally), etc.
+
+ In the discussion below, 'BST' (best) means optionally RAW or CNV
+
+ * detection : RAW & CHISQ (of RAW)
+ * moments : used by psf analysis & classification (BST)
+ * rough class : uses moments, not pixels
+ * image quality : uses moments as well
+ * generate PSF : (BST)
+ * guess models (BST)
+ * linear fit (BST)
+ * find ApResid (BST) -- uses sources not pixels
+ * extended analysis (BST)
+ * extended fits (BST)
+ * detection efficiency (BST)
+ 
+ * somehow need to copy the sources so they point at the pixels on the
+ * OUT image 
+
+ * foreach target PSF
+   * radial aperture
+   * convolve to next target PSF
+
+   * somehow need to organize the output file to have the values from
+   * the different PSFs in separate tables (with header info to
+   * specify the size of that PSF)
+
+   */