Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 9730) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 9770) @@ -1,37 +1,52 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - params->data.F32[PM_PAR_SKY] = So; - params->data.F32[PM_PAR_I0] = Zo; - params->data.F32[PM_PAR_XPOS] = Xo; - params->data.F32[PM_PAR_YPOS] = Yo; - params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX; - params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY; - params->data.F32[PM_PAR_SXY] = Sxy; -*****************************************************************************/ - -psF32 pmModelFunc_GAUSS(psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 X = x->data.F32[0] - params->data.F32[PM_PAR_XPOS]; - psF32 Y = x->data.F32[1] - params->data.F32[PM_PAR_YPOS]; - psF32 px = params->data.F32[PM_PAR_SXX]*X; - psF32 py = params->data.F32[PM_PAR_SYY]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[PM_PAR_SXY]*X*Y; + * 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 + * models use a psVector to represent the set of parameters, with the sequence used to specify + * the meaning of the parameter. The meaning of the parameters may thus vary depending on the + * specifics of the model. All models which are used a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_GAUSS +# define PM_MODEL_FLUX pmModelFlux_GAUSS +# define PM_MODEL_GUESS pmModelGuess_GAUSS +# define PM_MODEL_LIMITS pmModelLimits_GAUSS +# define PM_MODEL_RADIUS pmModelRadius_GAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_GAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS + +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + // XXX this is fitting sqrt(2)/sigma_x, sqrt(2)/sigma_y + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; psF32 r = exp(-z); - psF32 q = params->data.F32[PM_PAR_I0]*r; - psF32 f = q + params->data.F32[PM_PAR_SKY]; + psF32 q = PAR[PM_PAR_I0]*r; + psF32 f = q + PAR[PM_PAR_SKY]; if (deriv != NULL) { - deriv->data.F32[PM_PAR_SKY] = +1.0; - deriv->data.F32[PM_PAR_I0] = +r; - deriv->data.F32[PM_PAR_XPOS] = q*(2*px*params->data.F32[PM_PAR_SXX] + params->data.F32[PM_PAR_SXY]*Y); - deriv->data.F32[PM_PAR_YPOS] = q*(2*py*params->data.F32[PM_PAR_SYY] + params->data.F32[PM_PAR_SXY]*X); - deriv->data.F32[PM_PAR_SXX] = -2.0*q*px*X; - deriv->data.F32[PM_PAR_SYY] = -2.0*q*py*Y; - deriv->data.F32[PM_PAR_SXY] = -q*X*Y; + psF32 *dPAR = deriv->data.F32; + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; } return(f); @@ -39,5 +54,5 @@ // define the parameter limits -bool pmModelLimits_GAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) { @@ -58,6 +73,6 @@ 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.01; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.01; + 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; @@ -66,6 +81,6 @@ 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] = 2.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 2.0; + 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; @@ -74,27 +89,34 @@ // make an initial guess for parameters -bool pmModelGuess_GAUSS (pmModel *model, pmSource *source) +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { pmMoments *moments = source->moments; - psF32 *params = model->params->data.F32; - - params[PM_PAR_SKY] = moments->Sky; - params[PM_PAR_I0] = moments->Peak - moments->Sky; - params[PM_PAR_XPOS] = moments->x; - params[PM_PAR_YPOS] = moments->y; - params[PM_PAR_SXX] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx); - params[PM_PAR_SYY] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy); - params[PM_PAR_SXY] = 0.0; + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = moments->x; + PAR[PM_PAR_YPOS] = moments->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = 0.0; return(true); } -psF64 pmModelFlux_GAUSS(const psVector *params) -{ - psF64 A1 = PS_SQR(params->data.F32[PM_PAR_SXX]); - psF64 A2 = PS_SQR(params->data.F32[PM_PAR_SYY]); - psF64 A3 = PS_SQR(params->data.F32[PM_PAR_SXY]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); +psF64 PM_MODEL_FLUX (const psVector *params) +{ + + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; psF64 Flux = params->data.F32[PM_PAR_I0] * Area; @@ -105,6 +127,7 @@ // return the radius which yields the requested flux // this function is never allowed to return <= 0 -psF64 pmModelRadius_GAUSS (const psVector *params, psF64 flux) -{ +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psEllipseShape shape; if (flux <= 0) @@ -115,11 +138,17 @@ return (1.0); - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[PM_PAR_SXX], 1.0 / params->data.F32[PM_PAR_SYY]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); return (radius); } // construct the PSF model from the FLT model and the psf -bool pmModelFromPSF_GAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { @@ -127,19 +156,28 @@ psF32 *in = modelFLT->params->data.F32; - out[PM_PAR_SKY] = in[PM_PAR_SKY]; - out[PM_PAR_I0] = in[PM_PAR_I0]; - out[PM_PAR_XPOS] = in[PM_PAR_XPOS]; - out[PM_PAR_YPOS] = in[PM_PAR_YPOS]; - - assert (PM_PAR_YPOS + 1 == 4); // so starting at 4 is correct - for (int i = 4; i < 7; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]); + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + return(true); } // check the status of the fitted model -bool pmModelFitStatus_GAUSS (pmModel *model) +// this test is invalid if the parameters are derived +// from the PSF model +bool PM_MODEL_FIT_STATUS (pmModel *model) { @@ -164,2 +202,10 @@ return false; } + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 9730) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 9770) @@ -1,46 +1,58 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - /****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2.0) / SigmaX; - params->data.F32[5] = sqrt(2.0) / SigmaY; - params->data.F32[6] = Sxy; -*****************************************************************************/ - -psF32 pmModelFunc_PGAUSS(psVector *deriv, - const psVector *params, - const psVector *x) + * 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 + * models use a psVector to represent the set of parameters, with the sequence used to specify + * the meaning of the parameter. The meaning of the parameters may thus vary depending on the + * specifics of the model. All models which are used a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + *****************************************************************************/ + +# define PM_MODEL_FUNC pmModelFunc_PGAUSS +# define PM_MODEL_FLUX pmModelFlux_PGAUSS +# define PM_MODEL_GUESS pmModelGuess_PGAUSS +# define PM_MODEL_LIMITS pmModelLimits_PGAUSS +# define PM_MODEL_RADIUS pmModelRadius_PGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_PGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *pixcoord) { psF32 *PAR = params->data.F32; - 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 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; psF32 t = 1 + z + z*z/2.0; psF32 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ - psF32 f = PAR[1]*r + PAR[0]; + psF32 f = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY]; if (deriv != NULL) { psF32 *dPAR = deriv->data.F32; - psF32 q = PAR[1]*r*r*t; - dPAR[0] = +1.0; - dPAR[1] = +r; - dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - dPAR[4] = -2.0*q*px*X; - dPAR[5] = -2.0*q*py*Y; - dPAR[6] = -q*X*Y; + psF32 q = PAR[PM_PAR_I0]*r*r*t; + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; } return(f); } -bool pmModelLimits_PGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) +bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) { @@ -49,27 +61,27 @@ *params_max = psVectorAlloc (7, PS_TYPE_F32); - beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 3e6; - 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; - - 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_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e8; - 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] = +5.0; + beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; + beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; + beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; + beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; + beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; + + params_min[0][0].data.F32[PM_PAR_SKY] = -1000; + params_min[0][0].data.F32[PM_PAR_I0] = 0; + params_min[0][0].data.F32[PM_PAR_XPOS] = -100; + params_min[0][0].data.F32[PM_PAR_YPOS] = -100; + params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; + params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; + params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; + + params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; + params_max[0][0].data.F32[PM_PAR_I0] = 1e8; + params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; + params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; + params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; return (TRUE); @@ -77,32 +89,35 @@ // make an initial guess for parameters -bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source) +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { pmMoments *moments = source->moments; - psF32 *params = model->params->data.F32; - - params[0] = moments->Sky; - params[1] = moments->Peak - moments->Sky; - params[2] = moments->x; - params[3] = moments->y; - params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx); - params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy); - params[6] = 0.0; + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + PAR[PM_PAR_XPOS] = moments->x; + PAR[PM_PAR_YPOS] = moments->y; + PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy); + PAR[PM_PAR_SXY] = 0.0; return(true); } -psF64 pmModelFlux_PGAUSS(const psVector *params) +psF64 PM_MODEL_FLUX(const psVector *params) { float norm, z; + psEllipseShape shape; 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); + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; // the area needs to be multiplied by the integral of f(z) @@ -122,5 +137,5 @@ norm *= DZ / 3.0; - psF64 Flux = PAR[1] * Area * norm; + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; return(Flux); @@ -129,22 +144,32 @@ // define this function so it never returns Inf or NaN // return the radius which yields the requested flux -psF64 pmModelRadius_PGAUSS (const psVector *params, psF64 flux) +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { if (flux <= 0) return (1.0); - if (params->data.F32[1] <= 0) + if (params->data.F32[PM_PAR_I0] <= 0) return (1.0); - if (flux >= params->data.F32[1]) + if (flux >= params->data.F32[PM_PAR_I0]) 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"); - } + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); + shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0); + + // this estimates the radius assuming f(z) is roughly exp(-z) + psEllipseAxes axes = psEllipseShapeToAxes (shape); + psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux)); + + if (isnan(radius)) + psAbort ("psphot.model", "error in code: never return invalid radius"); + if (radius < 0) + psAbort ("psphot.model", "error in code: never return invalid radius"); + return (radius); } -bool pmModelFromPSF_PGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { @@ -152,17 +177,24 @@ 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 < 7; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = psPolynomial2DEval(poly, out[2], out[3]); + // we require these two parameters to exist + assert (psf->params_NEW->n > PM_PAR_YPOS); + assert (psf->params_NEW->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) { + out[i] = in[i]; + } else { + psPolynomial2D *poly = psf->params_NEW->data[i]; + out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } } + + // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here + out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + return(true); } -bool pmModelFitStatus_PGAUSS (pmModel *model) +bool PM_MODEL_FIT_STATUS (pmModel *model) { @@ -174,12 +206,12 @@ dP = 0; - dP += PS_SQR(dPAR[4] / PAR[4]); - dP += PS_SQR(dPAR[5] / PAR[5]); + 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[1] > 0); - status &= ((dPAR[1]/PAR[1]) < 0.5); + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); if (status) @@ -187,2 +219,10 @@ return false; } + +# undef PM_MODEL_FUNC +# undef PM_MODEL_FLUX +# undef PM_MODEL_GUESS +# undef PM_MODEL_LIMITS +# undef PM_MODEL_RADIUS +# undef PM_MODEL_FROM_PSF +# undef PM_MODEL_FIT_STATUS