Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 32666) +++ trunk/psphot/src/psphot.h (revision 32695) @@ -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); @@ -453,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: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 32666) +++ trunk/psphot/src/psphotBlendFit.c (revision 32695) @@ -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: trunk/psphot/src/psphotEfficiency.c =================================================================== --- trunk/psphot/src/psphotEfficiency.c (revision 32666) +++ trunk/psphot/src/psphotEfficiency.c (revision 32695) @@ -255,5 +255,5 @@ // remove all sources, adding noise for subtracted sources - psphotRemoveAllSources(realSources, recipe); + psphotRemoveAllSourcesByArray(realSources, recipe); #if TESTING Index: trunk/psphot/src/psphotFake.c =================================================================== --- trunk/psphot/src/psphotFake.c (revision 32666) +++ trunk/psphot/src/psphotFake.c (revision 32695) @@ -166,5 +166,5 @@ // remove all sources, adding noise for subtracted sources - psphotRemoveAllSources(realSources, recipe); + psphotRemoveAllSourcesByArray(realSources, recipe); psphotAddNoise(readout, realSources, recipe); Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 32666) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 32695) @@ -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: trunk/psphot/src/psphotMergeSources.c =================================================================== --- trunk/psphot/src/psphotMergeSources.c (revision 32666) +++ trunk/psphot/src/psphotMergeSources.c (revision 32695) @@ -520,4 +520,176 @@ } +// 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; + + 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 (!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; + } + } + return true; +} + +// 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 + 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 *detectionsSrc = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS"); + psAssert (detectionsSrc, "missing detections?"); + + psArray *sourcesSrc = detectionsSrc->allSources; + psAssert (sourcesSrc, "missing sources?"); + + // 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?"); + + // loop over the sources, redefine their pixels to point at the new filerule image, + // copy the source data, and add a reference back to the original source + + // copy the sources from sourceSrcs to the new detection structure + for (int i = 0; i < sourcesSrc->n; i++) { + pmSource *sourceSrc = sourcesSrc->data[i]; + if (sourceSrc->parent) { + nNonOrphans ++; + continue; // Not an orphan + } + + pmSource *sourceOut = pmSourceCopy(sourceSrc); + sourceOut->parent = sourceSrc; + + // keep the original source flags + sourceOut->seq = sourceSrc->seq; + sourceOut->type = sourceSrc->type; + sourceOut->mode = sourceSrc->mode; + sourceOut->mode2 = sourceSrc->mode2; + sourceOut->tmpFlags = sourceSrc->tmpFlags; + + // does this copy all model data? (NO) + sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF); + sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT); + + if (sourceSrc->modelFits) { + sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n); + for (int j = 0; j < sourceSrc->modelFits->n; j++) { + sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]); + } + } + + // drop the references to the original image pixels: + pmSourceFreePixels (sourceOut); + + // 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); + } + + // child sources have not been subtracted in this image, but this flag may be raised if + // they were subtracted in the parent's image + sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + + nParents ++; + psArrayAdd (detectionsOut->allSources, 100, sourceOut); + psFree (sourceOut); + } + psLogMsg ("psphot", 3, "%d parents created, %d unorphaned children, %ld input vs %ld output", nParents, nNonOrphans, sourcesSrc->n, detectionsOut->allSources->n); + + return true; +} + // create source children from ruleSrc for ruleOut bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) @@ -542,7 +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) +// repeatedly). bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { @@ -569,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; } @@ -593,4 +760,5 @@ // keep the original source flags + sourceOut->seq = sourceSrc->seq; sourceOut->type = sourceSrc->type; sourceOut->mode = sourceSrc->mode; @@ -612,8 +780,8 @@ 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) { - 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 @@ -624,13 +792,14 @@ psFree (sourceOut); } - psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n); - - - return true; -} - -// create source children associated with 'filerule' from the objectsSrc. XXX currently, this -// is only used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be -// called repeatedly) + psLogMsg ("psphot", 3, "created %ld children", detectionsOut->allSources->n); + + psFree(detectionsOut); // a copy remains on the analysis metadata + + return true; +} + +// create source children associated with 'filerule' from the objectsSrc. returns a new object +// array containing the child sources. XXX currently, this is only used by psphotStackReadout +// (sources go on allSources so that psphotChoosePSF can be called repeatedly) psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) { @@ -652,4 +821,5 @@ psAssert (readout, "missing readout?"); + // create DETECTIONS containers for each image, in case one lacks it pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); if (!detections) { @@ -665,4 +835,6 @@ psAssert (detections, "missing detections?"); } + + // we need to save the new sources on the detection arrays of the appropriate image detArrays->data[i] = psMemIncrRefCounter(detections); readouts->data[i] = psMemIncrRefCounter(readout); @@ -695,7 +867,8 @@ pmSource *sourceOut = pmSourceCopy(sourceSrc); + sourceOut->parent = sourceSrc; + + // save on the output object array at the same location objectOut->sources->data[i] = sourceOut; - - sourceOut->parent = sourceSrc; // keep the original source flags and sequence ID (if set) @@ -720,5 +893,5 @@ pmSourceFreePixels (sourceOut); - // set the output readotu + // set the output readout int index = sourceOut->imageID; if (index >= readouts->n) continue; // skip the sources generated by the chisq image Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 32666) +++ trunk/psphot/src/psphotOutput.c (revision 32695) @@ -407,4 +407,71 @@ 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 dump tests are disabled unless this is commented out: + return 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; + 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); + } + } + fclose (f); + npass ++; + + return true; +} + +# if (0) bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule, char *filename) { @@ -449,3 +516,3 @@ return true; } - +# endif Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 32666) +++ trunk/psphot/src/psphotReadout.c (revision 32695) @@ -9,55 +9,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) { @@ -188,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) - psphotDumpChisqs (config, view, filerule); // measure the radial profiles to the sky @@ -208,13 +156,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 @@ -253,5 +199,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 @@ -262,5 +208,4 @@ // NOTE: apply to ALL sources psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) - psphotDumpChisqs (config, view, filerule); } @@ -295,5 +240,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: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 32666) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 32695) @@ -9,9 +9,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 +23,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,5 +35,5 @@ // 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; @@ -43,8 +43,7 @@ } -bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { - - bool status; - pmSource *source; +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) { + + bool status; psTimerStart ("psphot.replace"); @@ -75,10 +74,24 @@ 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, do not replace unsubtracted sources + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + } + + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } @@ -88,5 +101,84 @@ } -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; + } + } + 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"); + 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, only remove unsubtracted sources + if ((source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; + } + + pmSourceSub (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; +} + +bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) { bool status; @@ -100,10 +192,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 +250,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 +321,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: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 32666) +++ trunk/psphot/src/psphotSourceFits.c (revision 32695) @@ -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: trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- trunk/psphot/src/psphotStackChisqImage.c (revision 32666) +++ trunk/psphot/src/psphotStackChisqImage.c (revision 32695) @@ -13,4 +13,6 @@ psAssert (chisqFile, "missing chisq image FPA?"); + // the readout containing the chisq image is generated in the first pass of this loop and + // used by the successive passes pmReadout *chiReadout = NULL; Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 32666) +++ trunk/psphot/src/psphotStackReadout.c (revision 32695) @@ -82,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)) { @@ -91,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); } @@ -114,5 +113,7 @@ } - // if DET and SRC are different images, copy the detections from DET to SRC + // If DET and SRC are different images, copy the detections from DET to SRC. This 'copy' + // is just a copy of the container pointer; the sources on both DET and SRC are the same + // memory objects if (strcmp(STACK_SRC, STACK_DET)) { if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) { @@ -130,4 +131,5 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } + // psphotDumpTest (config, view, STACK_SRC); psMemDump("sourcestats"); @@ -166,6 +168,4 @@ // 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); @@ -186,6 +186,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 +194,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), not DET psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources) @@ -200,4 +201,16 @@ // NOTE: this block performs the 2nd pass low-significance PSF detection stage { + // 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 psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources) @@ -212,7 +225,11 @@ // 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)) { + // 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); @@ -237,5 +254,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 @@ -312,5 +329,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