Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 9770) +++ /trunk/psphot/src/Makefile.am (revision 9771) @@ -3,11 +3,14 @@ libpsphot_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) -bin_PROGRAMS = psphot +bin_PROGRAMS = psphot polyfitTest + psphot_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) psphot_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) psphot_LDADD = libpsphot.la +psphot_SOURCES = psphot.c -psphot_SOURCES = \ - psphot.c +polyfitTest_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) +polyfitTest_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) +polyfitTest_SOURCES = polyfitTest.c libpsphot_la_SOURCES = \ @@ -44,4 +47,5 @@ psphotWeightBias.c \ psphotSourceFreePixels.c \ + psphotTestPSF.c \ psphotCleanup.c Index: /trunk/psphot/src/models/pmModel_TEST1.c =================================================================== --- /trunk/psphot/src/models/pmModel_TEST1.c (revision 9771) +++ /trunk/psphot/src/models/pmModel_TEST1.c (revision 9771) @@ -0,0 +1,236 @@ +/****************************************************************************** + * this file defines the TEST1 source shape model. Note that these model functions are loaded + * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own. The + * models use a psVector to represent the set of parameters, with the sequence used to specify + * the meaning of the parameter. The meaning of the parameters may thus vary depending on the + * specifics of the model. All models which are used a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) * SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) * SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_TEST1 +# define PM_MODEL_FLUX pmModelFlux_TEST1 +# define PM_MODEL_GUESS pmModelGuess_TEST1 +# define PM_MODEL_LIMITS pmModelLimits_TEST1 +# define PM_MODEL_RADIUS pmModelRadius_TEST1 +# define PM_MODEL_FROM_PSF pmModelFromPSF_TEST1 +# define PM_MODEL_FIT_STATUS pmModelFitStatus_TEST1 + +// XXX consider changing to form PAR[SXY]*(1/PAR[SXX]^2 + 1/PAR[SYY]^2)*X*Y +// this would provide natural limits on PAR[SXY] of -0.25 : +0.25 + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + // XXX this is fitting sigma_x/sqrt(2), sigma_y/sqrt(2) + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; + psF32 t = 1 + z + z*z/2.0; + psF32 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + psF32 f = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + psF32 q = PAR[PM_PAR_I0]*r*r*t; + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; + } + return(f); +} + +// define the parameter limits +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) +{ + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; + beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; + + params_min[0][0].data.F32[PM_PAR_SKY] = -1000; + params_min[0][0].data.F32[PM_PAR_I0] = 0; + params_min[0][0].data.F32[PM_PAR_XPOS] = -100; + params_min[0][0].data.F32[PM_PAR_YPOS] = -100; + params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; + + params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; + params_max[0][0].data.F32[PM_PAR_I0] = 1e8; + params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; + params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; // XXX if we change the fitted term as above, this is naturally +/-1 + + return (TRUE); +} + +// make an initial guess for parameters +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +{ + + pmMoments *moments = source->moments; + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = moments->x; + PAR[PM_PAR_YPOS] = moments->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = 0.0; // XXX we can get this right if we do the integral + + return(true); +} + +psF64 PM_MODEL_FLUX(const psVector *params) +{ + float norm, z; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + z += DZ; + f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psEllipseShape shape; + + if (flux <= 0) + return (1.0); + if (params->data.F32[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= params->data.F32[PM_PAR_I0]) + return (1.0); + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + + // this estimates the radius assuming f(z) is roughly exp(-z) + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + + if (isnan(radius)) psAbort ("psphot.model", "error in code: never return invalid radius"); + if (radius < 0) psAbort ("psphot.model", "error in code: never return invalid radius"); + + return (radius); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + + return(true); +} + +// XXX double-check these definitions below +// this test is invalid if the parameters are derived +// from the PSF model +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); + dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + if (status) + return true; + return false; +} + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 9770) +++ /trunk/psphot/src/psphot.h (revision 9771) @@ -15,4 +15,8 @@ // top-level psphot functions const char *psphotCVSName(void); + +// XXX test functions +bool psphotTestPSF (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool pmPSFtestModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask); pmConfig *psphotArguments (int argc, char **argv); Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 9770) +++ /trunk/psphot/src/psphotApResid.c (revision 9771) @@ -1,3 +1,5 @@ # include "psphot.h" +// XXXX this code fails if there are too few sources to measure the aperture residual +// the larger problem is that the rules for accepting more polynomial terms are weak. static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_SKYBIAS; @@ -12,4 +14,10 @@ pmModel *model; pmSource *source; + + PS_ASSERT_PTR_NON_NULL(psf, false); + PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false); + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); psTimerStart ("psphot"); @@ -115,170 +123,151 @@ break; case PM_PSF_APTREND_CONSTANT: - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); - return false; - } - break; + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction"); + return false; + } + break; case PM_PSF_APTREND_SKYBIAS: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + break; case PM_PSF_APTREND_SKYSAT: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat"); + return false; + } + break; case PM_PSF_APTREND_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN"); + return false; + } + break; case PM_PSF_APTREND_XY_QUAD: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD"); + return false; + } + break; case PM_PSF_APTREND_SKY_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin"); + return false; + } + break; case PM_PSF_APTREND_SKYSAT_XY_LIN: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); - return false; - } - // apply the fit - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias"); + return false; + } + // apply the fit + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin"); + return false; + } + break; case PM_PSF_APTREND_ALL: - // first clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); - psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); - return false; - } - // fit just SkyBias and clip out objects which are too far from the median - stats->clipIter = 2; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); - return false; - } - // finally, fit x, y, SkyBias and clip out objects which are too far from the median - stats->clipIter = 3; - pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); - psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); - if (psf->ApTrend == NULL) { - psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); - return false; - } - break; + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT); + if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend"); + return false; + } + // fit just SkyBias and clip out objects which are too far from the median + stats->clipIter = 2; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias"); + return false; + } + // finally, fit x, y, SkyBias and clip out objects which are too far from the median + stats->clipIter = 3; + pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL); + if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) { + psError(PSPHOT_ERR_PHOTOM, false, "fitting all"); + return false; + } + break; default: psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 9770) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 9771) @@ -163,8 +163,8 @@ // make a model with unit central intensity at the image center - modelEXT->params->data.F32[0] = 0; - modelEXT->params->data.F32[1] = 1; - modelEXT->params->data.F32[2] = 0.5*image->numCols; - modelEXT->params->data.F32[3] = 0.5*image->numRows; + modelEXT->params->data.F32[PM_PAR_SKY] = 0; + modelEXT->params->data.F32[PM_PAR_I0] = 1; + modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*image->numCols; + modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*image->numRows; // construct a PSF model at this coordinate @@ -178,7 +178,8 @@ psF64 FWHM_X = 2*modelRadiusFunc (modelPSF->params, 0.5); - shape.sx = modelPSF->params->data.F32[4]; - shape.sy = modelPSF->params->data.F32[5]; - shape.sxy = modelPSF->params->data.F32[6]; + // XXX make sure this is consistent with the re-definition of PM_PAR_SXX + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; axes = psEllipseShapeToAxes (shape); @@ -195,4 +196,5 @@ } +// XXX fix this wrt re-definition of PM_PAR_SXX bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) { Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 9770) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 9771) @@ -126,4 +126,10 @@ pmModelAdd (flux, mask, model, false, false); + // XXX test output line + if ((fabs(x - 298) < 2) && (fabs(y - 583) < 2)) { + psphotSaveImage (NULL, inSource->pixels, "weird.obj.fits"); + psphotSaveImage (NULL, flux, "weird.mod.fits"); + } + // calculate nDOF (nPix - 1) // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); Index: /trunk/psphot/src/psphotEvalPSF.c =================================================================== --- /trunk/psphot/src/psphotEvalPSF.c (revision 9770) +++ /trunk/psphot/src/psphotEvalPSF.c (revision 9771) @@ -6,15 +6,12 @@ // identify objects consistent with PSF shape/magnitude distribution -// we expect dparams[4],dparams[5] to have a scatter of: -// sigma_x / (S/N) * sqrt(2) -// 1 / (params[4],params[5])*(S/N) - -// sigma_x : 1 / SX -// dsx : 1 / (SX * SN) -// dsx_o = hypot (1/(SX*SN), MIN_DSX) +// we expect dparams[PM_PAR_SXX],dparams[PM_PAR_SXX] to have a scatter of: +// sigma_x / (S/N) // any objects which is consistent with the PSF should have -// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) +// abs(dparams[PM_PAR_SXX]) < N * dsxLine(mag) & abs(dparams[PM_PAR_SYY]) < N * dsyLine(mag) // this includes a minimum buffer (DS) for the brighter objects + +// EAM : 2006.10.20 : I have re-defined params[PM_PAR_SXX] : it is now sqrt(2)*sigma_x, not 1/sigma_x // saturated stars should fall outside (larger), but have peaks above SATURATION @@ -51,4 +48,105 @@ float dSX, dSY, SX, SY, SN; float nSx, nSy, Chi; + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_MODE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_TYPE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_MODE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[PM_PAR_I0] >= SATURATION) { + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_MODE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[PM_PAR_I0] <= 0) { + source->mode |= PM_SOURCE_MODE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_MODE_SATSTAR; + } + + SN = model->params->data.F32[PM_PAR_I0]/model->dparams->data.F32[1]; + SX = model->params->data.F32[PM_PAR_SXX]/M_SQRT2; + SY = model->params->data.F32[PM_PAR_SYY]/M_SQRT2; + dSX = model->dparams->data.F32[PM_PAR_SXX]; + dSY = model->dparams->data.F32[PM_PAR_SYY]; + Chi = model->chisqNorm / model->nDOF; + + // swing of sigma_x,y in sigmas + nSx = dSX / hypot (MIN_DS, (SX / SN)); + nSy = dSY / hypot (MIN_DS, (SY / SN)); + + // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[] + keep = TRUE; + keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); + keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); + keep &= (SN > PSF_MIN_SN); + keep &= (Chi < PSF_MAX_CHI); + if (keep) return true; + + // this source is not a star, warn if it was a PSFSTAR + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality (%f, %f : %f %f %f %f)\n", + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], nSx, nSy, SN, Chi); + } else { + psLogMsg ("psphot", 5, "fails PSF fit (%f, %f : %f %f %f %f)\n", + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], nSx, nSy, SN, Chi); + } + + // object appears to be small, suspected defect + if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_TYPE_DEFECT; + return false; + } + + // object appears to be large, suspected extended source + if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + return false; + } + + // poor-quality fit; only keep if nothing else works... + source->mode |= PM_SOURCE_MODE_POOR; + return false; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalDBL (pmSource *source, pmModel *model) { // do we actually have a valid PSF model? @@ -103,105 +201,4 @@ source->mode &= ~PM_SOURCE_MODE_SATSTAR; } - - SN = model->params->data.F32[1]/model->dparams->data.F32[1]; - SX = model->params->data.F32[4]; - SY = model->params->data.F32[5]; - dSX = model->dparams->data.F32[4]; - dSY = model->dparams->data.F32[5]; - Chi = model->chisqNorm / model->nDOF; - - // swing of sigma_x,y in sigmas - nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); - nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); - - // assign PM_SOURCE_MODE_GOODSTAR to bright objects within PSF region of dparams[] - keep = TRUE; - keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); - keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); - keep &= (SN > PSF_MIN_SN); - keep &= (Chi < PSF_MAX_CHI); - if (keep) return true; - - // this source is not a star, warn if it was a PSFSTAR - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality (%f, %f : %f %f %f %f)\n", - model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi); - } else { - psLogMsg ("psphot", 5, "fails PSF fit (%f, %f : %f %f %f %f)\n", - model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi); - } - - // object appears to be small, suspected defect - if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) { - source->type = PM_SOURCE_TYPE_DEFECT; - return false; - } - - // object appears to be large, suspected extended source - if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { - source->type = PM_SOURCE_TYPE_EXTENDED; - return false; - } - - // poor-quality fit; only keep if nothing else works... - source->mode |= PM_SOURCE_MODE_POOR; - return false; -} - -// examine the model->status, fit parameters, etc and decide if the model succeeded -// set the source->type and source->mode appropriately -bool psphotEvalDBL (pmSource *source, pmModel *model) { - - // do we actually have a valid PSF model? - if (model == NULL) { - source->mode &= ~PM_SOURCE_MODE_FITTED; - return false; - } - - // did the model fit fail for one or another reason? - switch (model->status) { - case PM_MODEL_SUCCESS: - break; - case PM_MODEL_UNTRIED: - source->mode &= ~PM_SOURCE_MODE_FITTED; - return false; - case PM_MODEL_BADARGS: - case PM_MODEL_NONCONVERGE: - case PM_MODEL_OFFIMAGE: - default: - source->mode |= PM_SOURCE_MODE_FAIL; - return false; - } - - // unless we prove otherwise, this object is a star. - source->type = PM_SOURCE_TYPE_STAR; - - // the following source->mode information pertains to modelPSF: - source->mode |= PM_SOURCE_MODE_PSFMODEL; - - // if the object has fitted peak above saturation, label as SATSTAR - // this is a valid PSF object, but ignore the other quality tests - // remember: fit does not use saturated pixels (masked) - // XXX no extended object can saturate and stay extended... - if (model->params->data.F32[1] >= SATURATION) { - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); - } - source->mode |= PM_SOURCE_MODE_SATSTAR; - return true; - } - - // if the object has a fitted peak below 0, the fit did not converge cleanly - if (model->params->data.F32[1] <= 0) { - source->mode |= PM_SOURCE_MODE_FAIL; - return false; - } - - // if the source was predicted to be a SATSTAR, but it fitted below saturation, - // make a note to the user - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); - source->mode &= ~PM_SOURCE_MODE_SATSTAR; - } return true; } Index: /trunk/psphot/src/psphotModelGroupInit.c =================================================================== --- /trunk/psphot/src/psphotModelGroupInit.c (revision 9770) +++ /trunk/psphot/src/psphotModelGroupInit.c (revision 9771) @@ -4,8 +4,10 @@ // psModule/src/objects/models -# include "models/pmModel_TAUSS.c" +// # include "models/pmModel_TAUSS.c" +# include "models/pmModel_TEST1.c" static pmModelGroup userModels[] = { - {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, + {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1, pmModelFlux_TEST1, pmModelRadius_TEST1, pmModelLimits_TEST1, pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelFitStatus_TEST1}, +// {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, }; Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 9770) +++ /trunk/psphot/src/psphotReadout.c (revision 9771) @@ -62,4 +62,7 @@ return true; } + + // psphotTestPSF (readout, sources, recipe); + // exit (1); // use bright stellar objects to measure PSF Index: /trunk/psphot/src/psphotTestPSF.c =================================================================== --- /trunk/psphot/src/psphotTestPSF.c (revision 9771) +++ /trunk/psphot/src/psphotTestPSF.c (revision 9771) @@ -0,0 +1,132 @@ +# include "psphot.h" + +bool psphotTestPSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + char *modelName; + pmPSF *psf = NULL; + psArray *stars = NULL; + + psTimerStart ("psphot"); + + psphotSaveImage (NULL, readout->image, "image.fits"); + + // check if a PSF model is supplied by the user + psf = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.PSF"); + if (psf != NULL) return psf; + + // examine PSF sources in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // array to store candidate PSF stars + int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS"); + if (!status) { + NSTARS = PS_MIN (sources->n, 200); + psWarning("PSF_MAX_NSTARS is not set in the recipe --- defaulting to %d\n", NSTARS); + } + + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS"); + if (!status) { + POISSON_ERRORS = true; + psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n"); + } + + // how to model the PSF variations across the field + // XXX make a default value? or not? + psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK"); + psPolynomial2D *psfTrendMask; + if (!status || !md) { + psWarning("PSF.TREND.MASK is not set in the recipe --- defaulting to use zeroth order.\n"); + psfTrendMask = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0); + } else { + psfTrendMask = psPolynomial2DfromMetadata (md); + if (!psfTrendMask) { + psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); + return false; + } + } + + pmSourceFitModelInit (15, 0.1, POISSON_ERRORS); + + stars = psArrayAllocEmpty (sources->n); + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + // if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source); + psArrayAdd (stars, 200, source); + } + psLogMsg ("psphot.pspsf", 4, "selected candidate %ld PSF objects\n", stars->n); + + if (stars->n == 0) { + psError(PSPHOT_ERR_PSF, true, "Failed to find any PSF candidates"); + return NULL; + } + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + if (!status) { + psWarning("PSF_FIT_RADIUS is not set in the recipe --- defaulting to 20.0\n"); + RADIUS = 20.0; + } + + // for this test, require a single model + psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + if (mdi->type != PS_DATA_STRING) psAbort ("psphotChoosePSF", "choose a single PSF_MODEL"); + modelName = mdi->data.V; + + pmPSFtestModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask); + + psphotSaveImage (NULL, readout->image, "resid.fits"); + psphotSaveImage (NULL, readout->mask, "mask.fits"); + psphotSaveImage (NULL, readout->weight, "weight.fits"); + + return true; +} + +bool pmPSFtestModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask) +{ + bool status; + float x; + float y; + + pmModelType type = pmModelSetType (modelName); + pmPSF *psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); + if (psf == NULL) psAbort ("psphotTestPSF", "unknown model"); + + FILE *f = fopen ("params.dat", "w"); + + // stage 1: fit an independent model (freeModel) to all sources + psTimerStart ("fit"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmModel *model = pmSourceModelGuess (source, psf->type); + x = source->peak->x; + y = source->peak->y; + + // set temporary object mask and fit object + // fit model as EXT, not PSF + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); + status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // write fitted parameters to file + fprintf (f, "%f ", model->params->data.F32[PM_PAR_XPOS]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_YPOS]); + + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SXX]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SYY]); + fprintf (f, "%f ", model->params->data.F32[PM_PAR_SXY]); + + fprintf (f, "%f %d\n", model->chisq, model->nIter); + + // subtract model flux + pmModelSub (source->pixels, source->mask, model, false, false); + } + fclose (f); + psLogMsg ("psphot.psftest", 4, "fit ext: %f sec for %ld sources\n", psTimerMark ("fit"), sources->n); + return true; +}