Index: /trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileIO.c (revision 14999) +++ /trunk/psModules/src/camera/pmFPAfileIO.c (revision 15000) @@ -774,4 +774,5 @@ } pmFPAview *thisView = pmFPAAddSourceFromHeader (file->fpa, phu, file->format); + assert (thisView); // XXX we are having some trouble with input psf files not having the Cell and fpa names matching. psFree (thisView); psFree (phu); Index: /trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 14999) +++ /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 15000) @@ -274,6 +274,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -331,7 +331,6 @@ if (i == PM_PAR_XPOS) continue; if (i == PM_PAR_YPOS) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 14999) +++ /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 15000) @@ -328,6 +328,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -384,7 +384,6 @@ if (i == PM_PAR_XPOS) continue; if (i == PM_PAR_YPOS) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 14999) +++ /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 15000) @@ -358,6 +358,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -412,7 +412,6 @@ if (i == PM_PAR_XPOS) continue; if (i == PM_PAR_YPOS) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 14999) +++ /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 15000) @@ -351,6 +351,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -405,7 +405,6 @@ if (i == PM_PAR_XPOS) continue; if (i == PM_PAR_YPOS) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 14999) +++ /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 15000) @@ -343,6 +343,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -397,7 +397,6 @@ if (i == PM_PAR_XPOS) continue; if (i == PM_PAR_YPOS) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psModules/src/objects/pmPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSF.c (revision 14999) +++ /trunk/psModules/src/objects/pmPSF.c (revision 15000) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 02:46:25 $ + * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -41,26 +41,38 @@ /*****************************************************************************/ -/* DEFINE STATEMENTS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* TYPE DEFINITIONS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ - -/*****************************************************************************/ /* FUNCTION IMPLEMENTATION - PUBLIC */ /*****************************************************************************/ + +static void pmPSFOptionsFree (pmPSFOptions *options) { + + if (!options) return; + + psFree (options->stats); + return; +} + +pmPSFOptions *pmPSFOptionsAlloc () { + + pmPSFOptions *options = (pmPSFOptions *) psAlloc(sizeof(pmPSFOptions)); + psMemSetDeallocator(options, (psFreeFunc) pmPSFOptionsFree); + + options->type = 0; + + options->stats = NULL; + + options->psfTrendMode = PM_TREND_NONE; + options->psfTrendNx = 0; + options->psfTrendNy = 0; + options->psfFieldNx = 0; + options->psfFieldNy = 0; + options->psfFieldXo = 0; + options->psfFieldYo = 0; + + options->poissonErrorsPhotLMM = true; + options->poissonErrorsPhotLin = false; + options->poissonErrorsParams = true; + + return options; +} /***************************************************************************** @@ -74,4 +86,5 @@ psFree (psf->ChiTrend); + psFree (psf->psfTrendStats); psFree (psf->ApTrend); psFree (psf->FluxScale); @@ -94,11 +107,12 @@ Object Normalization *****************************************************************************/ -pmPSF *pmPSFAlloc (pmModelType type, bool poissonErrors, psPolynomial2D *psfTrendMask) +pmPSF *pmPSFAlloc (pmPSFOptions *options) { int Nparams; pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF)); - - psf->type = type; + psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); + + psf->type = options->type; psf->chisq = 0.0; psf->ApResid = 0.0; @@ -108,5 +122,8 @@ psf->nPSFstars = 0; psf->nApResid = 0; - psf->poissonErrors = poissonErrors; + + psf->poissonErrorsPhotLMM = options->poissonErrorsPhotLMM; + psf->poissonErrorsPhotLin = options->poissonErrorsPhotLin; + psf->poissonErrorsParams = options->poissonErrorsParams; // the ApTrend components are (x, y). It may be represented with a polynomial or with a @@ -121,5 +138,5 @@ psf->FluxScale = NULL; - if (psf->poissonErrors) { + if (psf->poissonErrorsPhotLMM) { psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); } else { @@ -133,5 +150,5 @@ psf->residuals = NULL; - Nparams = pmModelClassParameterCount (type); + Nparams = pmModelClassParameterCount (options->type); if (!Nparams) { psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); @@ -139,4 +156,7 @@ } psf->params = psArrayAlloc(Nparams); + + // save the trend stats on the psf for use in pmPSFFromPSFtry + psf->psfTrendStats = psMemIncrRefCounter (options->stats); // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial @@ -147,27 +167,30 @@ // PM_PAR_XPOS, etc) - if (psfTrendMask) { + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = options->psfTrendNx; + binning->nYruff = options->psfTrendNy; + binning->nXfine = options->psfFieldNx; + binning->nYfine = options->psfFieldNy; + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo); + + psf->fieldNx = options->psfFieldNx; + psf->fieldNy = options->psfFieldNy; + psf->fieldXo = options->psfFieldXo; + psf->fieldYo = options->psfFieldYo; + + // define the parameter trends + if (options->psfTrendMode != PM_TREND_NONE) { for (int i = 0; i < psf->params->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; - - psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY); - for (int nx = 0; nx < param->nX + 1; nx++) { - for (int ny = 0; ny < param->nY + 1; ny++) { - param->mask[nx][ny] = psfTrendMask->mask[nx][ny]; - } - } - psf->params->data[i] = param; + if (i == PM_PAR_SKY) continue; + if (i == PM_PAR_I0) continue; + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + + psf->params->data[i] = pmTrend2DNoImageAlloc (options->psfTrendMode, binning, options->stats); } } - - psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); - return(psf); + psFree (binning); + return psf; } @@ -300,7 +323,11 @@ va_start(ap, sxy); - pmModelType type = pmModelClassGetType (typeName); - psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0); - pmPSF *psf = pmPSFAlloc (type, true, psfTrend); + pmPSFOptions *options = pmPSFOptionsAlloc (); + options->type = pmModelClassGetType (typeName); + options->psfTrendMode = PM_TREND_POLY_ORD; + options->psfTrendNx = 0; + options->psfTrendNy = 0; + + pmPSF *psf = pmPSFAlloc (options); psVector *par = psVectorAlloc (psf->params->n, PS_TYPE_F32); @@ -311,23 +338,24 @@ psEllipsePol pol = pmPSF_ModelToFit(par->data.F32); + pmTrend2D *trend = NULL; + // set the psf shape parameters - psPolynomial2D *poly = NULL; - poly = psf->params->data[PM_PAR_E0]; - poly->coeff[0][0] = pol.e0; - - poly = psf->params->data[PM_PAR_E1]; - poly->coeff[0][0] = pol.e1; - - poly = psf->params->data[PM_PAR_E2]; - poly->coeff[0][0] = pol.e2; + trend = psf->params->data[PM_PAR_E0]; + trend->poly->coeff[0][0] = pol.e0; + + trend = psf->params->data[PM_PAR_E1]; + trend->poly->coeff[0][0] = pol.e1; + + trend = psf->params->data[PM_PAR_E2]; + trend->poly->coeff[0][0] = pol.e2; for (int i = PM_PAR_SXY + 1; i < psf->params->n; i++) { - poly = psf->params->data[i]; - poly->coeff[0][0] = (psF32)va_arg(ap, psF64); + trend = psf->params->data[i]; + trend->poly->coeff[0][0] = (psF32)va_arg(ap, psF64); } va_end(ap); psFree (par); - psFree (psfTrend); + psFree (options); return psf; } Index: /trunk/psModules/src/objects/pmPSF.h =================================================================== --- /trunk/psModules/src/objects/pmPSF.h (revision 14999) +++ /trunk/psModules/src/objects/pmPSF.h (revision 15000) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:05:35 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -19,17 +19,4 @@ // type of model carried by the pmModel structure typedef int pmModelType; - -typedef enum { - PM_PSF_APTREND_ERROR = 0, - PM_PSF_APTREND_NONE, - PM_PSF_APTREND_CONSTANT, - PM_PSF_APTREND_SKYBIAS, - PM_PSF_APTREND_SKYSAT, - PM_PSF_APTREND_XY_LIN, - PM_PSF_APTREND_XY_QUAD, - PM_PSF_APTREND_SKY_XY_LIN, - PM_PSF_APTREND_SKYSAT_XY_LIN, - PM_PSF_APTREND_ALL -} pmPSFApTrendOptions; /** pmPSF data structure @@ -49,4 +36,5 @@ pmModelType type; ///< PSF Model in use psArray *params; ///< Model parameters (psPolynomial2D) + psStats *psfTrendStats; ///< psf parameter trend clipping stats psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) pmTrend2D *ApTrend; ///< ApResid vs (x,y) @@ -59,9 +47,31 @@ int nPSFstars; ///< number of stars used to measure PSF int nApResid; ///< number of stars used to measure ApResid - bool poissonErrors; + int fieldNx; + int fieldNy; + int fieldXo; + int fieldYo; + bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting + bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting + bool poissonErrorsParams; ///< use poission errors for model parameter fitting pmGrowthCurve *growth; ///< apMag vs Radius pmResiduals *residuals; ///< normalized residual image (no spatial variation) } pmPSF; + +typedef struct { + pmModelType type; + psStats *stats; // psfTrend clipping stats + pmTrend2DMode psfTrendMode; + int psfTrendNx; + int psfTrendNy; + int psfFieldNx; + int psfFieldNy; + int psfFieldXo; + int psfFieldYo; + bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting + bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting + bool poissonErrorsParams; ///< use poission errors for model parameter fitting + float radius; +} pmPSFOptions; # define PM_PAR_E0 PM_PAR_SXX @@ -74,13 +84,6 @@ * */ -pmPSF *pmPSFAlloc( - pmModelType type, // type of model for PSF - bool poissonErrors, - psPolynomial2D *psfTrendMask -); - -bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option); - -pmPSFApTrendOptions pmPSFApTrendOptionFromName (char *name); +pmPSF *pmPSFAlloc (pmPSFOptions *options); +pmPSFOptions *pmPSFOptionsAlloc (); double pmPSF_SXYfromModel (psF32 *modelPar); @@ -91,7 +94,8 @@ bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes); +bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); + +psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR); -bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); -psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); /// @} Index: /trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- /trunk/psModules/src/objects/pmPSF_IO.c (revision 14999) +++ /trunk/psModules/src/objects/pmPSF_IO.c (revision 15000) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:05:57 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -148,8 +148,11 @@ // for a pmPSF supplied on the analysis metadata, we write out +// if needed: +// - a PHU blank header // - image header : FITS Image NAXIS = 0 -// - psf table (+header) : FITS Table -// - psf resid (+header) : FITS Image -// if needed, we also write out a PHU blank header +// if (trendMode == MAP) +// - psf resid (+header) : FITS Image +// else +// - psf table (+header) : FITS Table bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { @@ -159,4 +162,11 @@ if (!analysis) return false; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT"); + if (!recipe) { + psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT"); + return false; + } // write a PHU? (only if input image is MEF) @@ -215,5 +225,5 @@ // write out the IMAGE header segment - // this header block is new, write it to disk + // if this header block is new, write it to disk if (hdu->header != file->header) { // add EXTNAME, EXTHEAD, EXTTYPE to header @@ -250,25 +260,59 @@ psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName); - psMetadataAddBool (header, PS_LIST_TAIL, "POISSON", 0, "Use Poisson errors in fits?", psf->poissonErrors); + psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM", 0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM); + psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN", 0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin); + psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR", 0, "Use Poisson errors in fits?", psf->poissonErrorsParams); int nPar = pmModelClassParameterCount (psf->type) ; psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); + + psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx); + psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy); + psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1", 0, "Image X Ref", psf->fieldXo); + psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2", 0, "Image Y Ref", psf->fieldYo); + + // extract PSF Clump info + pmPSFClump psfClump; + + psfClump.X = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X"); assert (status); + psfClump.Y = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y"); assert (status); + psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX"); assert (status); + psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY"); assert (status); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X); + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y); + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX); + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY); // save the dimensions of each parameter for (int i = 0; i < nPar; i++) { char name[9]; - psPolynomial2D *poly = psf->params->data[i]; - if (poly == NULL) continue; + int nX, nY; + + pmTrend2D *trend = psf->params->data[i]; + if (trend == NULL) continue; + + if (trend->mode == PM_TREND_MAP) { + nX = trend->map->map->numCols; + nY = trend->map->map->numRows; + } else { + nX = trend->poly->nX; + nY = trend->poly->nY; + } snprintf (name, 9, "PAR%02d_NX", i); - psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX); + psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX); snprintf (name, 9, "PAR%02d_NY", i); - psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); + psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY); + snprintf (name, 9, "PAR%02d_MD", i); + char *modeName = pmTrend2DModeToString (trend->mode); + psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName); + psFree (modeName); } // other required information describing the PSF - psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); - psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); - psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); - psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); + psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); + psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); + psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); // XXX can we drop this now? @@ -278,19 +322,40 @@ psArray *psfTable = psArrayAllocEmpty (100); for (int i = 0; i < nPar; i++) { - psPolynomial2D *poly = psf->params->data[i]; - if (poly == NULL) continue; // skip unset parameters (eg, XPOS) - for (int ix = 0; ix <= poly->nX; ix++) { - for (int iy = 0; iy <= poly->nY; iy++) { - - psMetadata *row = psMetadataAlloc (); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); - psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); - psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); - psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); - psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); - - psArrayAdd (psfTable, 100, row); - psFree (row); + pmTrend2D *trend = psf->params->data[i]; + if (trend == NULL) continue; // skip unset parameters (eg, XPOS) + + if (trend->mode == PM_TREND_MAP) { + // write the image components into a table: this is needed because they may each be a different size + psImageMap *map = trend->map; + for (int ix = 0; ix < map->map->numCols; ix++) { + for (int iy = 0; iy < map->map->numRows; iy++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); + psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); + psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); + psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", map->map->data.F32[iy][ix]); + psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", map->error->data.F32[iy][ix]); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", 0); // no cells are masked + + psArrayAdd (psfTable, 100, row); + psFree (row); + } + } + } else { + // write the polynomial components into a table + psPolynomial2D *poly = trend->poly; + for (int ix = 0; ix <= poly->nX; ix++) { + for (int iy = 0; iy <= poly->nY; iy++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); + psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); + psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); + psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); + psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); + + psArrayAdd (psfTable, 100, row); + psFree (row); + } } } @@ -356,4 +421,5 @@ // XXX save the growth curve // XXX save ApTrend (as image?) + // XXX write the ApTrend with the same API as will be used for the PSF parameters above # if (0) @@ -474,4 +540,11 @@ psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename); + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT"); + if (!recipe) { + psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT"); + return false; + } + // Menu of EXTNAME rules psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); @@ -506,16 +579,42 @@ if (!header) psAbort("cannot read table header"); + pmPSFOptions *options = pmPSFOptionsAlloc(); + // load the PSF model parameters from the FITS table char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME"); - pmModelType type = pmModelClassGetType (modelName); - if (type == -1) { + options->type = pmModelClassGetType (modelName); + if (options->type == -1) { psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename); return false; } - bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON"); + // psf clump data + pmPSFClump psfClump; + + psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); + psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); + psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); + psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); + + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X" , 0, "psf clump center", psfClump.X); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y" , 0, "psf clump center", psfClump.Y); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump size", psfClump.dX); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump size", psfClump.dY); + + options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, header, "ERR_LMM"); + options->poissonErrorsPhotLin = psMetadataLookupBool (&status, header, "ERR_LIN"); + options->poissonErrorsParams = psMetadataLookupBool (&status, header, "ERR_PAR"); + + options->psfFieldNx = psMetadataLookupS32 (&status, header, "IMAXIS1"); + options->psfFieldNy = psMetadataLookupS32 (&status, header, "IMAXIS2"); + options->psfFieldXo = psMetadataLookupS32 (&status, header, "IMREF1"); + options->psfFieldYo = psMetadataLookupS32 (&status, header, "IMREF2"); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = options->psfFieldNx; + binning->nYfine = options->psfFieldNy; // we determine the PSF parameter polynomials from the MD-defined polynomials - pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL); + pmPSF *psf = pmPSFAlloc (options); // check the number of expected parameters @@ -524,18 +623,34 @@ psAbort("mismatch model par count"); - // load the dimensions of each parameter + // load the trend mode and dimensions of each parameter for (int i = 0; i < nPar; i++) { char name[9]; snprintf (name, 9, "PAR%02d_NX", i); - int nXorder = psMetadataLookupS32 (&status, header, name); + binning->nXruff = psMetadataLookupS32 (&status, header, name); if (!status) continue; // not all parameters are defined + snprintf (name, 9, "PAR%02d_NY", i); - int nYorder = psMetadataLookupS32 (&status, header, name); + binning->nYruff = psMetadataLookupS32 (&status, header, name); if (!status) { - psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i); + psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i); return false; } - psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); - } + + snprintf (name, 9, "PAR%02d_MD", i); + char *modeName = psMetadataLookupStr (&status, header, name); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i); + return false; + } + pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName); + if (psfTrendMode == PM_TREND_NONE) { + psfTrendMode = PM_TREND_POLY_ORD; + } + + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo); + psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL); + } + psFree (binning); // other required information describing the PSF @@ -557,15 +672,33 @@ int xPow = psMetadataLookupS32 (&status, row, "X_POWER"); int yPow = psMetadataLookupS32 (&status, row, "Y_POWER"); - // XXX sanity check here - - psPolynomial2D *poly = psf->params->data[iPar]; - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar); + + pmTrend2D *trend = psf->params->data[iPar]; + if (trend == NULL) { + psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar); return false; } - poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); - poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); - poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); + if (trend->mode == PM_TREND_MAP) { + psImageMap *map = trend->map; + assert (map); + assert (map->map); + assert (map->error); + assert (xPow >= 0); + assert (yPow >= 0); + assert (xPow < map->map->numCols); + assert (yPow < map->map->numRows); + map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); + map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); + } else { + psPolynomial2D *poly = trend->poly; + assert (poly); + assert (xPow >= 0); + assert (yPow >= 0); + assert (xPow <= poly->nX); + assert (yPow <= poly->nY); + poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); + poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); + poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); + } } psFree (header); @@ -618,91 +751,3 @@ // create a psMetadata representation (human-readable) of a psf model -psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf) -{ - - if (metadata == NULL) { - metadata = psMetadataAlloc (); - } - - char *modelName = pmModelClassGetName (psf->type); - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); - - int nPar = pmModelClassParameterCount (psf->type) ; - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); - - for (int i = 0; i < nPar; i++) { - psPolynomial2D *poly = psf->params->data[i]; - if (poly == NULL) - continue; - psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i); - } - - // XXX fix this - psWarning ("APTREND is currently missing"); - // psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND"); - - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); - - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); - psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors); - - return metadata; -} - -// parse a psMetadata representation (human-readable) of a psf model -pmPSF *pmPSFfromMetadata (psMetadata *metadata) -{ - - bool status; - char keyword[80]; - - char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME"); - pmModelType type = pmModelClassGetType (modelName); - - bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS"); - if (!status) - poissonErrors = true; - - // we determine the PSF parameter polynomials from the MD-defined polynomials - pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL); - - int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR"); - if (nPar != pmModelClassParameterCount (psf->type)) - psAbort("mismatch model par count"); - - // un-fitted terms, not in the Metadata, are left NULL - // XXX add a double-check of the expected number? - for (int i = 0; i < nPar; i++) { - sprintf (keyword, "PSF_PAR%02d", i); - psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); - if (!status) - continue; - psPolynomial2D *poly = psPolynomial2DfromMetadata (folder); - psFree (psf->params->data[i]); - psf->params->data[i] = poly; - } - - // load the APTREND data - // XXX fix this to work with pmTrend2D - psWarning ("APTREND is not being read"); - # if (0) - sprintf (keyword, "APTREND"); - psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); - psPolynomial4D *poly = psPolynomial4DfromMetadata (folder); - psFree (psf->ApTrend); - psf->ApTrend = poly; - # endif - - psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID"); - psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID"); - psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS"); - - psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ"); - psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS"); - - psFree (metadata); - return (psf); -} +// XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped Index: /trunk/psModules/src/objects/pmPSF_IO.h =================================================================== --- /trunk/psModules/src/objects/pmPSF_IO.h (revision 14999) +++ /trunk/psModules/src/objects/pmPSF_IO.h (revision 15000) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-24 00:11:02 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -33,7 +33,4 @@ bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip); -psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf); -pmPSF *pmPSFfromMetadata (psMetadata *metadata); - /// @} # endif Index: /trunk/psModules/src/objects/pmPSFtoMetadata.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtoMetadata.c (revision 15000) +++ /trunk/psModules/src/objects/pmPSFtoMetadata.c (revision 15000) @@ -0,0 +1,91 @@ + +// create a psMetadata representation (human-readable) of a psf model +psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf) +{ + + if (metadata == NULL) { + metadata = psMetadataAlloc (); + } + + char *modelName = pmModelClassGetName (psf->type); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); + + int nPar = pmModelClassParameterCount (psf->type) ; + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); + + for (int i = 0; i < nPar; i++) { + psPolynomial2D *poly = psf->params->data[i]; + if (poly == NULL) + continue; + psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i); + } + + // XXX fix this + psWarning ("APTREND is currently missing"); + // psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND"); + + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); + + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors); + + return metadata; +} + +// parse a psMetadata representation (human-readable) of a psf model +pmPSF *pmPSFfromMetadata (psMetadata *metadata) +{ + + bool status; + char keyword[80]; + + char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME"); + pmModelType type = pmModelClassGetType (modelName); + + bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS"); + if (!status) + poissonErrors = true; + + // we determine the PSF parameter polynomials from the MD-defined polynomials + pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL); + + int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR"); + if (nPar != pmModelClassParameterCount (psf->type)) + psAbort("mismatch model par count"); + + // un-fitted terms, not in the Metadata, are left NULL + // XXX add a double-check of the expected number? + for (int i = 0; i < nPar; i++) { + sprintf (keyword, "PSF_PAR%02d", i); + psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); + if (!status) + continue; + psPolynomial2D *poly = psPolynomial2DfromMetadata (folder); + psFree (psf->params->data[i]); + psf->params->data[i] = poly; + } + + // load the APTREND data + // XXX fix this to work with pmTrend2D + psWarning ("APTREND is not being read"); + # if (0) + sprintf (keyword, "APTREND"); + psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); + psPolynomial4D *poly = psPolynomial4DfromMetadata (folder); + psFree (psf->ApTrend); + psf->ApTrend = poly; + # endif + + psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID"); + psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID"); + psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS"); + + psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ"); + psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS"); + + psFree (metadata); + return (psf); +} Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 14999) +++ /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; } Index: /trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.h (revision 14999) +++ /trunk/psModules/src/objects/pmPSFtry.h (revision 15000) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-15 20:21:18 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -78,5 +78,5 @@ * */ -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask); +pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options); /** pmPSFtryModel() @@ -87,5 +87,5 @@ * */ -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); /** pmPSFtryMetric() @@ -123,5 +123,5 @@ * */ -bool pmPSFFromPSFtry (pmPSFtry *psfTry, bool applyWeights); +bool pmPSFFromPSFtry (pmPSFtry *psfTry); /// @} Index: /trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 14999) +++ /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 15000) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 02:46:46 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -81,6 +81,8 @@ // we must have a valid model model = pmSourceGetModel (&isPSF, source); - if (model == NULL) - return false; + if (model == NULL) { + psTrace ("psModules.objects", 3, "fail mag : no valid model"); + return false; + } if (model->dparams->data.F32[PM_PAR_I0] > 0) { @@ -102,5 +104,4 @@ source->psfMag = NAN; } - fprintf (stderr, "."); } else { status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); @@ -109,12 +110,14 @@ // measure EXT model photometry status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS - if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf) { - // convert to the equivalent 2D model? + if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) { source->psfMag += pmTrend2DEval (psf->ApTrend, x, y); } - if (!isfinite(SN) || (SN < AP_MIN_SN)) - return false; + if (!isfinite(SN) || (SN < AP_MIN_SN)) { + psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN); + return false; + } // replace source flux @@ -209,4 +212,5 @@ status = pmSourcePhotometryAper (&source->apMag, model, flux, source->maskObj, maskVal); if (!status) { + psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag"); psErrorCode last = psErrorCodeLast(); if (last == PM_ERR_PHOTOM) { Index: /trunk/psModules/src/objects/pmTrend2D.c =================================================================== --- /trunk/psModules/src/objects/pmTrend2D.c (revision 14999) +++ /trunk/psModules/src/objects/pmTrend2D.c (revision 15000) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:09:18 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * * Copyright 2004 Institute for Astronomy, University of Hawaii @@ -14,4 +14,5 @@ #endif +# include # include # include "pmTrend2D.h" @@ -31,5 +32,4 @@ { assert (image); - assert (stats); pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); @@ -44,4 +44,14 @@ case PM_TREND_POLY_ORD: trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend); + // set masking somehow + for (int nx = 0; nx < trend->poly->nX + 1; nx++) { + for (int ny = 0; ny < trend->poly->nY + 1; ny++) { + if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { + trend->poly->mask[nx][ny] = 1; + } else { + trend->poly->mask[nx][ny] = 0; + } + } + } break; @@ -69,5 +79,5 @@ } -pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats) +pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats) { pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); @@ -81,5 +91,56 @@ switch (mode) { case PM_TREND_POLY_ORD: + trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff); + // set masking somehow + for (int nx = 0; nx < trend->poly->nX + 1; nx++) { + for (int ny = 0; ny < trend->poly->nY + 1; ny++) { + if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { + trend->poly->mask[nx][ny] = 1; + } else { + trend->poly->mask[nx][ny] = 0; + } + } + } + break; + + case PM_TREND_POLY_CHEB: + trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff); + break; + + case PM_TREND_MAP: { + // binning defines the map scale relationship + trend->map = psImageMapNoImageAlloc (binning, stats); + break; + } + + default: + psAbort ("error"); + } + return (trend); +} + +pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats) +{ + pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); + psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree); + + trend->map = NULL; + trend->poly = NULL; + trend->stats = psMemIncrRefCounter (stats); + trend->mode = mode; + + switch (mode) { + case PM_TREND_POLY_ORD: trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend); + // set masking somehow + for (int nx = 0; nx < trend->poly->nX + 1; nx++) { + for (int ny = 0; ny < trend->poly->nY + 1; ny++) { + if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { + trend->poly->mask[nx][ny] = 1; + } else { + trend->poly->mask[nx][ny] = 0; + } + } + } break; @@ -111,4 +172,9 @@ bool status; + assert (trend); + assert (x); + assert (y); + assert (f); + switch (trend->mode) { case PM_TREND_POLY_ORD: @@ -136,5 +202,5 @@ double result; - if (!trend) return 0.0; + assert (trend); switch (trend->mode) { @@ -158,4 +224,8 @@ psVector *result; + assert (trend); + assert (x); + assert (y); + switch (trend->mode) { case PM_TREND_POLY_ORD: @@ -173,2 +243,45 @@ return result; } + +psString pmTrend2DModeToString (pmTrend2DMode mode) { + + psString name; + + switch (mode) { + case PM_TREND_NONE: + name = psStringCopy ("NONE"); + break; + case PM_TREND_POLY_ORD: + name = psStringCopy ("POLY_ORD"); + break; + case PM_TREND_POLY_CHEB: + name = psStringCopy ("POLY_CHEB"); + break; + case PM_TREND_MAP: + name = psStringCopy ("MAP"); + break; + default: + psAbort ("invalid mode %d", mode); + } + return name; +} + +pmTrend2DMode pmTrend2DModeFromString (psString name) { + + if (!name) return PM_TREND_NONE; + + if (!strcasecmp (name, "NONE")) { + return PM_TREND_NONE; + } + if (!strcasecmp (name, "POLY_ORD")) { + return PM_TREND_POLY_ORD; + } + if (!strcasecmp (name, "POLY_CHEB")) { + return PM_TREND_POLY_CHEB; + } + if (!strcasecmp (name, "MAP")) { + return PM_TREND_MAP; + } + psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name); + return PM_TREND_NONE; +} Index: /trunk/psModules/src/objects/pmTrend2D.h =================================================================== --- /trunk/psModules/src/objects/pmTrend2D.h (revision 14999) +++ /trunk/psModules/src/objects/pmTrend2D.h (revision 15000) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:09:18 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-24 21:27:58 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -17,4 +17,5 @@ typedef enum { + PM_TREND_NONE, PM_TREND_POLY_ORD, PM_TREND_POLY_CHEB, @@ -33,4 +34,6 @@ pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats); +pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats); + // allocate a pmTrend2D tied to an abstract field with size nXfield,nYfield pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats); @@ -41,4 +44,7 @@ psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y); +psString pmTrend2DModeToString (pmTrend2DMode mode); +pmTrend2DMode pmTrend2DModeFromString (psString name); + /// @} # endif Index: /trunk/psphot/doc/psfmodel.txt =================================================================== --- /trunk/psphot/doc/psfmodel.txt (revision 14999) +++ /trunk/psphot/doc/psfmodel.txt (revision 15000) @@ -1,2 +1,31 @@ + +2007.09.21 + + there are three places where we can choose to use errors in the fits or not: + + * non-linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLMM) + * linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLin) + * fitting of the 2D variations in the psf parameters (poissonErrorsParams) + * fitting of the 2D variations in the aperture residuals + +2007.09.20 + + I am upgrading the PSF model to allow the parameter variation to be + modeled with pmTrend2D instead of just polynomials. I am making a + list of places to modify the code: + +pmPSFAlloc : need a method beyong psfTrendMask to carry in the psf +options + +pmPSF_ModelToFit : no need to change these + +update pmPSFBuildSimple to set the parameters of the pmTrend, which +ever is used. + +pmPSFtry.c: some significant re-work! + + + +pmPSF_IO : need new functions to save / load the trend (psImages) 2006.10.27 Index: /trunk/psphot/src/models/pmModel_STRAIL.c =================================================================== --- /trunk/psphot/src/models/pmModel_STRAIL.c (revision 14999) +++ /trunk/psphot/src/models/pmModel_STRAIL.c (revision 15000) @@ -530,6 +530,6 @@ for (int i = 4; i < 7; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = psPolynomial2DEval (poly, out[2], out[3]); + pmTrend2D *trend = psf->params->data[i-4]; + out[i] = pmTrend2DEval (trend, out[2], out[3]); } return(true); @@ -554,7 +554,6 @@ for (int i = 0; i < psf->params->n; i++) { if (i == PM_PAR_SKY) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psphot/src/models/pmModel_TEST1.c =================================================================== --- /trunk/psphot/src/models/pmModel_TEST1.c (revision 14999) +++ /trunk/psphot/src/models/pmModel_TEST1.c (revision 15000) @@ -217,6 +217,6 @@ out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params->data[i]; - out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } } @@ -246,7 +246,6 @@ for (int i = 0; i < psf->params->n; i++) { if (i == PM_PAR_SKY) continue; - psPolynomial2D *poly = psf->params->data[i]; - assert (poly); - PAR[i] = psPolynomial2DEval(poly, Xo, Yo); + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); } Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 14999) +++ /trunk/psphot/src/psphot.h (revision 15000) @@ -47,5 +47,5 @@ bool psphotReplaceAll (psArray *sources, psMaskType maskVal); bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark); -bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup); +bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup); bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag); bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark); Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 14999) +++ /trunk/psphot/src/psphotApResid.c (revision 15000) @@ -1,4 +1,5 @@ # include "psphotInternal.h" +# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; } // measure the aperture residual statistics and 2D variations bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark) @@ -63,9 +64,9 @@ model = source->modelPSF; - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->mode & PM_SOURCE_MODE_FAIL) continue; - if (source->mode & PM_SOURCE_MODE_POOR) continue; + if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR"); + if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR"); + if (source->mode & PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND"); + if (source->mode & PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR"); + if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures @@ -73,9 +74,11 @@ if (!pmSourceMagnitudes (source, psf, photMode, maskVal, mark)) { Nskip ++; + psTrace ("psphot", 3, "skip : bad source mag"); continue; } - if (!isfinite(source->apMag) || !isfinite(source->apMag)) { + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { Nfail ++; + psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); continue; } @@ -87,4 +90,5 @@ if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) { Nfail ++; + psTrace ("psphot", 3, "fail : bad dap %f %f", dap, MAX_AP_OFFSET); continue; } @@ -99,4 +103,5 @@ dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + psVectorExtend (mag, 100, 1); psVectorExtend (mask, 100, 1); psVectorExtend (xPos, 100, 1); @@ -146,4 +151,5 @@ } + // XXX catch error condition psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag); @@ -152,6 +158,4 @@ psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit); psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32)); - - psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits"); if (psTraceGetLevel("psphot") >= 5) { @@ -170,4 +174,5 @@ float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + psf->ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center psf->dApResid = errorFloor; @@ -206,5 +211,5 @@ */ -bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup) { +bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) { psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); @@ -225,4 +230,5 @@ psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); + psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_U8); int n = 0; @@ -233,13 +239,15 @@ dMSubset->data.F32[j] = dMag->data.F32[N]; dASubset->data.F32[j] = dap->data.F32[N]; + mkSubset->data.U8[j] = mask->data.U8[N]; } dMSubset->n = j; dASubset->n = j; + mkSubset->n = j; psStatsInit (statsS); psStatsInit (statsM); - psVectorStats (statsS, dASubset, NULL, NULL, 0xff); - psVectorStats (statsM, dMSubset, NULL, NULL, 0xff); + psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff); + psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff); dSo->data.F32[i] = statsS->robustStdev; @@ -250,4 +258,5 @@ psFree (dMSubset); psFree (dASubset); + psFree (mkSubset); psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); @@ -286,4 +295,7 @@ } + // the mask marks the values not used to calculate the ApTrend + psVectorInit (mask, 0); + // XXX stats structure for use by ApTrend : make parameters user setable psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); @@ -305,5 +317,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); - psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, nGroup); + psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup); psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 14999) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 15000) @@ -11,39 +11,55 @@ sources = psArraySort (sources, pmSourceSortBySN); + // structure to store user options defining the psf + pmPSFOptions *options = pmPSFOptionsAlloc (); + + // load user options from the recipe. no need to check existence -- they are // array to store candidate PSF stars int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS"); - PS_ASSERT (status, NULL); + assert (status); float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS"); - PS_ASSERT (status, NULL); + assert (status); // supply the measured sky variance for optional constant errors (non-poissonian) float SKY_SIG = psMetadataLookupF32 (&status, recipe, "SKY_SIG"); - PS_ASSERT (status, NULL); + assert (status); // use poissonian errors or local-sky errors - bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); - PS_ASSERT (status, NULL); + options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LMM"); + assert (status); // use poissonian errors or local-sky errors - bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS"); - PS_ASSERT (status, NULL); + options->poissonErrorsPhotLin = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LIN"); + assert (status); + + // use poissonian errors or local-sky errors + options->poissonErrorsParams = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PARAMS"); + assert (status); // how to model the PSF variations across the field - psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK"); - PS_ASSERT (status, NULL); + options->psfTrendMode = pmTrend2DModeFromString (psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE")); + assert (status); + + options->psfTrendNx = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NX"); + assert (status); + + options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY"); + assert (status); // get the fixed PSF fit radius // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS - float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); - PS_ASSERT (status, NULL); - - pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), POISSON_ERRORS); - - psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md); - if (!psfTrendMask) { - psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); - return NULL; - } + options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + assert (status); + + options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + + // dimensions of the field for which the PSF is defined + options->psfFieldNx = readout->image->numCols; + options->psfFieldNy = readout->image->numRows; + options->psfFieldXo = readout->image->col0; + options->psfFieldYo = readout->image->row0; + + pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM); psArray *stars = psArrayAllocEmpty (sources->n); @@ -66,5 +82,4 @@ psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates"); psFree (stars); - psFree (psfTrendMask); return NULL; } @@ -91,5 +106,5 @@ psMetadataItem *item = psListGetAndIncrement (iter); char *modelName = item->data.V; - models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark); + models->data[i] = pmPSFtryModel (stars, modelName, options, maskVal, mark); } @@ -97,5 +112,4 @@ psFree (list); psFree (stars); - psFree (psfTrendMask); // select the best of the models @@ -128,11 +142,5 @@ if (psTraceGetLevel("psphot") >= 5) { for (int i = PM_PAR_SXX; i < try->psf->params->n; i++) { - psPolynomial2D *poly = try->psf->params->data[i]; - for (int nx = 0; nx <= poly->nX; nx++) { - for (int ny = 0; ny <= poly->nY; ny++) { - if (poly->mask[nx][ny]) continue; - fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny); - } - } + // XXX re-write the output or some other status info } } @@ -236,7 +244,7 @@ // use pmModelSub because modelFlux has not been generated assert (source->maskObj); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); + psImageKeepCircle (source->maskObj, x, y, options->radius, "OR", PM_MASK_MARK); pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + psImageKeepCircle (source->maskObj, x, y, options->radius, "AND", PS_NOT_U8(PM_MASK_MARK)); } @@ -283,7 +291,7 @@ pmSourcesWritePSFs (try->sources, "psfstars.dat"); pmSourcesWriteEXTs (try->sources, "extstars.dat", false); - psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); - psMetadataConfigWrite (psfData, "psfmodel.dat"); - psFree (psfData); + // XXX need alternative output function + // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); + // psMetadataConfigWrite (psfData, "psfmodel.dat"); psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n"); exit (0); @@ -306,4 +314,5 @@ psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + psFree (options); return (psf); } Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 14999) +++ /trunk/psphot/src/psphotRoughClass.c (revision 15000) @@ -1,8 +1,14 @@ # include "psphotInternal.h" +# define CHECK_STATUS(S,MSG) { \ + if (!status) { \ + psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \ + return false; \ + } } // 2006.02.02 : no leaks bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool havePSF, psMaskType maskSat) { - static pmPSFClump psfClump; + bool status; + static pmPSFClump psfClump; psTimerStart ("psphot"); @@ -11,27 +17,14 @@ // determine the PSF parameters from the source moment values psfClump = pmSourcePSFClump (sources, recipe); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X", 0, "psf clump center", psfClump.X); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y", 0, "psf clump center", psfClump.Y); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump center", psfClump.dX); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump center", psfClump.dY); } else { // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y - bool status_x, status_y, status_a; - float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); - float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); - float ANGLE = psMetadataLookupF32 (&status_a, recipe, "ANGLE"); - if (!status_x | !status_y | !status_a) { - psError(PSPHOT_ERR_CONFIG, false, "FWHM_X, FWHM_Y, or ANGLE not defined"); - return false; - } - - psEllipseAxes axes; - axes.major = FWHM_X / (2.0*sqrt(2.0*log(2.0))); - axes.minor = FWHM_Y / (2.0*sqrt(2.0*log(2.0))); - axes.theta = ANGLE; - psEllipseShape shape = psEllipseAxesToShape (axes); - - psfClump.X = shape.sx; - psfClump.Y = shape.sy; - psfClump.dX = 0.1*shape.sx; - psfClump.dY = 0.1*shape.sy; - // dX,dY are somewhat crudely defined, but only used to select PSF candidates. - // if we already have a PSF, this is not actually used... + psfClump.X = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X"); CHECK_STATUS (status, "PSF.CLUMP.X"); + psfClump.Y = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y"); CHECK_STATUS (status, "PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX"); CHECK_STATUS (status, "PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY"); CHECK_STATUS (status, "PSF.CLUMP.DY"); } Index: /trunk/psphot/src/psphotSourcePlots.c =================================================================== --- /trunk/psphot/src/psphotSourcePlots.c (revision 14999) +++ /trunk/psphot/src/psphotSourcePlots.c (revision 15000) @@ -20,8 +20,9 @@ int dY = 0; // height of row so far int NX = 20*DX; // full width of output image - int NY = 0; // total height of output image + int NY = DY; // total height of output image so far - // first, examine the PSF and SAT stars: - // - determine bounding boxes for summary image + // first, examine the PSF and SAT stars to set output image size: + // - add stamp widths until we exceed output image width, + // - then start a new row offset by max height for (int i = 0; i < sources->n; i++) { @@ -31,5 +32,8 @@ keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); - if (!keep) continue; + if (!keep) { + psTrace ("psphot", 4, "not plotting star %d: %x", i, source->mode); + continue; + } // how does this subimage get placed into the output image? @@ -57,7 +61,15 @@ } + if (NY == 0) { + psWarning ("no PSF or SAT stars to plot? skipping.\n"); + return false; + } + // allocate output image psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32); psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32); + + psImageInit (outpos, 0.0); + psImageInit (outsub, 0.0); int Xo = 0; // starting corner of next box @@ -142,5 +154,5 @@ psFree (outpos); psFree (outsub); - return (sources); + return true; }