Index: /branches/eam_branches/ipp-20111110/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphot.h (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphot.h (revision 32689) @@ -455,4 +455,10 @@ psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc); +bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); +bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index); + +bool psphotCopyPeaks (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); +bool psphotCopyPeaksReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index); + bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source); void psphotSersicModelClassInit (); Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotBlendFit.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotBlendFit.c (revision 32689) @@ -236,13 +236,15 @@ pmSource *source = sources->data[i]; -# define TEST_X -420.0 -# define TEST_Y 300.0 +# if (0) +# define TEST_X 34 +# define TEST_Y 28 if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { - fprintf (stderr, "test galaxy\n"); + fprintf (stderr, "test object\n"); } # undef TEST_X # undef TEST_Y +# endif // skip non-astronomical objects (very likely defects) Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotFitSourcesLinear.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotFitSourcesLinear.c (revision 32689) @@ -342,9 +342,8 @@ psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); - // XXXX **** philosophical question: - // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit: - // should retain the chisq and errors from the intermediate non-linear fit? - // the non-linear fit provides better values for the position errors, and for - // extended sources, the shape errors + // Philosophical question: we measure bright objects in three passes: 1) linear fit; 2) + // non-linear fit; 3) linear fit: should we retain the chisq and errors from the + // intermediate non-linear fit? the non-linear fit provides better values for the position + // errors, and for extended sources, the shape errors // adjust I0 for fitSources and subtract Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotMergeSources.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotMergeSources.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotMergeSources.c (revision 32689) @@ -520,6 +520,6 @@ } -// create source parents children from ruleSrc for ruleOut for orphans -bool psphotSourceParent (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) +// copy the newPeaks from the detections of one pmFPAfile to another +bool psphotCopyPeaks (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) { bool status = true; @@ -534,5 +534,71 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotSourceParentReadout (config, view, ruleOut, ruleSrc, i)) { + if (!psphotCopyPeaksReadout (config, view, ruleOut, ruleSrc, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to copy peaks from %s to %s entry %d", ruleSrc, ruleOut, i); + return false; + } + } + return true; +} + +// add newly detected peaks to the existing list of sources +bool psphotCopyPeaksReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest + psAssert (fileSrc, "missing file?"); + + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + // find the currently selected readout + pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest + psAssert (fileOut, "missing file?"); + + pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa); + psAssert (readoutOut, "missing readout?"); + + // generate a new detection structure for the output filerule + pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS"); + psAssert (detectionsOut, "missing PSPHOT.DETECTIONS?"); + + psAssert (detectionsOut->peaks, "programming error"); + psAssert (!detectionsOut->oldPeaks, "programming error"); + psAssert (detections->peaks, "programming error"); + + // save the OUT existing peaks on oldPeaks + detectionsOut->oldPeaks = detectionsOut->peaks; + detectionsOut->peaks = psArrayAllocEmpty(detections->peaks->n); + + for (int i = 0; i < detections->peaks->n; i++) { + psAssert (detections->peaks->data[i], "programming error"); + pmPeak *peak = pmPeakAlloc (0, 0, 0.0, PM_PEAK_LONE); + pmPeakCopy(peak, detections->peaks->data[i]); + psArrayAdd (detectionsOut->peaks, 100, peak); + psFree (peak); + } + return true; +} + +// create source parents children from ruleSrc for ruleOut for orphans +bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) +{ + bool status = true; + + int num = psphotFileruleCount(config, ruleSrc); + + // skip the chisq image because it is a duplicate of the detection version + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotSourceParentsReadout (config, view, ruleOut, ruleSrc, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i); return false; @@ -542,14 +608,10 @@ } -// create source parents from ruleSrc for ruleOut for orphaned children for this readout. XXX currently, this is only -// used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called -// repeatedly) - -// XXX this function will do bad things if called repeatedly. It needs to either (1) delete -// the existing DETECTION container on the output images, or (2) intelligently supplement the -// existing DETECTION container with only the new detections. -bool psphotSourceParentReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { +// create source parents from ruleSrc for ruleOut for orphaned children for this readout. +bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { bool status; + int nParents = 0; + int nNonOrphans = 0; // find the currently selected readout @@ -583,5 +645,8 @@ for (int i = 0; i < sourcesSrc->n; i++) { pmSource *sourceSrc = sourcesSrc->data[i]; - if (sourceSrc->parent) continue; // I am not an orphan + if (sourceSrc->parent) { + nNonOrphans ++; + continue; // Not an orphan + } pmSource *sourceOut = pmSourceCopy(sourceSrc); @@ -609,6 +674,4 @@ pmSourceFreePixels (sourceOut); - // XXX do we need to skip the Chisq image sources? - // allocate image, weight, mask for the new image for each peak if (sourceOut->modelPSF) { @@ -620,9 +683,9 @@ sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + nParents ++; psArrayAdd (detectionsOut->allSources, 100, sourceOut); psFree (sourceOut); } - psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n); - + psLogMsg ("psphot", 3, "%d parents created, %d unorphaned children, %ld input vs %ld output", nParents, nNonOrphans, sourcesSrc->n, detectionsOut->allSources->n); return true; @@ -651,11 +714,7 @@ } -// create source children from ruleSrc for ruleOut for this entry. XXX currently, this is only +// Create source children from ruleSrc for ruleOut for this entry. Currently, this is only // used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called -// repeatedly) - -// XXX this function will do bad things if called repeatedly. It needs to either (1) delete -// the existing DETECTION container on the output images, or (2) intelligently supplement the -// existing DETECTION container with only the new detections. +// repeatedly). bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { @@ -682,15 +741,10 @@ psAssert (readoutOut, "missing readout?"); - // generate a new detection structure for the output filerule - pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS"); - if (!detectionsOut) { - detectionsOut = pmDetectionsAlloc(); - detectionsOut->allSources = psArrayAllocEmpty (100); - // save detections on the readout->analysis - if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) { - psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); - return false; - } - psFree(detectionsOut); // a copy remains on the analysis metadata + // replace any existing DETECTION container on readoutOut->analysis with the new one + pmDetections *detectionsOut = pmDetectionsAlloc(); + detectionsOut->allSources = psArrayAllocEmpty (100); + if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; } @@ -729,7 +783,5 @@ // allocate image, weight, mask for the new image for each peak - if (sourceOut->modelPSF) { - pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); - } + pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->windowRadius); // child sources have not been subtracted in this image, but this flag may be raised if @@ -740,6 +792,7 @@ psFree (sourceOut); } - psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n); - + psLogMsg ("psphot", 3, "created %ld children", detectionsOut->allSources->n); + + psFree(detectionsOut); // a copy remains on the analysis metadata return true; Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotOutput.c (revision 32689) @@ -414,5 +414,6 @@ char filename[64]; - // XXX uncomment to disreturn true; + // XXX dump tests are disabled unless this is commented out: + return true; bool status = true; @@ -453,14 +454,15 @@ // 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; + // 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); + // float Npix = source->moments ? source->moments->nPixels : NAN; + float Io = source->modelPSF ? source->modelPSF->params->data.F32[PM_PAR_I0] : NAN; + fprintf (f, "%d %f %f : %f %f : %f %f\n", source->imageID, peak->xf, peak->yf, Mx, My, Msum, Io); } } Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotReadout.c (revision 32689) @@ -1,4 +1,3 @@ # include "psphotInternal.h" -bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule); // this should be called by every program that links against libpsphot @@ -83,5 +82,4 @@ return psphotReadoutCleanup (config, view, filerule); } - psphotDumpTest (config, view, filerule); // find blended neighbors of very saturated stars (detections->newSources) @@ -139,5 +137,4 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) - psphotDumpTest (config, view, filerule); // measure the radial profiles to the sky @@ -188,5 +185,4 @@ // NOTE: new sources are saved on detections->newSources psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) - psphotDumpTest (config, view, filerule); // set source type @@ -209,5 +205,4 @@ // 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 Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c (revision 32689) @@ -1,6 +1,3 @@ # include "psphotInternal.h" - -static int replace_pass = 0; -static int remove_pass = 0; // replace the flux for sources which failed @@ -43,5 +40,4 @@ } } - replace_pass ++; return true; } @@ -60,7 +56,82 @@ 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?"); + + 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]; + + 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, do not replace unsubtracted sources + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + } + + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + psphotVisualShowImage(readout); + psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); + return true; +} + +// 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; + } + } + 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?"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -98,95 +169,9 @@ } 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); - psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); - return true; -} - -// 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; - + if ((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); Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotSourceFits.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotSourceFits.c (revision 32689) @@ -257,4 +257,5 @@ okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); + okDBL = false; // XXX this is failing badly... // XXX should I keep / save the flags set in the eval functions? Index: /branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c =================================================================== --- /branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c (revision 32688) +++ /branches/eam_branches/ipp-20111110/psphot/src/psphotStackReadout.c (revision 32689) @@ -1,4 +1,3 @@ # include "psphotInternal.h" -bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule); // we have 3 possible real filesets: @@ -83,8 +82,4 @@ } - // 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) if (!psphotModelBackground (config, view, STACK_DET)) { @@ -92,4 +87,7 @@ } if (!psphotSubtractBackground (config, view, STACK_DET)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } + if (!psphotSubtractBackground (config, view, STACK_SRC)) { return psphotReadoutCleanup (config, view, STACK_SRC); } @@ -133,5 +131,5 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } - psphotDumpTest (config, view, STACK_SRC); + // psphotDumpTest (config, view, STACK_SRC); psMemDump("sourcestats"); @@ -170,9 +168,6 @@ // 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); - psphotDumpTest (config, view, STACK_SRC); // measure the radial profiles to the sky @@ -206,14 +201,15 @@ // NOTE: this block performs the 2nd pass low-significance PSF detection stage { - // XXX the steps below to subtract the sources from DET should be skipped if SRC == DET! - - // generate children sources for all sources in the SRC image - // XXX (DROP existing detections list) - psphotSourceChildren (config, view, STACK_DET, STACK_SRC); - - // subtract all sources from DET (this will subtract using the psf model for SRC, which - // will somewhat oversubtract the sources -- this is OK - // *** this fails because the source->pixels are still pointing at SRC, not DET - psphotRemoveAllSources (config, view, STACK_DET, true); // ignore subtraction state for sources + // if DET and SRC are different images, generate children sources for all sources in + // the SRC image. This operation replaces the existing DETECTION container on DET + // which is currently a view to the one on SRC). children sources go to + // det->allSources + if (strcmp(STACK_SRC, STACK_DET)) { + psphotSourceChildren (config, view, STACK_DET, STACK_SRC); + + // 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, false); // ignore subtraction state for sources + } // add noise for subtracted objects @@ -231,12 +227,9 @@ // (this operation just ensures the metadata container has a view on SRC as well if (strcmp(STACK_SRC, STACK_DET)) { - // XXX this operation now needs to create source parents for the new detections - // XXX if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) { - // XXX psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - // XXX return psphotReadoutCleanup (config, view, STACK_SRC); - // XXX } - - // the old detection here are children of SRC; generate parents for the new detections - if (!psphotSourceParents (config, view, STACK_SRC, STACK_DET)) { + // replace all sources in DET + psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources + + // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them + if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) { psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); return psphotReadoutCleanup (config, view, STACK_SRC); @@ -244,11 +237,7 @@ } - // 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 @@ -271,5 +260,4 @@ // XXX check on free of sources... psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources) - psphotDumpTest (config, view, STACK_SRC); } @@ -285,5 +273,4 @@ psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); - psphotDumpTest (config, view, STACK_SRC); // NOTE: apply to ALL sources @@ -324,6 +311,4 @@ int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); - - psphotDumpTest (config, view, STACK_SRC); for (int entry = 1; entry < nRadialEntries; entry++) {