Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 14345) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 14652) @@ -19,12 +19,14 @@ *****************************************************************************/ -# define PM_MODEL_FUNC pmModelFunc_GAUSS -# define PM_MODEL_FLUX pmModelFlux_GAUSS -# define PM_MODEL_GUESS pmModelGuess_GAUSS -# define PM_MODEL_LIMITS pmModelLimits_GAUSS -# define PM_MODEL_RADIUS pmModelRadius_GAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_GAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS - +# define PM_MODEL_FUNC pmModelFunc_GAUSS +# define PM_MODEL_FLUX pmModelFlux_GAUSS +# define PM_MODEL_GUESS pmModelGuess_GAUSS +# define PM_MODEL_LIMITS pmModelLimits_GAUSS +# define PM_MODEL_RADIUS pmModelRadius_GAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_GAUSS +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_GAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) psF32 PM_MODEL_FUNC(psVector *deriv, const psVector *params, @@ -38,4 +40,6 @@ psF32 py = Y / PAR[PM_PAR_SYY]; psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; + assert (z >= 0.0); + psF32 r = exp(-z); psF32 q = PAR[PM_PAR_I0]*r; @@ -61,4 +65,5 @@ # define AR_MAX 20.0 # define AR_RATIO 0.99 + bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) { @@ -93,8 +98,8 @@ break; case PM_PAR_SXX: - beta_lim = 0.5; + beta_lim = 2.0; break; case PM_PAR_SYY: - beta_lim = 0.5; + beta_lim = 2.0; break; case PM_PAR_SXY: @@ -257,18 +262,17 @@ bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { - psF32 *out = modelPSF->params->data.F32; psF32 *in = modelFLT->params->data.F32; // we require these two parameters to exist - assert (psf->params_NEW->n > PM_PAR_YPOS); - assert (psf->params_NEW->n > PM_PAR_XPOS); + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); // supply the model-fitted parameters, or copy from the input - for (int i = 0; i < psf->params_NEW->n; i++) { - if (psf->params_NEW->data[i] == NULL) { + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { out[i] = in[i]; } else { - psPolynomial2D *poly = psf->params_NEW->data[i]; + psPolynomial2D *poly = psf->params->data[i]; out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } @@ -289,7 +293,7 @@ // apply the model limits here: this truncates excessive extrapolation // XXX do we need to do this still? should we put in asserts to test? - for (int i = 0; i < psf->params_NEW->n; i++) { + for (int i = 0; i < psf->params->n; i++) { // apply the limits to all components or just the psf-model parameters? - if (psf->params_NEW->data[i] == NULL) + if (psf->params->data[i] == NULL) continue; @@ -306,4 +310,56 @@ } +// 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, 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; + psPolynomial2D *poly = psf->params->data[i]; + assert (poly); + PAR[i] = psPolynomial2DEval(poly, 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)) { + psError(PM_ERR_PSF, false, "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); +} + // check the status of the fitted model // this test is invalid if the parameters are derived @@ -311,5 +367,4 @@ bool PM_MODEL_FIT_STATUS (pmModel *model) { - psF32 dP; bool status; @@ -339,3 +394,4 @@ # undef PM_MODEL_RADIUS # undef PM_MODEL_FROM_PSF +# undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 14345) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 14652) @@ -19,11 +19,12 @@ *****************************************************************************/ -# define PM_MODEL_FUNC pmModelFunc_PGAUSS -# define PM_MODEL_FLUX pmModelFlux_PGAUSS -# define PM_MODEL_GUESS pmModelGuess_PGAUSS -# define PM_MODEL_LIMITS pmModelLimits_PGAUSS -# define PM_MODEL_RADIUS pmModelRadius_PGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_PGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS +# define PM_MODEL_FUNC pmModelFunc_PGAUSS +# define PM_MODEL_FLUX pmModelFlux_PGAUSS +# define PM_MODEL_GUESS pmModelGuess_PGAUSS +# define PM_MODEL_LIMITS pmModelLimits_PGAUSS +# define PM_MODEL_RADIUS pmModelRadius_PGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_PGAUSS +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) @@ -66,12 +67,4 @@ # define AR_RATIO 0.99 -// we constraint limits by the min valid minor axis: -# define MIN_MINOR_AXIS 0.5 - -// f3 = (s_b^-2 - s_a^-2); F3_SQ_MAX is MIN_MINOR_AXIS^-4 -# define F3_SQ_MAX 16.0 - -static float saveParams[8]; - bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) { @@ -150,5 +143,4 @@ psAbort("invalid parameter %d for param min test", nParam); } - saveParams[nParam] = params[nParam]; if (params[nParam] < params_min) { params[nParam] = params_min; @@ -184,5 +176,4 @@ psAbort("invalid parameter %d for param max test", nParam); } - saveParams[nParam] = params[nParam]; if (params[nParam] > params_max) { params[nParam] = params_max; @@ -263,4 +254,6 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { + psF64 z, f; + int Nstep = 0; psEllipseShape shape; @@ -280,5 +273,36 @@ // this estimates the radius assuming f(z) is roughly exp(-z) psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); + psF64 sigma = axes.major; + + psF64 limit = flux / PAR[PM_PAR_I0]; + + // use the fact that f is monotonically decreasing + z = 0; + Nstep = 0; + + // choose a z value guaranteed to be beyond our limit + float z0 = pow((1.0 / limit), (1.0 / 3.0)); + float z1 = (1.0 / limit); + z1 = PS_MAX (z0, z1); + z0 = 0.0; + + // perform a type of bisection to find the value + float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0); + float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0); + while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) { + z = 0.5*(z0 + z1); + f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + if (f > limit) { + z0 = z; + f0 = f; + } else { + z1 = z; + f1 = f; + } + Nstep ++; + } + psF64 radius = sigma * sqrt (2.0 * z); + + // psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); if (isnan(radius)) @@ -297,12 +321,12 @@ // we require these two parameters to exist - assert (psf->params_NEW->n > PM_PAR_YPOS); - assert (psf->params_NEW->n > PM_PAR_XPOS); - - for (int i = 0; i < psf->params_NEW->n; i++) { - if (psf->params_NEW->data[i] == NULL) { + 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 { - psPolynomial2D *poly = psf->params_NEW->data[i]; + psPolynomial2D *poly = psf->params->data[i]; out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } @@ -322,7 +346,7 @@ // apply the model limits here: this truncates excessive extrapolation // XXX do we need to do this still? should we put in asserts to test? - for (int i = 0; i < psf->params_NEW->n; i++) { + for (int i = 0; i < psf->params->n; i++) { // apply the limits to all components or just the psf-model parameters? - if (psf->params_NEW->data[i] == NULL) + if (psf->params->data[i] == NULL) continue; @@ -339,4 +363,56 @@ } +// 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, 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; + psPolynomial2D *poly = psf->params->data[i]; + assert (poly); + PAR[i] = psPolynomial2DEval(poly, 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)) { + psError(PM_ERR_PSF, false, "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) { @@ -367,3 +443,4 @@ # undef PM_MODEL_RADIUS # undef PM_MODEL_FROM_PSF +# undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 14345) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 14652) @@ -20,11 +20,12 @@ *****************************************************************************/ -# define PM_MODEL_FUNC pmModelFunc_QGAUSS -# define PM_MODEL_FLUX pmModelFlux_QGAUSS -# define PM_MODEL_GUESS pmModelGuess_QGAUSS -# define PM_MODEL_LIMITS pmModelLimits_QGAUSS -# define PM_MODEL_RADIUS pmModelRadius_QGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_QGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS +# define PM_MODEL_FUNC pmModelFunc_QGAUSS +# define PM_MODEL_FLUX pmModelFlux_QGAUSS +# define PM_MODEL_GUESS pmModelGuess_QGAUSS +# define PM_MODEL_LIMITS pmModelLimits_QGAUSS +# define PM_MODEL_RADIUS pmModelRadius_QGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_QGAUSS +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS psF32 PM_MODEL_FUNC (psVector *deriv, @@ -350,12 +351,12 @@ // we require these two parameters to exist - assert (psf->params_NEW->n > PM_PAR_YPOS); - assert (psf->params_NEW->n > PM_PAR_XPOS); - - for (int i = 0; i < psf->params_NEW->n; i++) { - if (psf->params_NEW->data[i] == NULL) { + 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 { - psPolynomial2D *poly = psf->params_NEW->data[i]; + psPolynomial2D *poly = psf->params->data[i]; out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } @@ -372,7 +373,7 @@ // apply the model limits here: this truncates excessive extrapolation // XXX do we need to do this still? should we put in asserts to test? - for (int i = 0; i < psf->params_NEW->n; i++) { + for (int i = 0; i < psf->params->n; i++) { // apply the limits to all components or just the psf-model parameters? - if (psf->params_NEW->data[i] == NULL) + if (psf->params->data[i] == NULL) continue; @@ -390,4 +391,56 @@ } +// 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, 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; + psPolynomial2D *poly = psf->params->data[i]; + assert (poly); + PAR[i] = psPolynomial2DEval(poly, 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)) { + psError(PM_ERR_PSF, false, "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) { @@ -418,3 +471,4 @@ # undef PM_MODEL_RADIUS # undef PM_MODEL_FROM_PSF +# undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 14345) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 14652) @@ -20,11 +20,12 @@ *****************************************************************************/ -# define PM_MODEL_FUNC pmModelFunc_RGAUSS -# define PM_MODEL_FLUX pmModelFlux_RGAUSS -# define PM_MODEL_GUESS pmModelGuess_RGAUSS -# define PM_MODEL_LIMITS pmModelLimits_RGAUSS -# define PM_MODEL_RADIUS pmModelRadius_RGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_RGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS +# define PM_MODEL_FUNC pmModelFunc_RGAUSS +# define PM_MODEL_FLUX pmModelFlux_RGAUSS +# define PM_MODEL_GUESS pmModelGuess_RGAUSS +# define PM_MODEL_LIMITS pmModelLimits_RGAUSS +# define PM_MODEL_RADIUS pmModelRadius_RGAUSS +# define PM_MODEL_FROM_PSF pmModelFromPSF_RGAUSS +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS +# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS psF32 PM_MODEL_FUNC (psVector *deriv, @@ -343,12 +344,12 @@ // we require these two parameters to exist - assert (psf->params_NEW->n > PM_PAR_YPOS); - assert (psf->params_NEW->n > PM_PAR_XPOS); - - for (int i = 0; i < psf->params_NEW->n; i++) { - if (psf->params_NEW->data[i] == NULL) { + 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 { - psPolynomial2D *poly = psf->params_NEW->data[i]; + psPolynomial2D *poly = psf->params->data[i]; out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } @@ -365,7 +366,7 @@ // apply the model limits here: this truncates excessive extrapolation // XXX do we need to do this still? should we put in asserts to test? - for (int i = 0; i < psf->params_NEW->n; i++) { + for (int i = 0; i < psf->params->n; i++) { // apply the limits to all components or just the psf-model parameters? - if (psf->params_NEW->data[i] == NULL) + if (psf->params->data[i] == NULL) continue; @@ -383,4 +384,56 @@ } +// 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, 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; + psPolynomial2D *poly = psf->params->data[i]; + assert (poly); + PAR[i] = psPolynomial2DEval(poly, 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)) { + psError(PM_ERR_PSF, false, "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) { @@ -411,3 +464,4 @@ # undef PM_MODEL_RADIUS # undef PM_MODEL_FROM_PSF +# undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 14345) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 14652) @@ -23,11 +23,12 @@ *****************************************************************************/ -# define PM_MODEL_FUNC pmModelFunc_SERSIC -# define PM_MODEL_FLUX pmModelFlux_SERSIC -# define PM_MODEL_GUESS pmModelGuess_SERSIC -# define PM_MODEL_LIMITS pmModelLimits_SERSIC -# define PM_MODEL_RADIUS pmModelRadius_SERSIC -# define PM_MODEL_FROM_PSF pmModelFromPSF_SERSIC -# define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC +# define PM_MODEL_FUNC pmModelFunc_SERSIC +# define PM_MODEL_FLUX pmModelFlux_SERSIC +# define PM_MODEL_GUESS pmModelGuess_SERSIC +# define PM_MODEL_LIMITS pmModelLimits_SERSIC +# define PM_MODEL_RADIUS pmModelRadius_SERSIC +# define PM_MODEL_FROM_PSF pmModelFromPSF_SERSIC +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC +# define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC psF32 PM_MODEL_FUNC (psVector *deriv, @@ -360,12 +361,12 @@ // we require these two parameters to exist - assert (psf->params_NEW->n > PM_PAR_YPOS); - assert (psf->params_NEW->n > PM_PAR_XPOS); - - for (int i = 0; i < psf->params_NEW->n; i++) { - if (psf->params_NEW->data[i] == NULL) { + 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 { - psPolynomial2D *poly = psf->params_NEW->data[i]; + psPolynomial2D *poly = psf->params->data[i]; out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); } @@ -382,7 +383,7 @@ // apply the model limits here: this truncates excessive extrapolation // XXX do we need to do this still? should we put in asserts to test? - for (int i = 0; i < psf->params_NEW->n; i++) { + for (int i = 0; i < psf->params->n; i++) { // apply the limits to all components or just the psf-model parameters? - if (psf->params_NEW->data[i] == NULL) + if (psf->params->data[i] == NULL) continue; @@ -400,4 +401,56 @@ } +// 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, 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; + psPolynomial2D *poly = psf->params->data[i]; + assert (poly); + PAR[i] = psPolynomial2DEval(poly, 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)) { + psError(PM_ERR_PSF, false, "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) { @@ -428,3 +481,4 @@ # undef PM_MODEL_RADIUS # undef PM_MODEL_FROM_PSF +# undef PM_MODEL_PARAMS_FROM_PSF # undef PM_MODEL_FIT_STATUS