Index: trunk/psphot/doc/psfmodel.txt =================================================================== --- trunk/psphot/doc/psfmodel.txt (revision 14967) +++ 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 14967) +++ 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 14967) +++ 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 14967) +++ 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 14967) +++ 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 14967) +++ 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 14967) +++ 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 14967) +++ 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; }