IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 28, 2006, 10:23:51 AM (20 years ago)
Author:
magnier
Message:

This is an important change to the PSF shape model.

  • I have changed the definition of the elliptical contour parameters

(PM_PAR_SXX,SYY are now sigma_x,y/sqrt(2), rather than
sqrt(2)/sigma_x.

  • I have changed the way in which the SXY term is modelled. rather

than fit SXY to a polynomial, I now fit SXY / (1/SXX2 + 1/SYY2)2.
this representation is more likely to be well-described by a
polynomial. the other form is likely to vary as xy/rN rather than xy
rN.

these changes imply that the models all need to treat the SXX,SYY,SXY
terms in the same fashion (ie, fitting a contour of the form (x/SXX)2
+ (y/SYY)2 + SXY*x*y). If an arbitrary model uses some other meaning
for SXX, SYY, SXY, then it will have to define its own exceptions in
the functions pmPSFFromModels and pmModelFromPSF (the latter has an
implementation for each model already).

  • I have cleaned up the pmModel_NAME.c functions to use the #define names

for the parameters (PM_PAR_XXX) everywhere, and to use #def names for
the functions. this latter change makes it easy to specify the model
function names in a single location.

Location:
trunk/psModules/src/objects/models
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/models/pmModel_GAUSS.c

    r9730 r9770  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    51/******************************************************************************
    6     params->data.F32[PM_PAR_SKY] = So;
    7     params->data.F32[PM_PAR_I0] = Zo;
    8     params->data.F32[PM_PAR_XPOS] = Xo;
    9     params->data.F32[PM_PAR_YPOS] = Yo;
    10     params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX;
    11     params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY;
    12     params->data.F32[PM_PAR_SXY] = Sxy;
    13 *****************************************************************************/
    14 
    15 psF32 pmModelFunc_GAUSS(psVector *deriv,
    16                         const psVector *params,
    17                         const psVector *x)
    18 {
    19     psF32 X  = x->data.F32[0] - params->data.F32[PM_PAR_XPOS];
    20     psF32 Y  = x->data.F32[1] - params->data.F32[PM_PAR_YPOS];
    21     psF32 px = params->data.F32[PM_PAR_SXX]*X;
    22     psF32 py = params->data.F32[PM_PAR_SYY]*Y;
    23     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[PM_PAR_SXY]*X*Y;
     2 * this file defines the PGAUSS source shape model.  Note that these model functions are loaded
     3 * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
     4 * models use a psVector to represent the set of parameters, with the sequence used to specify
     5 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
     6 * specifics of the model.  All models which are used a PSF representations share a few
     7 * parameters, for which # define names are listed in pmModel.h:
     8 
     9 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
     10 * PM_PAR_I0 1    - central intensity
     11 * PM_PAR_XPOS 2  - X center of object
     12 * PM_PAR_YPOS 3  - Y center of object
     13 * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
     14 * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     15 * PM_PAR_SXY 6   - X*Y term of elliptical contour
     16 *****************************************************************************/
     17
     18# define PM_MODEL_FUNC       pmModelFunc_GAUSS
     19# define PM_MODEL_FLUX       pmModelFlux_GAUSS
     20# define PM_MODEL_GUESS      pmModelGuess_GAUSS
     21# define PM_MODEL_LIMITS     pmModelLimits_GAUSS
     22# define PM_MODEL_RADIUS     pmModelRadius_GAUSS
     23# define PM_MODEL_FROM_PSF   pmModelFromPSF_GAUSS
     24# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS
     25
     26psF32 PM_MODEL_FUNC(psVector *deriv,
     27                    const psVector *params,
     28                    const psVector *pixcoord)
     29{
     30    psF32 *PAR = params->data.F32;
     31
     32    // XXX this is fitting sqrt(2)/sigma_x, sqrt(2)/sigma_y
     33    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
     34    psF32 Y  = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS];
     35    psF32 px = X / PAR[PM_PAR_SXX];
     36    psF32 py = Y / PAR[PM_PAR_SYY];
     37    psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    2438    psF32 r  = exp(-z);
    25     psF32 q  = params->data.F32[PM_PAR_I0]*r;
    26     psF32 f  = q + params->data.F32[PM_PAR_SKY];
     39    psF32 q  = PAR[PM_PAR_I0]*r;
     40    psF32 f  = q + PAR[PM_PAR_SKY];
    2741
    2842    if (deriv != NULL) {
    29         deriv->data.F32[PM_PAR_SKY] = +1.0;
    30         deriv->data.F32[PM_PAR_I0] = +r;
    31         deriv->data.F32[PM_PAR_XPOS] = q*(2*px*params->data.F32[PM_PAR_SXX] + params->data.F32[PM_PAR_SXY]*Y);
    32         deriv->data.F32[PM_PAR_YPOS] = q*(2*py*params->data.F32[PM_PAR_SYY] + params->data.F32[PM_PAR_SXY]*X);
    33         deriv->data.F32[PM_PAR_SXX] = -2.0*q*px*X;
    34         deriv->data.F32[PM_PAR_SYY] = -2.0*q*py*Y;
    35         deriv->data.F32[PM_PAR_SXY] = -q*X*Y;
     43        psF32 *dPAR = deriv->data.F32;
     44        dPAR[PM_PAR_SKY]  = +1.0;
     45        dPAR[PM_PAR_I0]   = +r;
     46        dPAR[PM_PAR_XPOS] = q*(2*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
     47        dPAR[PM_PAR_YPOS] = q*(2*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     48        dPAR[PM_PAR_SXX]  = +2.0*q*px*px/PAR[PM_PAR_SXX];
     49        dPAR[PM_PAR_SYY]  = +2.0*q*py*py/PAR[PM_PAR_SYY];
     50        dPAR[PM_PAR_SXY]  = -q*X*Y;
    3651    }
    3752    return(f);
    … …  
    3954
    4055// define the parameter limits
    41 bool pmModelLimits_GAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
     56bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max)
    4257{
    4358
    … …  
    5873    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
    5974    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
    60     params_min[0][0].data.F32[PM_PAR_SXX] = 0.01;
    61     params_min[0][0].data.F32[PM_PAR_SYY] = 0.01;
     75    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     76    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
    6277    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    6378
    … …  
    6681    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
    6782    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
    68     params_max[0][0].data.F32[PM_PAR_SXX] = 2.0;
    69     params_max[0][0].data.F32[PM_PAR_SYY] = 2.0;
     83    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     84    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
    7085    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    7186
    … …  
    7489
    7590// make an initial guess for parameters
    76 bool pmModelGuess_GAUSS (pmModel *model, pmSource *source)
     91bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    7792{
    7893
    7994    pmMoments *moments = source->moments;
    80     psF32     *params  = model->params->data.F32;
    81 
    82     params[PM_PAR_SKY] = moments->Sky;
    83     params[PM_PAR_I0] = moments->Peak - moments->Sky;
    84     params[PM_PAR_XPOS] = moments->x;
    85     params[PM_PAR_YPOS] = moments->y;
    86     params[PM_PAR_SXX] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
    87     params[PM_PAR_SYY] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
    88     params[PM_PAR_SXY] = 0.0;
     95    psF32     *PAR  = model->params->data.F32;
     96
     97    PAR[PM_PAR_SKY] = moments->Sky;
     98    PAR[PM_PAR_I0] = moments->Peak - moments->Sky;
     99    PAR[PM_PAR_XPOS] = moments->x;
     100    PAR[PM_PAR_YPOS] = moments->y;
     101    PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx);
     102    PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy);
     103    PAR[PM_PAR_SXY] = 0.0;
    89104    return(true);
    90105}
    91106
    92 psF64 pmModelFlux_GAUSS(const psVector *params)
    93 {
    94     psF64 A1   = PS_SQR(params->data.F32[PM_PAR_SXX]);
    95     psF64 A2   = PS_SQR(params->data.F32[PM_PAR_SYY]);
    96     psF64 A3   = PS_SQR(params->data.F32[PM_PAR_SXY]);
    97     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     107psF64 PM_MODEL_FLUX (const psVector *params)
     108{
     109
     110    psEllipseShape shape;
     111
     112    psF32 *PAR = params->data.F32;
     113
     114    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     115    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     116    shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     117
    98118    // Area is equivalent to 2 pi sigma^2
     119    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     120    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    99121
    100122    psF64 Flux = params->data.F32[PM_PAR_I0] * Area;
    … …  
    105127// return the radius which yields the requested flux
    106128// this function is never allowed to return <= 0
    107 psF64 pmModelRadius_GAUSS  (const psVector *params, psF64 flux)
    108 {
     129psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
     130{
     131    psEllipseShape shape;
    109132
    110133    if (flux <= 0)
    … …  
    115138        return (1.0);
    116139
    117     psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[PM_PAR_SXX], 1.0 / params->data.F32[PM_PAR_SYY]);
    118     psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
     140    psF32 *PAR = params->data.F32;
     141
     142    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     143    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     144    shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     145
     146    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     147    psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
    119148    return (radius);
    120149}
    121150
    122151// construct the PSF model from the FLT model and the psf
    123 bool pmModelFromPSF_GAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
     152bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
    124153{
    125154
    … …  
    127156    psF32 *in  = modelFLT->params->data.F32;
    128157
    129     out[PM_PAR_SKY] = in[PM_PAR_SKY];
    130     out[PM_PAR_I0] = in[PM_PAR_I0];
    131     out[PM_PAR_XPOS] = in[PM_PAR_XPOS];
    132     out[PM_PAR_YPOS] = in[PM_PAR_YPOS];
    133 
    134     assert (PM_PAR_YPOS + 1 == 4);  // so starting at 4 is correct
    135     for (int i = 4; i < 7; i++) {
    136         psPolynomial2D *poly = psf->params->data[i-4];
    137         out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]);
     158    // we require these two parameters to exist
     159    assert (psf->params_NEW->n > PM_PAR_YPOS);
     160    assert (psf->params_NEW->n > PM_PAR_XPOS);
     161
     162    // supply the model-fitted parameters, or copy from the input
     163    for (int i = 0; i < psf->params_NEW->n; i++) {
     164        if (psf->params_NEW->data[i] == NULL) {
     165            out[i] = in[i];
     166        } else {
     167            psPolynomial2D *poly = psf->params_NEW->data[i];
     168            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
     169        }
    138170    }
     171
     172    // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here
     173    out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
     174
    139175    return(true);
    140176}
    141177
    142178// check the status of the fitted model
    143 bool pmModelFitStatus_GAUSS (pmModel *model)
     179// this test is invalid if the parameters are derived
     180// from the PSF model
     181bool PM_MODEL_FIT_STATUS (pmModel *model)
    144182{
    145183
    … …  
    164202    return false;
    165203}
     204
     205# undef PM_MODEL_FUNC
     206# undef PM_MODEL_FLUX
     207# undef PM_MODEL_GUESS
     208# undef PM_MODEL_LIMITS
     209# undef PM_MODEL_RADIUS
     210# undef PM_MODEL_FROM_PSF
     211# undef PM_MODEL_FIT_STATUS

  • trunk/psModules/src/objects/models/pmModel_PGAUSS.c

    r9730 r9770  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    51/******************************************************************************
    6     params->data.F32[0] = So;
    7     params->data.F32[1] = Zo;
    8     params->data.F32[2] = Xo;
    9     params->data.F32[3] = Yo;
    10     params->data.F32[4] = sqrt(2.0) / SigmaX;
    11     params->data.F32[5] = sqrt(2.0) / SigmaY;
    12     params->data.F32[6] = Sxy;
    13 *****************************************************************************/
    14 
    15 psF32 pmModelFunc_PGAUSS(psVector *deriv,
    16                          const psVector *params,
    17                          const psVector *x)
     2 * this file defines the PGAUSS source shape model.  Note that these model functions are loaded
     3 * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
     4 * models use a psVector to represent the set of parameters, with the sequence used to specify
     5 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
     6 * specifics of the model.  All models which are used a PSF representations share a few
     7 * parameters, for which # define names are listed in pmModel.h:
     8 
     9 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
     10 * PM_PAR_I0 1    - central intensity
     11 * PM_PAR_XPOS 2  - X center of object
     12 * PM_PAR_YPOS 3  - Y center of object
     13 * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
     14 * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     15 * PM_PAR_SXY 6   - X*Y term of elliptical contour
     16 *****************************************************************************/
     17
     18# define PM_MODEL_FUNC       pmModelFunc_PGAUSS
     19# define PM_MODEL_FLUX       pmModelFlux_PGAUSS
     20# define PM_MODEL_GUESS      pmModelGuess_PGAUSS
     21# define PM_MODEL_LIMITS     pmModelLimits_PGAUSS
     22# define PM_MODEL_RADIUS     pmModelRadius_PGAUSS
     23# define PM_MODEL_FROM_PSF   pmModelFromPSF_PGAUSS
     24# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS
     25
     26// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
     27psF32 PM_MODEL_FUNC(psVector *deriv,
     28                    const psVector *params,
     29                    const psVector *pixcoord)
    1830{
    1931    psF32 *PAR = params->data.F32;
    2032
    21     psF32 X  = x->data.F32[0] - PAR[2];
    22     psF32 Y  = x->data.F32[1] - PAR[3];
    23     psF32 px = PAR[4]*X;
    24     psF32 py = PAR[5]*Y;
    25     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;
     33    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
     34    psF32 Y  = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS];
     35    psF32 px = X / PAR[PM_PAR_SXX];
     36    psF32 py = Y / PAR[PM_PAR_SYY];
     37    psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    2638    psF32 t  = 1 + z + z*z/2.0;
    2739    psF32 r  = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */
    28     psF32 f  = PAR[1]*r + PAR[0];
     40    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
    2941
    3042    if (deriv != NULL) {
    3143        psF32 *dPAR = deriv->data.F32;
    32         psF32 q = PAR[1]*r*r*t;
    33         dPAR[0] = +1.0;
    34         dPAR[1] = +r;
    35         dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    36         dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    37         dPAR[4] = -2.0*q*px*X;
    38         dPAR[5] = -2.0*q*py*Y;
    39         dPAR[6] = -q*X*Y;
     44        psF32 q = PAR[PM_PAR_I0]*r*r*t;
     45        dPAR[PM_PAR_SKY] = +1.0;
     46        dPAR[PM_PAR_I0] = +r;
     47        dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
     48        dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     49        dPAR[PM_PAR_SXX] =  +2.0*q*px*px/PAR[PM_PAR_SXX];
     50        dPAR[PM_PAR_SYY] =  +2.0*q*py*py/PAR[PM_PAR_SYY];
     51        dPAR[PM_PAR_SXY] = -q*X*Y;
    4052    }
    4153    return(f);
    4254}
    4355
    44 bool pmModelLimits_PGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
     56bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max)
    4557{
    4658
    … …  
    4961    *params_max = psVectorAlloc (7, PS_TYPE_F32);
    5062
    51     beta_lim[0][0].data.F32[0] = 1000;
    52     beta_lim[0][0].data.F32[1] = 3e6;
    53     beta_lim[0][0].data.F32[2] = 5;
    54     beta_lim[0][0].data.F32[3] = 5;
    55     beta_lim[0][0].data.F32[4] = 0.5;
    56     beta_lim[0][0].data.F32[5] = 0.5;
    57     beta_lim[0][0].data.F32[6] = 0.5;
    58 
    59     params_min[0][0].data.F32[0] = -1000;
    60     params_min[0][0].data.F32[1] = 0;
    61     params_min[0][0].data.F32[2] = -100;
    62     params_min[0][0].data.F32[3] = -100;
    63     params_min[0][0].data.F32[4] = 0.01;
    64     params_min[0][0].data.F32[5] = 0.01;
    65     params_min[0][0].data.F32[6] = -5.0;
    66 
    67     params_max[0][0].data.F32[0] = 1e5;
    68     params_max[0][0].data.F32[1] = 1e8;
    69     params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    70     params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    71     params_max[0][0].data.F32[4] = 2.0;
    72     params_max[0][0].data.F32[5] = 2.0;
    73     params_max[0][0].data.F32[6] = +5.0;
     63    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     64    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     65    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     66    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     67    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     68    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     69    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     70
     71    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     72    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     73    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     74    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     75    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     76    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
     77    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     78
     79    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     80    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     81    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     82    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     83    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     84    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
     85    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    7486
    7587    return (TRUE);
    … …  
    7789
    7890// make an initial guess for parameters
    79 bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source)
     91bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    8092{
    8193
    8294    pmMoments *moments = source->moments;
    83     psF32     *params  = model->params->data.F32;
    84 
    85     params[0] = moments->Sky;
    86     params[1] = moments->Peak - moments->Sky;
    87     params[2] = moments->x;
    88     params[3] = moments->y;
    89     params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
    90     params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
    91     params[6] = 0.0;
     95    psF32     *PAR     = model->params->data.F32;
     96
     97    PAR[PM_PAR_SKY] = moments->Sky;
     98    PAR[PM_PAR_I0]  = moments->Peak - moments->Sky;
     99    PAR[PM_PAR_XPOS] = moments->x;
     100    PAR[PM_PAR_YPOS] = moments->y;
     101    PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx);
     102    PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy);
     103    PAR[PM_PAR_SXY] = 0.0;
    92104
    93105    return(true);
    94106}
    95107
    96 psF64 pmModelFlux_PGAUSS(const psVector *params)
     108psF64 PM_MODEL_FLUX(const psVector *params)
    97109{
    98110    float norm, z;
     111    psEllipseShape shape;
    99112
    100113    psF32 *PAR = params->data.F32;
    101114
    102     psF64 A1   = PS_SQR(PAR[4]);
    103     psF64 A2   = PS_SQR(PAR[5]);
    104     psF64 A3   = PS_SQR(PAR[6]);
    105     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     115    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     116    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     117    shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     118
    106119    // Area is equivalent to 2 pi sigma^2
     120    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     121    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    107122
    108123    // the area needs to be multiplied by the integral of f(z)
    … …  
    122137    norm *= DZ / 3.0;
    123138
    124     psF64 Flux = PAR[1] * Area * norm;
     139    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
    125140
    126141    return(Flux);
    … …  
    129144// define this function so it never returns Inf or NaN
    130145// return the radius which yields the requested flux
    131 psF64 pmModelRadius_PGAUSS (const psVector *params, psF64 flux)
     146psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    132147{
    133148    if (flux <= 0)
    134149        return (1.0);
    135     if (params->data.F32[1] <= 0)
     150    if (params->data.F32[PM_PAR_I0] <= 0)
    136151        return (1.0);
    137     if (flux >= params->data.F32[1])
     152    if (flux >= params->data.F32[PM_PAR_I0])
    138153        return (1.0);
    139154
    140     psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
    141     psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
    142     if (isnan(radius)) {
    143         fprintf (stderr, "error in code\n");
    144     }
     155    psF32 *PAR = params->data.F32;
     156
     157    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     158    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     159    shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     160
     161    // this estimates the radius assuming f(z) is roughly exp(-z)
     162    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     163    psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
     164
     165    if (isnan(radius))
     166        psAbort ("psphot.model", "error in code: never return invalid radius");
     167    if (radius < 0)
     168        psAbort ("psphot.model", "error in code: never return invalid radius");
     169
    145170    return (radius);
    146171}
    147172
    148 bool pmModelFromPSF_PGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
     173bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
    149174{
    150175
    … …  
    152177    psF32 *in  = modelFLT->params->data.F32;
    153178
    154     out[0] = in[0];
    155     out[1] = in[1];
    156     out[2] = in[2];
    157     out[3] = in[3];
    158 
    159     for (int i = 4; i < 7; i++) {
    160         psPolynomial2D *poly = psf->params->data[i-4];
    161         out[i] = psPolynomial2DEval(poly, out[2], out[3]);
     179    // we require these two parameters to exist
     180    assert (psf->params_NEW->n > PM_PAR_YPOS);
     181    assert (psf->params_NEW->n > PM_PAR_XPOS);
     182
     183    for (int i = 0; i < psf->params_NEW->n; i++) {
     184        if (psf->params_NEW->data[i] == NULL) {
     185            out[i] = in[i];
     186        } else {
     187            psPolynomial2D *poly = psf->params_NEW->data[i];
     188            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
     189        }
    162190    }
     191
     192    // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here
     193    out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
     194
    163195    return(true);
    164196}
    165197
    166 bool pmModelFitStatus_PGAUSS (pmModel *model)
     198bool PM_MODEL_FIT_STATUS (pmModel *model)
    167199{
    168200
    … …  
    174206
    175207    dP = 0;
    176     dP += PS_SQR(dPAR[4] / PAR[4]);
    177     dP += PS_SQR(dPAR[5] / PAR[5]);
     208    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     209    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    178210    dP = sqrt (dP);
    179211
    180212    status = true;
    181213    status &= (dP < 0.5);
    182     status &= (PAR[1] > 0);
    183     status &= ((dPAR[1]/PAR[1]) < 0.5);
     214    status &= (PAR[PM_PAR_I0] > 0);
     215    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    184216
    185217    if (status)
    … …  
    187219    return false;
    188220}
     221
     222# undef PM_MODEL_FUNC
     223# undef PM_MODEL_FLUX
     224# undef PM_MODEL_GUESS
     225# undef PM_MODEL_LIMITS
     226# undef PM_MODEL_RADIUS
     227# undef PM_MODEL_FROM_PSF
     228# undef PM_MODEL_FIT_STATUS
Note: See TracChangeset for help on using the changeset viewer.