Index: branches/eam_branches/ipp-20120627/psModules/src/objects/Makefile.am =================================================================== --- branches/eam_branches/ipp-20120627/psModules/src/objects/Makefile.am (revision 34099) +++ branches/eam_branches/ipp-20120627/psModules/src/objects/Makefile.am (revision 34133) @@ -80,5 +80,6 @@ models/pmModel_SERSIC.c \ models/pmModel_EXP.c \ - models/pmModel_DEV.c + models/pmModel_DEV.c \ + models/pmModel_TRAIL.c pkginclude_HEADERS = \ @@ -126,5 +127,6 @@ models/pmModel_SERSIC.h \ models/pmModel_EXP.h \ - models/pmModel_DEV.h + models/pmModel_DEV.h \ + models/pmModel_TRAIL.h CLEANFILES = *~ Index: branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c (revision 34133) +++ branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c (revision 34133) @@ -0,0 +1,357 @@ +/****************************************************************************** + * this file defines the TRAIL source shape model. This represents an infinitely thin + * line of length convolved with a Gaussian PSF. 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: + + pure Gaussian: + exp(-z) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - flux normalization + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_LENGTH 4 - trail length + * PM_PAR_THETA 5 - position angle + * PM_PAR_SIGMA 6 - PSF Gaussian sigma (not fitted?) + *****************************************************************************/ + +#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_TRAIL.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_TRAIL +# define PM_MODEL_FLUX pmModelFlux_TRAIL +# define PM_MODEL_GUESS pmModelGuess_TRAIL +# define PM_MODEL_LIMITS pmModelLimits_TRAIL +# define PM_MODEL_RADIUS pmModelRadius_TRAIL +# define PM_MODEL_FROM_PSF pmModelFromPSF_TRAIL +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_TRAIL +# define PM_MODEL_FIT_STATUS pmModelFitStatus_TRAIL +# define PM_MODEL_SET_LIMITS pmModelSetLimits_TRAIL + +// 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 }; + +// 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, 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) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords +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 ST = sin(PAR[PM_PAR_THETA]); + psF32 CT = cos(PAR[PM_PAR_THETA]); + + psF32 S2 = 2.0 * PS_SQR(PAR[PM_PAR_SIGMA]); + + psF32 Zp = (X*CT + Y*ST + 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + psF32 Zm = (X*CT + Y*ST - 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + + // psF32 Zp = (X*CT + Y*ST + 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + // psF32 Zm = (X*CT + Y*ST - 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + + psF32 Ep = erf(Zp); + psF32 Em = erf(Zm); + + psF32 Rxy = Y*CT - X*ST; + psF32 Gxy = exp(-Rxy*Rxy/S2); + + psF32 Pxy = Gxy * (Ep - Em); + psF32 f = Pxy * PAR[PM_PAR_I0] + PAR[PM_PAR_SKY]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + dPAR[PM_PAR_SKY] = 1.0; + dPAR[PM_PAR_I0] = Pxy; + + float dGdR = -2.0 * Rxy * Gxy / S2; // -R Gxy / (2 Sigma^2) + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dRdX = +ST; + float dRdY = -CT; + float dRdT = -Y*ST - X*CT; + + float dGdX = dGdR * dRdX; + float dGdY = dGdR * dRdY; + float dGdT = dGdR * dRdT; + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dZpdX = -CT / sqrt(S2); + float dZmdX = dZpdX; // float dZmdX = -CT / sqrt(S2); dZmdX = dZpdX + + float dZpdY = -ST / sqrt(S2); // float dZmdY = -ST / sqrt(S2); dZmdY = dZpdY + float dZmdY = dZpdY; + + float dZpdL = +0.5 / sqrt(S2); + float dZmdL = -0.5 / sqrt(S2); + + float dZpdT = (-X*ST + Y*CT) / sqrt(S2); + float dZmdT = dZpdT; // dZpdT = dZmdT + + float dEdZp = exp (-Zp*Zp) * M_2_SQRTPI; + float dEdZm = exp (-Zm*Zm) * M_2_SQRTPI; + + float dEpdX = dEdZp * dZpdX; + float dEmdX = dEdZm * dZmdX; + + float dEpdY = dEdZp * dZpdY; + float dEmdY = dEdZm * dZmdY; + + float dEpdL = dEdZp * dZpdL; + float dEmdL = dEdZm * dZmdL; + + float dEpdT = dEdZp * dZpdT; + float dEmdT = dEdZm * dZmdT; + + float dPdX = dGdX * (Ep - Em) + Gxy * (dEpdX - dEmdX); + float dPdY = dGdY * (Ep - Em) + Gxy * (dEpdY - dEmdY); + + // dGdL is 0.0 because dRdL is 0.0 + float dPdL = Gxy * (dEpdL - dEmdL); + + float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT); + + dPAR[PM_PAR_XPOS] = PAR[PM_PAR_I0] * dPdX; + dPAR[PM_PAR_YPOS] = PAR[PM_PAR_I0] * dPdY; + + dPAR[PM_PAR_LENGTH] = PAR[PM_PAR_I0] * dPdL; + dPAR[PM_PAR_THETA] = PAR[PM_PAR_I0] * dPdT; + dPAR[PM_PAR_SIGMA] = NAN; // we don't actually allow this to vary, so we do not need to calculate it + } + return(f); +} + +// 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) +{ + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 + PAR[PM_PAR_SKY] = 0.0; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + psF32 *PAR = params->data.F32; + psF64 Flux = PAR[PM_PAR_I0] * PAR[PM_PAR_LENGTH] * PAR[PM_PAR_SIGMA] * 2.0 * sqrt(2.0 * M_PI); + return(Flux); +} + +// return the radius which yields the requested flux +// this function is never allowed to return <= 0 +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF32 *PAR = params->data.F32; + + return (0.5*PAR[PM_PAR_LENGTH]); +} + +// construct the PSF model from the FLT model and the psf +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) +{ + psWarning ("do you really want to use a trail as a PSF model??"); + return false; +} + +// generate a model based on a the psf model +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psWarning ("do you really want to use a trail as a PSF model??"); + return false; +} + +// check the status of the fitted model +// 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) +{ + 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: branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.h =================================================================== --- branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.h (revision 34133) +++ branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.h (revision 34133) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_TRAIL_H + +#include "pmModel.h" + +psF32 pmModelFunc_TRAIL(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_TRAIL(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_TRAIL(pmModel *model, pmSource *source); +psF64 pmModelFlux_TRAIL(const psVector *params); +psF64 pmModelRadius_TRAIL(const psVector *params, psF64 flux); +bool pmModelFromPSF_TRAIL(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_TRAIL(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_TRAIL(pmModel *model); +void pmModelSetLimits_TRAIL(pmModelLimitsType type); + +#endif Index: branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelClass.c =================================================================== --- branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelClass.c (revision 34099) +++ branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelClass.c (revision 34133) @@ -51,4 +51,5 @@ # include "models/pmModel_EXP.h" # include "models/pmModel_DEV.h" +# include "models/pmModel_TRAIL.h" static pmModelClass defaultModels[] = { @@ -59,6 +60,7 @@ {"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 }, - {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, - {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, + {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, + {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, + {"PS_MODEL_TRAIL", 7, (pmModelFunc)pmModelFunc_TRAIL, (pmModelFlux)pmModelFlux_TRAIL, (pmModelRadius)pmModelRadius_TRAIL, (pmModelLimits)pmModelLimits_TRAIL, (pmModelGuessFunc)pmModelGuess_TRAIL, (pmModelFromPSFFunc)pmModelFromPSF_TRAIL, (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL, (pmModelFitStatusFunc)pmModelFitStatus_TRAIL, (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL }, }; Index: branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelFuncs.h =================================================================== --- branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelFuncs.h (revision 34099) +++ branches/eam_branches/ipp-20120627/psModules/src/objects/pmModelFuncs.h (revision 34133) @@ -81,4 +81,9 @@ #define PM_PAR_8 8 ///< Model-dependent parameter +// these are used by pmModel_TRAIL, with refers to L and Theta explicitly +#define PM_PAR_LENGTH 4 ///< trail length +#define PM_PAR_THETA 5 ///< position angle +#define PM_PAR_SIGMA 6 ///< position angle + /*** these prototype classes are used to define elements of the pmModelClass structure below ***/