Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 15002) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 15016) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-25 04:17:29 $ + * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-25 22:05:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,15 +35,4 @@ #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 @@ -363,10 +352,4 @@ 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; - if (psf->poissonErrorsParams) { - dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - } // construct the x,y terms @@ -378,117 +361,15 @@ 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] / source->modelEXT->params->data.F32[PM_PAR_I0]; - dz->data.F32[i] = 1.0; - } - } - - // 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? - - // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very - // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for - // each source and fit this set of parameters. These values are less tightly coupled, but - // are still inter-related. The fitted values do a good job of constraining the major axis - // and the position angle, but the minor axis is weakly measured. When we apply the PSF - // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape - // parameters, with the constraint that the minor axis must be greater than a minimum - // threshold. - - // convert the measured source shape paramters to polarization terms - psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *e1 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *e2 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelEXT == NULL) continue; - - psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); - - e0->data.F32[i] = pol.e0; - e1->data.F32[i] = pol.e1; - e2->data.F32[i] = pol.e2; - } - - 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); - } - } - - // 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 %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]), - 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); - } - - psFree (e0); - psFree (e1); - psFree (e2); + } + + if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { + psFree(x); + psFree(y); + psFree(z); + return false; + } // 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 + // 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) { @@ -508,27 +389,22 @@ for (int j = 0; j < psfTry->sources->n; j++) { pmSource *source = psfTry->sources->data[j]; - if (source->modelEXT == NULL) - continue; + if (source->modelEXT == NULL) continue; z->data.F32[j] = source->modelEXT->params->data.F32[i]; } - 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->nXruff = psf->trendNx; + binning->nYruff = psf->trendNy; binning->nXfine = psf->fieldNx; binning->nYfine = psf->fieldNy; - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - + + if (psf->psfTrendMode == PM_TREND_MAP) { + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); + } + + // free existing trend, re-alloc psFree (psf->params->data[i]); - psf->params->data[i] = pmTrend2DNoImageAlloc (trendE0->mode, binning, psf->psfTrendStats); + psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); psFree (binning); @@ -541,5 +417,4 @@ psFree(y); psFree(z); - psFree(dz); return false; } @@ -552,20 +427,16 @@ 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]); + if (source == NULL) continue; + if (source->modelEXT == NULL) continue; + + pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); + + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); for (int i = 0; i < psf->params->n; i++) { if (psf->params->data[i] == NULL) continue; - fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]); + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); } - fprintf (f, "%f %d\n", model->chisq, model->nIter); + fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); psFree(modelPSF); @@ -577,11 +448,156 @@ psFree (y); 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) { + +bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { + + // 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. + + // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very + // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for + // each source and fit this set of parameters. These values are less tightly coupled, but + // are still inter-related. The fitted values do a good job of constraining the major axis + // and the position angle, but the minor axis is weakly measured. When we apply the PSF + // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape + // parameters, with the constraint that the minor axis must be greater than a minimum + // threshold. + + // convert the measured source shape paramters to polarization terms + psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *mag = psVectorAlloc (sources->n, PS_TYPE_F32); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelEXT == NULL) continue; + // XXX I am relying on the fact that none of the masked sources + // have modelEXT set here. perhaps use the value of psfTry->mask instead? + + psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); + + e0->data.F32[i] = pol.e0; + e1->data.F32[i] = pol.e1; + e2->data.F32[i] = pol.e2; + + float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; + mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; + } + + if (psf->psfTrendMode == PM_TREND_MAP) { + float errorFloor = 0.0; + float errorTotal = 0.0; + float errorTotalMin = FLT_MAX; + int entryMin = -1; + + psVector *dz = NULL; + psVector *mask = psVectorAlloc (sources->n, PS_TYPE_U8); + + // check the fit residuals and increase Nx,Ny until the error is minimized + // pmPSFParamTrend increases the number along the longer of x or y + for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) { + + psVectorInit (mask, 0); + if (!pmPSFFitShapeParamsMap (psf, i, &errorFloor, &errorTotal, 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) { + // 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; + entryMin = i; + } + } + if (entryMin == -1) { + psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params"); + return false; + } + + // 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)) { + psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); + } + + pmTrend2D *trend = psf->params->data[PM_PAR_E0]; + psf->trendNx = trend->map->map->numCols; + psf->trendNy = trend->map->map->numRows; + + psFree (mask); + psFree (dz); + } else { + + // XXX iterate Nx, Ny based on scatter? + // XXX we force the x & y order to be the same + // modify the order to correspond to the actual number of matched stars: + int order = PS_MAX (psf->trendNx, psf->trendNy); + if ((sources->n < 15) && (order >= 3)) order = 2; + if ((sources->n < 11) && (order >= 2)) order = 1; + if ((sources->n < 8) && (order >= 1)) order = 0; + if ((sources->n < 3)) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + psf->trendNx = order; + psf->trendNy = order; + + // 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? + bool status = true; + status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); + status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); + status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL); + psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); + + if (!status) { + psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); + return false; + } + } + } + + // 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 : %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], + 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]), + srcMask->data.U8[i]); + } + fclose (f); + } + + psFree (e0); + psFree (e1); + psFree (e2); + psFree (mag); + return true; +} + +// 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) { int Nx, Ny; @@ -639,17 +655,4 @@ 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); @@ -700,5 +703,5 @@ // 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); + pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup); *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); @@ -718,18 +721,19 @@ } -bool pmPSFErrorFloor (float *errorFloor, psVector *dMag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) { +// calculate the minimum scatter of the parameters +bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, 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); + int nBin = PS_MAX (mag->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 + // use mag to group parameters in sequence + psVector *index = psVectorSortIndex (NULL, mag); + + // subset vectors for mag and dap values within the given range psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); @@ -740,5 +744,5 @@ for (int i = 0; i < dSo->n; i++) { int j; - for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) { + for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { int N = index->data.U32[n]; dE0subset->data.F32[j] = e0res->data.F32[N];