Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 6511) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 6873) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,8 +13,16 @@ # include -# include "pmObjects.h" -# include "pmPSF.h" -# include "pmPSFtry.h" -# include "pmModelGroup.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmGrowthCurve.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmModelGroup.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" // ******** pmPSFtry functions ************************************************** @@ -33,5 +41,5 @@ psFree (test->psf); psFree (test->sources); - psFree (test->modelFLT); + psFree (test->modelEXT); psFree (test->modelPSF); psFree (test->metric); @@ -42,5 +50,5 @@ // allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) +pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors) { @@ -51,7 +59,7 @@ // XXX probably need to increment ref counter type = pmModelSetType (modelName); - test->psf = pmPSFAlloc (type); + test->psf = pmPSFAlloc (type, poissonErrors); test->sources = psMemIncrRefCounter(sources); - test->modelFLT = psArrayAlloc (sources->n); + test->modelEXT = psArrayAlloc (sources->n); test->modelPSF = psArrayAlloc (sources->n); test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); @@ -59,13 +67,7 @@ test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); - test->modelFLT->n = test->modelFLT->nalloc; - test->modelPSF->n = test->modelPSF->nalloc; - test->metric->n = test->metric->nalloc; - test->fitMag->n = test->fitMag->nalloc; - test->mask->n = test->mask->nalloc; - - for (int i = 0; i < test->modelFLT->n; i++) { + for (int i = 0; i < test->modelEXT->n; i++) { test->mask->data.U8[i] = 0; - test->modelFLT->data[i] = NULL; + test->modelEXT->data[i] = NULL; test->modelPSF->data[i] = NULL; test->metric->data.F64[i] = 0; @@ -86,5 +88,5 @@ // mask values indicate the reason the source was rejected: -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) +pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors) { bool status; @@ -93,8 +95,8 @@ float x; float y; - int Nflt = 0; + int Next = 0; int Npsf = 0; - pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName); + pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors); // stage 1: fit an independent model (freeModel) to all sources @@ -108,26 +110,24 @@ // set temporary object mask and fit object - // fit model as FLT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + // fit model as EXT, not PSF + + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); + status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED); // exclude the poor fits if (!status) { - psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; + psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL; psFree (model); continue; } - psfTry->modelFLT->data[i] = model; - Nflt ++; - } - psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); - psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); - - // make this optional? - // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat"); + psfTry->modelEXT->data[i] = model; + Next ++; + } + psLogMsg ("psphot.psftry", 4, "fit ext: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); + psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n); // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask); + pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask); // stage 3: refit with fixed shape parameters @@ -139,13 +139,13 @@ pmSource *source = psfTry->sources->data[i]; - pmModel *modelFLT = psfTry->modelFLT->data[i]; + pmModel *modelEXT = psfTry->modelEXT->data[i]; // set shape for this model based on PSF - pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf); + pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf); x = source->peak->x; y = source->peak->y; - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, modelPSF, true); + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); + status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF); // skip poor fits @@ -162,5 +162,9 @@ // XXX : use a different estimator for the local sky? // XXX : pass 'source' as input? - if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { + if (!pmSourcePhotometryModel (&fitMag, modelPSF)) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; + goto next_source; + } + if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) { psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; goto next_source; @@ -172,31 +176,55 @@ next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - } + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED); + + } + psfTry->psf->nPSFstars = Npsf; + psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); - // make this optional - // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat"); + // measure the chi-square trend as a function of flux (PAR[1]) + // this should be linear for Poisson errors and quadratic for constant sky errors + psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK); + + // write sources with models first + for (int i = 0; i < psfTry->sources->n; i++) { + pmModel *model = psfTry->modelPSF->data[i]; + if (model == NULL) { + flux->data.F64[i] = 0.0; + chisq->data.F64[i] = 0.0; + mask->data.U8[i] = 1; + } else { + flux->data.F64[i] = model->params->data.F32[1]; + chisq->data.F64[i] = model[0].chisq/model[0].nDOF; + mask->data.U8[i] = 0; + } + } + // use 3hi/3lo sigma clipping on the r2rflux vs metric fit + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + stats->clipSigma = 3.0; + stats->clipIter = 3; + // linear clipped fit of ApResid to r2rflux + psVectorClipFitPolynomial1D (psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux); + 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]); + } + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + psFree (stats); + psFree (flux); + psFree (chisq); // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric_Alt (psfTry, RADIUS); - 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); + pmPSFtryMetric (psfTry, RADIUS); + 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, float RADIUS) { - - float dBin; - int nKeep, nSkip; - // the measured (aperture - fit) magnitudes (dA == psfTry->metric) // depend on both the true ap-fit (dAo) and the bias in the sky measurement: @@ -208,197 +236,43 @@ // we use an outlier rejection to avoid this bias - // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - rflux->n = rflux->nalloc; + // r2rflux = radius^2 * ten(0.4*fitMag); + psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); for (int i = 0; i < psfTry->sources->n; i++) { if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); - } - - // find min and max of (1/flux): - psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL); - - // build binned versions of rflux, metric - dBin = (stats->max - stats->min) / 10.0; - psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64); - psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64); - psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8); - daBin->n = daBin->nalloc; - maskB->n = maskB->nalloc; - psFree (stats); - - psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin); - - // group data in daBin bins, measure lower 50% mean - for (int i = 0; i < daBin->n; i++) { - - psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - tmp->n = 0; - - // accumulate data within bin range - for (int j = 0; j < psfTry->sources->n; j++) { - // masked for: bad model fit, outlier in parameters - if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL) - continue; - - // skip points with extreme dA values - if (fabs(psfTry->metric->data.F64[j]) > 0.5) - continue; - - // skip points outside of this bin - if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) - continue; - if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) - continue; - - tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j]; - tmp->n ++; - } - - // is this a valid point? - maskB->data.U8[i] = 0; - if (tmp->n < 2) { - maskB->data.U8[i] = 1; - psFree (tmp); - continue; - } - - // dA values are contaminated with low outliers - // measure statistics only on upper 50% of points - // this would be easier if we could sort in reverse: - // - // psVectorSort (tmp, tmp); - // tmp->n = 0.5*tmp->n; - // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - // psVectorStats (stats, tmp, NULL, NULL, 0); - // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); - - psVectorSort (tmp, tmp); - nKeep = 0.5*tmp->n; - nSkip = tmp->n - nKeep; - - psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64); - tmp2->n = tmp2->nalloc; - for (int j = 0; j < tmp2->n; j++) { - tmp2->data.F64[j] = tmp->data.F64[j + nSkip]; - } - - stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, tmp2, NULL, NULL, 0); - psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); - - daBin->data.F64[i] = stats->sampleMedian; - - psFree (stats); - psFree (tmp); - psFree (tmp2); - } - - // linear fit to rfBin, daBin - psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD); - psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); - - // XXX EAM : this is the intended API (cycle 7? cycle 8?) - // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin); - - // XXX EAM : replace this when the above version is implemented - // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL); - - psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); - psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); - - stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); - stats = psVectorStats (stats, daResid, NULL, maskB, 1); - - psfTry->psf->ApResid = poly->coeff[0]; - psfTry->psf->dApResid = stats->clippedStdev; - psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); - - psFree (rflux); - psFree (rfBin); - psFree (daBin); - psFree (maskB); - psFree (daBinFit); - psFree (daResid); - psFree (poly); - psFree (stats); - psFree (fitstat); - - return true; -} - -bool pmPSFtryMetric_Alt (pmPSFtry *try - , float RADIUS) -{ - - // the measured (aperture - fit) magnitudes (dA == try->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 - - // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (try - ->sources->n, PS_TYPE_F64); - rflux->n = rflux->nalloc; - for (int i = 0; i < try - ->sources->n; i++) { - if (try - ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*try - ->fitMag->data.F64[i]); - } - - // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit + r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); + } + + // use 3hi/1lo sigma clipping on the r2rflux vs metric fit psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - // XXX EAM stats->min = 1.0; stats->max = 3.0; stats->clipIter = 3; - // linear clipped fit to rfBin, daBin - psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); - poly = psVectorClipFitPolynomial1D (poly, stats, try - ->mask, PSFTRY_MASK_ALL, try - ->metric, NULL, rflux); - fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - try - ->psf->ApResid = poly->coeff[0]; - try - ->psf->dApResid = stats->sampleStdev; - try - ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); - - fprintf (stderr, "*******************************************************************************\n"); - - FILE *f; - f = fopen ("apresid.dat", "w"); - if (f == NULL) - psAbort ("pmPSFtry", "can't open output file"); - - for (int i = 0; i < try - ->sources->n; i++) { - fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try - ->fitMag->data.F64[i], rflux->data.F64[i], try - ->metric->data.F64[i], try - ->mask->data.U8[i]); - } - fclose (f); - - psFree (rflux); + // 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 = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux); + psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS); + psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0]; + psfTry->psf->ApTrend->coeff[0][0][1][0] = 0; + + psfTry->psf->skyBias = poly->coeff[1]; + psfTry->psf->ApResid = poly->coeff[0]; + psfTry->psf->dApResid = stats->sampleStdev; + + psFree (r2rflux); psFree (poly); psFree (stats); - // psFree (daFit); - // psFree (daResid); - return true; } + +/* + (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) + (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) + (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) +*/ +