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