Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 14937) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 15000) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:06:57 $ + * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -44,7 +44,5 @@ static void pmPSFtryFree (pmPSFtry *test) { - - if (test == NULL) - return; + if (test == NULL) return; psFree (test->psf); @@ -58,17 +56,10 @@ // allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask) +pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options) { - - // validate the requested model name - pmModelType type = pmModelClassGetType (modelName); - if (type == -1) { - psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName); - return NULL; - } - pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); - - test->psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); + psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); + + test->psf = pmPSFAlloc (options); test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); @@ -76,14 +67,15 @@ test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); + psVectorInit (test->mask, 0); + psVectorInit (test->metric, 0.0); + psVectorInit (test->metricErr, 0.0); + psVectorInit (test->fitMag, 0.0); + test->sources = psArrayAlloc (sources->n); + for (int i = 0; i < sources->n; i++) { test->sources->data[i] = pmSourceCopy (sources->data[i]); - test->mask->data.U8[i] = 0; - test->metric->data.F32[i] = 0; - test->metricErr->data.F32[i] = 0; - test->fitMag->data.F32[i] = 0; - } - - psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); + } + return (test); } @@ -99,5 +91,5 @@ // generate a pmPSFtry with a copy of the test PSF sources -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark) +pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark) { bool status; @@ -105,5 +97,12 @@ int Npsf = 0; - pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask); + // 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"); @@ -111,5 +110,5 @@ } - // stage 1: fit an independent model (freeModel) to all sources + // stage 1: fit an EXT model to all candidates PSF sources psTimerStart ("fit"); for (int i = 0; i < psfTry->sources->n; i++) { @@ -123,6 +122,6 @@ } - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark); + // set object mask to define valid pixels -- XXX not unset? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); // fit model as EXT, not PSF @@ -140,5 +139,5 @@ // stage 2: construct a psf (pmPSF) from this collection of model fits - if (!pmPSFFromPSFtry (psfTry, applyWeights)) { + if (!pmPSFFromPSFtry (psfTry)) { psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); psFree(psfTry); @@ -162,8 +161,8 @@ continue; } - source->modelPSF->radiusFit = RADIUS; - - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark); + source->modelPSF->radiusFit = options->radius; + + // set object mask to define valid pixels -- XXX not unset? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); // fit the PSF model to the source @@ -241,5 +240,5 @@ // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric (psfTry, RADIUS); + pmPSFtryMetric (psfTry, options->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); @@ -339,5 +338,5 @@ /***************************************************************************** pmPSFFromPSFtry (psfTry): build a PSF model from a collection of -source->modelEXT entires. The PSF ignores the first 4 (independent) model +source->modelEXT entries. The PSF ignores the first 4 (independent) model parameters and constructs a polynomial fit to the remaining as a function of image coordinate. @@ -345,5 +344,5 @@ Note: some of the array entries may be NULL (failed fits); ignore them. *****************************************************************************/ -bool pmPSFFromPSFtry (pmPSFtry *psfTry, bool applyWeight) +bool pmPSFFromPSFtry (pmPSFtry *psfTry) { pmPSF *psf = psfTry->psf; @@ -355,5 +354,5 @@ psVector *dz = NULL; - if (applyWeight) { + if (psf->poissonErrorsParams) { dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); } @@ -377,6 +376,4 @@ // more deviant than three sigma. // mask is currently updated for each pass. this is inconsistent? - - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very @@ -389,4 +386,6 @@ // threshold. + // XXX this function needs to check the fit residuals and modify Nx,Ny as needed + // convert the measured source shape paramters to polarization terms psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); @@ -406,20 +405,27 @@ // 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 < stats->clipIter; i++) { - psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y); - psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n); - psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y); - psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n); - psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y); - psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n); - } - - // XXX temporary dump of the psf parameters + 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? + pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz); + psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); + pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz); + psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); + pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz); + psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); + } + + // 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 : %d\n", + 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], psfTry->mask->data.U8[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]), + psfTry->mask->data.U8[i]); } fclose (f); @@ -455,8 +461,7 @@ // 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 'applyWeights' - if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) { + // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' + if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) { psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); - psFree(stats); psFree(x); psFree(y); @@ -485,6 +490,5 @@ for (int i = 0; i < psf->params->n; i++) { - if (psf->params->data[i] == NULL) - continue; + if (psf->params->data[i] == NULL) continue; fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]); } @@ -496,10 +500,9 @@ } - psFree (stats); psFree (x); psFree (y); psFree (z); psFree (dz); - return (true); + return true; }