IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 28, 2006, 4:55:24 PM (20 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_GAUSS.c

    r10181 r10262  
    3333    psF32 *PAR = params->data.F32;
    3434
    35     // XXX this is fitting sqrt(2)/sigma_x, sqrt(2)/sigma_y
    3635    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
    3736    psF32 Y  = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS];
    … …  
    5857
    5958// define the parameter limits
    60 bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max)
    61 {
    62 
    63     *beta_lim   = psVectorAlloc (7, PS_TYPE_F32);
    64     *params_min = psVectorAlloc (7, PS_TYPE_F32);
    65     *params_max = psVectorAlloc (7, PS_TYPE_F32);
    66 
    67     beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
    68     beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
    69     beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
    70     beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
    71     beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
    72     beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
    73     beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
    74 
    75     params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
    76     params_min[0][0].data.F32[PM_PAR_I0] = 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.5;
    80     params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
    81     params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    82 
    83     params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
    84     params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
    85     params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
    86     params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
    87     params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
    88     params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
    89     params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    90 
    91     return (TRUE);
     59// AR_MAX is the maximum allowed axis ratio
     60// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2)
     61# define AR_MAX 20.0
     62# define AR_RATIO 0.99
     63bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
     64{
     65    float beta_lim, params_min, params_max;
     66    float f1, f2, q1, q2;
     67
     68    // we need to calculate the limits for SXY specially
     69    if (nParam == PM_PAR_SXY) {
     70        f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     71        f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     72        q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     73        assert (q1 > 0);
     74        q2  = 0.5*sqrt (q1);
     75    }
     76
     77    switch (mode) {
     78    case PS_MINIMIZE_BETA_LIMIT:
     79        switch (nParam) {
     80        case PM_PAR_SKY:
     81            beta_lim = 1000;
     82            break;
     83        case PM_PAR_I0:
     84            beta_lim = 3e6;
     85            break;
     86        case PM_PAR_XPOS:
     87            beta_lim = 5;
     88            break;
     89        case PM_PAR_YPOS:
     90            beta_lim = 5;
     91            break;
     92        case PM_PAR_SXX:
     93            beta_lim = 0.5;
     94            break;
     95        case PM_PAR_SYY:
     96            beta_lim = 0.5;
     97            break;
     98        case PM_PAR_SXY:
     99            beta_lim =  q2;
     100            break;
     101        default:
     102            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam);
     103        }
     104        if (fabs(beta[nParam]) > fabs(beta_lim)) {
     105            beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
     106            return false;
     107        }
     108        return true;
     109    case PS_MINIMIZE_PARAM_MIN:
     110        switch (nParam) {
     111        case PM_PAR_SKY:
     112            params_min = -1000;
     113            break;
     114        case PM_PAR_I0:
     115            params_min =     0;
     116            break;
     117        case PM_PAR_XPOS:
     118            params_min =  -100;
     119            break;
     120        case PM_PAR_YPOS:
     121            params_min =  -100;
     122            break;
     123        case PM_PAR_SXX:
     124            params_min =   0.5;
     125            break;
     126        case PM_PAR_SYY:
     127            params_min =   0.5;
     128            break;
     129        case PM_PAR_SXY:
     130            params_min =   -q2;
     131            break;
     132        default:
     133            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam);
     134        }
     135        if (params[nParam] < params_min) {
     136            params[nParam] = params_min;
     137            return false;
     138        }
     139        return true;
     140    case PS_MINIMIZE_PARAM_MAX:
     141        switch (nParam) {
     142        case PM_PAR_SKY:
     143            params_max =   1e5;
     144            break;
     145        case PM_PAR_I0:
     146            params_max =   1e8;
     147            break;
     148        case PM_PAR_XPOS:
     149            params_max =   1e4;
     150            break;
     151        case PM_PAR_YPOS:
     152            params_max =   1e4;
     153            break;
     154        case PM_PAR_SXX:
     155            params_max =   100;
     156            break;
     157        case PM_PAR_SYY:
     158            params_max =   100;
     159            break;
     160        case PM_PAR_SXY:
     161            params_max =   +q2;
     162            break;
     163        default:
     164            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam);
     165        }
     166        if (params[nParam] > params_max) {
     167            params[nParam] = params_max;
     168            return false;
     169        }
     170        return true;
     171    default:
     172        psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits");
     173    }
     174    psAbort ("psModules.pmModel_GAUSS", "should not reach here");
     175    return false;
    92176}
    93177
    … …  
    98182    psF32     *PAR  = model->params->data.F32;
    99183
    100     # if (0)
    101 
    102         psEllipseMoments emoments;
     184    psEllipseMoments emoments;
    103185    emoments.x2 = moments->Sx;
    104     emoments.y2 = moments->Sx;
     186    emoments.y2 = moments->Sy;
    105187    emoments.xy = moments->Sxy;
    106188
    107     psEllipseAxes axes = psEllipseMomentsToAxes (emoments);
     189    // force the axis ratio to be < 20.0
     190    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
    108191    psEllipseShape shape = psEllipseAxesToShape (axes);
    109     # endif
    110 
    111     PAR[PM_PAR_SKY] = moments->Sky;
    112     PAR[PM_PAR_I0] = moments->Peak - moments->Sky;
     192
     193    PAR[PM_PAR_SKY]  = moments->Sky;
     194    PAR[PM_PAR_I0]   = moments->Peak - moments->Sky;
    113195    PAR[PM_PAR_XPOS] = moments->x;
    114196    PAR[PM_PAR_YPOS] = moments->y;
    115     PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*moments->Sx);
    116     PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*moments->Sy);
    117     PAR[PM_PAR_SXY] = 0.0;
    118     # if (0)
    119 
    120         PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
     197    PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
    121198    PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy);
    122199    PAR[PM_PAR_SXY] = shape.sxy;
    123     # endif
    124 
    125200    return(true);
    126201}
    … …  
    138213
    139214    // Area is equivalent to 2 pi sigma^2
    140     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     215    // axes ratio < 20
     216    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    141217    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    142218
    … …  
    165241    shape.sxy = PAR[PM_PAR_SXY];
    166242
    167     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     243    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    168244    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
    169245    return (radius);
    … …  
    194270    out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
    195271
     272    // apply the model limits here: this truncates excessive extrapolation
     273    for (int i = 0; i < psf->params_NEW->n; i++) {
     274        // apply the limits to all components or just the psf-model parameters?
     275        if (psf->params_NEW->data[i] == NULL)
     276            continue;
     277        PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, out, NULL);
     278        PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, out, NULL);
     279    }
    196280    return(true);
    197281}
Note: See TracChangeset for help on using the changeset viewer.