Index: /branches/eam_branches/ipp-20101205/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphot.h (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphot.h (revision 30141) @@ -75,6 +75,6 @@ bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources); +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources); bool psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -404,6 +404,6 @@ bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); -bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int nextSize); -bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM); +bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize); +bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize); // psphotStackMatchPSFsUtils @@ -443,3 +443,10 @@ bool psphotCleanInputs (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + #endif + Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotChoosePSF.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotChoosePSF.c (revision 30141) @@ -2,5 +2,5 @@ // generate a PSF model for inputs without PSF models already loaded -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule) +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources) { bool status = true; @@ -19,5 +19,5 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) { + if (!psphotChoosePSFReadout (config, view, filerule, i, recipe, newSources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i); return false; @@ -28,5 +28,5 @@ // try PSF models and select best option -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources) { bool status; @@ -50,5 +50,5 @@ psAssert (detections, "missing detections?"); - psArray *sources = detections->newSources; + psArray *sources = newSources ? detections->newSources : detections->allSources; psAssert (sources, "missing sources?"); Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotCleanup.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotCleanup.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotCleanup.c (revision 30141) @@ -22,6 +22,6 @@ pmVisualCleanup (); psLibFinalize(); - // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); + // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); return; } Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotReadoutKnownSources.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotReadoutKnownSources.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotReadoutKnownSources.c (revision 30141) @@ -43,5 +43,5 @@ } - if (!psphotChoosePSF (config, view, filerule)) { + if (!psphotChoosePSF (config, view, filerule, true)) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); return psphotReadoutCleanup (config, view, filerule); Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotReplaceUnfit.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotReplaceUnfit.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotReplaceUnfit.c (revision 30141) @@ -75,4 +75,6 @@ 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; @@ -101,2 +103,151 @@ return true; } + +// modify the sources to point at the corresponding pixels for the given filerule +bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + 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 (!psphotRedefinePixelsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + pmSource *source; + + psTimerStart ("psphot.replace"); + + // 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?"); + + // XXX the sources have already been copied (merge into here?) + 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?"); + + 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; + + 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; + if (radius < 1) { + fprintf (stderr, "!"); + } + + pmSourceRedefinePixels (source, readout, Xo, Yo, radius, true); + } + return true; +} + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + 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 (!psphotResetModelsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + pmSource *source; + + psTimerStart ("psphot.replace"); + + // 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?"); + + // XXX the sources have already been copied (merge into here?) + 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?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + 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; + + // the guess central intensity comes from the peak: + float Io = source->peak->flux; + 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; + + // 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; + + pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) + } + + // psphotSaveImage(NULL, readoutIn->image, "image.in.fits"); + // psphotSaveImage(NULL, readoutOut->image, "image.out.sub.fits"); + // psphotVisualShowImage (readout); + + psLogMsg ("psphot.replace", PS_LOG_INFO, "subtracted models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); + return true; +} + Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotSourceMatch.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotSourceMatch.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotSourceMatch.c (revision 30141) @@ -243,4 +243,9 @@ peak->footprint = pmFootprintCopyData(footprint, readout->image); + // the peak does not claim ownership of the footprint (it does not free it). save a copy of this + // footprint on detections->footprints so we can free it later + psArrayAdd(detections->footprints, 100, peak->footprint); + psFree (peak->footprint); + // create a new source pmSource *source = pmSourceAlloc(); @@ -249,5 +254,5 @@ // add the peak - source->peak = psMemIncrRefCounter(peak); + source->peak = peak; // allocate space for moments @@ -261,5 +266,4 @@ psArrayAdd (detections->newSources, 100, source); psFree (source); - psFree (peak); } psFree (footprint); Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c (revision 30141) @@ -26,7 +26,4 @@ } - float targetFWHM = fwhmValues->data.F32[lastSize + 1]; - float currentFWHM = fwhmValues->data.F32[lastSize]; - int num = psphotFileruleCount(config, filerule); @@ -36,6 +33,6 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, currentFWHM, targetFWHM)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, targetFWHM); + if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, fwhmValues, lastSize)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, fwhmValues->data.F32[lastSize+1]); psImageConvolveSetThreads(oldThreads); return false; @@ -47,5 +44,5 @@ } -bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM) { +bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) { bool status = false; @@ -70,4 +67,7 @@ minGauss = 0.5; } + + float targetFWHM = fwhmValues->data.F32[lastSize + 1]; + float currentFWHM = fwhmValues->data.F32[lastSize]; if (targetFWHM <= currentFWHM) { @@ -131,7 +131,20 @@ // so we save this in a vector. if the vector is not yet defined, create it - psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); - psAssert(fwhmValues, "should already exist.."); - psVectorAppend(fwhmValues, targetFWHM); + psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + psAssert(fwhmValuesOut, "should already exist.."); + psVectorAppend(fwhmValuesOut, targetFWHM); + + // do not generate a PSF if we already were supplied one + pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + if (psfOld) { + // save PSF on readout->analysis + char psfEntry[64]; + snprintf (psfEntry, 64, "PSPHOT.PSF.V%d", lastSize); + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, psfEntry, PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psfOld)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF"); + } return true; Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c =================================================================== --- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c (revision 30140) +++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c (revision 30141) @@ -10,4 +10,24 @@ // SRC (source analysis image) : nominally CNV (optionally RAW) // OUT (psf-matched images) : always OUT + +bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) { + + int num = psphotFileruleCount(config, filerule); + + // 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?"); + + psphotVisualShowImage (readout); + psphotVisualShowResidualImage (readout); + } + return true; +} bool psphotStackReadout (pmConfig *config, const pmFPAview *view) { @@ -92,4 +112,8 @@ } + if (!strcasecmp (breakPt, "TEST1")) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } + // generate the objects (object unify the sources from the different images) // XXX this could just match the detections for the chisq image, and not bother measuring the @@ -97,6 +121,12 @@ 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_SRC, false)) { // pass 1 + psFree(objects); psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); return psphotReadoutCleanup (config, view, STACK_SRC); @@ -118,4 +148,5 @@ // 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); @@ -124,8 +155,10 @@ // 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); } @@ -133,9 +166,11 @@ // 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)) { // 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_SRC); } if (!strcasecmp (breakPt, "PSFMODEL")) { + psFree(objects); return psphotReadoutCleanup (config, view, STACK_SRC); } @@ -150,4 +185,5 @@ // 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 @@ -174,13 +210,24 @@ // copy the detections from SRC to OUT (for radial aperture photometry) if (!psphotCopySources (config, view, STACK_OUT, STACK_SRC)) { + psFree(objects); psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); return psphotReadoutCleanup (config, view, STACK_SRC); } - // XXX need to do something to reassign the source pixels here bool smoothAgain = true; for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) { + psphotRedefinePixels (config, view, STACK_OUT); + + psphotChoosePSF (config, view, STACK_OUT, false); + + psphotResetModels (config, view, STACK_OUT); + + psphotFitSourcesLinear (config, view, STACK_OUT, false); + // measure circular, radial apertures (objects sorted by S/N) psphotRadialAperturesByObject (config, objects, view, STACK_OUT, nMatchedPSF); + + psphotReplaceAllSources (config, view, STACK_OUT); + psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF); }