IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9775


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

Location:
trunk/psModules/src/objects/models
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r9770 r9775  
    11/******************************************************************************
    2  * this file defines the PGAUSS source shape model.  Note that these model functions are loaded
     2 * this file defines the GAUSS source shape model.  Note that these model functions are loaded
    33 * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
    44 * models use a psVector to represent the set of parameters, with the sequence used to specify
    … …  
    66 * specifics of the model.  All models which are used a PSF representations share a few
    77 * parameters, for which # define names are listed in pmModel.h:
     8 
     9   pure Gaussian:
     10   exp(-z)
    811 
    912 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
    … …  
    114117    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
    115118    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
    116     shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     119    shape.sxy = PAR[PM_PAR_SXY];
    117120
    118121    // Area is equivalent to 2 pi sigma^2
    … …  
    131134    psEllipseShape shape;
    132135
     136    psF32 *PAR = params->data.F32;
     137
    133138    if (flux <= 0)
    134139        return (1.0);
    135     if (params->data.F32[PM_PAR_I0] <= 0)
     140    if (PAR[PM_PAR_I0] <= 0)
    136141        return (1.0);
    137     if (flux >= params->data.F32[PM_PAR_I0])
     142    if (flux >= PAR[PM_PAR_I0])
    138143        return (1.0);
    139 
    140     psF32 *PAR = params->data.F32;
    141144
    142145    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
    143146    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
    144     shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     147    shape.sxy = PAR[PM_PAR_SXY];
    145148
    146149    psEllipseAxes axes = psEllipseShapeToAxes (shape);
    147     psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
     150    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
    148151    return (radius);
    149152}

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

    r9770 r9775  
    66 * specifics of the model.  All models which are used a PSF representations share a few
    77 * parameters, for which # define names are listed in pmModel.h:
     8 
     9   Gaussian taylor expansion
     10   1 / (1 + z + z^2/2 + z^3/6)
    811 
    912 * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
    … …  
    115118    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
    116119    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
    117     shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     120    shape.sxy = PAR[PM_PAR_SXY];
    118121
    119122    // Area is equivalent to 2 pi sigma^2
    … …  
    146149psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    147150{
     151    psEllipseShape shape;
     152
     153    psF32 *PAR = params->data.F32;
     154
    148155    if (flux <= 0)
    149156        return (1.0);
    150     if (params->data.F32[PM_PAR_I0] <= 0)
     157    if (PAR[PM_PAR_I0] <= 0)
    151158        return (1.0);
    152     if (flux >= params->data.F32[PM_PAR_I0])
     159    if (flux >= PAR[PM_PAR_I0])
    153160        return (1.0);
    154 
    155     psF32 *PAR = params->data.F32;
    156161
    157162    shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
    158163    shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
    159     shape.sxy = PAR[PM_PAR_SXY] / sqrt(2.0);
     164    shape.sxy = PAR[PM_PAR_SXY];
    160165
    161166    // this estimates the radius assuming f(z) is roughly exp(-z)
    162167    psEllipseAxes axes = psEllipseShapeToAxes (shape);
    163     psF64 radius = axes.major * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
     168    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
    164169
    165170    if (isnan(radius))
    166         psAbort ("psphot.model", "error in code: never return invalid radius");
     171        psAbort ("psphot.model", "error in code: radius is NaN");
    167172    if (radius < 0)
    168         psAbort ("psphot.model", "error in code: never return invalid radius");
     173        psAbort ("psphot.model", "error in code: radius is negative");
    169174
    170175    return (radius);
    … …  
    215220    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    216221
    217     if (status)
    218         return true;
    219     return false;
     222    return status;
    220223}
    221224

  • 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

  • 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

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

    r9730 r9775  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 
    61/******************************************************************************
    7     one component, two slopes:
    8     1 / (1 + z^Npow + St*z^Ntide)
     2 * this file defines the SGAUSS 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:
    98 
    10     params->data.F32[0] = So;
    11     params->data.F32[1] = Zo;
    12     params->data.F32[2] = Xo;
    13     params->data.F32[3] = Yo;
    14     params->data.F32[4] = 1 / SigmaX;
    15     params->data.F32[5] = 1 / SigmaY;
    16     params->data.F32[6] = Sxy;
    17     params->data.F32[7] = Npow
    18     params->data.F32[8] = St
    19 *****************************************************************************/
    20 
    21 # define SQ(A)((A)*(A))
     9   power-law with fitted slope and outer tidal radius
     10   1 / (1 + z^N + kz^4)
     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   - slope of power-law component (N)
     20   * PM_PAR_8   8   - amplitude of the tidal cutoff (k)
     21   *****************************************************************************/
     22
     23/***
     24    XXXX the model in this file needs to be tested more carefully.
     25    the code for guessing the power-law slope based on the radial profile
     26    is either too slow or does not work well.
     27    fix up the code to follow conventions in the other model function files.
     28***/
     29
     30XXX broken code
     31
     32# define PM_MODEL_FUNC       pmModelFunc_SGAUSS
     33# define PM_MODEL_FLUX       pmModelFlux_SGAUSS
     34# define PM_MODEL_GUESS      pmModelGuess_SGAUSS
     35# define PM_MODEL_LIMITS     pmModelLimits_SGAUSS
     36# define PM_MODEL_RADIUS     pmModelRadius_SGAUSS
     37# define PM_MODEL_FROM_PSF   pmModelFromPSF_SGAUSS
     38# define PM_MODEL_FIT_STATUS pmModelFitStatus_SGAUSS
     39
    2240psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image);
    23 psF64 p_psImageGetElementF64(psImage *a, int i, int j);
    24 
    25 psF32 pmModelFunc_SGAUSS(psVector *deriv,
    26                          const psVector *params,
    27                          const psVector *x)
     41
     42psF32 PM_MODEL_FUNC (psVector *deriv,
     43                     const psVector *params,
     44                     const psVector *x)
    2845{
    2946    psF32 *PAR = params->data.F32;
    … …  
    6178}
    6279
    63 bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
     80bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max)
    6481{
    6582
    … …  
    99116
    100117    return (TRUE);
     118}
     119
     120bool PM_MODEL_GUESS  (pmModel *model, pmSource *source)
     121{
     122
     123    pmMoments *sMoments = source->moments;
     124    pmPeak    *peak     = source->peak;
     125    psF32     *params   = model->params->data.F32;
     126
     127    psEllipseAxes axes;
     128    psEllipseShape shape;
     129    psEllipseMoments moments;
     130
     131    moments.x2 = PS_SQR(sMoments->Sx);
     132    moments.y2 = PS_SQR(sMoments->Sy);
     133    moments.xy = sMoments->Sxy;
     134
     135    // solve the math to go from Moments To Shape
     136    axes = psEllipseMomentsToAxes(moments);
     137    shape = psEllipseAxesToShape(axes);
     138
     139    params[0] = sMoments->Sky;
     140    params[1] = sMoments->Peak - sMoments->Sky;
     141    params[2] = peak->x;
     142    params[3] = peak->y;
     143    params[4] = 1.0 / shape.sx;
     144    params[5] = 1.0 / shape.sy;
     145    params[6] = shape.sxy;
     146
     147    # if (0)
     148
     149        f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
     150    axes.major *= 2.0;
     151    axes.minor *= 2.0;
     152    f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
     153
     154    if (f1 > f2) {
     155        params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
     156    } else {
     157        params[7] = 0.5;
     158    }
     159    # endif
     160
     161    params[7] = 1.8;
     162    params[8] = 0.1;
     163
     164
     165    return(true);
     166}
     167
     168psF64 PM_MODEL_FLUX (const psVector *params)
     169{
     170    float f, norm, z;
     171
     172    psF32 *PAR = params->data.F32;
     173
     174    psF64 A1   = PS_SQR(PAR[4]);
     175    psF64 A2   = PS_SQR(PAR[5]);
     176    psF64 A3   = PS_SQR(PAR[6]);
     177    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     178    // Area is equivalent to 2 pi sigma^2
     179
     180    // the area needs to be multiplied by the integral of f(z)
     181    norm = 0.0;
     182    for (z = 0.005; z < 50; z += 0.01) {
     183        psF32 pr = PAR[8]*z;
     184        f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr)));
     185        norm += f;
     186    }
     187    norm *= 0.01;
     188
     189    psF64 Flux = PAR[1] * Area * norm;
     190
     191    return(Flux);
     192}
     193
     194// XXX need to define the radius along the major axis
     195// define this function so it never returns Inf or NaN
     196// return the radius which yields the requested flux
     197psF64 PM_MODEL_RADIUS   (const psVector *params, psF64 flux)
     198{
     199    psF64 r, z = 0.0, pr, f;
     200    psF32 *PAR = params->data.F32;
     201
     202    psEllipseAxes axes;
     203    psEllipseShape shape;
     204
     205    if (flux <= 0)
     206        return (1.0);
     207    if (PAR[1] <= 0)
     208        return (1.0);
     209    if (flux >= PAR[1])
     210        return (1.0);
     211
     212    // convert Sx,Sy,Sxy to major/minor axes
     213    shape.sx = 1.0 / PAR[4];
     214    shape.sy = 1.0 / PAR[5];
     215    shape.sxy = PAR[6];
     216
     217    axes = psEllipseShapeToAxes (shape);
     218    psF64 dr = 1.0 / axes.major;
     219    psF64 limit = flux / PAR[1];
     220
     221    // XXX : we can do this faster with an intelligent starting choice
     222    for (r = 0.0; r < 20.0; r += dr) {
     223        z = PS_SQR(r);
     224        pr = PAR[8]*z;
     225        f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr)));
     226        if (f < limit)
     227            break;
     228    }
     229    psF64 radius = 2.0 * axes.major * sqrt (z);
     230    if (isnan(radius)) {
     231        fprintf (stderr, "error in code\n");
     232    }
     233    return (radius);
     234}
     235
     236bool PM_MODEL_FROM_PSF  (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
     237{
     238
     239    psF32 *out = modelPSF->params->data.F32;
     240    psF32 *in  = modelFLT->params->data.F32;
     241
     242    out[0] = in[0];
     243    out[1] = in[1];
     244    out[2] = in[2];
     245    out[3] = in[3];
     246
     247    for (int i = 4; i < 9; i++) {
     248        psPolynomial2D *poly = psf->params->data[i-4];
     249        // XXX: Verify this (from EAM change)
     250        //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]);
     251        out[i] = psPolynomial2DEval(poly, out[2], out[3]);
     252    }
     253    return(true);
     254}
     255
     256bool PM_MODEL_FIT_STATUS  (pmModel *model)
     257{
     258
     259    psF32 dP;
     260    bool  status;
     261    psEllipseAxes axes;
     262    psEllipseShape shape;
     263
     264    psF32 *PAR  = model->params->data.F32;
     265    psF32 *dPAR = model->dparams->data.F32;
     266
     267    shape.sx = 1.0 / PAR[4];
     268    shape.sy = 1.0 / PAR[5];
     269    shape.sxy = PAR[6];
     270
     271    axes = psEllipseShapeToAxes (shape);
     272
     273    dP = 0;
     274    dP += PS_SQR(dPAR[4] / PAR[4]);
     275    dP += PS_SQR(dPAR[5] / PAR[5]);
     276    dP += PS_SQR(dPAR[7] / PAR[7]);
     277    dP = sqrt (dP);
     278
     279    status = true;
     280    status &= (dP < 0.5);
     281    status &= (PAR[1] > 0);
     282    status &= ((dPAR[1]/PAR[1]) < 0.5);
     283    status &= (fabs(PAR[8]) < 0.5);
     284    status &= (dPAR[8] < 0.1);
     285    status &= (axes.major > 1.41);
     286    status &= (axes.minor > 1.41);
     287    status &= ((axes.major / axes.minor) < 5.0);
     288    status &= (PAR[7] > 0.5);
     289
     290    if (status)
     291        return true;
     292    return false;
    101293}
    102294
    … …  
    142334}
    143335
    144 bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source)
     336// XXX EAM : test version using flux contours to guess slope
     337bool PM_MODEL_GUESS_HARD (pmModel *model, pmSource *source)
    145338{
    146339
    … …  
    148341    pmPeak    *peak     = source->peak;
    149342    psF32     *params   = model->params->data.F32;
     343    float f1, f2;
    150344
    151345    psEllipseAxes axes;
    … …  
    169363    params[6] = shape.sxy;
    170364
    171     # if (0)
    172 
    173         f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
     365    f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
    174366    axes.major *= 2.0;
    175367    axes.minor *= 2.0;
    … …  
    181373        params[7] = 0.5;
    182374    }
    183     # endif
    184 
    185     params[7] = 1.8;
    186375    params[8] = 0.1;
    187376
    188 
    189377    return(true);
    190378}
    191379
    192 // XXX EAM : test version using flux contours to guess slope
    193 bool pmModelGuess_SGAUSS_HARD (pmModel *model, pmSource *source)
    194 {
    195 
    196     pmMoments *sMoments = source->moments;
    197     pmPeak    *peak     = source->peak;
    198     psF32     *params   = model->params->data.F32;
    199     float f1, f2;
    200 
    201     psEllipseAxes axes;
    202     psEllipseShape shape;
    203     psEllipseMoments moments;
    204 
    205     moments.x2 = PS_SQR(sMoments->Sx);
    206     moments.y2 = PS_SQR(sMoments->Sy);
    207     moments.xy = sMoments->Sxy;
    208 
    209     // solve the math to go from Moments To Shape
    210     axes = psEllipseMomentsToAxes(moments);
    211     shape = psEllipseAxesToShape(axes);
    212 
    213     params[0] = sMoments->Sky;
    214     params[1] = sMoments->Peak - sMoments->Sky;
    215     params[2] = peak->x;
    216     params[3] = peak->y;
    217     params[4] = 1.0 / shape.sx;
    218     params[5] = 1.0 / shape.sy;
    219     params[6] = shape.sxy;
    220 
    221     f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
    222     axes.major *= 2.0;
    223     axes.minor *= 2.0;
    224     f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
    225 
    226     if (f1 > f2) {
    227         params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
    228     } else {
    229         params[7] = 0.5;
    230     }
    231     params[8] = 0.1;
    232 
    233     return(true);
    234 }
    235 
    236 psF64 pmModelFlux_SGAUSS(const psVector *params)
    237 {
    238     float f, norm, z;
    239 
    240     psF32 *PAR = params->data.F32;
    241 
    242     psF64 A1   = PS_SQR(PAR[4]);
    243     psF64 A2   = PS_SQR(PAR[5]);
    244     psF64 A3   = PS_SQR(PAR[6]);
    245     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    246     // Area is equivalent to 2 pi sigma^2
    247 
    248     // the area needs to be multiplied by the integral of f(z)
    249     norm = 0.0;
    250     for (z = 0.005; z < 50; z += 0.01) {
    251         psF32 pr = PAR[8]*z;
    252         f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
    253         norm += f;
    254     }
    255     norm *= 0.01;
    256 
    257     psF64 Flux = PAR[1] * Area * norm;
    258 
    259     return(Flux);
    260 }
    261 
    262 // XXX need to define the radius along the major axis
    263 // define this function so it never returns Inf or NaN
    264 // return the radius which yields the requested flux
    265 psF64 pmModelRadius_SGAUSS  (const psVector *params, psF64 flux)
    266 {
    267     psF64 r, z = 0.0, pr, f;
    268     psF32 *PAR = params->data.F32;
    269 
    270     psEllipseAxes axes;
    271     psEllipseShape shape;
    272 
    273     if (flux <= 0)
    274         return (1.0);
    275     if (PAR[1] <= 0)
    276         return (1.0);
    277     if (flux >= PAR[1])
    278         return (1.0);
    279 
    280     // convert Sx,Sy,Sxy to major/minor axes
    281     shape.sx = 1.0 / PAR[4];
    282     shape.sy = 1.0 / PAR[5];
    283     shape.sxy = PAR[6];
    284 
    285     axes = psEllipseShapeToAxes (shape);
    286     psF64 dr = 1.0 / axes.major;
    287     psF64 limit = flux / PAR[1];
    288 
    289     // XXX : we can do this faster with an intelligent starting choice
    290     for (r = 0.0; r < 20.0; r += dr) {
    291         z = SQ(r);
    292         pr = PAR[8]*z;
    293         f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
    294         if (f < limit)
    295             break;
    296     }
    297     psF64 radius = 2.0 * axes.major * sqrt (z);
    298     if (isnan(radius)) {
    299         fprintf (stderr, "error in code\n");
    300     }
    301     return (radius);
    302 }
    303 
    304 bool pmModelFromPSF_SGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
    305 {
    306 
    307     psF32 *out = modelPSF->params->data.F32;
    308     psF32 *in  = modelFLT->params->data.F32;
    309 
    310     out[0] = in[0];
    311     out[1] = in[1];
    312     out[2] = in[2];
    313     out[3] = in[3];
    314 
    315     for (int i = 4; i < 9; i++) {
    316         psPolynomial2D *poly = psf->params->data[i-4];
    317         // XXX: Verify this (from EAM change)
    318         //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]);
    319         out[i] = psPolynomial2DEval(poly, out[2], out[3]);
    320     }
    321     return(true);
    322 }
    323 
    324 bool pmModelFitStatus_SGAUSS (pmModel *model)
    325 {
    326 
    327     psF32 dP;
    328     bool  status;
    329     psEllipseAxes axes;
    330     psEllipseShape shape;
    331 
    332     psF32 *PAR  = model->params->data.F32;
    333     psF32 *dPAR = model->dparams->data.F32;
    334 
    335     shape.sx = 1.0 / PAR[4];
    336     shape.sy = 1.0 / PAR[5];
    337     shape.sxy = PAR[6];
    338 
    339     axes = psEllipseShapeToAxes (shape);
    340 
    341     dP = 0;
    342     dP += PS_SQR(dPAR[4] / PAR[4]);
    343     dP += PS_SQR(dPAR[5] / PAR[5]);
    344     dP += PS_SQR(dPAR[7] / PAR[7]);
    345     dP = sqrt (dP);
    346 
    347     status = true;
    348     status &= (dP < 0.5);
    349     status &= (PAR[1] > 0);
    350     status &= ((dPAR[1]/PAR[1]) < 0.5);
    351     status &= (fabs(PAR[8]) < 0.5);
    352     status &= (dPAR[8] < 0.1);
    353     status &= (axes.major > 1.41);
    354     status &= (axes.minor > 1.41);
    355     status &= ((axes.major / axes.minor) < 5.0);
    356     status &= (PAR[7] > 0.5);
    357 
    358     if (status)
    359         return true;
    360     return false;
    361 }
     380# undef PM_MODEL_FUNC
     381# undef PM_MODEL_FLUX
     382# undef PM_MODEL_GUESS
     383# undef PM_MODEL_LIMITS
     384# undef PM_MODEL_RADIUS
     385# undef PM_MODEL_FROM_PSF
     386# undef PM_MODEL_FIT_STATUS

  • trunk/psModules/src/objects/models/pmModel_TGAUSS.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^M + z^N)
     2 * this file defines the TGAUSS 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] = sqrt(2.0) / SigmaX;
    14     params->data.F32[5] = sqrt(2.0) / SigmaY;
    15     params->data.F32[6] = Sxy;
    16     params->data.F32[7] =
    17 *****************************************************************************/
    18 
    19 psF64 psModelFunc_TGAUSS(psVector *deriv,
    20                          const psVector *params,
    21                          const psVector *x)
     9   power-law with fixed slope and fitted amplitude
     10   1 / (1 + z + kz^2.2)
     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 high-order component (k)
     20   *****************************************************************************/
     21
     22/***
     23    XXXX the model in this file needs to be tested more carefully.
     24    fix up the code to follow conventions in the other model function files.
     25***/
     26
     27XXX broken code
     28
     29# define PM_MODEL_FUNC       pmModelFunc_TGAUSS
     30# define PM_MODEL_FLUX       pmModelFlux_TGAUSS
     31# define PM_MODEL_GUESS      pmModelGuess_TGAUSS
     32# define PM_MODEL_LIMITS     pmModelLimits_TGAUSS
     33# define PM_MODEL_RADIUS     pmModelRadius_TGAUSS
     34# define PM_MODEL_FROM_PSF   pmModelFromPSF_TGAUSS
     35# define PM_MODEL_FIT_STATUS pmModelFitStatus_TGAUSS
     36
     37psF64 PS_MODEL_FUNC (psVector *deriv,
     38                     const psVector *params,
     39                     const psVector *x)
    2240{
    2341    psF32 *PAR = params->data.F32;
    … …  
    5068}
    5169
    52 psF64 psModelFlux_TGAUSS(const psVector *params)
    53 {
    54     psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
    55     psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
    56     psF64 A3   = params->data.F32[6];
    57     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3));
    58     // Area is equivalent to 2 pi sigma^2
    59 
    60     psF64 Flux = params->data.F32[1] * Area;
    61 
    62     return(Flux);
    63 }
    64 
    65 // define this function so it never returns Inf or NaN
    66 // return the radius which yields the requested flux
    67 psF64 psModelRadius_TGAUSS  (const psVector *params, psF64 flux)
    68 {
    69     if (flux <= 0)
    70         return (1.0);
    71     if (params->data.F32[1] <= 0)
    72         return (1.0);
    73     if (flux >= params->data.F32[1])
    74         return (1.0);
    75 
    76     psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
    77     psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
    78     if (isnan(radius)) {
    79         fprintf (stderr, "error in code\n");
    80     }
    81     return (radius);
    82 }
    83 
    84 bool psModelGuess_TGAUSS (psModel *model, psSource *source)
     70bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
     71{
     72
     73    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
     74    *params_min = psVectorAlloc (8, PS_TYPE_F32);
     75    *params_max = psVectorAlloc (8, PS_TYPE_F32);
     76
     77    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     78    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     79    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     80    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     81    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     82    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     83    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     84    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
     85
     86    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     87    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     88    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     89    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     90    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     91    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
     92    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     93    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
     94
     95    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     96    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     97    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     98    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     99    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     100    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
     101    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
     102    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
     103
     104    return (TRUE);
     105}
     106
     107bool PS_MODEL_GUESS  (psModel *model, psSource *source)
    85108{
    86109
    … …  
    99122}
    100123
    101 bool psModelFromPSF_TGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
     124psF64 PS_MODEL_FLUX (const psVector *params)
     125{
     126    psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
     127    psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
     128    psF64 A3   = params->data.F32[6];
     129    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3));
     130    // Area is equivalent to 2 pi sigma^2
     131
     132    psF64 Flux = params->data.F32[1] * Area;
     133
     134    return(Flux);
     135}
     136
     137// define this function so it never returns Inf or NaN
     138// return the radius which yields the requested flux
     139psF64 PS_MODEL_RADIUS   (const psVector *params, psF64 flux)
     140{
     141    if (flux <= 0)
     142        return (1.0);
     143    if (params->data.F32[1] <= 0)
     144        return (1.0);
     145    if (flux >= params->data.F32[1])
     146        return (1.0);
     147
     148    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
     149    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
     150    if (isnan(radius)) {
     151        fprintf (stderr, "error in code\n");
     152    }
     153    return (radius);
     154}
     155
     156bool PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
    102157{
    103158
    … …  
    116171    return(true);
    117172}
     173
     174bool PM_MODEL_FIT_STATUS (pmModel *model)
     175{
     176
     177    psF32 dP;
     178    bool  status;
     179
     180    psF32 *PAR  = model->params->data.F32;
     181    psF32 *dPAR = model->dparams->data.F32;
     182
     183    dP = 0;
     184    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     185    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
     186    dP = sqrt (dP);
     187
     188    status = true;
     189    status &= (dP < 0.5);
     190    status &= (PAR[PM_PAR_I0] > 0);
     191    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     192
     193    return status;
     194}
     195
     196# undef PM_MODEL_FUNC
     197# undef PM_MODEL_FLUX
     198# undef PM_MODEL_GUESS
     199# undef PM_MODEL_LIMITS
     200# undef PM_MODEL_RADIUS
     201# undef PM_MODEL_FROM_PSF
     202# undef PM_MODEL_FIT_STATUS

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

    r9621 r9775  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
     1/******************************************************************************
     2 * this file defines the WAUSS 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:
     8 
     9   power-law with fitted linear term
     10   1 / (1 + Az + Bz^2 + z^3/6)
     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 (SigmaX / sqrt(2))
     17   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
     18   * PM_PAR_SXY 6   - X*Y term of elliptical contour
     19   * PM_PAR_7   7   - amplitude of the linear component (A)
     20   * PM_PAR_8   8   - amplitude of the quadratic component (B)
     21   *****************************************************************************/
    422
    5 /******************************************************************************
    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 *****************************************************************************/
     23/***
     24    XXXX the model in this file needs to be tested more carefully.
     25    fix up the code to follow conventions in the other model function files.
     26***/
    1427
    15 psF64 psModelFunc_WAUSS(psVector *deriv,
    16                         const psVector *params,
    17                         const psVector *x)
     28XXX broken code
     29
     30# define PM_MODEL_FUNC       pmModelFunc_WAUSS
     31# define PM_MODEL_FLUX       pmModelFlux_WAUSS
     32# define PM_MODEL_GUESS      pmModelGuess_WAUSS
     33# define PM_MODEL_LIMITS     pmModelLimits_WAUSS
     34# define PM_MODEL_RADIUS     pmModelRadius_WAUSS
     35# define PM_MODEL_FROM_PSF   pmModelFromPSF_WAUSS
     36# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS
     37
     38psF64 PS_MODEL_FUNC (psVector *deriv,
     39                     const psVector *params,
     40                     const psVector *x)
    1841{
    1942    psF32 X = x->data.F32[0] - params->data.F32[2];
    … …  
    4366}
    4467
    45 // this is probably wrong since it uses the gauss integral 2 pi sigma^2
    46 psF64 psModelFlux_WAUSS(const psVector *params)
     68bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
    4769{
    48     psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
    49     psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
    50     psF64 A3   = params->data.F32[6];
    51     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3));
    52     // Area is equivalent to 2 pi sigma^2
    5370
    54     psF64 Flux = params->data.F32[1] * Area;
     71    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
     72    *params_min = psVectorAlloc (8, PS_TYPE_F32);
     73    *params_max = psVectorAlloc (8, PS_TYPE_F32);
    5574
    56     return(Flux);
     75    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     76    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     77    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     78    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     79    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     80    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     81    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     82    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
     83
     84    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     85    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     86    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     87    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     88    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     89    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
     90    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     91    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
     92
     93    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     94    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     95    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     96    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     97    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     98    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
     99    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
     100    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
     101
     102    return (TRUE);
    57103}
    58104
    59 // return the radius which yields the requested flux
    60 psF64 psModelRadius_WAUSS  (const psVector *params, psF64 flux)
    61 {
    62     if (flux <= 0)
    63         return (1.0);
    64     if (params->data.F32[1] <= 0)
    65         return (1.0);
    66     if (flux >= params->data.F32[1] <= 0)
    67         return (1.0);
    68 
    69     psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
    70     psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
    71     return (radius);
    72 }
    73 
    74 bool psModelGuess_WAUSS (psModel *model, psSource *source)
     105bool PS_MODEL_GUESS  (psModel *model, psSource *source)
    75106{
    76107
    … …  
    90121}
    91122
    92 bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
     123// this is probably wrong since it uses the gauss integral 2 pi sigma^2
     124psF64 PS_MODEL_FLUX (const psVector *params)
     125{
     126    psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
     127    psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
     128    psF64 A3   = params->data.F32[6];
     129    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3));
     130    // Area is equivalent to 2 pi sigma^2
     131
     132    psF64 Flux = params->data.F32[1] * Area;
     133
     134    return(Flux);
     135}
     136
     137// return the radius which yields the requested flux
     138psF64 PS_MODEL_RADIUS   (const psVector *params, psF64 flux)
     139{
     140    if (flux <= 0)
     141        return (1.0);
     142    if (params->data.F32[1] <= 0)
     143        return (1.0);
     144    if (flux >= params->data.F32[1] <= 0)
     145        return (1.0);
     146
     147    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
     148    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
     149    return (radius);
     150}
     151
     152bool PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
    93153{
    94154
    … …  
    107167    return(true);
    108168}
     169
     170bool PM_MODEL_FIT_STATUS (pmModel *model)
     171{
     172
     173    psF32 dP;
     174    bool  status;
     175
     176    psF32 *PAR  = model->params->data.F32;
     177    psF32 *dPAR = model->dparams->data.F32;
     178
     179    dP = 0;
     180    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     181    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
     182    dP = sqrt (dP);
     183
     184    status = true;
     185    status &= (dP < 0.5);
     186    status &= (PAR[PM_PAR_I0] > 0);
     187    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     188
     189    return status;
     190}
     191
     192# undef PM_MODEL_FUNC
     193# undef PM_MODEL_FLUX
     194# undef PM_MODEL_GUESS
     195# undef PM_MODEL_LIMITS
     196# undef PM_MODEL_RADIUS
     197# undef PM_MODEL_FROM_PSF
     198# undef PM_MODEL_FIT_STATUS

  • trunk/psModules/src/objects/models/pmModel_ZGAUSS.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^Npow + PAR8*z^4)
     2 * this file defines the ZGAUSS 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 # define SQ(A)((A)*(A))
     9   power-law with fitted slope and tidal cutoff
     10   1 / (1 + z^N + (Az)^4)
     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   - slope of power-law component (N)
     20   *****************************************************************************/
     21
     22/***
     23    XXXX the model in this file needs to be tested more carefully.
     24    fix up the code to follow conventions in the other model function files.
     25***/
     26
     27XXX broken code
     28
     29# define PM_MODEL_FUNC       pmModelFunc_ZGAUSS
     30# define PM_MODEL_FLUX       pmModelFlux_ZGAUSS
     31# define PM_MODEL_GUESS      pmModelGuess_ZGAUSS
     32# define PM_MODEL_LIMITS     pmModelLimits_ZGAUSS
     33# define PM_MODEL_RADIUS     pmModelRadius_ZGAUSS
     34# define PM_MODEL_FROM_PSF   pmModelFromPSF_ZGAUSS
     35# define PM_MODEL_FIT_STATUS pmModelFitStatus_ZGAUSS
     36
    2037# define PAR8 0.1
    2138
    22 psF64 psModelFunc_ZGAUSS(psVector *deriv,
    23                          const psVector *params,
    24                          const psVector *x)
     39psF64 PS_MODEL_FUNC (psVector *deriv,
     40                     const psVector *params,
     41                     const psVector *x)
    2542{
    2643    psF32 *PAR = params->data.F32;
    … …  
    5471}
    5572
    56 psF64 psModelFlux_ZGAUSS(const psVector *params)
    57 {
    58     float f, norm, z;
    59 
    60     psF32 *PAR = params->data.F32;
    61 
    62     psF64 A1   = PS_SQR(PAR[4]);
    63     psF64 A2   = PS_SQR(PAR[5]);
    64     psF64 A3   = PS_SQR(PAR[6]);
    65     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    66     // Area is equivalent to 2 pi sigma^2
    67 
    68     // the area needs to be multiplied by the integral of f(z)
    69     norm = 0.0;
    70     psF32 pr = PAR8*z;
    71     for (z = 0.005; z < 50; z += 0.01) {
    72         f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
    73         norm += f;
    74     }
    75     norm *= 0.01;
    76 
    77     psF64 Flux = PAR[1] * Area * norm;
    78 
    79     return(Flux);
    80 }
    81 
    82 // XXX need to define the radius along the major axis
    83 // define this function so it never returns Inf or NaN
    84 // return the radius which yields the requested flux
    85 psF64 psModelRadius_ZGAUSS  (const psVector *params, psF64 flux)
    86 {
    87     psF64 r, z, pr, f;
    88     psF32 *PAR = params->data.F32;
    89 
    90     psEllipseAxes axes;
    91     psEllipseShape shape;
    92 
    93     if (flux <= 0)
    94         return (1.0);
    95     if (PAR[1] <= 0)
    96         return (1.0);
    97     if (flux >= PAR[1])
    98         return (1.0);
    99 
    100     // convert Sx,Sy,Sxy to major/minor axes
    101     shape.sx = 1.0 / PAR[4];
    102     shape.sy = 1.0 / PAR[5];
    103     shape.sxy = PAR[6];
    104 
    105     axes = psEllipseShapeToAxes (shape);
    106     psF64 dr = 1.0 / axes.major;
    107     psF64 limit = flux / PAR[1];
    108 
    109     // XXX : we can do this faster with an intelligent starting choice
    110     for (r = 0.0; r < 20.0; r += dr) {
    111         z = SQ(r);
    112         pr = PAR8*z;
    113         f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
    114         if (f < limit)
    115             break;
    116     }
    117     psF64 radius = 2.0 * axes.major * sqrt (z);
    118     if (isnan(radius)) {
    119         fprintf (stderr, "error in code\n");
    120     }
    121     return (radius);
    122 }
    123 
    124 bool psModelGuess_ZGAUSS (psModel *model, psSource *source)
     73bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
     74{
     75
     76    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
     77    *params_min = psVectorAlloc (8, PS_TYPE_F32);
     78    *params_max = psVectorAlloc (8, PS_TYPE_F32);
     79
     80    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     81    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     82    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     83    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     84    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     85    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     86    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     87    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
     88
     89    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     90    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     91    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     92    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     93    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
     94    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
     95    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     96    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
     97
     98    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     99    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     100    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     101    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     102    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
     103    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
     104    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
     105    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
     106
     107    return (TRUE);
     108}
     109
     110bool PS_MODEL_GUESS  (psModel *model, psSource *source)
    125111{
    126112
    … …  
    149135}
    150136
    151 bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
     137psF64 PS_MODEL_FLUX (const psVector *params)
     138{
     139    float f, norm, z;
     140
     141    psF32 *PAR = params->data.F32;
     142
     143    psF64 A1   = PS_SQR(PAR[4]);
     144    psF64 A2   = PS_SQR(PAR[5]);
     145    psF64 A3   = PS_SQR(PAR[6]);
     146    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     147    // Area is equivalent to 2 pi sigma^2
     148
     149    // the area needs to be multiplied by the integral of f(z)
     150    norm = 0.0;
     151    psF32 pr = PAR8*z;
     152    for (z = 0.005; z < 50; z += 0.01) {
     153        f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
     154        norm += f;
     155    }
     156    norm *= 0.01;
     157
     158    psF64 Flux = PAR[1] * Area * norm;
     159
     160    return(Flux);
     161}
     162
     163// XXX need to define the radius along the major axis
     164// define this function so it never returns Inf or NaN
     165// return the radius which yields the requested flux
     166psF64 PS_MODEL_RADIUS   (const psVector *params, psF64 flux)
     167{
     168    psF64 r, z, pr, f;
     169    psF32 *PAR = params->data.F32;
     170
     171    psEllipseAxes axes;
     172    psEllipseShape shape;
     173
     174    if (flux <= 0)
     175        return (1.0);
     176    if (PAR[1] <= 0)
     177        return (1.0);
     178    if (flux >= PAR[1])
     179        return (1.0);
     180
     181    // convert Sx,Sy,Sxy to major/minor axes
     182    shape.sx = 1.0 / PAR[4];
     183    shape.sy = 1.0 / PAR[5];
     184    shape.sxy = PAR[6];
     185
     186    axes = psEllipseShapeToAxes (shape);
     187    psF64 dr = 1.0 / axes.major;
     188    psF64 limit = flux / PAR[1];
     189
     190    // XXX : we can do this faster with an intelligent starting choice
     191    for (r = 0.0; r < 20.0; r += dr) {
     192        z = SQ(r);
     193        pr = PAR8*z;
     194        f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr)));
     195        if (f < limit)
     196            break;
     197    }
     198    psF64 radius = 2.0 * axes.major * sqrt (z);
     199    if (isnan(radius)) {
     200        fprintf (stderr, "error in code\n");
     201    }
     202    return (radius);
     203}
     204
     205bool PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
    152206{
    153207
    … …  
    166220    return(true);
    167221}
     222
     223bool PM_MODEL_FIT_STATUS (pmModel *model)
     224{
     225
     226    psF32 dP;
     227    bool  status;
     228
     229    psF32 *PAR  = model->params->data.F32;
     230    psF32 *dPAR = model->dparams->data.F32;
     231
     232    dP = 0;
     233    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     234    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
     235    dP = sqrt (dP);
     236
     237    status = true;
     238    status &= (dP < 0.5);
     239    status &= (PAR[PM_PAR_I0] > 0);
     240    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     241
     242    return status;
     243}
     244
     245# undef PM_MODEL_FUNC
     246# undef PM_MODEL_FLUX
     247# undef PM_MODEL_GUESS
     248# undef PM_MODEL_LIMITS
     249# undef PM_MODEL_RADIUS
     250# undef PM_MODEL_FROM_PSF
     251# undef PM_MODEL_FIT_STATUS
Note: See TracChangeset for help on using the changeset viewer.