IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10262


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

Location:
trunk/psModules/src/objects/models
Files:
3 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}

  • 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

  • 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.