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_QGAUSS.c

    r9730 r9775  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    51/******************************************************************************
    6     one component, two slopes:
    7     1 / (1 + z^M + z^N)
     2 * this file defines the QGAUSS 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[PM_PAR_SKY] = So;
    10     params->data.F32[PM_PAR_I0] = Zo;
    11     params->data.F32[PM_PAR_XPOS] = Xo;
    12     params->data.F32[PM_PAR_YPOS] = Yo;
    13     params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX;
    14     params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY;
    15     params->data.F32[PM_PAR_SXY] = Sxy;
    16     params->data.F32[PM_PAR_7] =
    17     params->data.F32[PM_PAR_8] =
    18 *****************************************************************************/
    19 
    20 psF32 pmModelFunc_QGAUSS(psVector *deriv,
    21                          const psVector *params,
    22                          const psVector *x)
     9   power-law with fitted linear term
     10   1 / (1 + kz + z^2.25)
     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   - amplitude of the linear component (k)
     20   *****************************************************************************/
     21
     22# define PM_MODEL_FUNC       pmModelFunc_QGAUSS
     23# define PM_MODEL_FLUX       pmModelFlux_QGAUSS
     24# define PM_MODEL_GUESS      pmModelGuess_QGAUSS
     25# define PM_MODEL_LIMITS     pmModelLimits_QGAUSS
     26# define PM_MODEL_RADIUS     pmModelRadius_QGAUSS
     27# define PM_MODEL_FROM_PSF   pmModelFromPSF_QGAUSS
     28# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS
     29
     30psF32 PM_MODEL_FUNC (psVector *deriv,
     31                     const psVector *params,
     32                     const psVector *pixcoord)
    2333{
    2434    psF32 *PAR = params->data.F32;
    2535
    26     psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
    27     psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
    28     psF32 px = PAR[PM_PAR_SXX]*X;
    29     psF32 py = PAR[PM_PAR_SYY]*Y;
    30     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
     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) + PAR[PM_PAR_SXY]*X*Y;
     41
     42    // XXX if the elliptical contour is defined in valid way, this step should not be required.
     43    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
     44    // for negative values of z
    3145    if (z < 0)
    3246        z = 0;
    … …  
    4559        psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
    4660
    47         dPAR[PM_PAR_SKY] = +1.0;
    48         dPAR[PM_PAR_I0] = +r;
    49         dPAR[PM_PAR_XPOS] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
    50         dPAR[PM_PAR_YPOS] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
    51         dPAR[PM_PAR_SXX] = -2.0*q*px*X;
    52         dPAR[PM_PAR_SYY] = -2.0*q*py*Y;
    53         dPAR[PM_PAR_SXY] = -q*X*Y;
    54         dPAR[PM_PAR_7] = -t*z;
     61        dPAR[PM_PAR_SKY]  = +1.0;
     62        dPAR[PM_PAR_I0]   = +r;
     63        dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
     64        dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     65        dPAR[PM_PAR_SXX]  = +2.0*q*px*px/PAR[PM_PAR_SXX];
     66        dPAR[PM_PAR_SYY]  = +2.0*q*py*py/PAR[PM_PAR_SYY];
     67        dPAR[PM_PAR_SXY]  = -q*X*Y;
     68        dPAR[PM_PAR_7]    = -t*z;
    5569    }
    5670    return(f);
    5771}
    5872
    59 bool pmModelLimits_QGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
     73bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max)
    6074{
    6175
    … …  
    7791    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
    7892    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
    79     params_min[0][0].data.F32[PM_PAR_SXX] = 0.01;
    80     params_min[0][0].data.F32[PM_PAR_SYY] = 0.01;
     93    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     94    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
    8195    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    8296    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
    … …  
    86100    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
    87101    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
    88     params_max[0][0].data.F32[PM_PAR_SXX] = 2.0;
    89     params_max[0][0].data.F32[PM_PAR_SYY] = 2.0;
     102    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     103    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
    90104    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    91105    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
    … …  
    95109
    96110// make an initial guess for parameters
    97 bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source)
     111bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    98112{
    99113
    100114    pmMoments *moments = source->moments;
    101115    pmPeak    *peak    = source->peak;
    102     psF32     *params  = model->params->data.F32;
    103 
    104     params[PM_PAR_SKY] = moments->Sky;
    105     params[PM_PAR_I0] = moments->Peak - moments->Sky;
    106     params[PM_PAR_XPOS] = peak->x;
    107     params[PM_PAR_YPOS] = peak->y;
    108     params[PM_PAR_SXX] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
    109     params[PM_PAR_SYY] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
    110     params[PM_PAR_SXY] = 0.0;
    111     params[PM_PAR_7] = 1.0;
     116    psF32     *PAR  = model->params->data.F32;
     117
     118    PAR[PM_PAR_SKY] = moments->Sky;
     119    PAR[PM_PAR_I0]  = moments->Peak - moments->Sky;
     120    PAR[PM_PAR_XPOS] = peak->x;
     121    PAR[PM_PAR_YPOS] = peak->y;
     122    PAR[PM_PAR_SXX]  = PS_MAX(0.5, moments->Sx);
     123    PAR[PM_PAR_SYY]  = PS_MAX(0.5, moments->Sy);
     124    PAR[PM_PAR_SXY] = 0.0;
     125    PAR[PM_PAR_7]    = 1.0;
    112126
    113127    return(true);
    114128}
    115129
    116 psF64 pmModelFlux_QGAUSS(const psVector *params)
    117 {
    118     float z;
    119     float norm;
     130psF64 PM_MODEL_FLUX (const psVector *params)
     131{
     132    float z, norm;
     133    psEllipseShape shape;
    120134
    121135    psF32 *PAR = params->data.F32;
    122136
    123     psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
    124     psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
    125     psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
    126     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     137    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     138    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     139    shape.sxy = PAR[PM_PAR_SXY];
     140
    127141    // Area is equivalent to 2 pi sigma^2
     142    psEllipseAxes axes = psEllipseShapeToAxes (shape);
     143    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    128144
    129145    // the area needs to be multiplied by the integral of f(z)
    … …  
    150166// define this function so it never returns Inf or NaN
    151167// return the radius which yields the requested flux
    152 psF64 pmModelRadius_QGAUSS (const psVector *params, psF64 flux)
     168psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    153169{
    154170    psF64 z, f;
     171    int Nstep = 0;
     172    psEllipseShape shape;
     173
    155174    psF32 *PAR = params->data.F32;
    156     int Nstep = 0;
    157175
    158176    if (flux <= 0)
    … …  
    163181        return (1.0);
    164182
    165     // if Sx == Sy, sigma = Sx == Sy
    166     psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / 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
    167190    psF64 limit = flux / PAR[PM_PAR_I0];
    168191
    169     # if 0
    170     /* test example will just use both, printing both */
    171     psF64 dz = 1.0 / (2.0 * sigma*sigma);
    172 
    173     // we can do this much better with intelligent choices here
    174     for (z = 0.0; z < 30.0; z += dz) {
    175         Nstep ++;
    176         f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    177         // test: f = 1.0 / (1 + PAR[PM_PAR_7]*z + PS_SQR(z));
    178         if (f < limit)
    179             break;
    180     }
    181     // fprintf (stderr, "rad 1: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep);
    182 
    183     # else
    184 
    185         /* use the fact that f is monotonically decreasing */
    186         z = 0;
     192    // use the fact that f is monotonically decreasing
     193    z = 0;
    187194    Nstep = 0;
    188195
    … …  
    193200    z0 = 0.0;
    194201
     202    // perform a type of bisection to find the value
    195203    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
    196204    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
    … …  
    198206        z = 0.5*(z0 + z1);
    199207        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    200         // fprintf (stderr, "%f  %f  %f   :   %f  %f  %f\n", f0, f, f1, z0, z, z1);
    201208        if (f > limit) {
    202209            z0 = z;
    … …  
    208215        Nstep ++;
    209216    }
    210     // fprintf (stderr, "rad 2: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep);
    211     # endif
    212 
    213217    psF64 radius = sigma * sqrt (2.0 * z);
    214     if (isnan(radius)) {
    215         fprintf (stderr, "error in code\n");
    216     }
     218
     219    if (isnan(radius))
     220        psAbort ("psphot.model", "error in code: radius is NaN");
     221
    217222    return (radius);
    218223}
    219224
    220 bool pmModelFromPSF_QGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
     225bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
    221226{
    222227
    … …  
    224229    psF32 *in  = modelFLT->params->data.F32;
    225230
    226     out[PM_PAR_SKY] = in[PM_PAR_SKY];
    227     out[PM_PAR_I0] = in[PM_PAR_I0];
    228     out[PM_PAR_XPOS] = in[PM_PAR_XPOS];
    229     out[PM_PAR_YPOS] = in[PM_PAR_YPOS];
    230 
    231     assert (PM_PAR_YPOS + 1 == 4);  // so starting at 4 is correct
    232     for (int i = 4; i < 8; i++) {
    233         psPolynomial2D *poly = psf->params->data[i-4];
    234         // XXX: Verify this (from EAM change)
    235         //out[i] = Polynomial2DEval_EAM(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]);
    236         out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]);
    237     }
     231    // we require these two parameters to exist
     232    assert (psf->params_NEW->n > PM_PAR_YPOS);
     233    assert (psf->params_NEW->n > PM_PAR_XPOS);
     234
     235    for (int i = 0; i < psf->params_NEW->n; i++) {
     236        if (psf->params_NEW->data[i] == NULL) {
     237            out[i] = in[i];
     238        } else {
     239            psPolynomial2D *poly = psf->params_NEW->data[i];
     240            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
     241        }
     242    }
     243
     244    // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here
     245    out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
     246
    238247    return(true);
    239248}
    240249
    241 bool pmModelFitStatus_QGAUSS (pmModel *model)
     250bool PM_MODEL_FIT_STATUS (pmModel *model)
    242251{
    243252
    … …  
    258267    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    259268
    260     if (!status)
    261         return false;
    262     return true;
    263 }
     269    return status;
     270}
     271
     272# undef PM_MODEL_FUNC
     273# undef PM_MODEL_FLUX
     274# undef PM_MODEL_GUESS
     275# undef PM_MODEL_LIMITS
     276# undef PM_MODEL_RADIUS
     277# undef PM_MODEL_FROM_PSF
     278# undef PM_MODEL_FIT_STATUS
Note: See TracChangeset for help on using the changeset viewer.