Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFFromPSFtry.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFFromPSFtry.c (revision 25455) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFFromPSFtry.c (revision 25455) @@ -0,0 +1,155 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +/***************************************************************************** +pmPSFFromPSFtry (psfTry): build a PSF model from a collection of source->modelEXT entries +using the specified order in X,Y. The PSF ignores the first 4 (independent) model +parameters and constructs a polynomial fit to the remaining as a function of image +coordinate. Input: psfTry with fitted source->modelEXT collection, pre-allocated psf +Note: some of the array entries may be NULL (failed fits); ignore them. + *****************************************************************************/ +bool pmPSFFromPSFtry (pmPSFtry *psfTry, int Nx, int Ny) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + pmPSF *psf = psfTry->psf; + + // construct the fit vectors from the collection of objects + psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + // construct the x,y terms + for (int i = 0; i < psfTry->sources->n; i++) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelEXT == NULL) + continue; + + x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + } + + if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { + psFree(x); + psFree(y); + psFree(z); + return false; + } + + // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) + // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters + for (int i = 0; i < psf->params->n; i++) { + switch (i) { + case PM_PAR_SKY: + case PM_PAR_I0: + case PM_PAR_XPOS: + case PM_PAR_YPOS: + case PM_PAR_SXX: + case PM_PAR_SYY: + case PM_PAR_SXY: + continue; + default: + break; + } + + // select the per-object fitted data for this parameter + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[j]; + if (source->modelEXT == NULL) continue; + z->data.F32[j] = source->modelEXT->params->data.F32[i]; + } + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = psf->trendNx; + binning->nYruff = psf->trendNy; + binning->nXfine = psf->fieldNx; + binning->nYfine = psf->fieldNy; + + if (psf->psfTrendMode == PM_TREND_MAP) { + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); + } + + // free existing trend, re-alloc + psFree (psf->params->data[i]); + psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); + psFree (binning); + + // fit the collection of measured parameters to the PSF 2D model + // 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 (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) { + psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); + psFree(x); + psFree(y); + psFree(z); + return false; + } + } + + // test dump of star parameters vs position (compare with fitted values) + if (psTraceGetLevel("psModules.objects") >= 4) { + FILE *f = fopen ("params.dat", "w"); + + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[j]; + if (source == NULL) continue; + if (source->modelEXT == NULL) continue; + + pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); + + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) continue; + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); + } + fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); + + psFree(modelPSF); + } + fclose (f); + } + + psFree (x); + psFree (y); + psFree (z); + return true; +} + Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (revision 25454) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (revision 25455) @@ -37,51 +37,4 @@ #include "pmSourceVisual.h" -bool printTrendMap (pmTrend2D *trend) { - - if (!trend->map) return false; - if (!trend->map->map) return false; - - for (int j = 0; j < trend->map->map->numRows; j++) { - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); - } - fprintf (stderr, "\t\t\t"); - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); - } - fprintf (stderr, "\n"); - } - return true; -} - -bool psImageMapCleanup (psImageMap *map) { - - if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true; - - // find the weighted average of all pixels - float Sum = 0.0; - float Wt = 0.0; - for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (!isfinite(map->map->data.F32[j][i])) continue; - Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i]; - Wt += map->error->data.F32[j][i]; - } - } - - float Mean = Sum / Wt; - - // do any of the pixels in the map need to be repaired? - // XXX for now, we are just replacing bad pixels with the Mean - for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (isfinite(map->map->data.F32[j][i]) && - (map->error->data.F32[j][i] > 0.0)) continue; - map->map->data.F32[j][i] = Mean; - } - } - return true; -} - // ******** pmPSFtry functions ************************************************** // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely @@ -170,123 +123,27 @@ // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { - - pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - if (!source->moments->nPixels) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - - source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); - if (source->modelEXT == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); - continue; - } - - // set object mask to define valid pixels - // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); - - // fit model as EXT, not PSF - status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); - - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // exclude the poor fits - if (!status) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); - continue; - } - Next ++; - } - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); - - if (Next == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); + if (!pmPSFtryFitEXT(psfTry, options, maskVal, markVal)) { + psError(PS_ERR_UNKNOWN, false, "failed to fit EXT models to sources for psf model"); psFree(psfTry); return NULL; - } - - // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation - if (!pmPSFFromPSFtry (psfTry)) { - psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); - psFree(psfTry); - return NULL; - } - - // stage 3: refit with fixed shape parameters - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { - - pmSource *source = psfTry->sources->data[i]; - - // masked for: bad model fit, outlier in parameters - if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); - continue; - } - - // set shape for this model based on PSF - source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); - if (source->modelPSF == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; - abort(); - continue; - } - source->modelPSF->radiusFit = options->radius; - - // XXXX use a different radius for the aperture magnitude than for the PSF fit? - - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); - - // fit the PSF model to the source - status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal); - - // skip poor fits - if (!status) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); - continue; - } - - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); - if (!status || isnan(source->apMag)) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); - continue; - } - - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - psfTry->fitMag->data.F32[i] = source->psfMag; - psfTry->metric->data.F32[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F32[i] = source->errMag; - - psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); - Npsf ++; - } - psfTry->psf->nPSFstars = Npsf; - - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); - - if (Npsf == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); - psFree(psfTry); - return NULL; - } + } + + for (int i = 0; i < Norder; i++) { + // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation + if (!pmPSFFromPSFtry (psfTry, Nx, Ny)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + + // stage 3: refit with fixed shape parameters, measure pmPSFtryMetric + if (!pmPSFtryFitPSF (psfTry, Nx, Ny)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + } + // XXXXX this is probably not used any more. Are the chisq of the fits so bad? can we + // fix them by softening the errors on the brightest pixels? // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) @@ -485,122 +342,4 @@ (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) */ - -/***************************************************************************** -pmPSFFromPSFtry (psfTry): build a PSF model from a collection of -source->modelEXT entries. The PSF ignores the first 4 (independent) model -parameters and constructs a polynomial fit to the remaining as a function of -image coordinate. - input: psfTry with fitted source->modelEXT collection, pre-allocated psf -Note: some of the array entries may be NULL (failed fits); ignore them. - *****************************************************************************/ -bool pmPSFFromPSFtry (pmPSFtry *psfTry) -{ - PS_ASSERT_PTR_NON_NULL(psfTry, false); - PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); - - pmPSF *psf = psfTry->psf; - - // construct the fit vectors from the collection of objects - psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - - // construct the x,y terms - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelEXT == NULL) - continue; - - x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; - y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; - } - - if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { - psFree(x); - psFree(y); - psFree(z); - return false; - } - - // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) - // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters - for (int i = 0; i < psf->params->n; i++) { - switch (i) { - case PM_PAR_SKY: - case PM_PAR_I0: - case PM_PAR_XPOS: - case PM_PAR_YPOS: - case PM_PAR_SXX: - case PM_PAR_SYY: - case PM_PAR_SXY: - continue; - default: - break; - } - - // select the per-object fitted data for this parameter - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source->modelEXT == NULL) continue; - z->data.F32[j] = source->modelEXT->params->data.F32[i]; - } - - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = psf->trendNx; - binning->nYruff = psf->trendNy; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - - if (psf->psfTrendMode == PM_TREND_MAP) { - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - } - - // free existing trend, re-alloc - psFree (psf->params->data[i]); - psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); - - // fit the collection of measured parameters to the PSF 2D model - // 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 (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) { - psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); - psFree(x); - psFree(y); - psFree(z); - return false; - } - } - - // test dump of star parameters vs position (compare with fitted values) - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("params.dat", "w"); - - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source == NULL) continue; - if (source->modelEXT == NULL) continue; - - pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); - - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); - - for (int i = 0; i < psf->params->n; i++) { - if (psf->params->data[i] == NULL) continue; - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); - } - fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); - - psFree(modelPSF); - } - fclose (f); - } - - psFree (x); - psFree (y); - psFree (z); - return true; -} bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { @@ -900,7 +639,7 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); psImageMapCleanup (trend->map); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); @@ -910,7 +649,7 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); psImageMapCleanup (trend->map); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); @@ -920,7 +659,7 @@ stdev = psStatsGetValue (trend->stats, stdevOption); psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); psImageMapCleanup (trend->map); - // printTrendMap (trend); + // pmTrend2DPrintMap (trend); pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); } Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.old.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.old.c (revision 25455) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.old.c (revision 25455) @@ -0,0 +1,1142 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// ******** pmPSFtry functions ************************************************** +// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely +// * fitted version of the model, the pmPSF fit to the fitted model parameters, +// * and the PSF fits to the source. It also includes the statistics from the +// * fits, both the individual sources, and the collection + +// free a pmPSFtry structure +static void pmPSFtryFree (pmPSFtry *test) +{ + if (test == NULL) return; + + psFree (test->psf); + psFree (test->sources); + psFree (test->metric); + psFree (test->metricErr); + psFree (test->fitMag); + psFree (test->mask); + return; +} + +// allocate a pmPSFtry based on the desired sources and the model (identified by name) +pmPSFtry *pmPSFtryAlloc (const psArray *sources, const pmPSFOptions *options) +{ + pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); + psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); + + test->psf = pmPSFAlloc (options); + test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); + test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); + test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F32); + test->mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK); + + psVectorInit (test->mask, 0); + psVectorInit (test->metric, 0.0); + psVectorInit (test->metricErr, 0.0); + psVectorInit (test->fitMag, 0.0); + + test->sources = psArrayAlloc (sources->n); + + for (int i = 0; i < sources->n; i++) { + test->sources->data[i] = pmSourceCopy (sources->data[i]); + } + + return (test); +} + +bool psMemCheckPSFtry(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmPSFtryFree); +} + + +// build a pmPSFtry for the given model: +// - fit each source with the free-floating model +// - construct the pmPSF from the collection of models +// - fit each source with the PSF-parameter models +// - measure the pmPSF quality metric (dApResid) + +// sources used in for pmPSFtry may be masked by the analysis +// mask values indicate the reason the source was rejected: + +// generate a pmPSFtry with a copy of the test PSF sources +pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) +{ + bool status; + int Next = 0; + int Npsf = 0; + + // validate the requested model name + options->type = pmModelClassGetType (modelName); + if (options->type == -1) { + psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName); + return NULL; + } + + pmPSFtry *psfTry = pmPSFtryAlloc (sources, options); + if (psfTry == NULL) { + psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model"); + return NULL; + } + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF + psTimerStart ("psf.fit"); + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + if (!source->moments) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + if (!source->moments->nPixels) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + + source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); + if (source->modelEXT == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); + continue; + } + + // set object mask to define valid pixels + // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + + // fit model as EXT, not PSF + status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); + + // clear object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // exclude the poor fits + if (!status) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); + continue; + } + Next ++; + } + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); + + if (Next == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); + psFree(psfTry); + return NULL; + } + + // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation + if (!pmPSFFromPSFtry (psfTry)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + + // stage 3: refit with fixed shape parameters + psTimerStart ("psf.fit"); + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); + continue; + } + + // set shape for this model based on PSF + source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); + if (source->modelPSF == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; + abort(); + continue; + } + source->modelPSF->radiusFit = options->radius; + + // XXXX use a different radius for the aperture magnitude than for the PSF fit? + + // set object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + + // fit the PSF model to the source + status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal); + + // skip poor fits + if (!status) { + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); + continue; + } + + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); + if (!status || isnan(source->apMag)) { + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); + continue; + } + + // clear object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + + psfTry->fitMag->data.F32[i] = source->psfMag; + psfTry->metric->data.F32[i] = source->apMag - source->psfMag; + psfTry->metricErr->data.F32[i] = source->errMag; + + psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); + Npsf ++; + } + psfTry->psf->nPSFstars = Npsf; + + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); + + if (Npsf == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); + psFree(psfTry); + return NULL; + } + + // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) + // this should be linear for Poisson errors and quadratic for constant sky errors + psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK); + + // generate the x and y vectors, and mask missing models + for (int i = 0; i < psfTry->sources->n; i++) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelPSF == NULL) { + flux->data.F32[i] = 0.0; + chisq->data.F32[i] = 0.0; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; + } else { + flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; + chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + } + } + + // use 3hi/3lo sigma clipping on the chisq fit + psStats *stats = options->stats; + + // linear clipped fit of chisq trend vs flux + if (options->chiFluxTrend) { + bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, + mask, 0xff, chisq, NULL, flux); + psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean + psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev + + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", + psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); + + psFree(flux); + psFree(mask); + psFree(chisq); + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); + psFree(psfTry); + return NULL; + } + } + + for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) { + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, + psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), + psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i)); + } + + // XXX this function wants aperture radius for pmSourcePhotometry + if (!pmPSFtryMetric (psfTry, options)) { + psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName); + psFree (psfTry); + return NULL; + } + + psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", + modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); + + return (psfTry); +} + +bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(options, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + float RADIUS = options->radius; + + // the measured (aperture - fit) magnitudes (dA == psfTry->metric) + // depend on both the true ap-fit (dAo) and the bias in the sky measurement: + // dA = dAo + dsky/flux + // where flux is the flux of the star + // we fit this trend to find the infinite flux aperture correction (dAo), + // the nominal sky bias (dsky), and the error on dAo + // the values of dA are contaminated by stars with close neighbors in the aperture + // we use an outlier rejection to avoid this bias + + // r2rflux = radius^2 * ten(0.4*fitMag); + psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) + continue; + r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]); + } + + // XXX test dump of aperture residual data + if (psTraceGetLevel("psModules.objects") >= 5) { + FILE *f = fopen ("apresid.dat", "w"); + for (int i = 0; i < psfTry->sources->n; i++) { + int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); + + pmSource *source = psfTry->sources->data[i]; + float x = source->peak->x; + float y = source->peak->y; + + fprintf (f, "%d %d, %f %f %f %f %f %f \n", + i, keep, x, y, + psfTry->fitMag->data.F32[i], + r2rflux->data.F32[i], + psfTry->metric->data.F32[i], + psfTry->metricErr->data.F32[i]); + } + fclose (f); + } + + // This analysis of the apResid statistics is only approximate. The fitted magnitudes + // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a + // function of magnitude. We re-measure the apResid statistics later in psphot using the + // linear, constant-error fitting. Do not reject outliers with excessive vigor here. + + // fit ApTrend only to r2rflux, ignore x,y,flux variations for now + // linear clipped fit of ApResid to r2rflux + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); + poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS) + + // XXX replace this with a psVectorStats call? since we are not fitting the trend + bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL, + psfTry->metric, psfTry->metricErr, r2rflux); + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly"); + + psFree(poly); + psFree(r2rflux); + + return false; + } + psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev + psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0], + psStatsGetValue(options->stats, stdevStat), psfTry->sources->n); + + float dSys = psVectorSystematicError (psfTry->metric, psfTry->metricErr, 0.1); + fprintf (stderr, "systematic error: %f\n", dSys); + + int n = 0; + psVector *bright = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + psVector *brightErr = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + for (int i = 0; i < psfTry->metric->n; i++) { + if (!isfinite(psfTry->metric->data.F32[i])) continue; + if (!isfinite(psfTry->metricErr->data.F32[i])) continue; + if (psfTry->metricErr->data.F32[i] <= 0.0) continue; + if (psfTry->metricErr->data.F32[i] > 0.005) continue; + bright->data.F32[n] = psfTry->metric->data.F32[i]; + brightErr->data.F32[n] = psfTry->metricErr->data.F32[i]; + n++; + } + bright->n = brightErr->n = n; + + float dSysBright = psVectorSystematicError (bright, brightErr, 0.1); + fprintf (stderr, "bright systematic error: %f\n", dSysBright); + psFree(bright); + psFree(brightErr); + + // XXX test dump of fitted model (dump when tracing?) + if (psTraceGetLevel("psModules.objects") >= 4) { + FILE *f = fopen ("resid.dat", "w"); + psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux); + for (int i = 0; i < psfTry->sources->n; i++) { + int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); + + pmSource *source = psfTry->sources->data[i]; + float x = source->peak->x; + float y = source->peak->y; + + fprintf (f, "%d %d, %f %f %f %f %f %f %f\n", + i, keep, x, y, + psfTry->fitMag->data.F32[i], + r2rflux->data.F32[i], + psfTry->metric->data.F32[i], + psfTry->metricErr->data.F32[i], + apfit->data.F32[i]); + } + fclose (f); + psFree (apfit); + } + + // XXX drop the skyBias value, or include above?? + psfTry->psf->skyBias = poly->coeff[1]; + psfTry->psf->ApResid = poly->coeff[0]; + psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat); + + psFree (r2rflux); + psFree (poly); + + return true; +} + +/* + (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) + (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) + (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) +*/ + +/***************************************************************************** +pmPSFFromPSFtry (psfTry): build a PSF model from a collection of +source->modelEXT entries. The PSF ignores the first 4 (independent) model +parameters and constructs a polynomial fit to the remaining as a function of +image coordinate. + input: psfTry with fitted source->modelEXT collection, pre-allocated psf +Note: some of the array entries may be NULL (failed fits); ignore them. + *****************************************************************************/ +bool pmPSFFromPSFtry (pmPSFtry *psfTry) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + pmPSF *psf = psfTry->psf; + + // construct the fit vectors from the collection of objects + psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + // construct the x,y terms + for (int i = 0; i < psfTry->sources->n; i++) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelEXT == NULL) + continue; + + x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + } + + if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { + psFree(x); + psFree(y); + psFree(z); + return false; + } + + // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) + // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters + for (int i = 0; i < psf->params->n; i++) { + switch (i) { + case PM_PAR_SKY: + case PM_PAR_I0: + case PM_PAR_XPOS: + case PM_PAR_YPOS: + case PM_PAR_SXX: + case PM_PAR_SYY: + case PM_PAR_SXY: + continue; + default: + break; + } + + // select the per-object fitted data for this parameter + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[j]; + if (source->modelEXT == NULL) continue; + z->data.F32[j] = source->modelEXT->params->data.F32[i]; + } + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = psf->trendNx; + binning->nYruff = psf->trendNy; + binning->nXfine = psf->fieldNx; + binning->nYfine = psf->fieldNy; + + if (psf->psfTrendMode == PM_TREND_MAP) { + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); + } + + // free existing trend, re-alloc + psFree (psf->params->data[i]); + psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); + psFree (binning); + + // fit the collection of measured parameters to the PSF 2D model + // 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 (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) { + psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); + psFree(x); + psFree(y); + psFree(z); + return false; + } + } + + // test dump of star parameters vs position (compare with fitted values) + if (psTraceGetLevel("psModules.objects") >= 4) { + FILE *f = fopen ("params.dat", "w"); + + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[j]; + if (source == NULL) continue; + if (source->modelEXT == NULL) continue; + + pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); + + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) continue; + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); + } + fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); + + psFree(modelPSF); + } + fclose (f); + } + + psFree (x); + psFree (y); + psFree (z); + return true; +} + +bool pmPSFFitShapeParams (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 + // three sigma. mask is currently updated for each pass. + + // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very + // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for + // each source and fit this set of parameters. These values are less tightly coupled, but + // are still inter-related. The fitted values do a good job of constraining the major axis + // and the position angle, but the minor axis is weakly measured. When we apply the PSF + // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape + // parameters, with the constraint that the minor axis must be greater than a minimum + // threshold. + + // XXX re-read the sextractor manual on handling 'infinitely thin' sources... + + // convert the measured source shape paramters to polarization terms + psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *mag = psVectorAlloc (sources->n, PS_TYPE_F32); + + for (int i = 0; i < sources->n; i++) { + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + pmSource *source = sources->data[i]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); + + e0->data.F32[i] = pol.e0; + e1->data.F32[i] = pol.e1; + e2->data.F32[i] = pol.e2; + + float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; + mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; + } + + if (psf->psfTrendMode == PM_TREND_MAP) { + float scatterTotal = 0.0; + float scatterTotalMin = FLT_MAX; + int entryMin = -1; + + psVector *dz = NULL; + psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK); + + // check the fit residuals and increase Nx,Ny until the error is minimized + // pmPSFParamTrend increases the number along the longer of x or y + for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) { + + // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) + for (int i = 0; i < mask->n; i++) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } + if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { + break; + } + + // store the resulting scatterTotal values and the scales, redo the best + if (scatterTotal < scatterTotalMin) { + scatterTotalMin = scatterTotal; + entryMin = i; + } + } + if (entryMin == -1) { + psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params"); + return false; + } + + // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) + for (int i = 0; i < mask->n; i++) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } + if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { + psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); + } + + pmTrend2D *trend = psf->params->data[PM_PAR_E0]; + psf->trendNx = trend->map->map->numCols; + psf->trendNy = trend->map->map->numRows; + + // copy mask back to srcMask + for (int i = 0; i < mask->n; i++) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } + + psFree (mask); + psFree (dz); + } else { + + // XXX iterate Nx, Ny based on scatter? + // XXX we force the x & y order to be the same + // modify the order to correspond to the actual number of matched stars: + int order = PS_MAX (psf->trendNx, psf->trendNy); + if ((sources->n < 15) && (order >= 3)) order = 2; + if ((sources->n < 11) && (order >= 2)) order = 1; + if ((sources->n < 8) && (order >= 1)) order = 0; + if ((sources->n < 3)) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + psf->trendNx = order; + psf->trendNy = order; + + // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. + // This way, the parameters masked by one of the fits will be applied to the others + for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { + + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + + trend = psf->params->data[PM_PAR_E0]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); + pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); + + trend = psf->params->data[PM_PAR_E1]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); + pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); + + trend = psf->params->data[PM_PAR_E2]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); + + if (!status) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + } + } + + // test dump of the psf parameters + if (psTraceGetLevel("psModules.objects") >= 4) { + FILE *f = fopen ("pol.dat", "w"); + fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); + for (int i = 0; i < e0->n; i++) { + fprintf (f, "%f %f : %f %f %f : %f %f %f : %d\n", + x->data.F32[i], y->data.F32[i], + e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], + pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), + pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), + pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); + } + fclose (f); + } + + psFree (e0); + psFree (e1); + psFree (e2); + psFree (mag); + return true; +} + +// fit the shape variations as a psImageMap for the given scale factor +bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { + + int Nx, Ny; + + // set the map scale to match the aspect ratio : for a scale of 1, we guarantee + // that we have a single cell + if (psf->fieldNx > psf->fieldNy) { + Nx = scale; + float AR = psf->fieldNy / (float) psf->fieldNx; + Ny = (int) (Nx * AR + 0.5); + Ny = PS_MAX (1, Ny); + } else { + Ny = scale; + float AR = psf->fieldNx / (float) psf->fieldNy; + Nx = (int) (Ny * AR + 0.5); + Nx = PS_MAX (1, Nx); + } + + // do we have enough sources for this fine of a grid? + if (x->n < 10*Nx*Ny) { + return false; + } + + // XXX check this against the exising type + pmTrend2DMode psfTrendMode = PM_TREND_MAP; + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = Nx; + binning->nYruff = Ny; + binning->nXfine = psf->fieldNx; + binning->nYfine = psf->fieldNy; + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); + + psFree (psf->params->data[PM_PAR_E0]); + psFree (psf->params->data[PM_PAR_E1]); + psFree (psf->params->data[PM_PAR_E2]); + + int nIter = psf->psfTrendStats->clipIter; + psf->psfTrendStats->clipIter = 1; + psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); + psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); + psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); + psFree (binning); + + // if the map is 1x1 (a single value), we measure the resulting ensemble scatter + + // if the map is more finely sampled, divide the values into two sets: measure the fit from + // one set and the scatter from the other set. + psVector *x_fit = NULL; + psVector *y_fit = NULL; + psVector *x_tst = NULL; + psVector *y_tst = NULL; + + psVector *e0obs_fit = NULL; + psVector *e1obs_fit = NULL; + psVector *e2obs_fit = NULL; + psVector *e0obs_tst = NULL; + psVector *e1obs_tst = NULL; + psVector *e2obs_tst = NULL; + + if (scale == 1) { + x_fit = psMemIncrRefCounter (x); + y_fit = psMemIncrRefCounter (y); + x_tst = psMemIncrRefCounter (x); + y_tst = psMemIncrRefCounter (y); + e0obs_fit = psMemIncrRefCounter (e0obs); + e1obs_fit = psMemIncrRefCounter (e1obs); + e2obs_fit = psMemIncrRefCounter (e2obs); + e0obs_tst = psMemIncrRefCounter (e0obs); + e1obs_tst = psMemIncrRefCounter (e1obs); + e2obs_tst = psMemIncrRefCounter (e2obs); + } else { + x_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + y_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + x_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + y_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); + for (int i = 0; i < e0obs_fit->n; i++) { + // e0obs->n == 8 or 9: + // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6 + // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7 + x_fit->data.F32[i] = x->data.F32[2*i+0]; + x_tst->data.F32[i] = x->data.F32[2*i+1]; + y_fit->data.F32[i] = y->data.F32[2*i+0]; + y_tst->data.F32[i] = y->data.F32[2*i+1]; + + e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; + e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; + e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0]; + e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1]; + e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0]; + e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1]; + } + } + + // the mask marks the values not used to calculate the ApTrend + psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK); + // copy mask values to fitMask as a starting point + for (int i = 0; i < fitMask->n; i++) { + fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } + + // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. + // This way, the parameters masked by one of the fits will be applied to the others + for (int i = 0; i < nIter; i++) { + // XXX we are using the same stats structure on each pass: do we need to re-init it? + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + + trend = psf->params->data[PM_PAR_E0]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); + // pmTrend2DPrintMap (trend); + psImageMapCleanup (trend->map); + // pmTrend2DPrintMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); + + trend = psf->params->data[PM_PAR_E1]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); + // pmTrend2DPrintMap (trend); + psImageMapCleanup (trend->map); + // pmTrend2DPrintMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); + + trend = psf->params->data[PM_PAR_E2]; + status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); + // pmTrend2DPrintMap (trend); + psImageMapCleanup (trend->map); + // pmTrend2DPrintMap (trend); + pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); + } + psf->psfTrendStats->clipIter = nIter; // restore default setting + + // construct the fitted values and the residuals + psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst); + psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst); + psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst); + + psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit); + psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit); + psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit); + + // measure scatter for the unfitted points + // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10); + // psTraceSetLevel ("psLib.math.vectorClippedStats", 10); + pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options)); + // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0); + // psTraceSetLevel ("psLib.math.vectorClippedStats", 0); + + psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny); + psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal); + + psFree (x_fit); + psFree (y_fit); + psFree (x_tst); + psFree (y_tst); + + psFree (e0obs_fit); + psFree (e1obs_fit); + psFree (e2obs_fit); + psFree (e0obs_tst); + psFree (e1obs_tst); + psFree (e2obs_tst); + + psFree (e0fit); + psFree (e1fit); + psFree (e2fit); + + psFree (e0res); + psFree (e1res); + psFree (e2res); + + // XXX copy fitMask values back to mask + for (int i = 0; i < fitMask->n; i++) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; + } + psFree (fitMask); + + return true; +} + +// calculate the scatter of the parameters +bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt) +{ + + // psStats *stats = psStatsAlloc(stdevOpt); + psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV); + + // calculate the root-mean-square of E0, E1, E2 + float dEsquare = 0.0; + psStatsInit (stats); + if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + psStatsInit (stats); + if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + psStatsInit (stats); + if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + *scatterTotal = sqrtf(dEsquare); + + psFree(stats); + return true; +} + +// calculate the minimum scatter of the parameters +bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, + psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt) +{ + + psStats *statsS = psStatsAlloc(stdevOpt); + + // measure the trend in bins with 10 values each; if < 10 total, use them all + int nBin = PS_MAX (mag->n / nGroup, 1); + + // use mag to group parameters in sequence + psVector *index = psVectorSortIndex (NULL, mag); + + // subset vectors for mag and dap values within the given range + psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); + psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); + psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); + psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK); + + int n = 0; + float min = INFINITY; // Minimum error + for (int i = 0; i < nBin; i++) { + int j; + int nValid = 0; + for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { + int N = index->data.U32[n]; + dE0subset->data.F32[j] = e0res->data.F32[N]; + dE1subset->data.F32[j] = e1res->data.F32[N]; + dE2subset->data.F32[j] = e2res->data.F32[N]; + + mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N]; + if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; + } + if (nValid < 3) continue; + + dE0subset->n = j; + dE1subset->n = j; + dE2subset->n = j; + mkSubset->n = j; + + // calculate the root-mean-square of E0, E1, E2 + float dEsquare = 0.0; + psStatsInit (statsS); + if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) { + } + dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); + + psStatsInit (statsS); + if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); + + psStatsInit (statsS); + if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); + + if (isfinite(dEsquare)) { + float err = sqrtf(dEsquare); + if (err < min) { + min = err; + } + } + } + *errorFloor = min; + + psFree (dE0subset); + psFree (dE1subset); + psFree (dE2subset); + psFree (mkSubset); + + psFree(index); + + psFree(statsS); + + return true; +} + +float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) { + + psAssert(residuals, "residuals cannot be NULL"); + psAssert(errors, "errors cannot be NULL"); + psAssert(residuals->n == errors->n, "residuals and errors must be the same length"); + + // given a vector of residuals and their formal errors, calculated the necessary systematic + // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction + // highest chi-square contributors (allowed outliers) + + psVector *mask = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK); + psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32); + + // calculate the chisq vector: + int Ngood = 0; + for (int i = 0; i < residuals->n; i++) { + chisq->data.F32[i] = PS_MAX_F32; + if (!isfinite(residuals->data.F32[i])) continue; + if (!isfinite(errors->data.F32[i])) continue; + if (errors->data.F32[i] <= 0.0) continue; + chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]); + Ngood ++; + } + + psVector *index = psVectorSortIndex(NULL, chisq); + + // toss out the clipFraction highest chisq values + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (i < (1.0 - clipFraction)*Ngood) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0; + } else { + mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1; + } + } + + // Ndof ~= Ngood + // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof + // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0 + + // use Newton-Raphson to solve for S2: + + // use median sigma to calculate the initial guess for S2: + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, errors, NULL, mask, 1); + float errorMedian = stats->sampleMedian; + + float nPts = 0.0; + float res2mean = 0.0; + float ChiSq = 0.0; + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; + res2mean += PS_SQR(residuals->data.F32[n]); + ChiSq += PS_SQR(residuals->data.F32[n] / errors->data.F32[n]); + nPts += 1.0; + } + res2mean /= nPts; + ChiSq /= nPts; + + float S2guess = res2mean - PS_SQR(errorMedian); + + psLogMsg ("psModules", 3, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", + ChiSq, residuals->n, Ngood, nPts, S2guess); + + for (int iter = 0; iter < 10; iter++) { + + ChiSq = 0.0; + float dRdS = 0.0; + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; + float error2 = PS_SQR(errors->data.F32[n]) + S2guess; + ChiSq += PS_SQR(residuals->data.F32[n]) / error2; + dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2); + } + ChiSq /= nPts; + dRdS /= nPts; + + // Note the sign on dS: dRdS above is -1 * dR/dS formally + float dS = (ChiSq - 1.0) / dRdS; + S2guess += dS; + + psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess); + } + + // free local allocations + psFree (mask); + psFree (chisq); + psFree (stats); + psFree (index); + + return (sqrt(S2guess)); +} + Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25455) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25455) @@ -0,0 +1,97 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// Fit an EXT model to all candidates PSF sources. +// Note: this is independent of the modeled 2D variations in the PSF. +bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { + + bool status; + + psTimerStart ("psf.fit"); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + int Next = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + if (!source->moments) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + if (!source->moments->nPixels) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + + source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); + if (source->modelEXT == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); + continue; + } + + // set object mask to define valid pixels + // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + + // fit model as EXT, not PSF + status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); + + // clear object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // exclude the poor fits + if (!status) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); + continue; + } + Next ++; + } + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); + + if (Next == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); + return false; + } + + return true; +} Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25455) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25455) @@ -0,0 +1,108 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// stage 3: Refit with fixed shape parameters. This function uses the LMM fitting, but could +// be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c) +bool pmPSFtryFitPSF (pmPSFtry *psfTry) { + + psTimerStart ("psf.fit"); + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + psAssert (source->moments, "how can a psf source not have moments?"); + + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); + continue; + } + + // set shape for this model based on PSF + source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); + if (source->modelPSF == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; + abort(); + continue; + } + source->modelPSF->radiusFit = options->radius; + // XXXX use a different radius for the aperture magnitude than for the PSF fit? + + // set object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + + // fit the PSF model to the source + status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_NORM, maskVal); + + // skip poor fits + if (!status) { + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); + continue; + } + + // XXX this function calculates the aperture magnitude, but we now use the moments->Sum as the flux + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); + if (!status || isnan(source->apMag)) { + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); + continue; + } + + // clear object mask to define valid pixels + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + + psfTry->fitMag->data.F32[i] = source->psfMag; + psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag; + psfTry->metricErr->data.F32[i] = source->errMag; + + psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); + Npsf ++; + } + psfTry->psf->nPSFstars = Npsf; + + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); + + if (Npsf == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); + psFree(psfTry); + return NULL; + } + Index: /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.c (revision 25454) +++ /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.c (revision 25455) @@ -298,2 +298,21 @@ return PM_TREND_NONE; } + +bool pmTrend2DPrintMap (pmTrend2D *trend) { + + if (!trend->map) return false; + if (!trend->map->map) return false; + + for (int j = 0; j < trend->map->map->numRows; j++) { + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); + } + fprintf (stderr, "\t\t\t"); + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); + } + fprintf (stderr, "\n"); + } + return true; +} + Index: /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.h (revision 25454) +++ /branches/eam_branches/20090715/psModules/src/objects/pmTrend2D.h (revision 25455) @@ -97,4 +97,6 @@ pmTrend2DMode pmTrend2DModeFromString(psString name); +bool pmTrend2DPrintMap (pmTrend2D *trend); + /// @} # endif