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_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.