Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 14936) +++ 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