Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 25477) +++ trunk/psModules/src/objects/Makefile.am (revision 25738) @@ -54,11 +54,9 @@ pmGrowthCurve.c \ pmSourceMatch.c \ - pmDetEff.c - -EXTRA_DIST = \ + pmDetEff.c \ models/pmModel_GAUSS.c \ models/pmModel_PGAUSS.c \ + models/pmModel_PS1_V1.c \ models/pmModel_QGAUSS.c \ - models/pmModel_SGAUSS.c \ models/pmModel_RGAUSS.c \ models/pmModel_SERSIC.c @@ -92,5 +90,11 @@ pmGrowthCurve.h \ pmSourceMatch.h \ - pmDetEff.h + pmDetEff.h \ + models/pmModel_GAUSS.h \ + models/pmModel_PGAUSS.h \ + models/pmModel_PS1_V1.h \ + models/pmModel_QGAUSS.h \ + models/pmModel_RGAUSS.h \ + models/pmModel_SERSIC.h CLEANFILES = *~ Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 25738) @@ -19,4 +19,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_GAUSS.h" + # define PM_MODEL_FUNC pmModelFunc_GAUSS # define PM_MODEL_FLUX pmModelFlux_GAUSS @@ -27,4 +36,20 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_GAUSS # define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS +# define PM_MODEL_SET_LIMITS pmModelSetLimits_GAUSS + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// 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, 2.0, 2.0, 0.5 }; + +static bool limitsApply = true; // Apply limits? // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) @@ -68,118 +93,66 @@ 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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 = PS_MAX (0.0, q1); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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 = 2.0; - break; - case PM_PAR_SYY: - beta_lim = 2.0; - break; - case PM_PAR_SXY: - beta_lim = 0.5*q2; - 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.01; - 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("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; - 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; + 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"); @@ -388,4 +361,33 @@ } +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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} + # undef PM_MODEL_FUNC # undef PM_MODEL_FLUX @@ -396,2 +398,3 @@ # undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS +# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_GAUSS.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_GAUSS_H + +#include "pmModel.h" + +psF32 pmModelFunc_GAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_GAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_GAUSS(pmModel *model, pmSource *source); +psF64 pmModelFlux_GAUSS(const psVector *params); +psF64 pmModelRadius_GAUSS(const psVector *params, psF64 flux); +bool pmModelFromPSF_GAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_GAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_GAUSS(pmModel *model); +void pmModelSetLimits_GAUSS(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 25738) @@ -19,4 +19,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_PGAUSS.h" + # define PM_MODEL_FUNC pmModelFunc_PGAUSS # define PM_MODEL_FLUX pmModelFlux_PGAUSS @@ -27,4 +36,20 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS # define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS +# define PM_MODEL_SET_LIMITS pmModelSetLimits_PGAUSS + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// 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, 2.0, 2.0, 0.5 }; + +static bool limitsApply = true; // Apply limits? // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) @@ -69,119 +94,67 @@ 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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 = PS_MAX (0.0, q1); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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 = 2.0; - break; - case PM_PAR_SYY: - beta_lim = 2.0; - break; - case PM_PAR_SXY: - beta_lim = 0.5*q2; - 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.01; - 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("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; - 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: + 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"); } @@ -189,4 +162,5 @@ return false; } + // make an initial guess for parameters @@ -434,4 +408,34 @@ } + +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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} + # undef PM_MODEL_FUNC # undef PM_MODEL_FLUX @@ -442,2 +446,3 @@ # undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS +# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_PGAUSS_H + +#include "pmModel.h" + +psF32 pmModelFunc_PGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_PGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_PGAUSS(pmModel *model, pmSource *source); +psF64 pmModelFlux_PGAUSS(const psVector *params); +psF64 pmModelRadius_PGAUSS(const psVector *params, psF64 flux); +bool pmModelFromPSF_PGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_PGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_PGAUSS(pmModel *model); +void pmModelSetLimits_PGAUSS(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 25738) @@ -20,4 +20,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_PS1_V1.h" + # define PM_MODEL_FUNC pmModelFunc_PS1_V1 # define PM_MODEL_FLUX pmModelFlux_PS1_V1 @@ -28,7 +37,25 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PS1_V1 # define PM_MODEL_FIT_STATUS pmModelFitStatus_PS1_V1 +# define PM_MODEL_SET_LIMITS pmModelSetLimits_PS1_V1 # define ALPHA 1.666 # define ALPHA_M 0.666 + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; + +// Strict parameter limits +// k = PAR_7 < 0 is very undesirable (big divot in the middle) +static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; + +// 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, 2.0 }; + +static bool limitsApply = true; // Apply limits? + psF32 PM_MODEL_FUNC (psVector *deriv, @@ -84,128 +111,67 @@ 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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 = 1.0; - break; - case PM_PAR_SYY: - beta_lim = 1.0; - break; - case PM_PAR_SXY: - beta_lim = 0.5*q2; - break; - case PM_PAR_7: - beta_lim = 2.0; - 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.01; - 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.0; - 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 = 20.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: + 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"); } @@ -299,10 +265,16 @@ 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); + if (flux <= 0) { + return 1.0; + } + if (PAR[PM_PAR_I0] <= 0) { + return 1.0; + } + if (flux >= PAR[PM_PAR_I0]) { + return 1.0; + } + if (PAR[PM_PAR_7] == 0.0) { + return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); + } shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; @@ -468,4 +440,35 @@ } + +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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} + + # undef PM_MODEL_FUNC # undef PM_MODEL_FLUX @@ -476,4 +479,5 @@ # 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_PS1_V1.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_PS1_V1_H + +#include "pmModel.h" + +psF32 pmModelFunc_PS1_V1(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_PS1_V1(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_PS1_V1(pmModel *model, pmSource *source); +psF64 pmModelFlux_PS1_V1(const psVector *params); +psF64 pmModelRadius_PS1_V1(const psVector *params, psF64 flux); +bool pmModelFromPSF_PS1_V1(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_PS1_V1(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_PS1_V1(pmModel *model); +void pmModelSetLimits_PS1_V1(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 25738) @@ -20,4 +20,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_QGAUSS.h" + # define PM_MODEL_FUNC pmModelFunc_QGAUSS # define PM_MODEL_FLUX pmModelFlux_QGAUSS @@ -28,4 +37,20 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS # define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS +# define PM_MODEL_SET_LIMITS pmModelSetLimits_QGAUSS + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; + +// 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, @@ -79,129 +104,69 @@ # 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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 = 1.0; - break; - case PM_PAR_SYY: - beta_lim = 1.0; - break; - case PM_PAR_SXY: - beta_lim = 0.5*q2; - break; - case PM_PAR_7: - beta_lim = 2.0; - 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.01; - 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("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 = 20.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; + 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"); @@ -464,4 +429,34 @@ } + +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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} + # undef PM_MODEL_FUNC # undef PM_MODEL_FLUX @@ -472,2 +467,3 @@ # undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS +# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_QGAUSS_H + +#include "pmModel.h" + +psF32 pmModelFunc_QGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_QGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_QGAUSS(pmModel *model, pmSource *source); +psF64 pmModelFlux_QGAUSS(const psVector *params); +psF64 pmModelRadius_QGAUSS(const psVector *params, psF64 flux); +bool pmModelFromPSF_QGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_QGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_QGAUSS(pmModel *model); +void pmModelSetLimits_QGAUSS(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 25738) @@ -20,4 +20,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_RGAUSS.h" + # define PM_MODEL_FUNC pmModelFunc_RGAUSS # define PM_MODEL_FLUX pmModelFlux_RGAUSS @@ -28,4 +37,20 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS # define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS +# define PM_MODEL_SET_LIMITS pmModelSetLimits_RGAUSS + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; + +// 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, 0.5, 0.5, 0.5, 0.5 }; + +static bool limitsApply = true; // Apply limits? psF32 PM_MODEL_FUNC (psVector *deriv, @@ -73,130 +98,70 @@ # 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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.01; - 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: + 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"); } @@ -204,4 +169,5 @@ return false; } + // make an initial guess for parameters @@ -456,4 +422,34 @@ } + +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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} + # undef PM_MODEL_FUNC # undef PM_MODEL_FLUX @@ -464,2 +460,3 @@ # undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS +# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_RGAUSS_H + +#include "pmModel.h" + +psF32 pmModelFunc_RGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_RGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_RGAUSS(pmModel *model, pmSource *source); +psF64 pmModelFlux_RGAUSS(const psVector *params); +psF64 pmModelRadius_RGAUSS(const psVector *params, psF64 flux); +bool pmModelFromPSF_RGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_RGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_RGAUSS(pmModel *model); +void pmModelSetLimits_RGAUSS(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 25477) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 25738) @@ -23,4 +23,13 @@ *****************************************************************************/ +#include +#include + +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmModel.h" +#include "pmModel_SERSIC.h" + # define PM_MODEL_FUNC pmModelFunc_SERSIC # define PM_MODEL_FLUX pmModelFlux_SERSIC @@ -31,4 +40,20 @@ # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC # define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC +# define PM_MODEL_SET_LIMITS pmModelSetLimits_SERSIC + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; + +// 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, 2.0 }; + +static bool limitsApply = true; // Apply limits? psF32 PM_MODEL_FUNC (psVector *deriv, @@ -91,128 +116,67 @@ 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; + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially + float q2 = NAN; 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); + 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*sqrt (q1); + q2 = 0.5*sqrtf(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 = 1.0; - break; - case PM_PAR_SYY: - beta_lim = 1.0; - break; - case PM_PAR_SXY: - beta_lim = 0.5*q2; - break; - case PM_PAR_7: - beta_lim = 2.0; - 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.01; - break; - case PM_PAR_XPOS: - params_min = -100; - break; - case PM_PAR_YPOS: - params_min = -100; - break; - case PM_PAR_SXX: - params_min = 0.05; - break; - case PM_PAR_SYY: - params_min = 0.05; - break; - case PM_PAR_SXY: - params_min = -q2; - break; - case PM_PAR_7: - params_min = 0.05; - 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: + 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"); } @@ -220,5 +184,4 @@ return false; } - // make an initial guess for parameters @@ -447,7 +410,37 @@ fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n", - dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); + dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); 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_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; } @@ -460,2 +453,3 @@ # undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS +# undef PM_MODEL_SET_LIMITS Index: trunk/psModules/src/objects/models/pmModel_SERSIC.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.h (revision 25738) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.h (revision 25738) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_SERSIC_H + +#include "pmModel.h" + +psF32 pmModelFunc_SERSIC(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_SERSIC(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_SERSIC(pmModel *model, pmSource *source); +psF64 pmModelFlux_SERSIC(const psVector *params); +psF64 pmModelRadius_SERSIC(const psVector *params, psF64 flux); +bool pmModelFromPSF_SERSIC(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_SERSIC(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_SERSIC(pmModel *model); +void pmModelSetLimits_SERSIC(pmModelLimitsType type); + +#endif Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 25477) +++ trunk/psModules/src/objects/pmModel.c (revision 25738) @@ -86,4 +86,5 @@ tmp->modelParamsFromPSF = class->modelParamsFromPSF; tmp->modelFitStatus = class->modelFitStatus; + tmp->modelSetLimits = class->modelSetLimits; psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); @@ -235,14 +236,14 @@ if (model->residuals) { - DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5; - DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5; - Ro = (model->residuals->Ro) ? model->residuals->Ro->data.F32 : NULL; - Rx = (model->residuals->Rx) ? model->residuals->Rx->data.F32 : NULL; - Ry = (model->residuals->Ry) ? model->residuals->Ry->data.F32 : NULL; - Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL; - if (Ro) { - NX = model->residuals->Ro->numCols; - NY = model->residuals->Ro->numRows; - } + DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5; + DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5; + Ro = (model->residuals->Ro) ? model->residuals->Ro->data.F32 : NULL; + Rx = (model->residuals->Rx) ? model->residuals->Rx->data.F32 : NULL; + Ry = (model->residuals->Ry) ? model->residuals->Ry->data.F32 : NULL; + Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL; + if (Ro) { + NX = model->residuals->Ro->numCols; + NY = model->residuals->Ro->numRows; + } } @@ -256,5 +257,5 @@ // XXX should we use using 0.5 pixel offset? - // Convert to coordinate in parent image, with offset (dx,dy) + // Convert to coordinate in parent image, with offset (dx,dy) imageCol = ix + image->col0 - dx; imageRow = iy + image->row0 - dy; @@ -276,56 +277,56 @@ float rx = xBin*ix + DX; - int rx0 = rx - 0.5; - int rx1 = rx + 0.5; - int ry0 = ry - 0.5; - int ry1 = ry + 0.5; - - if (rx0 < 0) goto skip; - if (ry0 < 0) goto skip; - if (rx1 >= NX) goto skip; - if (ry1 >= NY) goto skip; - - // these go from 0.0 to 1.0 between the centers of the pixels - float fx = rx - 0.5 - rx0; - float Fx = 1.0 - fx; - float fy = ry - 0.5 - ry0; - float Fy = 1.0 - fy; - - // check the residual image mask (if set). give up if any of the 4 pixels are masked. - if (Rm) { - if (Rm[ry0][rx0]) goto skip; - if (Rm[ry0][rx1]) goto skip; - if (Rm[ry1][rx0]) goto skip; - if (Rm[ry1][rx1]) goto skip; - } - - // a possible further optimization if we re-use these values - // XXX allow for masked pixels, and add pixel weights - float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx); - float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx); - float Vo = V0*Fy + V1*fy; - if (!isfinite(Vo)) goto skip; - - float Vx = 0.0; - float Vy = 0.0; - - // skip Rx,Ry if Ro is masked - if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) { - V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx); - V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx); - Vx = V0*Fy + V1*fy; - - V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx); - V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx); - Vy = V0*Fy + V1*fy; - } - if (!isfinite(Vx)) goto skip; - if (!isfinite(Vy)) goto skip; - - // 2D residual variations are set for the true source position - pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy); + int rx0 = rx - 0.5; + int rx1 = rx + 0.5; + int ry0 = ry - 0.5; + int ry1 = ry + 0.5; + + if (rx0 < 0) goto skip; + if (ry0 < 0) goto skip; + if (rx1 >= NX) goto skip; + if (ry1 >= NY) goto skip; + + // these go from 0.0 to 1.0 between the centers of the pixels + float fx = rx - 0.5 - rx0; + float Fx = 1.0 - fx; + float fy = ry - 0.5 - ry0; + float Fy = 1.0 - fy; + + // check the residual image mask (if set). give up if any of the 4 pixels are masked. + if (Rm) { + if (Rm[ry0][rx0]) goto skip; + if (Rm[ry0][rx1]) goto skip; + if (Rm[ry1][rx0]) goto skip; + if (Rm[ry1][rx1]) goto skip; + } + + // a possible further optimization if we re-use these values + // XXX allow for masked pixels, and add pixel weights + float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx); + float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx); + float Vo = V0*Fy + V1*fy; + if (!isfinite(Vo)) goto skip; + + float Vx = 0.0; + float Vy = 0.0; + + // skip Rx,Ry if Ro is masked + if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) { + V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx); + V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx); + Vx = V0*Fy + V1*fy; + + V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx); + V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx); + Vy = V0*Fy + V1*fy; + } + if (!isfinite(Vx)) goto skip; + if (!isfinite(Vy)) goto skip; + + // 2D residual variations are set for the true source position + pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy); } - skip: + skip: // add or subtract the value if (add) { Index: trunk/psModules/src/objects/pmModel.h =================================================================== --- trunk/psModules/src/objects/pmModel.h (revision 25477) +++ trunk/psModules/src/objects/pmModel.h (revision 25738) @@ -44,4 +44,12 @@ } pmModelOpMode; +/// Parameter limit types +typedef enum { + PM_MODEL_LIMITS_NONE, ///< Apply no limits: suitable for debugging + PM_MODEL_LIMITS_IGNORE, ///< Ignore all limits: fit can go to town + PM_MODEL_LIMITS_LAX, ///< Lax limits: attempting to reproduce mildly bad data + PM_MODEL_LIMITS_STRICT, ///< Strict limits: good quality data +} pmModelLimitsType; + typedef struct pmModel pmModel; typedef struct pmSource pmSource; @@ -74,4 +82,7 @@ // This function returns the success / failure status of the given model fit typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +// This function sets the parameter limits for the given model +typedef bool (*pmModelSetLimitsFunc)(pmModelLimitsType type); /** pmModel data structure @@ -108,4 +119,5 @@ pmModelParamsFromPSF modelParamsFromPSF; pmModelFitStatusFunc modelFitStatus; + pmModelSetLimitsFunc modelSetLimits; }; @@ -151,5 +163,5 @@ pmModel *model, ///< The input pmModel pmModelOpMode mode, ///< mode to control how the model is added into the image - psImageMaskType maskVal ///< Value to mask + psImageMaskType maskVal ///< Value to mask ); @@ -169,5 +181,5 @@ pmModel *model, ///< The input pmModel pmModelOpMode mode, ///< mode to control how the model is added into the image - psImageMaskType maskVal ///< Value to mask + psImageMaskType maskVal ///< Value to mask ); @@ -202,4 +214,14 @@ ); + +/// Set the model parameter limits for the given model +/// +/// Wraps the model-specific pmModelSetLimitsFunc function. +bool pmModelSetLimits( + const pmModel *model, ///< Model of interest + pmModelLimits type ///< Type of limits + ); + + /// @} # endif /* PM_MODEL_H */ Index: trunk/psModules/src/objects/pmModelClass.c =================================================================== --- trunk/psModules/src/objects/pmModelClass.c (revision 25477) +++ trunk/psModules/src/objects/pmModelClass.c (revision 25738) @@ -40,18 +40,18 @@ double sqrt (double x); -# include "models/pmModel_GAUSS.c" -# include "models/pmModel_PGAUSS.c" -# include "models/pmModel_QGAUSS.c" -# include "models/pmModel_PS1_V1.c" -# include "models/pmModel_RGAUSS.c" -# include "models/pmModel_SERSIC.c" +# include "models/pmModel_GAUSS.h" +# include "models/pmModel_PGAUSS.h" +# include "models/pmModel_QGAUSS.h" +# include "models/pmModel_PS1_V1.h" +# include "models/pmModel_RGAUSS.h" +# include "models/pmModel_SERSIC.h" static pmModelClass defaultModels[] = { - {"PS_MODEL_GAUSS", 7, pmModelFunc_GAUSS, pmModelFlux_GAUSS, pmModelRadius_GAUSS, pmModelLimits_GAUSS, pmModelGuess_GAUSS, pmModelFromPSF_GAUSS, pmModelParamsFromPSF_GAUSS, pmModelFitStatus_GAUSS}, - {"PS_MODEL_PGAUSS", 7, pmModelFunc_PGAUSS, pmModelFlux_PGAUSS, pmModelRadius_PGAUSS, pmModelLimits_PGAUSS, pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelParamsFromPSF_PGAUSS, pmModelFitStatus_PGAUSS}, - {"PS_MODEL_QGAUSS", 8, pmModelFunc_QGAUSS, pmModelFlux_QGAUSS, pmModelRadius_QGAUSS, pmModelLimits_QGAUSS, pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelParamsFromPSF_QGAUSS, pmModelFitStatus_QGAUSS}, - {"PS_MODEL_PS1_V1", 8, pmModelFunc_PS1_V1, pmModelFlux_PS1_V1, pmModelRadius_PS1_V1, pmModelLimits_PS1_V1, pmModelGuess_PS1_V1, pmModelFromPSF_PS1_V1, pmModelParamsFromPSF_PS1_V1, pmModelFitStatus_PS1_V1}, - {"PS_MODEL_RGAUSS", 8, pmModelFunc_RGAUSS, pmModelFlux_RGAUSS, pmModelRadius_RGAUSS, pmModelLimits_RGAUSS, pmModelGuess_RGAUSS, pmModelFromPSF_RGAUSS, pmModelParamsFromPSF_RGAUSS, pmModelFitStatus_RGAUSS}, - {"PS_MODEL_SERSIC", 8, pmModelFunc_SERSIC, pmModelFlux_SERSIC, pmModelRadius_SERSIC, pmModelLimits_SERSIC, pmModelGuess_SERSIC, pmModelFromPSF_SERSIC, pmModelParamsFromPSF_SERSIC, pmModelFitStatus_SERSIC} + {"PS_MODEL_GAUSS", 7, (pmModelFunc)pmModelFunc_GAUSS, (pmModelFlux)pmModelFlux_GAUSS, (pmModelRadius)pmModelRadius_GAUSS, (pmModelLimits)pmModelLimits_GAUSS, (pmModelGuessFunc)pmModelGuess_GAUSS, (pmModelFromPSFFunc)pmModelFromPSF_GAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_GAUSS, (pmModelFitStatusFunc)pmModelFitStatus_GAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_GAUSS }, + {"PS_MODEL_PGAUSS", 7, (pmModelFunc)pmModelFunc_PGAUSS, (pmModelFlux)pmModelFlux_PGAUSS, (pmModelRadius)pmModelRadius_PGAUSS, (pmModelLimits)pmModelLimits_PGAUSS, (pmModelGuessFunc)pmModelGuess_PGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_PGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_PGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_PGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_PGAUSS }, + {"PS_MODEL_QGAUSS", 8, (pmModelFunc)pmModelFunc_QGAUSS, (pmModelFlux)pmModelFlux_QGAUSS, (pmModelRadius)pmModelRadius_QGAUSS, (pmModelLimits)pmModelLimits_QGAUSS, (pmModelGuessFunc)pmModelGuess_QGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_QGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_QGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_QGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_QGAUSS }, + {"PS_MODEL_PS1_V1", 8, (pmModelFunc)pmModelFunc_PS1_V1, (pmModelFlux)pmModelFlux_PS1_V1, (pmModelRadius)pmModelRadius_PS1_V1, (pmModelLimits)pmModelLimits_PS1_V1, (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 }, + {"PS_MODEL_RGAUSS", 8, (pmModelFunc)pmModelFunc_RGAUSS, (pmModelFlux)pmModelFlux_RGAUSS, (pmModelRadius)pmModelRadius_RGAUSS, (pmModelLimits)pmModelLimits_RGAUSS, (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS }, + {"PS_MODEL_SERSIC", 8, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC } }; @@ -168,2 +168,18 @@ return (models[type].name); } + + +void pmModelClassSetLimits(pmModelLimitsType type) +{ + if (!models) { + pmModelClassInit(); + } + + for (int i = 0; i < Nmodels; i++) { + if (models[i].modelSetLimits) { + models[i].modelSetLimits(type); + } + } + +} + Index: trunk/psModules/src/objects/pmModelClass.h =================================================================== --- trunk/psModules/src/objects/pmModelClass.h (revision 25477) +++ trunk/psModules/src/objects/pmModelClass.h (revision 25738) @@ -29,9 +29,10 @@ # define PM_MODEL_CLASS_H +#include + /// @addtogroup Objects Object Detection / Analysis Functions /// @{ -typedef struct -{ +typedef struct { char *name; int nParams; @@ -44,4 +45,5 @@ pmModelParamsFromPSF modelParamsFromPSF; pmModelFitStatusFunc modelFitStatus; + pmModelSetLimitsFunc modelSetLimits; } pmModelClass; @@ -73,4 +75,8 @@ pmModelType pmModelClassGetType (const char *name); +/// Set parameter limits for all models +void pmModelClassSetLimits(pmModelLimitsType type); + + /// @} # endif /* PM_MODEL_CLASS_H */