IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9770


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
Files:
8 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

  • trunk/psModules/src/objects/pmModelGroup.c

    r8815 r9770  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-09-15 09:49:01 $
     8 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-10-28 20:23:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3535double sqrt (double x);
    3636
    37 #include "models/pmModel_GAUSS.c"
    38 #include "models/pmModel_PGAUSS.c"
    39 #include "models/pmModel_QGAUSS.c"
    40 #include "models/pmModel_SGAUSS.c"
     37# include "models/pmModel_GAUSS.c"
     38// # include "models/pmModel_PGAUSS.c"
     39// # include "models/pmModel_QGAUSS.c"
     40// # include "models/pmModel_SGAUSS.c"
    4141
    4242static pmModelGroup defaultModels[] = {
    4343                                          {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
    44                                           {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
    45                                           {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
    46                                           {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
     44                                          //                                          {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
     45                                          //                                          {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
     46                                          //                                          {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
    4747                                      };
    4848
    … …  
    208208{
    209209    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
    210     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    211     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     210    PS_ASSERT_PTR_NON_NULL(source, NULL);
     211    PS_ASSERT_PTR_NON_NULL(source->moments, NULL);
     212    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
    212213
    213214    pmModel *model = pmModelAlloc(modelType);

  • trunk/psModules/src/objects/pmPSF.c

    r9730 r9770  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-10-24 22:55:05 $
     8 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-10-28 20:23:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6969    psFree (psf->ApTrend);
    7070    psFree (psf->growth);
    71     psFree (psf->params);
     71    psFree (psf->params_NEW);
    7272    return;
    7373}
    7474
    75 
    76 
    7775/*****************************************************************************
    78 pmPSFAlloc (type): allocate a pmPSF.
    79     NOTE: a PSF always has 4 parameters fewer than the equivalent model.
    80       This is because those 4 parameters are
     76 pmPSFAlloc (type): allocate a pmPSF.
     77 
     78 NOTE: PSF model parameters which are not modeled on an image are set to NULL in psf->params.
     79 
     80 These are normally:
     81 
    8182 X-center
    8283 Y-center
    … …  
    119120        return(NULL);
    120121    }
    121 
    122     psf->params = psArrayAlloc(Nparams - 4);
    123 
    124     // the order of the PSF parameter (X,Y) fits is determined by the
    125     // psfTrendMask polynomial (user-specified as in the recipe). the
    126     // masks of psfTrendMask are applied to each parameter.
    127     // if psfTrendMask is NULL, these polynomials are not allocated.
    128     // in this case, the user must set them by hand (as in pmPSFfromMD)
    129     // XXX should we drop the hard-wired '4' above and use NULL to identify
    130     // the parameters which are not fitted.  these could be selected by
    131     // testing for the value of PM_PAR_XPOS, etc.
     122    psf->params_NEW = psArrayAlloc(Nparams);
     123
     124    // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
     125    // (user-specified as in the recipe). the masks of psfTrendMask are applied to each
     126    // parameter.  if psfTrendMask is NULL, these polynomials are not allocated.  in this case,
     127    // the user must set them by hand (as in pmPSFfromMD).  the parameters which are not fitted
     128    // are left as NULL.  these are selected by testing for them by the named values (
     129    // PM_PAR_XPOS, etc)
     130
    132131    if (psfTrendMask) {
    133         for (int i = 0; i < psf->params->n; i++) {
     132        for (int i = 0; i < psf->params_NEW->n; i++) {
     133            if (i == PM_PAR_SKY)
     134                continue;
     135            if (i == PM_PAR_I0)
     136                continue;
     137            if (i == PM_PAR_XPOS)
     138                continue;
     139            if (i == PM_PAR_YPOS)
     140                continue;
     141
    134142            psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY);
    135143            for (int nx = 0; nx < param->nX + 1; nx++) {
    … …  
    138146                }
    139147            }
    140             psf->params->data[i] = param;
     148            psf->params_NEW->data[i] = param;
    141149        }
    142150    }
    … …  
    166174    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
    167175
     176    // construct the x,y terms
    168177    for (int i = 0; i < models->n; i++) {
    169178        pmModel *model = models->data[i];
    … …  
    171180            continue;
    172181
    173         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
    174         // note that the data is provided in the F64 array
    175         x->data.F64[i] = model->params->data.F32[2];
    176         y->data.F64[i] = model->params->data.F32[3];
     182        // use F64 for polynomial fitting
     183        x->data.F64[i] = model->params->data.F32[PM_PAR_XPOS];
     184        y->data.F64[i] = model->params->data.F32[PM_PAR_YPOS];
    177185    }
    178186
    … …  
    183191    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    184192
    185     for (int i = 0; i < psf->params->n; i++) {
     193    // skip the unfitted parameters (X, Y, Io, Sky)
     194    for (int i = 0; i < psf->params_NEW->n; i++) {
     195        if (i == PM_PAR_SKY)
     196            continue;
     197        if (i == PM_PAR_I0)
     198            continue;
     199        if (i == PM_PAR_XPOS)
     200            continue;
     201        if (i == PM_PAR_YPOS)
     202            continue;
     203
     204        // select the per-object fitted data for this parameter
    186205        for (int j = 0; j < models->n; j++) {
    187206            pmModel *model = models->data[j];
    188207            if (model == NULL)
    189208                continue;
    190             z->data.F64[j] = model->params->data.F32[i + 4];
    191             dz->data.F64[j] = 1;
    192             // XXX EAM : need to use actual errors?
    193             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
     209            z->data.F64[j] = model->params->data.F32[i];
     210            dz->data.F64[j] = 1; // use the model-fitted error? or S/N?
     211
     212            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
     213            if (i == PM_PAR_SXY) {
     214                z->data.F64[j] = pmPSF_SXYfromModel (model->params->data.F32);
     215            }
    194216        }
    195 
    196         // XXX EAM : this is the expected API (cycle 7? cycle 8?)
    197         psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y);
    198 
    199         // XXX EAM : drop this when above is implemented...
    200         // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
    201 
     217        psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, mask, 0xff, z, dz, x, y);
    202218        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    203         // psPolynomial2DDump (psf->params->data[i]);
     219
     220        // XXX Test output
     221        psPolynomial2D *poly = psf->params_NEW->data[i];
     222        fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
     223        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][0], poly->coeff[1][0], poly->coeff[0][1]);
     224        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][2], poly->coeff[1][1], poly->coeff[2][0]);
     225    }
     226
     227    // XXX test dump of star parameters vs position (compare with fitted values)
     228    {
     229        FILE *f = fopen ("params.dat", "w");
     230
     231        for (int j = 0; j < models->n; j++)
     232        {
     233            pmModel *model = models->data[j];
     234            if (model == NULL)
     235                continue;
     236
     237            pmModel *modelPSF = pmModelFromPSF (model, psf);
     238
     239            fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     240
     241            for (int i = 0; i < psf->params_NEW->n; i++) {
     242                if (psf->params_NEW->data[i] == NULL)
     243                    continue;
     244                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
     245            }
     246            fprintf (f, "%f %d\n", model->chisq, model->nIter);
     247        }
     248        fclose (f);
    204249    }
    205250
    … …  
    211256    return (true);
    212257}
    213 
    214 
    215258
    216259/*****************************************************************************
    … …  
    250293    }
    251294    return true;
     295}
     296
     297// the PSF models the \sigma_{xy} variation of the elliptical contour as a function of position in the image with a
     298// polynomial.  an individual object has a contour of the form (x^2/2sx^2) + (y^2/2sy^2) + sxy*x*y
     299// these are the values of the model->params.  the psf->params term for sxy is actually fitted
     300// to sxy/(sxx^-2 + syy^-2)^2
     301
     302// input: model->param, output: psf->param[PM_PAR_SXY]
     303double pmPSF_SXYfromModel (psF32 *modelPar)
     304{
     305
     306    double SXX = modelPar[PM_PAR_SXX];
     307    double SYY = modelPar[PM_PAR_SYY];
     308    double SXY = modelPar[PM_PAR_SXY];
     309
     310    double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
     311    return (par);
     312}
     313
     314// input: fitted psf->param, output: model->param[PM_PAR_SXY]
     315double pmPSF_SXYtoModel (psF32 *fittedPar)
     316{
     317
     318    double SXX = fittedPar[PM_PAR_SXX];
     319    double SYY = fittedPar[PM_PAR_SYY];
     320    double fit = fittedPar[PM_PAR_SXY];
     321
     322    double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
     323    return SXY;
    252324}
    253325

  • trunk/psModules/src/objects/pmPSF.h

    r9562 r9770  
    4242{
    4343    pmModelType type;   ///< PSF Model in use
    44     psArray *params;   ///< Model parameters (psPolynomial2D)
     44    psArray *params_NEW;   ///< Model parameters (psPolynomial2D)
     45    // XXXXX I am changing params: we will allocate elements for the
     46    // unfitted elements (So, Io, Xo, Yo) and leave them as NULL
     47    // I am using a new name to catch all refs to params with gcc
    4548    psPolynomial1D *ChiTrend;  ///< Chisq vs flux fit (correction for systematic errors)
    4649    psPolynomial4D *ApTrend;  ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst))

  • trunk/psModules/src/objects/pmPSF_IO.c

    r9563 r9770  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-10-14 00:56:13 $
     8 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-10-28 20:23:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    5050    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
    5151
    52     for (int i = 0; i < nPar - 4; i++) {
    53         psPolynomial2D *poly = psf->params->data[i];
     52    for (int i = 0; i < nPar; i++) {
     53        psPolynomial2D *poly = psf->params_NEW->data[i];
     54        if (poly == NULL)
     55            continue;
    5456        psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i);
    5557    }
    … …  
    8789        psAbort ("read PSF" , "mismatch model par count");
    8890
    89     for (int i = 0; i < nPar - 4; i++) {
     91    // un-fitted terms, not in the Metadata, are left NULL
     92    // XXX add a double-check of the expected number?
     93    for (int i = 0; i < nPar; i++) {
    9094        sprintf (keyword, "PSF_PAR%02d", i);
    9195        psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
     96        if (!status)
     97            continue;
    9298        psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
    93         psFree (psf->params->data[i]);
    94         psf->params->data[i] = poly;
     99        psFree (psf->params_NEW->data[i]);
     100        psf->params_NEW->data[i] = poly;
    95101    }
    96102    sprintf (keyword, "APTREND");

  • trunk/psModules/src/objects/pmPSFtry.c

    r9730 r9770  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-10-24 22:55:05 $
     7 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-10-28 20:23:51 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    111111        pmSource *source = psfTry->sources->data[i];
    112112        pmModel  *model  = pmSourceModelGuess (source, psfTry->psf->type);
     113        if (model == NULL) {
     114            psError(PS_ERR_UNKNOWN, false, "failed to build model");
     115            return NULL;
     116        }
    113117        x = source->peak->x;
    114118        y = source->peak->y;

  • trunk/psModules/src/objects/pmSourceIO_CMP.c

    r9653 r9770  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-10-19 21:16:49 $
     5 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-10-28 20:23:51 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3535// XXX make sure the angle in correctly translated to/from degrees
    3636// XXX we lose all information from the 'type' field
     37
     38// XXX update this file is we convert to PAR[4] : SigmaX*sqrt(2) (not 1/SigmaX)
    3739
    3840// elixir-style pseudo FITS table (header + ascii list)
    … …  
    267269            shape = psEllipseAxesToShape (axes);
    268270
    269             PAR[4] = shape.sx;
    270             PAR[5] = shape.sy;
    271             PAR[6] = shape.sxy;
     271            PAR[PM_PAR_SXX] = shape.sx;
     272            PAR[PM_PAR_SYY] = shape.sy;
     273            PAR[PM_PAR_SXY] = shape.sxy;
    272274
    273275            source->modelPSF = model;
Note: See TracChangeset for help on using the changeset viewer.