Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 29935) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 30031) @@ -62,4 +62,6 @@ #include "pmSourceIO.h" +bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header); + bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file) { @@ -851,65 +853,13 @@ // read the psf clump data for each region + status = false; if (roAnalysis) { - int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN"); - if (!status) { - // read old-style psf clump data - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); - psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); - - psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); - if (!regionMD) { - regionMD = psMetadataAlloc(); - psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); - psFree (regionMD); - } - - // psf clump data - pmPSFClump psfClump; - - psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); - psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); - psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); - psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); - - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - } else { - psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); - - for (int i = 0; i < nRegions; i++) { - char key[10]; - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - - psMetadata *regionMD = psMetadataLookupPtr (&status, roAnalysis, regionName); - if (!regionMD) { - regionMD = psMetadataAlloc(); - psMetadataAddMetadata (roAnalysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); - psFree (regionMD); - } - - // psf clump data - pmPSFClump psfClump; - - snprintf (key, 9, "CLX_%03d", i); - psfClump.X = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLY_%03d", i); - psfClump.Y = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLDX_%03d", i); - psfClump.dX = psMetadataLookupF32 (&status, header, key); assert(status); - snprintf (key, 9, "CLDY_%03d", i); - psfClump.dY = psMetadataLookupF32 (&status, header, key); assert(status); - - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - } - } + status = pmPSFmodelReadPSFClump (roAnalysis, header); + if (!status) { + psMetadataAddS32 (roAnalysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 0); + } + } + if (!roAnalysis || !status) { + psWarning ("no PSF.CLUMP data available for PSF model"); } @@ -1123,3 +1073,66 @@ } -// XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped +bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header) { + + bool status = false;; + + int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN"); + if (!status) { + // read old-style psf clump data + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + + // psf clump data + pmPSFClump psfClump; + psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); if (!status) return false; + psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); if (!status) return false; + psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); if (!status) return false; + psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); if (!status) return false; + + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + } else { + for (int i = 0; i < nRegions; i++) { + char key[10]; + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + + // psf clump data + pmPSFClump psfClump; + + snprintf (key, 9, "CLX_%03d", i); + psfClump.X = psMetadataLookupF32 (&status, header, key); if (!status) return false; + snprintf (key, 9, "CLY_%03d", i); + psfClump.Y = psMetadataLookupF32 (&status, header, key); if (!status) return false; + snprintf (key, 9, "CLDX_%03d", i); + psfClump.dX = psMetadataLookupF32 (&status, header, key); if (!status) return false; + snprintf (key, 9, "CLDY_%03d", i); + psfClump.dY = psMetadataLookupF32 (&status, header, key); if (!status) return false; + + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + } + psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); + } + return true; +} Index: trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- trunk/psModules/src/objects/pmPSFtry.h (revision 29935) +++ trunk/psModules/src/objects/pmPSFtry.h (revision 30031) @@ -100,5 +100,5 @@ bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); -bool pmPSFtryMakePSF (pmPSFtry *psfTry); +bool pmPSFtryMakePSF (bool *goodFit, pmPSFtry *psfTry); bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); @@ -123,5 +123,5 @@ ); -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); +bool pmPSFFitShapeParams (bool *goodFit, pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); Index: trunk/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 29935) +++ trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 30031) @@ -50,5 +50,5 @@ Note: some of the array entries may be NULL (failed fits); ignore them. *****************************************************************************/ -bool pmPSFtryMakePSF (pmPSFtry *psfTry) +bool pmPSFtryMakePSF (bool *goodFit, pmPSFtry *psfTry) { PS_ASSERT_PTR_NON_NULL(psfTry, false); @@ -74,8 +74,14 @@ // fit the shape parameters (SXX, SYY, SXY) as a function of position - if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, srcMask)) { + if (!pmPSFFitShapeParams (goodFit, psf, psfTry->sources, x, y, srcMask)) { psFree(x); psFree(y); return false; + } + if (!goodFit) { + psWarning ("poor fit to PSF shape parameters for trend order %d, %d, skipping\n", psf->trendNx, psf->trendNy); + psFree(x); + psFree(y); + return true; } @@ -115,5 +121,5 @@ // the mask is carried from previous steps and updated with this operation // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' - if (!pmTrend2DFit (trend, srcMask, 0xff, x, y, z, NULL)) { + if (!pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, z, NULL)) { psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); psFree(x); @@ -122,4 +128,13 @@ return false; } + if (!goodFit) { + // if we do not get a good fit (but do not actually hit an error), + // tell the calling program to try something else + psWarning ("poor fit to PSF parameter %d for trend order %d, %d, skipping\n", i, psf->trendNx, psf->trendNy); + psFree(x); + psFree(y); + psFree(z); + return true; + } if (trend->mode == PM_TREND_MAP) { // p_psImagePrint (2, trend->map->map, "param N Before"); // XXX TEST: @@ -163,5 +178,5 @@ // fit the shape parameters using the supplied order (pmPSF->trendNx,trendNy) -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { +bool pmPSFFitShapeParams (bool *goodFit, pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { // we are doing a robust fit. after each pass, we drop points which are more deviant than @@ -219,5 +234,11 @@ trend = psf->params->data[PM_PAR_E0]; trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); + status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e0, NULL); + if (!goodFit) { + psFree (e0); + psFree (e1); + psFree (e2); + return true; + } mean = psStatsGetValue (trend->stats, meanOption); stdev = psStatsGetValue (trend->stats, stdevOption); @@ -228,5 +249,11 @@ trend = psf->params->data[PM_PAR_E1]; trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); + status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e1, NULL); + if (!goodFit) { + psFree (e0); + psFree (e1); + psFree (e2); + return true; + } mean = psStatsGetValue (trend->stats, meanOption); stdev = psStatsGetValue (trend->stats, stdevOption); @@ -237,5 +264,11 @@ trend = psf->params->data[PM_PAR_E2]; trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); + status &= pmTrend2DFit (goodFit, trend, srcMask, 0xff, x, y, e2, NULL); + if (!goodFit) { + psFree (e0); + psFree (e1); + psFree (e2); + return true; + } mean = psStatsGetValue (trend->stats, meanOption); stdev = psStatsGetValue (trend->stats, stdevOption); @@ -246,4 +279,7 @@ if (!status) { psError (PS_ERR_UNKNOWN, true, "failed to fit PSF shape params"); + psFree (e0); + psFree (e1); + psFree (e2); return false; } Index: trunk/psModules/src/objects/pmPSFtryModel.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryModel.c (revision 29935) +++ trunk/psModules/src/objects/pmPSFtryModel.c (revision 30031) @@ -136,9 +136,14 @@ // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation - if (!pmPSFtryMakePSF (psfTry)) { + bool goodFit = false; + if (!pmPSFtryMakePSF (&goodFit, psfTry)) { psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); psFree(psfTry); return NULL; } + if (!goodFit) { + psWarning ("poor psf fit for order %d, skipping\n", i); + continue; + } // stage 3: refit with fixed shape parameters, measure pmPSFtry->metric @@ -169,4 +174,10 @@ } psFree (srcMask); + + if (!minPSF) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a valid psf model from the sources"); + psFree(psfTry); + return NULL; + } // keep the ones matching the min systematic error: Index: trunk/psModules/src/objects/pmTrend2D.c =================================================================== --- trunk/psModules/src/objects/pmTrend2D.c (revision 29935) +++ trunk/psModules/src/objects/pmTrend2D.c (revision 30031) @@ -179,5 +179,5 @@ } -bool pmTrend2DFit(pmTrend2D *trend, psVector *mask, psVectorMaskType maskVal, const psVector *x, +bool pmTrend2DFit(bool *goodFit, pmTrend2D *trend, psVector *mask, psVectorMaskType maskVal, const psVector *x, const psVector *y, const psVector *f, const psVector *df) { @@ -188,5 +188,13 @@ PS_ASSERT_VECTOR_NON_NULL(f, false); - bool status; + bool status = false; + *goodFit = false; + // for the psImageMap fit, it is possible to have valid data but no valid solution for + // example, an isolated cell may not be reached from other cells, making the solution + // degenerate. psImageMapFit should probably handle this case, but until it does, we allow + // it to fail on the result, but not yield an error (goodFit = false). + // psVectorClipFitPolynomial2D can not fail in this way (really?), so goodFit is always + // true + switch (trend->mode) { case PM_TREND_POLY_ORD: @@ -196,4 +204,5 @@ // of points in the image, and potentially based on the fractional range of the // data? + *goodFit = true; break; @@ -201,5 +210,5 @@ // XXX supply fraction from trend elements // XXX need to add the API which adjusts the scale - status = psImageMapClipFit(trend->map, trend->stats, mask, maskVal, x, y, f, df); + status = psImageMapClipFit(goodFit, trend->map, trend->stats, mask, maskVal, x, y, f, df); break; Index: trunk/psModules/src/objects/pmTrend2D.h =================================================================== --- trunk/psModules/src/objects/pmTrend2D.h (revision 29935) +++ trunk/psModules/src/objects/pmTrend2D.h (revision 30031) @@ -76,5 +76,6 @@ ); -bool pmTrend2DFit(pmTrend2D *trend, +bool pmTrend2DFit(bool *goodFit, + pmTrend2D *trend, psVector *mask, // Warning: mask is modified! psVectorMaskType maskVal,