Index: /branches/eam_branches/20090715/psModules/src/extras/pmVisual.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/extras/pmVisual.c (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/extras/pmVisual.c (revision 25476) @@ -22,4 +22,6 @@ #include "pmSubtractionStamps.h" #include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmPSFtry.h" #include "pmFPAExtent.h" Index: /branches/eam_branches/20090715/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/Makefile.am (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/Makefile.am (revision 25476) @@ -50,4 +50,9 @@ pmPSF_IO.c \ pmPSFtry.c \ + pmPSFtryModel.c \ + pmPSFtryFitEXT.c \ + pmPSFtryMakePSF.c \ + pmPSFtryFitPSF.c \ + pmPSFtryMetric.c \ pmTrend2D.c \ pmGrowthCurveGenerate.c \ Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSF.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSF.h (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSF.h (revision 25476) @@ -38,6 +38,5 @@ * */ -typedef struct -{ +typedef struct { pmModelType type; ///< PSF Model in use psArray *params; ///< Model parameters (psPolynomial2D) @@ -65,6 +64,5 @@ pmGrowthCurve *growth; ///< apMag vs Radius pmResiduals *residuals; ///< normalized residual image (no spatial variation) -} -pmPSF; +} pmPSF; typedef struct { Index: anches/eam_branches/20090715/psModules/src/objects/pmPSFFromPSFtry.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFFromPSFtry.c (revision 25475) +++ (revision ) @@ -1,155 +1,0 @@ -/** @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 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (revision 25476) @@ -63,5 +63,5 @@ psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); - test->psf = pmPSFAlloc (options); + test->psf = NULL; test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); @@ -87,747 +87,4 @@ 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 - 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; - } - - 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]) - // 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) -*/ - -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; } @@ -886,5 +143,5 @@ 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]); + ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]); nPts += 1.0; } @@ -914,4 +171,5 @@ float dS = (ChiSq - 1.0) / dRdS; S2guess += dS; + S2guess = PS_MAX(0.0, S2guess); psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess); Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.h (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.h (revision 25476) @@ -89,5 +89,18 @@ * */ -pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType mark); +pmPSFtry *pmPSFtryModel ( + const psArray *sources, ///< PSF sources to use in the pmPSF model analysis + const char *modelName, ///< human-readable name of desired model + pmPSFOptions *options, + psImageMaskType maskVal, + psImageMaskType mark + ); + +/** fit EXT models to all possible psf sources */ +bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); + +bool pmPSFtryMakePSF (pmPSFtry *psfTry); + +bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); /** pmPSFtryMetric() @@ -97,8 +110,5 @@ * */ -bool pmPSFtryMetric( - pmPSFtry *psfTry, ///< Add comment. - pmPSFOptions *options ///< PSF fitting options -); +bool pmPSFtryMetric(pmPSFtry *psfTry); /** pmPSFtryMetric_Alt() @@ -112,4 +122,11 @@ float RADIUS ///< Add comment. ); + +bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); + +float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); + +/// @} +# endif /** @@ -125,13 +142,3 @@ * */ -bool pmPSFFromPSFtry (pmPSFtry *psfTry); -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); -bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz); -bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt); -bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt); - -float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); - -/// @} -# endif Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25476) @@ -61,5 +61,5 @@ } - source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); + source->modelEXT = pmSourceModelGuess (source, options->type); if (source->modelEXT == NULL) { psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; @@ -86,6 +86,6 @@ 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); + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n); if (Next == 0) { Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25476) @@ -39,7 +39,14 @@ // 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) { +bool pmPSFtryFitPSF (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 Npsf = 0; for (int i = 0; i < psfTry->sources->n; i++) { @@ -98,11 +105,12 @@ 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); + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, psfTry->sources->n); if (Npsf == 0) { psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); - psFree(psfTry); - return NULL; + return false; } + return true; +} Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c (revision 25476) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c (revision 25476) @@ -0,0 +1,252 @@ +/** @file pmPSFtry.c + * @brief generate a pmPSF from a collection of EXT measurments of likely PSF stars. + * + * @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 pmPSFtryMakePSF (pmPSFtry *psfTry) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + pmPSF *psf = psfTry->psf; + psVector *srcMask = psfTry->mask; + + // 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); + + // construct the x,y terms + for (int i = 0; i < psfTry->sources->n; i++) { + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + pmSource *source = psfTry->sources->data[i]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + } + + // fit the shape parameters (SXX, SYY, SXY) as a function of position + if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, srcMask)) { + psFree(x); + psFree(y); + return false; + } + + // vector to store the other parameter values + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY); + // fit 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++) { + // 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 = psfTry->sources->data[j]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + z->data.F32[j] = source->modelEXT->params->data.F32[i]; + } + + // 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], srcMask, 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; +} + +// fit the shape parameters using the supplied order (pmPSF->trendNx,trendNy) +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. the source mask (srcMask) is 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... + + // storage vectors for the polarization terms & mags + psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); + + // convert the measured source shape paramters to polarization terms + 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; + } + + // 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++) { + // XXX we are using the same stats structure on each pass: do we need to re-init it? + // XXX we hardwire this to SAMPLE stats above (psphotChoosePSF.c), hardwire here instead? + 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]; + 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); + 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); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + // pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); + + 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); + 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); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + // pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); + + 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); + 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); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + // pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); + + if (!status) { + psError (PS_ERR_UNKNOWN, true, "failed to fit PSF 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); + return true; +} + Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMetric.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMetric.c (revision 25476) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMetric.c (revision 25476) @@ -0,0 +1,84 @@ +/** @file pmPSFtry.c + * @brief: measure the systematic error in the aperture-psf magnitude + * + * @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" + +// The quality of the PSF model is determined by the systematic scatter in the fit - aperture +// magnitude. We use the moments->Sum value, calculated with a Gaussian window, as a proxy for +// the aperture magnitude. We measure the systmatic scatter with the function +// psVectorSystematicError which solves for the value of SysErr that needs to be added in +// quadrature to the errors of a set of residual measurements in order to yield a ChiSq of 1.0. +bool pmPSFtryMetric (pmPSFtry *psfTry) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + if (!psVectorStats (stats, psfTry->metric, NULL, psfTry->mask, PSFTRY_MASK_ALL)) { + psError(PS_ERR_UNKNOWN, false, "Failed to measured clipped mean"); + psFree (stats); + return false; + } + + // generate a residual vector + psVector *resid = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + psVector *error = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + int n = 0; + for (int i = 0; i < psfTry->metric->n; i++) { + if (psfTry->mask && (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] && PSFTRY_MASK_ALL)) continue; + // if (psfTry->metricErr->data.F32[i] > 0.005) continue; + resid->data.F32[n] = psfTry->metric->data.F32[i] - stats->sampleMedian; + error->data.F32[n] = psfTry->metricErr->data.F32[i]; + n++; + } + resid->n = error->n = n; + + float psfSysErr = psVectorSystematicError (resid, error, 0.05); + + psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f (systamtic error) from statistics of %ld psf stars (%d used)\n", + stats->sampleMedian, psfSysErr, psfTry->sources->n, n); + + psfTry->psf->ApResid = stats->sampleMedian; + psfTry->psf->dApResid = psfSysErr; + + pmSourceVisualPlotPSFMetric (psfTry); + + psFree (stats); + psFree (resid); + psFree (error); + + return true; +} Index: /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryModel.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryModel.c (revision 25476) +++ /branches/eam_branches/20090715/psModules/src/objects/pmPSFtryModel.c (revision 25476) @@ -0,0 +1,214 @@ +/** @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" + +// 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) +{ + // 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 + 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; + } + + int orderMin = (options->psfTrendMode == PM_TREND_MAP) ? 1 : 0; + int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy); + + // XXX set the min number of needed source more carefully (depends on psfTrendMode?) + if ((sources->n < 15) && (orderMax >= 3)) orderMax = 2; + if ((sources->n < 11) && (orderMax >= 2)) orderMax = 1; + if ((sources->n < 8) && (orderMax >= 1)) orderMax = 0; + if ((sources->n < 3)) { + psError (PS_ERR_UNKNOWN, true, "failed to determine PSF parameters"); + return NULL; + } + + // save the raw source mask (generated by pmPSFtryFitEXT) + psVector *srcMask = psVectorCopy (NULL, psfTry->mask, PS_TYPE_VECTOR_MASK); + + // we will save the PSF with the best fit (min systematic error) + pmPSF *minPSF = NULL; + psVector *minMask = NULL; + + // as we loop over orders, we need to refer to the initial selection, but we modify the + // option values to match the current guess: save the max values here: + int Nx = options->psfTrendNx; + int Ny = options->psfTrendNy; + for (int i = orderMin; i <= orderMax; i++) { + + if (Nx > Ny) { + options->psfTrendNx = i; + options->psfTrendNy = PS_MAX (orderMin, (int)(i * (Ny / Nx) + 0.5)); + } else { + options->psfTrendNy = i; + options->psfTrendNx = PS_MAX (orderMin, (int)(i * (Nx / Ny) + 0.5)); + } + + // free existing data, if any + psFree(psfTry->psf); + psFree(psfTry->mask); + + // allocate a mask and a psf model using the current Nx,Ny order values; + psfTry->psf = pmPSFAlloc (options); + psfTry->mask = psVectorCopy (NULL, srcMask, PS_TYPE_VECTOR_MASK); + + // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation + if (!pmPSFtryMakePSF (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, measure pmPSFtry->metric + if (!pmPSFtryFitPSF (psfTry, options, maskVal, markVal)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + + // stage 4: measure systematic error in pmPSFtry->metric + if (!pmPSFtryMetric (psfTry)) { + psError(PS_ERR_UNKNOWN, false, "failed to measure systematic error of metric"); + psFree(psfTry); + return NULL; + } + + if (!minPSF) { + minPSF = psMemIncrRefCounter(psfTry->psf); + minMask = psMemIncrRefCounter(psfTry->mask); + } + + if (psfTry->psf->dApResid < minPSF->dApResid) { + psFree (minPSF); + minPSF = psMemIncrRefCounter(psfTry->psf); + psFree (minMask); + minMask = psMemIncrRefCounter(psfTry->mask); + } + } + psFree (srcMask); + + // keep the ones matching the min systematic error: + psFree (psfTry->psf); + psfTry->psf = psMemIncrRefCounter(minPSF); + psFree (psfTry->mask); + psfTry->mask = psMemIncrRefCounter(minMask); + + // 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]) + // 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)); + } + + 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); +} + Index: /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c (revision 25476) @@ -5,4 +5,6 @@ #include #include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmPSFtry.h" #include "pmSourceVisual.h" @@ -27,4 +29,86 @@ bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi); +bool pmSourceVisualPlotPSFMetric (pmPSFtry *psfTry) { + + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa1); + KapaInitGraph (&graphdata); + + psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32); + + graphdata.xmin = +32.0; + graphdata.xmax = -32.0; + graphdata.ymin = +32.0; + graphdata.ymax = -32.0; + + // construct the plot vectors + int n = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + x->data.F32[n] = psfTry->fitMag->data.F32[i]; + y->data.F32[n] = psfTry->metric->data.F32[i]; + dy->data.F32[n] = psfTry->metricErr->data.F32[i]; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + n++; + } + x->n = y->n = dy->n = n; + + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + + KapaSetLimits (kapa1, &graphdata); + + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + graphdata.etype |= 0x01; + + KapaPrepPlot (kapa1, n, &graphdata); + KapaPlotVector (kapa1, n, x->data.F32, "x"); + KapaPlotVector (kapa1, n, y->data.F32, "y"); + KapaPlotVector (kapa1, n, dy->data.F32, "dym"); + KapaPlotVector (kapa1, n, dy->data.F32, "dyp"); + + psFree (x); + psFree (y); + psFree (dy); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) { Index: /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h (revision 25475) +++ /branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h (revision 25476) @@ -18,4 +18,5 @@ bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask); +bool pmSourceVisualPlotPSFMetric (pmPSFtry *try); /// @}