Index: /trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 10261) +++ /trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 10262) @@ -33,5 +33,4 @@ 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]; @@ -58,36 +57,121 @@ // define the parameter limits -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (7, PS_TYPE_F32); - *params_min = psVectorAlloc (7, PS_TYPE_F32); - *params_max = psVectorAlloc (7, PS_TYPE_F32); - - beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; - beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; - beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; - - params_min[0][0].data.F32[PM_PAR_SKY] = -1000; - params_min[0][0].data.F32[PM_PAR_I0] = 0; - params_min[0][0].data.F32[PM_PAR_XPOS] = -100; - params_min[0][0].data.F32[PM_PAR_YPOS] = -100; - params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; - params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; - - params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; - params_max[0][0].data.F32[PM_PAR_I0] = 1e8; - params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; - params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; - - return (TRUE); +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + float beta_lim, params_min, params_max; + float f1, f2, q1, q2; + + // we need to calculate the limits for SXY specially + if (nParam == PM_PAR_SXY) { + f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + assert (q1 > 0); + q2 = 0.5*sqrt (q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: + switch (nParam) { + case PM_PAR_SKY: + beta_lim = 1000; + break; + case PM_PAR_I0: + beta_lim = 3e6; + break; + case PM_PAR_XPOS: + beta_lim = 5; + break; + case PM_PAR_YPOS: + beta_lim = 5; + break; + case PM_PAR_SXX: + beta_lim = 0.5; + break; + case PM_PAR_SYY: + beta_lim = 0.5; + break; + case PM_PAR_SXY: + beta_lim = q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam); + } + if (fabs(beta[nParam]) > fabs(beta_lim)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); + return false; + } + return true; + case PS_MINIMIZE_PARAM_MIN: + switch (nParam) { + case PM_PAR_SKY: + params_min = -1000; + break; + case PM_PAR_I0: + params_min = 0; + break; + case PM_PAR_XPOS: + params_min = -100; + break; + case PM_PAR_YPOS: + params_min = -100; + break; + case PM_PAR_SXX: + params_min = 0.5; + break; + case PM_PAR_SYY: + params_min = 0.5; + break; + case PM_PAR_SXY: + params_min = -q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam); + } + if (params[nParam] < params_min) { + params[nParam] = params_min; + return false; + } + return true; + case PS_MINIMIZE_PARAM_MAX: + switch (nParam) { + case PM_PAR_SKY: + params_max = 1e5; + break; + case PM_PAR_I0: + params_max = 1e8; + break; + case PM_PAR_XPOS: + params_max = 1e4; + break; + case PM_PAR_YPOS: + params_max = 1e4; + break; + case PM_PAR_SXX: + params_max = 100; + break; + case PM_PAR_SYY: + params_max = 100; + break; + case PM_PAR_SXY: + params_max = +q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam); + } + if (params[nParam] > params_max) { + params[nParam] = params_max; + return false; + } + return true; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits"); + } + psAbort ("psModules.pmModel_GAUSS", "should not reach here"); + return false; } @@ -98,29 +182,20 @@ psF32 *PAR = model->params->data.F32; - # if (0) - - psEllipseMoments emoments; + psEllipseMoments emoments; emoments.x2 = moments->Sx; - emoments.y2 = moments->Sx; + emoments.y2 = moments->Sy; emoments.xy = moments->Sxy; - psEllipseAxes axes = psEllipseMomentsToAxes (emoments); + // force the axis ratio to be < 20.0 + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); psEllipseShape shape = psEllipseAxesToShape (axes); - # endif - - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + + 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, M_SQRT2*moments->Sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*moments->Sy); - PAR[PM_PAR_SXY] = 0.0; - # if (0) - - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); + PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - # endif - return(true); } @@ -138,5 +213,6 @@ // Area is equivalent to 2 pi sigma^2 - psEllipseAxes axes = psEllipseShapeToAxes (shape); + // axes ratio < 20 + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 Area = 2.0 * M_PI * axes.major * axes.minor; @@ -165,5 +241,5 @@ shape.sxy = PAR[PM_PAR_SXY]; - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); return (radius); @@ -194,4 +270,12 @@ out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); + // apply the model limits here: this truncates excessive extrapolation + for (int i = 0; i < psf->params_NEW->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params_NEW->data[i] == NULL) + continue; + PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, out, NULL); + PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, out, NULL); + } return(true); } Index: /trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 10261) +++ /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 10262) @@ -58,36 +58,122 @@ } -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (7, PS_TYPE_F32); - *params_min = psVectorAlloc (7, PS_TYPE_F32); - *params_max = psVectorAlloc (7, PS_TYPE_F32); - - beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; - beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; - beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; - - params_min[0][0].data.F32[PM_PAR_SKY] = -1000; - params_min[0][0].data.F32[PM_PAR_I0] = 0; - params_min[0][0].data.F32[PM_PAR_XPOS] = -100; - params_min[0][0].data.F32[PM_PAR_YPOS] = -100; - params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; - params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; - - params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; - params_max[0][0].data.F32[PM_PAR_I0] = 1e8; - params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; - params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; - - return (TRUE); +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + float beta_lim, params_min, params_max; + float f1, f2, q1, q2; + + // we need to calculate the limits for SXY specially + if (nParam == PM_PAR_SXY) { + f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + assert (q1 > 0); + q2 = 0.5*sqrt (q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: + switch (nParam) { + case PM_PAR_SKY: + beta_lim = 1000; + break; + case PM_PAR_I0: + beta_lim = 3e6; + break; + case PM_PAR_XPOS: + beta_lim = 5; + break; + case PM_PAR_YPOS: + beta_lim = 5; + break; + case PM_PAR_SXX: + beta_lim = 0.5; + break; + case PM_PAR_SYY: + beta_lim = 0.5; + break; + case PM_PAR_SXY: + beta_lim = q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam); + } + if (fabs(beta[nParam]) > fabs(beta_lim)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); + return false; + } + return true; + case PS_MINIMIZE_PARAM_MIN: + switch (nParam) { + case PM_PAR_SKY: + params_min = -1000; + break; + case PM_PAR_I0: + params_min = 0; + break; + case PM_PAR_XPOS: + params_min = -100; + break; + case PM_PAR_YPOS: + params_min = -100; + break; + case PM_PAR_SXX: + params_min = 0.5; + break; + case PM_PAR_SYY: + params_min = 0.5; + break; + case PM_PAR_SXY: + params_min = -q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam); + } + if (params[nParam] < params_min) { + params[nParam] = params_min; + return false; + } + return true; + case PS_MINIMIZE_PARAM_MAX: + switch (nParam) { + case PM_PAR_SKY: + params_max = 1e5; + break; + case PM_PAR_I0: + params_max = 1e8; + break; + case PM_PAR_XPOS: + params_max = 1e4; + break; + case PM_PAR_YPOS: + params_max = 1e4; + break; + case PM_PAR_SXX: + params_max = 100; + break; + case PM_PAR_SYY: + params_max = 100; + break; + case PM_PAR_SXY: + params_max = +q2; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam); + } + if (params[nParam] > params_max) { + params[nParam] = params_max; + return false; + } + return true; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits"); + } + psAbort ("psModules.pmModel_GAUSS", "should not reach here"); + return false; } @@ -98,14 +184,11 @@ psF32 *PAR = model->params->data.F32; - # if (0) - - psEllipseMoments emoments; + psEllipseMoments emoments; emoments.x2 = moments->Sx; emoments.y2 = moments->Sx; emoments.xy = moments->Sxy; - psEllipseAxes axes = psEllipseMomentsToAxes (emoments); + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); psEllipseShape shape = psEllipseAxesToShape (axes); - # endif PAR[PM_PAR_SKY] = moments->Sky; @@ -113,15 +196,7 @@ PAR[PM_PAR_XPOS] = moments->x; // XXX use peak->xf, peak->yf? PAR[PM_PAR_YPOS] = moments->y; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*moments->Sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*moments->Sy); - PAR[PM_PAR_SXY] = 0.0; - - # if (0) - - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); + PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - # endif - return(true); } @@ -139,5 +214,5 @@ // Area is equivalent to 2 pi sigma^2 - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 Area = 2.0 * M_PI * axes.major * axes.minor; @@ -183,5 +258,5 @@ // this estimates the radius assuming f(z) is roughly exp(-z) - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); Index: /trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 10261) +++ /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 10262) @@ -43,6 +43,7 @@ // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values // for negative values of z - if (z < 0) - z = 0; + // XXX use an assert here to force the elliptical parameters to be correctly determined + // if (z < 0) z = 0; + assert (z > 0); psF32 zp = pow(z,1.25); @@ -72,48 +73,149 @@ } -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); -} +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + float beta_lim, params_min, params_max; + float f1, f2, q1, q2; + + // we need to calculate the limits for SXY specially + if (nParam == PM_PAR_SXY) { + f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + assert (q1 > 0); + q2 = 0.5*sqrt (q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: + switch (nParam) { + case PM_PAR_SKY: + beta_lim = 1000; + break; + case PM_PAR_I0: + beta_lim = 3e6; + break; + case PM_PAR_XPOS: + beta_lim = 5; + break; + case PM_PAR_YPOS: + beta_lim = 5; + break; + case PM_PAR_SXX: + beta_lim = 0.5; + break; + case PM_PAR_SYY: + beta_lim = 0.5; + break; + case PM_PAR_SXY: + beta_lim = q2; + break; + case PM_PAR_7: + beta_lim = 0.5; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam); + } + if (fabs(beta[nParam]) > fabs(beta_lim)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); + return false; + } + return true; + case PS_MINIMIZE_PARAM_MIN: + switch (nParam) { + case PM_PAR_SKY: + params_min = -1000; + break; + case PM_PAR_I0: + params_min = 0; + break; + case PM_PAR_XPOS: + params_min = -100; + break; + case PM_PAR_YPOS: + params_min = -100; + break; + case PM_PAR_SXX: + params_min = 0.5; + break; + case PM_PAR_SYY: + params_min = 0.5; + break; + case PM_PAR_SXY: + params_min = -q2; + break; + case PM_PAR_7: + params_min = 0.1; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam); + } + if (params[nParam] < params_min) { + params[nParam] = params_min; + return false; + } + return true; + case PS_MINIMIZE_PARAM_MAX: + switch (nParam) { + case PM_PAR_SKY: + params_max = 1e5; + break; + case PM_PAR_I0: + params_max = 1e8; + break; + case PM_PAR_XPOS: + params_max = 1e4; + break; + case PM_PAR_YPOS: + params_max = 1e4; + break; + case PM_PAR_SXX: + params_max = 100; + break; + case PM_PAR_SYY: + params_max = 100; + break; + case PM_PAR_SXY: + params_max = +q2; + break; + case PM_PAR_7: + params_max = 10.0; + break; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam); + } + if (params[nParam] > params_max) { + params[nParam] = params_max; + return false; + } + return true; + default: + psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits"); + } + psAbort ("psModules.pmModel_GAUSS", "should not reach here"); + return false; +} + // make an initial guess for parameters bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; pmPeak *peak = source->peak; psF32 *PAR = model->params->data.F32; + + psEllipseMoments emoments; + emoments.x2 = moments->Sx; + emoments.y2 = moments->Sy; + emoments.xy = moments->Sxy; + + // force the axis ratio to be < 20.0 + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); + psEllipseShape shape = psEllipseAxesToShape (axes); PAR[PM_PAR_SKY] = moments->Sky; @@ -121,7 +223,7 @@ PAR[PM_PAR_XPOS] = peak->x; PAR[PM_PAR_YPOS] = peak->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; + PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); + PAR[PM_PAR_SXY] = shape.sxy; PAR[PM_PAR_7] = 1.0; @@ -141,5 +243,5 @@ // Area is equivalent to 2 pi sigma^2 - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 Area = 2.0 * M_PI * axes.major * axes.minor; @@ -186,5 +288,5 @@ shape.sxy = PAR[PM_PAR_SXY]; - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 sigma = axes.major;