Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 15000) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 15002) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-24 21:27:58 $ + * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-25 04:17:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,4 +35,15 @@ #include "pmSourcePhotometry.h" +int testSaveImage (psMetadata *header, psImage *image, char *filename) { + + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + return (TRUE); +} + +bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *dz); +bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup); + // ******** pmPSFtry functions ************************************************** // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely @@ -172,5 +183,5 @@ if (!status) { psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; - psTrace ("psModules.pmPSFtry", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); continue; } @@ -179,5 +190,5 @@ if (!pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, mark)) { psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; - psTrace ("psModules.pmPSFtry", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); continue; } @@ -352,4 +363,5 @@ psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); psVector *dz = NULL; @@ -366,8 +378,10 @@ x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + flux->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_I0]; // weight by the error on the source flux if (dz) { - dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0]; + // dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0] / source->modelEXT->params->data.F32[PM_PAR_I0]; + dz->data.F32[i] = 1.0; } } @@ -386,6 +400,4 @@ // 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); @@ -403,15 +415,56 @@ } - // 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? - 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); + pmTrend2D *trend = psf->params->data[PM_PAR_E0]; + if (trend->mode == PM_TREND_MAP) { + float errorFloor = 0.0; + float errorTotal = 0.0; + float errorTotalMin = FLT_MAX; + int entryMin = -1; + + // XXX do I need to store the initial value of psfTry->mask and reset? + + // check the fit residuals and increase Nx,Ny until the error is minimized + for (int i = 1; i < 10; i++) { + + if (!pmPSFParamTrend (psf, i, &errorFloor, &errorTotal, psfTry->mask, x, y, 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 + if (i == 1) { + for (int i = 0; i < psfTry->sources->n; i++) { + if (dz) { + dz->data.F32[i] = errorFloor; + } + } + } + + // store the resulting errorTotal values and the scales, redo the best + if (errorTotal < errorTotalMin) { + errorTotalMin = errorTotal; + entryMin = i; + } + } + if (entryMin == -1) { + psAbort ("failed on ApResid Trend"); + } + + // XXX catch error condition + pmPSFParamTrend (psf, entryMin, &errorFloor, &errorTotal, psfTry->mask, x, y, e0, e1, e2, dz); + } else { + // 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? + 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); + } } @@ -421,6 +474,6 @@ 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], + fprintf (f, "%f %f %f : %f %f %f : %f %f %f : %d\n", + x->data.F32[i], y->data.F32[i], flux->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]), @@ -437,4 +490,5 @@ // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) + // XXX 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) { @@ -459,8 +513,28 @@ } + pmTrend2D *trendE0 = psf->params->data[PM_PAR_E0]; + + // re-create the pmTrend2D using the scales defined for E0 + psImageBinning *binning = psImageBinningAlloc(); + if (trendE0->mode == PM_TREND_MAP) { + binning->nXruff = trend->map->map->numCols; + binning->nYruff = trend->map->map->numRows; + } else { + binning->nXruff = trend->poly->nX; + binning->nYruff = trend->poly->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[i]); + psf->params->data[i] = pmTrend2DNoImageAlloc (trendE0->mode, 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, e0, NULL)) { + 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); @@ -504,5 +578,212 @@ psFree (z); psFree (dz); + psFree (flux); return true; } +// XXX this function only works for MAP type -- extend to poly... +bool pmPSFParamTrend (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { + + int Nx, Ny; + + // set the map scale to match the aspect ratio + if (psf->fieldNx > psf->fieldNy) { + Nx = scale; + Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx)); + } else { + Ny = scale; + Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy)); + } + + // do we have enough sources for this fine of a grid? + if (x->n < 3*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 + 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); + + // 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); + } + psf->psfTrendStats->clipIter = nIter; // restore default setting + + # if (0) + pmTrend2D *trend; + trend = psf->params->data[PM_PAR_E0]; + testSaveImage (NULL, trend->map->map, "e0.fits"); + testSaveImage (NULL, trend->map->error, "e0d.fits"); + trend = psf->params->data[PM_PAR_E1]; + testSaveImage (NULL, trend->map->map, "e1.fits"); + testSaveImage (NULL, trend->map->error, "e1d.fits"); + trend = psf->params->data[PM_PAR_E2]; + testSaveImage (NULL, trend->map->map, "e2.fits"); + testSaveImage (NULL, trend->map->error, "e2d.fits"); + # endif + + // 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); + + // measure errors on the map pixels + pmTrend2D *trend; + + 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)); + + // 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); + pmPSFErrorFloor (errorFloor, dz, e0res, e1res, e2res, mask, nGroup); + + *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); + + 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); + + psFree (e0fit); + psFree (e1fit); + psFree (e2fit); + + psFree (e0res); + psFree (e1res); + psFree (e2res); + + return true; +} + +bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) { + + psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + + // measure the trend in bins with 10 values each; if < 10 total, use them all + int nBin = PS_MAX (dMag->n / nGroup, 1); + + // output vectors for ApResid trend + psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32); + + // use dMag to group the dMag and dap vectors + psVector *index = psVectorSortIndex (NULL, dMag); + + // subset vectors for dMag 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_U8); + + int n = 0; + for (int i = 0; i < dSo->n; i++) { + int j; + for (j = 0; (j < nGroup) && (n < dMag->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.U8[j] = mask->data.U8[N]; + } + 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); + psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + psStatsInit (statsS); + psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + psStatsInit (statsS); + psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff); + dEsquare += PS_SQR(statsS->robustStdev); + + dSo->data.F32[i] = sqrt(dEsquare); + } + psFree (dE0subset); + psFree (dE1subset); + psFree (dE2subset); + psFree (mkSubset); + + psStats *stats = psStatsAlloc (PS_STAT_MIN); + psVectorStats (stats, dSo, NULL, NULL, 0); + + *errorFloor = stats->min; + + psFree (stats); + psFree (index); + + psFree (dSo); + + psFree (statsS); + + return true; +}