Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 9885) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 9886) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-11-03 15:12:24 $ + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-11-07 09:14:00 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -46,6 +46,6 @@ psFree (test->psf); psFree (test->sources); - psFree (test->modelEXT); - psFree (test->modelPSF); + // psFree (test->modelEXT); + // psFree (test->modelPSF); psFree (test->metric); psFree (test->metricErr); @@ -69,19 +69,20 @@ test->psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); - test->sources = psMemIncrRefCounter(sources); - test->modelEXT = psArrayAlloc (sources->n); - test->modelPSF = psArrayAlloc (sources->n); test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F64); test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); - - for (int i = 0; i < test->modelEXT->n; i++) { + // test->modelEXT = psArrayAlloc (sources->n); + // test->modelPSF = psArrayAlloc (sources->n); + + test->sources = psArrayAlloc (sources->n); + for (int i = 0; i < sources->n; i++) { + test->sources->data[i] = pmSourceCopy (sources->data[i]); test->mask->data.U8[i] = 0; - test->modelEXT->data[i] = NULL; - test->modelPSF->data[i] = NULL; test->metric->data.F64[i] = 0; test->metricErr->data.F64[i] = 0; test->fitMag->data.F64[i] = 0; + // test->modelEXT->data[i] = NULL; + // test->modelPSF->data[i] = NULL; } @@ -99,9 +100,8 @@ // mask values indicate the reason the source was rejected: +// generate a pmPSFtry with a copy of the test PSF sources pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask) { bool status; - float obsMag; - float fitMag; float x; float y; @@ -111,5 +111,5 @@ pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask); if (psfTry == NULL) { - psError (PS_ERR_UNKNOWN, false, "failed to create allocate psf model"); + psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model"); return NULL; } @@ -120,6 +120,6 @@ pmSource *source = psfTry->sources->data[i]; - pmModel *model = pmSourceModelGuess (source, psfTry->psf->type); - if (model == NULL) { + source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); + if (source->modelEXT == NULL) { psError(PS_ERR_UNKNOWN, false, "failed to build model"); return NULL; @@ -132,5 +132,5 @@ psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); - status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT); + status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT); psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); @@ -138,8 +138,8 @@ if (!status) { psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL; - psFree (model); - continue; - } - psfTry->modelEXT->data[i] = model; + continue; + } + // XXX is this still needed? + // psfTry->modelEXT->data[i] = psMemIncrRefCounter (source->modelEXT); Next ++; } @@ -148,5 +148,5 @@ // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask); + pmPSFFromPSFtry (psfTry); // stage 3: refit with fixed shape parameters @@ -158,46 +158,35 @@ pmSource *source = psfTry->sources->data[i]; - pmModel *modelEXT = psfTry->modelEXT->data[i]; // set shape for this model based on PSF - pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf); + source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); + source->modelPSF->radiusFit = RADIUS; x = source->peak->x; y = source->peak->y; + // set the mask and fit the PSF model psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); - status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF); + status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); // skip poor fits if (!status) { psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; - psFree (modelPSF); - goto next_source; + continue; } // otherwise, save the resulting model - psfTry->modelPSF->data[i] = modelPSF; - - // I'm not using the pmSourceMagnitudes API, why? + // psfTry->modelPSF->data[i] = psMemIncrRefCounter (source->modelPSF); + // XXX : use a different aperture radius from the fit radius? - // XXX : use a different estimator for the local sky? - // XXX : pass 'source' as input? - if (!pmSourcePhotometryModel (&fitMag, modelPSF)) { + if (!pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP)) { 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; - } - - psfTry->metric->data.F64[i] = obsMag - fitMag; - psfTry->metricErr->data.F64[i] = modelPSF->dparams->data.F32[PM_PAR_I0] / modelPSF->params->data.F32[PM_PAR_I0]; - - psfTry->fitMag->data.F64[i] = fitMag; + continue; + } + + psfTry->fitMag->data.F64[i] = source->psfMag; + psfTry->metric->data.F64[i] = source->apMag - source->psfMag; + psfTry->metricErr->data.F64[i] = source->errMag; Npsf ++; - -next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); - } psfTry->psf->nPSFstars = Npsf; @@ -212,24 +201,25 @@ psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK); - // write sources with models first + // generate the x and y vectors, and mask missing models for (int i = 0; i < psfTry->sources->n; i++) { - pmModel *model = psfTry->modelPSF->data[i]; - if (model == NULL) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelPSF == 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[PM_PAR_I0]; - chisq->data.F64[i] = model[0].chisq/model[0].nDOF; + flux->data.F64[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; + chisq->data.F64[i] = source->modelPSF->chisq / source->modelPSF->nDOF; mask->data.U8[i] = 0; } } - // use 3hi/3lo sigma clipping on the r2rflux vs metric fit + + // use 3hi/3lo sigma clipping on the chisq 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 - psfTry->psf->ChiTrend = - psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux); + + // linear clipped fit of chisq trend vs flux + psfTry->psf->ChiTrend = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux); psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); @@ -246,5 +236,7 @@ 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 ("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)); } @@ -342,2 +334,101 @@ */ +/***************************************************************************** +pmPSFFromPSFtry (psfTry): build a PSF model from a collection of +source->modelEXT entires. 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) +{ + pmPSF *psf = psfTry->psf; + + // construct the fit vectors from the collection of objects + psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + + // 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; + + // use F64 for polynomial fitting + x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + } + + // 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. this is inconsistent? + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + // skip the unfitted parameters (X, Y, Io, Sky) + for (int i = 0; i < psf->params_NEW->n; i++) { + if (i == PM_PAR_SKY) + continue; + if (i == PM_PAR_I0) + continue; + if (i == PM_PAR_XPOS) + continue; + if (i == PM_PAR_YPOS) + continue; + + // select the per-object fitted data for this parameter + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelEXT == NULL) + continue; + + z->data.F64[j] = source->modelEXT->params->data.F32[i]; + dz->data.F64[j] = 1; // use the model-fitted error? or S/N? + + // for SXY, we actually fit SXY / (SXX^-2 + SYY^-2) + if (i == PM_PAR_SXY) { + z->data.F64[j] = pmPSF_SXYfromModel (source->modelEXT->params->data.F32); + } + } + psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, dz, x, y); + // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, psfTry->mask->n, i); + } + + // XXX test dump of star parameters vs position (compare with fitted values) + if (0) { + 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; + + pmModel *model = source->modelEXT; + if (model == NULL) + continue; + + pmModel *modelPSF = pmModelFromPSF (model, psf); + + fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) + continue; + fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]); + } + fprintf (f, "%f %d\n", model->chisq, model->nIter); + } + fclose (f); + } + + psFree (stats); + psFree (x); + psFree (y); + psFree (z); + psFree (dz); + return (true); +} +