Index: /trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 9769) +++ /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 9769) +++ /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 Index: /trunk/psModules/src/objects/pmModelGroup.c =================================================================== --- /trunk/psModules/src/objects/pmModelGroup.c (revision 9769) +++ /trunk/psModules/src/objects/pmModelGroup.c (revision 9770) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-15 09:49:01 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-28 20:23:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,14 +35,14 @@ double sqrt (double x); -#include "models/pmModel_GAUSS.c" -#include "models/pmModel_PGAUSS.c" -#include "models/pmModel_QGAUSS.c" -#include "models/pmModel_SGAUSS.c" +# include "models/pmModel_GAUSS.c" +// # include "models/pmModel_PGAUSS.c" +// # include "models/pmModel_QGAUSS.c" +// # include "models/pmModel_SGAUSS.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_SGAUSS", 9, pmModelFunc_SGAUSS, pmModelFlux_SGAUSS, pmModelRadius_SGAUSS, pmModelLimits_SGAUSS, pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS}, + // {"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_SGAUSS", 9, pmModelFunc_SGAUSS, pmModelFlux_SGAUSS, pmModelRadius_SGAUSS, pmModelLimits_SGAUSS, pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS}, }; @@ -208,6 +208,7 @@ { psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source, NULL); + PS_ASSERT_PTR_NON_NULL(source->moments, NULL); + PS_ASSERT_PTR_NON_NULL(source->peak, NULL); pmModel *model = pmModelAlloc(modelType); Index: /trunk/psModules/src/objects/pmPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSF.c (revision 9769) +++ /trunk/psModules/src/objects/pmPSF.c (revision 9770) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-24 22:55:05 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-28 20:23:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -69,14 +69,15 @@ psFree (psf->ApTrend); psFree (psf->growth); - psFree (psf->params); + psFree (psf->params_NEW); return; } - - /***************************************************************************** -pmPSFAlloc (type): allocate a pmPSF. - NOTE: a PSF always has 4 parameters fewer than the equivalent model. - This is because those 4 parameters are + pmPSFAlloc (type): allocate a pmPSF. + + NOTE: PSF model parameters which are not modeled on an image are set to NULL in psf->params. + + These are normally: + X-center Y-center @@ -119,17 +120,24 @@ return(NULL); } - - psf->params = psArrayAlloc(Nparams - 4); - - // the order of the PSF parameter (X,Y) fits is determined by the - // psfTrendMask polynomial (user-specified as in the recipe). the - // masks of psfTrendMask are applied to each parameter. - // if psfTrendMask is NULL, these polynomials are not allocated. - // in this case, the user must set them by hand (as in pmPSFfromMD) - // XXX should we drop the hard-wired '4' above and use NULL to identify - // the parameters which are not fitted. these could be selected by - // testing for the value of PM_PAR_XPOS, etc. + psf->params_NEW = psArrayAlloc(Nparams); + + // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial + // (user-specified as in the recipe). the masks of psfTrendMask are applied to each + // parameter. if psfTrendMask is NULL, these polynomials are not allocated. in this case, + // the user must set them by hand (as in pmPSFfromMD). the parameters which are not fitted + // are left as NULL. these are selected by testing for them by the named values ( + // PM_PAR_XPOS, etc) + if (psfTrendMask) { - for (int i = 0; i < psf->params->n; i++) { + for (int i = 0; i < psf->params_NEW->n; i++) { + if (i == PM_PAR_SKY) + continue; + if (i == PM_PAR_I0) + continue; + if (i == PM_PAR_XPOS) + continue; + if (i == PM_PAR_YPOS) + continue; + psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY); for (int nx = 0; nx < param->nX + 1; nx++) { @@ -138,5 +146,5 @@ } } - psf->params->data[i] = param; + psf->params_NEW->data[i] = param; } } @@ -166,4 +174,5 @@ psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64); + // construct the x,y terms for (int i = 0; i < models->n; i++) { pmModel *model = models->data[i]; @@ -171,8 +180,7 @@ continue; - // XXX EAM : this is fragile: x and y must be stored consistently at 2,3 - // note that the data is provided in the F64 array - x->data.F64[i] = model->params->data.F32[2]; - y->data.F64[i] = model->params->data.F32[3]; + // use F64 for polynomial fitting + x->data.F64[i] = model->params->data.F32[PM_PAR_XPOS]; + y->data.F64[i] = model->params->data.F32[PM_PAR_YPOS]; } @@ -183,23 +191,60 @@ psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - for (int i = 0; i < psf->params->n; i++) { + // skip the unfitted parameters (X, Y, Io, Sky) + for (int i = 0; i < psf->params_NEW->n; i++) { + if (i == PM_PAR_SKY) + continue; + if (i == PM_PAR_I0) + continue; + if (i == PM_PAR_XPOS) + continue; + if (i == PM_PAR_YPOS) + continue; + + // select the per-object fitted data for this parameter for (int j = 0; j < models->n; j++) { pmModel *model = models->data[j]; if (model == NULL) continue; - z->data.F64[j] = model->params->data.F32[i + 4]; - dz->data.F64[j] = 1; - // XXX EAM : need to use actual errors? - // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4 + z->data.F64[j] = model->params->data.F32[i]; + dz->data.F64[j] = 1; // use the model-fitted error? or S/N? + + // for SXY, we actually fit SXY / (SXX^-2 + SYY^-2) + if (i == PM_PAR_SXY) { + z->data.F64[j] = pmPSF_SXYfromModel (model->params->data.F32); + } } - - // XXX EAM : this is the expected API (cycle 7? cycle 8?) - psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y); - - // XXX EAM : drop this when above is implemented... - // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz); - + psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, mask, 0xff, z, dz, x, y); // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i); - // psPolynomial2DDump (psf->params->data[i]); + + // XXX Test output + psPolynomial2D *poly = psf->params_NEW->data[i]; + fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev); + fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][0], poly->coeff[1][0], poly->coeff[0][1]); + fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][2], poly->coeff[1][1], poly->coeff[2][0]); + } + + // XXX test dump of star parameters vs position (compare with fitted values) + { + FILE *f = fopen ("params.dat", "w"); + + for (int j = 0; j < models->n; j++) + { + pmModel *model = models->data[j]; + if (model == NULL) + continue; + + pmModel *modelPSF = pmModelFromPSF (model, psf); + + fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + + for (int i = 0; i < psf->params_NEW->n; i++) { + if (psf->params_NEW->data[i] == NULL) + continue; + fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]); + } + fprintf (f, "%f %d\n", model->chisq, model->nIter); + } + fclose (f); } @@ -211,6 +256,4 @@ return (true); } - - /***************************************************************************** @@ -250,4 +293,33 @@ } return true; +} + +// the PSF models the \sigma_{xy} variation of the elliptical contour as a function of position in the image with a +// polynomial. an individual object has a contour of the form (x^2/2sx^2) + (y^2/2sy^2) + sxy*x*y +// these are the values of the model->params. the psf->params term for sxy is actually fitted +// to sxy/(sxx^-2 + syy^-2)^2 + +// input: model->param, output: psf->param[PM_PAR_SXY] +double pmPSF_SXYfromModel (psF32 *modelPar) +{ + + double SXX = modelPar[PM_PAR_SXX]; + double SYY = modelPar[PM_PAR_SYY]; + double SXY = modelPar[PM_PAR_SXY]; + + double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); + return (par); +} + +// input: fitted psf->param, output: model->param[PM_PAR_SXY] +double pmPSF_SXYtoModel (psF32 *fittedPar) +{ + + double SXX = fittedPar[PM_PAR_SXX]; + double SYY = fittedPar[PM_PAR_SYY]; + double fit = fittedPar[PM_PAR_SXY]; + + double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); + return SXY; } Index: /trunk/psModules/src/objects/pmPSF.h =================================================================== --- /trunk/psModules/src/objects/pmPSF.h (revision 9769) +++ /trunk/psModules/src/objects/pmPSF.h (revision 9770) @@ -42,5 +42,8 @@ { pmModelType type; ///< PSF Model in use - psArray *params; ///< Model parameters (psPolynomial2D) + psArray *params_NEW; ///< Model parameters (psPolynomial2D) + // XXXXX I am changing params: we will allocate elements for the + // unfitted elements (So, Io, Xo, Yo) and leave them as NULL + // I am using a new name to catch all refs to params with gcc psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) psPolynomial4D *ApTrend; ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)) Index: /trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- /trunk/psModules/src/objects/pmPSF_IO.c (revision 9769) +++ /trunk/psModules/src/objects/pmPSF_IO.c (revision 9770) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-14 00:56:13 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-28 20:23:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -50,6 +50,8 @@ psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); - for (int i = 0; i < nPar - 4; i++) { - psPolynomial2D *poly = psf->params->data[i]; + for (int i = 0; i < nPar; i++) { + psPolynomial2D *poly = psf->params_NEW->data[i]; + if (poly == NULL) + continue; psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i); } @@ -87,10 +89,14 @@ psAbort ("read PSF" , "mismatch model par count"); - for (int i = 0; i < nPar - 4; i++) { + // un-fitted terms, not in the Metadata, are left NULL + // XXX add a double-check of the expected number? + for (int i = 0; i < nPar; i++) { sprintf (keyword, "PSF_PAR%02d", i); psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); + if (!status) + continue; psPolynomial2D *poly = psPolynomial2DfromMetadata (folder); - psFree (psf->params->data[i]); - psf->params->data[i] = poly; + psFree (psf->params_NEW->data[i]); + psf->params_NEW->data[i] = poly; } sprintf (keyword, "APTREND"); Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 9769) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 9770) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-24 22:55:05 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-28 20:23:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -111,4 +111,8 @@ pmSource *source = psfTry->sources->data[i]; pmModel *model = pmSourceModelGuess (source, psfTry->psf->type); + if (model == NULL) { + psError(PS_ERR_UNKNOWN, false, "failed to build model"); + return NULL; + } x = source->peak->x; y = source->peak->y; Index: /trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 9769) +++ /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 9770) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-19 21:16:49 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-28 20:23:51 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -35,4 +35,6 @@ // XXX make sure the angle in correctly translated to/from degrees // XXX we lose all information from the 'type' field + +// XXX update this file is we convert to PAR[4] : SigmaX*sqrt(2) (not 1/SigmaX) // elixir-style pseudo FITS table (header + ascii list) @@ -267,7 +269,7 @@ shape = psEllipseAxesToShape (axes); - PAR[4] = shape.sx; - PAR[5] = shape.sy; - PAR[6] = shape.sxy; + PAR[PM_PAR_SXX] = shape.sx; + PAR[PM_PAR_SYY] = shape.sy; + PAR[PM_PAR_SXY] = shape.sxy; source->modelPSF = model;