Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 29936) +++ trunk/psphot/src/psphotStackReadout.c (revision 30624) @@ -1,6 +1,44 @@ # 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,8 +58,13 @@ 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,9 +73,9 @@ // 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,57 +83,63 @@ // 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,5 +147,5 @@ // 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,66 +153,102 @@ // 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,16 +261,116 @@ // 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) + + */