Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 29004) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 29004) @@ -0,0 +1,472 @@ +/****************************************************************************** + * this file defines the DEV 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: + + f = exp(-z^n) + + * 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 - normalized dev parameter + + note that a standard dev model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, + the -1 offset) are absorbed in this model by the normalization, the exponent, and the + radial scale. We fit the elements in this form, then re-normalize them on output. + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + +#include "pmModel_DEV.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_DEV +# define PM_MODEL_FLUX pmModelFlux_DEV +# define PM_MODEL_GUESS pmModelGuess_DEV +# define PM_MODEL_LIMITS pmModelLimits_DEV +# define PM_MODEL_RADIUS pmModelRadius_DEV +# define PM_MODEL_FROM_PSF pmModelFromPSF_DEV +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_DEV +# define PM_MODEL_FIT_STATUS pmModelFitStatus_DEV +# define PM_MODEL_SET_LIMITS pmModelSetLimits_DEV + +// f = exp(-z^0.125) +# define ALPHA 0.125 +// # define ALPHA 0.25 + +// 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.001, 0.001, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// Moderate parameter limits +static float *paramsMinModerate = paramsMinLax; +static float *paramsMaxModerate = paramsMaxLax; + +// Strict parameter limits +static float *paramsMinStrict = paramsMinLax; +static float *paramsMaxStrict = paramsMaxLax; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 }; + +static bool limitsApply = true; // Apply limits? + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + float index = 0.5 / ALPHA; + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + 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); + + assert (z >= 0); + + psF32 f2 = bn*pow(z,ALPHA); + psF32 f1 = Io*exp(-f2); + psF32 z0 = PAR[PM_PAR_I0]*f1; + psF32 f0 = PAR[PM_PAR_SKY] + z0; + + assert (isfinite(f2)); + assert (isfinite(f1)); + assert (isfinite(z0)); + assert (isfinite(f0)); + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + // gradient is infinite for z = 0; saturate at z = 0.01 + psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0); + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +2.0*f1; + + assert (isfinite(z1)); + + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } + return (f0); +} + +// 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) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + // we need to calculate the limits for SXY specially + float q2 = NAN; + if (nParam == PM_PAR_SXY) { + float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + float 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*sqrtf(q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + +// 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) +{ + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *PAR = model->params->data.F32; + + psEllipseMoments emoments; + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + // 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; + + // the other parameters depend on the guess for PAR_7 + float index = 0.5 / ALPHA; + float bn = 1.9992*index - 0.3271; + // float fR = 1.0 / (sqrt(2.0) * pow (bn, index)); + float Io = exp(0.5*bn); + + float Sxx = PS_MAX(0.5, shape.sx); + float Syy = PS_MAX(0.5, shape.sy); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = peak->flux / Io; + PAR[PM_PAR_XPOS] = peak->xf; + PAR[PM_PAR_YPOS] = peak->yf; + // PAR[PM_PAR_SXX] = Sxx * fR; + // PAR[PM_PAR_SYY] = Syy * fR; + PAR[PM_PAR_SXX] = Sxx; + PAR[PM_PAR_SYY] = Syy; + PAR[PM_PAR_SXY] = shape.sxy; + + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + 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, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = exp(-pow(z,ALPHA)); + z += DZ; + f2 = exp(-pow(z,ALPHA)); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// 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) +{ + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + if (flux <= 0) + return (1.0); + if (PAR[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= PAR[PM_PAR_I0]) + return (1.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, 20.0); + + // f = Io exp(-z^n) -> z^n = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); + + psAssert (isfinite(radius), "fix this code: z should not be nan"); + 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; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + out[i] = in[i]; + } else { + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + if (!pmPSF_FitToModel (out, 0.1)) { + psTrace("psModules.objects", 5, "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->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->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; +} + +// construct the PSF model from the FLT model and the psf +// XXX is this sufficiently general do be a global function, not a pmModelClass function? +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psF32 *PAR = model->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = Io; + PAR[PM_PAR_XPOS] = Xo; + PAR[PM_PAR_YPOS] = Yo; + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params->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; + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + // XXX user-defined value for limit? + if (!pmPSF_FitToModel (PAR, 0.1)) { + psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); + 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->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); + model->flags |= PM_MODEL_STATUS_LIMITS; + } + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: trunk/psModules/src/objects/models/pmModel_DEV.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.h (revision 29004) +++ trunk/psModules/src/objects/models/pmModel_DEV.h (revision 29004) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_DEV_H + +#include "pmModel.h" + +psF32 pmModelFunc_DEV(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_DEV(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_DEV(pmModel *model, pmSource *source); +psF64 pmModelFlux_DEV(const psVector *params); +psF64 pmModelRadius_DEV(const psVector *params, psF64 flux); +bool pmModelFromPSF_DEV(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_DEV(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_DEV(pmModel *model); +void pmModelSetLimits_DEV(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 29004) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 29004) @@ -0,0 +1,453 @@ +/****************************************************************************** + * this file defines the Exponential (EXP) 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: + + f = exp(-sqrt(z)) (since z is r^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 + + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + +#include "pmModel_EXP.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_EXP +# define PM_MODEL_FLUX pmModelFlux_EXP +# define PM_MODEL_GUESS pmModelGuess_EXP +# define PM_MODEL_LIMITS pmModelLimits_EXP +# define PM_MODEL_RADIUS pmModelRadius_EXP +# define PM_MODEL_FROM_PSF pmModelFromPSF_EXP +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_EXP +# define PM_MODEL_FIT_STATUS pmModelFitStatus_EXP +# define PM_MODEL_SET_LIMITS pmModelSetLimits_EXP + +// 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.05, 0.05, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// Moderate parameter limits +static float *paramsMinModerate = paramsMinLax; +static float *paramsMaxModerate = paramsMaxLax; + +// Strict parameter limits +static float *paramsMinStrict = paramsMinLax; +static float *paramsMaxStrict = paramsMaxLax; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5}; + +static bool limitsApply = true; // Apply limits? + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + 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; + + // XXX if the elliptical contour is defined in valid way, this step should not be required. + // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values + // for negative values of z + // XXX use an assert here to force the elliptical parameters to be correctly determined + // if (z < 0) z = 0; + assert (z >= 0); + + psF32 f2 = sqrt(z); + psF32 f1 = exp(-f2); + psF32 z0 = PAR[PM_PAR_I0]*f1; + psF32 f0 = PAR[PM_PAR_SKY] + z0; + + assert (isfinite(f2)); + assert (isfinite(f1)); + assert (isfinite(z0)); + assert (isfinite(f0)); + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +f1; + + // gradient is infinite for z = 0; saturate at z = 0.01 + // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i] + psF32 z1 = (z < 0.01) ? 0.5*z0/sqrt(0.01) : 0.5*z0/sqrt(z); + + // XXX dampen SXX and SYY as in GAUSS? + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } + return (f0); +} + +// 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) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + // we need to calculate the limits for SXY specially + float q2 = NAN; + if (nParam == PM_PAR_SXY) { + float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + float 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*sqrtf(q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + +// 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) +{ + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *PAR = model->params->data.F32; + + psEllipseMoments emoments; + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + // 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] = 0.0; + PAR[PM_PAR_I0] = peak->flux; + PAR[PM_PAR_XPOS] = peak->xf; + PAR[PM_PAR_YPOS] = peak->yf; + 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; + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + 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, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity] + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = exp(-sqrt(z)); + z += DZ; + f2 = exp(-sqrt(z)); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// 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) +{ + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + if (flux <= 0) + return (1.0); + if (PAR[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= PAR[PM_PAR_I0]) + return (1.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, 20.0); + + // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + psF64 radius = axes.major * sqrt (2.0) * zn; + + psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); + 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; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + out[i] = in[i]; + } else { + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + if (!pmPSF_FitToModel (out, 0.1)) { + psTrace("psModules.objects", 5, "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->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->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; +} + +// construct the PSF model from the FLT model and the psf +// XXX is this sufficiently general do be a global function, not a pmModelClass function? +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psF32 *PAR = model->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = Io; + PAR[PM_PAR_XPOS] = Xo; + PAR[PM_PAR_YPOS] = Yo; + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params->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; + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + // XXX user-defined value for limit? + if (!pmPSF_FitToModel (PAR, 0.1)) { + psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); + 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->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); + model->flags |= PM_MODEL_STATUS_LIMITS; + } + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: trunk/psModules/src/objects/models/pmModel_EXP.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.h (revision 29004) +++ trunk/psModules/src/objects/models/pmModel_EXP.h (revision 29004) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_EXP_H + +#include "pmModel.h" + +psF32 pmModelFunc_EXP(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_EXP(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_EXP(pmModel *model, pmSource *source); +psF64 pmModelFlux_EXP(const psVector *params); +psF64 pmModelRadius_EXP(const psVector *params, psF64 flux); +bool pmModelFromPSF_EXP(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_EXP(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_EXP(pmModel *model); +void pmModelSetLimits_EXP(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 29004) @@ -21,11 +21,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_GAUSS.h" +# define PM_MODEL_NPARAM 7 # define PM_MODEL_FUNC pmModelFunc_GAUSS # define PM_MODEL_FLUX pmModelFlux_GAUSS @@ -83,4 +103,5 @@ 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]); + // the extra factor of 2 below is needed to avoid excessive swings dPAR[PM_PAR_SXX] = +4.0*q*px*px/PAR[PM_PAR_SXX]; @@ -102,5 +123,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -347,4 +368,5 @@ // this test is invalid if the parameters are derived // from the PSF model +// XXX how is this used? it prevents forced photometry from ever being 'successful' bool PM_MODEL_FIT_STATUS (pmModel *model) { @@ -394,12 +416,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 29004) @@ -21,11 +21,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_PGAUSS.h" +# define PM_MODEL_NPARAM 7 # define PM_MODEL_FUNC pmModelFunc_PGAUSS # define PM_MODEL_FLUX pmModelFlux_PGAUSS @@ -103,5 +123,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -448,12 +468,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 29004) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_PS1_V1.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_PS1_V1 # define PM_MODEL_FLUX pmModelFlux_PS1_V1 @@ -122,5 +142,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -222,5 +242,5 @@ PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 1.0; + PAR[PM_PAR_7] = 0.5; return(true); @@ -468,14 +488,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS -# undef ALPHA -# undef ALPHA_M Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 29004) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_QGAUSS.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_QGAUSS # define PM_MODEL_FLUX pmModelFlux_QGAUSS @@ -123,5 +143,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -469,14 +489,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS -# undef ALPHA -# undef ALPHA_M Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 29004) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_RGAUSS.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_RGAUSS # define PM_MODEL_FLUX pmModelFlux_RGAUSS @@ -112,5 +132,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -209,8 +229,8 @@ PAR[PM_PAR_XPOS] = peak->xf; PAR[PM_PAR_YPOS] = peak->yf; - 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_SXX] = PS_MAX(0.5, shape.sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 2.25; + PAR[PM_PAR_7] = 1.5; return(true); @@ -463,12 +483,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 27565) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 29004) @@ -18,5 +18,5 @@ * PM_PAR_7 7 - normalized sersic parameter - note that a standard sersic model uses exp(-K*(z^(1/n) - 1). the additional elements (K, + note that a standard sersic model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, the -1 offset) are absorbed in this model by the normalization, the exponent, and the radial scale. We fit the elements in this form, then re-normalize them on output. @@ -25,11 +25,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" #include "pmMoments.h" -#include "pmPeaks.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_SERSIC.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_SERSIC # define PM_MODEL_FLUX pmModelFlux_SERSIC @@ -47,5 +67,5 @@ // Lax parameter limits -static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 }; +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 }; static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; @@ -84,6 +104,10 @@ assert (z >= 0); - psF32 f2 = pow(z,PAR[PM_PAR_7]); - psF32 f1 = exp(-f2); + float index = 0.5 / PAR[PM_PAR_7]; + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + psF32 f2 = bn*pow(z,PAR[PM_PAR_7]); + psF32 f1 = Io*exp(-f2); psF32 z0 = PAR[PM_PAR_I0]*f1; psF32 f0 = PAR[PM_PAR_SKY] + z0; @@ -98,9 +122,10 @@ // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0); + psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0); dPAR[PM_PAR_SKY] = +1.0; dPAR[PM_PAR_I0] = +f1; - dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -z0*f2*log(z); + dPAR[PM_PAR_7] = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z); + dPAR[PM_PAR_7] *= 3.0; assert (isfinite(z1)); @@ -109,7 +134,6 @@ dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag? dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; } @@ -127,5 +151,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -201,4 +225,15 @@ psF32 *PAR = model->params->data.F32; + // the other parameters depend on the guess for PAR_7 + if (!isfinite(PAR[PM_PAR_7])) { + PAR[PM_PAR_7] = 0.25; + } + float index = 0.5 / PAR[PM_PAR_7]; + + // the scale-length is a function of the moments and the index: + // Rmajor_guess = Rmajor_moments * Scale * Zero + float Scale = 0.70 + 0.053 * PAR[PM_PAR_7]; + float Zero = 1.16 - 0.615 * PAR[PM_PAR_7]; + psEllipseMoments emoments; emoments.x2 = moments->Mxx; @@ -213,4 +248,12 @@ if (!isfinite(axes.theta)) return false; + // set a lower limit to avoid absurd solutions.. + float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero); + float Rminor = Rmajor * (axes.minor / axes.major); + + fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor); + + axes.major = Rmajor; + axes.minor = Rminor; psEllipseShape shape = psEllipseAxesToShape (axes); @@ -219,12 +262,22 @@ if (!isfinite(shape.sxy)) return false; + float bn = 1.9992*index - 0.3271; + // float fR = 1.0 / (sqrt(2.0) * pow (bn, index)); + float Io = exp(0.5*bn); + + // XXX do we need this factor of sqrt(2)? + // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); + // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + + float Sxx = PS_MAX(0.5, shape.sx); + float Syy = PS_MAX(0.5, shape.sy); + PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; + PAR[PM_PAR_I0] = peak->flux / Io; PAR[PM_PAR_XPOS] = peak->xf; PAR[PM_PAR_YPOS] = peak->yf; - 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_SXX] = Sxx; + PAR[PM_PAR_SYY] = Syy; PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 0.5; return(true); @@ -254,7 +307,8 @@ float f1, f2; for (z = DZ; z < 50; z += DZ) { - f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + f1 = exp(-pow(z,PAR[PM_PAR_7])); z += DZ; - f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + f2 = exp(-pow(z,PAR[PM_PAR_7])); norm += f0 + 4*f1 + f2; f0 = f2; @@ -287,17 +341,12 @@ psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 sigma = axes.major; - - psF64 limit = flux / PAR[PM_PAR_I0]; - - 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)) - psAbort("error in code: radius is NaN"); - + + // f = Io exp(-z^n) -> z^n = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]); + + fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major); + + psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); return (radius); } @@ -448,12 +497,2 @@ return; } - -# 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_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS