Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 14322) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 14323) @@ -28,5 +28,5 @@ # define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS -psF64 PS_MODEL_FUNC (psVector *deriv, +psF32 PM_MODEL_FUNC (psVector *deriv, const psVector *params, const psVector *pixcoord) @@ -39,4 +39,6 @@ psF32 py = Y / PAR[PM_PAR_SYY]; psF32 z = PS_SQR(px) + PS_SQR(py) + X*Y*PAR[PM_PAR_SXY]; + + assert (z >= 0); psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); @@ -58,76 +60,186 @@ dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; dPAR[PM_PAR_SXY] = -q*X*Y; - dPAR[PM_PAR_7] = -5.0*t*log(z)*p*z; + + // this model derivative is undefined at z = 0.0, but is actually 0.0 + dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -5.0*t*log(z)*p*z; } 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] = 1.25; - - 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] = 4.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) +// 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 = 0, params_min = 0, params_max = 0; + float f1 = 0, f2 = 0, q1 = 0, q2 = 0; + + // 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); + q1 = (q1 < 0.0) ? 0.0 : q1; + // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that + // angle and let f2,f1 fight it out + 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 = 0.5*q2; + break; + case PM_PAR_7: + beta_lim = 0.5; + break; + default: + psAbort("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); + psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], 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 = 1.25; + break; + default: + psAbort("invalid parameter %d for param min test", nParam); + } + if (params[nParam] < params_min) { + params[nParam] = params_min; + psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, 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 = 4.0; + break; + default: + psAbort("invalid parameter %d for param max test", nParam); + } + if (params[nParam] > params_max) { + params[nParam] = params_max; + psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], params_max); + return false; + } + return true; + default: + psAbort("invalid choice for limits"); + } + psAbort("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; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments tmpMoments; - - // XXX fix this stuff : should be using correct ellipse relationships... - tmpMoments.x2 = PS_SQR(source->moments->Sx); - tmpMoments.y2 = PS_SQR(source->moments->Sy); - tmpMoments.xy = source->moments->Sxy; - - axes = psEllipseMomentsToAxes(tmpMoments); - shape = psEllipseAxesToShape(axes); - - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + 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); + + if (!isfinite(axes.major)) return false; + if (!isfinite(axes.minor)) return false; + if (!isfinite(axes.theta)) return false; + + psEllipseShape shape = psEllipseAxesToShape (axes); + + if (!isfinite(shape.sx)) return false; + if (!isfinite(shape.sy)) return false; + if (!isfinite(shape.sxy)) return false; + + PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_I0] = moments->Peak - moments->Sky; 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_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] = 2.0; + PAR[PM_PAR_7] = 2.25; + return(true); } -psF64 PS_MODEL_FLUX (const psVector *params) +psF64 PM_MODEL_FLUX (const psVector *params) { float norm, z; @@ -136,10 +248,10 @@ psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); - shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; shape.sxy = PAR[PM_PAR_SXY]; // 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; @@ -167,7 +279,8 @@ // 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 z, f, p; +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF64 z, f; + int Nstep = 0; psEllipseShape shape; @@ -181,12 +294,11 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0); - shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0); + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; shape.sxy = PAR[PM_PAR_SXY]; - psEllipseAxes axes = psEllipseShapeToAxes (shape); + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 sigma = axes.major; - psF64 dz = 1.0 / (2.0 * sigma*sigma); psF64 limit = flux / PAR[PM_PAR_I0]; @@ -196,6 +308,6 @@ // choose a z value guaranteed to be beyond our limit - float z0 = pow((1.0 / limit), (1.0 / 2.25)); - float z1 = (1.0 / limit) / PAR[PM_PAR_7]; + float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7])); + float z1 = (1.0 / limit); z1 = PS_MAX (z0, z1); z0 = 0.0; @@ -224,5 +336,5 @@ } -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { @@ -243,8 +355,30 @@ } - // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here - out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); - - return(true); + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + if (!pmPSF_FitToModel (out, 0.1)) { + psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", + in[PM_PAR_YPOS], in[PM_PAR_XPOS]); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + 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; + + bool status = true; + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL); + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", + in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + modelPSF->flags |= PM_MODEL_STATUS_LIMITS; + } + } + + return true; }