Index: branches/eam_branches/ipp-20111110/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphot.h (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphot.h (revision 32685) @@ -108,6 +108,6 @@ bool psphotBlendFit_Threaded (psThreadJob *job); -bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState); +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool ignoreState); bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -176,5 +176,7 @@ // in psphotReplaceUnfit.c: -bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe); +bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe); +bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState); +bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState); bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); Index: branches/eam_branches/ipp-20111110/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotEfficiency.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotEfficiency.c (revision 32685) @@ -255,5 +255,5 @@ // remove all sources, adding noise for subtracted sources - psphotRemoveAllSources(realSources, recipe); + psphotRemoveAllSourcesByArray(realSources, recipe); #if TESTING Index: branches/eam_branches/ipp-20111110/psphot/src/psphotFake.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotFake.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotFake.c (revision 32685) @@ -166,5 +166,5 @@ // remove all sources, adding noise for subtracted sources - psphotRemoveAllSources(realSources, recipe); + psphotRemoveAllSourcesByArray(realSources, recipe); psphotAddNoise(readout, realSources, recipe); Index: branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c (revision 32685) @@ -407,4 +407,69 @@ return true; } + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + static int npass = 0; + char filename[64]; + + // XXX uncomment to disreturn true; + + bool status = true; + + int num = psphotFileruleCount(config, filerule); + + snprintf (filename, 64, "testdump.%02d.dat", npass); + FILE *f = fopen (filename, "w"); + + // 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->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + if (detections->newSources) { + fprintf (f, "## --- from new sources ---\n"); + } else { + fprintf (f, "## --- from all sources ---\n"); + } + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + + pmPeak *peak = source->peak; + if (!peak) continue; + + // XXX only dump a given region + if (peak->xf < 20) continue; + if (peak->yf < 20) continue; + if (peak->xf > 40) continue; + if (peak->yf > 70) continue; + + float Msum = source->moments ? source->moments->Sum : NAN; + float Mx = source->moments ? source->moments->Mx : NAN; + float My = source->moments ? source->moments->My : NAN; + float Npix = source->moments ? source->moments->nPixels : NAN; + fprintf (f, "%f %f : %f %f : %f %f\n", peak->xf, peak->yf, Mx, My, Msum, Npix); + } + } + fclose (f); + npass ++; + + return true; +} + +# if (0) bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule, char *filename) { @@ -449,3 +514,3 @@ return true; } - +# endif Index: branches/eam_branches/ipp-20111110/psphot/src/psphotRadialApertures.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotRadialApertures.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotRadialApertures.c (revision 32685) @@ -156,5 +156,5 @@ // set this to 0 to run without threading -# if (1) +# if (0) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); Index: branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c (revision 32685) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule); // this should be called by every program that links against libpsphot @@ -9,55 +10,4 @@ } -// for now, let's store the detections on the readout->analysis for each readout -bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule) -{ - static int npass = 0; - char filename[64]; - - return true; - - bool status = true; - - int num = psphotFileruleCount(config, filerule); - - snprintf (filename, 64, "chisq.%02d.dat", npass); - FILE *f = fopen (filename, "w"); - - // 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?"); - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) continue; - - pmModel *model = pmSourceGetModel (NULL, source); - if (!model) continue; - - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - fprintf (f, "%f %f %f %d %d %f 1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm); - } else { - fprintf (f, "%f %f %f %d %d %f 0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm); - } - } - } - fclose (f); - npass ++; - - return true; -} - bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) { @@ -133,4 +83,5 @@ return psphotReadoutCleanup (config, view, filerule); } + psphotDumpTest (config, view, filerule); // find blended neighbors of very saturated stars (detections->newSources) @@ -188,5 +139,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) - psphotDumpChisqs (config, view, filerule); + psphotDumpTest (config, view, filerule); // measure the radial profiles to the sky @@ -208,13 +159,11 @@ // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources) - psphotDumpChisqs (config, view, filerule); // replace all sources - psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources) + psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources) // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources) - psphotDumpChisqs (config, view, filerule); // if we only do one pass, skip to extended source analysis @@ -239,4 +188,5 @@ // NOTE: new sources are saved on detections->newSources psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) + psphotDumpTest (config, view, filerule); // set source type @@ -253,5 +203,5 @@ // 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 + psphotReplaceAllSources (config, view, filerule, false); // pass 2 // merge the newly selected sources into the existing list @@ -259,8 +209,8 @@ // XXX check on free of sources... psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) + psphotDumpTest (config, view, filerule); // NOTE: apply to ALL sources psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) - psphotDumpChisqs (config, view, filerule); } @@ -295,5 +245,5 @@ // 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 + psphotReplaceAllSources (config, view, filerule, false); // pass 2 // merge the newly selected sources into the existing list Index: branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c (revision 32685) @@ -1,3 +1,6 @@ # include "psphotInternal.h" + +static int replace_pass = 0; +static int remove_pass = 0; // replace the flux for sources which failed @@ -9,9 +12,9 @@ for (int i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // replace other sources? - if (source->mode & PM_SOURCE_MODE_FAIL) goto replace; - continue; + source = sources->data[i]; + + // replace other sources? + if (source->mode & PM_SOURCE_MODE_FAIL) goto replace; + continue; replace: @@ -23,5 +26,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule) +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState) { bool status = true; @@ -35,16 +38,16 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) { + if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) { psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); return false; } } - return true; -} - -bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { + replace_pass ++; + return true; +} + +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) { bool status; - pmSource *source; psTimerStart ("psphot.replace"); @@ -57,4 +60,8 @@ psAssert (readout, "missing readout?"); + char name[128]; + snprintf (name, 128, "testadd.neg.%02d.%02d.fits", replace_pass, index); + psphotSaveImage (NULL, readout->image, name); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -75,11 +82,28 @@ for (int i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // replace other sources? - if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; - - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } + pmSource *source = sources->data[i]; + + if (ignoreState) { + // rely on the type of source to decide if we subtract it or not + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // do not include CRs in the full ensemble fit + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + + // do not include MOMENTS_FAILURES in the fit + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + } else { + // if we respect the state, only remove unsubtracted sources + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + } + + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + snprintf (name, 128, "testadd.pos.%02d.%02d.fits", replace_pass, index); + psphotSaveImage (NULL, readout->image, name); psphotVisualShowImage(readout); @@ -88,5 +112,88 @@ } -bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState) +{ + bool status = true; + + // 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 (!psphotRemoveAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); + return false; + } + } + remove_pass ++; + return true; +} + +bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) { + + bool status; + + psTimerStart ("psphot.replace"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + char name[128]; + snprintf (name, 128, "testsub.pos.%02d.%02d.fits", remove_pass, index); + psphotSaveImage (NULL, readout->image, name); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // 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 + psAssert (maskVal, "missing mask value?"); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // do not include CRs in the full ensemble fit + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + + // do not include MOMENTS_FAILURES in the fit + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // if we respect the state, only remove unsubtracted sources + if (!ignoreState && (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + snprintf (name, 128, "testsub.neg.%02d.%02d.fits", remove_pass, index); + psphotSaveImage (NULL, readout->image, name); + + psphotVisualShowImage(readout); + psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); + return true; +} + +bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) { bool status; @@ -100,10 +207,10 @@ for (int i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // replace other sources? - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source = sources->data[i]; + + // replace other sources? + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); @@ -158,17 +265,17 @@ for (int i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // sources have not yet been subtracted in this image (but this flag may be raised) - source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; - if (!source->modelPSF) continue; - - float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - float radius = source->modelPSF->fitRadius; - - // force a redefine to this image - pmSourceFreePixels(source); - pmSourceRedefinePixels (source, readout, Xo, Yo, radius); + source = sources->data[i]; + + // sources have not yet been subtracted in this image (but this flag may be raised) + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + if (!source->modelPSF) continue; + + float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + float radius = source->modelPSF->fitRadius; + + // force a redefine to this image + pmSourceFreePixels(source); + pmSourceRedefinePixels (source, readout, Xo, Yo, radius); } return true; @@ -229,62 +336,62 @@ for (int i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // *** we need to cache the 'best' model, and we have 3 cases: - // 1) model is the psf model --> generate from the new psf - // 2) model is an unconvolved extended model --> just cache the copy (not perfect) - // 3) model is a convolved extended model --> re-generate - - // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits - bool isPSF = false; - pmModel *model = pmSourceGetModel(&isPSF, source); - if (!model) continue; - - float radius = model->fitRadius; // save for future use below - - // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM - if (isPSF || model->isPCM) { - // the guess central intensity comes from the peak: - float Io = source->peak->rawFlux; - float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - - // generate a model for this object with Io = 1.0 - pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); - if (modelPSF == NULL) { - psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo); - continue; - } - - // set the source PSF model - psFree (source->modelPSF); - source->modelPSF = modelPSF; - source->modelPSF->fitRadius = radius; - } - - if (model->isPCM) { - psAssert(false, "this section is not complete"); - - pmSourceCachePSF (source, maskVal); - - psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize); - if (!psfKernel) { - psWarning ("no psf kernel"); - } - - // generate an image of the right size - psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - psImageInit (rawModelFlux, 0.0); + source = sources->data[i]; + + // *** we need to cache the 'best' model, and we have 3 cases: + // 1) model is the psf model --> generate from the new psf + // 2) model is an unconvolved extended model --> just cache the copy (not perfect) + // 3) model is a convolved extended model --> re-generate + + // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits + bool isPSF = false; + pmModel *model = pmSourceGetModel(&isPSF, source); + if (!model) continue; + + float radius = model->fitRadius; // save for future use below + + // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM + if (isPSF || model->isPCM) { + // the guess central intensity comes from the peak: + float Io = source->peak->rawFlux; + float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // generate a model for this object with Io = 1.0 + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); + if (modelPSF == NULL) { + psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo); + continue; + } + + // set the source PSF model + psFree (source->modelPSF); + source->modelPSF = modelPSF; + source->modelPSF->fitRadius = radius; + } + + if (model->isPCM) { + psAssert(false, "this section is not complete"); + + pmSourceCachePSF (source, maskVal); + + psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize); + if (!psfKernel) { + psWarning ("no psf kernel"); + } + + // generate an image of the right size + psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + psImageInit (rawModelFlux, 0.0); - // insert the model image normalized to 1.0 - pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); - - psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel); + // insert the model image normalized to 1.0 + pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + + psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel); - psFree (psfKernel); - psFree (rawModelFlux); - } - - pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) + psFree (psfKernel); + psFree (rawModelFlux); + } + + pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) } Index: branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c =================================================================== --- branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c (revision 32645) +++ branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c (revision 32685) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule); // we have 3 possible real filesets: @@ -130,4 +131,5 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } + psphotDumpTest (config, view, STACK_SRC); psMemDump("sourcestats"); @@ -170,4 +172,5 @@ psphotFitSourcesLinear (config, view, STACK_SRC, false); psphotStackVisualFilerule(config, view, STACK_SRC); + psphotDumpTest (config, view, STACK_SRC); // measure the radial profiles to the sky @@ -186,6 +189,6 @@ psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources) - // replace all sources - psphotReplaceAllSources (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 @@ -194,4 +197,5 @@ // 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) psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources) @@ -200,4 +204,8 @@ // NOTE: this block performs the 2nd pass low-significance PSF detection stage { + // 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, true); // ignore subtraction state for sources + // add noise for subtracted objects psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources) @@ -212,6 +220,6 @@ // 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)) { - // 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"); @@ -220,7 +228,11 @@ } + // replace all sources in DET + psphotReplaceAllSources (config, view, STACK_DET, true); // ignore subtraction state for sources + // 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) + psphotDumpTest (config, view, STACK_SRC); // set source type @@ -237,5 +249,5 @@ // 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 + psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2 // merge the newly selected sources into the existing list @@ -243,4 +255,5 @@ // XXX check on free of sources... psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources) + psphotDumpTest (config, view, STACK_SRC); } @@ -256,4 +269,5 @@ psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); + psphotDumpTest (config, view, STACK_SRC); // NOTE: apply to ALL sources @@ -295,4 +309,6 @@ int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); + psphotDumpTest (config, view, STACK_SRC); + for (int entry = 1; entry < nRadialEntries; entry++) { // NOTE: entry 0 is the unmatched image set @@ -312,5 +328,5 @@ // 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