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.

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.