Index: /branches/eam_branches/ipp-20110213/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/Makefile.am (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/Makefile.am (revision 30749) @@ -129,4 +129,5 @@ psphotReadoutMinimal.c \ psphotModelBackground.c \ + psphotMaskBackground.c \ psphotSubtractBackground.c \ psphotFindDetections.c \ Index: /branches/eam_branches/ipp-20110213/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30749) @@ -12,4 +12,6 @@ #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs + +# define READOUT_OR_INTERNAL(VIEW,FILE)((FILE)->mode == PM_FPA_MODE_INTERNAL) ? (FILE)->readout : pmFPAviewThisReadout((VIEW), (FILE)->fpa) // top-level psphot functions @@ -52,4 +54,7 @@ bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotMaskBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotMaskBackground.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotMaskBackground.c (revision 30749) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotMaskBackground.c (revision 30749) @@ -0,0 +1,102 @@ +# include "psphotInternal.h" +# define NSIGMA 2.0 + +// mask pixel in the input image above model + N*stdev (using mark) +bool psphotMaskBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) +{ + bool status = true; + + psTimerStart ("psphot.background"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + psAssert (readout, "missing readout?"); + + // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename' + pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest + assert (modelFile); + + pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL.STDEV", index); + assert (stdevFile); + + float skyMean = psMetadataLookupF32(&status, readout->analysis, "SKY_MEAN"); + float skyStdv = psMetadataLookupF32(&status, readout->analysis, "SKY_STDEV"); + + pmReadout *model = READOUT_OR_INTERNAL(view, modelFile); + psAssert (model, "this must exist"); + + pmReadout *stdev = READOUT_OR_INTERNAL(view, stdevFile); + psAssert (stdev, "this must exist"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (markVal); + + psImageBinning *binning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING"); + assert (binning); + + // create a binned version of the threshold image (= model + NSIGMA*stdev) + psImage *threshBinned = psImageCopy (NULL, model->image, PS_TYPE_F32); + for (int iy = 0; iy < model->image->numRows; iy++) { + for (int ix = 0; ix < model->image->numCols; ix++) { + // threshBinned->data.F32[iy][ix] = model->image->data.F32[iy][ix] + NSIGMA*stdev->image->data.F32[iy][ix]; + threshBinned->data.F32[iy][ix] = skyMean + NSIGMA*skyStdv; + } + } + psphotSaveImage (NULL, threshBinned, "threshbin.fits"); + + // create a threshold image + psImage *threshold = psImageCopy (NULL, readout->image, PS_TYPE_F32); + + // linear interpolation to full-scale + if (!psImageUnbin (threshold, threshBinned, binning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + return false; + } + psphotSaveImage (NULL, threshold, "threshold.fits"); + + psLogMsg ("psphot", PS_LOG_MINUTIA, "build threshold image: %f sec\n", psTimerMark ("psphot.background")); + + // raise the 'markVal' mask bit for pixels above threshold + for (int iy = 0; iy < readout->image->numRows; iy++) { + for (int ix = 0; ix < readout->image->numCols; ix++) { + if (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + if (readout->image->data.F32[iy][ix] < threshold->data.F32[iy][ix]) continue; + readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= markVal; + } + } + psphotSaveImage (NULL, readout->mask, "newmask.fits"); + psLogMsg ("psphot", PS_LOG_INFO, "masked image based on 1st pass background model: %f sec\n", psTimerMark ("psphot.background")); + + psFree (threshold); + psFree (threshBinned); + + return true; +} + +bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotMaskBackgroundReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for %s entry %d", filerule, i); + return false; + } + } + return true; +} Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotModelBackground.c (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotModelBackground.c (revision 30749) @@ -58,4 +58,11 @@ psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (markVal); + + // apply both MASK and MARK to make background model + maskVal |= markVal; psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); @@ -308,4 +315,12 @@ } + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "backfill.%02d.fits", npass); + psphotSaveImage (NULL, model, name); + sprintf (name, "backfist.%02d.fits", npass); + psphotSaveImage (NULL, modelStdev, name); + } + psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background")); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotOutput.c (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotOutput.c (revision 30749) @@ -29,33 +29,24 @@ } -pmReadout *psphotSelectBackground (pmConfig *config, - const pmFPAview *view) { +pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view) { bool status; - pmReadout *background; + pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL"); if (!file) return NULL; - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + + pmReadout *background = READOUT_OR_INTERNAL(view, file); return background; } -pmReadout *psphotSelectBackgroundStdev (pmConfig *config, - const pmFPAview *view) { +pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view) { bool status; - pmReadout *background; pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV"); if (!file) return NULL; - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + + pmReadout *background = READOUT_OR_INTERNAL(view, file); return background; } Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30749) @@ -52,4 +52,17 @@ } +# if (0) + // XXX test to mask outliers, not very helpful + // mask the high values in the image (with MARK) + if (!psphotMaskBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } +# endif + if (!psphotSubtractBackground (config, view, filerule)) { return psphotReadoutCleanup (config, view, filerule); @@ -157,41 +170,85 @@ // NOTE: possibly re-measure background model here with objects subtracted / or masked - // add noise for subtracted objects - psphotAddNoise (config, view, filerule); // 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, filerule, 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, filerule); // pass 1 (detections->allSources) - - // define new sources based on only the new peaks - // NOTE: new sources are saved on detections->newSources - psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) - - // set source type - // NOTE: apply only to detections->newSources - if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, view, filerule); - } - - // create full input models, set the radius to fitRadius, set circular fit mask - // NOTE: apply only to detections->newSources - psphotGuessModels (config, view, filerule); // 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, filerule); // 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, filerule); // (detections->newSources + detections->allSources -> detections->allSources) - - // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + // NOTE: this block performs the 2nd pass low-significance PSF detection stage + { + // add noise for subtracted objects + psphotAddNoise (config, view, filerule); // 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, filerule, 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, filerule); // pass 1 (detections->allSources) + + // define new sources based on only the new peaks + // NOTE: new sources are saved on detections->newSources + psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) + + // set source type + // NOTE: apply only to detections->newSources + if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view, filerule); + } + + // create full input models, set the radius to fitRadius, set circular fit mask + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view, filerule); // 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, filerule); // 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, filerule); // (detections->newSources + detections->allSources -> detections->allSources) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + } + + // NOTE: this block performs the 2nd pass low-significance EXT detection stage (smooth or rebin by NxN times PSF size) + if (0) { + // add noise for subtracted objects + psphotAddNoise (config, view, filerule); // 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, filerule, 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, filerule); // pass 1 (detections->allSources) + + // define new sources based on only the new peaks + // NOTE: new sources are saved on detections->newSources + psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) + + // set source type + // NOTE: apply only to detections->newSources + if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view, filerule); + } + + // create full input models, set the radius to fitRadius, set circular fit mask + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view, filerule); // 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, filerule); // 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, filerule); // (detections->newSources + detections->allSources -> detections->allSources) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + } pass1finish: Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotSubtractBackground.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotSubtractBackground.c (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotSubtractBackground.c (revision 30749) @@ -24,10 +24,5 @@ assert (modelFile); - pmReadout *model = NULL; - if (modelFile->mode == PM_FPA_MODE_INTERNAL) { - model = modelFile->readout; - } else { - model = pmFPAviewThisReadout (view, modelFile->fpa); - } + pmReadout *model = READOUT_OR_INTERNAL(view, modelFile); assert (model); @@ -44,9 +39,5 @@ if (file) { // we are using PSPHOT.BACKGND as an I/O file: select readout or create - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + background = READOUT_OR_INTERNAL(view, file); if (background == NULL) { // readout does not yet exist: create from input Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotVisual.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotVisual.c (revision 30748) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotVisual.c (revision 30749) @@ -293,6 +293,4 @@ bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) { - pmReadout *backgnd; - if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true; @@ -303,9 +301,5 @@ pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); - if (file->mode == PM_FPA_MODE_INTERNAL) { - backgnd = file->readout; - } else { - backgnd = pmFPAviewThisReadout (view, file->fpa); - } + pmReadout *backgnd = READOUT_OR_INTERNAL(view, file); psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2);