Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 13899) +++ /trunk/psphot/src/psphot.h (revision 13900) @@ -16,5 +16,5 @@ psString psphotVersionLong(void); -bool psphotModelTest (pmReadout *readout, psMetadata *recipe); +bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark); bool psphotReadout (pmConfig *config, const pmFPAview *view); bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources); @@ -26,14 +26,14 @@ // psphotReadout functions -bool psphotImageMedian (pmConfig *config, const pmFPAview *view); +bool psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal); psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, - const bool returnFootprints, const int pass); + 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); -bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump); +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *allpeaks, psMaskType maskVal, psMaskType mark); +bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump, psMaskType maskSat); bool psphotBasicDeblend (psArray *sources, psMetadata *recipe); -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe); +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark); bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); @@ -41,11 +41,11 @@ 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); -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); -bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); -bool psphotReplaceUnfit (psArray *sources); -bool psphotReplaceAll (psArray *sources); -bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); -bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background); +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 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 psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark); bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); @@ -56,6 +56,6 @@ int pmSourceSortBySN (const void **a, const void **b); int pmSourceSortByY (const void **a, const void **b); -bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore); -bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, pmConfig *config); +bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psMaskType maskVal); +bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal); void psphotSourceFreePixels (psArray *sources); @@ -83,12 +83,12 @@ bool psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout); bool psphotInitLimitsEXT (psMetadata *recipe); -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf); -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf); -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source); -psArray *psphotFitDBL (pmReadout *readout, pmSource *source); +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal); +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal); +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal); +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal); // functions to support simultaneous multi-source fitting -bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode); +bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode, psMaskType maskVal); // plotting functions (available if libkapa is installed) @@ -101,16 +101,16 @@ 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); +bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal); bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source); -bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal); bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY); -bool psphotAddWithTest (pmSource *source, bool useState); -bool psphotSubWithTest (pmSource *source, bool useState); -bool psphotSetState (pmSource *source, bool curState); +bool psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal); +bool psphotSubWithTest (pmSource *source, bool useState, psMaskType maskVal); +bool psphotSetState (pmSource *source, bool curState, psMaskType maskVal); bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe); bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf); +bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); #endif Index: /trunk/psphot/src/psphotAddNoise.c =================================================================== --- /trunk/psphot/src/psphotAddNoise.c (revision 13899) +++ /trunk/psphot/src/psphotAddNoise.c (revision 13900) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) { +bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add, psMaskType maskVal) { bool status = false; @@ -34,50 +34,50 @@ // loop over all source for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + pmSource *source = sources->data[i]; - // skip sources which were not subtracted - if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; + // skip sources which were not subtracted + if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; // model must be defined - - if (add) { - psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } else { - psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } + // select appropriate model + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) continue; // model must be defined - psF32 *PAR = model->params->data.F32; + if (add) { + psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + } else { + psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + } - // save original values - float oldI0 = PAR[PM_PAR_I0]; - oldshape.sx = PAR[PM_PAR_SXX]; - oldshape.sy = PAR[PM_PAR_SYY]; - oldshape.sxy = PAR[PM_PAR_SXY]; + psF32 *PAR = model->params->data.F32; - // increase size and height of source - axes = psEllipseShapeToAxes (oldshape, 20.0); - axes.major *= SIZE; - axes.minor *= SIZE; - newshape = psEllipseAxesToShape (axes); - PAR[PM_PAR_I0] = FACTOR*PS_SQR(oldI0); - PAR[PM_PAR_SXX] = newshape.sx; - PAR[PM_PAR_SYY] = newshape.sy; - PAR[PM_PAR_SXY] = newshape.sxy; + // save original values + float oldI0 = PAR[PM_PAR_I0]; + oldshape.sx = PAR[PM_PAR_SXX]; + oldshape.sy = PAR[PM_PAR_SYY]; + oldshape.sxy = PAR[PM_PAR_SXY]; - // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected - pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add); - - // restore original values - PAR[PM_PAR_I0] = oldI0; - PAR[PM_PAR_SXX] = oldshape.sx; - PAR[PM_PAR_SYY] = oldshape.sy; - PAR[PM_PAR_SXY] = oldshape.sxy; + // increase size and height of source + axes = psEllipseShapeToAxes (oldshape, 20.0); + axes.major *= SIZE; + axes.minor *= SIZE; + newshape = psEllipseAxesToShape (axes); + PAR[PM_PAR_I0] = FACTOR*PS_SQR(oldI0); + PAR[PM_PAR_SXX] = newshape.sx; + PAR[PM_PAR_SYY] = newshape.sy; + PAR[PM_PAR_SXY] = newshape.sxy; + + // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected + pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal); + + // restore original values + PAR[PM_PAR_I0] = oldI0; + PAR[PM_PAR_SXX] = oldshape.sx; + PAR[PM_PAR_SYY] = oldshape.sy; + PAR[PM_PAR_SXY] = oldshape.sxy; } if (add) { - psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); } else { - psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); } return true; Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 13899) +++ /trunk/psphot/src/psphotApResid.c (revision 13900) @@ -6,5 +6,5 @@ // measure the aperture residual statistics -bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) +bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark) { int Nfail = 0; @@ -31,7 +31,7 @@ bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); int NSTAR_APERTURE_CORRECTION_MIN = - psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN"); + psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN"); if (!status) { - NSTAR_APERTURE_CORRECTION_MIN = 5; + NSTAR_APERTURE_CORRECTION_MIN = 5; } @@ -47,8 +47,8 @@ psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS); - if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH)) { - psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections"); - psFree(psf->growth); psf->growth = NULL; - return false; + if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, mark)) { + psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections"); + psFree(psf->growth); psf->growth = NULL; + return false; } @@ -80,5 +80,5 @@ // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures // will fail if below S/N threshold or model is missing - if (!pmSourceMagnitudes (source, psf, photMode)) { + if (!pmSourceMagnitudes (source, psf, photMode, maskVal, mark)) { Nskip ++; continue; @@ -135,5 +135,5 @@ if (Npsf < NSTAR_APERTURE_CORRECTION_MIN) { psError(PSPHOT_ERR_APERTURE, true, "Only %d valid aperture residual sources (need %d), giving up", - Npsf, NSTAR_APERTURE_CORRECTION_MIN); + Npsf, NSTAR_APERTURE_CORRECTION_MIN); return false; } @@ -154,5 +154,5 @@ stats->max = 3.0; -#define P_APTREND_SWITCH_CLEANUP /* Cleanup memory on error in ApTrendOption switch */ \ +#define P_APTREND_SWITCH_CLEANUP /* Cleanup memory on error in ApTrendOption switch */ \ psFree(psf->growth); psf->growth = NULL; \ psFree(mask); \ @@ -172,170 +172,170 @@ break; case PM_PSF_APTREND_CONSTANT: - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_SKYBIAS: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_SKYSAT: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_XY_QUAD: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_SKY_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_SKYSAT_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; case PM_PSF_APTREND_ALL: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // fit just SkyBias and clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - // finally, fit x, y, SkyBias and clip out objects which are too far from the median - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); - if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); - P_APTREND_SWITCH_CLEANUP; - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // fit just SkyBias and clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + // finally, fit x, y, SkyBias and clip out objects which are too far from the median + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); + P_APTREND_SWITCH_CLEANUP; + return false; + } + break; default: psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName); @@ -343,5 +343,5 @@ } #undef P_APTREND_SWITCH_CLEANUP - + // construct the fitted values and the residuals psVector *fit = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux); @@ -376,6 +376,6 @@ psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid); - psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", - Nkeep, Npsf, psTimerMark ("psphot")); + psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", + Nkeep, Npsf, psTimerMark ("psphot")); psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n", psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat); Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 13899) +++ /trunk/psphot/src/psphotBlendFit.c (revision 13900) @@ -2,5 +2,5 @@ // XXX I don't like this name -bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { +bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) { int Nfit = 0; @@ -14,5 +14,5 @@ // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, pmSourceSortBySN); - + // S/N limit to perform full non-linear fits float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM"); @@ -30,66 +30,66 @@ for (int i = 0; i < sources->n; i++) { - // if (i%100 == 0) psphotFitSummary (); + // if (i%100 == 0) psphotFitSummary (); - pmSource *source = sources->data[i]; + pmSource *source = sources->data[i]; - // skip non-astronomical objects (very likely defects) - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - // skip DBL second sources (ie, added by psphotFitBlob) - if (source->mode & PM_SOURCE_MODE_PAIR) continue; + // skip DBL second sources (ie, added by psphotFitBlob) + if (source->mode & PM_SOURCE_MODE_PAIR) continue; - // limit selection to some SN limit - // XXX this should use peak? - if (source->moments == NULL) continue; - if (source->moments->SN < FIT_SN_LIM) continue; + // limit selection to some SN limit + // XXX this should use peak? + if (source->moments == NULL) continue; + if (source->moments->SN < FIT_SN_LIM) continue; - // XXX this should use peak? - if (source->moments->x < AnalysisRegion.x0) continue; - if (source->moments->y < AnalysisRegion.y0) continue; - if (source->moments->x > AnalysisRegion.x1) continue; - if (source->moments->y > AnalysisRegion.y1) continue; + // XXX this should use peak? + if (source->moments->x < AnalysisRegion.x0) continue; + if (source->moments->y < AnalysisRegion.y0) continue; + if (source->moments->x > AnalysisRegion.x1) continue; + if (source->moments->y > AnalysisRegion.y1) continue; - // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; - // skip sources which are insignificant flux? - if (source->modelPSF->params->data.F32[1] < 0.1) { - psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", - source->modelPSF->params->data.F32[1], - source->modelPSF->params->data.F32[2], - source->modelPSF->params->data.F32[3]); - continue; - } + // skip sources which are insignificant flux? + if (source->modelPSF->params->data.F32[1] < 0.1) { + psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", + source->modelPSF->params->data.F32[1], + source->modelPSF->params->data.F32[2], + source->modelPSF->params->data.F32[3]); + continue; + } - // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL); - } - Nfit ++; + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Nfit ++; - // try fitting PSFs, then try extended sources - // these functions subtract the resulting fitted source (XXX and update the modelFlux?) - if (psphotFitBlend (readout, source, psf)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); - Npsf ++; - continue; - } - if (psphotFitBlob (readout, source, sources, psf)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); - Next ++; - continue; - } + // try fitting PSFs, then try extended sources + // these functions subtract the resulting fitted source (XXX and update the modelFlux?) + if (psphotFitBlend (readout, source, psf, maskVal)) { + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); + Npsf ++; + continue; + } + if (psphotFitBlob (readout, source, sources, psf, maskVal)) { + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); + Next ++; + continue; + } - psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); - Nfail ++; + psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); + Nfail ++; - // re-subtract the object, leave local sky - pmSourceCacheModel (source); - pmSourceSub (source, PM_MODEL_OP_FULL); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode |= PM_SOURCE_MODE_TEMPSUB; + // re-subtract the object, leave local sky + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->mode |= PM_SOURCE_MODE_TEMPSUB; } Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 13899) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 13900) @@ -2,5 +2,5 @@ // try PSF models and select best option -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { bool status; @@ -43,6 +43,6 @@ psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md); if (!psfTrendMask) { - psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); - return NULL; + psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); + return NULL; } @@ -53,11 +53,11 @@ pmSource *source = sources->data[i]; if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - // keep NSTARS PSF stars, unmark the rest - if (stars->n < NSTARS) { - psArrayAdd (stars, 200, source); - } else { - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - } - } + // keep NSTARS PSF stars, unmark the rest + if (stars->n < NSTARS) { + psArrayAdd (stars, 200, source); + } else { + source->mode &= ~PM_SOURCE_MODE_PSFSTAR; + } + } } psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n); @@ -65,6 +65,6 @@ if (stars->n == 0) { psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates"); - psFree (stars); - psFree (psfTrendMask); + psFree (stars); + psFree (psfTrendMask); return NULL; } @@ -91,5 +91,5 @@ psMetadataItem *item = psListGetAndIncrement (iter); char *modelName = item->data.V; - models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS); + models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark); } @@ -127,13 +127,13 @@ // print/dump psf parameters if (psTraceGetLevel("psphot") >= 5) { - for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) { - psPolynomial2D *poly = try->psf->params_NEW->data[i]; - for (int nx = 0; nx <= poly->nX; nx++) { - for (int ny = 0; ny <= poly->nY; ny++) { - if (poly->mask[nx][ny]) continue; - fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny); - } - } - } + for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) { + psPolynomial2D *poly = try->psf->params_NEW->data[i]; + for (int nx = 0; nx <= poly->nX; nx++) { + for (int ny = 0; ny <= poly->nY; ny++) { + if (poly->mask[nx][ny]) continue; + fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny); + } + } + } } @@ -147,19 +147,19 @@ psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32); for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) - continue; - - pmSource *source = try->sources->data[i]; - Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX]; - Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY]; - dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]; - dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX]; - dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY]; - Sx->n ++; - Sy->n ++; - dSN->n ++; - dSx->n ++; - dSy->n ++; + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX]; + Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY]; + dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]; + dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX]; + dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY]; + Sx->n ++; + Sy->n ++; + dSN->n ++; + dSx->n ++; + dSy->n ++; } psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); @@ -174,6 +174,6 @@ for (int i = 0; i < dSN->n; i++) { - dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]); - dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]); + dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]); + dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]); } @@ -206,79 +206,79 @@ // build a PSF residual image - if (!psphotMakeResiduals (try->sources, recipe, try->psf)) { - psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF"); - psFree (models); - return NULL; + if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) { + psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF"); + psFree (models); + return NULL; } // XXX test dump of psf star data and psf-subtracted image - if (psTraceGetLevel("psphot.psfstars") > 5) { - psphotSaveImage (NULL, readout->image, "rawstars.fits"); - - for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) - continue; - - pmSource *source = try->sources->data[i]; - float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - - // set the mask and subtract the PSF model - // XXX should we be using maskObj? should we be unsetting the mask? - // use pmModelSub because modelFlux has not been generated - assert (source->maskObj); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); - pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); - } - - FILE *f = fopen ("shapes.dat", "w"); - for (int i = 0; i < try->sources->n; i++) { - psF32 inPar[10]; - - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - - pmSource *source = try->sources->data[i]; - psF32 *outPar = source->modelEXT->params->data.F32; - - psEllipseShape shape; - - shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = outPar[PM_PAR_SXY]; - - psEllipsePol pol = pmPSF_ModelToFit (outPar); - inPar[PM_PAR_E0] = pol.e0; - inPar[PM_PAR_E1] = pol.e1; - inPar[PM_PAR_E2] = pol.e2; - pmPSF_FitToModel (inPar, 0.1); - - psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); - psEllipseAxes axes2; - (void)psEllipsePolToAxes(pol, 0.1, &axes2); - psEllipsePol pol2 = psEllipseAxesToPol (axes1); - - fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", - i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], - outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], - pol.e0, pol.e1, pol.e2, - pol2.e0, pol2.e1, pol2.e2, - inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], - axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, - axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD - ); - } - fclose (f); - - psphotSaveImage (NULL, readout->image, "psfstars.fits"); - pmSourcesWritePSFs (try->sources, "psfstars.dat"); - pmSourcesWriteEXTs (try->sources, "extstars.dat", false); - psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); - psMetadataConfigWrite (psfData, "psfmodel.dat"); - psFree (psfData); - psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n"); - exit (0); + if (psTraceGetLevel("psphot.psfstars") > 5) { + psphotSaveImage (NULL, readout->image, "rawstars.fits"); + + for (int i = 0; i < try->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // set the mask and subtract the PSF model + // XXX should we be using maskObj? should we be unsetting the mask? + // use pmModelSub because modelFlux has not been generated + assert (source->maskObj); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); + pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + FILE *f = fopen ("shapes.dat", "w"); + for (int i = 0; i < try->sources->n; i++) { + psF32 inPar[10]; + + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = try->sources->data[i]; + psF32 *outPar = source->modelEXT->params->data.F32; + + psEllipseShape shape; + + shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; + shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; + shape.sxy = outPar[PM_PAR_SXY]; + + psEllipsePol pol = pmPSF_ModelToFit (outPar); + inPar[PM_PAR_E0] = pol.e0; + inPar[PM_PAR_E1] = pol.e1; + inPar[PM_PAR_E2] = pol.e2; + pmPSF_FitToModel (inPar, 0.1); + + psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes2; + (void)psEllipsePolToAxes(pol, 0.1, &axes2); + psEllipsePol pol2 = psEllipseAxesToPol (axes1); + + fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", + i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], + outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], + pol.e0, pol.e1, pol.e2, + pol2.e0, pol2.e1, pol2.e2, + inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], + axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, + axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD + ); + } + fclose (f); + + psphotSaveImage (NULL, readout->image, "psfstars.fits"); + pmSourcesWritePSFs (try->sources, "psfstars.dat"); + pmSourcesWriteEXTs (try->sources, "extstars.dat", false); + psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); + psMetadataConfigWrite (psfData, "psfmodel.dat"); + psFree (psfData); + psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n"); + exit (0); } @@ -290,7 +290,7 @@ if (!psphotPSFstats (readout, recipe, psf)) { - psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); - psFree(psf); - return NULL; + psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); + psFree(psf); + return NULL; } @@ -327,7 +327,7 @@ pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre"); - psFree(modelEXT); - return false; + psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre"); + psFree(modelEXT); + return false; } @@ -346,7 +346,7 @@ psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X", PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y", PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X", PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y", PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true); @@ -386,5 +386,5 @@ moments.xy = source->moments->Sxy; - // limit axis ratio < 20.0 + // limit axis ratio < 20.0 axes = psEllipseMomentsToAxes (moments, 20.0); Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 13899) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 13900) @@ -1,8 +1,8 @@ # include "psphotInternal.h" -// In this function, we smooth the image, then search for the peaks +// In this function, we smooth the image, then search for the peaks psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, - bool returnFootprints, - const int pass) { + bool returnFootprints, + const int pass, psMaskType maskVal) { float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK; @@ -13,57 +13,57 @@ if (pass == 1) { - SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); } else { - bool status_x, status_y; - float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); - float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); - if (!status_x | !status_y) { - psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined"); - return false; - } - SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + bool status_x, status_y; + float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); + float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); + if (!status_x | !status_y) { + psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined"); + return false; + } + SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); } // smooth the image, applying the mask as we go psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_im, readout->mask, 0xff, SIGMA_SMTH, NSIGMA_SMTH); + psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot")); // smooth the weight, applying the mask as we go psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_wt, readout->mask, 0xff, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); + psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot")); psImage *mask = readout->mask; - // optionally save example images under trace + // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - char name[64]; - sprintf (name, "imsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); - sprintf (name, "wtsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_wt, name); + char name[64]; + sprintf (name, "imsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); + sprintf (name, "wtsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_wt, name); } // build the significance image on top of smooth_im for (int j = 0; j < smooth_im->numRows; j++) { - for (int i = 0; i < smooth_im->numCols; i++) { - float value = smooth_im->data.F32[j][i]; - if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || mask->data.U8[j][i]) { - smooth_im->data.F32[j][i] = 0.0; - } else { - smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; - } - } + for (int i = 0; i < smooth_im->numCols; i++) { + float value = smooth_im->data.F32[j][i]; + if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) { + smooth_im->data.F32[j][i] = 0.0; + } else { + smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; + } + } } psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot")); - // optionally save example images under trace + // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - char name[64]; - sprintf (name, "snsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); + char name[64]; + sprintf (name, "snsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); } @@ -73,9 +73,9 @@ // signal/noise limit for the detected peaks if (pass == 1) { - NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); } else { - NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); - } - + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + } + // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing // gaussian SIGMA. a peak in the significance image has an effective S/N for faint sources @@ -87,8 +87,8 @@ // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to // S/N = sqrt(Io)*4*pi*sigma_sm^2 - // thus, the threshold is: + // thus, the threshold is: float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); if (effArea < 1) { - effArea = 1; // never less than a pixel + effArea = 1; // never less than a pixel } float threshold = PS_SQR(NSIGMA_PEAK) / effArea; @@ -97,11 +97,11 @@ psArray *peaks = pmFindImagePeaks (smooth_im, threshold); if (peaks == NULL) { - // XXX this may also be due to a programming or config error - // XXX do we need to set something in the readout->analysis to indicate that - // we tried and failed to find peaks (something in the header data) - psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); - return false; + // XXX this may also be due to a programming or config error + // XXX do we need to set something in the readout->analysis to indicate that + // we tried and failed to find peaks (something in the header data) + psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); + return false; } - + // correct the peak values to S/N = sqrt(value*effArea) // get the peak flux from the unsmoothed image @@ -110,27 +110,27 @@ int col0 = readout->image->col0; for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - peak->SN = sqrt(peak->value*effArea); - peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; + pmPeak *peak = peaks->data[i]; + peak->SN = sqrt(peak->value*effArea); + peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; } // limit the total number of returned peak as specified if (pass == 1) { - psArraySort (peaks, pmPeakSortBySN); - int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); - if (NMAX && (peaks->n > NMAX)) { - psArray *tmpPeaks = psArrayAllocEmpty (NMAX); - for (int i = 0; i < NMAX; i++) { - psArrayAdd (tmpPeaks, 100, peaks->data[i]); - } - psFree (peaks); - peaks = tmpPeaks; - } - } + psArraySort (peaks, pmPeakSortBySN); + int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); + if (NMAX && (peaks->n > NMAX)) { + psArray *tmpPeaks = psArrayAllocEmpty (NMAX); + for (int i = 0; i < NMAX; i++) { + psArrayAdd (tmpPeaks, 100, peaks->data[i]); + } + psFree (peaks); + peaks = tmpPeaks; + } + } // optional dump of all peak data char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE"); if (status && (output != NULL) && (output[0])) { - pmPeaksWriteText (peaks, output); + pmPeaksWriteText (peaks, output); } psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); @@ -142,22 +142,22 @@ // want to do that // - if (returnFootprints) { // We want an array of pmFootprint, not pmPeak - int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); - if (!status) { - npixMin = 1; - } - float FOOTPRINT_NSIGMA_LIMIT = - psMetadataLookupS32(&status, recipe, - (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2"); - if (!status) { - FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK; - } - threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; + if (returnFootprints) { // We want an array of pmFootprint, not pmPeak + int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); + if (!status) { + npixMin = 1; + } + float FOOTPRINT_NSIGMA_LIMIT = + psMetadataLookupS32(&status, recipe, + (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2"); + if (!status) { + FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK; + } + threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; - psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin); - pmPeaksAssignToFootprints(footprints, peaks); + psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin); + pmPeaksAssignToFootprints(footprints, peaks); - psFree(peaks); - peaks = footprints; // well, you know what I mean + psFree(peaks); + peaks = footprints; // well, you know what I mean } @@ -174,20 +174,20 @@ */ 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 +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. + nsigma_delta = 0; // min. } float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); if (!status) { - nsigma_min = 0; + nsigma_min = 0; } const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); return pmFootprintArrayCullPeaks(image, weight, footprints, - nsigma_delta, nsigma_min*skyStdev); + nsigma_delta, nsigma_min*skyStdev); } Index: /trunk/psphot/src/psphotFitSet.c =================================================================== --- /trunk/psphot/src/psphotFitSet.c (revision 13899) +++ /trunk/psphot/src/psphotFitSet.c (revision 13900) @@ -2,5 +2,5 @@ // This is not used in main psphot code (only in psphotModelTest.c) -bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) { +bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode, psMaskType maskVal) { double x, y, Io; @@ -25,5 +25,5 @@ // XXX pmSourceFitSet must cache the modelFlux? - pmSourceFitSet (source, modelSet, mode); + pmSourceFitSet (source, modelSet, mode, maskVal); // write out positive object @@ -33,5 +33,5 @@ for (int i = 0; i < modelSet->n; i++) { pmModel *model = modelSet->data[i]; - pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL); + pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); fprintf (stderr, "output parameters (obj %d):\n", i); Index: /trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinear.c (revision 13899) +++ /trunk/psphot/src/psphotFitSourcesLinear.c (revision 13900) @@ -5,12 +5,12 @@ // model with selected pixels, and the fit radius must be defined -// given the set of sources, each of which points to the pixels in the -// science image, we construct a set of simulated sources with their own pixels. -// these are used to determine the simultaneous linear fit of fluxes. +// given the set of sources, each of which points to the pixels in the +// science image, we construct a set of simulated sources with their own pixels. +// these are used to determine the simultaneous linear fit of fluxes. // the analysis is performed wrt the simulated pixel values 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) { +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal) { bool status; @@ -41,9 +41,9 @@ int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); if (!status) { - SKY_FIT_ORDER = 0; + SKY_FIT_ORDER = 0; } bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR"); if (!status) { - SKY_FIT_LINEAR = false; + SKY_FIT_LINEAR = false; } @@ -55,5 +55,5 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - if (source->type == PM_SOURCE_TYPE_STAR && + if (source->type == PM_SOURCE_TYPE_STAR && source->mode & PM_SOURCE_MODE_SATSTAR) continue; if (final) { @@ -67,5 +67,5 @@ y = source->peak->yf; - // is the source in the region of interest? + // is the source in the region of interest? if (x < AnalysisRegion.x0) continue; if (y < AnalysisRegion.y0) continue; @@ -97,11 +97,11 @@ psSparseMatrixElement (sparse, i, i, f); - // the formal error depends on the weighting scheme - if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false); - errors->data.F32[i] = 1.0 / sqrt(var); - } else { - errors->data.F32[i] = 1.0 / sqrt(f); - } + // the formal error depends on the weighting scheme + if (CONSTANT_PHOTOMETRIC_WEIGHTS) { + float var = pmSourceModelDotModel (SRCi, SRCi, false); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(f); + } @@ -110,21 +110,21 @@ psSparseVectorElement (sparse, i, f); - // add the per-source weights (border region) - switch (SKY_FIT_ORDER) { - case 1: - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 2, f); - - case 0: - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 0, f); - break; - - default: - psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER); - break; - } + // add the per-source weights (border region) + switch (SKY_FIT_ORDER) { + case 1: + f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 1, f); + f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 2, f); + + case 0: + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 0, f); + break; + + default: + psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER); + break; + } // loop over all other stars following this one @@ -159,9 +159,9 @@ psVector *skyfit = NULL; if (SKY_FIT_LINEAR) { - psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); - fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); + psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); + fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); } else { - norm = psSparseSolve (NULL, constraint, sparse, 5); - skyfit = NULL; + norm = psSparseSolve (NULL, constraint, sparse, 5); + skyfit = NULL; } @@ -169,5 +169,5 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); + pmModel *model = pmSourceGetModel (NULL, source); // assign linearly-fitted normalization @@ -177,11 +177,11 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - // XXX is the value of 'errors' modified by the sky fit? + // XXX is the value of 'errors' modified by the sky fit? // subtract object - pmSourceSub (source, PM_MODEL_OP_FULL); + 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 + // XXX not sure about the use of TEMPSUB } @@ -190,5 +190,5 @@ pmSource *source = fitSources->data[i]; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->weight); + pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal); } @@ -219,10 +219,10 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) continue; float x = model->params->data.F32[PM_PAR_XPOS]; float y = model->params->data.F32[PM_PAR_YPOS]; - psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); - } + psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + } // accumulate the image statistics from the masked regions @@ -233,31 +233,31 @@ double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; w = x = y = x2 = xy = y2 = fo = fx = fy = 0; - + int col0 = readout->image->col0; int row0 = readout->image->row0; for (int j = 0; j < readout->image->numRows; j++) { - for (int i = 0; i < readout->image->numCols; i++) { - if (mask[j][i]) continue; - if (constant_weights) { - wt = 1.0; - } else { - wt = weight[j][i]; - } - f = image[j][i]; - w += 1/wt; - fo += f/wt; - if (SKY_FIT_ORDER == 0) continue; - - xc = i + col0; - yc = j + row0; - x += xc/wt; - y += yc/wt; - x2 += xc*xc/wt; - xy += xc*yc/wt; - y2 += yc*yc/wt; - fx += f*xc/wt; - fy += f*yc/wt; - } + for (int i = 0; i < readout->image->numCols; i++) { + if (mask[j][i]) continue; + if (constant_weights) { + wt = 1.0; + } else { + wt = weight[j][i]; + } + f = image[j][i]; + w += 1/wt; + fo += f/wt; + if (SKY_FIT_ORDER == 0) continue; + + xc = i + col0; + yc = j + row0; + x += xc/wt; + y += yc/wt; + x2 += xc*xc/wt; + xy += xc*yc/wt; + y2 += yc*yc/wt; + fx += f*xc/wt; + fy += f*yc/wt; + } } @@ -269,15 +269,15 @@ psSparseBorderElementT (border, 0, 0, w); if (SKY_FIT_ORDER > 0) { - psSparseBorderElementG (border, 0, fx); - psSparseBorderElementG (border, 0, fy); - psSparseBorderElementT (border, 1, 0, x); - psSparseBorderElementT (border, 2, 0, y); - psSparseBorderElementT (border, 0, 1, x); - psSparseBorderElementT (border, 1, 1, x2); - psSparseBorderElementT (border, 2, 1, xy); - psSparseBorderElementT (border, 0, 2, y); - psSparseBorderElementT (border, 1, 2, xy); - psSparseBorderElementT (border, 2, 2, y2); - } + psSparseBorderElementG (border, 0, fx); + psSparseBorderElementG (border, 0, fy); + psSparseBorderElementT (border, 1, 0, x); + psSparseBorderElementT (border, 2, 0, y); + psSparseBorderElementT (border, 0, 1, x); + psSparseBorderElementT (border, 1, 1, x2); + psSparseBorderElementT (border, 2, 1, xy); + psSparseBorderElementT (border, 0, 2, y); + psSparseBorderElementT (border, 1, 2, xy); + psSparseBorderElementT (border, 2, 2, y2); + } return true; } Index: /trunk/psphot/src/psphotGrowthCurve.c =================================================================== --- /trunk/psphot/src/psphotGrowthCurve.c (revision 13899) +++ /trunk/psphot/src/psphotGrowthCurve.c (revision 13900) @@ -7,5 +7,5 @@ // XXX add a option to turn off the curve-of-growth (ie, make the apMag = fitMag everywhere); -bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore) { +bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psMaskType maskVal) { // bool status; @@ -45,5 +45,5 @@ // loop over a range of source fluxes - // no need to interpolate since we have force the object center + // no need to interpolate since we have force the object center // to 0.5, 0.5 above for (int i = 0; i < psf->growth->radius->n; i++) { @@ -53,16 +53,16 @@ radius = psf->growth->radius->data.F32[i]; - // NOTE: we use pmModelAdd not pmSourceAdd because we are not working with a normal source + // NOTE: we use pmModelAdd not pmSourceAdd because we are not working with a normal source psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); - pmModelAdd (image, mask, model, PM_MODEL_OP_FULL); - pmSourcePhotometryAper (&apMag, model, image, mask); + pmModelAdd (image, mask, model, PM_MODEL_OP_FULL, maskVal); + pmSourcePhotometryAper (&apMag, model, image, mask, maskVal); psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); - // the 'ignore' mode is for testing - if (ignore) { - psf->growth->apMag->data.F32[i] = fitMag; - } else { - psf->growth->apMag->data.F32[i] = apMag; - } + // the 'ignore' mode is for testing + if (ignore) { + psf->growth->apMag->data.F32[i] = fitMag; + } else { + psf->growth->apMag->data.F32[i] = apMag; + } } Index: /trunk/psphot/src/psphotGuessModels.c =================================================================== --- /trunk/psphot/src/psphotGuessModels.c (revision 13899) +++ /trunk/psphot/src/psphotGuessModels.c (revision 13900) @@ -17,5 +17,5 @@ // construct an initial PSF model for each object -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) { psTimerStart ("psphot"); @@ -57,21 +57,21 @@ 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; - } + 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; + source->mode |= PM_SOURCE_MODE_BADPSF; } psFree (modelEXT); @@ -85,5 +85,5 @@ source->modelPSF->residuals = psf->residuals; - pmSourceCacheModel (source); + pmSourceCacheModel (source, maskVal); } psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); Index: /trunk/psphot/src/psphotImageMedian.c =================================================================== --- /trunk/psphot/src/psphotImageMedian.c (revision 13899) +++ /trunk/psphot/src/psphotImageMedian.c (revision 13900) @@ -40,5 +40,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotImageMedian (pmConfig *config, const pmFPAview *view) +bool psphotImageMedian (pmConfig *config, const pmFPAview *view, psMaskType maskVal) { bool status = true; @@ -195,5 +195,5 @@ // XXX don't bother trying if there are no valid pixels... - if (psImageBackground(stats, subset, submask, 0xff, rng)) { + if (psImageBackground(stats, subset, submask, maskVal, rng)) { if (stats->options & PS_STAT_ROBUST_QUARTILE) { modelData[iy][ix] = stats->robustMedian; @@ -205,5 +205,5 @@ psStatsOptions currentOptions = stats->options; stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; - if (!psImageBackground(stats, subset, submask, 0xff, rng)) { + if (!psImageBackground(stats, subset, submask, maskVal, rng)) { psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " "(%dx%d, (row0,col0) = (%d,%d)", Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 13899) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 13900) @@ -2,7 +2,9 @@ bool psphotMagnitudes(psArray *sources, - psMetadata *recipe, - pmPSF *psf, - pmReadout *background) + psMetadata *recipe, + pmPSF *psf, + pmReadout *background, + psMaskType maskVal, + psMaskType mark) { bool status = false; @@ -13,9 +15,9 @@ pmSourceMagnitudesInit (recipe); - // XXX require that we have a background model, or + // XXX require that we have a background model, or // allow it to be missing, setting local sky to 0.0? PS_ASSERT (background, false); - // the binning details are saved on the analysis metadata + // the binning details are saved on the analysis metadata psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING"); PS_ASSERT (status, false); @@ -29,15 +31,15 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - status = pmSourceMagnitudes (source, psf, photMode); - if (status) Nap ++; + pmSource *source = (pmSource *) sources->data[i]; + status = pmSourceMagnitudes (source, psf, photMode, maskVal, mark); + if (status) Nap ++; - source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, background->image, binning); - if (isnan(source->sky) && false) { - psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky"); - psErrorStackPrint(NULL, " "); - psErrorClear(); - } - } + source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, background->image, binning); + if (isnan(source->sky) && false) { + psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky"); + psErrorStackPrint(NULL, " "); + psErrorClear(); + } + } psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot"), sources->n, Nap); Index: /trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- /trunk/psphot/src/psphotMakeResiduals.c (revision 13899) +++ /trunk/psphot/src/psphotMakeResiduals.c (revision 13900) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) { +bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) { bool status, isPSF; @@ -94,5 +94,5 @@ psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32); psImage *mask = psImageCopy (NULL, source->maskView, PS_TYPE_U8); - pmModelSub (image, mask, model, PM_MODEL_OP_FUNC); + pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); // re-normalize image and weight Index: /trunk/psphot/src/psphotMaskReadout.c =================================================================== --- /trunk/psphot/src/psphotMaskReadout.c (revision 13899) +++ /trunk/psphot/src/psphotMaskReadout.c (revision 13900) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, pmConfig *config) { +bool psphotMaskReadout (pmReadout *readout, psMetadata *recipe, psMaskType maskVal) { bool status; @@ -19,5 +19,5 @@ // psImageKeepRegion assumes the region refers to the parent coordinates - psImageKeepRegion (readout->mask, keep, "OR", pmConfigMask("BAD", config)); + psImageKeepRegion (readout->mask, keep, "OR", maskVal); return true; Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 13899) +++ /trunk/psphot/src/psphotModelTest.c (revision 13900) @@ -3,5 +3,5 @@ // XXX consider this function : add more test information? -bool psphotModelTest (pmReadout *readout, psMetadata *recipe) { +bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { bool status; @@ -94,5 +94,5 @@ // find the local sky - status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); if (!status) psAbort("pmSourceLocalSky error"); @@ -159,17 +159,17 @@ // define the pixels used for the fit - psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", PM_MASK_MARK); + psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", mark); char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); if (status) { - status = psphotFitSet (source, model, fitset, fitMode); + status = psphotFitSet (source, model, fitset, fitMode, maskVal); exit (0); } - status = pmSourceFitModel (source, model, fitMode); + status = pmSourceFitModel (source, model, fitMode, maskVal); // measure the source mags pmSourcePhotometryModel (&fitMag, model); - pmSourcePhotometryAper (&obsMag, model, source->pixels, source->maskObj); + pmSourcePhotometryAper (&obsMag, model, source->pixels, source->maskObj, maskVal); fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag); @@ -178,5 +178,5 @@ // subtract object, leave local sky - pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL); + pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); fprintf (stderr, "output parameters: \n"); Index: /trunk/psphot/src/psphotMosaicChip.c =================================================================== --- /trunk/psphot/src/psphotMosaicChip.c (revision 13899) +++ /trunk/psphot/src/psphotMosaicChip.c (revision 13900) @@ -28,5 +28,5 @@ psTrace("pmChipMosaic", 5, "mosaic chip %s to %s (xbin,ybin: %d,%d to %d,%d)\n", in->name, out->name, in->xBin, in->yBin, out->xBin, out->yBin); - status = pmChipMosaic(outChip, inChip, true); + status = pmChipMosaic(outChip, inChip, true, pmConfigMask("BLANK", config)); return status; } Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 13899) +++ /trunk/psphot/src/psphotReadout.c (revision 13900) @@ -12,4 +12,15 @@ } + // 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"); @@ -24,17 +35,15 @@ // generate mask & weight images if they don't already exit - psMaskType satMask = pmConfigMask("SAT", config); - psMaskType badMask = pmConfigMask("BAD", config); if (!readout->mask) { - if (!pmReadoutGenerateMask(readout, satMask, badMask)) { - psError (PSPHOT_ERR_CONFIG, false, "trouble creating mask"); - return false; - } + 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; - } + if (!pmReadoutGenerateWeight(readout, true)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight"); + return false; + } } @@ -46,8 +55,8 @@ // I have a valid mask, now mask in the analysis region of interest - psphotMaskReadout (readout, recipe, config); + psphotMaskReadout (readout, recipe, maskBad); // run a single-model test if desired - psphotModelTest (readout, recipe); + psphotModelTest (readout, recipe, maskVal, maskMark); if (psTraceGetLevel("psphot") >= 5) { @@ -62,5 +71,5 @@ // generate a background model (median, smoothed image) - if (!psphotImageMedian (config, view)) { + if (!psphotImageMedian (config, view, maskVal)) { return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); } @@ -74,5 +83,5 @@ psArray *footprints = NULL; if (useFootprints) { - footprints = psphotFindPeaks (readout, recipe, useFootprints, 1); + footprints = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal); int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS"); if (growRadius > 0) { @@ -85,5 +94,5 @@ peaks = pmFootprintArrayToPeaks(footprints); } else { - peaks = psphotFindPeaks (readout, recipe, useFootprints, 1); + peaks = psphotFindPeaks (readout, recipe, useFootprints, 1, maskVal); } @@ -97,5 +106,5 @@ // construct sources and measure basic stats - psArray *sources = psphotSourceStats (readout, recipe, peaks); + psArray *sources = psphotSourceStats (readout, recipe, peaks, maskVal, maskMark); if (!sources) return false; psFree (peaks); @@ -114,5 +123,5 @@ // classify sources based on moments, brightness - if (!psphotRoughClass (sources, recipe, true)) { + if (!psphotRoughClass (sources, recipe, true, maskSat)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); return psphotReadoutCleanup (config, readout, recipe, NULL, sources); @@ -126,5 +135,5 @@ if (!psf) { // use bright stellar objects to measure PSF - psf = psphotChoosePSF (readout, sources, recipe); + psf = psphotChoosePSF (readout, sources, recipe, maskVal, maskMark); if (psf == NULL) { psLogMsg ("psphot", 3, "failure to construct a psf model"); @@ -143,5 +152,5 @@ // construct an initial model for each object - psphotGuessModels (readout, sources, recipe, psf); + psphotGuessModels (readout, sources, recipe, psf, maskVal); // XXX test output of models @@ -151,5 +160,5 @@ // linear PSF fit to peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE, maskVal); if (dump) psphotSaveImage (NULL, readout->image, "image.v1.fits"); @@ -162,13 +171,13 @@ // non-linear PSF and EXT fit to brighter sources - psphotBlendFit (readout, sources, recipe, psf); + psphotBlendFit (readout, sources, recipe, psf, maskVal); if (dump) psphotSaveImage (NULL, readout->image, "image.v2.fits"); // replace all sources - psphotReplaceAll (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); + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal); if (dump) psphotSaveImage (NULL, readout->image, "image.v4.fits"); if (!strcasecmp (breakPt, "PASS1")) { @@ -183,10 +192,10 @@ // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe, true); + psphotAddNoise (readout, sources, recipe, true, maskVal); // find the peaks in the image psArray *newPeaks; if (useFootprints) { - psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2); + psArray *newFootprints = psphotFindPeaks (readout, recipe, useFootprints, 2, maskVal); int growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2"); @@ -209,22 +218,22 @@ psFree(mergedFootprints); } else { - newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2); + newPeaks = psphotFindPeaks(readout, recipe, useFootprints, 2, maskVal); } // remove noise for subtracted objects - psphotAddNoise (readout, sources, recipe, false); + psphotAddNoise (readout, sources, recipe, false, maskVal); // define new sources based on the new peaks - psArray *newSources = psphotSourceStats (readout, recipe, newPeaks); + psArray *newSources = psphotSourceStats (readout, recipe, newPeaks, maskVal, maskMark); psFree (newPeaks); // set source type - psphotRoughClass (newSources, recipe, false); + psphotRoughClass (newSources, recipe, false, maskSat); // create full input models - psphotGuessModels (readout, newSources, recipe, psf); + psphotGuessModels (readout, newSources, recipe, psf, maskVal); // replace all sources - psphotReplaceAll (sources); + psphotReplaceAll (sources, maskVal); if (dump) psphotSaveImage (NULL, readout->image, "image.v5.fits"); @@ -234,5 +243,5 @@ // linear PSF fit to remaining peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE, maskVal); if (dump) psphotSaveImage (NULL, readout->image, "image.v6.fits"); @@ -243,5 +252,5 @@ // measure aperture photometry corrections - if (!psphotApResid (readout, sources, recipe, psf)) { + if (!psphotApResid (readout, sources, recipe, psf, maskVal, maskMark)) { psTrace ("psphot", 4, "failure on psphotApResid"); psError(PSPHOT_ERR_PHOTOM, false, "Measure aperture photometry corrections"); @@ -251,5 +260,5 @@ // calculate source magnitudes pmReadout *background = psphotSelectBackground (config, view, false); - psphotMagnitudes(sources, recipe, psf, background); + psphotMagnitudes(sources, recipe, psf, background, maskVal, maskMark); // replace failed sources? Index: /trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- /trunk/psphot/src/psphotReplaceUnfit.c (revision 13899) +++ /trunk/psphot/src/psphotReplaceUnfit.c (revision 13900) @@ -2,5 +2,5 @@ // replace the flux for sources which failed -bool psphotReplaceUnfit (psArray *sources) { +bool psphotReplaceUnfit (psArray *sources, psMaskType maskVal) { pmSource *source; @@ -16,7 +16,7 @@ replace: - pmSourceAdd (source, PM_MODEL_OP_FULL); - source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; - source->mode &= ~PM_SOURCE_MODE_TEMPSUB; + 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,5 +24,5 @@ } -bool psphotReplaceAll (psArray *sources) { +bool psphotReplaceAll (psArray *sources, psMaskType maskVal) { pmSource *source; @@ -36,5 +36,5 @@ if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - pmSourceAdd (source, PM_MODEL_OP_FULL); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; } @@ -44,13 +44,13 @@ // add source, if the source has been subtracted (or if we ignore the state) -bool psphotAddWithTest (pmSource *source, bool useState) { +bool psphotAddWithTest (pmSource *source, bool useState, psMaskType maskVal) { // what is current state? (true : add; false : sub) bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); if (state && useState) return true; - + // replace the model if 1) state says it is missing or 2) useState is false (just do it) if (!state || !useState) { - pmSourceAdd (source, PM_MODEL_OP_FULL); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } return true; @@ -58,5 +58,5 @@ // sub source, if the source has been added (or if we ignore the state) -bool psphotSubWithTest (pmSource *source, bool useState) { +bool psphotSubWithTest (pmSource *source, bool useState, psMaskType maskVal) { // what is current state? (true : sub; false : add) @@ -66,21 +66,21 @@ // replace the model if 1) state says it is missing or 2) useState is false (just do it) if (!state || !useState) { - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } return true; } -// add or sub source replace or if the source has -bool psphotSetState (pmSource *source, bool curState) { +// add or sub source replace or if the source has +bool psphotSetState (pmSource *source, bool curState, psMaskType maskVal) { // what is desired state? (true : add; false : sub) bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); if (curState == newState) return true; - + if (curState && !newState) { - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } if (newState && !curState) { - pmSourceAdd (source, PM_MODEL_OP_FULL); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } return true; Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 13899) +++ /trunk/psphot/src/psphotRoughClass.c (revision 13900) @@ -2,5 +2,5 @@ // 2006.02.02 : no leaks -bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump) { +bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool findPsfClump, psMaskType maskSat) { static pmPSFClump psfClump; @@ -10,19 +10,19 @@ if (findPsfClump) { - psfClump = pmSourcePSFClump (sources, recipe); + psfClump = pmSourcePSFClump (sources, recipe); } else if (!havePsfClump) { - psError(PSPHOT_ERR_PROG, false, "You must find the PSF clump before reusing it"); + psError(PSPHOT_ERR_PROG, false, "You must find the PSF clump before reusing it"); } else { - ; + ; } - - havePsfClump = false; // we don't know if it's valid + + havePsfClump = false; // we don't know if it's valid if (psfClump.X < 0) { - psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump"); - return false; + psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump"); + return false; } if (!psfClump.X || !psfClump.Y) { - psError(PSPHOT_ERR_DATA, true, "Failed to find a valid PSF clump"); - return false; + psError(PSPHOT_ERR_DATA, true, "Failed to find a valid PSF clump"); + return false; } psLogMsg ("psphot", 3, "psf clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); @@ -30,7 +30,7 @@ // group into STAR, COSMIC, EXTENDED, SATURATED, etc. - if (!pmSourceRoughClass (sources, recipe, psfClump)) { - psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass"); - return false; + if (!pmSourceRoughClass (sources, recipe, psfClump, maskSat)) { + psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass"); + return false; } @@ -40,6 +40,6 @@ psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot")); - havePsfClump = true; // we have set psfClump - + havePsfClump = true; // we have set psfClump + return true; } Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 13899) +++ /trunk/psphot/src/psphotSourceFits.c (revision 13900) @@ -16,10 +16,10 @@ bool psphotFitSummary () { - fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", - NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); - return true; -} - -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { + fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", + NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); + return true; +} + +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal) { float x, y, dR; @@ -27,5 +27,5 @@ // if this source is not a possible blend, just fit as PSF if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { - bool status = psphotFitPSF (readout, source, psf); + bool status = psphotFitPSF (readout, source, psf, maskVal); return status; } @@ -65,6 +65,6 @@ model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf; - // these should never be invalid values - // XXX drop these tests eventually + // these should never be invalid values + // XXX drop these tests eventually if (isnan(model->params->data.F32[PM_PAR_I0])) psAbort("nan in blend fit"); if (isnan(model->params->data.F32[PM_PAR_XPOS])) psAbort("nan in blend fit"); @@ -83,5 +83,5 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); + pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); // correct model chisq for flux trend @@ -109,7 +109,7 @@ psTrace ("psphot", 5, "fitted blend as PSF\n"); - // build cached model and subtract - pmSourceCacheModel (blend); - pmSourceSub (blend, PM_MODEL_OP_FULL); + // build cached model and subtract + pmSourceCacheModel (blend, maskVal); + pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal); blend->mode |= PM_SOURCE_MODE_SUBTRACTED; blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -121,6 +121,6 @@ psTrace ("psphot", 4, "failed on blend 0 of %ld\n", modelSet->n); psFree (PSF); - psFree (modelSet); - psFree (sourceSet); + psFree (modelSet); + psFree (sourceSet); return false; } @@ -134,6 +134,6 @@ // build cached model and subtract - pmSourceCacheModel (source); - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -141,5 +141,5 @@ } -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal) { double chiTrend; @@ -155,5 +155,5 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF); + pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal); // correct model chisq for flux trend @@ -173,6 +173,6 @@ // build cached model and subtract - pmSourceCacheModel (source); - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -201,5 +201,5 @@ } -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal) { bool okEXT, okDBL; @@ -222,9 +222,9 @@ pmSource *tmpSrc = pmSourceAlloc (); - pmModel *EXT = psphotFitEXT (readout, source); + pmModel *EXT = psphotFitEXT (readout, source, maskVal); okEXT = psphotEvalEXT (tmpSrc, EXT); chiEXT = EXT->chisq / EXT->nDOF; - psArray *DBL = psphotFitDBL (readout, source); + psArray *DBL = psphotFitDBL (readout, source, maskVal); okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); @@ -269,6 +269,6 @@ // build cached model and subtract - pmSourceCacheModel (source); - pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); @@ -294,8 +294,8 @@ // build cached models and subtract - pmSourceCacheModel (source); - pmSourceSub (source, PM_MODEL_OP_FULL); - pmSourceCacheModel (newSrc); - pmSourceSub (newSrc, PM_MODEL_OP_FULL); + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceCacheModel (newSrc, maskVal); + pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal); psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]); @@ -307,5 +307,5 @@ // fit a double PSF source to an extended blob -psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal) { float dx, dy; @@ -348,9 +348,9 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); + pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); return (modelSet); } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal) { NfitEXT ++; @@ -359,5 +359,5 @@ pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); PS_ASSERT (EXT, NULL); - + // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit"); @@ -369,5 +369,5 @@ // fit EXT (not PSF) model (set/unset the pixel mask) - pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT); + pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal); return (EXT); } Index: /trunk/psphot/src/psphotSourcePlots.c =================================================================== --- /trunk/psphot/src/psphotSourcePlots.c (revision 13899) +++ /trunk/psphot/src/psphotSourcePlots.c (revision 13900) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal) { // bool status = false; @@ -17,8 +17,8 @@ // counters to track the size of the image and area used in a row - int dX = 0; // starting corner of next box - int dY = 0; // height of row so far - int NX = 20*DX; // full width of output image - int NY = 0; // total height of output image + int dX = 0; // starting corner of next box + int dY = 0; // height of row so far + int NX = 20*DX; // full width of output image + int NY = 0; // total height of output image // first, examine the PSF and SAT stars: @@ -28,31 +28,31 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); - if (!keep) continue; + if (!keep) continue; - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } @@ -61,11 +61,11 @@ psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32); - int Xo = 0; // starting corner of next box - int Yo = 0; // starting corner of next box - dY = 0; // height of row so far + int Xo = 0; // starting corner of next box + int Yo = 0; // starting corner of next box + dY = 0; // height of row so far int nPSF = 0; int nSAT = 0; - int kapa = 0; // file descriptor for plotting routine + int kapa = 0; // file descriptor for plotting routine // first, examine the PSF and SAT stars: @@ -76,62 +76,62 @@ pmSource *source = sources->data[i]; - bool keep = false; + bool keep = false; if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - keep = true; - } + nPSF ++; + keep = true; + } if (source->mode & PM_SOURCE_MODE_SATSTAR) { - nSAT ++; - keep = true; - } - if (!keep) continue; + nSAT ++; + keep = true; + } + if (!keep) continue; - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - psphotAddWithTest (source, true); // replace source if subtracted - psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotRadialPlot (&kapa, "radial.plots.ps", source); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); - psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); - psphotSetState (source, false); // replace source (has been subtracted) - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - psphotAddWithTest (source, true); // replace source if subtracted - psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotRadialPlot (&kapa, "radial.plots.ps", source); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); - psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); - psphotSetState (source, false); // replace source (has been subtracted) + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) - Xo = DX; - dY = DY; - } - } else { - // extend this row - psphotAddWithTest (source, true); // replace source if subtracted - psphotRadialPlot (&kapa, "radial.plots.ps", source); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); + Xo = DX; + dY = DY; + } + } else { + // extend this row + psphotAddWithTest (source, true, maskVal); // replace source if subtracted + psphotRadialPlot (&kapa, "radial.plots.ps", source); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY); - psphotSubWithTest (source, false); // remove source (force) - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); - psphotSetState (source, false); // replace source (has been subtracted) + psphotSubWithTest (source, false, maskVal); // remove source (force) + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY); + psphotSetState (source, false, maskVal); // replace source (has been subtracted) - Xo += DX; - dY = PS_MAX (dY, DY); - } + Xo += DX; + dY = PS_MAX (dY, DY); + } } Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 13899) +++ /trunk/psphot/src/psphotSourceStats.c (revision 13900) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks) +psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks, psMaskType maskVal, psMaskType mark) { bool status = false; @@ -40,16 +40,16 @@ // skip faint sources - if (source->peak->SN < MIN_SN) { + if (source->peak->SN < MIN_SN) { psArrayAdd (sources, 100, source); psFree (source); - continue; - } + continue; + } // measure a local sky value // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); if (!status) { psFree (source); - Nfail ++; + Nfail ++; continue; } @@ -57,8 +57,8 @@ // measure the local sky variance (needed if noise is not sqrt(signal)) // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER); + status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); if (!status) { psFree (source); - Nfail ++; + Nfail ++; continue; } @@ -70,5 +70,5 @@ psArrayAdd (sources, 100, source); psFree (source); - Nmoments ++; + Nmoments ++; continue; } @@ -83,10 +83,10 @@ psArrayAdd (sources, 100, source); psFree (source); - Nmoments ++; + Nmoments ++; continue; } psFree (source); - Nfail ++; + Nfail ++; continue; }