Changeset 27840 for branches/simtest_nebulous_branches/psphot/src/psphotReadout.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotReadout.c (modified) (6 diffs)

branches/simtest_nebulous_branches/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/simtest_nebulous_branches/psphot/src/psphotReadout.c

@@ -3 +3 @@
 // this should be called by every program that links against libpsphot
 bool psphotInit (void) {
-
     psphotErrorRegister();              // register our error codes/messages
     psphotModelClassInit ();            // load implementation-specific models
@@ -12 +11 @@
 bool psphotReadout(pmConfig *config, const pmFPAview *view) {
 
-    // measure the total elapsed time in psphotReadout.  XXX the current threading plan
-    // for psphot envisions threading within psphotReadout, not multiple threads calling
-    // the same psphotReadout.  In the current plan, this dtime is the elapsed time used
-    // jointly by the multiple threads, not the total time used by all threads.
+    // 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");
+
+    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
 
     // select the current recipe
@@ -24 +23 @@
         return false;
     }
+    // optional break-point for processing
+    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
+    psAssert (breakPt, "configuration error: set BREAK_POINT");
 
     // set the photcode for this image
-    if (!psphotAddPhotcode (recipe, config, view, "PSPHOT.INPUT")) {
+    if (!psphotAddPhotcode (config, view)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
     }
 
-    // find the currently selected readout
-    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
-    PS_ASSERT_PTR_NON_NULL (readout, false);
-
-    // optional break-point for processing
-    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
-    PS_ASSERT_PTR_NON_NULL (breakPt, false);
-
     // Generate the mask and weight images, including the user-defined analysis region of interest
-    psphotSetMaskAndVariance (config, readout, recipe);
+    if (!psphotSetMaskAndVariance (config, view)) {
+        return psphotReadoutCleanup(config, view);
+    }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
-    }
-
-    // display the image, weight, mask (ch 1,2,3)
-    psphotVisualShowImage (readout);
+        return psphotReadoutCleanup(config, view);
+    }
 
     // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
-    }
-    if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
+    if (!psphotModelBackground (config, view)) {
+        return psphotReadoutCleanup (config, view);
+    }
+    if (!psphotSubtractBackground (config, view)) {
+        return psphotReadoutCleanup (config, view);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
-    }
-
-    // display the backsub and backgnd images
-    psphotVisualShowBackground (config, view, readout);
-    
-    // run a single-model test if desired (exits from here if test is run)
-    psphotModelTest (config, view, recipe);
-
-    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
-    pmPSF *psf = psphotLoadPSF (config, view, recipe);
-
-    // several functions below behave differently if we have a PSF model already
-    bool havePSF = (psf != NULL);
+        return psphotReadoutCleanup (config, view);
+    }
+
+    // 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)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
+        return psphotReadoutCleanup (config, view);
+    }
 
     // find the detections (by peak and/or footprint) in the image.
-    pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
-    if (!detections) {
-        psLogMsg ("psphot", 3, "unable to find detections in this image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
-    }
-
-    // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
-    if (!sources) return false;
+    if (!psphotFindDetections (config, view, 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);
+    }
+
+    // construct sources and measure basic stats (saved on detections->newSources)
+    if (!psphotSourceStats (config, view, true)) { // pass 1
+        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
+        return psphotReadoutCleanup (config, view);
+    }
     if (!strcasecmp (breakPt, "PEAKS")) {
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
-    }
-
-    // find blended neighbors of very saturated stars
-    // XXX merge this with Basic Deblend?
-    psphotDeblendSatstars (readout, sources, recipe);
-
-    // mark blended peaks PS_SOURCE_BLEND
-    if (!psphotBasicDeblend (sources, recipe)) {
-        psLogMsg ("psphot", 3, "failed on deblend analysis");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-    }
-
-    // classify sources based on moments, brightness
-    if (!psphotRoughClass (readout, sources, recipe, havePSF)) {
-        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
+        return psphotReadoutCleanup(config, view);
+    }
+
+    // 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);
+    }
+
+    // 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);
+    }
+
+    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
+    // clump defined for it is used not measured (detections->newSources)
+    if (!psphotRoughClass (config, view)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
+        return psphotReadoutCleanup (config, view);
+    }
+    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
+    if (!psphotImageQuality (config, view)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
+        return psphotReadoutCleanup(config, view);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
-    }
-
-    if (!havePSF && !psphotImageQuality (recipe, sources)) {
-        psLogMsg("psphot", 3, "failed to measure image quality");
-        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
-    }
-
-    // if we were not supplied a PSF, choose one here
-    if (psf == NULL) {
-        // use bright stellar objects to measure PSF
-        // XXX if we do not have enough stars to generate the PSF, build one
-        // from the SEEING guess and model class
-        psf = psphotChoosePSF (readout, sources, recipe);
-        if (psf == NULL) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model");
-            return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-        }
-        havePSF = true;
+        return psphotReadoutCleanup(config, view);
+    }
+
+    // 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)) { // pass 1
+        psLogMsg ("psphot", 3, "failure to construct a psf model");
+        return psphotReadoutCleanup (config, view);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-    }
-
-    psphotVisualShowPSFModel (readout, psf);
+        return psphotReadoutCleanup (config, view);
+    }
 
     // include externally-supplied sources
-    psphotLoadExtSources (config, view, sources);
-
-    // construct an initial model for each object
-    psphotGuessModels (config, readout, sources, psf);
-
-    // XXX test output of models
-    // psphotTestSourceOutput (readout, sources, recipe, psf);
-
-    // linear PSF fit to source peaks
-    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
-
-    // We have to place these visualizations here because the models are not realized until
-    // psphotGuessModels or fitted until psphotFitSourcesLinear.
-    psphotVisualShowPSFStars (recipe, psf, sources);
-    psphotVisualShowSatStars (recipe, psf, sources);
-
-    // identify CRs and extended sources
-    psphotSourceSize (config, readout, sources, recipe, 0);
+    // XXX fix this in the new multi-input context
+    // psphotLoadExtSources (config, view); // pass 1
+
+    // construct an initial model for each object, set the radius to fitRadius, set circular
+    // fit mask (detections->newSources)
+    psphotGuessModels (config, view); // pass 1
+
+    // merge the newly selected sources into the existing list
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view);
+
+    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
+    psphotFitSourcesLinear (config, view, false); // pass 1 (detections->allSources)
+
+    // identify CRs and extended sources (only unmeasured sources are measured)
+    psphotSourceSize (config, view, true); // pass 1 (detections->allSources)
     if (!strcasecmp (breakPt, "ENSEMBLE")) {
         goto finish;
@@ -151 +134 @@
 
     // non-linear PSF and EXT fit to brighter sources
-    psphotBlendFit (config, readout, sources, psf);
+    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
+    psphotBlendFit (config, view); // pass 1 (detections->allSources)
 
     // replace all sources
-    psphotReplaceAllSources (sources, recipe);
-
-    // linear fit to include all sources
-    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
+    psphotReplaceAllSources (config, view); // pass 1 (detections->allSources)
+
+    // linear fit to include all sources (subtract again)
+    // NOTE : apply to ALL sources (extended + psf)
+    psphotFitSourcesLinear (config, view, true); // pass 2 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
-    if (!strcasecmp (breakPt, "PASS1")) {
-        goto pass1finish;
-    }
-
-    // XXX for the moment, drop the re-calc of the background (prove this works)
-    // replace background in residual image
-    // psphotSkyReplace (config, view);
-    // re-measure background model (median, smoothed image)
-    // psphotImageMedian (config, view);
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
+    // NOTE: possibly re-measure background model here with objects subtracted / or masked
 
     // add noise for subtracted objects
-    psphotAddNoise (readout, sources, recipe);
-
-    // find fainter sources (pass 2)
-    detections = psphotFindDetections (detections, readout, recipe);
+    psphotAddNoise (config, view); // 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, false); // pass 2 (detections->peaks, detections->footprints)
 
     // remove noise for subtracted objects (ie, return to normal noise level)
-    psphotSubNoise (readout, sources, recipe);
+    // NOTE: this needs to operate only on the OLD sources
+    psphotSubNoise (config, view); // pass 1 (detections->allSources)
 
     // define new sources based on only the new peaks
-    psArray *newSources = psphotSourceStats (config, readout, detections);
+    // NOTE: new sources are saved on detections->newSources
+    psphotSourceStats (config, view, false); // pass 2 (detections->newSources)
 
     // set source type
-    if (!psphotRoughClass (readout, newSources, recipe, havePSF)) {
+    // NOTE: apply only to detections->newSources
+    if (!psphotRoughClass (config, view)) { // pass 2 (detections->newSources)
         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
-    }
-
-    // create full input models
-    psphotGuessModels (config, readout, newSources, psf);
+        return psphotReadoutCleanup (config, view);
+    }
+
+    // create full input models, set the radius to fitRadius, set circular fit mask
+    // NOTE: apply only to detections->newSources
+    psphotGuessModels (config, view); // pass 2 (detections->newSources)
 
     // replace all sources so fit below applies to all at once
-    psphotReplaceAllSources (sources, recipe);
+    // NOTE: apply only to OLD sources (which have been subtracted)
+    psphotReplaceAllSources (config, view); // pass 2
 
     // merge the newly selected sources into the existing list
-    psphotMergeSources (sources, newSources);
-    psFree (newSources);
-
-    // linear fit to all sources
-    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
+    // NOTE: merge OLD and NEW
+    // XXX check on free of sources...
+    psphotMergeSources (config, view); // (detections->newSources + detections->allSources -> detections->allSources)
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, true); // pass 3 (detections->allSources)
 
 pass1finish:
 
     // measure source size for the remaining sources
-    psphotSourceSize (config, readout, sources, recipe, 0);
-
-    psphotExtendedSourceAnalysis (readout, sources, recipe);
-
-    psphotExtendedSourceFits (readout, sources, recipe);
+    // NOTE: applies only to NEW (unmeasured) sources
+    psphotSourceSize (config, view, false); // pass 2 (detections->allSources)
+
+    psphotExtendedSourceAnalysis (config, view); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view); // pass 1 (detections->allSources)
 
 finish:
@@ -216 +202 @@
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, readout, sources, psf)) {
+    if (!psphotApResid (config, view)) { // pass 1 (detections->allSources)
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
+        return psphotReadoutCleanup (config, view);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, readout, view, sources, psf);
+    psphotMagnitudes(config, view); // pass 1 (detections->allSources)
+
+    if (!psphotEfficiency(config, view)) { // pass 1
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
 
     // replace failed sources?
@@ -228 +219 @@
 
     // replace background in residual image
-    psphotSkyReplace (config, view);
+    psphotSkyReplace (config, view); // pass 1
 
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (sources);
+    psphotSourceFreePixels (config, view); // pass 1
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
+    return psphotReadoutCleanup(config, view);
 }
-

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack