Index: /branches/eam_branch_20070921/psModules/src/objects/pmPSF.c =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmPSF.c (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmPSF.c (revision 14980) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.28.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 21:45:00 $ + * @version $Revision: 1.28.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -45,4 +45,9 @@ static void pmPSFOptionsFree (pmPSFOptions *options) { + + if (!options) return; + + psFree (options->image); + psFree (options->stats); return; } @@ -54,8 +59,13 @@ options->type = 0; - options->poissonErrors = false; + + options->image = NULL; + options->stats = NULL; + options->psfTrendMode = PM_TREND_NONE; options->psfTrendNx = 0; options->psfTrendNy = 0; + + options->poissonErrors = false; } @@ -136,4 +146,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 @@ -144,8 +157,5 @@ // PM_PAR_XPOS, etc) - // XXX make this a user option - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - - // XXX define new structure to carry in the psf options + // define the parameter trends if (options->psfTrendMode != PM_TREND_NONE) { for (int i = 0; i < psf->params->n; i++) { @@ -155,11 +165,7 @@ if (i == PM_PAR_YPOS) continue; - // XXX need to set stats & image or else nXfine,nYfine - pmTrend2D *param = pmTrend2DAlloc (options->psfTrendMode, image, options->psfTrendNx, options->psfTrendNx, stats); - psf->params->data[i] = param; + psf->params->data[i] = pmTrend2DAlloc (options->psfTrendMode, options->image, options->psfTrendNx, options->psfTrendNx, options->stats); } } - - psFree (stats); return psf; } @@ -287,4 +293,5 @@ // generate a psf model of the requested type, with fixed shape +// XXX update this. pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...) { Index: /branches/eam_branch_20070921/psModules/src/objects/pmPSF.h =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmPSF.h (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmPSF.h (revision 14980) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.16.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 18:55:12 $ + * @version $Revision: 1.16.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -36,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) @@ -54,8 +55,10 @@ typedef struct { pmModelType type; - bool poissonErrors; + psImage *image; // image for which the PSF is defined + psStats *stats; // psfTrend clipping stats pmTrend2DMode psfTrendMode; int psfTrendNx; int psfTrendNy; + bool poissonErrors; } pmPSFOptions; @@ -75,8 +78,4 @@ ); -bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option); - -pmPSFApTrendOptions pmPSFApTrendOptionFromName (char *name); - double pmPSF_SXYfromModel (psF32 *modelPar); double pmPSF_SXYtoModel (psF32 *fittedPar); @@ -86,7 +85,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: /branches/eam_branch_20070921/psModules/src/objects/pmPSF_IO.c =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmPSF_IO.c (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmPSF_IO.c (revision 14980) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.22.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 18:55:12 $ + * @version $Revision: 1.22.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -148,9 +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 -// XXX if we are using psImageMap, this is also a FITS Image, with header data -// - 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) { @@ -216,5 +218,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 @@ -244,4 +246,5 @@ // write the PSF model parameters in a FITS table + assert (psf->psfTrendMode != PM_TREND_NONE); { // we need to write a header for the table, @@ -267,9 +270,13 @@ } + char *modeName = pmTrend2DModeToString(trend->mode); + // 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); + psMetadataAddStr (header, PS_LIST_TAIL, "MODE", 0, "", modeName); + psFree (modeName); // XXX can we drop this now? @@ -279,19 +286,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[ix][iy]); + 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); + } } } @@ -535,7 +563,14 @@ int nYorder = 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; } + 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; + } + // XXX allocate pmTrend2d based on the values here psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); } @@ -549,4 +584,7 @@ // XXX can we drop this now? psf->skyBias = psMetadataLookupF32 (&status, header, "SKY_BIAS"); + + char *modeName = psMetadataLookupStr (&status, header, "MODE"); + trend->mode = pmTrend2DModeToString(modeName); // read the raw table data @@ -561,13 +599,21 @@ // 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; + map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); + map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); + // poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); + } else { + psPolynomial2D *poly = trend->poly; + 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); Index: /branches/eam_branch_20070921/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmPSFtry.c (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmPSFtry.c (revision 14980) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.45.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 21:45:00 $ + * @version $Revision: 1.45.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * * 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,15 +56,8 @@ // allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, pmPSFOptions *options) +pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options) { - - // 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 *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); + psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); test->psf = pmPSFAlloc (options); @@ -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); } @@ -105,7 +97,12 @@ int Npsf = 0; - // XXX use the options structure here instead? - - pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, options); + // 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"); @@ -113,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++) { @@ -125,5 +122,5 @@ } - // set object mask to define valid pixels + // set object mask to define valid pixels -- XXX not unset? psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); @@ -142,4 +139,5 @@ // stage 2: construct a psf (pmPSF) from this collection of model fits + // XXX take 'applyWeights' from the psf options? if (!pmPSFFromPSFtry (psfTry, options->applyWeights)) { psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); @@ -166,5 +164,5 @@ source->modelPSF->radiusFit = options->radius; - // set object mask to define valid pixels + // set object mask to define valid pixels -- XXX not unset? psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); @@ -380,6 +378,4 @@ // 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 // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for @@ -391,4 +387,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); @@ -408,20 +406,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); @@ -458,5 +463,5 @@ // 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)) { + 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); @@ -487,6 +492,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]); } Index: /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.c =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.c (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.c (revision 14980) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 18:55:12 $ + * @version $Revision: 1.2.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * * Copyright 2004 Institute for Astronomy, University of Hawaii @@ -179,2 +179,43 @@ return result; } + +psString pmTrend2DModeToString (pmTrend2DMode mode) { + + switch (mode) { + case PM_TREND_NONE: + psString name = psStringCopy ("NONE"); + break; + case PM_TREND_POLY_ORD: + psString name = psStringCopy ("POLY_ORD"); + break; + case PM_TREND_POLY_CHEB: + psString name = psStringCopy ("POLY_CHEB"); + break; + case PM_TREND_MAP: + psString 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: /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.h =================================================================== --- /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.h (revision 14979) +++ /branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.h (revision 14980) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 21:45:00 $ + * @version $Revision: 1.2.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-22 02:20:11 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -42,4 +42,7 @@ psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y); +psString pmTrend2DModeToString (pmTrend2DMode mode); +pmTrend2DMode pmTrend2DModeFromString (psString name); + /// @} # endif Index: /branches/eam_branch_20070921/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branch_20070921/psphot/src/psphotChoosePSF.c (revision 14979) +++ /branches/eam_branch_20070921/psphot/src/psphotChoosePSF.c (revision 14980) @@ -35,7 +35,12 @@ // how to model the PSF variations across the field - char *psfTrendMode = psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE"); - assert (status); - options->psfTrendMode = pmTrend2DModeFromString (psfTrendMode); + 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 @@ -44,11 +49,8 @@ assert (status); - 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->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + options->image = psMemIncrRefCounter (readout->image); + + pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotom); psArray *stars = psArrayAllocEmpty (sources->n); @@ -71,5 +73,4 @@ psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates"); psFree (stars); - psFree (psfTrendMask); return NULL; } @@ -96,5 +97,4 @@ 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); } @@ -103,5 +103,4 @@ psFree (list); psFree (stars); - psFree (psfTrendMask); // select the best of the models @@ -134,11 +133,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 } } @@ -289,7 +282,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);