Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.c =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.c (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.c (revision 35646) @@ -118,5 +118,66 @@ } -bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) { +bool pmModelUseReff (pmModelType type) { + bool useReff = false; + useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + return useReff; +} + +// this function and the one below handle the two cases, where the model shape is uses R_eff or Sigma +bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff) { + + psEllipseShape shape = psEllipseAxesToShape (axes); + + if (!isfinite(shape.sx)) return false; + if (!isfinite(shape.sy)) return false; + if (!isfinite(shape.sxy)) return false; + + // set the shape parameters + if (useReff) { + *Sxx = PS_MAX(0.5, shape.sx); + *Syy = PS_MAX(0.5, shape.sy); + *Sxy = shape.sxy * 2.0; + } else { + *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); + *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + *Sxy = shape.sxy; + } + + return true; +} + +bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff) { + + psEllipseShape shape; + + // set the shape parameters + if (useReff) { + shape.sx = Sxx; + shape.sy = Syy; + shape.sxy = Sxy / 2.0; + } else { + shape.sx = Sxx / M_SQRT2; + shape.sy = Syy / M_SQRT2; + shape.sxy = Sxy; + } + + if ((shape.sx == 0) || (shape.sy == 0)) { + axes->major = 0.0; + axes->minor = 0.0; + axes->theta = 0.0; + } else { + // axes ratio < 20 + // replace with maxAR argument? + *axes = psEllipseShapeToAxes (shape, 20.0); + } + + return true; +} + +// Reff says if this is a model which uses R_eff (like exp or dev) instead of Sigma +// set the parameter values SXX, SXY, SYY +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff) { psEllipseMoments emoments; @@ -137,14 +198,5 @@ axes.minor *= scale; - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - - // set the shape parameters - *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); - *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); - *Sxy = shape.sxy; + pmModelAxesToParams (Sxx, Sxy, Syy, axes, useReff); return true; Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.h =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.h (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmModelUtils.h (revision 35646) @@ -44,5 +44,9 @@ bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source); bool pmModelSetNorm (float *Io, pmSource *source); -bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments); +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff); + +bool pmModelUseReff (pmModelType type); +bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff); +bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff); // XXX void pmModelSetModelVarOption (bool option); Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmPCMdata.c =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmPCMdata.c (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmPCMdata.c (revision 35646) @@ -291,13 +291,10 @@ psAssert (modelPSF, "psf model must be defined"); - psEllipseShape shape; psEllipseAxes axes; - - shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; - shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; - shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; - axes = psEllipseShapeToAxes (shape, 20.0); + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); @@ -451,13 +448,10 @@ psAssert (modelPSF, "psf model must be defined"); - psEllipseShape shape; psEllipseAxes axes; - - shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; - shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; - shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; - axes = psEllipseShapeToAxes (shape, 20.0); + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.c =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.c (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.c (revision 35646) @@ -249,39 +249,43 @@ // Mxy = SXY * (SXX^-4 + SYY^-4 - 2 SXY ^2) +// XXX deprecated // input: model->param, output: psf->param[PM_PAR_SXY] -double pmPSF_SXYfromModel (psF32 *modelPar) -{ - PS_ASSERT_PTR_NON_NULL(modelPar, NAN); - - double SXX = modelPar[PM_PAR_SXX]; - double SYY = modelPar[PM_PAR_SYY]; - double SXY = modelPar[PM_PAR_SXY]; - - double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); - return (par); -} - +// XXX double pmPSF_SXYfromModel (psF32 *modelPar) +// XXX { +// XXX PS_ASSERT_PTR_NON_NULL(modelPar, NAN); +// XXX +// XXX double SXX = modelPar[PM_PAR_SXX]; +// XXX double SYY = modelPar[PM_PAR_SYY]; +// XXX double SXY = modelPar[PM_PAR_SXY]; +// XXX +// XXX double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); +// XXX return (par); +// XXX } + +// XXX deprecated // input: fitted psf->param, output: model->param[PM_PAR_SXY] -double pmPSF_SXYtoModel (psF32 *fittedPar) -{ - PS_ASSERT_PTR_NON_NULL(fittedPar, NAN); - - double SXX = fittedPar[PM_PAR_SXX]; - double SYY = fittedPar[PM_PAR_SYY]; - double fit = fittedPar[PM_PAR_SXY]; - - double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); - - assert (!isnan(SXY)); - - return SXY; -} - -// New Concept: the PSF modelling function fits the polarization terms e0, e1, e2: - -// convert the parameters used in the fitted source model -// to the parameters used in the 2D PSF model -// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) -bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis) +// XXX double pmPSF_SXYtoModel (psF32 *fittedPar) +// XXX { +// XXX PS_ASSERT_PTR_NON_NULL(fittedPar, NAN); +// XXX +// XXX double SXX = fittedPar[PM_PAR_SXX]; +// XXX double SYY = fittedPar[PM_PAR_SYY]; +// XXX double fit = fittedPar[PM_PAR_SXY]; +// XXX +// XXX double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); +// XXX +// XXX assert (!isnan(SXY)); +// XXX +// XXX return SXY; +// XXX } + +// The PSF modelling function fits the polarization terms e0, e1, e2: + +// the FIT is the 2D representation of the shape using polarization parameters for the elliptical contour +// the MODEL is the realized psf model for a given location + +// convert the parameters (in situ) used in the fitted source model to the parameters used in +// the 2D PSF model +bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff) { PS_ASSERT_PTR_NON_NULL(fittedPar, false); @@ -298,17 +302,12 @@ return false; } - psEllipseShape shape = psEllipseAxesToShape (axes); - - fittedPar[PM_PAR_SXX] = shape.sx * M_SQRT2; - fittedPar[PM_PAR_SYY] = shape.sy * M_SQRT2; - fittedPar[PM_PAR_SXY] = shape.sxy; - + + pmModelAxesToParams (&fittedPar[PM_PAR_SXX], &fittedPar[PM_PAR_SXX], &fittedPar[PM_PAR_SXX], axes, useReff); return true; } -// convert the PSF parameters used in the 2D PSF model fit into the -// parameters used in the source model -// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) -psEllipsePol pmPSF_ModelToFit (psF32 *modelPar) +// convert the parameters (in situ) used in the 2D PSF model fit into the parameters used in +// the source model +psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff) { // must assert non-NULL input parameter @@ -319,11 +318,8 @@ PS_ASSERT_PTR_NON_NULL(modelPar, pol); - psEllipseShape shape; - - shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = modelPar[PM_PAR_SXY]; - - pol = psEllipseShapeToPol (shape); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff); + + pol = psEllipseAxesToPol (axes); return pol; @@ -332,40 +328,15 @@ // convert the parameters used in the fitted source model to the psEllipseAxes representation // (major,minor,theta) -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type) -{ - psEllipseShape shape; +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type) +{ psEllipseAxes axes; axes.major = NAN; axes.minor = NAN; axes.theta = NAN; - // XXX: must assert non-NULL input parameter + PS_ASSERT_PTR_NON_NULL(modelPar, axes); - bool useReff = false; - useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); - - if (useReff) { - shape.sx = modelPar[PM_PAR_SXX]; - shape.sy = modelPar[PM_PAR_SYY]; - shape.sxy = modelPar[PM_PAR_SXY] / 2.0; - // XXX I *think* dividing by 2.0 is the right direction, but this - // needs to be checked with a real test - } else { - shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = modelPar[PM_PAR_SXY]; - } - - if ((shape.sx == 0) || (shape.sy == 0)) { - axes.major = 0.0; - axes.minor = 0.0; - axes.theta = 0.0; - } else { - // XXX this is not really consistent with the model fit range above - axes = psEllipseShapeToAxes (shape, maxAR); - } - + bool useReff = pmModelUseReff (type); + pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff); return axes; } @@ -377,27 +348,14 @@ PS_ASSERT_PTR_NON_NULL(modelPar, false); + modelPar[PM_PAR_SXX] = 0.0; + modelPar[PM_PAR_SYY] = 0.0; + modelPar[PM_PAR_SXY] = 0.0; + if ((axes.major <= 0) || (axes.minor <= 0)) { - modelPar[PM_PAR_SXX] = 0.0; - modelPar[PM_PAR_SYY] = 0.0; - modelPar[PM_PAR_SXY] = 0.0; return true; } - - psEllipseShape shape = psEllipseAxesToShape (axes); - - bool useReff = false; - useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC")); - useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV")); - useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP")); - - if (useReff) { - modelPar[PM_PAR_SXX] = shape.sx; - modelPar[PM_PAR_SYY] = shape.sy; - modelPar[PM_PAR_SXY] = shape.sxy * 2.0; // XXX NEED factor of 2 here for correct angle conversion - } else { - modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2; - modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2; - modelPar[PM_PAR_SXY] = shape.sxy; - } + + bool useReff = pmModelUseReff (type); + pmModelAxesToParams (&modelPar[PM_PAR_SXX], &modelPar[PM_PAR_SXY], &modelPar[PM_PAR_SYY], axes, useReff); return true; } @@ -423,5 +381,6 @@ par->data.F32[PM_PAR_SXY] = sxy; - psEllipsePol pol = pmPSF_ModelToFit(par->data.F32); + bool useReff = pmModelUseReff (options->type); + psEllipsePol pol = pmPSF_ModelToFit(par->data.F32, useReff); pmTrend2D *trend = NULL; Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.h =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.h (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSF.h (revision 35646) @@ -107,8 +107,8 @@ bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type); -bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); +bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff); -psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type); +psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff); +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type); /// Calculate FWHM value from a PSF Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryMakePSF.c (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryMakePSF.c (revision 35646) @@ -214,5 +214,6 @@ assert (source->modelEXT); // all unmasked sources should have modelEXT - psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); + bool useReff = pmModelUseReff (source->modelEXT->type); + psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32, useReff); e0->data.F32[i] = pol.e0; Index: branches/eam_branches/ipp-20130509/psModules/src/objects/pmSource.c =================================================================== --- branches/eam_branches/ipp-20130509/psModules/src/objects/pmSource.c (revision 35645) +++ branches/eam_branches/ipp-20130509/psModules/src/objects/pmSource.c (revision 35646) @@ -1145,5 +1145,4 @@ bool status; psEllipseShape oldshape; - psEllipseShape newshape; psEllipseAxes axes; @@ -1166,13 +1165,13 @@ if (!isfinite(oldI0)) return false; + bool useReff = pmModelUseReff (model->type); + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); + // increase size and height of source - axes = psEllipseShapeToAxes (oldshape, 20.0); axes.major *= SIZE; axes.minor *= SIZE; - newshape = psEllipseAxesToShape (axes); + + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, useReff); PAR[PM_PAR_I0] = FACTOR*oldI0; - PAR[PM_PAR_SXX] = newshape.sx; - PAR[PM_PAR_SYY] = newshape.sy; - PAR[PM_PAR_SXY] = newshape.sxy; psImage *target = source->variance;