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

    r10195 r10262  
    5858}
    5959
    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);
     60// define the parameter limits
     61// AR_MAX is the maximum allowed axis ratio
     62// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2)
     63# define AR_MAX 20.0
     64# define AR_RATIO 0.99
     65bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
     66{
     67    float beta_lim, params_min, params_max;
     68    float f1, f2, q1, q2;
     69
     70    // we need to calculate the limits for SXY specially
     71    if (nParam == PM_PAR_SXY) {
     72        f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     73        f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     74        q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     75        assert (q1 > 0);
     76        q2  = 0.5*sqrt (q1);
     77    }
     78
     79    switch (mode) {
     80    case PS_MINIMIZE_BETA_LIMIT:
     81        switch (nParam) {
     82        case PM_PAR_SKY:
     83            beta_lim = 1000;
     84            break;
     85        case PM_PAR_I0:
     86            beta_lim = 3e6;
     87            break;
     88        case PM_PAR_XPOS:
     89            beta_lim = 5;
     90            break;
     91        case PM_PAR_YPOS:
     92            beta_lim = 5;
     93            break;
     94        case PM_PAR_SXX:
     95            beta_lim = 0.5;
     96            break;
     97        case PM_PAR_SYY:
     98            beta_lim = 0.5;
     99            break;
     100        case PM_PAR_SXY:
     101            beta_lim =  q2;
     102            break;
     103        default:
     104            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for beta test", nParam);
     105        }
     106        if (fabs(beta[nParam]) > fabs(beta_lim)) {
     107            beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
     108            return false;
     109        }
     110        return true;
     111    case PS_MINIMIZE_PARAM_MIN:
     112        switch (nParam) {
     113        case PM_PAR_SKY:
     114            params_min = -1000;
     115            break;
     116        case PM_PAR_I0:
     117            params_min =     0;
     118            break;
     119        case PM_PAR_XPOS:
     120            params_min =  -100;
     121            break;
     122        case PM_PAR_YPOS:
     123            params_min =  -100;
     124            break;
     125        case PM_PAR_SXX:
     126            params_min =   0.5;
     127            break;
     128        case PM_PAR_SYY:
     129            params_min =   0.5;
     130            break;
     131        case PM_PAR_SXY:
     132            params_min =   -q2;
     133            break;
     134        default:
     135            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param min test", nParam);
     136        }
     137        if (params[nParam] < params_min) {
     138            params[nParam] = params_min;
     139            return false;
     140        }
     141        return true;
     142    case PS_MINIMIZE_PARAM_MAX:
     143        switch (nParam) {
     144        case PM_PAR_SKY:
     145            params_max =   1e5;
     146            break;
     147        case PM_PAR_I0:
     148            params_max =   1e8;
     149            break;
     150        case PM_PAR_XPOS:
     151            params_max =   1e4;
     152            break;
     153        case PM_PAR_YPOS:
     154            params_max =   1e4;
     155            break;
     156        case PM_PAR_SXX:
     157            params_max =   100;
     158            break;
     159        case PM_PAR_SYY:
     160            params_max =   100;
     161            break;
     162        case PM_PAR_SXY:
     163            params_max =   +q2;
     164            break;
     165        default:
     166            psAbort ("psModules.pmModel_GAUSS", "invalid parameter %d for param max test", nParam);
     167        }
     168        if (params[nParam] > params_max) {
     169            params[nParam] = params_max;
     170            return false;
     171        }
     172        return true;
     173    default:
     174        psAbort ("psModules.pmModel_GAUSS", "invalid choice for limits");
     175    }
     176    psAbort ("psModules.pmModel_GAUSS", "should not reach here");
     177    return false;
    92178}
    93179
    … …  
    98184    psF32     *PAR     = model->params->data.F32;
    99185
    100     # if (0)
    101 
    102         psEllipseMoments emoments;
     186    psEllipseMoments emoments;
    103187    emoments.x2 = moments->Sx;
    104188    emoments.y2 = moments->Sx;
    105189    emoments.xy = moments->Sxy;
    106190
    107     psEllipseAxes axes = psEllipseMomentsToAxes (emoments);
     191    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
    108192    psEllipseShape shape = psEllipseAxesToShape (axes);
    109     # endif
    110193
    111194    PAR[PM_PAR_SKY]  = moments->Sky;
    … …  
    113196    PAR[PM_PAR_XPOS] = moments->x; // XXX use peak->xf, peak->yf?
    114197    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 
    119     # if (0)
    120 
    121         PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
     198    PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
    122199    PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy);
    123200    PAR[PM_PAR_SXY] = shape.sxy;
    124     # endif
    125 
    126201    return(true);
    127202}
    … …  
    139214
    140215    // Area is equivalent to 2 pi sigma^2
    141     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     216    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    142217    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    143218
    … …  
    183258
    184259    // this estimates the radius assuming f(z) is roughly exp(-z)
    185     psEllipseAxes axes = psEllipseShapeToAxes (shape);
     260    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    186261    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
    187262
Note: See TracChangeset for help on using the changeset viewer.