Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 19906) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 20078) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-06 13:05:13 $ + * @version $Revision: 1.63 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-13 02:00:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,4 +35,51 @@ #include "pmSourcePhotometry.h" +bool printTrendMap (pmTrend2D *trend) { + + if (!trend->map) return false; + if (!trend->map->map) return false; + + for (int j = 0; j < trend->map->map->numRows; j++) { + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); + } + fprintf (stderr, "\t\t\t"); + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); + } + fprintf (stderr, "\n"); + } + return true; +} + +bool psImageMapCleanup (psImageMap *map) { + + if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true; + + // find the weighted average of all pixels + float Sum = 0.0; + float Wt = 0.0; + for (int j = 0; j < map->map->numRows; j++) { + for (int i = 0; i < map->map->numCols; i++) { + if (!isfinite(map->map->data.F32[j][i])) continue; + Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i]; + Wt += map->error->data.F32[j][i]; + } + } + + float Mean = Sum / Wt; + + // do any of the pixels in the map need to be repaired? + // XXX for now, we are just replacing bad pixels with the Mean + for (int j = 0; j < map->map->numRows; j++) { + for (int i = 0; i < map->map->numCols; i++) { + if (isfinite(map->map->data.F32[j][i]) && + (map->error->data.F32[j][i] > 0.0)) continue; + map->map->data.F32[j][i] = Mean; + } + } + return true; +} + // ******** pmPSFtry functions ************************************************** // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely @@ -120,5 +167,5 @@ maskVal |= markVal; - // stage 1: fit an EXT model to all candidates PSF sources + // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF psTimerStart ("fit"); for (int i = 0; i < psfTry->sources->n; i++) { @@ -158,5 +205,5 @@ } - // stage 2: construct a psf (pmPSF) from this collection of model fits + // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation if (!pmPSFFromPSFtry (psfTry)) { psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); @@ -558,7 +605,6 @@ if (psf->psfTrendMode == PM_TREND_MAP) { - float errorFloor = 0.0; - float errorTotal = 0.0; - float errorTotalMin = FLT_MAX; + float scatterTotal = 0.0; + float scatterTotalMin = FLT_MAX; int entryMin = -1; @@ -571,24 +617,11 @@ psVectorInit (mask, 0); - if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { + if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { break; } - // we do not have a good model for the error distributions on e0, e1, e2. - // use the bright end scatter from the constant fit to set the scale for the - // versions with 2D variation. potentially scale by poisson errors... - // if (i == 1) { - // XXX let's drop this for the moment: relies on valid result from errorFloor - if (0) { - // if (psf->poissonErrorsParams) do something else.. - dz = psVectorAlloc (sources->n, PS_TYPE_F32); - for (int i = 0; i < sources->n; i++) { - dz->data.F32[i] = errorFloor; - } - } - - // store the resulting errorTotal values and the scales, redo the best - if (errorTotal < errorTotalMin) { - errorTotalMin = errorTotal; + // store the resulting scatterTotal values and the scales, redo the best + if (scatterTotal < scatterTotalMin) { + scatterTotalMin = scatterTotal; entryMin = i; } @@ -601,5 +634,5 @@ // XXX supply the resulting mask values back to srcMask psVectorInit (mask, 0); - if (!pmPSFFitShapeParamsMap (psf, entryMin, &errorFloor, &errorTotal, mask, x, y, mag, e0, e1, e2, dz)) { + if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); } @@ -689,9 +722,10 @@ // fit the shape variations as a psImageMap for the given scale factor -bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { +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 + // 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; @@ -707,10 +741,7 @@ // do we have enough sources for this fine of a grid? - if (x->n < 3*Nx*Ny) { + if (x->n < 10*Nx*Ny) { return false; } - - // the mask marks the values not used to calculate the ApTrend - psVectorInit (mask, 0); // XXX check this against the exising type @@ -736,15 +767,68 @@ 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_U8); + psVectorInit (fitMask, 0); + // 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? - // pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz); - // psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n); - // pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz); - // psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n); - // pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz); - // psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n); - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); @@ -757,108 +841,62 @@ trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, mask, 0xff, x, y, e0obs, dz); + 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->n); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, mask, 0xff, x, y, e1obs, dz); + 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->n); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, mask, 0xff, x, y, e2obs, dz); + 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); - + // printTrendMap (trend); + psImageMapCleanup (trend->map); + // printTrendMap (trend); } psf->psfTrendStats->clipIter = nIter; // restore default setting // construct the fitted values and the residuals - psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x, y); - psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x, y); - psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], x, y); - - psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs, "-", (void *) e0fit); - psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs, "-", (void *) e1fit); - psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit); - -// XXX this code determines the formal error on the map, not the scatter -# if (0) - // measure errors on the map pixels - pmTrend2D *trend; - psStatsOptions meanOption; - psStatsOptions stdevOption; - float mapErrorSum = 0.0; - float mapError = 0.0; - - trend = psf->params->data[PM_PAR_E0]; - meanOption = psStatsMeanOption (trend->stats->options); - stdevOption = psStatsStdevOption (trend->stats->options); - mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); - mapErrorSum += mapError; - psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError)); - - trend = psf->params->data[PM_PAR_E1]; - meanOption = psStatsMeanOption (trend->stats->options); - stdevOption = psStatsStdevOption (trend->stats->options); - mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); - mapErrorSum += mapError; - psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError)); - - trend = psf->params->data[PM_PAR_E2]; - meanOption = psStatsMeanOption (trend->stats->options); - stdevOption = psStatsStdevOption (trend->stats->options); - mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); - mapErrorSum += mapError; - psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError)); - - trend = psf->params->data[PM_PAR_E0]; - float mapErrorSum = 0.0; - float mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } - } - psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError)); - mapErrorSum += mapError; - - trend = psf->params->data[PM_PAR_E1]; - mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } - } - psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError)); - mapErrorSum += mapError; - - trend = psf->params->data[PM_PAR_E2]; - mapError = 0.0; - for (int iy = 0; iy < trend->map->error->numRows; iy++) { - for (int ix = 0; ix < trend->map->error->numCols; ix++) { - mapError += PS_SQR (trend->map->error->data.F32[iy][ix]); - } - } - psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError)); - mapErrorSum += mapError; - - psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum)); -# endif - - // measure systematic errorFloor & systematic / photon scale factor - // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin - int nGroup = PS_MAX (3*Nx*Ny, 10); - pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup, - psStatsStdevOption(psf->psfTrendStats->options)); - - // XXX Bogus: *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); - *errorTotal = *errorFloor; - - psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup); - psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter floor: %f, error map total: %f\n", *errorFloor, *errorTotal); + psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], x_tst, y_tst); + psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], x_tst, y_tst); + psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], 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, 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); @@ -870,4 +908,33 @@ psFree (e2res); + psFree (fitMask); + + return true; +} + +// calculate the scatter of the parameters +bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, 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); + psVectorStats (stats, e0res, NULL, NULL, 0xff); + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + psStatsInit (stats); + psVectorStats (stats, e1res, NULL, NULL, 0xff); + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + psStatsInit (stats); + psVectorStats (stats, e2res, NULL, NULL, 0xff); + dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); + + *scatterTotal = sqrtf(dEsquare); + + psFree(stats); return true; }