Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 16611) +++ trunk/psphot/src/Makefile.am (revision 16820) @@ -27,4 +27,5 @@ psphotModelBackground.c \ psphotSubtractBackground.c \ + psphotFindDetections.c \ psphotFindPeaks.c \ psphotSourceStats.c \ Index: trunk/psphot/src/pmFootprint.c =================================================================== --- trunk/psphot/src/pmFootprint.c (revision 16611) +++ trunk/psphot/src/pmFootprint.c (revision 16820) @@ -1241,2 +1241,51 @@ return peaks; } + +void pmDetectionsFree (pmDetections *detections) { + + if (!detections) return; + + psFree (detections->footprints); + psFree (detections->peaks); + psFree (detections->oldPeaks); + return; +} + +// generate a pmDetections container with empty (allocated) footprints and peaks containers +pmDetections *pmDetectionsAlloc() { + + pmDetections *detections = (pmDetections *)psAlloc(sizeof(pmDetections)); + psMemSetDeallocator(detections, (psFreeFunc) pmDetectionsFree); + + detections->footprints = NULL; + detections->peaks = NULL; + detections->oldPeaks = NULL; + detections->last = 0; + + return (detections); +} + +/************************************************************************************************************/ +/* + * Cull a set of peaks contained in a psArray of pmFootprints + */ +psErrorCode +psphotCullPeaks(const psImage *image, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + const psMetadata *recipe, + psArray *footprints) { // array of pmFootprints + bool status = false; + float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA"); + if (!status) { + nsigma_delta = 0; // min. + } + float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); + if (!status) { + nsigma_min = 0; + } + const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); + + return pmFootprintArrayCullPeaks(image, weight, footprints, + nsigma_delta, nsigma_min*skyStdev); +} + Index: trunk/psphot/src/pmFootprint.h =================================================================== --- trunk/psphot/src/pmFootprint.h (revision 16611) +++ trunk/psphot/src/pmFootprint.h (revision 16820) @@ -52,3 +52,13 @@ psArray *pmFootprintArrayToPeaks(const psArray *footprints); +typedef struct { + psArray *footprints; // collection of footprints in the image + psArray *peaks; // collection of all peaks contained by the footprints + psArray *oldPeaks; // collection of all peaks previously found + int last; +} pmDetections; + +pmDetections *pmDetectionsAlloc (); + + #endif Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 16611) +++ trunk/psphot/src/psphot.h (revision 16820) @@ -9,4 +9,6 @@ #include "psphotErrorCodes.h" +#include "pmFootprint.h" + #define PSPHOT_RECIPE "PSPHOT" // Name of the recipe to use @@ -16,9 +18,8 @@ psString psphotVersionLong(void); -bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark); +bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe); bool psphotReadout (pmConfig *config, const pmFPAview *view); -bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources); +bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources); bool psphotDefineFiles (pmConfig *config, pmFPAfile *input); - // XXX test functions @@ -26,35 +27,36 @@ // psphotReadout functions -bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename, psMaskType maskVal); -bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename, psMaskType maskVal) ; -psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, - const bool returnFootprints, const int pass, psMaskType maskVal); -#include "pmFootprint.h" -psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, - const psMetadata *recipe, psArray *footprints); -psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *allpeaks, psMaskType maskVal, psMaskType mark); -bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump, psMaskType maskSat); +bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename); +bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename) ; +pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe); + +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections); +bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump); bool psphotBasicDeblend (psArray *sources, psMetadata *recipe); -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark); +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); -#if NOT_IN_LIBPSPHOT -bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); -#endif -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal); -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); -bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); bool psphotReplaceUnfit (psArray *sources, psMaskType maskVal); -bool psphotReplaceAll (psArray *sources, psMaskType maskVal); -bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark); +bool psphotReplaceAll (psArray *sources, psMetadata *recipe); +bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background); +bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); +bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe); + +// used by psphotFindDetections +psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, const bool returnFootprints, const int pass, psMaskType maskVal); +psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints); + +// used by ApResid bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup); bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag); -bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark); -bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); // basic support functions void psphotModelClassInit (void); bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psMaskType maskVal); -bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal); +bool psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe); void psphotSourceFreePixels (psArray *sources); @@ -100,7 +102,9 @@ pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe); bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources); -bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal); +bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add); bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source); -bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal); +bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY); @@ -109,9 +113,8 @@ bool psphotSetState (pmSource *source, bool curState, psMaskType maskVal); bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe); -bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, long first); bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); -bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal); bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal); psKernel *psphotKernelFromPSF (pmSource *source, int nPix); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 16611) +++ trunk/psphot/src/psphotAddNoise.c (revision 16820) @@ -1,5 +1,13 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal) { +bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) { + return psphotAddOrSubNoise (readout, sources, recipe, true); +} + +bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) { + return psphotAddOrSubNoise (readout, sources, recipe, false); +} + +bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) { bool status = false; @@ -13,4 +21,8 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // increase variance by factor*(object noise): Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 16611) +++ trunk/psphot/src/psphotApResid.c (revision 16820) @@ -3,5 +3,6 @@ # define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; } // measure the aperture residual statistics and 2D variations -bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark) + +bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { int Nfail = 0; @@ -19,4 +20,12 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType mark = psMetadataLookupU8(&status, recipe, "MASK.MARK"); // Mask value for bad pixels + assert (mark); // S/N limit to perform full non-linear fits Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 16611) +++ trunk/psphot/src/psphotBlendFit.c (revision 16820) @@ -2,5 +2,5 @@ // XXX I don't like this name -bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) { +bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { int Nfit = 0; @@ -11,4 +11,8 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // source analysis is done in S/N order (brightest first) @@ -73,5 +77,6 @@ // try fitting PSFs, then try extended sources // these functions subtract the resulting fitted source (XXX and update the modelFlux?) - // XXX consider conditions under which the source has EXT fit + // XXX re-consider conditions under which the source has EXT fit: + // I should try EXT if the source size measurement says it is large if (psphotFitBlend (readout, source, psf, maskVal)) { psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); @@ -92,5 +97,8 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode |= PM_SOURCE_MODE_TEMPSUB; + } + + if (psTraceGetLevel("psphot") >= 6) { + psphotSaveImage (NULL, readout->image, "image.v2.fits"); } Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 16611) +++ trunk/psphot/src/psphotChoosePSF.c (revision 16820) @@ -2,9 +2,17 @@ // try PSF models and select best option -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { bool status; psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType mark = psMetadataLookupU8(&status, recipe, "MASK.MARK"); // Mask value for bad pixels + assert (mark); // examine PSF sources in S/N order (brightest first) Index: trunk/psphot/src/psphotEvalPSF.c =================================================================== --- trunk/psphot/src/psphotEvalPSF.c (revision 16611) +++ trunk/psphot/src/psphotEvalPSF.c (revision 16820) @@ -134,6 +134,6 @@ // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[] keep = TRUE; - keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); - keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); + // keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); + // keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); keep &= (SN > PSF_MIN_SN); keep &= (Chi < PSF_MAX_CHI); Index: trunk/psphot/src/psphotExtendedSources.c =================================================================== --- trunk/psphot/src/psphotExtendedSources.c (revision 16611) +++ trunk/psphot/src/psphotExtendedSources.c (revision 16820) @@ -1,9 +1,13 @@ # include "psphot.h" -bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal) { +bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe) { bool status; int Next = 0; int Npsf = 0; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // S/N limit to perform full non-linear fits @@ -117,5 +121,4 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode |= PM_SOURCE_MODE_TEMPSUB; } Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 16820) +++ trunk/psphot/src/psphotFindDetections.c (revision 16820) @@ -0,0 +1,78 @@ +# include "psphotInternal.h" + +// XXX let's make a better distinction between 'pass' (used for counting intput / output +// entries) and needing/having a PSF. In this function, we smooth the image, then search for +// the peaks +pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe) { + + bool status; + int pass; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // Use the new pmFootprints approach? + const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS"); + + // on first pass, detections have not yet been allocated + if (!detections) { + detections = pmDetectionsAlloc(); + pass = 1; + } else { + pass = 2; + } + + // move the old peaks array (if it exists) to oldPeaks + assert (detections->oldPeaks == NULL); + detections->oldPeaks = detections->peaks; + detections->peaks = NULL; + + // if footprints are not requested, find the peaks and return them + if (!useFootprints) { + detections->peaks = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal); + return detections; + } + + // psphotFindPeaks returns a different set of objects if useFootprints is true! + // XXX fix this by splitting psphotFindPeaks into smooth / peaks / footprint stages? + psArray *footprints = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal); + + int growRadius = 0; + if (pass == 1) { + growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS"); + } else { + growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2"); + } + + if (growRadius > 0) { + psArray *tmp = pmGrowFootprintArray(footprints, growRadius); + psFree(footprints); + footprints = tmp; + } + + if (pass == 2) { + // merge in old peaks; + const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints + + psLogMsg ("psphot", PS_LOG_INFO, "merging %ld new footprints into %ld existing ones\n", footprints->n, detections->footprints->n); + psArray *mergedFootprints = pmMergeFootprintArrays(detections->footprints, footprints, includePeaks); + psFree(footprints); + psFree(detections->footprints); + detections->footprints = mergedFootprints; + } else { + detections->footprints = footprints; + } + + psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); + detections->peaks = pmFootprintArrayToPeaks(detections->footprints); + psLogMsg ("psphot", PS_LOG_INFO, "peaks: %ld, total footprints: %ld\n", detections->peaks->n, detections->footprints->n); + + return detections; +} + +// if we use the footprints, the output peaks list contains both old and new peaks, +// otherwise it only contains the new peaks. + +// if psf is defined, we should treat this as pass "2", not "1". re-define this boolean to +// be "have PSF" Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 16611) +++ trunk/psphot/src/psphotFindPeaks.c (revision 16820) @@ -166,4 +166,6 @@ psFree(peaks); peaks = footprints; // well, you know what I mean + + psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("psphot")); } @@ -174,26 +176,2 @@ } - -/************************************************************************************************************/ -/* - * Cull a set of peaks contained in a psArray of pmFootprints - */ -psErrorCode -psphotCullPeaks(const psImage *image, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - const psMetadata *recipe, - psArray *footprints) { // array of pmFootprints - bool status = false; - float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA"); - if (!status) { - nsigma_delta = 0; // min. - } - float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); - if (!status) { - nsigma_min = 0; - } - const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); - - return pmFootprintArrayCullPeaks(image, weight, footprints, - nsigma_delta, nsigma_min*skyStdev); -} Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 16611) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 16820) @@ -12,5 +12,5 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal) { +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { bool status; @@ -21,4 +21,8 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // source analysis is done in spatial order @@ -187,6 +191,4 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; - if (!final) source->mode |= PM_SOURCE_MODE_TEMPSUB; - // XXX not sure about the use of TEMPSUB } Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 16611) +++ trunk/psphot/src/psphotGuessModels.c (revision 16820) @@ -17,76 +17,82 @@ // construct an initial PSF model for each object -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) { +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { - psTimerStart ("psphot"); + bool status; - // setup the PSF fit radius details - psphotInitRadiusPSF (recipe, psf->type); + psTimerStart ("psphot"); - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); - // skip non-astronomical objects (very likely defects) - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + // setup the PSF fit radius details + psphotInitRadiusPSF (recipe, psf->type); - // XXX if a source is faint, it will not have moments measured. - // it must be modelled as a PSF. In this case, we need to use - // the peak centroid to get the coordinates and get the peak flux - // from the image? - pmModel *modelEXT; - if (!source->moments) { - modelEXT = wildGuess(source, psf); - } else { - // use the source moments, etc to guess basic model parameters - modelEXT = pmSourceModelGuess (source, psf->type); - if (!modelEXT) { - modelEXT = wildGuess(source, psf); - } - // these valuse are set in pmSourceModelGuess, should this rule be in there as well? - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->x; - modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y; - } else { - modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf; - modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf; - } + 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; + + // XXX if a source is faint, it will not have moments measured. + // it must be modelled as a PSF. In this case, we need to use + // the peak centroid to get the coordinates and get the peak flux + // from the image? + pmModel *modelEXT; + if (!source->moments) { + modelEXT = wildGuess(source, psf); + } else { + // use the source moments, etc to guess basic model parameters + modelEXT = pmSourceModelGuess (source, psf->type); + if (!modelEXT) { + modelEXT = wildGuess(source, psf); + } + // these valuse are set in pmSourceModelGuess, should this rule be in there as well? + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->x; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y; + } else { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf; + } + } + + // set PSF parameters for this model (apply 2D shape model) + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + if (modelPSF == NULL) { + psError(PSPHOT_ERR_PSF, false, + "Failed to determine PSF model at r,c = (%d,%d); trying centre of image", + source->peak->y, source->peak->x); + // + // Try the centre of the image + // + modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols; + modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows; + modelPSF = pmModelFromPSF (modelEXT, psf); + if (modelPSF == NULL) { + psError(PSPHOT_ERR_PSF, false, + "Failed to determine PSF model at centre of image"); + psFree(modelEXT); + return false; + } + + source->mode |= PM_SOURCE_MODE_BADPSF; + } + psFree (modelEXT); + + // XXX need to define the guess flux? + // set the fit radius based on the object flux limit and the model + psphotCheckRadiusPSF (readout, source, modelPSF); + + // set the source PSF model + source->modelPSF = modelPSF; + source->modelPSF->residuals = psf->residuals; + + pmSourceCacheModel (source, maskVal); } - - // set PSF parameters for this model (apply 2D shape model) - pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); - if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at r,c = (%d,%d); trying centre of image", - source->peak->y, source->peak->x); - // - // Try the centre of the image - // - modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols; - modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows; - modelPSF = pmModelFromPSF (modelEXT, psf); - if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at centre of image"); - psFree(modelEXT); - return false; - } - - source->mode |= PM_SOURCE_MODE_BADPSF; - } - psFree (modelEXT); - - // XXX need to define the guess flux? - // set the fit radius based on the object flux limit and the model - psphotCheckRadiusPSF (readout, source, modelPSF); - - // set the source PSF model - source->modelPSF = modelPSF; - source->modelPSF->residuals = psf->residuals; - - pmSourceCacheModel (source, maskVal); - } - psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); - return true; + psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + return true; } Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 16611) +++ trunk/psphot/src/psphotMagnitudes.c (revision 16820) @@ -1,15 +1,18 @@ # include "psphotInternal.h" -bool psphotMagnitudes(psArray *sources, - psMetadata *recipe, - pmPSF *psf, - pmReadout *background, - psMaskType maskVal, - psMaskType mark) -{ +bool psphotMagnitudes(psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background) { + bool status = false; int Nap = 0; psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType mark = psMetadataLookupU8(&status, recipe, "MASK.MARK"); // Mask value for bad pixels + assert (mark); pmSourceMagnitudesInit (recipe); Index: trunk/psphot/src/psphotMaskReadout.c =================================================================== --- trunk/psphot/src/psphotMaskReadout.c (revision 16611) +++ trunk/psphot/src/psphotMaskReadout.c (revision 16820) @@ -1,7 +1,47 @@ # include "psphotInternal.h" -bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal) { +// generate mask and weight if not defined, additional mask for restricted subregion +bool psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe) { bool status; + + // ** Interpret the mask values: + + // single-bit named masks + psMaskType maskMark = pmConfigMask("MARK", config); // Mask value for marking + psMetadataAddU8 (recipe, PS_LIST_TAIL, "MASK.MARK", PS_META_REPLACE, "user-defined mask", maskMark); + + psMaskType maskSat = pmConfigMask("SAT", config); // Mask value for saturated pixels + psMetadataAddU8 (recipe, PS_LIST_TAIL, "MASK.SAT", PS_META_REPLACE, "user-defined mask", maskSat); + + psMaskType maskBad = pmConfigMask("BAD", config); // Mask value for bad pixels + psMetadataAddU8 (recipe, PS_LIST_TAIL, "MASK.BAD", PS_META_REPLACE, "user-defined mask", maskBad); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "PSPHOT"); // Mask value for bad pixels + if (!maskVal) { + const char *maskValStr = psMetadataLookupStr(NULL, recipe, "MASKVAL"); // String with mask names + if (!maskValStr) { + psError(PSPHOT_ERR_CONFIG, false, "Missing recipe item: MASKVAL(STR)"); + return false; + } + maskVal = pmConfigMask(maskValStr, config); // Mask values to mask against + psMetadataAddU8 (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal); + assert (maskVal); + } + + // generate mask & weight images if they don't already exit + if (!readout->mask) { + if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask"); + return false; + } + } + if (!readout->weight) { + if (!pmReadoutGenerateWeight(readout, true)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight"); + return false; + } + } // mask the excluded outer pixels @@ -19,5 +59,12 @@ // psImageKeepRegion assumes the region refers to the parent coordinates - psImageKeepRegion (readout->mask, keep, "OR", maskVal); + psImageKeepRegion (readout->mask, keep, "OR", maskBad); + + // test output of files at this stage + if (psTraceGetLevel("psphot") >= 5) { + psphotSaveImage (NULL, readout->image, "image.fits"); + psphotSaveImage (NULL, readout->mask, "mask.fits"); + psphotSaveImage (NULL, readout->weight, "weight.fits"); + } return true; Index: trunk/psphot/src/psphotModelBackground.c =================================================================== --- trunk/psphot/src/psphotModelBackground.c (revision 16611) +++ trunk/psphot/src/psphotModelBackground.c (revision 16820) @@ -4,5 +4,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename, psMaskType maskVal) +bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename) { bool status = true; @@ -20,4 +20,9 @@ // select the appropriate recipe information psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // user supplied seed, if available Index: trunk/psphot/src/psphotModelTest.c =================================================================== --- trunk/psphot/src/psphotModelTest.c (revision 16611) +++ trunk/psphot/src/psphotModelTest.c (revision 16820) @@ -3,5 +3,5 @@ // XXX add more test information? -bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { +bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { bool status; @@ -11,4 +11,10 @@ pmPSF *psf = NULL; pmSourceFitMode fitMode; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + psMaskType mark = psMetadataLookupU8(&status, recipe, "MASK.MARK"); // Mask value for bad pixels + assert (maskVal); + assert (mark); // run model fitting tests on a single source? Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 16611) +++ trunk/psphot/src/psphotReadout.c (revision 16820) @@ -9,6 +9,4 @@ psTimerStart ("psphotReadout"); - bool dump = (psTraceGetLevel("psphot") >= 6); - // select the current recipe psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); @@ -18,40 +16,4 @@ } - // Interpret the mask values - const char *maskValStr = psMetadataLookupStr(NULL, recipe, "MASKVAL"); // String with mask names - if (!maskValStr) { - psError(PSPHOT_ERR_CONFIG, false, "Missing recipe item: MASKVAL(STR)"); - return false; - } - psMaskType maskVal = pmConfigMask(maskValStr, config); // Mask values to mask against - psMaskType maskMark = pmConfigMask("MARK", config); // Mask value for marking - psMaskType maskSat = pmConfigMask("SAT", config); // Mask value for saturated pixels - psMaskType maskBad = pmConfigMask("BAD", config); // Mask value for bad pixels - - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); - - // optional break-point for processing - char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); - PS_ASSERT_PTR_NON_NULL (breakPt, false); - - // Use the new pmFootprints approach? - const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS"); - - // generate mask & weight images if they don't already exit - if (!readout->mask) { - if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) { - psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask"); - return false; - } - } - if (!readout->weight) { - if (!pmReadoutGenerateWeight(readout, true)) { - psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight"); - return false; - } - } - // set the photcode for this image if (!psphotAddPhotcode (recipe, config, view)) { @@ -60,31 +22,31 @@ } - // I have a valid mask, now mask in the analysis region of interest - psphotMaskReadout (readout, recipe, maskBad); - - if (psTraceGetLevel("psphot") >= 5) { - psphotSaveImage (NULL, readout->image, "image.fits"); - psphotSaveImage (NULL, readout->mask, "mask.fits"); - psphotSaveImage (NULL, readout->weight, "weight.fits"); - } - + // find the currently selected readout + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + PS_ASSERT_PTR_NON_NULL (breakPt, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL); + return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); } // generate a background model (median, smoothed image) - if (!psphotModelBackground (config, view, "PSPHOT.INPUT", maskVal)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); - } - if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT", maskVal)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); - } - + if (!psphotModelBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); - } - - // run a single-model test if desired - psphotModelTest (config, view, recipe, maskVal, maskMark); + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // run a single-model test if desired (exits from here if test is run) + psphotModelTest (config, view, recipe); // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe @@ -94,39 +56,20 @@ bool havePSF = (psf != NULL); - // find the peaks in the image. - - // XXX clean this up into a single function. if psf is defined, we should treat this as - // pass "2", not "1". re-define this boolean to be "have PSF" - psArray *peaks; - psArray *footprints = NULL; - if (useFootprints) { - footprints = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal); - int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS"); - if (growRadius > 0) { - psArray *tmp = pmGrowFootprintArray(footprints, growRadius); - psFree(footprints); - footprints = tmp; - } - psphotCullPeaks(readout->image, readout->weight, recipe, footprints); - - peaks = pmFootprintArrayToPeaks(footprints); - } else { - peaks = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal); - } - - if (!peaks) { - psLogMsg ("psphot", 3, "unable to find peaks in this image"); - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); + // find the detections (by peak and/or footprint) in the image. + pmDetections *detections = psphotFindDetections (NULL, readout, recipe); + if (!detections) { + psLogMsg ("psphot", 3, "unable to find detections in this image"); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); } // construct sources and measure basic stats - psArray *sources = psphotSourceStats (readout, recipe, peaks, maskVal, maskMark); + psArray *sources = psphotSourceStats (readout, recipe, detections); if (!sources) return false; - psFree (peaks); - if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, sources); - } - + return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + } + + // find blended neighbors of very saturated stars + // XXX merge this with Basic Deblend? psphotDeblendSatstars (sources, recipe); @@ -134,14 +77,14 @@ if (!psphotBasicDeblend (sources, recipe)) { psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, NULL, sources); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } // classify sources based on moments, brightness - if (!psphotRoughClass (sources, recipe, havePSF, maskSat)) { + if (!psphotRoughClass (sources, recipe, havePSF)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, NULL, sources); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, sources); + return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); } @@ -151,16 +94,13 @@ // XXX if we do not have enough stars to generate the PSF, build one // from the SEEING guess and model class - psf = psphotChoosePSF (readout, sources, recipe, maskVal, maskMark); + psf = psphotChoosePSF (readout, sources, recipe); if (psf == NULL) { psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, psf, sources); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } havePSF = true; } - - // XXX Test dump of PSF parameters across image field. - if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, readout, recipe, psf, sources); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } @@ -169,19 +109,15 @@ // construct an initial model for each object - psphotGuessModels (readout, sources, recipe, psf, maskVal); + psphotGuessModels (readout, sources, recipe, psf); // XXX test output of models // psphotTestSourceOutput (readout, sources, recipe, psf); - if (dump) psphotSaveImage (NULL, readout->image, "image.v0.fits"); - - // linear PSF fit to peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v1.fits"); + // linear PSF fit to source peaks + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX test the CR/EXT measurement // XXX we need to call this here and after the second pass: option for starting number? - psphotSourceSize (readout, sources, recipe); - + psphotSourceSize (readout, sources, recipe, 0); if (!strcasecmp (breakPt, "ENSEMBLE")) { goto finish; @@ -189,14 +125,11 @@ // non-linear PSF and EXT fit to brighter sources - psphotBlendFit (readout, sources, recipe, psf, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v2.fits"); - - // replace all sources - psphotReplaceAll (sources, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v3.fits"); - - // linear PSF fit to remaining peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v4.fits"); + psphotBlendFit (readout, sources, recipe, psf); + + // replace all sources (make this part of psphotFitSourcesLinear?) + psphotReplaceAll (sources, recipe); + + // linear fit to include all sources + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); if (!strcasecmp (breakPt, "PASS1")) { goto finish; @@ -210,62 +143,37 @@ // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe, true, maskVal); - - // find the peaks in the image - psArray *newPeaks; - if (useFootprints) { - psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2, maskVal); - - int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2"); - if (growRadius > 0) { - psArray *tmp = pmGrowFootprintArray(newFootprints, growRadius); - psFree(newFootprints); - newFootprints = tmp; - } - - // merge in old peaks - const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints - psArray *mergedFootprints = pmMergeFootprintArrays(footprints, newFootprints, includePeaks); - psFree(footprints); - psFree(newFootprints); - - psphotCullPeaks(readout->image, readout->weight, recipe, mergedFootprints); - - newPeaks = pmFootprintArrayToPeaks(mergedFootprints); - - psFree(mergedFootprints); - } else { - newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2, maskVal); - } + psphotAddNoise (readout, sources, recipe); + + detections = psphotFindDetections (detections, readout, recipe); // remove noise for subtracted objects - psphotAddNoise (readout, sources, recipe, false, maskVal); - - // define new sources based on the new peaks - psArray *newSources = psphotSourceStats (readout, recipe, newPeaks, maskVal, maskMark); - psFree (newPeaks); + psphotSubNoise (readout, sources, recipe); + + // define new sources based on only the new peaks + psArray *newSources = psphotSourceStats (readout, recipe, detections); // set source type - if (!psphotRoughClass (newSources, recipe, havePSF, maskSat)) { + if (!psphotRoughClass (newSources, recipe, havePSF)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, psf, sources); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } // create full input models - psphotGuessModels (readout, newSources, recipe, psf, maskVal); - - // replace all sources - psphotReplaceAll (sources, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v5.fits"); - - // merge the newly selected peaks into the existing list + psphotGuessModels (readout, newSources, recipe, psf); + + // replace all sources so fit below applies to all at once + psphotReplaceAll (sources, recipe); + + // merge the newly selected sources into the existing list psphotMergeSources (sources, newSources); psFree (newSources); - // linear PSF fit to remaining peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal); - if (dump) psphotSaveImage (NULL, readout->image, "image.v6.fits"); - - psphotExtendedSources (readout, sources, recipe, maskVal); + // linear fit to all sources + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + + // measure source size for the remaining sources + psphotSourceSize (readout, sources, recipe, 0); + + psphotExtendedSources (readout, sources, recipe); finish: @@ -273,17 +181,16 @@ // plot positive sources if (!havePSF) { - // psphotSourcePlots (readout, sources, recipe, maskVal); + // psphotSourcePlots (readout, sources, recipe); } // measure aperture photometry corrections - if (!havePSF && !psphotApResid (readout, sources, recipe, psf, maskVal, maskMark)) { - psTrace ("psphot", 4, "failure on psphotApResid"); - psError(PSPHOT_ERR_PHOTOM, false, "Measure aperture photometry corrections"); - return false; + if (!havePSF && !psphotApResid (readout, sources, recipe, psf)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } // calculate source magnitudes pmReadout *background = psphotSelectBackground (config, view, false); - psphotMagnitudes(sources, recipe, psf, background, maskVal, maskMark); + psphotMagnitudes(sources, recipe, psf, background); // replace failed sources? @@ -297,4 +204,4 @@ // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, psf, sources); + return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); } Index: trunk/psphot/src/psphotReadoutCleanup.c =================================================================== --- trunk/psphot/src/psphotReadoutCleanup.c (revision 16611) +++ trunk/psphot/src/psphotReadoutCleanup.c (revision 16820) @@ -1,8 +1,8 @@ # include "psphotInternal.h" -// psphotReadoutCleanup is called on exit from psphotReadout if the last raised -// error is not a DATA error, then there was a serious problem. only in this -// case, or if the fail on the stats measurement, do we return false -bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) { +// psphotReadoutCleanup is called on exit from psphotReadout. If the last raised error is +// not a DATA error, then there was a serious problem. Only in this case, or if the fail +// on the stats measurement, do we return false +bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources) { // remove internal pmFPAfiles, if created @@ -73,4 +73,5 @@ pmKapaClose (); + psFree (detections); psFree (psf); psFree (header); Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 16611) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 16820) @@ -18,5 +18,4 @@ pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; } psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot"), sources->n); @@ -24,9 +23,14 @@ } -bool psphotReplaceAll (psArray *sources, psMaskType maskVal) { +bool psphotReplaceAll (psArray *sources, psMetadata *recipe) { + bool status; pmSource *source; psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); for (int i = 0; i < sources->n; i++) { @@ -39,5 +43,5 @@ source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; } - psLogMsg ("psphot.replace", PS_LOG_INFO, "replace models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); return true; } Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 16611) +++ trunk/psphot/src/psphotRoughClass.c (revision 16820) @@ -1,3 +1,4 @@ # include "psphotInternal.h" + # define CHECK_STATUS(S,MSG) { \ if (!status) { \ @@ -7,5 +8,5 @@ // 2006.02.02 : no leaks -bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool havePSF, psMaskType maskSat) { +bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool havePSF) { bool status; @@ -13,4 +14,8 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskSat = psMetadataLookupU8(&status, recipe, "MASK.SAT"); // Mask value for bad pixels + assert (maskSat); if (!havePSF) { Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 16611) +++ trunk/psphot/src/psphotSourceFits.c (revision 16820) @@ -113,5 +113,4 @@ pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal); blend->mode |= PM_SOURCE_MODE_SUBTRACTED; - blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; } NfitBlend += modelSet->n; @@ -137,5 +136,4 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; return true; } @@ -176,6 +174,4 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; - return true; } @@ -275,8 +271,6 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); - - source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; return true; @@ -290,5 +284,4 @@ source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode |= PM_SOURCE_MODE_PAIR; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; // copy most data from the primary source (modelEXT, blends stay NULL) Index: trunk/psphot/src/psphotSourcePlots.c =================================================================== --- trunk/psphot/src/psphotSourcePlots.c (revision 16611) +++ trunk/psphot/src/psphotSourcePlots.c (revision 16820) @@ -1,8 +1,12 @@ # include "psphotInternal.h" -bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal) { +bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) { - // bool status = false; + bool status = false; psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); // the source images are written to an image 10x the size of a PSF object Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 16611) +++ trunk/psphot/src/psphotSourceSize.c (revision 16820) @@ -6,11 +6,25 @@ // compares the observed flux to the model. -bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, long first) { + + bool status; + + psTimerStart ("psphot"); + + float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT"); + assert (status); // loop over all source - for (int i = 0; i < sources->n; i++) { + for (int i = first; i < sources->n; i++) { pmSource *source = sources->data[i]; + // skip source if it was already measured + if (isfinite(source->crNsigma)) continue; + + source->crNsigma = -1.0; + source->extNsigma = -1.0; + // source must have been subtracted + source->crNsigma = -3.0; if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; @@ -23,4 +37,5 @@ // skip sources which are too close to a boundary + source->crNsigma = -4.0; if (xPeak < 1) continue; if (xPeak > source->pixels->numCols - 2) continue; @@ -29,4 +44,5 @@ // XXX for now, just skip any sources with masked pixels + source->crNsigma = -5.0; bool keep = true; for (int iy = -1; (iy <= +1) && keep; iy++) { @@ -97,5 +113,17 @@ source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0; // NOTE: abs needed to make the Nsigma value positive + + if (!isfinite(source->crNsigma)) { + fprintf (stderr, "."); + source->crNsigma = -6.0; + } + + if (source->crNsigma > CR_NSIGMA_LIMIT) { + source->mode |= PM_SOURCE_MODE_CRLIMIT; + } } + + psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot")); + return true; } Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 16611) +++ trunk/psphot/src/psphotSourceStats.c (revision 16820) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks, psMaskType maskVal, psMaskType mark) +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections) { bool status = false; @@ -8,4 +8,12 @@ psTimerStart ("psphot"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType mark = psMetadataLookupU8(&status, recipe, "MASK.MARK"); // Mask value for bad pixels + assert (mark); // determine properties (sky, moments) of initial sources @@ -20,4 +28,6 @@ char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); if (!status) return NULL; + + psArray *peaks = detections->peaks; sources = psArrayAllocEmpty (peaks->n); Index: trunk/psphot/src/psphotSubtractBackground.c =================================================================== --- trunk/psphot/src/psphotSubtractBackground.c (revision 16611) +++ trunk/psphot/src/psphotSubtractBackground.c (revision 16820) @@ -4,5 +4,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename, psMaskType maskVal) +bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename) { bool status = true; @@ -25,4 +25,9 @@ // select the appropriate recipe information psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");