Changeset 33415 for branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotStackReadout.c (modified) (10 diffs)

branches/meh_branches/ppstack_test/psphot/src

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

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

branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c

@@ -44 +44 @@
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
 
+    psArray *objects = NULL; // used below after 'pass1finish' label
+
+    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
+    // by the multiple threads, not the total time used by all threads.
     psTimerStart ("psphotReadout");
 
@@ -56 +60 @@
     // optional break-point for processing
     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
-    PS_ASSERT_PTR_NON_NULL (breakPt, false);
-
-    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    psAssert (breakPt, "configuration error: set BREAK_POINT");
+
+    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
+    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
     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 
+    char *STACK_DET = STACK_RAW;
 
     // set the photcode for each image
@@ -71 +74 @@
     }
 
-    // Generate the mask and weight images
-    // XXX this should be done before we perform the convolutions
+    // Generate the mask and weight images (if not supplied) and set mask bits
     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -81 +83 @@
 
     // generate a background model (median, smoothed image)
-    // 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_DET)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -89 +89 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+    if (!psphotSubtractBackground (config, view, STACK_SRC)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
     if (!strcasecmp (breakPt, "BACKMDL")) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
+    // also make the chisq detection image
     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
@@ -109 +113 @@
     }
 
-    // copy the detections from DET to SRC 
+    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
+    // is just a copy of the container pointer; the sources on both DET and SRC are the same
+    // memory objects
     if (strcmp(STACK_SRC, STACK_DET)) {
         if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
@@ -122 +128 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
-
-    if (!strcasecmp (breakPt, "TEST1")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
+    if (!strcasecmp (breakPt, "PEAKS")) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
-
-    // generate the objects (object unify the sources from the different images)
-    // 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);
-
-    psMemDump("matchsources");
-
-    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_SRC, false)) { // pass 1
-        psFree(objects);
-        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    psMemDump("sourcestats");
-
-    // find blended neighbors of very saturated stars (detections->newSources)
-    // if (!psphotDeblendSatstars (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_SRC);
-    // }
-
-    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
-    // if (!psphotBasicDeblend (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     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_SRC)) {
-        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
         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_SRC)) { // pass 1
-        psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        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
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // use bright stellar objects to measure PSF
     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_SRC);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        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_SRC);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
 
     // merge the newly selected sources into the existing list
@@ -200 +164 @@
     psphotMergeSources (config, view, STACK_SRC);
 
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, STACK_SRC);
+
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinearStack (config, objects, FALSE);
+    psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
 
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
+    // re-measure the kron mags with models subtracted.  this pass starts with a circular
+    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
+    // but iterates to an appropriately larger size
+    psphotKronIterate(config, view, STACK_SRC);
+
     // identify CRs and extended sources
-    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_SRC)) {
-        psFree (objects);
-        psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
+    psphotSourceSize (config, view, STACK_SRC, true);
+
+    // non-linear PSF and EXT fit to brighter sources
+    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
+    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+
+    // replace all sources (do NOT ignore subtraction state)
+    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
+    // linear fit to include all sources (subtract again)
+    // NOTE : apply to ALL sources (extended + psf)
+    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
+
+    // NOTE: possibly re-measure background model here with objects subtracted / or masked
+
+    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
+    { 
+        // if DET and SRC are different images, generate children sources for all sources in
+        // the SRC image.  This operation replaces the existing DETECTION container on DET
+        // which is currently a view to the one on SRC).  children sources go to
+        // det->allSources
+        if (strcmp(STACK_SRC, STACK_DET)) {
+            psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
+
+            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
+            //  will somewhat oversubtract the sources -- this is OK
+            psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+        }
+
+        // add noise for subtracted objects
+        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+
+        // find fainter sources
+        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
+        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+
+        // remove noise for subtracted objects (ie, return to normal noise level)
+        // NOTE: this needs to operate only on the OLD sources
+        psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+
+        // if DET and SRC are different images, copy the detections from DET to SRC 
+        // (this operation just ensures the metadata container has a view on SRC as well
+        if (strcmp(STACK_SRC, STACK_DET)) {
+            // replace all sources in DET
+            psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+
+            // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
+            if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
+                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+                return psphotReadoutCleanup (config, view, STACK_SRC);
+            }
+        }
+
+        // define new sources based on only the new peaks & measure moments
+        // NOTE: new sources are saved on detections->newSources
+        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+
+        // set source type
+        // NOTE: apply only to detections->newSources
+        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
+            return psphotReadoutCleanup (config, view, STACK_SRC);
+        }
+
+        // replace all sources so fit below applies to all at once
+        // NOTE: apply only to OLD sources (which have been subtracted)
+        psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+
+        // merge the newly selected sources into the existing list
+        // NOTE: merge OLD and NEW
+        // XXX check on free of sources...
+        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+
+        // Construct an initial model for each object, set the radius to fitRadius, set circular
+        // fit mask.  NOTE: only applied to sources without guess models
+        psphotGuessModels (config, view, STACK_SRC);
+    }
+
+pass1finish:
+
+    // generate the objects (objects unify the sources from the different images) NOTE: could
+    // this just match the detections for the chisq image, and not bother measuring the source
+    // stats in that case...?
+    objects = psphotMatchSources (config, view, STACK_SRC);
+    psMemDump("matchsources");
+
+    // Construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask.  NOTE: only applied to sources without guess models
+    psphotGuessModels (config, view, STACK_SRC);
+
+    psphotStackObjectsUnifyPosition (objects);
+
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
+
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
+    // re-measure the kron mags with models subtracted
+    // psphotKronMasked(config, view, STACK_SRC);
+    psphotKronIterate(config, view, STACK_SRC);
+
+    // measure source size for the remaining sources
+    // NOTE: applies only to NEW (unmeasured) sources
+    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
 
-    psphotStackObjectsUnifyPosition (objects);
-
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, 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_SRC); // pass 1 (detections->allSources)
-
-    // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
 
     // create source children for the OUT filerule (for radial aperture photometry) 
@@ -238 +310 @@
     }
 
+    // measure circular, radial apertures (objects sorted by S/N)
+    // this forces photometry on the undetected sources from other images
+    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");
 
-    bool smoothAgain = true;
-    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
+
+    for (int entry = 1; entry < nRadialEntries; entry++) {
+        // NOTE: entry 0 is the unmatched image set
 
         // re-measure the PSF for the smoothed image (using entries in 'allSources')
@@ -253 +330 @@
 
         // measure circular, radial apertures (objects sorted by S/N)
-        psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); 
+        // entry 0 == unmatched? pass entry + 1?
+        psphotRadialApertures (config, view, STACK_OUT, entry); 
 
         // replace the flux in the image so it is returned to its original state
-        psphotReplaceAllSources (config, view, STACK_OUT);
+        psphotReplaceAllSources (config, view, STACK_OUT, false);
 
         // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-        psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+        psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
         psMemDump("matched");
     }
 
-    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
+    // measure aperture photometry corrections
+    if (!psphotApResid (config, view, STACK_SRC)) {
+        psFree (objects);
+        psLogMsg ("psphot", 3, "failed on psphotApResid");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // calculate source magnitudes
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
+    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
+    if (false && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();

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