IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 12, 2006, 4:23:55 PM (20 years ago)
Author:
rhl
Message:

Provide symbolic names for model parameters

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6947 r9526  
    44    1 / (1 + z^M + z^N)
    55 
    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     params->data.F32[7] =
    14     params->data.F32[8] =
     6    params->data.F32[PM_PAR_SKY] = So;
     7    params->data.F32[PM_PAR_FLUX] = Zo;
     8    params->data.F32[PM_PAR_XPOS] = Xo;
     9    params->data.F32[PM_PAR_YPOS] = Yo;
     10    params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX;
     11    params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY;
     12    params->data.F32[PM_PAR_SXY] = Sxy;
     13    params->data.F32[PM_PAR_7] =
     14    params->data.F32[PM_PAR_8] =
    1515*****************************************************************************/
    1616
    … …  
    2121    psF32 *PAR = params->data.F32;
    2222
    23     psF32 X  = x->data.F32[0] - PAR[2];
    24     psF32 Y  = x->data.F32[1] - PAR[3];
    25     psF32 px = PAR[4]*X;
    26     psF32 py = PAR[5]*Y;
    27     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;
     23    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
     24    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
     25    psF32 px = PAR[PM_PAR_SXX]*X;
     26    psF32 py = PAR[PM_PAR_SYY]*Y;
     27    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    2828    if (z < 0)
    2929        z = 0;
    3030
    3131    psF32 zp = pow(z,1.25);
    32     psF32 r  = 1.0 / (1 + PAR[7]*z + z*zp);
    33 
    34     psF32 r1 = PAR[1]*r;
    35     psF32 f  = r1 + PAR[0];
     32    psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
     33
     34    psF32 r1 = PAR[PM_PAR_FLUX]*r;
     35    psF32 f  = r1 + PAR[PM_PAR_SKY];
    3636
    3737    if (deriv != NULL) {
    3838        psF32 *dPAR = deriv->data.F32;
    3939
    40         // note difference from a pure gaussian: q = params->data.F32[1]*r
     40        // note difference from a pure gaussian: q = params->data.F32[PM_PAR_FLUX]*r
    4141        psF32 t = r1*r;
    42         psF32 q = t*(PAR[7] + 2.25*zp);
    43 
    44         dPAR[0] = +1.0;
    45         dPAR[1] = +r;
    46         dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    47         dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    48         dPAR[4] = -2.0*q*px*X;
    49         dPAR[5] = -2.0*q*py*Y;
    50         dPAR[6] = -q*X*Y;
    51         dPAR[7] = -t*z;
     42        psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
     43
     44        dPAR[PM_PAR_SKY] = +1.0;
     45        dPAR[PM_PAR_FLUX] = +r;
     46        dPAR[PM_PAR_XPOS] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
     47        dPAR[PM_PAR_YPOS] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
     48        dPAR[PM_PAR_SXX] = -2.0*q*px*X;
     49        dPAR[PM_PAR_SYY] = -2.0*q*py*Y;
     50        dPAR[PM_PAR_SXY] = -q*X*Y;
     51        dPAR[PM_PAR_7] = -t*z;
    5252    }
    5353    return(f);
    … …  
    6464    (*params_max)->n = (*params_max)->nalloc;
    6565
    66     beta_lim[0][0].data.F32[0] = 1000;
    67     beta_lim[0][0].data.F32[1] = 3e6;
    68     beta_lim[0][0].data.F32[2] = 5;
    69     beta_lim[0][0].data.F32[3] = 5;
    70     beta_lim[0][0].data.F32[4] = 0.5;
    71     beta_lim[0][0].data.F32[5] = 0.5;
    72     beta_lim[0][0].data.F32[6] = 0.5;
    73     beta_lim[0][0].data.F32[7] = 0.5;
    74 
    75     params_min[0][0].data.F32[0] = -1000;
    76     params_min[0][0].data.F32[1] = 0;
    77     params_min[0][0].data.F32[2] = -100;
    78     params_min[0][0].data.F32[3] = -100;
    79     params_min[0][0].data.F32[4] = 0.01;
    80     params_min[0][0].data.F32[5] = 0.01;
    81     params_min[0][0].data.F32[6] = -5.0;
    82     params_min[0][0].data.F32[7] = 0.1;
    83 
    84     params_max[0][0].data.F32[0] = 1e5;
    85     params_max[0][0].data.F32[1] = 1e8;
    86     params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    87     params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    88     params_max[0][0].data.F32[4] = 2.0;
    89     params_max[0][0].data.F32[5] = 2.0;
    90     params_max[0][0].data.F32[6] = +5.0;
    91     params_max[0][0].data.F32[7] = 10.0;
     66    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     67    beta_lim[0][0].data.F32[PM_PAR_FLUX] = 3e6;
     68    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     69    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     70    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     71    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     72    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
     73    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
     74
     75    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     76    params_min[0][0].data.F32[PM_PAR_FLUX] = 0;
     77    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     78    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;
     81    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
     82    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
     83
     84    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     85    params_max[0][0].data.F32[PM_PAR_FLUX] = 1e8;
     86    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     87    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;
     90    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
     91    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
    9292
    9393    return (TRUE);
    … …  
    102102    psF32     *params  = model->params->data.F32;
    103103
    104     params[0] = moments->Sky;
    105     params[1] = moments->Peak - moments->Sky;
    106     params[2] = peak->x;
    107     params[3] = peak->y;
    108     params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
    109     params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
    110     params[6] = 0.0;
    111     params[7] = 1.0;
     104    params[PM_PAR_SKY] = moments->Sky;
     105    params[PM_PAR_FLUX] = 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;
    112112
    113113    return(true);
    … …  
    121121    psF32 *PAR = params->data.F32;
    122122
    123     psF64 A1   = PS_SQR(PAR[4]);
    124     psF64 A2   = PS_SQR(PAR[5]);
    125     psF64 A3   = PS_SQR(PAR[6]);
     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]);
    126126    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    127127    // Area is equivalent to 2 pi sigma^2
    … …  
    135135    float f1, f2;
    136136    for (z = DZ; z < 50; z += DZ) {
    137         f1 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
     137        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    138138        z += DZ;
    139         f2 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
     139        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    140140        norm += f0 + 4*f1 + f2;
    141141        f0 = f2;
    … …  
    143143    norm *= DZ / 3.0;
    144144
    145     psF64 Flux = PAR[1] * Area * norm;
     145    psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm;
    146146
    147147    return(Flux);
    … …  
    158158    if (flux <= 0)
    159159        return (1.0);
    160     if (PAR[1] <= 0)
     160    if (PAR[PM_PAR_FLUX] <= 0)
    161161        return (1.0);
    162     if (flux >= PAR[1])
     162    if (flux >= PAR[PM_PAR_FLUX])
    163163        return (1.0);
    164164
    165165    // if Sx == Sy, sigma = Sx == Sy
    166     psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
    167     psF64 limit = flux / PAR[1];
     166    psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0);
     167    psF64 limit = flux / PAR[PM_PAR_FLUX];
    168168
    169169    # if 0
    … …  
    174174    for (z = 0.0; z < 30.0; z += dz) {
    175175        Nstep ++;
    176         f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
    177         // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z));
     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));
    178178        if (f < limit)
    179179            break;
    … …  
    189189    // choose a z value guaranteed to be beyond our limit
    190190    float z0 = pow((1.0 / limit), (1.0 / 2.25));
    191     float z1 = (1.0 / limit) / PAR[7];
     191    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
    192192    z1 = PS_MAX (z0, z1);
    193193    z0 = 0.0;
    194194
    195     float f0 = 1.0 / (1 + PAR[7]*z0 + pow(z0, 2.25));
    196     float f1 = 1.0 / (1 + PAR[7]*z1 + pow(z1, 2.25));
     195    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
     196    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
    197197    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    198198        z = 0.5*(z0 + z1);
    199         f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
     199        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    200200        // fprintf (stderr, "%f  %f  %f   :   %f  %f  %f\n", f0, f, f1, z0, z, z1);
    201201        if (f > limit) {
    … …  
    224224    psF32 *in  = modelFLT->params->data.F32;
    225225
    226     out[0] = in[0];
    227     out[1] = in[1];
    228     out[2] = in[2];
    229     out[3] = in[3];
    230 
     226    out[PM_PAR_SKY] = in[PM_PAR_SKY];
     227    out[PM_PAR_FLUX] = in[PM_PAR_FLUX];
     228    out[PM_PAR_XPOS] = in[PM_PAR_XPOS];
     229    out[PM_PAR_YPOS] = in[PM_PAR_YPOS];
     230
     231    assert (PM_PAR_YPOS == 4);  // so that we copy the rest of the params
    231232    for (int i = 4; i < 8; i++) {
    232233        psPolynomial2D *poly = psf->params->data[i-4];
    233234        // XXX: Verify this (from EAM change)
    234         //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]);
    235         out[i] = psPolynomial2DEval(poly, out[2], out[3]);
     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]);
    236237    }
    237238    return(true);
    … …  
    248249
    249250    dP = 0;
    250     dP += PS_SQR(dPAR[4] / PAR[4]);
    251     dP += PS_SQR(dPAR[5] / PAR[5]);
     251    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     252    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    252253    dP = sqrt (dP);
    253254
    254255    status = true;
    255256    status &= (dP < 0.5);
    256     status &= (PAR[1] > 0);
    257     status &= ((dPAR[1]/PAR[1]) < 0.5);
     257    status &= (PAR[PM_PAR_FLUX] > 0);
     258    status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5);
    258259
    259260    if (!status)
Note: See TracChangeset for help on using the changeset viewer.