IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 1, 2009, 4:03:58 PM (17 years ago)
Author:
Paul Price
Message:

Merging stack development branch (implement lax and strict PSF model limits) from branches/pap

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk

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

    r20001 r25738  
    1919 *****************************************************************************/
    2020
     21#include <stdio.h>
     22#include <pslib.h>
     23
     24#include "pmMoments.h"
     25#include "pmPeaks.h"
     26#include "pmSource.h"
     27#include "pmModel.h"
     28#include "pmModel_PGAUSS.h"
     29
    2130# define PM_MODEL_FUNC            pmModelFunc_PGAUSS
    2231# define PM_MODEL_FLUX            pmModelFlux_PGAUSS
    … …  
    2736# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS
    2837# define PM_MODEL_FIT_STATUS      pmModelFitStatus_PGAUSS
     38# define PM_MODEL_SET_LIMITS      pmModelSetLimits_PGAUSS
     39
     40// Lax parameter limits
     41static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 };
     42static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     43
     44// Strict parameter limits
     45static float *paramsMinStrict = paramsMinLax;
     46static float *paramsMaxStrict = paramsMaxLax;
     47
     48// Parameter limits to use
     49static float *paramsMinUse = paramsMinLax;
     50static float *paramsMaxUse = paramsMaxLax;
     51static float betaUse[] = { 1000, 3e6, 5, 5, 2.0, 2.0, 0.5 };
     52
     53static bool limitsApply = true;         // Apply limits?
    2954
    3055// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    … …  
    6994bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    7095{
    71     float beta_lim = 0, params_min = 0, params_max = 0;
    72     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     96    if (!limitsApply) {
     97        return true;
     98    }
     99    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    73100
    74101    // we need to calculate the limits for SXY specially
     102    float q2 = NAN;
    75103    if (nParam == PM_PAR_SXY) {
    76         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    77         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    78         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    79         q1 = PS_MAX (0.0, q1);
     104        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     105        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     106        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     107        q1 = (q1 < 0.0) ? 0.0 : q1;
    80108        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    81109        // angle and let f2,f1 fight it out
    82         q2  = 0.5*sqrt (q1);
     110        q2 = 0.5*sqrtf(q1);
    83111    }
    84112
    85113    switch (mode) {
    86     case PS_MINIMIZE_BETA_LIMIT:
    87         switch (nParam) {
    88         case PM_PAR_SKY:
    89             beta_lim = 1000;
    90             break;
    91         case PM_PAR_I0:
    92             beta_lim = 3e6;
    93             break;
    94         case PM_PAR_XPOS:
    95             beta_lim = 5;
    96             break;
    97         case PM_PAR_YPOS:
    98             beta_lim = 5;
    99             break;
    100         case PM_PAR_SXX:
    101             beta_lim = 2.0;
    102             break;
    103         case PM_PAR_SYY:
    104             beta_lim = 2.0;
    105             break;
    106         case PM_PAR_SXY:
    107             beta_lim =  0.5*q2;
    108             break;
    109         default:
    110             psAbort("invalid parameter %d for beta test", nParam);
    111         }
    112         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    113             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    114             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    115                      nParam, beta[nParam], beta_lim);
    116             return false;
    117         }
    118         return true;
    119     case PS_MINIMIZE_PARAM_MIN:
    120         switch (nParam) {
    121         case PM_PAR_SKY:
    122             params_min = -1000;
    123             break;
    124         case PM_PAR_I0:
    125             params_min =  0.01;
    126             break;
    127         case PM_PAR_XPOS:
    128             params_min =  -100;
    129             break;
    130         case PM_PAR_YPOS:
    131             params_min =  -100;
    132             break;
    133         case PM_PAR_SXX:
    134             params_min =   0.5;
    135             break;
    136         case PM_PAR_SYY:
    137             params_min =   0.5;
    138             break;
    139         case PM_PAR_SXY:
    140             params_min =   -q2;
    141             break;
    142         default:
    143             psAbort("invalid parameter %d for param min test", nParam);
    144         }
    145         if (params[nParam] < params_min) {
    146             params[nParam] = params_min;
    147             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    148                      nParam, params[nParam], params_min);
    149             return false;
    150         }
    151         return true;
    152     case PS_MINIMIZE_PARAM_MAX:
    153         switch (nParam) {
    154         case PM_PAR_SKY:
    155             params_max =   1e5;
    156             break;
    157         case PM_PAR_I0:
    158             params_max =   1e8;
    159             break;
    160         case PM_PAR_XPOS:
    161             params_max =   1e4;
    162             break;
    163         case PM_PAR_YPOS:
    164             params_max =   1e4;
    165             break;
    166         case PM_PAR_SXX:
    167             params_max =   100;
    168             break;
    169         case PM_PAR_SYY:
    170             params_max =   100;
    171             break;
    172         case PM_PAR_SXY:
    173             params_max =   +q2;
    174             break;
    175         default:
    176             psAbort("invalid parameter %d for param max test", nParam);
    177         }
    178         if (params[nParam] > params_max) {
    179             params[nParam] = params_max;
    180             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    181                      nParam, params[nParam], params_max);
    182             return false;
    183         }
    184         return true;
    185     default:
     114      case PS_MINIMIZE_BETA_LIMIT: {
     115          psAssert(beta, "Require beta to limit beta");
     116          float limit = betaUse[nParam];
     117          if (nParam == PM_PAR_SXY) {
     118              limit *= q2;
     119          }
     120          if (fabs(beta[nParam]) > fabs(limit)) {
     121              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     122              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     123                      nParam, beta[nParam], limit);
     124              return false;
     125          }
     126          return true;
     127      }
     128      case PS_MINIMIZE_PARAM_MIN: {
     129          psAssert(params, "Require parameters to limit parameters");
     130          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     131          float limit = paramsMinUse[nParam];
     132          if (nParam == PM_PAR_SXY) {
     133              limit *= q2;
     134          }
     135          if (params[nParam] < limit) {
     136              params[nParam] = limit;
     137              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     138                      nParam, params[nParam], limit);
     139              return false;
     140          }
     141          return true;
     142      }
     143      case PS_MINIMIZE_PARAM_MAX: {
     144          psAssert(params, "Require parameters to limit parameters");
     145          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     146          float limit = paramsMaxUse[nParam];
     147          if (nParam == PM_PAR_SXY) {
     148              limit *= q2;
     149          }
     150          if (params[nParam] > limit) {
     151              params[nParam] = limit;
     152              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     153                      nParam, params[nParam], limit);
     154              return false;
     155          }
     156          return true;
     157      }
     158      default:
    186159        psAbort("invalid choice for limits");
    187160    }
    … …  
    189162    return false;
    190163}
     164
    191165
    192166// make an initial guess for parameters
    … …  
    434408}
    435409
     410
     411void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     412{
     413    switch (type) {
     414      case PM_MODEL_LIMITS_NONE:
     415        paramsMinUse = NULL;
     416        paramsMaxUse = NULL;
     417        limitsApply = true;
     418        break;
     419      case PM_MODEL_LIMITS_IGNORE:
     420        paramsMinUse = NULL;
     421        paramsMaxUse = NULL;
     422        limitsApply = false;
     423        break;
     424      case PM_MODEL_LIMITS_LAX:
     425        paramsMinUse = paramsMinLax;
     426        paramsMaxUse = paramsMaxLax;
     427        limitsApply = true;
     428        break;
     429      case PM_MODEL_LIMITS_STRICT:
     430        paramsMinUse = paramsMinStrict;
     431        paramsMaxUse = paramsMaxStrict;
     432        limitsApply = true;
     433        break;
     434      default:
     435        psAbort("Unrecognised model limits type: %x", type);
     436    }
     437    return;
     438}
     439
    436440# undef PM_MODEL_FUNC
    437441# undef PM_MODEL_FLUX
    … …  
    442446# undef PM_MODEL_PARAMS_FROM_PSF
    443447# undef PM_MODEL_FIT_STATUS
     448# undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.