Index: trunk/psModules/src/camera/pmFPARead.c =================================================================== --- trunk/psModules/src/camera/pmFPARead.c (revision 25738) +++ trunk/psModules/src/camera/pmFPARead.c (revision 25754) @@ -387,4 +387,5 @@ psFree(regionString); psFree(readout); + psFree(iter); return false; } Index: trunk/psModules/src/camera/pmHDUGenerate.c =================================================================== --- trunk/psModules/src/camera/pmHDUGenerate.c (revision 25738) +++ trunk/psModules/src/camera/pmHDUGenerate.c (revision 25754) @@ -611,4 +611,5 @@ if (cells->n == 0) { // Nothing to do + psFree (cells); return true; } @@ -660,4 +661,5 @@ if (cells->n == 0) { // Nothing to do + psFree (cells); return true; } @@ -707,8 +709,11 @@ if (cells->n == 0) { // Nothing to do + psFree (cells); return true; } - - return generateHDU(hdu, cells); + + bool status = generateHDU(hdu, cells); + psFree (cells); + return status; } default: Index: trunk/psModules/src/camera/pmReadoutFake.c =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.c (revision 25738) +++ trunk/psModules/src/camera/pmReadoutFake.c (revision 25754) @@ -97,4 +97,17 @@ continue; } + // check that all params are valid: + bool validParams = true; + for (int n = 0; validParams && (n < normModel->params->n); n++) { + if (n == PM_PAR_SKY) continue; + if (n == PM_PAR_I0) continue; + if (n == PM_PAR_XPOS) continue; + if (n == PM_PAR_YPOS) continue; + if (!isfinite(normModel->params->data.F32[n])) validParams = false; + } + if (!validParams) { + psFree(normModel); + continue; + } if (circularise && !circulariseModel(normModel)) { psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); @@ -112,4 +125,17 @@ continue; } + // check that all params are valid: + bool validParams = true; + for (int n = 0; validParams && (n < fakeModel->params->n); n++) { + if (n == PM_PAR_SKY) continue; + if (n == PM_PAR_I0) continue; + if (n == PM_PAR_XPOS) continue; + if (n == PM_PAR_YPOS) continue; + if (!isfinite(fakeModel->params->data.F32[n])) validParams = false; + } + if (!validParams) { + psFree(fakeModel); + continue; + } if (circularise && !circulariseModel(fakeModel)) { psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); Index: trunk/psModules/src/extras/pmVisual.c =================================================================== --- trunk/psModules/src/extras/pmVisual.c (revision 25738) +++ trunk/psModules/src/extras/pmVisual.c (revision 25754) @@ -22,4 +22,7 @@ #include "pmSubtractionStamps.h" #include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmSource.h" #include "pmFPAExtent.h" @@ -86,9 +89,13 @@ { char key[10]; - fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?"); + if (plotFlag) { + fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) "); + } else { + fprintf (stdout, "[c]ontinue? [a]bort all visual plots? (c) "); + } if (!fgets(key, 8, stdin)) { psWarning("Unable to read option"); } - if (key[0] == 's') { + if (plotFlag && (key[0] == 's')) { *plotFlag = false; } Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 25738) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 25754) @@ -381,5 +381,6 @@ options->poissonErrorsParams = true; options->stats = psStatsAlloc(PSF_STATS); - options->radius = maxRadius; + options->fitRadius = maxRadius; + options->apRadius = maxRadius; // XXX need to decide if aperture mags need a different radius options->psfTrendMode = PM_TREND_MAP; options->psfTrendNx = xOrder; Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 25738) +++ trunk/psModules/src/objects/Makefile.am (revision 25754) @@ -50,4 +50,9 @@ pmPSF_IO.c \ pmPSFtry.c \ + pmPSFtryModel.c \ + pmPSFtryFitEXT.c \ + pmPSFtryMakePSF.c \ + pmPSFtryFitPSF.c \ + pmPSFtryMetric.c \ pmTrend2D.c \ pmGrowthCurveGenerate.c \ Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 25754) @@ -1,8 +1,8 @@ /****************************************************************************** * this file defines the GAUSS 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 + * by pmModelClass.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 + * specifics of the model. All models which are used as a PSF representations share a few * parameters, for which # define names are listed in pmModel.h: @@ -54,4 +54,6 @@ // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords psF32 PM_MODEL_FUNC(psVector *deriv, const psVector *params, @@ -163,4 +165,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -178,5 +181,5 @@ psEllipseShape shape = psEllipseAxesToShape (axes); - PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; PAR[PM_PAR_XPOS] = peak->xf; @@ -230,4 +233,6 @@ psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); + psAssert (isfinite(radius), "fix this code: radius should not be nan for %f", PAR[PM_PAR_I0]); + return (radius); } @@ -340,5 +345,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - psF32 dP; bool status; @@ -346,17 +350,9 @@ 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; + return status; } Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 25754) @@ -1,8 +1,8 @@ /****************************************************************************** * this file defines the PGAUSS 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 + * by pmModelClass.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 + * specifics of the model. All models which are used as a PSF representations share a few * parameters, for which # define names are listed in pmModel.h: @@ -54,4 +54,6 @@ // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords psF32 PM_MODEL_FUNC(psVector *deriv, const psVector *params, @@ -165,4 +167,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -179,5 +182,5 @@ psEllipseShape shape = psEllipseAxesToShape (axes); - PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; PAR[PM_PAR_XPOS] = peak->xf; @@ -258,5 +261,7 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / 3.0)); + psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_I0]); float z1 = (1.0 / limit); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_I0]); z1 = PS_MAX (z0, z1); z0 = 0.0; @@ -389,5 +394,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - psF32 dP; bool status; @@ -395,11 +399,5 @@ 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); Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 25754) @@ -1,9 +1,9 @@ /****************************************************************************** - * this file defines the PS1_V1 source shape model (XXX need a better name!). 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: + * this file defines the PS1_V1 source shape model. Note that these model functions are loaded + * by pmModelClass.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 as a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: power-law with fitted linear term @@ -42,4 +42,8 @@ # define ALPHA_M 0.666 +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords + // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; @@ -57,5 +61,4 @@ static bool limitsApply = true; // Apply limits? - psF32 PM_MODEL_FUNC (psVector *deriv, @@ -182,4 +185,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -206,5 +210,4 @@ if (!isfinite(shape.sxy)) return false; - // XXX turn this off here for now PAR[PM_PAR_SKY] = moments->Sky; PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; @@ -292,8 +295,8 @@ // choose a z value guaranteed to be beyond our limit - float z0 = pow((1.0 / limit), (1.0 / ALPHA)); - float z1 = (1.0 / limit) / PAR[PM_PAR_7]; - z1 = PS_MAX (z0, z1); - z0 = 0.0; + float z0 = 0.0; + float z1 = pow((1.0 / limit), (1.0 / ALPHA)); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); + if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0; // perform a type of bisection to find the value @@ -420,6 +423,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - - psF32 dP; bool status; @@ -427,11 +428,5 @@ 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); Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 25754) @@ -1,5 +1,5 @@ /****************************************************************************** * this file defines the QGAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * model functions are loaded by pmModelClass.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 @@ -38,4 +38,8 @@ # define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS # define PM_MODEL_SET_LIMITS pmModelSetLimits_QGAUSS + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords // Lax parameter limits @@ -178,4 +182,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -202,5 +207,4 @@ if (!isfinite(shape.sxy)) return false; - // XXX turn this off here for now PAR[PM_PAR_SKY] = moments->Sky; PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; @@ -283,5 +287,7 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / 2.25)); + psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]); float z1 = (1.0 / limit) / PAR[PM_PAR_7]; + psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); z1 = PS_MAX (z0, z1); z0 = 0.0; @@ -409,6 +415,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - - psF32 dP; bool status; @@ -416,11 +420,5 @@ 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); Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 25754) @@ -1,8 +1,8 @@ /****************************************************************************** * this file defines the RGAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' + * model functions are loaded by pmModelClass.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 + * may thus vary depending on the specifics of the model. All models which are used as a PSF * representations share a few parameters, for which # define names are listed in pmModel.h: @@ -39,4 +39,8 @@ # define PM_MODEL_SET_LIMITS pmModelSetLimits_RGAUSS +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords + // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 }; @@ -87,5 +91,5 @@ dPAR[PM_PAR_SXY] = -q*X*Y; - // this model derivative is undefined at z = 0.0, but is actually 0.0 + // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0 dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -5.0*t*log(z)*p*z; } @@ -172,4 +176,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -196,5 +201,5 @@ if (!isfinite(shape.sxy)) return false; - PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; PAR[PM_PAR_XPOS] = peak->xf; @@ -276,5 +281,7 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7])); + psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]); float z1 = (1.0 / limit); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); z1 = PS_MAX (z0, z1); z0 = 0.0; @@ -402,6 +409,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - - psF32 dP; bool status; @@ -409,11 +414,5 @@ 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); Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 25754) @@ -1,8 +1,8 @@ /****************************************************************************** * this file defines the SERSIC 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 + * by pmModelClass.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 + * specifics of the model. All models which are used as a PSF representations share a few * parameters, for which # define names are listed in pmModel.h: @@ -41,4 +41,8 @@ # define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC # define PM_MODEL_SET_LIMITS pmModelSetLimits_SERSIC + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords // Lax parameter limits @@ -186,4 +190,5 @@ // make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { @@ -210,5 +215,4 @@ if (!isfinite(shape.sxy)) return false; - // XXX PAR[PM_PAR_SKY] = moments->Sky; PAR[PM_PAR_SKY] = 0.0; PAR[PM_PAR_I0] = peak->flux; @@ -284,6 +288,8 @@ psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7])); + psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); psF64 radius = sigma * sqrt (2.0 * z); + psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma); if (isnan(radius)) @@ -392,6 +398,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - - psF32 dP; bool status; @@ -399,17 +403,14 @@ 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); +<<<<<<< .working fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n", dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); +======= +>>>>>>> .merge-right.r25750 return status; } Index: trunk/psModules/src/objects/models/pmModel_SGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SGAUSS.c (revision 25738) +++ (revision ) @@ -1,389 +1,0 @@ -/****************************************************************************** - * this file defines the SGAUSS source shape model (XXX need a better name!). 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: - - power-law with fitted slope and outer tidal radius - 1 / (1 + z^N + kz^4) - - * 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 - * PM_PAR_7 7 - slope of power-law component (N) - * PM_PAR_8 8 - amplitude of the tidal cutoff (k) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - the code for guessing the power-law slope based on the radial profile - is either too slow or does not work well. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_SGAUSS -# define PM_MODEL_FLUX pmModelFlux_SGAUSS -# define PM_MODEL_GUESS pmModelGuess_SGAUSS -# define PM_MODEL_LIMITS pmModelLimits_SGAUSS -# define PM_MODEL_RADIUS pmModelRadius_SGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_SGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_SGAUSS - -psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image); - -psF32 PM_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8); - // note that if z -> 0, dPAR[7] -> -inf - // also z^(PAR[7]-1) -> Inf - - psF32 pr = z*PAR[8]; - psF32 pr3 = pr*pr*pr; - psF32 p = pow(z, PAR[7] - 1.0); - psF32 r = 1.0 / (1 + z*p + pr*pr3); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3); - - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -2*t*log(z)*z*p; - deriv->data.F32[8] = -2*t*4*z*pr3; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (9, PS_TYPE_F32); - *params_min = psVectorAlloc (9, PS_TYPE_F32); - *params_max = psVectorAlloc (9, PS_TYPE_F32); - - beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 10000; - beta_lim[0][0].data.F32[2] = 5; - beta_lim[0][0].data.F32[3] = 5; - beta_lim[0][0].data.F32[4] = 0.5; - beta_lim[0][0].data.F32[5] = 0.5; - beta_lim[0][0].data.F32[6] = 0.5; - beta_lim[0][0].data.F32[7] = 0.5; - beta_lim[0][0].data.F32[8] = 0.05; - - params_min[0][0].data.F32[0] = -1000; - params_min[0][0].data.F32[1] = 0; - params_min[0][0].data.F32[2] = -100; - params_min[0][0].data.F32[3] = -100; - params_min[0][0].data.F32[4] = 0.01; - params_min[0][0].data.F32[5] = 0.01; - params_min[0][0].data.F32[6] = -5.0; - params_min[0][0].data.F32[7] = 0.5; - params_min[0][0].data.F32[8] = 0.001; - - params_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e6; - params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[4] = 2.0; - params_max[0][0].data.F32[5] = 2.0; - params_max[0][0].data.F32[6] = +3.0; - params_max[0][0].data.F32[7] = 5.0; - params_max[0][0].data.F32[8] = 0.5; - - return (TRUE); -} - -bool PM_MODEL_GUESS (pmModel *model, pmSource *source) -{ - - pmMoments *sMoments = source->moments; - pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(sMoments->Sx); - moments.y2 = PS_SQR(sMoments->Sy); - moments.xy = sMoments->Sxy; - - // solve the math to go from Moments To Shape - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params[0] = sMoments->Sky; - params[1] = sMoments->Peak - sMoments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = 1.0 / shape.sx; - params[5] = 1.0 / shape.sy; - params[6] = shape.sxy; - - # if (0) - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - axes.major *= 2.0; - axes.minor *= 2.0; - f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - - if (f1 > f2) { - params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); - } else { - params[7] = 0.5; - } - # endif - - params[7] = 1.8; - params[8] = 0.1; - - - return(true); -} - -psF64 PM_MODEL_FLUX (const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - for (z = 0.005; z < 50; z += 0.01) { - psF32 pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// 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) -{ - psF64 r, z = 0.0, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - return (1.0); - - // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // XXX : we can do this faster with an intelligent starting choice - for (r = 0.0; r < 20.0; r += dr) { - z = PS_SQR(r); - pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 9; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - // XXX: Verify this (from EAM change) - //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]); - out[i] = psPolynomial2DEval(poly, out[2], out[3]); - } - return(true); -} - -bool PM_MODEL_FIT_STATUS (pmModel *model) -{ - - psF32 dP; - bool status; - psEllipseAxes axes; - psEllipseShape shape; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - - dP = 0; - dP += PS_SQR(dPAR[4] / PAR[4]); - dP += PS_SQR(dPAR[5] / PAR[5]); - dP += PS_SQR(dPAR[7] / PAR[7]); - dP = sqrt (dP); - - status = true; - status &= (dP < 0.5); - status &= (PAR[1] > 0); - status &= ((dPAR[1]/PAR[1]) < 0.5); - status &= (fabs(PAR[8]) < 0.5); - status &= (dPAR[8] < 0.1); - status &= (axes.major > 1.41); - status &= (axes.minor > 1.41); - status &= ((axes.major / axes.minor) < 5.0); - status &= (PAR[7] > 0.5); - - if (status) - return true; - return false; -} - -// measure the flux for the elliptical contour -psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image) -{ - - double t, dt, ct, st, xo, yo, value; - int N, Nt, x, y; - - // choose dt to uniformly divide contour, with ~1 pix spacing at most - dt = asin (1 / axes.minor); - Nt = (int)(2*M_PI / dt) + 1; - dt = 2*M_PI / Nt; - - ct = cos(axes.theta); - st = sin(axes.theta); - xo = xc - image->col0; - yo = yc - image->row0; - - psVector *contour = psVectorAlloc (Nt, PS_TYPE_F32); - contour->n = contour->nalloc; - for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) { - x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo; - y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo; - value = p_psImageGetElementF64(image, x, y); - if (isfinite(value)) { - contour->data.F32[N] = value; - N++; - } - } - contour->n = N; - // accept every pixel: double counting is not so problematic here... - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - if (!psVectorStats (stats, contour, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - value = stats->sampleMedian; - - psFree (stats); - psFree (contour); - - return (value); -} - -// XXX EAM : test version using flux contours to guess slope -bool PM_MODEL_GUESS_HARD (pmModel *model, pmSource *source) -{ - - pmMoments *sMoments = source->moments; - pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - float f1, f2; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(sMoments->Sx); - moments.y2 = PS_SQR(sMoments->Sy); - moments.xy = sMoments->Sxy; - - // solve the math to go from Moments To Shape - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params[0] = sMoments->Sky; - params[1] = sMoments->Peak - sMoments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = 1.0 / shape.sx; - params[5] = 1.0 / shape.sy; - params[6] = shape.sxy; - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - axes.major *= 2.0; - axes.minor *= 2.0; - f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - - if (f1 > f2) { - params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); - } else { - params[7] = 0.5; - } - params[8] = 0.1; - - return(true); -} - -# 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/psModules/src/objects/models/pmModel_TGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_TGAUSS.c (revision 25738) +++ (revision ) @@ -1,202 +1,0 @@ -/****************************************************************************** - * this file defines the TGAUSS source shape model (XXX need a better name!). 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: - - power-law with fixed slope and fitted amplitude - 1 / (1 + z + kz^2.2) - - * 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 - * PM_PAR_7 7 - amplitude of high-order component (k) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_TGAUSS -# define PM_MODEL_FLUX pmModelFlux_TGAUSS -# define PM_MODEL_GUESS pmModelGuess_TGAUSS -# define PM_MODEL_LIMITS pmModelLimits_TGAUSS -# define PM_MODEL_RADIUS pmModelRadius_TGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_TGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_TGAUSS - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; - - psF32 p = pow(z, 1.2); - psF32 r = 1.0 / (1 + z + PAR[7]*z*p); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(1 + PAR[7]*2.2*p); - - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -t*z*p; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, 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; - beta_lim[0][0].data.F32[PM_PAR_7] = 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_min[0][0].data.F32[PM_PAR_7] = 0.1; - - 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; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = 1.0/source->moments->Sx; - params->data.F32[5] = 1.0/source->moments->Sy; - // params->data.F32[6] = source->moments->Sxy; - params->data.F32[6] = 0.0; - params->data.F32[7] = 5.0; - return(true); -} - -psF64 PS_MODEL_FLUX (const psVector *params) -{ - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - - psF64 Flux = params->data.F32[1] * Area; - - return(Flux); -} - -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1]) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 8; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -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); - - return status; -} - -# 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/psModules/src/objects/models/pmModel_WAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_WAUSS.c (revision 25738) +++ (revision ) @@ -1,198 +1,0 @@ -/****************************************************************************** - * this file defines the WAUSS source shape model (XXX need a better name!). 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: - - power-law with fitted linear term - 1 / (1 + Az + Bz^2 + z^3/6) - - * 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 (SigmaX / sqrt(2)) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) - * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - amplitude of the linear component (A) - * PM_PAR_8 8 - amplitude of the quadratic component (B) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_WAUSS -# define PM_MODEL_FLUX pmModelFlux_WAUSS -# define PM_MODEL_GUESS pmModelGuess_WAUSS -# define PM_MODEL_LIMITS pmModelLimits_WAUSS -# define PM_MODEL_RADIUS pmModelRadius_WAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_WAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[2]; - psF32 px = params->data.F32[4]*X; - psF32 py = params->data.F32[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; - psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); - psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - if (deriv != NULL) { - // note difference from gaussian: q = params->data.F32[1]*r - psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; - deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; - // The values of 100 dampen the swing of params->data.F32[7,8] */ - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, 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; - beta_lim[0][0].data.F32[PM_PAR_7] = 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_min[0][0].data.F32[PM_PAR_7] = 0.1; - - 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; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - // XXX: What are these? - // params->data.F32[7] = B2; - // params->data.F32[8] = B3; - return(true); -} - -// this is probably wrong since it uses the gauss integral 2 pi sigma^2 -psF64 PS_MODEL_FLUX (const psVector *params) -{ - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - - psF64 Flux = params->data.F32[1] * Area; - - return(Flux); -} - -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1] <= 0) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 9; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -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); - - return status; -} - -# 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/psModules/src/objects/models/pmModel_ZGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 25738) +++ (revision ) @@ -1,251 +1,0 @@ -/****************************************************************************** - * this file defines the ZGAUSS source shape model (XXX need a better name!). 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: - - power-law with fitted slope and tidal cutoff - 1 / (1 + z^N + (Az)^4) - - * 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 - * PM_PAR_7 7 - slope of power-law component (N) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_ZGAUSS -# define PM_MODEL_FLUX pmModelFlux_ZGAUSS -# define PM_MODEL_GUESS pmModelGuess_ZGAUSS -# define PM_MODEL_LIMITS pmModelLimits_ZGAUSS -# define PM_MODEL_RADIUS pmModelRadius_ZGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_ZGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_ZGAUSS - -# define PAR8 0.1 - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; - - psF32 pr = PAR8*z; - psF32 p = pow(z, PAR[7] - 1.0); - psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR8*pr*pr*pr); - - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - deriv->data.F32[4] = -q*px*X; - deriv->data.F32[5] = -q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -t*log(z)*z*p; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, 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; - beta_lim[0][0].data.F32[PM_PAR_7] = 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_min[0][0].data.F32[PM_PAR_7] = 0.1; - - 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; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(source->moments->Sx); - moments.y2 = PS_SQR(source->moments->Sy); - moments.xy = source->moments->Sxy; - - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = 1.0 / shape.sx; - params->data.F32[5] = 1.0 / shape.sy; - params->data.F32[6] = shape.sxy; - params->data.F32[7] = 1.9; - return(true); -} - -psF64 PS_MODEL_FLUX (const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - psF32 pr = PAR8*z; - for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - psF64 r, z, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - return (1.0); - - // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // XXX : we can do this faster with an intelligent starting choice - for (r = 0.0; r < 20.0; r += dr) { - z = SQ(r); - pr = PAR8*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 8; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -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); - - return status; -} - -# 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/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 25738) +++ trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 25754) @@ -65,4 +65,5 @@ // use the center of the center pixel of the image + // 0.5 PIX: is this offset needed? probably -- the psf model uses 0.5,0.5 as the center, double check float xc = (int)(ix*readout->image->numCols + 0.5*readout->image->numCols) + readout->image->col0 + 0.5; float yc = (int)(iy*readout->image->numRows + 0.5*readout->image->numRows) + readout->image->row0 + 0.5; @@ -195,5 +196,5 @@ return NULL; } - psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageMaskPixels (mask, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask // the 'ignore' mode is for testing Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 25738) +++ trunk/psModules/src/objects/pmModel.c (revision 25754) @@ -56,10 +56,11 @@ tmp->type = type; - tmp->chisq = 0.0; - tmp->chisqNorm = 0.0; + tmp->mag = NAN; + tmp->chisq = NAN; + tmp->chisqNorm = NAN; tmp->nDOF = 0; tmp->nPix = 0; tmp->nIter = 0; - tmp->radiusFit = 0; + tmp->fitRadius = 0; tmp->flags = PM_MODEL_STATUS_NONE; tmp->residuals = NULL; // XXX should the model own this memory? @@ -109,5 +110,5 @@ new->nIter = model->nIter; new->flags = model->flags; - new->radiusFit = model->radiusFit; + new->fitRadius = model->fitRadius; for (int i = 0; i < new->params->n; i++) { @@ -190,7 +191,7 @@ psVector *params = model->params; - psS32 imageCol; - psS32 imageRow; - psF32 pixelValue; + float imageCol; + float imageRow; + float pixelValue; // save original values; restore before returning @@ -233,17 +234,17 @@ psF32 **Rx = NULL; psF32 **Ry = NULL; - psImageMaskType **Rm = NULL; + pmResidMaskType **Rm = NULL; if (model->residuals) { - DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5; - DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5; - Ro = (model->residuals->Ro) ? model->residuals->Ro->data.F32 : NULL; - Rx = (model->residuals->Rx) ? model->residuals->Rx->data.F32 : NULL; - Ry = (model->residuals->Ry) ? model->residuals->Ry->data.F32 : NULL; - Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL; - if (Ro) { - NX = model->residuals->Ro->numCols; - NY = model->residuals->Ro->numRows; - } + DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5; + DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5; + Ro = (model->residuals->Ro) ? model->residuals->Ro->data.F32 : NULL; + Rx = (model->residuals->Rx) ? model->residuals->Rx->data.F32 : NULL; + Ry = (model->residuals->Ry) ? model->residuals->Ry->data.F32 : NULL; + Rm = (model->residuals->mask) ? model->residuals->mask->data.PM_TYPE_RESID_MASK_DATA : NULL; + if (Ro) { + NX = model->residuals->Ro->numCols; + NY = model->residuals->Ro->numRows; + } } @@ -256,11 +257,11 @@ continue; - // XXX should we use using 0.5 pixel offset? // Convert to coordinate in parent image, with offset (dx,dy) - imageCol = ix + image->col0 - dx; - imageRow = iy + image->row0 - dy; - - x->data.F32[0] = (float) imageCol; - x->data.F32[1] = (float) imageRow; + // 0.5 PIX: the model take pixel coordinates so convert the pixel index here + imageCol = ix + 0.5 + image->col0 - dx; + imageRow = iy + 0.5 + image->row0 - dy; + + x->data.F32[0] = imageCol; + x->data.F32[1] = imageRow; pixelValue = 0.0; Index: trunk/psModules/src/objects/pmModel.h =================================================================== --- trunk/psModules/src/objects/pmModel.h (revision 25738) +++ trunk/psModules/src/objects/pmModel.h (revision 25754) @@ -107,5 +107,5 @@ int nIter; ///< number of iterations to reach min pmModelStatus flags; ///< model status flags - float radiusFit; ///< fit radius actually used + float fitRadius; ///< fit radius actually used pmResiduals *residuals; ///< normalized PSF residuals Index: trunk/psModules/src/objects/pmModelClass.c =================================================================== --- trunk/psModules/src/objects/pmModelClass.c (revision 25738) +++ trunk/psModules/src/objects/pmModelClass.c (revision 25754) @@ -36,5 +36,5 @@ #include "pmErrorCodes.h" -// XXX shouldn't these be defined for us in pslib.h ??? +// XXX shouldn't these be defined for us in math.h ??? double hypot(double x, double y); double sqrt (double x); Index: trunk/psModules/src/objects/pmModelGroup.c =================================================================== --- trunk/psModules/src/objects/pmModelGroup.c (revision 25738) +++ (revision ) @@ -1,204 +1,0 @@ -/** @file pmModelGroup.c - * - * Functions to define and manipulate object model attributes - * - * @author GLG, MHPCC - * @author EAM, IfA - * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-12-08 02:51:14 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmSpan.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmModel.h" -#include "pmModelGroup.h" -#include "pmErrorCodes.h" - -// XXX shouldn't these be defined for us in pslib.h ??? -double hypot(double x, double y); -double sqrt (double x); - -# include "models/pmModel_GAUSS.c" -# include "models/pmModel_PGAUSS.c" -# include "models/pmModel_QGAUSS.c" -# include "models/pmModel_RGAUSS.c" -# include "models/pmModel_SERSIC.c" - -static pmModelGroup defaultModels[] = { - {"PS_MODEL_GAUSS", 7, pmModelFunc_GAUSS, pmModelFlux_GAUSS, pmModelRadius_GAUSS, pmModelLimits_GAUSS, pmModelGuess_GAUSS, pmModelFromPSF_GAUSS, pmModelFitStatus_GAUSS}, - {"PS_MODEL_PGAUSS", 7, pmModelFunc_PGAUSS, pmModelFlux_PGAUSS, pmModelRadius_PGAUSS, pmModelLimits_PGAUSS, pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS}, - {"PS_MODEL_QGAUSS", 8, pmModelFunc_QGAUSS, pmModelFlux_QGAUSS, pmModelRadius_QGAUSS, pmModelLimits_QGAUSS, pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS}, - {"PS_MODEL_RGAUSS", 8, pmModelFunc_RGAUSS, pmModelFlux_RGAUSS, pmModelRadius_RGAUSS, pmModelLimits_RGAUSS, pmModelGuess_RGAUSS, pmModelFromPSF_RGAUSS, pmModelFitStatus_RGAUSS}, - {"PS_MODEL_SERSIC", 8, pmModelFunc_SERSIC, pmModelFlux_SERSIC, pmModelRadius_SERSIC, pmModelLimits_SERSIC, pmModelGuess_SERSIC, pmModelFromPSF_SERSIC, pmModelFitStatus_SERSIC} - }; - -static pmModelGroup *models = NULL; -static int Nmodels = 0; - -static void ModelGroupFree (pmModelGroup *modelGroup) -{ - - if (modelGroup == NULL) - return; - return; -} - -pmModelGroup *pmModelGroupAlloc (int nModels) -{ - - pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup)); - psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree); - return (modelGroup); -} - -bool psMemCheckModelGroup(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) ModelGroupFree); -} - -void pmModelGroupAdd (pmModelGroup *model) -{ - - if (models == NULL) { - pmModelGroupInit (); - } - - Nmodels ++; - models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup)); - models[Nmodels-1] = model[0]; - return; -} - -bool pmModelGroupInit (void) -{ - - // if we do not need to init, return false; - if (models != NULL) - return false; - - int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup); - - models = pmModelGroupAlloc (Nnew); - for (int i = 0; i < Nnew; i++) { - models[i] = defaultModels[i]; - } - Nmodels = Nnew; - return true; -} - -void pmModelGroupCleanup (void) -{ - psFree (models); - models = NULL; - return; -} - -pmModelFunc pmModelFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFunc); -} - -pmModelFlux pmModelFlux_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFlux); -} - -pmModelRadius pmModelRadius_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelRadius); -} - -pmModelLimits pmModelLimits_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelLimits); -} - -pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelGuessFunc); -} - -pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFitStatusFunc); -} - -pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFromPSFFunc); -} - -psS32 pmModelParameterCount (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (0); - } - return (models[type].nParams); -} - -psS32 pmModelSetType (char *name) -{ - for (int i = 0; i < Nmodels; i++) { - if (!strcmp(models[i].name, name)) { - return (i); - } - } - return (-1); -} - -char *pmModelGetType (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].name); -} Index: trunk/psModules/src/objects/pmModelGroup.h =================================================================== --- trunk/psModules/src/objects/pmModelGroup.h (revision 25738) +++ (revision ) @@ -1,201 +1,0 @@ -/* @file pmModelGroup.h - * - * The object model function types are desined to allow for the flexible addition of new object - * models. Every object model, with parameters represented by pmModel, has an associated set of - * functions which provide necessary support operations. A set of abstract functions allow the - * programmer to select the approriate function or property for a specific named object model. - * - * Every model instance belongs to a class of models, defined by the value of - * the pmModelType type entry. Various functions need access to information about - * each of the models. Some of this information varies from model to model, and - * may depend on the current parameter values or other data quantities. In order - * to keep the code from requiring the information about each model to be coded - * into the low-level fitting routines, we define a collection of functions which - * allow us to abstract this type of model-dependent information. These generic - * functions take the model type and return the corresponding function pointer - * for the specified model. Each model is defined by creating this collection of - * specific functions, and placing them in a single file for each model. We - * define the following structure to carry the collection of information about - * the models. - * - * @author EAM, IfA - * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-10 01:09:20 $ - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -# ifndef PM_MODEL_GROUP_H -# define PM_MODEL_GROUP_H - -/// @addtogroup Objects Object Detection / Analysis Functions -/// @{ - -// This function is the model chi-square minimization function for this model. -typedef psMinimizeLMChi2Func pmModelFunc; - -// This function sets the parameter limits for this model. -typedef psMinimizeLMLimitFunc pmModelLimits; - -// This function returns the integrated flux for the given model parameters. -typedef psF64 (*pmModelFlux)(const psVector *params); - -// This function returns the radius at which the given model and parameters -// achieves the given flux. -typedef psF64 (*pmModelRadius)(const psVector *params, double flux); - -// This function provides the model guess parameters based on the details of -// the given source. -typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); - -// This function constructs the PSF model for the given source based on the -// supplied psf and the EXT model for the object. -typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf); - -// This function sets the model parameters based on the PSF for a given coordinate and central -// intensity -typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io); - -// This function returns the success / failure status of the given model fit -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - -typedef struct -{ - char *name; - int nParams; - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuessFunc; - pmModelFromPSFFunc modelFromPSFFunc; - pmModelParamsFromPSF modelParamsFromPSF; - pmModelFitStatusFunc modelFitStatusFunc; -} -pmModelGroup; - -// allocate a pmModelGroup to hold nModels entries -pmModelGroup *pmModelGroupAlloc (int nModels); - -bool psMemCheckModelGroup(psPtr ptr); - -// initialize the internal (static) model group with the default models -bool pmModelGroupInit (void); - -// free the internal (static) model group -void pmModelGroupCleanup (void); - -// add a new model to the internal (static) model group -void pmModelGroupAdd (pmModelGroup *model); - -/* This function returns the number of parameters used by the listed function. - */ -int pmModelParameterCount( - pmModelType type ///< Add comment. -); - - -/* This function returns the user-space model names for the specified model type. - */ -char *pmModelGetType( - pmModelType type ///< Add comment. -); - - -/** - * - * This function returns the internal model type code for the user-space model names. - * - */ -pmModelType pmModelSetType( - char *name ///< Add comment. -); - -/** - * - * Each of the function types above has a corresponding function which returns - * the function given the model type: - * - */ - -/** - * - * pmModelFunc is the function used to determine the value of the model at a - * specific coordinate, and is the one used by psMinimizeLMChi2. - * - */ -pmModelFunc pmModelFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFlux returns the total integrated flux for the given input parameters. - * - */ -pmModelFlux pmModelFlux_GetFunction (pmModelType type); - - -/** - * - * pmModelRadius returns the scaling radius at which the flux of the model - * matches the specified flux. This presumes that the model is a function of an - * elliptical contour. - * - */ -pmModelRadius pmModelRadius_GetFunction (pmModelType type); - - -/** - * - * pmModelLimits sets the parameter limit vectors for the function. - * - */ -pmModelLimits pmModelLimits_GetFunction (pmModelType type); - - -/** - * - * pmModelGuessFunc generates an initial guess for the model based on the - * provided source statistics (moments and pixel values as needed). - * - */ -pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFromPSFFunc takes as input a representation of the psf and a value - * for the model and fills in the PSF parameters of the model. The input - * primarily relies upon the centroid coordinates of the input model, though the - * normalization may potentially be used. - * - */ -pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFitStatusFunc returns a true or false values based on the success or - * failure of a model fit. The success is determined by quantities such as the - * chisq or the signal-to-noise. - * - */ -pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); - - -/** pmSourceModelGuess() - * - * Convert available data to an initial guess for the given model. This - * function allocates a pmModel entry for the pmSource structure based on the - * provided model selection. The method of defining the model parameter guesses - * are specified for each model below. The guess values are placed in the model - * parameters. The function returns TRUE on success or FALSE on failure. - * - */ -pmModel *pmSourceModelGuess( - pmSource *source, ///< The input pmSource - pmModelType model ///< The type of model to be created. -); - -/// @} -# endif /* PM_MODEL_GROUP_H */ Index: trunk/psModules/src/objects/pmModelUtils.c =================================================================== --- trunk/psModules/src/objects/pmModelUtils.c (revision 25738) +++ trunk/psModules/src/objects/pmModelUtils.c (revision 25754) @@ -48,5 +48,5 @@ return NULL; } - // XXX note that model->residuals is just a reference + // note that model->residuals is just a reference modelPSF->residuals = psf->residuals; @@ -65,11 +65,9 @@ // set model parameters for this source based on PSF information if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) { - // XXX we do not want to raise an error here, just note that we failed - // psError(PM_ERR_PSF, false, "Failed to set model params from PSF"); psFree(modelPSF); return NULL; } - // XXX note that model->residuals is just a reference + // note that model->residuals is just a reference modelPSF->residuals = psf->residuals; Index: trunk/psModules/src/objects/pmPSF.h =================================================================== --- trunk/psModules/src/objects/pmPSF.h (revision 25738) +++ trunk/psModules/src/objects/pmPSF.h (revision 25754) @@ -38,6 +38,5 @@ * */ -typedef struct -{ +typedef struct { pmModelType type; ///< PSF Model in use psArray *params; ///< Model parameters (psPolynomial2D) @@ -65,6 +64,5 @@ pmGrowthCurve *growth; ///< apMag vs Radius pmResiduals *residuals; ///< normalized residual image (no spatial variation) -} -pmPSF; +} pmPSF; typedef struct { @@ -81,5 +79,6 @@ bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting bool poissonErrorsParams; ///< use poission errors for model parameter fitting - float radius; + float fitRadius; + float apRadius; bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? } pmPSFOptions; Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 25738) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 25754) @@ -37,51 +37,4 @@ #include "pmSourceVisual.h" -bool printTrendMap (pmTrend2D *trend) { - - if (!trend->map) return false; - if (!trend->map->map) return false; - - for (int j = 0; j < trend->map->map->numRows; j++) { - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); - } - fprintf (stderr, "\t\t\t"); - for (int i = 0; i < trend->map->map->numCols; i++) { - fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); - } - fprintf (stderr, "\n"); - } - return true; -} - -bool psImageMapCleanup (psImageMap *map) { - - if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true; - - // find the weighted average of all pixels - float Sum = 0.0; - float Wt = 0.0; - for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (!isfinite(map->map->data.F32[j][i])) continue; - Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i]; - Wt += map->error->data.F32[j][i]; - } - } - - float Mean = Sum / Wt; - - // do any of the pixels in the map need to be repaired? - // XXX for now, we are just replacing bad pixels with the Mean - for (int j = 0; j < map->map->numRows; j++) { - for (int i = 0; i < map->map->numCols; i++) { - if (isfinite(map->map->data.F32[j][i]) && - (map->error->data.F32[j][i] > 0.0)) continue; - map->map->data.F32[j][i] = Mean; - } - } - return true; -} - // ******** pmPSFtry functions ************************************************** // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely @@ -110,5 +63,5 @@ psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); - test->psf = pmPSFAlloc (options); + test->psf = NULL; test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); @@ -136,934 +89,98 @@ } +float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) { -// build a pmPSFtry for the given model: -// - fit each source with the free-floating model -// - construct the pmPSF from the collection of models -// - fit each source with the PSF-parameter models -// - measure the pmPSF quality metric (dApResid) + psAssert(residuals, "residuals cannot be NULL"); + psAssert(errors, "errors cannot be NULL"); + psAssert(residuals->n == errors->n, "residuals and errors must be the same length"); -// sources used in for pmPSFtry may be masked by the analysis -// mask values indicate the reason the source was rejected: + // given a vector of residuals and their formal errors, calculated the necessary systematic + // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction + // highest chi-square contributors (allowed outliers) -// generate a pmPSFtry with a copy of the test PSF sources -pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) -{ - bool status; - int Next = 0; - int Npsf = 0; + psVector *mask = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK); + psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32); - // 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; + // calculate the chisq vector: + int Ngood = 0; + for (int i = 0; i < residuals->n; i++) { + chisq->data.F32[i] = PS_MAX_F32; + if (!isfinite(residuals->data.F32[i])) continue; + if (!isfinite(errors->data.F32[i])) continue; + if (errors->data.F32[i] <= 0.0) continue; + chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]); + Ngood ++; } - pmPSFtry *psfTry = pmPSFtryAlloc (sources, options); - if (psfTry == NULL) { - psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model"); - return NULL; + psVector *index = psVectorSortIndex(NULL, chisq); + + // toss out the clipFraction highest chisq values + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (i < (1.0 - clipFraction)*Ngood) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0; + } else { + mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1; + } } - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; + // Ndof ~= Ngood + // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof + // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0 + + // use Newton-Raphson to solve for S2: - // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { + // use median sigma to calculate the initial guess for S2: + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, errors, NULL, mask, 1); + float errorMedian = stats->sampleMedian; + + float nPts = 0.0; + float res2mean = 0.0; + float ChiSq = 0.0; + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; + res2mean += PS_SQR(residuals->data.F32[n]); + ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]); + nPts += 1.0; + } + res2mean /= nPts; + ChiSq /= nPts; + + float S2guess = res2mean - PS_SQR(errorMedian); - pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - if (!source->moments->nPixels) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } + psLogMsg ("psModules", 10, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", + ChiSq, residuals->n, Ngood, nPts, S2guess); - source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); - if (source->modelEXT == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); - continue; - } + for (int iter = 0; iter < 10; iter++) { - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + ChiSq = 0.0; + float dRdS = 0.0; + for (int i = 0; i < residuals->n; i++) { + int n = index->data.S32[i]; + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; + float error2 = PS_SQR(errors->data.F32[n]) + S2guess; + ChiSq += PS_SQR(residuals->data.F32[n]) / error2; + dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2); + } + ChiSq /= nPts; + dRdS /= nPts; - // fit model as EXT, not PSF - status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); + // Note the sign on dS: dRdS above is -1 * dR/dS formally + float dS = (ChiSq - 1.0) / dRdS; + S2guess += dS; + S2guess = PS_MAX(0.0, S2guess); - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // exclude the poor fits - if (!status) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); - continue; - } - Next ++; - } - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); - - if (Next == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); - psFree(psfTry); - return NULL; + psLogMsg ("psModules", 10, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess); } - // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation - if (!pmPSFFromPSFtry (psfTry)) { - psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); - psFree(psfTry); - return NULL; - } + // free local allocations + psFree (mask); + psFree (chisq); + psFree (stats); + psFree (index); - // stage 3: refit with fixed shape parameters - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { - - pmSource *source = psfTry->sources->data[i]; - - // masked for: bad model fit, outlier in parameters - if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); - continue; - } - - // set shape for this model based on PSF - source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); - if (source->modelPSF == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; - abort(); - continue; - } - source->modelPSF->radiusFit = options->radius; - - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); - - // fit the PSF model to the source - status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal); - - // skip poor fits - if (!status) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); - continue; - } - - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); - if (!status || isnan(source->apMag)) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); - continue; - } - - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - psfTry->fitMag->data.F32[i] = source->psfMag; - psfTry->metric->data.F32[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F32[i] = source->errMag; - - psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); - Npsf ++; - } - psfTry->psf->nPSFstars = Npsf; - - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); - - if (Npsf == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); - psFree(psfTry); - return NULL; - } - - // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) - // this should be linear for Poisson errors and quadratic for constant sky errors - psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK); - - // generate the x and y vectors, and mask missing models - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelPSF == NULL) { - flux->data.F32[i] = 0.0; - chisq->data.F32[i] = 0.0; - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - } else { - flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; - chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF; - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; - } - } - - // use 3hi/3lo sigma clipping on the chisq fit - psStats *stats = options->stats; - - // linear clipped fit of chisq trend vs flux - if (options->chiFluxTrend) { - bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, - mask, 0xff, chisq, NULL, flux); - psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean - psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev - - psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", - psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); - - psFree(flux); - psFree(mask); - psFree(chisq); - - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); - psFree(psfTry); - return NULL; - } - } - - for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) { - psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, - psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), - psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i)); - } - - // XXX this function wants aperture radius for pmSourcePhotometry - if (!pmPSFtryMetric (psfTry, options)) { - psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName); - psFree (psfTry); - return NULL; - } - - psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", - modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); - - return (psfTry); + return (sqrt(S2guess)); } -bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options) -{ - PS_ASSERT_PTR_NON_NULL(psfTry, false); - PS_ASSERT_PTR_NON_NULL(options, false); - PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); - - float RADIUS = options->radius; - - // the measured (aperture - fit) magnitudes (dA == psfTry->metric) - // depend on both the true ap-fit (dAo) and the bias in the sky measurement: - // dA = dAo + dsky/flux - // where flux is the flux of the star - // we fit this trend to find the infinite flux aperture correction (dAo), - // the nominal sky bias (dsky), and the error on dAo - // the values of dA are contaminated by stars with close neighbors in the aperture - // we use an outlier rejection to avoid this bias - - // r2rflux = radius^2 * ten(0.4*fitMag); - psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - - for (int i = 0; i < psfTry->sources->n; i++) { - if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) - continue; - r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]); - } - - // XXX test dump of aperture residual data - if (psTraceGetLevel("psModules.objects") >= 5) { - FILE *f = fopen ("apresid.dat", "w"); - for (int i = 0; i < psfTry->sources->n; i++) { - int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); - - pmSource *source = psfTry->sources->data[i]; - float x = source->peak->x; - float y = source->peak->y; - - fprintf (f, "%d %d, %f %f %f %f %f %f \n", - i, keep, x, y, - psfTry->fitMag->data.F32[i], - r2rflux->data.F32[i], - psfTry->metric->data.F32[i], - psfTry->metricErr->data.F32[i]); - } - fclose (f); - } - - // This analysis of the apResid statistics is only approximate. The fitted magnitudes - // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a - // function of magnitude. We re-measure the apResid statistics later in psphot using the - // linear, constant-error fitting. Do not reject outliers with excessive vigor here. - - // fit ApTrend only to r2rflux, ignore x,y,flux variations for now - // linear clipped fit of ApResid to r2rflux - psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); - poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS) - - // XXX replace this with a psVectorStats call? since we are not fitting the trend - bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL, - psfTry->metric, psfTry->metricErr, r2rflux); - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly"); - - psFree(poly); - psFree(r2rflux); - - return false; - } - psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev - psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0], - psStatsGetValue(options->stats, stdevStat), psfTry->sources->n); - - - - // XXX test dump of fitted model (dump when tracing?) - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("resid.dat", "w"); - psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux); - for (int i = 0; i < psfTry->sources->n; i++) { - int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); - - pmSource *source = psfTry->sources->data[i]; - float x = source->peak->x; - float y = source->peak->y; - - fprintf (f, "%d %d, %f %f %f %f %f %f %f\n", - i, keep, x, y, - psfTry->fitMag->data.F32[i], - r2rflux->data.F32[i], - psfTry->metric->data.F32[i], - psfTry->metricErr->data.F32[i], - apfit->data.F32[i]); - } - fclose (f); - psFree (apfit); - } - - // XXX drop the skyBias value, or include above?? - psfTry->psf->skyBias = poly->coeff[1]; - psfTry->psf->ApResid = poly->coeff[0]; - psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat); - - psFree (r2rflux); - psFree (poly); - - return true; -} - -/* - (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) - (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) - (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) -*/ - -/***************************************************************************** -pmPSFFromPSFtry (psfTry): build a PSF model from a collection of -source->modelEXT entries. The PSF ignores the first 4 (independent) model -parameters and constructs a polynomial fit to the remaining as a function of -image coordinate. - input: psfTry with fitted source->modelEXT collection, pre-allocated psf -Note: some of the array entries may be NULL (failed fits); ignore them. - *****************************************************************************/ -bool pmPSFFromPSFtry (pmPSFtry *psfTry) -{ - PS_ASSERT_PTR_NON_NULL(psfTry, false); - PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); - - pmPSF *psf = psfTry->psf; - - // construct the fit vectors from the collection of objects - psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - - // construct the x,y terms - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelEXT == NULL) - continue; - - x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; - y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; - } - - if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { - psFree(x); - psFree(y); - psFree(z); - return false; - } - - // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) - // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters - for (int i = 0; i < psf->params->n; i++) { - switch (i) { - case PM_PAR_SKY: - case PM_PAR_I0: - case PM_PAR_XPOS: - case PM_PAR_YPOS: - case PM_PAR_SXX: - case PM_PAR_SYY: - case PM_PAR_SXY: - continue; - default: - break; - } - - // select the per-object fitted data for this parameter - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source->modelEXT == NULL) continue; - z->data.F32[j] = source->modelEXT->params->data.F32[i]; - } - - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = psf->trendNx; - binning->nYruff = psf->trendNy; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - - if (psf->psfTrendMode == PM_TREND_MAP) { - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - } - - // free existing trend, re-alloc - psFree (psf->params->data[i]); - psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); - - // fit the collection of measured parameters to the PSF 2D model - // the mask is carried from previous steps and updated with this operation - // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' - if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) { - psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); - psFree(x); - psFree(y); - psFree(z); - return false; - } - } - - // test dump of star parameters vs position (compare with fitted values) - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("params.dat", "w"); - - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source == NULL) continue; - if (source->modelEXT == NULL) continue; - - pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); - - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); - - for (int i = 0; i < psf->params->n; i++) { - if (psf->params->data[i] == NULL) continue; - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); - } - fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); - - psFree(modelPSF); - } - fclose (f); - } - - psFree (x); - psFree (y); - psFree (z); - return true; -} - - -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { - - // we are doing a robust fit. after each pass, we drop points which are more deviant than - // three sigma. mask is currently updated for each pass. - - // 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 - // each source and fit this set of parameters. These values are less tightly coupled, but - // are still inter-related. The fitted values do a good job of constraining the major axis - // and the position angle, but the minor axis is weakly measured. When we apply the PSF - // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape - // parameters, with the constraint that the minor axis must be greater than a minimum - // threshold. - - // convert the measured source shape paramters to polarization terms - psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *mag = psVectorAlloc (sources->n, PS_TYPE_F32); - - for (int i = 0; i < sources->n; i++) { - // skip any masked sources (failed to fit one of the model steps or get a magnitude) - if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - - pmSource *source = sources->data[i]; - assert (source->modelEXT); // all unmasked sources should have modelEXT - - psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); - - e0->data.F32[i] = pol.e0; - e1->data.F32[i] = pol.e1; - e2->data.F32[i] = pol.e2; - - float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; - mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; - } - - if (psf->psfTrendMode == PM_TREND_MAP) { - float scatterTotal = 0.0; - float scatterTotalMin = FLT_MAX; - int entryMin = -1; - - psVector *dz = NULL; - psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK); - - // check the fit residuals and increase Nx,Ny until the error is minimized - // pmPSFParamTrend increases the number along the longer of x or y - for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) { - - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { - break; - } - - // store the resulting scatterTotal values and the scales, redo the best - if (scatterTotal < scatterTotalMin) { - scatterTotalMin = scatterTotal; - entryMin = i; - } - } - if (entryMin == -1) { - psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params"); - return false; - } - - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { - psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); - } - - pmTrend2D *trend = psf->params->data[PM_PAR_E0]; - psf->trendNx = trend->map->map->numCols; - psf->trendNy = trend->map->map->numRows; - - // copy mask back to srcMask - for (int i = 0; i < mask->n; i++) { - srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - - psFree (mask); - psFree (dz); - } else { - - // XXX iterate Nx, Ny based on scatter? - // XXX we force the x & y order to be the same - // modify the order to correspond to the actual number of matched stars: - int order = PS_MAX (psf->trendNx, psf->trendNy); - if ((sources->n < 15) && (order >= 3)) order = 2; - if ((sources->n < 11) && (order >= 2)) order = 1; - if ((sources->n < 8) && (order >= 1)) order = 0; - if ((sources->n < 3)) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - psf->trendNx = order; - psf->trendNy = order; - - // 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 < psf->psfTrendStats->clipIter; i++) { - - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - - trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); - pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); - - trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); - pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); - - trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); - pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); - - if (!status) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - } - } - - // 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 : %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], - 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]), - srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); - } - fclose (f); - } - - psFree (e0); - psFree (e1); - psFree (e2); - psFree (mag); - return true; -} - -// fit the shape variations as a psImageMap for the given scale factor -bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { - - int Nx, Ny; - - // set the map scale to match the aspect ratio : for a scale of 1, we guarantee - // that we have a single cell - if (psf->fieldNx > psf->fieldNy) { - Nx = scale; - float AR = psf->fieldNy / (float) psf->fieldNx; - Ny = (int) (Nx * AR + 0.5); - Ny = PS_MAX (1, Ny); - } else { - Ny = scale; - float AR = psf->fieldNx / (float) psf->fieldNy; - Nx = (int) (Ny * AR + 0.5); - Nx = PS_MAX (1, Nx); - } - - // do we have enough sources for this fine of a grid? - if (x->n < 10*Nx*Ny) { - return false; - } - - // XXX check this against the exising type - pmTrend2DMode psfTrendMode = PM_TREND_MAP; - - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = Nx; - binning->nYruff = Ny; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - - psFree (psf->params->data[PM_PAR_E0]); - psFree (psf->params->data[PM_PAR_E1]); - psFree (psf->params->data[PM_PAR_E2]); - - int nIter = psf->psfTrendStats->clipIter; - psf->psfTrendStats->clipIter = 1; - psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); - - // if the map is 1x1 (a single value), we measure the resulting ensemble scatter - - // if the map is more finely sampled, divide the values into two sets: measure the fit from - // one set and the scatter from the other set. - psVector *x_fit = NULL; - psVector *y_fit = NULL; - psVector *x_tst = NULL; - psVector *y_tst = NULL; - - psVector *e0obs_fit = NULL; - psVector *e1obs_fit = NULL; - psVector *e2obs_fit = NULL; - psVector *e0obs_tst = NULL; - psVector *e1obs_tst = NULL; - psVector *e2obs_tst = NULL; - - if (scale == 1) { - x_fit = psMemIncrRefCounter (x); - y_fit = psMemIncrRefCounter (y); - x_tst = psMemIncrRefCounter (x); - y_tst = psMemIncrRefCounter (y); - e0obs_fit = psMemIncrRefCounter (e0obs); - e1obs_fit = psMemIncrRefCounter (e1obs); - e2obs_fit = psMemIncrRefCounter (e2obs); - e0obs_tst = psMemIncrRefCounter (e0obs); - e1obs_tst = psMemIncrRefCounter (e1obs); - e2obs_tst = psMemIncrRefCounter (e2obs); - } else { - x_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - x_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - for (int i = 0; i < e0obs_fit->n; i++) { - // e0obs->n == 8 or 9: - // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6 - // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7 - x_fit->data.F32[i] = x->data.F32[2*i+0]; - x_tst->data.F32[i] = x->data.F32[2*i+1]; - y_fit->data.F32[i] = y->data.F32[2*i+0]; - y_tst->data.F32[i] = y->data.F32[2*i+1]; - - e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; - e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; - e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0]; - e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1]; - e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0]; - e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1]; - } - } - - // the mask marks the values not used to calculate the ApTrend - psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK); - // copy mask values to fitMask as a starting point - for (int i = 0; i < fitMask->n; i++) { - fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - - // 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 < nIter; i++) { - // XXX we are using the same stats structure on each pass: do we need to re-init it? - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - - trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); - - trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); - - trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); - // printTrendMap (trend); - psImageMapCleanup (trend->map); - // printTrendMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); - } - psf->psfTrendStats->clipIter = nIter; // restore default setting - - // construct the fitted values and the residuals - psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst); - psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst); - psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst); - - psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit); - psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit); - psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit); - - // measure scatter for the unfitted points - // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10); - // psTraceSetLevel ("psLib.math.vectorClippedStats", 10); - pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options)); - // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0); - // psTraceSetLevel ("psLib.math.vectorClippedStats", 0); - - psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny); - psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal); - - psFree (x_fit); - psFree (y_fit); - psFree (x_tst); - psFree (y_tst); - - psFree (e0obs_fit); - psFree (e1obs_fit); - psFree (e2obs_fit); - psFree (e0obs_tst); - psFree (e1obs_tst); - psFree (e2obs_tst); - - psFree (e0fit); - psFree (e1fit); - psFree (e2fit); - - psFree (e0res); - psFree (e1res); - psFree (e2res); - - // XXX copy fitMask values back to mask - for (int i = 0; i < fitMask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - psFree (fitMask); - - return true; -} - -// calculate the scatter of the parameters -bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt) -{ - - // psStats *stats = psStatsAlloc(stdevOpt); - psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV); - - // calculate the root-mean-square of E0, E1, E2 - float dEsquare = 0.0; - psStatsInit (stats); - if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - psStatsInit (stats); - if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - psStatsInit (stats); - if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - *scatterTotal = sqrtf(dEsquare); - - psFree(stats); - return true; -} - -// calculate the minimum scatter of the parameters -bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, - psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt) -{ - - psStats *statsS = psStatsAlloc(stdevOpt); - - // measure the trend in bins with 10 values each; if < 10 total, use them all - int nBin = PS_MAX (mag->n / nGroup, 1); - - // use mag to group parameters in sequence - psVector *index = psVectorSortIndex (NULL, mag); - - // subset vectors for mag and dap values within the given range - psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK); - - int n = 0; - float min = INFINITY; // Minimum error - for (int i = 0; i < nBin; i++) { - int j; - int nValid = 0; - for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { - int N = index->data.U32[n]; - dE0subset->data.F32[j] = e0res->data.F32[N]; - dE1subset->data.F32[j] = e1res->data.F32[N]; - dE2subset->data.F32[j] = e2res->data.F32[N]; - - mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N]; - if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; - } - if (nValid < 3) continue; - - dE0subset->n = j; - dE1subset->n = j; - dE2subset->n = j; - mkSubset->n = j; - - // calculate the root-mean-square of E0, E1, E2 - float dEsquare = 0.0; - psStatsInit (statsS); - if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) { - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - psStatsInit (statsS); - if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - psStatsInit (statsS); - if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - if (isfinite(dEsquare)) { - float err = sqrtf(dEsquare); - if (err < min) { - min = err; - } - } - } - *errorFloor = min; - - psFree (dE0subset); - psFree (dE1subset); - psFree (dE2subset); - psFree (mkSubset); - - psFree(index); - - psFree(statsS); - - return true; -} Index: trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- trunk/psModules/src/objects/pmPSFtry.h (revision 25738) +++ trunk/psModules/src/objects/pmPSFtry.h (revision 25754) @@ -89,5 +89,18 @@ * */ -pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType mark); +pmPSFtry *pmPSFtryModel ( + const psArray *sources, ///< PSF sources to use in the pmPSF model analysis + const char *modelName, ///< human-readable name of desired model + pmPSFOptions *options, + psImageMaskType maskVal, + psImageMaskType mark + ); + +/** fit EXT models to all possible psf sources */ +bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); + +bool pmPSFtryMakePSF (pmPSFtry *psfTry); + +bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); /** pmPSFtryMetric() @@ -97,8 +110,5 @@ * */ -bool pmPSFtryMetric( - pmPSFtry *psfTry, ///< Add comment. - pmPSFOptions *options ///< PSF fitting options -); +bool pmPSFtryMetric(pmPSFtry *psfTry); /** pmPSFtryMetric_Alt() @@ -112,4 +122,11 @@ float RADIUS ///< Add comment. ); + +bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); + +float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); + +/// @} +# endif /** @@ -125,11 +142,3 @@ * */ -bool pmPSFFromPSFtry (pmPSFtry *psfTry); -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask); -bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz); -bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt); -bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt); - -/// @} -# endif Index: trunk/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryFitEXT.c (revision 25754) +++ trunk/psModules/src/objects/pmPSFtryFitEXT.c (revision 25754) @@ -0,0 +1,97 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// Fit an EXT model to all candidates PSF sources. +// Note: this is independent of the modeled 2D variations in the PSF. +bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { + + bool status; + + psTimerStart ("psf.fit"); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + int Next = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + if (!source->moments) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + if (!source->moments->nPixels) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + + source->modelEXT = pmSourceModelGuess (source, options->type); + if (source->modelEXT == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); + continue; + } + + // set object mask to define valid pixels + // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); + + // fit model as EXT, not PSF + status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); + + // clear object mask to define valid pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + + // exclude the poor fits + if (!status) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); + continue; + } + Next ++; + } + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n); + + if (Next == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); + return false; + } + + return true; +} Index: trunk/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 25754) +++ trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 25754) @@ -0,0 +1,149 @@ +/** @file pmPSFtryFitPSF.c + * @brief Fit the PSF model to the candidate PSF stars (part of PSF model generation) + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// stage 3: Refit with fixed shape parameters. This function uses the LMM fitting, but could +// be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c) +bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { + + bool status; + + psTimerStart ("psf.fit"); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // DEBUG code: save the PSF model fit data in detail +# ifdef DEBUG + char filename[64]; + snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy); + FILE *f = fopen (filename, "w"); + psAssert (f, "failed open"); +# endif + + int Npsf = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + psAssert (source->moments, "how can a psf source not have moments?"); + + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); + continue; + } + + // set shape for this model based on PSF + psFree (source->modelPSF); + source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); + if (source->modelPSF == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; + abort(); + continue; + } + // PSF fit and aperture mags use different radii + source->modelPSF->fitRadius = options->fitRadius; + source->apRadius = options->apRadius; + + // set object mask to define valid pixels for PSF model fit + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); + + // fit the PSF model to the source + status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_NORM, maskVal); + + // skip poor fits + if (!status) { + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); + continue; + } + + // set object mask to define valid pixels for APERTURE magnitude + if (options->fitRadius != options->apRadius) { + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal); + } + + // This function calculates the psf and aperture magnitudes + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); // raw PSF mag, AP mag + if (!status || isnan(source->apMag)) { + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); + continue; + } + + // clear object mask to define valid pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + + psfTry->fitMag->data.F32[i] = source->psfMag; + psfTry->metric->data.F32[i] = source->apMag - source->psfMag; + psfTry->metricErr->data.F32[i] = source->errMag; + + // XXX this did not work: modifies shape of psf too much + // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag; + +# ifdef DEBUG + fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", + source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], + source->psfMag, source->apMag, source->errMag, + source->modelPSF->params->data.F32[PM_PAR_I0], + source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], + source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]); +# endif + + psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); + Npsf ++; + } + psfTry->psf->nPSFstars = Npsf; + +# ifdef DEBUG + fclose (f); +# endif + + pmSourceVisualShowModelFits (psfTry->psf, psfTry->sources, maskVal); + + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, psfTry->sources->n); + + if (Npsf == 0) { + psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); + return false; + } + + return true; +} Index: trunk/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 25754) +++ trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 25754) @@ -0,0 +1,267 @@ +/** @file pmPSFtry.c + * @brief generate a pmPSF from a collection of EXT measurments of likely PSF stars. + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +/***************************************************************************** +pmPSFFromPSFtry (psfTry): build a PSF model from a collection of source->modelEXT entries +using the specified order in X,Y. The PSF ignores the first 4 (independent) model +parameters and constructs a polynomial fit to the remaining as a function of image +coordinate. Input: psfTry with fitted source->modelEXT collection, pre-allocated psf +Note: some of the array entries may be NULL (failed fits); ignore them. + *****************************************************************************/ +bool pmPSFtryMakePSF (pmPSFtry *psfTry) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + pmPSF *psf = psfTry->psf; + psVector *srcMask = psfTry->mask; + + // construct the fit vectors from the collection of objects + psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + // construct the x,y terms + for (int i = 0; i < psfTry->sources->n; i++) { + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + pmSource *source = psfTry->sources->data[i]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + } + + // fit the shape parameters (SXX, SYY, SXY) as a function of position + if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, srcMask)) { + psFree(x); + psFree(y); + return false; + } + + // vector to store the other parameter values + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + + // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY); + // fit the remaining parameters. + for (int i = 0; i < psf->params->n; i++) { + switch (i) { + case PM_PAR_SKY: + case PM_PAR_I0: + case PM_PAR_XPOS: + case PM_PAR_YPOS: + case PM_PAR_SXX: + case PM_PAR_SYY: + case PM_PAR_SXY: + continue; + default: + break; + } + + // select the per-object fitted data for this parameter + for (int j = 0; j < psfTry->sources->n; j++) { + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue; + + pmSource *source = psfTry->sources->data[j]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + z->data.F32[j] = source->modelEXT->params->data.F32[i]; + } + + pmTrend2D *trend = psf->params->data[i]; + + // fit the collection of measured parameters to the PSF 2D model + // the mask is carried from previous steps and updated with this operation + // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' + if (!pmTrend2DFit (trend, srcMask, 0xff, x, y, z, NULL)) { + psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); + psFree(x); + psFree(y); + psFree(z); + return false; + } + if (trend->mode == PM_TREND_MAP) { + // p_psImagePrint (2, trend->map->map, "param N Before"); // XXX TEST: + psImageMapRepair (trend->map->map); + // p_psImagePrint (2, trend->map->map, "param N After"); // XXX TEST: + } + } + + // test dump of star parameters vs position (compare with fitted values) + if (psTraceGetLevel("psModules.objects") >= 4) { + FILE *f = fopen ("params.dat", "w"); + + for (int j = 0; j < psfTry->sources->n; j++) { + pmSource *source = psfTry->sources->data[j]; + if (source == NULL) continue; + if (source->modelEXT == NULL) continue; + + pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); + + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + psFree(modelPSF); + continue; + } + fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); + } + fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); + + psFree(modelPSF); + } + fclose (f); + } + + psFree (x); + psFree (y); + psFree (z); + return true; +} + +// fit the shape parameters using the supplied order (pmPSF->trendNx,trendNy) +bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { + + // we are doing a robust fit. after each pass, we drop points which are more deviant than + // three sigma. the source mask (srcMask) is updated for each pass. + + // 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 + // each source and fit this set of parameters. These values are less tightly coupled, but + // are still inter-related. The fitted values do a good job of constraining the major axis + // and the position angle, but the minor axis is weakly measured. When we apply the PSF + // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape + // parameters, with the constraint that the minor axis must be greater than a minimum + // threshold. + + // XXX re-read the sextractor manual on handling 'infinitely thin' sources... + + // storage vectors for the polarization terms & mags + psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); + psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); + + // convert the measured source shape paramters to polarization terms + for (int i = 0; i < sources->n; i++) { + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + pmSource *source = sources->data[i]; + assert (source->modelEXT); // all unmasked sources should have modelEXT + + psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); + + e0->data.F32[i] = pol.e0; + e1->data.F32[i] = pol.e1; + e2->data.F32[i] = pol.e2; + } + + // 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 + // NOTE : trend->stats (below) points to the same data as psfTrendStats; set the value to 1 + // to limit the iteration within the loop + int nIter = psf->psfTrendStats->clipIter; + psf->psfTrendStats->clipIter = 1; + for (int i = 0; i < nIter; i++) { + // XXX we are using the same stats structure on each pass: do we need to re-init it? + // XXX we hardwire this to SAMPLE stats above (psphotChoosePSF.c), hardwire here instead? + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + + trend = psf->params->data[PM_PAR_E0]; + trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); + + trend = psf->params->data[PM_PAR_E1]; + trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); + + trend = psf->params->data[PM_PAR_E2]; + trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); + if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); + + if (!status) { + psError (PS_ERR_UNKNOWN, true, "failed to fit PSF shape params"); + return false; + } + } + psf->psfTrendStats->clipIter = nIter; + + // 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 : %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], + 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]), + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); + } + fclose (f); + } + + psFree (e0); + psFree (e1); + psFree (e2); + return true; +} + Index: trunk/psModules/src/objects/pmPSFtryMetric.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryMetric.c (revision 25754) +++ trunk/psModules/src/objects/pmPSFtryMetric.c (revision 25754) @@ -0,0 +1,84 @@ +/** @file pmPSFtry.c + * @brief: measure the systematic error in the aperture-psf magnitude + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// The quality of the PSF model is determined by the systematic scatter in the fit - aperture +// magnitude. We use the moments->Sum value, calculated with a Gaussian window, as a proxy for +// the aperture magnitude. We measure the systmatic scatter with the function +// psVectorSystematicError which solves for the value of SysErr that needs to be added in +// quadrature to the errors of a set of residual measurements in order to yield a ChiSq of 1.0. +bool pmPSFtryMetric (pmPSFtry *psfTry) +{ + PS_ASSERT_PTR_NON_NULL(psfTry, false); + PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + if (!psVectorStats (stats, psfTry->metric, NULL, psfTry->mask, PSFTRY_MASK_ALL)) { + psError(PS_ERR_UNKNOWN, false, "Failed to measured clipped mean"); + psFree (stats); + return false; + } + + // generate a residual vector + psVector *resid = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + psVector *error = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); + int n = 0; + for (int i = 0; i < psfTry->metric->n; i++) { + if (psfTry->mask && (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] && PSFTRY_MASK_ALL)) continue; + // if (psfTry->metricErr->data.F32[i] > 0.005) continue; + resid->data.F32[n] = psfTry->metric->data.F32[i] - stats->sampleMedian; + error->data.F32[n] = psfTry->metricErr->data.F32[i]; + n++; + } + resid->n = error->n = n; + + float psfSysErr = psVectorSystematicError (resid, error, 0.05); + + psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f (systematic error) from statistics of %ld psf stars (%d used)\n", + stats->sampleMedian, psfSysErr, psfTry->sources->n, n); + + psfTry->psf->ApResid = stats->sampleMedian; + psfTry->psf->dApResid = psfSysErr; + + pmSourceVisualPlotPSFMetric (psfTry); + + psFree (stats); + psFree (resid); + psFree (error); + + return true; +} Index: trunk/psModules/src/objects/pmPSFtryModel.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryModel.c (revision 25754) +++ trunk/psModules/src/objects/pmPSFtryModel.c (revision 25754) @@ -0,0 +1,214 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +# include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmPSFtry.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" +#include "pmSourceVisual.h" + +// build a pmPSFtry for the given model: +// - fit each source with the free-floating model +// - construct the pmPSF from the collection of models +// - fit each source with the PSF-parameter models +// - measure the pmPSF quality metric (dApResid) + +// sources used in for pmPSFtry may be masked by the analysis +// mask values indicate the reason the source was rejected: + +// generate a pmPSFtry with a copy of the test PSF sources +pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) +{ + // 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"); + return NULL; + } + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF + if (!pmPSFtryFitEXT(psfTry, options, maskVal, markVal)) { + psError(PS_ERR_UNKNOWN, false, "failed to fit EXT models to sources for psf model"); + psFree(psfTry); + return NULL; + } + + int orderMin = (options->psfTrendMode == PM_TREND_MAP) ? 1 : 0; + int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy); + + // XXX set the min number of needed source more carefully (depends on psfTrendMode?) + if ((sources->n < 15) && (orderMax >= 3)) orderMax = 2; + if ((sources->n < 11) && (orderMax >= 2)) orderMax = 1; + if ((sources->n < 8) && (orderMax >= 1)) orderMax = 0; + if ((sources->n < 3)) { + psError (PS_ERR_UNKNOWN, true, "failed to determine PSF parameters"); + return NULL; + } + + // save the raw source mask (generated by pmPSFtryFitEXT) + psVector *srcMask = psVectorCopy (NULL, psfTry->mask, PS_TYPE_VECTOR_MASK); + + // we will save the PSF with the best fit (min systematic error) + pmPSF *minPSF = NULL; + psVector *minMask = NULL; + + // as we loop over orders, we need to refer to the initial selection, but we modify the + // option values to match the current guess: save the max values here: + int Nx = options->psfTrendNx; + int Ny = options->psfTrendNy; + for (int i = orderMin; i <= orderMax; i++) { + + if (Nx > Ny) { + options->psfTrendNx = i; + options->psfTrendNy = PS_MAX (orderMin, (int)(i * (Ny / Nx) + 0.5)); + } else { + options->psfTrendNy = i; + options->psfTrendNx = PS_MAX (orderMin, (int)(i * (Nx / Ny) + 0.5)); + } + + // free existing data, if any + psFree(psfTry->psf); + psFree(psfTry->mask); + + // allocate a mask and a psf model using the current Nx,Ny order values; + psfTry->psf = pmPSFAlloc (options); + psfTry->mask = psVectorCopy (NULL, srcMask, PS_TYPE_VECTOR_MASK); + + // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation + if (!pmPSFtryMakePSF (psfTry)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + + // stage 3: refit with fixed shape parameters, measure pmPSFtry->metric + if (!pmPSFtryFitPSF (psfTry, options, maskVal, markVal)) { + psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); + psFree(psfTry); + return NULL; + } + + // stage 4: measure systematic error in pmPSFtry->metric + if (!pmPSFtryMetric (psfTry)) { + psError(PS_ERR_UNKNOWN, false, "failed to measure systematic error of metric"); + psFree(psfTry); + return NULL; + } + + if (!minPSF) { + minPSF = psMemIncrRefCounter(psfTry->psf); + minMask = psMemIncrRefCounter(psfTry->mask); + } + + if (psfTry->psf->dApResid < minPSF->dApResid) { + psFree (minPSF); + psFree (minMask); + minPSF = psMemIncrRefCounter(psfTry->psf); + minMask = psMemIncrRefCounter(psfTry->mask); + } + } + psFree (srcMask); + + // keep the ones matching the min systematic error: + psFree (psfTry->psf); + psFree (psfTry->mask); + psfTry->psf = minPSF; + psfTry->mask = minMask; + + // XXXXX this is probably not used any more. Are the chisq of the fits so bad? can we + // fix them by softening the errors on the brightest pixels? + + // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) + // this should be linear for Poisson errors and quadratic for constant sky errors + psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK); + + // generate the x and y vectors, and mask missing models + for (int i = 0; i < psfTry->sources->n; i++) { + pmSource *source = psfTry->sources->data[i]; + if (source->modelPSF == NULL) { + flux->data.F32[i] = 0.0; + chisq->data.F32[i] = 0.0; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; + } else { + flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; + chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + } + } + + // use 3hi/3lo sigma clipping on the chisq fit + psStats *stats = options->stats; + + // linear clipped fit of chisq trend vs flux + if (options->chiFluxTrend) { + bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, + mask, 0xff, chisq, NULL, flux); + psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean + psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev + + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", + psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); + + psFree(flux); + psFree(mask); + psFree(chisq); + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); + psFree(psfTry); + return NULL; + } + } + + for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) { + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, + psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), + psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i)); + } + + psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", + modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); + + return (psfTry); +} + Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 25738) +++ trunk/psModules/src/objects/pmPeaks.c (revision 25754) @@ -60,7 +60,8 @@ // if min point is too deviant, use the peak value // XXX need to calculate dx, dy correctly + // 0.5 PIX: peaks are calculated using the pixel index and converted here to pixel coords if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { - peak->xf = min.x + ix + image->col0; - peak->yf = min.y + iy + image->row0; + peak->xf = min.x + ix + image->col0 + 0.5; + peak->yf = min.y + iy + image->row0 + 0.5; // These errors are fractional errors, and should be scaled by the @@ -73,6 +74,6 @@ peak->yf = PS_MAX (PS_MIN (peak->yf, image->numRows - 1), image->row0); } else { - peak->xf = ix; - peak->yf = iy; + peak->xf = ix + 0.5; + peak->yf = iy + 0.5; peak->dx = NAN; peak->dy = NAN; Index: trunk/psModules/src/objects/pmPetrosian.c =================================================================== --- trunk/psModules/src/objects/pmPetrosian.c (revision 25754) +++ trunk/psModules/src/objects/pmPetrosian.c (revision 25754) @@ -0,0 +1,107 @@ +/* @file pmPetrosian.c + * low-level petrosian functions + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pmPetrosian.h" + +static void pmPetrosianFree(pmPetrosian *petrosian) +{ + if (!petrosian) { + return; + } + psFree(petrosian->radii); + psFree(petrosian->fluxes); + psFree(petrosian->theta); + psFree(petrosian->isophotalRadii); + + psFree(petrosian->radiusElliptical); + psFree(petrosian->fluxElliptical); + + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); + psFree(petrosian->radialBins); + psFree(petrosian->area); +} + +pmPetrosian *pmPetrosianAlloc() +{ + pmPetrosian *petrosian = (pmPetrosian *)psAlloc(sizeof(pmPetrosian)); + psMemSetDeallocator(petrosian, (psFreeFunc) pmPetrosianFree); + + petrosian->radii = NULL; + petrosian->fluxes = NULL; + petrosian->theta = NULL; + petrosian->isophotalRadii = NULL; + + petrosian->radiusElliptical = NULL; + petrosian->fluxElliptical = NULL; + + petrosian->radialBins = NULL; + petrosian->area = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; + + petrosian->petrosianRadius = NAN; + petrosian->petrosianFlux = NAN; + + // petrosian->axes = {0.0, 0.0, 0.0}; + + return petrosian; +} + +bool pmPetrosianFreeVectors(pmPetrosian *petrosian) { + + psFree(petrosian->radii); + psFree(petrosian->fluxes); + psFree(petrosian->theta); + psFree(petrosian->isophotalRadii); + + psFree(petrosian->radiusElliptical); + psFree(petrosian->fluxElliptical); + + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); + psFree(petrosian->radialBins); + psFree(petrosian->area); + + petrosian->radii = NULL; + petrosian->fluxes = NULL; + petrosian->theta = NULL; + petrosian->isophotalRadii = NULL; + + petrosian->radiusElliptical = NULL; + petrosian->fluxElliptical = NULL; + + petrosian->radialBins = NULL; + petrosian->area = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; + + return true; +} + +# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B]) +# define SWAP_INDEX(TYPE,A,B) { \ + float tmp; \ + if (A != B) { \ + tmp = index->data.F32[A]; \ + index->data.F32[A] = index->data.F32[B]; \ + index->data.F32[B] = tmp; \ + tmp = extra->data.F32[A]; \ + extra->data.F32[A] = extra->data.F32[B]; \ + extra->data.F32[B] = tmp; \ + } \ +} + +bool pmPetrosianSortPair (psVector *index, psVector *extra) { + + // sort the vector set by the radius + PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE); + return true; +} Index: trunk/psModules/src/objects/pmPetrosian.h =================================================================== --- trunk/psModules/src/objects/pmPetrosian.h (revision 25754) +++ trunk/psModules/src/objects/pmPetrosian.h (revision 25754) @@ -0,0 +1,44 @@ +/* @file pmPetrosian.h + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_PETROSIAN_H +#define PM_PETROSIAN_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +typedef struct { + psArray *radii; // radii for raw radial profiles at evenly-spaced angles + psArray *fluxes; // fluxes measured at above radii + psVector *theta; // angles corresponding to above radial profiles + psVector *isophotalRadii; // isophotal radius for the above angles + + psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels + psVector *fluxElliptical; // flux for the above radial coordinates + + psVector *binSB; // mean surface brightness within radial bins + psVector *binSBstdev; // scatter of mean surface brightness within radial bins + psVector *radialBins; // radii corresponding to above binnedBlux + psVector *area; // differential area of the non-overlapping radial bins + + psEllipseAxes axes; // shape of elliptical contour + + float petrosianRadius; + float petrosianFlux; + +} pmPetrosian; + +pmPetrosian *pmPetrosianAlloc(); + +bool pmPetrosianFreeVectors(pmPetrosian *petrosian); +bool pmPetrosianSortPair (psVector *index, psVector *extra); + +/// @} + +# endif /* PM_PETROSIAN_H */ Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 25738) +++ trunk/psModules/src/objects/pmSource.c (revision 25754) @@ -282,6 +282,6 @@ psArray *peaks = NULL; - pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0}; - pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0}; + pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0, 0, 0.0}; + pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0, 0, 0.0}; pmPSFClump psfClump; @@ -319,5 +319,6 @@ // construct a sigma-plane image - int numCols = SX_MAX / PSF_CLUMP_GRID_SCALE, numRows = SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image + int numCols = 1 + SX_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image + int numRows = 1 + SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image psTrace("psModules.objects", 10, "sigma-plane dimensions: %dx%d\n", numCols, numRows); psImage *splane = psImageAlloc(numCols, numRows, PS_TYPE_F32); // sigma-plane image @@ -332,12 +333,16 @@ } - int x = src->peak->x, y = src->peak->y; // Coordinates of peak - if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) { - continue; - } + if (region) { + int x = src->peak->x, y = src->peak->y; // Coordinates of peak + if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) { + continue; + } + } if (src->mode & PM_SOURCE_MODE_BLEND) { continue; } + + if (!src->moments->nPixels) continue; if (src->moments->SN < PSF_SN_LIM) { @@ -384,5 +389,5 @@ // find the peak in this image - psStats *stats = psStatsAlloc (PS_STAT_MAX); + psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV); if (!psImageStats (stats, splane, NULL, 0)) { psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); @@ -393,4 +398,6 @@ peaks = pmPeaksInImage (splane, stats[0].max / 2); psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2); + + psfClump.nSigma = stats->sampleStdev; const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP"); @@ -430,4 +437,9 @@ psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value); + // XXX store the mean sigma? + float meanSigma = psfClump.nSigma; + psfClump.nStars = clump->value; + psfClump.nSigma = clump->value / meanSigma; + // define section window for clump minSx = clump->x * PSF_CLUMP_GRID_SCALE - 2.0*PSF_CLUMP_GRID_SCALE; @@ -452,8 +464,10 @@ continue; - if (tmpSrc->peak->x < region->x0) continue; - if (tmpSrc->peak->x > region->x1) continue; - if (tmpSrc->peak->y < region->y0) continue; - if (tmpSrc->peak->y > region->y1) continue; + if (region) { + if (tmpSrc->peak->x < region->x0) continue; + if (tmpSrc->peak->x > region->x1) continue; + if (tmpSrc->peak->y < region->y0) continue; + if (tmpSrc->peak->y > region->y1) continue; + } if (tmpSrc->moments->Mxx < minSx) @@ -919,4 +933,6 @@ if (model == NULL) return false; // model must be defined + bool addNoise = mode & PM_MODEL_OP_NOISE; + if (source->modelFlux) { // add in the pixels from the modelFlux image @@ -940,10 +956,9 @@ psF32 **target = source->pixels->data.F32; - if (mode & PM_MODEL_OP_NOISE) { - // XXX need to scale by the gain... + if (addNoise) { + // when adding noise, we assume the shape and Io have been modified target = source->variance->data.F32; } - // XXX need to respect the source and model masks? for (int iy = 0; iy < source->modelFlux->numRows; iy++) { int oy = iy + dY; @@ -959,16 +974,27 @@ } } + if (!addNoise) { + if (add) { + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + } else { + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + } + } return true; } psImage *target = source->pixels; - if (mode & PM_MODEL_OP_NOISE) { + if (addNoise) { target = source->variance; } - if (add) { - status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); - } else { - status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); + if (!addNoise) { + if (add) { + status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + } else { + status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + } } Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 25738) +++ trunk/psModules/src/objects/pmSource.h (revision 25754) @@ -38,6 +38,7 @@ typedef enum { - PM_SOURCE_TMPF_MODEL_GUESS = 0x0001, - PM_SOURCE_TMPF_SUBTRACTED = 0x0002, + PM_SOURCE_TMPF_MODEL_GUESS = 0x0001, + PM_SOURCE_TMPF_SUBTRACTED = 0x0002, + PM_SOURCE_TMPF_SIZE_MEASURED = 0x0004, } pmSourceTmpF; @@ -81,4 +82,5 @@ float extNsigma; ///< Nsigma deviation from PSF to EXT float sky, skyErr; ///< The sky and its error at the center of the object + float apRadius; psRegion region; ///< area on image covered by selected pixels pmSourceExtendedPars *extpars; ///< extended source parameters @@ -98,4 +100,6 @@ float Y; float dY; + int nStars; + float nSigma; } pmPSFClump; Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 25754) @@ -17,30 +17,140 @@ #endif -#include -#include -#include +// #include +// #include +// #include #include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPAMaskWeight.h" -#include "pmSpan.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" -#include "pmModel.h" -#include "pmSource.h" - +// #include "pmHDU.h" +// #include "pmFPA.h" +// #include "pmFPAMaskWeight.h" +// #include "pmSpan.h" +// #include "pmFootprint.h" +// #include "pmPeaks.h" +// #include "pmMoments.h" +// #include "pmResiduals.h" +// #include "pmGrowthCurve.h" +// #include "pmTrend2D.h" +// #include "pmPSF.h" +// #include "pmModel.h" +// #include "pmSource.h" +#include "pmSourceExtendedPars.h" + +// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and +// intermediate data used to measure the profile +static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile) +{ + if (!profile) { + return; + } + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + psFree(profile->isophotalRadii); + + psFree(profile->radiusElliptical); + psFree(profile->fluxElliptical); + + psFree(profile->binSB); + psFree(profile->binSBstdev); + psFree(profile->binSBerror); + + psFree(profile->radialBins); + psFree(profile->area); +} + +pmSourceRadialProfile *pmSourceRadialProfileAlloc() +{ + pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile)); + psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); + + profile->radii = NULL; + profile->fluxes = NULL; + profile->theta = NULL; + profile->isophotalRadii = NULL; + + profile->radiusElliptical = NULL; + profile->fluxElliptical = NULL; + + profile->binSB = NULL; + profile->binSBstdev = NULL; + profile->binSBerror = NULL; + + profile->radialBins = NULL; + profile->area = NULL; + + return profile; +} + +bool psMemCheckSourceRadialProfile(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree); +} + + +// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) { + + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + psFree(profile->isophotalRadii); + + psFree(profile->radiusElliptical); + psFree(profile->fluxElliptical); + + // psFree(profile->binSB); + // psFree(profile->binSBstdev); + // psFree(profile->binSBerror); + + // psFree(profile->radialBins); + psFree(profile->area); + + profile->radii = NULL; + profile->fluxes = NULL; + profile->theta = NULL; + profile->isophotalRadii = NULL; + + profile->radiusElliptical = NULL; + profile->fluxElliptical = NULL; + + // profile->binSB = NULL; + // profile->binSBstdev = NULL; + // profile->binSBerror = NULL; + + // profile->radialBins = NULL; + profile->area = NULL; + + return true; +} + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B]) +# define SWAP_INDEX(TYPE,A,B) { \ + float tmp; \ + if (A != B) { \ + tmp = index->data.F32[A]; \ + index->data.F32[A] = index->data.F32[B]; \ + index->data.F32[B] = tmp; \ + tmp = extra->data.F32[A]; \ + extra->data.F32[A] = extra->data.F32[B]; \ + extra->data.F32[B] = tmp; \ + } \ +} + +bool pmSourceRadialProfileSortPair (psVector *index, psVector *extra) { + + // sort the vector set by the radius + PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE); + return true; +} + +// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source static void pmSourceExtendedParsFree (pmSourceExtendedPars *pars) { if (!pars) return; psFree(pars->profile); - psFree(pars->annuli); - psFree(pars->isophot); - psFree(pars->petrosian); - psFree(pars->kron); + psFree(pars->petrosian_50); + psFree(pars->petrosian_80); return; } @@ -51,8 +161,6 @@ pars->profile = NULL; - pars->annuli = NULL; - pars->isophot = NULL; - pars->petrosian = NULL; - pars->kron = NULL; + pars->petrosian_50 = NULL; + pars->petrosian_80 = NULL; return pars; @@ -66,135 +174,27 @@ -static void pmSourceRadialProfileFree (pmSourceRadialProfile *profile) { - if (!profile) return; - - psFree(profile->radius); - psFree(profile->flux); - psFree(profile->variance); +// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind +static void pmSourceExtendedFluxFree (pmSourceExtendedFlux *flux) { + if (!flux) return; return; } -pmSourceRadialProfile *pmSourceRadialProfileAlloc (void) { - - pmSourceRadialProfile *profile = (pmSourceRadialProfile *) psAlloc(sizeof(pmSourceRadialProfile)); - psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); - - profile->radius = NULL; - profile->flux = NULL; - profile->variance = NULL; - - return profile; -} - -bool psMemCheckSourceRadialProfile(psPtr ptr) +pmSourceExtendedFlux *pmSourceExtendedFluxAlloc (void) { + + pmSourceExtendedFlux *flux = (pmSourceExtendedFlux *) psAlloc(sizeof(pmSourceExtendedFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceExtendedFluxFree); + + flux->flux = 0.0; + flux->fluxErr = 0.0; + flux->radius = 0.0; + flux->radiusErr = 0.0; + + return flux; +} + + +bool psMemCheckSourceExtendedFlux(psPtr ptr) { PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree); -} - - -static void pmSourceIsophotalValuesFree (pmSourceIsophotalValues *isophot) { - if (!isophot) return; - return; -} - -pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc (void) { - - pmSourceIsophotalValues *isophot = (pmSourceIsophotalValues *) psAlloc(sizeof(pmSourceIsophotalValues)); - psMemSetDeallocator(isophot, (psFreeFunc) pmSourceIsophotalValuesFree); - - isophot->mag = 0.0; - isophot->magErr = 0.0; - isophot->rad = 0.0; - isophot->radErr = 0.0; - - return isophot; -} - - -bool psMemCheckSourceIsophotalValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceIsophotalValuesFree); -} - - -static void pmSourcePetrosianValuesFree (pmSourcePetrosianValues *petrosian) { - if (!petrosian) return; - return; -} - -pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc (void) { - - pmSourcePetrosianValues *petrosian = (pmSourcePetrosianValues *) psAlloc(sizeof(pmSourcePetrosianValues)); - psMemSetDeallocator(petrosian, (psFreeFunc) pmSourcePetrosianValuesFree); - - petrosian->mag = 0.0; - petrosian->magErr = 0.0; - petrosian->rad = 0.0; - petrosian->radErr = 0.0; - - return petrosian; -} - - -bool psMemCheckSourcePetrosianValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourcePetrosianValuesFree); -} - -static void pmSourceKronValuesFree (pmSourceKronValues *kron) { - if (!kron) return; - return; -} - -pmSourceKronValues *pmSourceKronValuesAlloc (void) { - - pmSourceKronValues *kron = (pmSourceKronValues *) psAlloc(sizeof(pmSourceKronValues)); - psMemSetDeallocator(kron, (psFreeFunc) pmSourceKronValuesFree); - - kron->mag = 0.0; - kron->magErr = 0.0; - kron->rad = 0.0; - kron->radErr = 0.0; - - return kron; -} - - -bool psMemCheckSourceKronValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceKronValuesFree); -} - - -static void pmSourceAnnuliFree (pmSourceAnnuli *annuli) { - if (!annuli) return; - - psFree (annuli->flux); - psFree (annuli->fluxErr); - psFree (annuli->fluxVar); - - return; -} - -pmSourceAnnuli *pmSourceAnnuliAlloc (void) { - - pmSourceAnnuli *annuli = (pmSourceAnnuli *) psAlloc(sizeof(pmSourceAnnuli)); - psMemSetDeallocator(annuli, (psFreeFunc) pmSourceAnnuliFree); - - annuli->flux = NULL; - annuli->fluxErr = NULL; - annuli->fluxVar = NULL; - - return annuli; -} - - -bool psMemCheckSourceAnnuli(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceAnnuliFree); -} + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree); +} Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 25738) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 25754) @@ -15,56 +15,53 @@ typedef struct { - psVector *radius; - psVector *flux; - psVector *variance; + psArray *radii; // radii for raw radial profiles at evenly-spaced angles + psArray *fluxes; // fluxes measured at above radii + psVector *theta; // angles corresponding to above radial profiles + psVector *isophotalRadii; // isophotal radius for the above angles + + psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels + psVector *fluxElliptical; // flux for the above radial coordinates + + psVector *binSB; // mean surface brightness within radial bins + psVector *binSBstdev; // scatter of mean surface brightness within radial bins + psVector *binSBerror; // formal error on mean surface brightness within radial bins + + psVector *radialBins; // radii corresponding to above binnedBlux + psVector *area; // differential area of the non-overlapping radial bins + + psEllipseAxes axes; // shape of elliptical contour } pmSourceRadialProfile; typedef struct { - psVector *flux; - psVector *fluxErr; - psVector *fluxVar; -} pmSourceAnnuli; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourceIsophotalValues; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourcePetrosianValues; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourceKronValues; + float flux; + float fluxErr; + float radius; + float radiusErr; +} pmSourceExtendedFlux; typedef struct { pmSourceRadialProfile *profile; - pmSourceAnnuli *annuli; - pmSourceIsophotalValues *isophot; - pmSourcePetrosianValues *petrosian; - pmSourceKronValues *kron; + pmSourceExtendedFlux *petrosian_50; + pmSourceExtendedFlux *petrosian_80; } pmSourceExtendedPars; -pmSourceExtendedPars *pmSourceExtendedParsAlloc(void); +// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and +// intermediate data used to measure the profile +pmSourceRadialProfile *pmSourceRadialProfileAlloc(); +bool psMemCheckSourceRadialProfile(psPtr ptr); + +// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile); + +// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source +pmSourceExtendedPars *pmSourceExtendedParsAlloc (void); bool psMemCheckSourceExtendedPars(psPtr ptr); -pmSourceRadialProfile *pmSourceRadialProfileAlloc(void); -bool psMemCheckSourceRadialProfile(psPtr ptr); -pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc(void); -bool psMemCheckSourceIsophotalValues(psPtr ptr); -pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc(void); -bool psMemCheckSourcePetrosianValues(psPtr ptr); -pmSourceKronValues *pmSourceKronValuesAlloc(void); -bool psMemCheckSourceKronValues(psPtr ptr); -pmSourceAnnuli *pmSourceAnnuliAlloc(void); -bool psMemCheckSourceAnnuli(psPtr ptr); + +// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind +pmSourceExtendedFlux *pmSourceExtendedFluxAlloc(void); +bool psMemCheckSourceExtendedFlux(psPtr ptr); + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra); /// @} Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 25754) @@ -96,6 +96,7 @@ // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); + // 0.5 PIX: the coordinate values must be in pixel coords, not index + coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0); x->data[nPix] = (psPtr *) coord; y->data.F32[nPix] = source->pixels->data.F32[i][j]; @@ -175,7 +176,5 @@ continue; dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); - if (psTraceGetLevel("psModules.objects") >= 4) { - fprintf (stderr, "%f +/- %f\n", params->data.F32[i], dparams->data.F32[i]); - } + psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]); } psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); @@ -186,4 +185,5 @@ model->nPix = y->n; model->nDOF = y->n - nParams; + model->chisqNorm = model->chisq / model->nDOF; model->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 25754) @@ -321,9 +321,7 @@ for (int j = 0; j < model->params->n; j++, n++) { - if (psTraceGetLevel("psModules.objects") >= 4) { - fprintf (stderr, "%f ", param->data.F32[n]); - } model->params->data.F32[j] = param->data.F32[n]; model->dparams->data.F32[j] = dparam->data.F32[n]; + psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]); } psTrace ("psModules.objects", 4, " src %d", i); @@ -483,6 +481,7 @@ // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); + // 0.5 PIX: the coordinate values must be in pixel coords, not index + coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0); x->data[nPix] = (psPtr *) coord; y->data.F32[nPix] = source->pixels->data.F32[i][j]; Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 25754) @@ -572,4 +572,5 @@ psFree (xfitname); psFree (outhead); + psFree (deteffname); return false; } @@ -583,4 +584,5 @@ psFree (xfitname); psFree (outhead); + psFree (deteffname); break; @@ -938,4 +940,5 @@ psFree (headname); psFree (dataname); + psFree (deteffname); return true; } @@ -996,4 +999,5 @@ psFree (headname); psFree (dataname); + psFree (deteffname); psFree (tableHeader); break; Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 25754) @@ -128,5 +128,5 @@ nDOF = model->nDOF; nPix = model->nPix; - apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + apRadius = source->apRadius; errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { @@ -308,4 +308,5 @@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -313,5 +314,4 @@ model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); source->moments = pmMomentsAlloc (); @@ -353,8 +353,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -396,4 +396,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -476,4 +477,5 @@ } +# endif psArrayAdd (table, 100, row); psFree (row); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25754) @@ -128,5 +128,5 @@ nDOF = model->nDOF; nPix = model->nPix; - apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + apRadius = source->apRadius; errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { @@ -308,4 +308,5 @@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -313,5 +314,4 @@ model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); source->moments = pmMomentsAlloc (); @@ -353,8 +353,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -396,4 +396,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -476,4 +477,5 @@ } +# endif psArrayAdd (table, 100, row); psFree (row); Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 25754) @@ -349,8 +349,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -410,4 +410,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -501,4 +502,5 @@ } } +# endif psArrayAdd (table, 100, row); Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 25754) @@ -300,8 +300,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -343,4 +343,6 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +// XXX disable these outputs until we clean up the names +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -422,4 +424,5 @@ } } +# endif psArrayAdd (table, 100, row); Index: trunk/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 25754) @@ -119,5 +119,5 @@ logChi, logChiNorm, source[0].peak->SN, - model[0].radiusFit, + source[0].apRadius, source[0].pixWeight, model[0].nDOF, @@ -181,5 +181,5 @@ log10(model[0].chisqNorm/model[0].nDOF), source[0].peak->SN, - model[0].radiusFit, + source[0].apRadius, source[0].pixWeight, model[0].nDOF, Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 25754) @@ -36,6 +36,4 @@ #include "pmSource.h" -bool pmSourceMoments_Old(pmSource *source, psF32 radius); - /****************************************************************************** pmSourceMoments(source, radius): this function takes a subImage defined in the @@ -50,8 +48,7 @@ pmSource->mask (optional) -XXX: The peak calculations are done in image coords, not subImage coords. - -this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons +The peak calculations are done in image coords, not subImage coords. + +this version optionally clips input pixels on S/N *****************************************************************************/ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) @@ -64,5 +61,5 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // XXX supply sky in a different way? + // use sky from moments if defined, 0.0 otherwise psF32 sky = 0.0; if (source->moments == NULL) { @@ -91,10 +88,18 @@ // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in // this subimage. we subtract off the peak coordinates, adjusted to this subimage, to have - // minimal round-off error in the sums - - psF32 xOff = source->peak->x; - psF32 yOff = source->peak->y; - psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage - psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage + // minimal round-off error in the sums. since these values are subtracted just to minimize + // the dynamic range and are added back below, the exact value does not matter. these are + // (int) so they can be used in the image index below. + + int xOff = source->peak->x; + int yOff = source->peak->y; + int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage + int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage + + // we are guaranteed to have a valid pixel and variance at this location (right? right?) + // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -126,21 +131,20 @@ psF32 wDiff = *vWgt; - // skip pixels below specified significance level + // skip pixels below specified significance level. this is allowed, but should be + // avoided -- the over-weights the wings of bright stars compared to those of faint + // stars. if (PS_SQR(pDiff) < minSN2*wDiff) continue; if (pDiff < 0) continue; + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources if (sigma > 0.0) { - // apply a pseudo-gaussian weight - - // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - assert (z >= 0.0); - psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); - psF32 weight = 1.0 / t; - - // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); - - wDiff *= weight; - pDiff *= weight; + // XXX a lot of extra flops; can we pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; } @@ -154,9 +158,7 @@ Y1 += yWght; - // fprintf (stderr, " : %6.1f %6.1f %8.1f %8.1f\n", X1, Y1, Sum, Var); - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + } } @@ -164,6 +166,6 @@ // XXX EAM - the limit is a bit arbitrary. make it user defined? if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); - return (false); + psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); + return (false); } @@ -172,6 +174,6 @@ float My = Y1/Sum; if ((fabs(Mx) > radius) || (fabs(My) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); - return (false); + psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); + return (false); } @@ -179,6 +181,9 @@ // add back offset of peak in primary image - source->moments->Mx = Mx + xOff; - source->moments->My = My + yOff; + // also offset from pixel index to pixel coordinate + // (the calculation above uses pixel index instead of coordinate) + // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords + source->moments->Mx = Mx + xOff + 0.5; + source->moments->My = My + yOff + 0.5; source->moments->Sum = Sum; @@ -205,5 +210,6 @@ Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum - // center of mass in subimage + // center of mass in subimage. Note: the calculation below uses pixel index, so we do not + // correct xCM, yCM to pixel coordinate here. psF32 xCM = Mx + xPeak; // coord of peak in subimage psF32 yCM = My + yPeak; // coord of peak in subimage @@ -211,87 +217,85 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row - yCM; + psF32 yDiff = row - yCM; if (fabs(yDiff) > radius) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col - xCM; + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if (vMsk) { + if (*vMsk) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + psF32 xDiff = col - xCM; if (fabs(xDiff) > radius) continue; - // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - psF32 r = sqrt(r2); - if (r > radius) continue; - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : should this limit be user-defined? - - if (PS_SQR(pDiff) < minSN2*wDiff) continue; - if (pDiff < 0) continue; - + // radius is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + psF32 r = sqrt(r2); + if (r > radius) continue; + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // skip pixels below specified significance level. this is allowed, but should be + // avoided -- the over-weights the wings of bright stars compared to those of faint + // stars. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + if (pDiff < 0) continue; + + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources if (sigma > 0.0) { - // apply a pseudo-gaussian weight - - // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - assert (z >= 0.0); - psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); - psF32 weight = 1.0 / t; - - // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); - - wDiff *= weight; - pDiff *= weight; + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; } - Sum += pDiff; - - psF32 x = xDiff * pDiff; - psF32 y = yDiff * pDiff; - - psF32 xx = xDiff * x; - psF32 xy = xDiff * y; - psF32 yy = yDiff * y; - - psF32 xxx = xDiff * xx / r; - psF32 xxy = xDiff * xy / r; - psF32 xyy = xDiff * yy / r; - psF32 yyy = yDiff * yy / r; - - psF32 xxxx = xDiff * xxx / r2; - psF32 xxxy = xDiff * xxy / r2; - psF32 xxyy = xDiff * xyy / r2; - psF32 xyyy = xDiff * yyy / r2; - psF32 yyyy = yDiff * yyy / r2; - - XX += xx; - XY += xy; - YY += yy; - - XXX += xxx; - XXY += xxy; - XYY += xyy; - YYY += yyy; - - XXXX += xxxx; - XXXY += xxxy; - XXYY += xxyy; - XYYY += xyyy; - YYYY += yyyy; - } + Sum += pDiff; + + psF32 x = xDiff * pDiff; + psF32 y = yDiff * pDiff; + + psF32 xx = xDiff * x; + psF32 xy = xDiff * y; + psF32 yy = yDiff * y; + + psF32 xxx = xDiff * xx / r; + psF32 xxy = xDiff * xy / r; + psF32 xyy = xDiff * yy / r; + psF32 yyy = yDiff * yy / r; + + psF32 xxxx = xDiff * xxx / r2; + psF32 xxxy = xDiff * xxy / r2; + psF32 xxyy = xDiff * xyy / r2; + psF32 xyyy = xDiff * yyy / r2; + psF32 yyyy = yDiff * yyy / r2; + + XX += xx; + XY += xy; + YY += yy; + + XXX += xxx; + XXY += xxy; + XYY += xyy; + YYY += yyy; + + XXXX += xxxx; + XXXY += xxxy; + XXYY += xxyy; + XYYY += xyyy; + YYYY += yyyy; + } } @@ -312,151 +316,18 @@ if (source->moments->Mxx < 0) { - fprintf (stderr, "error: neg second moment??\n"); + fprintf (stderr, "error: neg second moment??\n"); } if (source->moments->Myy < 0) { - fprintf (stderr, "error: neg second moment??\n"); + fprintf (stderr, "error: neg second moment??\n"); } psTrace ("psModules.objects", 4, "Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", - source->moments->Mxx, source->moments->Mxy, source->moments->Myy, - source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, - source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, + source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", - source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, numPixels); - - if (source->moments->Mxx < 0) { - fprintf (stderr, "error: neg second moment??\n"); - } - if (source->moments->Myy < 0) { - fprintf (stderr, "error: neg second moment??\n"); - } - - // XXX TEST: - // pmSourceMoments_Old (source, radius); + source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, numPixels); + return(true); } - - -bool pmSourceMoments_Old(pmSource *source, psF32 radius) -{ - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - // a note about coordinates: coordinates of objects throughout psphot refer to the primary - // image coordinates. the source->pixels image has an offset relative to its parent of - // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in - // this subimage. we subtract off the peak coordinates, adjusted to this subimage, to have - // minimal round-off error in the sums - - psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage - psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col - xPeak; - psF32 yDiff = row - yPeak; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) continue; - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - X2 += xDiff * xWght; - XY += xDiff * yWght; - Y2 += yDiff * yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); - return (false); - } - - psTrace ("psModules.objects", 4, "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - return (false); - } - -# if (PS_TRACE_ON) - float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0); - float Sxy = XY / Sum; - float Syy = PS_MAX(Y2/Sum - PS_SQR(y), 0); - - psTrace ("psModules.objects", 3, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, x, y, Sxx, Sxy, Syy); -# endif - - return(true); -} - Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 25738) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 25754) @@ -84,4 +84,5 @@ // we must have a valid model + // XXX allow aperture magnitudes for sources without a model model = pmSourceGetModel (&isPSF, source); if (model == NULL) { @@ -90,4 +91,5 @@ } + // XXX handle negative flux, low-significance if (model->dparams->data.F32[PM_PAR_I0] > 0) { SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0]; @@ -101,37 +103,32 @@ // measure PSF model photometry - if (psf->FluxScale) { + // XXX TEST: do not use flux scale + if (0 && psf->FluxScale) { // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y); - if (!isfinite(fluxScale)) { - // XXX objects on the edge can be slightly outside -- if we get an - // error, use the full fit. - psErrorClear(); - status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); - } else { - if (isfinite(fluxScale) && (fluxScale > 0.0)) { - source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]); - } else { - source->psfMag = NAN; - } - } + psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); + psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); + source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]); } else { status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); } - // if we have a collection of model fits, one of them is a pointer to modelEXT? - // XXX not guaranteed + // if we have a collection of model fits, check if one of them is a pointer to modelEXT if (source->modelFits) { - for (int i = 0; i < source->modelFits->n; i++) { - pmModel *model = source->modelFits->data[i]; - status = pmSourcePhotometryModel (&model->mag, model); - } - if (source->modelEXT) { - source->extMag = source->modelEXT->mag; - } + bool foundEXT = false; + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + status = pmSourcePhotometryModel (&model->mag, model); + if (model == source->modelEXT) foundEXT = true; + } + if (foundEXT) { + source->extMag = source->modelEXT->mag; + } else { + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + } } else { - if (source->modelEXT) { - status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); - } + if (source->modelEXT) { + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + } } @@ -160,11 +157,6 @@ } - // replace source flux - // XXX need to be certain which model and size of mask for prior subtraction - // XXX full model or just analytical? - // XXX use pmSourceAdd instead? - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); - } + // if we measure aperture magnitudes, the source must not currently be subtracted! + psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!"); // if we are measuring aperture photometry and applying the growth correction, @@ -201,9 +193,10 @@ if (isfinite (source->apMag) && isPSF && psf) { if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) { - source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit); + source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius); } if (mode & PM_SOURCE_PHOT_APCORR) { + // XXX this should be removed -- we no longer fit for the 'sky bias' rflux = pow (10.0, 0.4*source->psfMag); - source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux; + source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux; } } @@ -211,10 +204,4 @@ psFree(flux); psFree(mask); - } - - // if source was originally subtracted, re-subtract object, leave local sky - // XXX replace with pmSourceSub... - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); } @@ -763,55 +750,2 @@ return flux; } - -// XXX this is test code to verify the shift is doing the right thing (seems to be) -# if (0) -// measure centroid of unshifted gaussian (should be 16.0,16.0) - { - psImage *image = source->pixels; - float xo = 0.0; - float yo = 0.0; - float xo2 = 0.0; - float yo2 = 0.0; - float no = 0.0; - for (int j = 0; j < image->numRows; j++) - { - for (int i = 0; i < image->numCols; i++) { - xo += i*image->data.F32[j][i]; - yo += j*image->data.F32[j][i]; - xo2 += i*i*image->data.F32[j][i]; - yo2 += j*j*image->data.F32[j][i]; - no += image->data.F32[j][i]; - } - } - xo /= no; - yo /= no; - xo2 = sqrt (xo2/no - xo*xo); - yo2 = sqrt (yo2/no - yo*yo); - fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no); - } - -// measure centroid of unshifted gaussian (should be 16.0,16.0) - { - psImage *image = flux; - float xo = 0.0; - float yo = 0.0; - float xo2 = 0.0; - float yo2 = 0.0; - float no = 0.0; - for (int j = 0; j < image->numRows; j++) - { - for (int i = 0; i < image->numCols; i++) { - xo += i*image->data.F32[j][i]; - yo += j*image->data.F32[j][i]; - xo2 += i*i*image->data.F32[j][i]; - yo2 += j*j*image->data.F32[j][i]; - no += image->data.F32[j][i]; - } - } - xo /= no; - yo /= no; - xo2 = sqrt (xo2/no - xo*xo); - yo2 = sqrt (yo2/no - yo*yo); - fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no); - } -# endif Index: trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.c (revision 25738) +++ trunk/psModules/src/objects/pmSourceVisual.c (revision 25754) @@ -5,4 +5,7 @@ #include #include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmSource.h" #include "pmSourceVisual.h" @@ -15,6 +18,6 @@ static int kapa1 = -1; +static int kapa2 = -1; static bool plotPSF = true; -// static int kapa2 = -1; // static int kapa3 = -1; @@ -27,12 +30,9 @@ bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi); - -bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) { - - KapaSection section; // put the positive profile in one and the residuals in another? +bool pmSourceVisualPlotPSFMetric (pmPSFtry *psfTry) { Graphdata graphdata; - if (!pmVisualIsVisual() || !plotPSF) return true; + if (!pmVisualIsVisual()) return true; if (kapa1 == -1) { @@ -45,41 +45,249 @@ } + KapaClearSections (kapa1); + KapaInitGraph (&graphdata); + + psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32); + + graphdata.xmin = +32.0; + graphdata.xmax = -32.0; + graphdata.ymin = +32.0; + graphdata.ymax = -32.0; + + // construct the plot vectors + int n = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + x->data.F32[n] = psfTry->fitMag->data.F32[i]; + y->data.F32[n] = psfTry->metric->data.F32[i]; + dy->data.F32[n] = psfTry->metricErr->data.F32[i]; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + n++; + } + x->n = y->n = dy->n = n; + + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + + // better choice for range? + // graphdata.xmin = -17.0; + // graphdata.xmax = -9.0; + graphdata.ymin = -0.51; + graphdata.ymax = +0.51; + + KapaSetLimits (kapa1, &graphdata); + + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + graphdata.etype |= 0x01; + + KapaPrepPlot (kapa1, n, &graphdata); + KapaPlotVector (kapa1, n, x->data.F32, "x"); + KapaPlotVector (kapa1, n, y->data.F32, "y"); + KapaPlotVector (kapa1, n, dy->data.F32, "dym"); + KapaPlotVector (kapa1, n, dy->data.F32, "dyp"); + + psFree (x); + psFree (y); + psFree (dy); + + pmVisualAskUser(NULL); + return true; +} + +// to see the structure of the psf model, place the sources in a fake image 1/10th the size +// at their appropriate relative location. later sources stomp on earlier sources +bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) { + + if (!pmVisualIsVisual()) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + // create images 1/10 scale: + psImage *image = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImage *model = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImage *resid = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImageInit (image, 0.0); + psImageInit (model, 0.0); + psImageInit (resid, 0.0); + + for (int i = sources->n - 1; i >= 0; i--) { + pmSource *source = sources->data[i]; + if (!source) continue; + if (!source->pixels) continue; + + pmSourceCacheModel (source, maskVal); + if (!source->modelFlux) continue; + + pmModel *srcModel = pmSourceGetModel (NULL, source); + if (!model) continue; + + float norm = srcModel->params->data.F32[PM_PAR_I0]; + + int Xo = 0.1*srcModel->params->data.F32[PM_PAR_XPOS]; + int Yo = 0.1*srcModel->params->data.F32[PM_PAR_YPOS]; + + // insert source pixels in the image at 1/10th offset + for (int iy = 0; iy < source->pixels->numRows; iy++) { + int jy = iy + Yo; + if (jy >= image->numRows) continue; + for (int ix = 0; ix < source->pixels->numCols; ix++) { + int jx = ix + Xo; + if (jx >= image->numCols) continue; + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue; + if (source->modelFlux->data.F32[iy][ix] < 0.001) continue; + image->data.F32[jy][jx] = source->pixels->data.F32[iy][ix]; + model->data.F32[jy][jx] = source->modelFlux->data.F32[iy][ix]; + resid->data.F32[jy][jx] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix]; + } + } + } + + // KapaClearSections (kapa2); + pmVisualScaleImage (kapa2, image, "image", 0, true); + pmVisualScaleImage (kapa2, model, "model", 1, true); + pmVisualScaleImage (kapa2, resid, "resid", 2, true); + +# ifdef DEBUG + { + psFits *fits = psFitsOpen ("image.fits", "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("model.fits", "w"); + psFitsWriteImage (fits, NULL, model, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("resid.fits", "w"); + psFitsWriteImage (fits, NULL, resid, 0, NULL); + psFitsClose (fits); + } +# endif + + psFree (image); + psFree (model); + psFree (resid); + + pmVisualAskUser(NULL); + return true; +} + +bool pmSourceVisualShowModelFit (pmSource *source) { + + if (!pmVisualIsVisual()) return true; + if (!source->pixels) return false; + if (!source->modelFlux) return false; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + // KapaClearSections (kapa2); + pmVisualScaleImage (kapa2, source->pixels, "source", 0, false); + pmVisualScaleImage (kapa2, source->modelFlux, "model", 1, false); + + pmModel *model = pmSourceGetModel (NULL, source); + float norm = model->params->data.F32[PM_PAR_I0]; + + psImage *resid = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + resid->data.F32[iy][ix] = NAN; + continue; + } + resid->data.F32[iy][ix] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix]; + } + } + pmVisualScaleImage (kapa2, resid, "resid", 2, false); + + psFree (resid); + + pmVisualAskUser(NULL); + return true; +} + +bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) { + + KapaSection section; // put the positive profile in one and the residuals in another? + + Graphdata graphdata; + + if (!pmVisualIsVisual() || !plotPSF) return true; + + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + KapaClearPlots (kapa1); KapaInitGraph (&graphdata); - float min = +1e32; - float max = -1e32; - float Min = +1e32; - float Max = -1e32; + float Xmin = +1e32; + float Xmax = -1e32; + float Ymin = +1e32; + float Ymax = -1e32; + float Fmin = +1e32; + float Fmax = -1e32; psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32); psVector *model = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *xm = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *ym = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *Fm = psVectorAlloc (x->n, PS_TYPE_F32); + + int n = 0; for (int i = 0; i < x->n; i++) { model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]); resid->data.F32[i] = param->data.F32[i] - model->data.F32[i]; if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - min = PS_MIN (min, resid->data.F32[i]); - max = PS_MAX (max, resid->data.F32[i]); - Min = PS_MIN (min, param->data.F32[i]); - Max = PS_MAX (max, param->data.F32[i]); - } - - psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32); - for (int i = 0; i < x->n; i++) { - xn->data.F32[i] = x->data.F32[i] / 5000.0; - yn->data.F32[i] = y->data.F32[i] / 5000.0; - zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min); - Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min); - Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min); - } + Xmin = PS_MIN (Xmin, x->data.F32[i]); + Xmax = PS_MAX (Xmax, x->data.F32[i]); + Ymin = PS_MIN (Ymin, y->data.F32[i]); + Ymax = PS_MAX (Ymax, y->data.F32[i]); + Fmin = PS_MIN (Fmin, param->data.F32[i]); + Fmax = PS_MAX (Fmax, param->data.F32[i]); + xm->data.F32[n] = x->data.F32[i]; + ym->data.F32[n] = y->data.F32[i]; + Fm->data.F32[n] = param->data.F32[i]; + n++; + } + xm->n = ym->n = Fm->n = n; // view 1 on resid - section.dx = 0.5; - section.dy = 0.33; + section.dx = 1.0; + section.dy = 0.5; section.x = 0.0; section.y = 0.0; @@ -88,67 +296,98 @@ KapaSetSection (kapa1, §ion); psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.xmax = Xmax; + graphdata.ymin = Fmin; + graphdata.ymax = Fmax; + + { + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + } + + KapaSetLimits (kapa1, &graphdata); + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "X (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM); + + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapa1, x->n, &graphdata); + KapaPlotVector (kapa1, x->n, x->data.F32, "x"); + KapaPlotVector (kapa1, x->n, param->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, xm->n, &graphdata); + KapaPlotVector (kapa1, xm->n, xm->data.F32, "x"); + KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, x->n, &graphdata); + KapaPlotVector (kapa1, x->n, x->data.F32, "x"); + KapaPlotVector (kapa1, x->n, model->data.F32, "y"); // view 2 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.5; - section.y = 0.0; + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.5; section.name = NULL; psStringAppend (§ion.name, "a2"); KapaSetSection (kapa1, §ion); psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); - - // view 3 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.0; - section.y = 0.33; - section.name = NULL; - psStringAppend (§ion.name, "a3"); - KapaSetSection (kapa1, §ion); - psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); - - // view 4 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.5; - section.y = 0.33; - section.name = NULL; - psStringAppend (§ion.name, "a4"); - KapaSetSection (kapa1, §ion); - psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); - - // view 5 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.0; - section.y = 0.66; - section.name = NULL; - psStringAppend (§ion.name, "a5"); - KapaSetSection (kapa1, §ion); - psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); - - // view 6 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.5; - section.y = 0.66; - section.name = NULL; - psStringAppend (§ion.name, "a6"); - KapaSetSection (kapa1, §ion); - psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Ymin; + graphdata.xmax = Ymax; + graphdata.ymin = Fmin; + graphdata.ymax = Fmax; + { + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + } + + KapaSetLimits (kapa1, &graphdata); + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "Y (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM); + + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapa1, y->n, &graphdata); + KapaPlotVector (kapa1, y->n, y->data.F32, "x"); + KapaPlotVector (kapa1, y->n, param->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, xm->n, &graphdata); + KapaPlotVector (kapa1, xm->n, ym->data.F32, "x"); + KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, y->n, &graphdata); + KapaPlotVector (kapa1, y->n, y->data.F32, "x"); + KapaPlotVector (kapa1, y->n, model->data.F32, "y"); psFree (resid); - - psFree (xn); - psFree (yn); - psFree (zn); - psFree (Zn); + psFree (model); // pause and wait for user input: @@ -159,6 +398,8 @@ } -// send in normalized points +// Somewhat broken 3D plotting function (was used by pmSourceVisualPSFModelResid, but not anymore) bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi) { + + return true; psVector *xv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32); @@ -192,9 +433,4 @@ KapaSetLimits (myKapa, graphdata); - // KapaSetFont (myKapa, "helvetica", 14); - // KapaBox (myKapa, graphdata); - // KapaSendLabel (myKapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM); - // KapaSendLabel (myKapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM); - graphdata->color = KapaColorByName ("black"); graphdata->ptype = 100; Index: trunk/psModules/src/objects/pmSourceVisual.h =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.h (revision 25738) +++ trunk/psModules/src/objects/pmSourceVisual.h (revision 25754) @@ -18,4 +18,7 @@ bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask); +bool pmSourceVisualPlotPSFMetric (pmPSFtry *try); +bool pmSourceVisualShowModelFit (pmSource *source); +bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal); /// @} Index: trunk/psModules/src/objects/pmTrend2D.c =================================================================== --- trunk/psModules/src/objects/pmTrend2D.c (revision 25738) +++ trunk/psModules/src/objects/pmTrend2D.c (revision 25754) @@ -298,2 +298,21 @@ return PM_TREND_NONE; } + +bool pmTrend2DPrintMap (pmTrend2D *trend) { + + if (!trend->map) return false; + if (!trend->map->map) return false; + + for (int j = 0; j < trend->map->map->numRows; j++) { + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->map->data.F32[j][i]); + } + fprintf (stderr, "\t\t\t"); + for (int i = 0; i < trend->map->map->numCols; i++) { + fprintf (stderr, "%5.2f ", trend->map->error->data.F32[j][i]); + } + fprintf (stderr, "\n"); + } + return true; +} + Index: trunk/psModules/src/objects/pmTrend2D.h =================================================================== --- trunk/psModules/src/objects/pmTrend2D.h (revision 25738) +++ trunk/psModules/src/objects/pmTrend2D.h (revision 25754) @@ -97,4 +97,6 @@ pmTrend2DMode pmTrend2DModeFromString(psString name); +bool pmTrend2DPrintMap (pmTrend2D *trend); + /// @} # endif Index: trunk/psModules/test/objects/tap_pmGrowthCurve.c =================================================================== --- trunk/psModules/test/objects/tap_pmGrowthCurve.c (revision 25738) +++ trunk/psModules/test/objects/tap_pmGrowthCurve.c (revision 25754) @@ -131,30 +131,30 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0003, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0020, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; + source->apRadius = 3.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 2.0; + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0075, 0.0001, "growth offset is is %f", refMag - source->apMag); @@ -234,29 +234,29 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0004, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0026, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; + source->apRadius = 3.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 2.0; + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0103, 0.0001, "growth offset is is %f", refMag - source->apMag); @@ -336,29 +336,29 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0006, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0038, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 2.0; + source->apRadius = 3.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0164, 0.0001, "growth offset is is %f", refMag - source->apMag); Index: trunk/psModules/test/objects/tap_pmModel.c =================================================================== --- trunk/psModules/test/objects/tap_pmModel.c (revision 25738) +++ trunk/psModules/test/objects/tap_pmModel.c (revision 25754) @@ -77,5 +77,5 @@ ok(model->nDOF == 0, "pmModelAlloc() set pmModel->nDOF correctly"); ok(model->nIter == 0, "pmModelAlloc() set pmModel->nIter correctly"); - ok(model->radiusFit == 0, "pmModelAlloc() set pmModel->radiusFit correctly"); + ok(model->fitRadius == 0, "pmModelAlloc() set pmModel->fitRadius correctly"); ok(model->flags == PM_MODEL_STATUS_NONE, "pmModelAlloc() set pmModel->flags correctly"); ok(model->residuals == NULL, "pmModelAlloc() set pmModel->residuals correctly"); @@ -132,5 +132,5 @@ modelSrc->nIter = 3; modelSrc->flags = PM_MODEL_STATUS_NONE; - modelSrc->radiusFit = 4; + modelSrc->fitRadius = 4; pmModel *modelDst = pmModelCopy(modelSrc); ok(modelDst != NULL && psMemCheckModel(modelDst), "pmModelCopy() returned a non-NULL pmModel"); @@ -139,5 +139,5 @@ ok(modelDst->nIter == 3, "pmModelCopy() set the pmModel->nIter member correctly"); ok(modelDst->flags == PM_MODEL_STATUS_NONE, "pmModelCopy() set the pmModel->flags member correctly"); - ok(modelDst->radiusFit == 4, "pmModelCopy() set the pmModel->radiusFit member correctly"); + ok(modelDst->fitRadius == 4, "pmModelCopy() set the pmModel->fitRadius member correctly"); psFree(modelSrc); Index: trunk/psModules/test/objects/tap_pmModelUtils.c =================================================================== --- trunk/psModules/test/objects/tap_pmModelUtils.c (revision 25738) +++ trunk/psModules/test/objects/tap_pmModelUtils.c (revision 25754) @@ -81,5 +81,5 @@ ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly"); ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly"); - ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly"); + ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly"); ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly"); ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly"); @@ -140,5 +140,5 @@ ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly"); ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly"); - ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly"); + ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly"); ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly"); ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly"); Index: trunk/psModules/test/objects/tap_pmSourcePhotometry.c =================================================================== --- trunk/psModules/test/objects/tap_pmSourcePhotometry.c (revision 25738) +++ trunk/psModules/test/objects/tap_pmSourcePhotometry.c (revision 25754) @@ -96,5 +96,5 @@ source->modelPSF = pmModelFromPSF (modelRef, psf); source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1; - source->modelPSF->radiusFit = radius; + source->apRadius = radius; // measure photometry for centered source (fractional pix : 0.5,0.5)