Index: branches/meh_branches/ppstack_test/psphot/src/Makefile.am =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/Makefile.am (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/Makefile.am (revision 34041) @@ -206,5 +206,6 @@ psphotPetrosianStats.c \ psphotPetrosianVisual.c \ - psphotEfficiency.c + psphotEfficiency.c \ + psphotSetNFrames.c # re-instate these Index: branches/meh_branches/ppstack_test/psphot/src/psphot.h =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphot.h (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphot.h (revision 34041) @@ -99,5 +99,10 @@ bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final); + +# if (HAVE_MODEL_VAR) +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode); +# else bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final); +# endif bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize); @@ -189,5 +194,5 @@ // used by psphotFindDetections pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); -psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks); bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); psErrorCode psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints); @@ -472,5 +477,5 @@ bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); +bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index); bool psphotKronIterate_Threaded (psThreadJob *job); @@ -481,3 +486,6 @@ bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); +bool psphotSetNFrames (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + #endif Index: branches/meh_branches/ppstack_test/psphot/src/psphotAddNoise.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotAddNoise.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotAddNoise.c (revision 34041) @@ -47,5 +47,7 @@ psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); + //psAssert (sources, "missing sources?"); + // if no work, should just return true + if (!sources) return true; psTimerStart ("psphot.noise"); Index: branches/meh_branches/ppstack_test/psphot/src/psphotArguments.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotArguments.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotArguments.c (revision 34041) @@ -217,4 +217,5 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "EXPNUM", "-expnum", "-expnumlist"); if (argc == 1) { Index: branches/meh_branches/ppstack_test/psphot/src/psphotDefineFiles.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotDefineFiles.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotDefineFiles.c (revision 34041) @@ -85,4 +85,26 @@ output->save = true; } +// // optionally save the smoothed variance model (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL"); +// return false; +// } +// output->save = true; +// } +// // optionally save the smoothed variance model's standard deviation (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.VARMDL.STDEV")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL.STDEV"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL.STDEV"); +// return false; +// } +// output->save = true; +// } // optionally save the PSF Model if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) { Index: branches/meh_branches/ppstack_test/psphot/src/psphotEfficiency.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotEfficiency.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotEfficiency.c (revision 34041) @@ -420,5 +420,9 @@ // psphotFitSourcesLinearReadout subtracts the model fits +# if (HAVE_MODEL_VAR) + if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) { +# else if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) { +# endif psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); Index: branches/meh_branches/ppstack_test/psphot/src/psphotFindDetections.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotFindDetections.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotFindDetections.c (revision 34041) @@ -48,4 +48,5 @@ // Use the new pmFootprints approach? const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS"); + const bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -64,4 +65,5 @@ } + bool replaceSourcesForFootprints = false; if (firstPass) { pass = 1; @@ -70,4 +72,5 @@ } else { pass = 2; + replaceSourcesForFootprints = footprintUseUnsubtracted; NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL); NMAX = 0; // unlimited number of peaks in final pass: allow a limit (PEAKS_NMAX_2) ? @@ -96,5 +99,6 @@ // detect the peaks in the significance image - detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); + int totalPeaks = 0; + detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks); psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); if (!detections->peaks) { @@ -105,8 +109,26 @@ return false; } + // hard limit on number of peaks we will accept. (To avoid memory overload in psphotStack) + int maxPeaks = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX_TOTAL"); PS_ASSERT (status, NULL); + if (maxPeaks && (totalPeaks > maxPeaks)) { + psFree (detections); + psError (PSPHOT_ERR_DATA, true, "Too many peaks %d found PEAKS_NMAX_TOTAL: %d", totalPeaks, maxPeaks); + return false; + } // optionally merge peaks into footprints if (useFootprints) { + if (replaceSourcesForFootprints) { + psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false); + psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false); + psFree (significance); + significance = psphotSignificanceImage (readout, recipe, maskVal); + } + psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal); + + if (replaceSourcesForFootprints) { + psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false); + } } Index: branches/meh_branches/ppstack_test/psphot/src/psphotFindPeaks.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotFindPeaks.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotFindPeaks.c (revision 34041) @@ -4,5 +4,5 @@ // image must be constructed to represent (S/N)^2. If nMax is non-zero, only return a maximum // of nMax peaks -psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) { +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) { bool status = false; @@ -18,4 +18,8 @@ psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); return NULL; + } + // return the total number of peaks found before the nMax limit is applied + if (totalPeaks) { + *totalPeaks = peaks->n; } Index: branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c (revision 34041) @@ -12,4 +12,10 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); +# if (HAVE_MODEL_VAR) +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources); +pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe); +bool psphotFreeModelVariance (pmReadout *readout, psArray *sources); +# endif + // for now, let's store the detections on the readout->analysis for each readout bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final) @@ -24,4 +30,16 @@ assert (recipe); +# if (HAVE_MODEL_VAR) + pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe); + if (!fitVarMode) { + psError (PSPHOT_ERR_CONFIG, true, "failed to get LINEAR_FIT_VARIANCE_MODE"); + return false; + } + // MODEL_VAR requires 2 passes -- in the first, we get the rough fluxes; in the second, we + // use the flux to define the model variance before fitting the objects. Other modes only + // do a single pass. + pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode; +# endif + int num = psphotFileruleCount(config, filerule); @@ -50,8 +68,41 @@ psAssert (psf, "missing psf?"); - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { +# if (HAVE_MODEL_VAR) + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) +# else + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) +# endif + { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; } + +# if (HAVE_MODEL_VAR) + // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance + // after we have determined the initial set of fits, then we can generate the variance image and + // re-run the fit against that variance. + if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) { + // generate the model variance image & source pointers + if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + + // replace all sources (use TMPF_SUBTRACTED as test) + psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, false); + + // rerun fit with correct fitVarMode + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + + // free the model variance image & source pointers + if (!psphotFreeModelVariance (readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + } +# endif psphotVisualShowResidualImage (readout, (num > 0)); @@ -62,6 +113,41 @@ } -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) { - +# if (HAVE_MODEL_VAR) +// look up the fit variance mode from the recipe; older recipes do not have the value +// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and +// set the value to either CONST or IMAGE_VAR +pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) { + + bool status = false; + + char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE"); + if (!status) { + bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + if (!status) { + psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS"); + } + pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR; + return fitVarMode; + } + if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) { + return PM_SOURCE_PHOTFIT_CONST; + } + if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) { + return PM_SOURCE_PHOTFIT_IMAGE_VAR; + } + if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) { + return PM_SOURCE_PHOTFIT_MODEL_VAR; + } + psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString); + return PM_SOURCE_PHOTFIT_NONE; +} +# endif + +# if (HAVE_MODEL_VAR) +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) +# else +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) +# endif +{ bool status; float x; @@ -99,4 +185,5 @@ if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); +# if (!HAVE_MODEL_VAR) bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); @@ -104,4 +191,5 @@ psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); } +# endif int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); if (!status) { @@ -231,4 +319,14 @@ psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder); +# if (HAVE_MODEL_VAR) + // if fitVarMode is MODEL_VAR, then we need to generate the model image variance + // XXX we have two possibilities here: + + // 1) do 2 passes, where in the first case we use the CONST weighting, and in the second + // use the fitted model values to define the model + + // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?) +# endif + // fill out the sparse matrix elements and border elements (B) // SRCi is the current source of interest @@ -238,7 +336,20 @@ // diagonal elements of the sparse matrix (auto-cross-product) +# if (HAVE_MODEL_VAR) + f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); +# else f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); +# endif psSparseMatrixElement (sparse, i, i, f); +# if (HAVE_MODEL_VAR) + // if we have used CONSTANT errors, then we need to calculate the value of the parameter error + if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) { + float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(f); + } +# else // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { @@ -248,8 +359,12 @@ errors->data.F32[i] = 1.0 / sqrt(f); } - +# endif // find the image x model value +# if (HAVE_MODEL_VAR) + f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); +# else f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); +# endif psSparseVectorElement (sparse, i, f); @@ -257,11 +372,22 @@ switch (SKY_FIT_ORDER) { case 1: +# if (HAVE_MODEL_VAR) + f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal); + psSparseBorderElementB (border, i, 1, f); + f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal); + psSparseBorderElementB (border, i, 2, f); +# else f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 1, f); f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 2, f); +# endif case 0: +# if (HAVE_MODEL_VAR) + f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal); +# else f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); +# endif psSparseBorderElementB (border, i, 0, f); break; @@ -283,5 +409,9 @@ // got an overlap; calculate cross-product and add to output array +# if (HAVE_MODEL_VAR) + f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal); +# else f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); +# endif psSparseMatrixElement (sparse, j, i, f); } @@ -321,5 +451,9 @@ // set the sky, sky_x, sky_y components of border matrix +# if (HAVE_MODEL_VAR) + SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal); +# else SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal); +# endif psSparseConstraint constraint; @@ -479,2 +613,114 @@ return true; } + +# if (HAVE_MODEL_VAR) +bool psphotModelBackgroundReadout(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psMetadata *analysis, // Analysis metadata for outputs + pmReadout *readout, // Readout for which to generate a background model + psImageBinning *binning, // Binning parameters + const pmConfig *config,// Configuration + bool useVarianceImage + ); + +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) { + + bool status = false; + psRegion fullRegion = psRegionSet (0, 0, 0, 0); + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // create a model variance image + psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32); + + // find the binning information + psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config); + assert (backBinning); + + psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + + if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + // linear interpolation to full-scale + if (!psImageUnbin (modelVar, varModel, backBinning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + psFree (varModel); + psFree (varModelStdev); + + // XXX for a test: + psphotSaveImage (NULL, modelVar, "model.bck.fits"); + + // insert all of the source models + for (int i = 0; i < sources->n; i++) { + + // source of interest + pmSource *source = sources->data[i]; + + // skip sources which were not fitted already + if (!(source->mode & PM_SOURCE_MODE_LINEAR_FIT)) continue; + + // pixel region appropriate for the source + psRegion region = psRegionForImage (source->pixels, fullRegion); + + // define the source->modelVar pixels (view on modelVar image) + psAssert (!source->modelVar, "programming error : modelVar should be NULL here"); + psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here"); + psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here"); + source->modelVar = psImageSubset(modelVar, region); + + // add the source model to the model variance image + pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal); + } + + // XXX for a test: + psphotSaveImage (NULL, modelVar, "model.var.fits"); + psphotSaveImage (NULL, readout->variance, "image.var.fits"); + + // we save the model variance for future reference + psMetadataAddImage(readout->analysis, PS_LIST_TAIL, "PSPHOT.MODEL.VAR", PS_META_REPLACE, "model variance", modelVar); + psFree (modelVar); + + return true; +} + +bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) { + + bool status = false; + + // find the binning information + psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR"); + assert (modelVar); + + psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR"); + + // clear modelVar pointers for all of the source models + for (int i = 0; i < sources->n; i++) { + + // source of interest + pmSource *source = sources->data[i]; + psFree (source->modelVar); + } + + return true; +} +# endif Index: branches/meh_branches/ppstack_test/psphot/src/psphotImageLoop.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotImageLoop.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotImageLoop.c (revision 34041) @@ -75,8 +75,11 @@ pmFPAfileActivate (config->files, false, "PSPHOT.MASK"); pmFPAfileActivate (config->files, false, "PSPHOT.VARIANCE"); + // do not load the exposure number file yet file until we are done with the psphot analysis + pmFPAfileActivate (config->files, false, "PSPHOT.EXPNUM"); if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot."); // re-activate files so they will be closed and freed below - pmFPAfileActivate (config->files, true, NULL); +// XXX: Defer this +// pmFPAfileActivate (config->files, true, NULL); // there is now only a single chip (multiple readouts?). loop over it and process @@ -150,8 +153,26 @@ } } - // save output which is saved at the chip level - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + // Defer output and closing of files until we've (possibly) done the NFrames analysis below + pmFPAfileActivate (config->files, false, NULL); + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphot."); } - // save output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot."); + + // activate the EXPNUM image, cause it to be read in (if it exists), use it to set N Frames for the sources + // and finally iterate up to trigger writing of the output files + pmFPAfileActivate (config->files, true, "PSPHOT.EXPNUM"); + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed attempting to load EXPNUM input for Chip in psphot."); + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + if (!psphotSetNFrames (config, view, "PSPHOT.INPUT")) ESCAPE ("failed to setNFrames."); + } + } + // now activate all files to trigger the output + pmFPAfileActivate (config->files, true, NULL); + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed to iterate up Chip in psphot."); + } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); Index: branches/meh_branches/ppstack_test/psphot/src/psphotKronIterate.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotKronIterate.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotKronIterate.c (revision 34041) @@ -1,9 +1,7 @@ # include "psphotInternal.h" -# ifndef ROUND -# define ROUND(X) ((int) ((X) + 0.5*SIGN(X))) -# endif - -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert); -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal); + +bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow); +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels); + bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -11,6 +9,4 @@ bool status = true; - // return true; - fprintf (stdout, "\n"); psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---"); @@ -46,5 +42,5 @@ // psAssert (psf, "missing psf?"); - if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) { + if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -58,5 +54,5 @@ bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max); -bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) { bool status = false; @@ -69,4 +65,5 @@ psTimerStart ("psphot.kron"); + // determine the number of allowed threads int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads @@ -84,4 +81,63 @@ MIN_KRON_RADIUS = 0.25*RADIUS; } + + int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS"); + if (!status) { + KRON_ITERATIONS = 1; + } + + bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); + if (!status) { + KRON_APPLY_WEIGHT = true; + } + + bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); + if (!status) { + KRON_APPLY_WINDOW = false; + } + bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH"); + if (!status) { + KRON_SMOOTH = false; + } + float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA"); + if (!status) { + KRON_SMOOTH_SIGMA = 1.7; + } + float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA"); + if (!status) { + KRON_SMOOTH_NSIGMA = 2; + } + /* + * Parameter for calculating maximum integration radius based on source's surface + * brightness + * Given minimum surface brightness SBmin and a flux the maximum radius is found from + * + * SBmin = source->flux / (pi * Rmax**2) + * Rmax = sqrt (source->flux ) / sqrt (SBmin * pi) + * + * Now what do we use for SBmin? + * + * SBmin = ( some flux ) / (some area) + * some flux ~ flux of Ns sigma PSF source + * some area ~ K times area of a PSF + * + * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA + * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2 + * (the 4 accounts for the fact that 1 sigma is not the total area, it is + * actually a larger region). + * + * SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2) + * = Ns * SKY_STDEV * 4 / K + * + * We combine the two parameters Ns and K, and the constant 4 into a single recipe value + * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value + * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10 + */ + float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR"); + if (!status) { + KRON_SB_MIN_FACTOR = 10; + } + float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV"); + float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV ); // bit-masks to test for good/bad pixels @@ -114,5 +170,38 @@ // set a window function for each source based on the moments // (this skips really bad sources (no peak, no moments, DEFECT) - psphotKronWindowSetSource (source, kronWindow, false, true); + psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW); + } + + // We measure the Kron Radius on a smoothed copy of the readout image + psImage *smoothedImage = NULL; + if (KRON_SMOOTH) { + // Build the smoothed source image + // Replace the subtracted sources + psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false); + // Copy the image and smooth + psTimerStart ("psphot.kron.smooth"); + smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32); + psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA); + psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth")); + + // remove the sources + psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false ); + // Now subtract smooth versions of the sources from the smoothed image + psTimerStart ("psphot.kron.smooth.sources"); + for (int i=0; i< sources->n; i++) { + pmSource *source = sources->data[i]; + // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + // cache copy of smoothedPixels in the source + // tmpPtr is for use by a single "module" and must be null otherwise + psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!"); + + psImage *smoothedPixels = psImageSubset(smoothedImage, source->region); + source->tmpPtr = (psPtr) smoothedPixels; + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + } + } + psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources")); + } @@ -136,8 +225,14 @@ psArrayAdd(job->args, 1, kronWindow); psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, smoothedImage); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32); // set this to 0 to run without threading @@ -150,5 +245,5 @@ if (!psphotKronIterate_Threaded(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); + // psFree(AnalysisRegion); return false; } @@ -174,4 +269,12 @@ psFree (cellGroups); psFree (kronWindow); + if (KRON_SMOOTH) { + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + psFree(source->tmpPtr); + source->tmpPtr = NULL; + } + } + psFree (smoothedImage); psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); @@ -184,11 +287,16 @@ psImage *kronWindow = job->args->data[1]; psArray *sources = job->args->data[2]; - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); - float RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); - float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); - - // XXX TEST : set iteration to 1 - for (int j = 0; j < 1; j++) { + psImage *smoothedImage = job->args->data[3]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + float RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); + float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[7],F32); + int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[8],S32); + bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[9],S32); + bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[10],S32); + float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[11],F32); + float KRON_SB_MIN_DIVISOR = PS_SCALAR_VALUE(job->args->data[12],F32); + + for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -196,54 +304,88 @@ if (!source->peak) continue; // XXX how can we have a peak-less source? - // allocate space for moments + // check status of this source's moments if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; // replace object in image bool reSubtract = false; + psImage *smoothedPixels = NULL; if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + smoothedPixels = (psImage *) source->tmpPtr; + if (smoothedPixels) { + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0); + } reSubtract = true; } - // use S/N to control max window size - // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr; - - // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux) - // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2) - // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0) - // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS; - - // iterate to the window radius - // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow); - - // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS - // if we lack the skyRadius (eg MATCHED sources), go to the default value - float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + // On first iteration set window radius to sky radius (if valid) on second iteration + // use a factor times the previous radial moment value up to a maximum value that + // depends on the surface brightness of the source + float maxWindow; + if (j == 0) { + maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + } else { + if (KRON_SB_MIN_DIVISOR) { + if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) { + // Limit window radius based on surface brightness + float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR; + + if (source->moments->Mrf > 0) { + maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax); + } else { + maxWindow = Rmax; + } + } else { + maxWindow = RADIUS; + } + } else { + // old code + maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + } + } float windowRadius = PS_MAX(RADIUS, maxWindow); // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); + if (extend && smoothedPixels) { + psFree(source->tmpPtr); + smoothedPixels = psImageSubset(smoothedImage, source->region); + psAssert (smoothedPixels, "WTF?"); + source->tmpPtr = (psPtr) smoothedPixels ; + } + // clear the window function for this source based on the moments - psphotKronWindowSetSource (source, kronWindow, (j > 0), false); - - // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal); - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set a window function for each source based on the moments - psphotKronWindowSetSource (source, kronWindow, true, true); + // Note: this function also applies cuts on the source and returns false if it + // does not meet the requirements for measuring the Kron Radius or Magnitude. + // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false + if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) { + + // this function populates moments->Mrf,KronFlux,KronFluxErr + psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // set a window function for each source based on the moments + psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW); + } // if we subtracted it above, re-subtract the object, leave local sky if (reSubtract) { - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (smoothedPixels) { + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + } } } } + return true; } -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) { +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, + bool applyWeight, psImage *smoothedPixels) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -253,5 +395,5 @@ psF32 R2 = PS_SQR(radius); - float rsigma2 = 0.5 / PS_SQR(radius/2.0); + float rsigma2 = applyWeight ? 0.5 / R2 : 0; // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -266,4 +408,5 @@ // Xn = SUM (x - xc)^n * (z - sky) + psF32 RF = 0.0; psF32 RS = 0.0; @@ -284,5 +427,10 @@ int Ywo = source->pixels->row0; - psF32 **vPix = source->pixels->data.F32; + psF32 **vPix; + if (smoothedPixels) { + vPix = smoothedPixels->data.F32; + } else { + vPix = source->pixels->data.F32; + } psF32 **vWin = kronWindow->data.F32; psF32 **vWgt = source->variance->data.F32; @@ -357,4 +505,9 @@ float Var = 0.0; float Win = 0.0; + + // set vPix to the source pixels (it may have been set to the + // smoothed image above) + vPix = source->pixels->data.F32; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -395,5 +548,5 @@ } -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) { +bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) { if (!source) return false; @@ -402,5 +555,13 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) return false; if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false; psAssert(kronWindow, "need a window"); + + if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false; + + if (!applyWindow) { + return true; + } // we have a source with moments Mx, My, Mxx, Mxy, Myy. we just need to define a Gaussian that has @@ -413,19 +574,29 @@ float Yo = source->moments->My; - float Mxx = source->moments->Mxx; - float Mxy = source->moments->Mxy; - float Myy = source->moments->Myy; - - float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - - // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width - float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor; - - float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2 - float Sxy = Mxy / (scale * Mmajor * Mminor); - float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2 - - float Smajor = sqrt(Mmajor); + psEllipseMoments moments; + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0); + if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) { + // Shall we log a proper warning? This happens often with matched sources (forced photometry) + // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta); + return false; + } + + // Why this factor of 0.5 ? + float scale = 0.5 * source->moments->Mrf / axes.major; + axes.major *= scale; + axes.minor *= scale; + + psEllipseShape shape = psEllipseAxesToShape(axes); + if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) { + // Shall we log a proper warning? This happens often with matched sources (forced photometry) + // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy); + return false; + } + + float Smajor = axes.major; int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1); @@ -434,7 +605,9 @@ int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1); - float rMxx = 0.5 / Sxx; - float rMyy = 0.5 / Syy; - + float rMxx = 0.5 / PS_SQR(shape.sx); + float rMyy = 0.5 / PS_SQR(shape.sy); + float Sxy = -1. * shape.sxy; // factor of -1 is included to match the previous window function + // implementation. XXX: Is this correct? + for (int iy = minY; iy < maxY; iy++) { for (int ix = minX; ix < maxX; ix++) { @@ -447,7 +620,8 @@ float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z)); - kronWindow->data.F32[iy][ix] *= f; + kronWindow->data.F32[iy][ix] *= f; } } + return true; } Index: branches/meh_branches/ppstack_test/psphot/src/psphotMergeSources.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotMergeSources.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotMergeSources.c (revision 34041) @@ -677,5 +677,8 @@ if (sourceOut->modelPSF) { pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); - } + } else { + // if we have no pixels we can't use it to determine the psf so make sure this bit is off + sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; + } // child sources have not been subtracted in this image, but this flag may be raised if @@ -741,6 +744,13 @@ psAssert (readoutOut, "missing readout?"); + pmDetections *detectionsOutOld = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS"); + psArray *oldFootprints = detectionsOutOld ? detectionsOutOld->footprints : NULL; + // replace any existing DETECTION container on readoutOut->analysis with the new one pmDetections *detectionsOut = pmDetectionsAlloc(); + if (oldFootprints) { + // ... but hang on to any existing footprints so that they can be merged with new footprints in pass 2 + detectionsOut->footprints = psMemIncrRefCounter(oldFootprints); + } detectionsOut->allSources = psArrayAllocEmpty (100); if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) { @@ -900,6 +910,10 @@ // allocate image, weight, mask for the new image for each peak if (sourceOut->modelPSF) { - pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); - } + pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, + sourceOut->modelPSF->fitRadius); + } else { + // if we have no pixels we can't use it to determine the psf so make sure this bit is off + sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; + } // child sources have not been subtracted in this image, but this flag may be raised if Index: branches/meh_branches/ppstack_test/psphot/src/psphotModelBackground.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotModelBackground.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotModelBackground.c (revision 34041) @@ -37,10 +37,11 @@ // corresponding to the model. Other information about the background model is saved on the // readout->analysis -static bool psphotModelBackgroundReadout(psImage *model, // Model image - psImage *modelStdev, // Model stdev image - psMetadata *analysis, // Analysis metadata for outputs - pmReadout *readout, // Readout for which to generate a background model - psImageBinning *binning, // Binning parameters - const pmConfig *config // Configuration +bool psphotModelBackgroundReadout(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psMetadata *analysis, // Analysis metadata for outputs + pmReadout *readout, // Readout for which to generate a background model + psImageBinning *binning, // Binning parameters + const pmConfig *config,// Configuration + bool useVarianceImage ) { @@ -49,5 +50,6 @@ bool status = true; - psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout + psImage *image = useVarianceImage ? readout->variance : readout->image; + psImage *mask = readout->mask; // Image and mask for readout // select the appropriate recipe information @@ -143,5 +145,7 @@ // we save the binning structure for use in psphotMagnitudes - psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + if (!useVarianceImage) { + psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + } psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -159,6 +163,6 @@ // measure clipped median for subimages - psRegion ruffRegion = {0,0,0,0}; - psRegion fineRegion = {0,0,0,0}; + psRegion ruffRegion = psRegionSet (0,0,0,0); + psRegion fineRegion = psRegionSet (0,0,0,0); for (int iy = 0; iy < model->numRows; iy++) { for (int ix = 0; ix < model->numCols; ix++) { @@ -330,25 +334,30 @@ psVectorStats (statsDQ, dQ, NULL, NULL, 0); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); - - // measure image and background stats and save for later output - psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | - PS_STAT_SAMPLE_STDEV | - PS_STAT_MIN | - PS_STAT_MAX); - psImageStats (statsBck, model, NULL, 0); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); - psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); - psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); - psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", - statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); + if (!useVarianceImage) { + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + + // measure image and background stats and save for later output + psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | + PS_STAT_SAMPLE_STDEV | + PS_STAT_MIN | + PS_STAT_MAX); + psImageStats (statsBck, model, NULL, 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); + psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", + statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); + psFree(statsBck); + } else { + psLogMsg ("psphot", PS_LOG_INFO, "variance data : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + } psFree(statsDQ); @@ -356,5 +365,4 @@ psFree(stats); - psFree(statsBck); psFree(statsDefaults); psFree(binning); @@ -375,5 +383,5 @@ psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation - if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config, false)) { psFree(model); psFree(modelStdev); @@ -400,5 +408,5 @@ pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); - if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); return false; Index: branches/meh_branches/ppstack_test/psphot/src/psphotParseCamera.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotParseCamera.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotParseCamera.c (revision 34041) @@ -2,4 +2,6 @@ // define the needed / desired I/O files + + bool psphotParseCamera (pmConfig *config) { @@ -40,4 +42,11 @@ } + pmFPAfileBindFromArgs (&status, input, config, "PSPHOT.EXPNUM", "EXPNUM"); + if (!status) { + psError (PS_ERR_UNKNOWN, false, "failed to load find definition"); + return NULL; + } + + // define the additional input/output files associated with psphot if (!psphotDefineFiles (config, input)) { Index: branches/meh_branches/ppstack_test/psphot/src/psphotRadialProfileWings.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotRadialProfileWings.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotRadialProfileWings.c (revision 34041) @@ -55,4 +55,5 @@ static float MIN_RADIUS = NAN; static float SKY_STDEV = NAN; +static float SKY_SLOPE_MIN = NAN; // static FILE *file = NULL; @@ -81,5 +82,5 @@ MIN_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); - if (!status) { + if (!status || (MIN_RADIUS > MAX_RADIUS)) { MIN_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); } @@ -89,4 +90,10 @@ if (!status) { SKY_STDEV = 1.0; // a crude default value (why would this not exist?) + } + + // SKY_SLOPE_MIN is the sigma of the sky model (ie, smoothed on large scales) + SKY_SLOPE_MIN = psMetadataLookupF32 (&status, recipe, "SKY_SLOPE_MIN"); + if (!status) { + SKY_SLOPE_MIN = 3.0; } @@ -302,10 +309,16 @@ limitSlope = slope; } - if (!limit && isfinite(slope) && (fabs(slope) < 3.0)) { + if (!limit && isfinite(slope) && (fabs(slope) < SKY_SLOPE_MIN)) { // SB no longer changing. limit = true; // linearly interpolate to the radius at which we hit the sky, using the last flux and the limiting slope if (isfinite(lastFlux)) { - limitRadius = lastRadius + lastFlux / 3.0; + float interpolatedRadius = lastRadius + lastFlux / SKY_SLOPE_MIN; + if (interpolatedRadius < MAX_RADIUS) { + limitRadius = interpolatedRadius; + } else { + // XXX should we keep going in this case? + limitRadius = meanRadius; + } } else { limitRadius = meanRadius; Index: branches/meh_branches/ppstack_test/psphot/src/psphotReadout.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotReadout.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotReadout.c (revision 34041) @@ -185,5 +185,9 @@ // remove noise for subtracted objects (ie, return to normal noise level) // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); + // Note: if footprintUseUnsubtracted is true the noise was already subtracted in psphotFindDetections() + if (!footprintUseUnsubtracted) { + psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + } // define new sources based on only the new peaks & measure moments Index: branches/meh_branches/ppstack_test/psphot/src/psphotReadoutKnownSources.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotReadoutKnownSources.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotReadoutKnownSources.c (revision 34041) @@ -71,4 +71,7 @@ psphotSourceFreePixels (config, view, filerule); + // set source->nFrames + psphotSetNFrames (config, view, filerule); + // create the exported-metadata and free local data return psphotReadoutCleanup (config, view, filerule); Index: branches/meh_branches/ppstack_test/psphot/src/psphotReadoutMinimal.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotReadoutMinimal.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotReadoutMinimal.c (revision 34041) @@ -64,5 +64,8 @@ // Construct an initial model for each object, set the radius to fitRadius, set circular // fit mask. NOTE: only applied to sources without guess models - psphotGuessModels (config, view, filerule); + if (!psphotGuessModels (config, view, filerule)) { + psLogMsg ("psphot", 3, "failure to Guess Model"); + return psphotReadoutCleanup (config, view, filerule); + } // linear PSF fit to source peaks Index: branches/meh_branches/ppstack_test/psphot/src/psphotReplaceUnfit.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotReplaceUnfit.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotReplaceUnfit.c (revision 34041) @@ -60,5 +60,6 @@ psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); + //psAssert (sources, "missing sources?"); + if (!sources) return true; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -139,6 +140,7 @@ psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - + //psAssert (sources, "missing sources?"); + if (!sources) return true; + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels Index: branches/meh_branches/ppstack_test/psphot/src/psphotSersicModelClass.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotSersicModelClass.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotSersicModelClass.c (revision 34041) @@ -546,4 +546,8 @@ // set the normalization by linear fit between model and data psphotSersicModelNorm (pcm, source); + if (!isfinite(PAR[PM_PAR_I0])) { + fprintf(stderr, "psphotSersicModelClassGuessPCM: psphotSerisicModelNorm set PM_PAR_I0 to NAN\n"); + return false; + } if (TIMING) { t5 = psTimerMark ("SersicGuess"); } Index: branches/meh_branches/ppstack_test/psphot/src/psphotSetNFrames.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotSetNFrames.c (revision 34041) +++ branches/meh_branches/ppstack_test/psphot/src/psphotSetNFrames.c (revision 34041) @@ -0,0 +1,69 @@ +# include "psphotInternal.h" + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + + +// Set source->nFrames based on the values in the EXPNUM image at the coordinates of each source's peak + +bool psphotSetNFrames (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Set N Frames ---"); + + int num = psphotFileruleCount(config, filerule); + + // loop over the inputs + for (int i = 0; i < num; i++) { + if (!psphotSetNFramesReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to set nFrames for %s entry %d", filerule, i); + return false; + } + } + + return true; +} + +bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) +{ + bool status; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psImage *expnum = psMetadataLookupPtr(&status, readout->analysis, "EXPNUM"); + if (!status || !expnum) { + psLogMsg ("psphot", PS_LOG_INFO, "No EXPNUM image for input %d", index); + return true; + } + psLogMsg ("psphot", PS_LOG_INFO, "Found EXPNUM image for input %d", index); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + psAssert (detections->allSources, "missing sources?"); + + psArray *sources = detections->allSources; + + int col0 = expnum->col0; + int row0 = expnum->row0; + int numCols = expnum->numCols; + int numRows = expnum->numRows; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + int x = source->peak->x; + int y = source->peak->y; + if (x >= col0 && x < numCols && y >= row0 && y < numRows) { + source->nFrames = expnum->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; + } else { + source->nFrames = 0; + } + } + + return true; +} Index: branches/meh_branches/ppstack_test/psphot/src/psphotSetThreads.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotSetThreads.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotSetThreads.c (revision 34041) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 7); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task); Index: branches/meh_branches/ppstack_test/psphot/src/psphotSourceMatch.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotSourceMatch.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotSourceMatch.c (revision 34041) @@ -3,4 +3,8 @@ bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); bool psphotMatchSourcesSetIDs (psArray *objects); + +static psArray *psphotMatchFootprintCacheAlloc (int nImages); +static pmFootprint *psphotMatchLookupFootprint (psArray *cache, int footprintID, int imageID); +static pmFootprint *psphotMatchCopyFootprint (psArray *cache, pmFootprint *footprint, int imageID, psImage *image, psArray *footprints); psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -193,4 +197,5 @@ // vector to track if source for an image is found psVector *found = psVectorAlloc(nImages, PS_TYPE_U8); + psArray *footprintCache = psphotMatchFootprintCacheAlloc(nImages); for (int i = 0; i < objects->n; i++) { @@ -219,5 +224,5 @@ // we make a copy of the largest footprint; this will be used for all new sources associated with this object - pmFootprint *footprint = NULL; + pmFootprint *largestFootprint = NULL; if (iSpansMax != -1) { // copy the footprint info pmSource *src = obj->sources->data[iSpansMax]; @@ -226,6 +231,5 @@ psAssert(src->peak->footprint->nspans == nSpansMax, "wrong footprint?"); - // we only care about the spans, do not worry about the image of this footprint - footprint = pmFootprintCopyData(src->peak->footprint, NULL); + largestFootprint = src->peak->footprint; } @@ -248,11 +252,16 @@ // assign to a footprint on this readout->image - if (footprint) { - peak->footprint = pmFootprintCopyData(footprint, readout->image); - - // the peak does not claim ownership of the footprint (it does not free it). save a copy of this - // footprint on detections->footprints so we can free it later - psArrayAdd(detections->footprints, 100, peak->footprint); - psFree (peak->footprint); + if (largestFootprint) { + // we save the copies that we make of the of the footprints in a hash so that we can reuse them + // for all sources that share the fooprint. Without this we had serious memory explosion in + // dense fields with lots of footprints (spans are small but when you have enough of them ...) + peak->footprint = psphotMatchLookupFootprint(footprintCache, largestFootprint->id, index); + if (!peak->footprint) { + // the peak does not claim ownership of the footprint (it does not free it). + // psphotMatchCopyFootprint saves a copy of this + // footprint on detections->footprints so we can free it later + peak->footprint = psphotMatchCopyFootprint(footprintCache, largestFootprint, + index, readout->image, detections->footprints); + } } @@ -272,4 +281,9 @@ pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); +#if (0) + fprintf(stderr, "Add mising source for obj: %5d %5d image: %d flux: %f size: %4d %4d\n", + i, obj->id, index, peakFlux, source->pixels->numRows, source->pixels->numCols); +#endif + peak->assigned = true; pmPhotObjAddSource(obj, source); @@ -277,6 +291,6 @@ psFree (source); } - psFree (footprint); - } + } + psFree(footprintCache); // how many sources do we have now? @@ -309,2 +323,54 @@ return true; } + +// Cache of footprints created for unmatched sources. +static psArray * psphotMatchFootprintCacheAlloc (int nImages) { + psArray *cache = psArrayAlloc(nImages); + for (int i = 0; i < nImages; i++) { + psHash *hash = psHashAlloc(5000); + cache->data[i] = hash; + } + return cache; +} + +// Find copy of footprint with given ID made for a given image +static pmFootprint *psphotMatchLookupFootprint (psArray *cache, int footprintID, int imageID) { + psHash *hash = (psHash *) cache->data[imageID]; + + psAssert(hash != NULL, "missing hash for image %d", imageID); + + // footprintID is the id of the original footprint. + char key[32]; + sprintf(key, "%d", footprintID); + + // The footprints in our hashes are the copies of that footprint for the respective images + pmFootprint *copy = psHashLookup(hash, key); + + return copy; +} + +// Create a copy of a given footprint for a given image +static pmFootprint *psphotMatchCopyFootprint (psArray *cache, pmFootprint *footprint, int imageID, psImage *image, psArray *footprints) { + + psAssert((imageID >= 0 && imageID < cache->n), "invalid imageID %d", imageID); + + psHash *hash = (psHash *) cache->data[imageID]; + + psAssert(hash != NULL, "missing hash for image %d", imageID); + + char key[32]; + sprintf(key, "%d", footprint->id); + + pmFootprint *copy = pmFootprintCopyData (footprint, image); + + // save a copy in this image's hash + psHashAdd(hash, key, copy); + + // save a copy of this footprint on the passed in footprints Array so we can free it later + psArrayAdd(footprints, 100, copy); + + // drop our ref + psFree(copy); + + return copy; +} Index: branches/meh_branches/ppstack_test/psphot/src/psphotSourceStats.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotSourceStats.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotSourceStats.c (revision 34041) @@ -605,4 +605,7 @@ } psAssert (isfinite(Sigma), "did we miss a case?"); + if (!isfinite(minKronRadius)) { + minKronRadius = Sigma; + } // choose a grid scale that is a fixed fraction of the psf sigma^2 Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackArguments.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackArguments.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackArguments.c (revision 34041) @@ -53,4 +53,8 @@ psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-ds9regions"))) { + psArgumentRemove (N, &argc, argv); + pmSubtractionRegions(true); } Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackImageLoop.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackImageLoop.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackImageLoop.c (revision 34041) @@ -30,7 +30,17 @@ } + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW"); + if (useRaw && inputRaw == NULL) { + psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input."); + useRaw = false; + } + pmFPAview *view = pmFPAviewAlloc (0); - // XXX for now, just load the full set of images up front + + // XXX for now, just load the full set of images up front except for EXPNUM which we defer + pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW"); + pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV"); if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); @@ -86,6 +96,9 @@ UpdateHeadersForChip(config, view); - // save output which is saved at the chip level - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + // Defer output until we have performed psphotSetNFrames() + pmFPAfileActivate (config->files, false, NULL); + + // Iterate up + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphot."); } psMemDump("doneloop"); @@ -93,5 +106,36 @@ UpdateHeadersForFPA(config, view); - // save output which is saved at the fpa level + // Iterate up output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput pmFPAfileIOChecks FPA in psphot."); + + // Load the appropriate EXPNUM image + pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV"); + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot."); + + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip EXPNUM in psphotStack."); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack."); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 6, "Readout %d: %x %x\n", view->readout, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames."); + } + } + // now activate all files to trigger final output + pmFPAfileActivate (config->files, true, NULL); + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsPrepare.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsPrepare.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsPrepare.c (revision 34041) @@ -32,7 +32,7 @@ pmDetections *detections = psMetadataLookupPtr(NULL, readout->analysis, "PSPHOT.DETECTIONS"); // Sources if (!detections || !detections->allSources) { - psWarning("No detections found for image %d --- rejecting.", index); - options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0x01; - return true; // XXX not an error: we continue processing other images + // An input must have sources. If we got here without any the input file was likely corrupt + psError(PSPHOT_ERR_IO, true, "Input %d has no sources.\n\tCheck all instances of: %s", index, inputSrc->origname); + return false; } psAssert (detections->allSources, "missing sources?"); Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsUtils.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackMatchPSFsUtils.c (revision 34041) @@ -105,9 +105,12 @@ psAssert(psphotRecipe, "Need PSPHOT recipe"); +#if (0) + // DON'T OVERWRITE PSPHOT recipe value psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); + psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); +#endif psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackObjects.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackObjects.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackObjects.c (revision 34041) @@ -119,6 +119,4 @@ } if (!skipSourcePetro) keepObjectPetro = true; - - keepObjectPetro = true; } Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackParseCamera.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackParseCamera.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackParseCamera.c (revision 34041) @@ -59,4 +59,11 @@ } } + psString expnum = psMetadataLookupStr(&status, input, "RAW:EXPNUM"); // Name of expnum image + if (expnum && strlen(expnum) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.EXPNUM.RAW", expnum, PM_FPA_FILE_EXPNUM)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from expnum %d (%s)", i, expnum); + return false; + } + } nRaw ++; } @@ -84,4 +91,11 @@ } } + psString expnum = psMetadataLookupStr(&status, input, "CNV:EXPNUM"); // Name of EXPNUM image + if (expnum && strlen(expnum) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.EXPNUM.CNV", expnum, PM_FPA_FILE_EXPNUM)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from expnum %d (%s)", i, expnum); + return false; + } + } psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask if (psf && strlen(psf) > 0) { @@ -115,4 +129,9 @@ } + psS64 stack_id = psMetadataLookupS64(&status, input, "STACK_ID"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to find STACK_ID from sources %d", i); + return false; + } // generate an pmFPAimage for the output convolved image // XXX output of these files should be optional @@ -125,5 +144,5 @@ } outputImage->save = true; - outputImage->fileID = i; // this is used to generate output names + outputImage->fileID = stack_id; // this is used to generate output names pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK"); @@ -137,5 +156,5 @@ } outputMask->save = true; - outputMask->fileID = i; // this is used to generate output names + outputMask->fileID = stack_id; // this is used to generate output names pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE"); @@ -149,5 +168,5 @@ } outputVariance->save = true; - outputVariance->fileID = i; // this is used to generate output names + outputVariance->fileID = stack_id; // this is used to generate output names // the output sources are carried on the outputImage->fpa structures @@ -158,5 +177,5 @@ } outsources->save = true; - outsources->fileID = i; // this is used to generate output names + outsources->fileID = stack_id; // this is used to generate output names } } Index: branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotStackReadout.c (revision 34041) @@ -89,6 +89,16 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } - if (!psphotSubtractBackground (config, view, STACK_SRC)) { - return psphotReadoutCleanup (config, view, STACK_SRC); + if (strcmp(STACK_SRC, STACK_DET)) { +#define MODEL_BACKGROUND_SRC 1 +#ifdef MODEL_BACKGROUND_SRC + // work around the fact that the background levels on the convolved + // and unconvolved stacks can be different + if (!psphotModelBackground (config, view, STACK_SRC)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } +#endif + if (!psphotSubtractBackground (config, view, STACK_SRC)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } } if (!strcasecmp (breakPt, "BACKMDL")) { @@ -190,4 +200,5 @@ psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources) + // if we only do one pass, skip to extended source analysis if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; @@ -211,5 +222,5 @@ // subtract all sources from DET (this will subtract using the psf model for SRC, which // will somewhat oversubtract the sources -- this is OK - psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources + psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources } @@ -223,5 +234,9 @@ // remove noise for subtracted objects (ie, return to normal noise level) // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources) + // NOTE: if fooprintsUseUnsubtracted, the noise has already been removed by psphotFindDetections + bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); + if (!footprintsUseUnsubtracted) { + psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources) + } // if DET and SRC are different images, copy the detections from DET to SRC @@ -302,7 +317,7 @@ psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources) - // create source children for the OUT filerule (for radial aperture photometry) - psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects); - if (!objectsRadial) { + // create source children for the OUT filerule (for radial aperture photometry and output) + psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects); + if (!objectsOut) { psFree(objects); psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); @@ -310,33 +325,51 @@ } - // measure circular, radial apertures (objects sorted by S/N) - // this forces photometry on the undetected sources from other images - psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched? - psMemDump("extmeas"); - - int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); - - for (int entry = 1; entry < nRadialEntries; entry++) { - // NOTE: entry 0 is the unmatched image set - - // re-measure the PSF for the smoothed image (using entries in 'allSources') - psphotChoosePSF (config, view, STACK_OUT, false); - - // this is necessary to update the models based on the new PSF - psphotResetModels (config, view, STACK_OUT); - - // this is necessary to get the right normalization for the new models - psphotFitSourcesLinear (config, view, STACK_OUT, false); - - // measure circular, radial apertures (objects sorted by S/N) - // entry 0 == unmatched? pass entry + 1? - psphotRadialApertures (config, view, STACK_OUT, entry); - - // replace the flux in the image so it is returned to its original state - psphotReplaceAllSources (config, view, STACK_OUT, false); - - // smooth to the next FWHM, or set 'smoothAgain' to false if no more - psphotStackMatchPSFsNext(config, view, STACK_OUT, entry); - psMemDump("matched"); + + bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES"); + if (radial_apertures) { + // measure circular, radial apertures (objects sorted by S/N) + // this forces photometry on the undetected sources from other images + + // NOTE: we always do the radial apertures analysis on the convolved image since + // those are the ones that are psf matched and the source of STACK_OUT's pixels + // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true. + // but in that case we don't get past the psf matching step because there is no + // target psf for the RAW inputs + + // If useRaw copy the sources to the convolved readout + if (strcmp(STACK_SRC, STACK_CNV)) { + if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + } + psphotRadialApertures (config, view, STACK_CNV, 0); // XXX entry 0 == unmatched? + psMemDump("extmeas"); + + int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); + + for (int entry = 1; entry < nRadialEntries; entry++) { + // NOTE: entry 0 is the unmatched image set + + // re-measure the PSF for the smoothed image (using entries in 'allSources') + psphotChoosePSF (config, view, STACK_OUT, false); + + // this is necessary to update the models based on the new PSF + psphotResetModels (config, view, STACK_OUT); + + // this is necessary to get the right normalization for the new models + psphotFitSourcesLinear (config, view, STACK_OUT, false); + + // measure circular, radial apertures (objects sorted by S/N) + // entry 0 == unmatched? pass entry + 1? + psphotRadialApertures (config, view, STACK_OUT, entry); + + // replace the flux in the image so it is returned to its original state + psphotReplaceAllSources (config, view, STACK_OUT, false); + + // smooth to the next FWHM, or set 'smoothAgain' to false if no more + psphotStackMatchPSFsNext(config, view, STACK_OUT, entry); + psMemDump("matched"); + } } @@ -344,4 +377,5 @@ if (!psphotApResid (config, view, STACK_SRC)) { psFree (objects); + psFree (objectsOut); psLogMsg ("psphot", 3, "failed on psphotApResid"); return psphotReadoutCleanup (config, view, STACK_SRC); @@ -375,5 +409,5 @@ psFree (objects); - psFree (objectsRadial); + psFree (objectsOut); // create the exported-metadata and free local data Index: branches/meh_branches/ppstack_test/psphot/src/psphotTest.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotTest.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotTest.c (revision 34041) @@ -103,5 +103,5 @@ # if (0) - psRegion region = {0,0,0,0}; // a region representing the entire array +psRegion region = psRegionSet (0,0,0,0); // a region representing the entire array psphotTestArguments (&argc, argv); Index: branches/meh_branches/ppstack_test/psphot/src/psphotTestSourceOutput.c =================================================================== --- branches/meh_branches/ppstack_test/psphot/src/psphotTestSourceOutput.c (revision 33415) +++ branches/meh_branches/ppstack_test/psphot/src/psphotTestSourceOutput.c (revision 34041) @@ -136,5 +136,5 @@ // generate working image for this source - psRegion region = {ix - dx, ix + dx, iy - dy, iy + dy}; + psRegion region = psRegionSet(ix - dx, ix + dx, iy - dy, iy + dy); psImage *vM = psImageSubset (imMo, region); Index: branches/meh_branches/ppstack_test/psphot/test/tap_psphot_varmodel.pro =================================================================== --- branches/meh_branches/ppstack_test/psphot/test/tap_psphot_varmodel.pro (revision 34041) +++ branches/meh_branches/ppstack_test/psphot/test/tap_psphot_varmodel.pro (revision 34041) @@ -0,0 +1,642 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# $KAPA = kapa -noX + +# PSF.CONVOLVE : if true, we insert delta functions (and optionally +# galaxies) and smooth the image with the psf model +# (uses a GAUSS regardless of the model). Note that +# PSF.CONVOLVE = T is faster than F, but (a) only +# allows Gauss models and (b) only yields quantized +# locations + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM F + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# create an image with fake sources and insert the resulting cmf file into a dvodb +$RefConfig = -camera SIMTEST +$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE +$RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RefConfig = $RefConfig -Db PSF.CONVOLVE F + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 +$RefOptions = $RefOptions -seeing 1.0 +$RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS +$RefOptions = $RefOptions -Df STARS.DENSITY 10.0 +$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5 + +# basic config for ppSim with randomly distributed stars and NO galaxies +$RealConfig = -camera SIMTEST +$RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN +$RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RealConfig = $RealConfig -Db STARS.FAKE F +$RealConfig = $RealConfig -Db STARS.REAL T +$RealConfig = $RealConfig -Db MATCH.DENSITY F +$RealConfig = $RealConfig -Db PSF.CONVOLVE F +$RealConfig = $RealConfig -Df STARS.DENSITY 10.0 +$RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5 +$RealConfig = $RealConfig -Db GALAXY.FAKE F +$RealConfig = $RealConfig -Db GALAXY.GRID F + +# options for the repeated images +$RealOptions = $BaseOptions -exptime 30.0 + +$ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS + +# sample alternate options: +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0 +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro go + + mkdir test + + $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS + mkexp test/image.00 1.0 + $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS + mkexp test/image.01 1.0 + $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1 + mkexp test/image.02 1.0 +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + exec rm -rf catdir.ref + exec rm -f refimage.fits + + exec time ppSim $RefConfig $RefOptions refimage -nx 3000 -ny 3000 + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR catdir.ref -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes refimage.cmf + exec relphot -averages -D CATDIR catdir.ref -update -region 260 280 -33 -13 +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 3) + echo "USAGE: mkexp basename fwhm" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + + # create the raw image + echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename + exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename +end + +macro runphot + if ($0 != 4) + echo "USAGE: runphot basename outname options" + break + end + + local basename + $basename = $1 + $outname = $2 + $options = $3 + + # create the chip output + echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options + exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options +end + +# compare two cmf files with extname Chip.psf +# things to compare: +# * completeness (which sources in (1) are not detected in (2) +# * positions (X_PSF, Y_PSF) +# * instrumental psf mags +# * position errors (no input errors; use a model?) +# * measured FWHM? +# * kron mags (fluxes) +# * etc, etc +macro ckchip + if ($0 != 5) + echo "USAGE: ckchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i NX NY nx ny N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end + + # plot (input - output) vs mag +end + +macro stchip + if ($0 != 5) + echo "USAGE: stchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +# compare chip to warp +macro ckwarp + if ($0 != 5) + echo "USAGE: ckwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.0 -index1 index1 -index2 index2 + + local i nx ny NX NY N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end +end + +macro stwarp + if ($0 != 5) + echo "USAGE: stwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +macro stats.pair + if ($0 != 3) + echo "USAGE: stats.pair (N) (output)" + break + end + + list word -split $spairs:$1 + if ($word:n != 8) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + reindex v2 = $word:1 using index2 + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7) + subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7) + + vstats -q d1 -sigma-clip 3.0 + $M1 = $MEAN + $S1 = $SIGMA + vstats -q d2 -sigma-clip 3.0 + $M2 = $MEAN + $S2 = $SIGMA + + output $2 + fprintf "%-18s %7.4f %7.4f %7.4f %7.4f" $word:0 $M1 $S1 $M2 $S2 + output stdout +end + +macro show.pair + if ($0 != 2) + echo "USAGE: show.pair (N)" + break + end + + list word -split $pairs:$1 + if ($word:n != 7) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + if ("$word:6" == "V") + reindex v2 = $word:1 using index2 + end + if ("$word:6" == "S") + set v2 = $word:1 + zero(index1) + end + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + if (("$word:4" == "def") || ("$word:5" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:4 $word:5; box; plot rv delta + end + label -y '$word:0' -x '$word:2' +end + +# This list is used to compare a pair of vectors (sans error) or a +# vector and an expected (constant) value. The last field defines a +# vector or constant for the comparison. It is assumed that the +# vector sets have been loaded and matched with match2d to generate +# index vectors 'index1' and index2'. The macro 'show.pair' generates +# a plot of the range vector vs (v1 - v2). The indices for v1, v2 are +# index1 and 2 respectively. The index for the range vector is defined +# by the integer following that vector. the y-limits of the plot are +# given by the last two numbers +list pairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_out X_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_SIG_out X_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_SIG_out Y_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_MAG_out PSF_INST_MAG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_SIG_out PSF_INST_MAG_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_FLUX_out PSF_INST_FLUX_raw PSF_INST_MAG_raw 2 def def V + #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2 def def V + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RAW_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RADIUS_out 0.0 PSF_INST_MAG_raw 2 -0.01 20.1 S + SKY_out 0.0 PSF_INST_MAG_raw 2 def def S + SKY_SIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_CHISQ_out 1.0 PSF_INST_MAG_raw 2 def def S + CR_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + EXT_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + PSF_MAJOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_MINOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_THETA_out 0.0 PSF_INST_MAG_raw 2 -1.61 1.61 S + PSF_QF_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_QF_PERFECT_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_NDOF_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_NPIX_out 0.0 PSF_INST_MAG_raw 2 def def S + MOMENTS_XX_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_XY_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_YY_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_M3C_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M3S_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M4C_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_M4S_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_R1_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + MOMENTS_RH_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + K_out M_raw PSF_INST_MAG_raw 2 def def V + KRON_FLUX_ERR_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_INNER_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_OUTER_out 0.0 PSF_INST_MAG_raw 2 def def S +# CAL_PSF_MAG CAL_PSF_MAG none Mraw 2 -1.01 1.01 V +# CAL_PSF_MAG_SIG CAL_PSF_MAG_SIG none Mraw 2 -1.01 1.01 V +# RA_PSF RA_PSF none Mraw 2 -1.01 1.01 V +# DEC_PSF DEC_PSF none Mraw 2 -1.01 1.01 V +# PEAK_FLUX_AS_MAG PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V +# FLAGS FLAGS 0.0 Mraw 2 -1.01 1.01 S +# FLAGS2 FLAGS2 0.0 Mraw 2 -1.01 1.01 S +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG + PSF_INST_FLUX + PSF_INST_FLUX_SIG + AP_MAG + AP_MAG_RAW + AP_MAG_RADIUS + SKY + SKY_SIGMA + PSF_CHISQ + CR_NSIGMA + EXT_NSIGMA + PSF_MAJOR + PSF_MINOR + PSF_THETA + PSF_QF + PSF_QF_PERFECT + PSF_NDOF + PSF_NPIX + MOMENTS_XX + MOMENTS_XY + MOMENTS_YY + MOMENTS_M3C + MOMENTS_M3S + MOMENTS_M4C + MOMENTS_M4S + MOMENTS_R1 + MOMENTS_RH + KRON_FLUX + KRON_FLUX_ERR + KRON_FLUX_INNER + KRON_FLUX_OUTER +# CAL_PSF_MAG +# CAL_PSF_MAG_SIG +# RA_PSF +# DEC_PSF +# PEAK_FLUX_AS_MAG +# FLAGS +# FLAGS2 +end + +# use these cmf entries to measure average stats of the given pairs +list spairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 +end + +# # XXX this is a hack: the cmf file created by ppSim is not +# # compatible with the image for warp (because one is a chip-mosaic +# # and the other is not) +# exec rm -f fix.hdr +# output fix.hdr +# echo "PSMOSAIC= 'CHIP ' / Mosaicked level" +# output stdout +# fix the header to be compatible with the chip file (to avoid running psphot) +# exec fits_insert $basename.cmf fix.hdr + +macro completeness + if ($0 != 4) + echo "USAGE: completeness (raw) (out) (output)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out 1.5 -index1 index1 -index2 index2 -closest + + histogram PSF_INST_MAG_raw nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (PSF_INST_MAG_raw > dMag[$i]) && (PSF_INST_MAG_raw <= dMag[$i+1]) + subset all = index1 if (inrange) + subset got = index1 if (inrange) && (index1 >= 0) + if (all[] == 0) + fMag[$i] = 0 + else + fMag[$i] = got[] / all[] + end + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + label -x mag_inst -y det_frac + + set found = (index1 >= 0) + plot PSF_INST_MAG_raw found -c red + + set ldmag = log(PSF_INST_MAG_SIG_raw) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot PSF_INST_MAG_raw ldmag + label +y log(S/N) + + png -name $3 +end + +macro show.dpair + if ($0 != 2) + echo "USAGE: show.dpair (N)" + break + end + + list word -split $dpairs:$1 + if ($word:n != 7) + echo "invalid dpair $1" + break + end + + $Nr = $word:4 + + reindex v1 = $word:0 using index1 + reindex dv = $word:1 using index1 + reindex v2 = $word:2 using index2 + reindex rv = $word:3 using index$Nr + + set delta = (v1 - v2) / dv + if (("$word:5" == "def") || ("$word:6" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:5 $word:6; box; plot rv delta + end + label -y '$word:0' -x '$word:3' +end + +# this list is used to compare a pair of vectors with an error it is +# assumed that the vector sets have been loaded and matched with +# match2d to generate index vectors 'index1' and index2'. the macro +# show.dpair generates a plot of the range vector vs (v1 - v2) / dv1. +# The indices for v1, dv1, and v2 are index1,1, and 2 respectively. The +# index for the range vector is defined by the integer following that +# vector. The y-limits of the plot are given by the last two numbers +# (use 'def') for the full default range of the delta vector +list dpairs + # v1 dv v2 range + X_PSF X_PSF_SIG Xraw Mraw 2 -10.0 10.0 + Y_PSF Y_PSF_SIG Yraw Mraw 2 -10.0 10.0 + PSF_INST_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + PSF_INST_FLUX PSF_INST_FLUX_SIG Fraw Mraw 2 -10.0 10.0 + AP_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + AP_MAG_RAW PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 +end + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +macro ckradialflux + data $1 + read -fits Chip.xrad X_APER Y_APER PSF_FWHM APER_FLUX APER_FLUX_ERR APER_FLUX_STDEV APER_FILL + read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG + set mag = PSF_INST_MAG + # XXX include nradii in header + for i 0 11 + set mA$i = -2.5*log(APER_FLUX:$i) + set dM$i = mA$i - PSF_INST_MAG + vstat dM$i + $DM$i = $MEDIAN + end + for i 1 11 + $j = $i - 1 + echo {$DM$i - $DM$j} + end + + lim mag -2 2; clear; box + for i 0 11 + plot mag dM$i -c black + end +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end