IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 29, 2006, 5:02:59 PM (20 years ago)
Author:
magnier
Message:

major cleanup of model code, fixes to use new param defs

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9621 r9775  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    51/******************************************************************************
    6     one component, two slopes:
    7     1 / (1 + z + z^Npow)
     2 * this file defines the RGAUSS source shape model (XXX need a better name!).  Note that these
     3 * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
     4 * statements of their own.  The models use a psVector to represent the set of parameters, with
     5 * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
     6 * may thus vary depending on the specifics of the model.  All models which are used a PSF
     7 * representations share a few parameters, for which # define names are listed in pmModel.h:
    88 
    9     params->data.F32[0] = So;
    10     params->data.F32[1] = Zo;
    11     params->data.F32[2] = Xo;
    12     params->data.F32[3] = Yo;
    13     params->data.F32[4] = 1 / SigmaX;
    14     params->data.F32[5] = 1 / SigmaY;
    15     params->data.F32[6] = Sxy;
    16     params->data.F32[7] = Npow
    17 *****************************************************************************/
    18 
    19 psF64 psModelFunc_RGAUSS(psVector *deriv,
    20                          const psVector *params,
    21                          const psVector *x)
     9   power-law with fitted slope
     10   1 / (1 + z + z^alpha)
     11 
     12 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
     13 * PM_PAR_I0 1    - central intensity
     14 * PM_PAR_XPOS 2  - X center of object
     15 * PM_PAR_YPOS 3  - Y center of object
     16 * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
     17 * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     18 * PM_PAR_SXY 6   - X*Y term of elliptical contour
     19 * PM_PAR_7   7   - power-law slope (alpha)
     20 *****************************************************************************/
     21
     22# define PM_MODEL_FUNC       pmModelFunc_RGAUSS
     23# define PM_MODEL_FLUX       pmModelFlux_RGAUSS
     24# define PM_MODEL_GUESS      pmModelGuess_RGAUSS
     25# define PM_MODEL_LIMITS     pmModelLimits_RGAUSS
     26# define PM_MODEL_RADIUS     pmModelRadius_RGAUSS
     27# define PM_MODEL_FROM_PSF   pmModelFromPSF_RGAUSS
     28# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS
     29
     30psF64 PS_MODEL_FUNC (psVector *deriv,
     31                     const psVector *params,
     32                     const psVector *pixcoord)
    2233{
    2334    psF32 *PAR = params->data.F32;
    2435
    25     psF32 X  = x->data.F32[0] - PAR[2];
    26     psF32 Y  = x->data.F32[1] - PAR[3];
    27     psF32 px = PAR[4]*X;
    28     psF32 py = PAR[5]*Y;
    29     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;
    30 
    31     // psF32 FACT = 1 + 5*exp(-5*PAR[7]);
    32 
    33     psF32 p  = pow(z, PAR[7] - 1.0);
     36    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
     37    psF32 Y  = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS];
     38    psF32 px = X / PAR[PM_PAR_SXX];
     39    psF32 py = Y / PAR[PM_PAR_SYY];
     40    psF32 z  = PS_SQR(px) + PS_SQR(py) + X*Y*PAR[PM_PAR_SXY];
     41
     42    psF32 p  = pow(z, PAR[PM_PAR_7] - 1.0);
    3443    psF32 r  = 1.0 / (1 + z + z*p);
    35     psF32 f  = PAR[1]*r + PAR[0];
     44    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
    3645
    3746    if (deriv != NULL) {
    38         // note difference from a pure gaussian: q = params->data.F32[1]*r
    39         psF32 t = PAR[1]*r*r;
    40         psF32 q = t*(1 + PAR[7]*p);
    41 
    42         deriv->data.F32[0] = +1.0;
    43         deriv->data.F32[1] = +r;
    44         deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    45         deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    46         deriv->data.F32[4] = -q*px*X*2;
    47         deriv->data.F32[5] = -q*py*Y*2;
    48         deriv->data.F32[6] = -q*X*Y;
    49         deriv->data.F32[7] = -5.0*t*log(z)*p*z;
    50 
    51         // deriv->data.F32[4] = -1.8*q*px*X*2;
    52         // deriv->data.F32[5] = -1.8*q*py*Y*2;
    53         // deriv->data.F32[6] = -1.5*q*X*Y;
    54         // deriv->data.F32[7] = -5.0*t*log(z)*p*z;
     47        psF32 *dPAR = deriv->data.F32;
     48
     49        // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
     50        psF32 t = PAR[PM_PAR_I0]*r*r;
     51        psF32 q = t*(1 + PAR[PM_PAR_7]*p);
     52
     53        dPAR[PM_PAR_SKY] = +1.0;
     54        dPAR[PM_PAR_I0] = +r;
     55        dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
     56        dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     57        dPAR[PM_PAR_SXX] = +2.0*q*px*px/PAR[PM_PAR_SXX];
     58        dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY];
     59        dPAR[PM_PAR_SXY] = -q*X*Y;
     60        dPAR[PM_PAR_7] = -5.0*t*log(z)*p*z;
    5561    }
    5662    return(f);
    5763}
    5864
    59 psF64 psModelFlux_RGAUSS(const psVector *params)
    60 {
    61     float f, norm, z;
     65bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
     66{
     67
     68    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
     69    *params_min = psVectorAlloc (8, PS_TYPE_F32);
     70    *params_max = psVectorAlloc (8, PS_TYPE_F32);
     71
     72    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     73    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     74    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     75    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     76    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     77    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     78    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     79    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
     80
     81    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     82    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     83    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     84    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     85    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     86    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
     87    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     88    params_min[0][0].data.F32[PM_PAR_7] = 1.25;
     89
     90    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     91    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     92    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     93    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     94    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     95    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
     96    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
     97    params_max[0][0].data.F32[PM_PAR_7] = 4.0;
     98
     99    return (TRUE);
     100}
     101
     102bool PS_MODEL_GUESS  (psModel *model, psSource *source)
     103{
     104    pmMoments *moments = source->moments;
     105    pmPeak    *peak    = source->peak;
     106    psF32     *PAR     = model->params->data.F32;
     107
     108    psEllipseAxes axes;
     109    psEllipseShape shape;
     110    psEllipseMoments tmpMoments;
     111
     112    // XXX fix this stuff : should be using correct ellipse relationships...
     113    tmpMoments.x2 = PS_SQR(source->moments->Sx);
     114    tmpMoments.y2 = PS_SQR(source->moments->Sy);
     115    tmpMoments.xy = source->moments->Sxy;
     116
     117    axes = psEllipseMomentsToAxes(tmpMoments);
     118    shape = psEllipseAxesToShape(axes);
     119
     120    PAR[PM_PAR_SKY] = moments->Sky;
     121    PAR[PM_PAR_I0] = moments->Peak - moments->Sky;
     122    PAR[PM_PAR_XPOS] = peak->x;
     123    PAR[PM_PAR_YPOS] = peak->y;
     124    PAR[PM_PAR_SXX]  = PS_MAX(0.5, moments->Sx);
     125    PAR[PM_PAR_SYY]  = PS_MAX(0.5, moments->Sy);
     126    PAR[PM_PAR_SXY]  = shape.sxy;
     127    PAR[PM_PAR_7]    = 2.0;
     128    return(true);
     129}
     130
     131psF64 PS_MODEL_FLUX (const psVector *params)
     132{
     133    float norm, z;
     134    psEllipseShape shape;
    62135
    63136    psF32 *PAR = params->data.F32;
    64137
    65     psF64 A1   = PS_SQR(PAR[4]);
    66     psF64 A2   = PS_SQR(PAR[5]);
    67     psF64 A3   = PS_SQR(PAR[6]);
    68     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     138    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     139    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     140    shape.sxy = PAR[PM_PAR_SXY];
     141
    69142    // Area is equivalent to 2 pi sigma^2
     143    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     144    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    70145
    71146    // the area needs to be multiplied by the integral of f(z)
    72147    norm = 0.0;
    73     for (z = 0.005; z < 50; z += 0.01) {
    74         f = 1.0 / (1 + z + pow(z, PAR[7]));
    75         norm += f;
    76     }
    77     norm *= 0.01;
    78 
    79     psF64 Flux = params->data.F32[1] * Area * norm;
     148
     149    # define DZ 0.25
     150
     151    float f0 = 1.0;
     152    float f1, f2;
     153    for (z = DZ; z < 50; z += DZ) {
     154        f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
     155        z += DZ;
     156        f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
     157        norm += f0 + 4*f1 + f2;
     158        f0 = f2;
     159    }
     160    norm *= DZ / 3.0;
     161
     162    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
    80163
    81164    return(Flux);
    … …  
    84167// define this function so it never returns Inf or NaN
    85168// return the radius which yields the requested flux
    86 psF64 psModelRadius_RGAUSS (const psVector *params, psF64 flux)
     169psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux)
    87170{
    88171    psF64 z, f, p;
     172    psEllipseShape shape;
     173
    89174    psF32 *PAR = params->data.F32;
    90175
    91176    if (flux <= 0)
    92177        return (1.0);
    93     if (PAR[1] <= 0)
     178    if (PAR[PM_PAR_I0] <= 0)
    94179        return (1.0);
    95     if (flux >= PAR[1])
     180    if (flux >= PAR[PM_PAR_I0])
    96181        return (1.0);
    97182
    98     // if Sx == Sy, sigma = Sx == Sy
    99     psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
     183    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     184    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     185    shape.sxy = PAR[PM_PAR_SXY];
     186
     187    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     188    psF64 sigma = axes.major;
     189
    100190    psF64 dz = 1.0 / (2.0 * sigma*sigma);
    101     psF64 limit = flux / PAR[1];
    102 
    103     // we can do this much better with intelligent choices here
    104     for (z = 0.0; z < 20.0; z += dz) {
    105         p = pow(z, PAR[7]);
    106         f = 1.0 / (1 + z + p);
    107         if (f < limit)
    108             break;
     191    psF64 limit = flux / PAR[PM_PAR_I0];
     192
     193    // use the fact that f is monotonically decreasing
     194    z = 0;
     195    Nstep = 0;
     196
     197    // choose a z value guaranteed to be beyond our limit
     198    float z0 = pow((1.0 / limit), (1.0 / 2.25));
     199    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
     200    z1 = PS_MAX (z0, z1);
     201    z0 = 0.0;
     202
     203    // perform a type of bisection to find the value
     204    float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7]));
     205    float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7]));
     206    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
     207        z = 0.5*(z0 + z1);
     208        f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
     209        if (f > limit) {
     210            z0 = z;
     211            f0 = f;
     212        } else {
     213            z1 = z;
     214            f1 = f;
     215        }
     216        Nstep ++;
    109217    }
    110218    psF64 radius = sigma * sqrt (2.0 * z);
    111     if (isnan(radius)) {
    112         fprintf (stderr, "error in code\n");
    113     }
     219
     220    if (isnan(radius))
     221        psAbort ("psphot.model", "error in code: radius is NaN");
     222
    114223    return (radius);
    115224}
    116225
    117 bool psModelGuess_RGAUSS (psModel *model, psSource *source)
    118 {
    119 
    120     psVector *params = model->params;
    121 
    122     psEllipseAxes axes;
    123     psEllipseShape shape;
    124     psEllipseMoments moments;
    125 
    126     moments.x2 = PS_SQR(source->moments->Sx);
    127     moments.y2 = PS_SQR(source->moments->Sy);
    128     moments.xy = source->moments->Sxy;
    129 
    130     axes = psEllipseMomentsToAxes(moments);
    131     shape = psEllipseAxesToShape(axes);
    132 
    133     params->data.F32[0] = source->moments->Sky;
    134     params->data.F32[1] = source->peak->counts - source->moments->Sky;
    135     params->data.F32[2] = source->moments->x;
    136     params->data.F32[3] = source->moments->y;
    137     params->data.F32[4] = 1.0 / shape.sx;
    138     params->data.F32[5] = 1.0 / shape.sy;
    139     params->data.F32[6] = shape.sxy;
    140     params->data.F32[7] = 2.0;
    141     return(true);
    142 }
    143 
    144 bool psModelFromPSF_RGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
     226bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
    145227{
    146228
    … …  
    148230    psF32 *in  = modelFLT->params->data.F32;
    149231
    150     out[0] = in[0];
    151     out[1] = in[1];
    152     out[2] = in[2];
    153     out[3] = in[3];
    154 
    155     for (int i = 4; i < 8; i++) {
    156         psPolynomial2D *poly = psf->params->data[i-4];
    157         out[i] = Polynomial2DEval (poly, out[2], out[3]);
    158     }
     232    // we require these two parameters to exist
     233    assert (psf->params_NEW->n > PM_PAR_YPOS);
     234    assert (psf->params_NEW->n > PM_PAR_XPOS);
     235
     236    for (int i = 0; i < psf->params_NEW->n; i++) {
     237        if (psf->params_NEW->data[i] == NULL) {
     238            out[i] = in[i];
     239        } else {
     240            psPolynomial2D *poly = psf->params_NEW->data[i];
     241            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
     242        }
     243    }
     244
     245    // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here
     246    out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
     247
    159248    return(true);
    160249}
     250
     251bool PM_MODEL_FIT_STATUS (pmModel *model)
     252{
     253
     254    psF32 dP;
     255    bool  status;
     256
     257    psF32 *PAR  = model->params->data.F32;
     258    psF32 *dPAR = model->dparams->data.F32;
     259
     260    dP = 0;
     261    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     262    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
     263    dP = sqrt (dP);
     264
     265    status = true;
     266    status &= (dP < 0.5);
     267    status &= (PAR[PM_PAR_I0] > 0);
     268    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     269
     270    return status;
     271}
     272
     273# undef PM_MODEL_FUNC
     274# undef PM_MODEL_FLUX
     275# undef PM_MODEL_GUESS
     276# undef PM_MODEL_LIMITS
     277# undef PM_MODEL_RADIUS
     278# undef PM_MODEL_FROM_PSF
     279# undef PM_MODEL_FIT_STATUS
Note: See TracChangeset for help on using the changeset viewer.