IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 28, 2006, 4:55:24 PM (19 years ago)
Author:
magnier
Message:

changed limits to new dynamic limits function, added ellipse parameter limits for sxy

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r10078 r10262  
    4343    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
    4444    // for negative values of z
    45     if (z < 0)
    46         z = 0;
     45    // XXX use an assert here to force the elliptical parameters to be correctly determined
     46    // if (z < 0) z = 0;
     47    assert (z > 0);
    4748
    4849    psF32 zp = pow(z,1.25);
    … …  
    7273}
    7374
    74 bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
    75 {
    76 
    77     *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
    78     *params_min = psVectorAlloc (8, PS_TYPE_F32);
    79     *params_max = psVectorAlloc (8, PS_TYPE_F32);
    80 
    81     beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
    82     beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
    83     beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
    84     beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
    85     beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
    86     beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
    87     beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
    88     beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
    89 
    90     params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
    91     params_min[0][0].data.F32[PM_PAR_I0] = 0;
    92     params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
    93     params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
    94     params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
    95     params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
    96     params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    97     params_min[0][0].data.F32[PM_PAR_7] = 0.1;
    98 
    99     params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
    100     params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
    101     params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
    102     params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
    103     params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
    104     params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
    105     params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    106     params_max[0][0].data.F32[PM_PAR_7] = 10.0;
    107 
    108     return (TRUE);
    109 }
     75// define the parameter limits
     76// AR_MAX is the maximum allowed axis ratio
     77// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2)
     78# define AR_MAX 20.0
     79# define AR_RATIO 0.99
     80bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
     81{
     82    float beta_lim, params_min, params_max;
     83    float f1, f2, q1, q2;
     84
     85    // we need to calculate the limits for SXY specially
     86    if (nParam == PM_PAR_SXY) {
     87        f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     88        f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     89        q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     90        assert (q1 > 0);
     91        q2  = 0.5*sqrt (q1);
     92    }
     93
     94    switch (mode) {
     95    case PS_MINIMIZE_BETA_LIMIT:
     96        switch (nParam) {
     97        case PM_PAR_SKY:
     98            beta_lim = 1000;
     99            break;
     100        case PM_PAR_I0:
     101            beta_lim = 3e6;
     102            break;
     103        case PM_PAR_XPOS:
     104            beta_lim = 5;
     105            break;
     106        case PM_PAR_YPOS:
     107            beta_lim = 5;
     108            break;
     109        case PM_PAR_SXX:
     110            beta_lim = 0.5;
     111            break;
     112        case PM_PAR_SYY:
     113            beta_lim = 0.5;
     114            break;
     115        case PM_PAR_SXY:
     116            beta_lim =  q2;
     117            break;
     118        case PM_PAR_7:
     119            beta_lim = 0.5;
     120            break;
     121        default:
     122            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam);
     123        }
     124        if (fabs(beta[nParam]) > fabs(beta_lim)) {
     125            beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
     126            return false;
     127        }
     128        return true;
     129    case PS_MINIMIZE_PARAM_MIN:
     130        switch (nParam) {
     131        case PM_PAR_SKY:
     132            params_min = -1000;
     133            break;
     134        case PM_PAR_I0:
     135            params_min =     0;
     136            break;
     137        case PM_PAR_XPOS:
     138            params_min =  -100;
     139            break;
     140        case PM_PAR_YPOS:
     141            params_min =  -100;
     142            break;
     143        case PM_PAR_SXX:
     144            params_min =   0.5;
     145            break;
     146        case PM_PAR_SYY:
     147            params_min =   0.5;
     148            break;
     149        case PM_PAR_SXY:
     150            params_min =  -q2;
     151            break;
     152        case PM_PAR_7:
     153            params_min =   0.1;
     154            break;
     155        default:
     156            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam);
     157        }
     158        if (params[nParam] < params_min) {
     159            params[nParam] = params_min;
     160            return false;
     161        }
     162        return true;
     163    case PS_MINIMIZE_PARAM_MAX:
     164        switch (nParam) {
     165        case PM_PAR_SKY:
     166            params_max =   1e5;
     167            break;
     168        case PM_PAR_I0:
     169            params_max =   1e8;
     170            break;
     171        case PM_PAR_XPOS:
     172            params_max =   1e4;
     173            break;
     174        case PM_PAR_YPOS:
     175            params_max =   1e4;
     176            break;
     177        case PM_PAR_SXX:
     178            params_max =   100;
     179            break;
     180        case PM_PAR_SYY:
     181            params_max =   100;
     182            break;
     183        case PM_PAR_SXY:
     184            params_max =  +q2;
     185            break;
     186        case PM_PAR_7:
     187            params_max =  10.0;
     188            break;
     189        default:
     190            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam);
     191        }
     192        if (params[nParam] > params_max) {
     193            params[nParam] = params_max;
     194            return false;
     195        }
     196        return true;
     197    default:
     198        psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits");
     199    }
     200    psAbort ("psModules.pmModel_GAUSS", "should not reach here");
     201    return false;
     202}
     203
    110204
    111205// make an initial guess for parameters
    112206bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    113207{
    114 
    115208    pmMoments *moments = source->moments;
    116209    pmPeak    *peak    = source->peak;
    117210    psF32     *PAR  = model->params->data.F32;
     211
     212    psEllipseMoments emoments;
     213    emoments.x2 = moments->Sx;
     214    emoments.y2 = moments->Sy;
     215    emoments.xy = moments->Sxy;
     216
     217    // force the axis ratio to be < 20.0
     218    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
     219    psEllipseShape shape = psEllipseAxesToShape (axes);
    118220
    119221    PAR[PM_PAR_SKY]  = moments->Sky;
    … …  
    121223    PAR[PM_PAR_XPOS] = peak->x;
    122224    PAR[PM_PAR_YPOS] = peak->y;
    123     PAR[PM_PAR_SXX]  = PS_MAX(0.5, moments->Sx);
    124     PAR[PM_PAR_SYY]  = PS_MAX(0.5, moments->Sy);
    125     PAR[PM_PAR_SXY]  = 0.0;
     225    PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
     226    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
     227    PAR[PM_PAR_SXY]  = shape.sxy;
    126228    PAR[PM_PAR_7]    = 1.0;
    127229
    … …  
    141243
    142244    // Area is equivalent to 2 pi sigma^2
    143     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     245    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    144246    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    145247
    … …  
    186288    shape.sxy = PAR[PM_PAR_SXY];
    187289
    188     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     290    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    189291    psF64 sigma = axes.major;
    190292
Note: See TracChangeset for help on using the changeset viewer.