Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 31154) +++ trunk/psphot/src/psphotStackReadout.c (revision 32348) @@ -44,4 +44,6 @@ bool psphotStackReadout (pmConfig *config, const pmFPAview *view) { + // 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,12 +58,11 @@ // 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,6 +72,5 @@ } - // 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); @@ -80,7 +80,9 @@ } + // XXX I think this is not defined correctly for an array of images. + // XXX I probably need to subtract the model (same model?) for both RAW and OUT. + // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level. + // 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); @@ -93,4 +95,5 @@ } + // 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,5 +112,5 @@ } - // copy the detections from DET to SRC + // if DET and SRC are different images, copy the detections from DET to SRC if (strcmp(STACK_SRC, STACK_DET)) { if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) { @@ -122,111 +125,150 @@ 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); + } psMemDump("sourcestats"); - // generate the objects (object unify the sources from the different images) + // classify sources based on moments, brightness + // only run this on detections from the input images, not chisq image + if (!psphotRoughClass (config, view, STACK_SRC)) { + 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 + psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + if (!strcasecmp (breakPt, "MOMENTS")) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } + + // use bright stellar objects to measure PSF + if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1 + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + if (!strcasecmp (breakPt, "PSFMODEL")) { + 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); + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (config, view, STACK_SRC); + + // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) + // XXX why do this as a stack operation? + // psphotFitSourcesLinearStack (config, objects, false); + psphotFitSourcesLinear (config, view, STACK_SRC, false); + psphotStackVisualFilerule(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); + + // 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 + psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources) + + // linear fit to include all sources (subtract again) + // NOTE : apply to ALL sources (extended + psf) + psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources) + + // if we only do one pass, skip to extended source analysis + if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; + + // NOTE: possibly re-measure background model here with objects subtracted / or masked + + // NOTE: this block performs the 2nd pass low-significance PSF detection stage + { + // 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 + if (strcmp(STACK_SRC, STACK_DET)) { + // XXX how does this handle 1st vs 2nd pass sources? + if (!psphotCopySources (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); + } + + // create full input models, set the radius to fitRadius, set circular fit mask + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view, STACK_SRC); // pass 2 (detections->newSources) + + // 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); // 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) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources) + } + +pass1finish: + + // 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"); + + // generate the objects (objects 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 - 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); - - // merge the newly selected sources into the existing list - // NOTE: merge OLD and NEW - 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_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); - } - - psMemDump("psfstats"); - psphotStackObjectsUnifyPosition (objects); + psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, 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,8 +280,15 @@ } + // measure circular, radial apertures (objects sorted by S/N) + // XXX can we just use psphotRadialApertures + // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none') + // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF); + 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,5 +302,6 @@ // 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 @@ -259,9 +309,19 @@ // 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); + + if (0 && !psphotEfficiency(config, view, STACK_DET)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear();