IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 12:11:34 PM (17 years ago)
Author:
eugene
Message:

merge changes from trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/objects/models/pmModel_GAUSS.c

    r25683 r25745  
    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_GAUSS.h"
     29
    2130# define PM_MODEL_FUNC            pmModelFunc_GAUSS
    2231# define PM_MODEL_FLUX            pmModelFlux_GAUSS
    … …  
    2736# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_GAUSS
    2837# define PM_MODEL_FIT_STATUS      pmModelFitStatus_GAUSS
     38# define PM_MODEL_SET_LIMITS      pmModelSetLimits_GAUSS
     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)
    … …  
    7095bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    7196{
    72     float beta_lim = 0, params_min = 0, params_max = 0;
    73     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     97    if (!limitsApply) {
     98        return true;
     99    }
     100    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    74101
    75102    // we need to calculate the limits for SXY specially
     103    float q2 = NAN;
    76104    if (nParam == PM_PAR_SXY) {
    77         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    78         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    79         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    80         q1 = PS_MAX (0.0, q1);
     105        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     106        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     107        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     108        q1 = (q1 < 0.0) ? 0.0 : q1;
    81109        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    82110        // angle and let f2,f1 fight it out
    83         q2  = 0.5*sqrt (q1);
     111        q2 = 0.5*sqrtf(q1);
    84112    }
    85113
    86114    switch (mode) {
    87     case PS_MINIMIZE_BETA_LIMIT:
    88         switch (nParam) {
    89         case PM_PAR_SKY:
    90             beta_lim = 1000;
    91             break;
    92         case PM_PAR_I0:
    93             beta_lim = 3e6;
    94             break;
    95         case PM_PAR_XPOS:
    96             beta_lim = 5;
    97             break;
    98         case PM_PAR_YPOS:
    99             beta_lim = 5;
    100             break;
    101         case PM_PAR_SXX:
    102             beta_lim = 2.0;
    103             break;
    104         case PM_PAR_SYY:
    105             beta_lim = 2.0;
    106             break;
    107         case PM_PAR_SXY:
    108             beta_lim =  0.5*q2;
    109             break;
    110         default:
    111             psAbort("invalid parameter %d for beta test", nParam);
    112         }
    113         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    114             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    115             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    116                      nParam, beta[nParam], beta_lim);
    117             return false;
    118         }
    119         return true;
    120     case PS_MINIMIZE_PARAM_MIN:
    121         switch (nParam) {
    122         case PM_PAR_SKY:
    123             params_min = -1000;
    124             break;
    125         case PM_PAR_I0:
    126             params_min =   0.01;
    127             break;
    128         case PM_PAR_XPOS:
    129             params_min =  -100;
    130             break;
    131         case PM_PAR_YPOS:
    132             params_min =  -100;
    133             break;
    134         case PM_PAR_SXX:
    135             params_min =   0.5;
    136             break;
    137         case PM_PAR_SYY:
    138             params_min =   0.5;
    139             break;
    140         case PM_PAR_SXY:
    141             params_min =   -q2;
    142             break;
    143         default:
    144             psAbort("invalid parameter %d for param min test", nParam);
    145         }
    146         if (params[nParam] < params_min) {
    147             params[nParam] = params_min;
    148             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    149                      nParam, params[nParam], params_min);
    150             return false;
    151         }
    152         return true;
    153     case PS_MINIMIZE_PARAM_MAX:
    154         switch (nParam) {
    155         case PM_PAR_SKY:
    156             params_max =   1e5;
    157             break;
    158         case PM_PAR_I0:
    159             params_max =   1e8;
    160             break;
    161         case PM_PAR_XPOS:
    162             params_max =   1e4;
    163             break;
    164         case PM_PAR_YPOS:
    165             params_max =   1e4;
    166             break;
    167         case PM_PAR_SXX:
    168             params_max =   100;
    169             break;
    170         case PM_PAR_SYY:
    171             params_max =   100;
    172             break;
    173         case PM_PAR_SXY:
    174             params_max =   +q2;
    175             break;
    176         default:
    177             psAbort("invalid parameter %d for param max test", nParam);
    178         }
    179         if (params[nParam] > params_max) {
    180             params[nParam] = params_max;
    181             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    182                      nParam, params[nParam], params_max);
    183             return false;
    184         }
    185         return true;
     115      case PS_MINIMIZE_BETA_LIMIT: {
     116          psAssert(beta, "Require beta to limit beta");
     117          float limit = betaUse[nParam];
     118          if (nParam == PM_PAR_SXY) {
     119              limit *= q2;
     120          }
     121          if (fabs(beta[nParam]) > fabs(limit)) {
     122              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     123              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     124                      nParam, beta[nParam], limit);
     125              return false;
     126          }
     127          return true;
     128      }
     129      case PS_MINIMIZE_PARAM_MIN: {
     130          psAssert(params, "Require parameters to limit parameters");
     131          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     132          float limit = paramsMinUse[nParam];
     133          if (nParam == PM_PAR_SXY) {
     134              limit *= q2;
     135          }
     136          if (params[nParam] < limit) {
     137              params[nParam] = limit;
     138              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     139                      nParam, params[nParam], limit);
     140              return false;
     141          }
     142          return true;
     143      }
     144      case PS_MINIMIZE_PARAM_MAX: {
     145          psAssert(params, "Require parameters to limit parameters");
     146          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     147          float limit = paramsMaxUse[nParam];
     148          if (nParam == PM_PAR_SXY) {
     149              limit *= q2;
     150          }
     151          if (params[nParam] > limit) {
     152              params[nParam] = limit;
     153              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     154                      nParam, params[nParam], limit);
     155              return false;
     156          }
     157          return true;
     158      }
    186159    default:
    187160        psAbort("invalid choice for limits");
    … …  
    384357}
    385358
     359void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     360{
     361    switch (type) {
     362      case PM_MODEL_LIMITS_NONE:
     363        paramsMinUse = NULL;
     364        paramsMaxUse = NULL;
     365        limitsApply = true;
     366        break;
     367      case PM_MODEL_LIMITS_IGNORE:
     368        paramsMinUse = NULL;
     369        paramsMaxUse = NULL;
     370        limitsApply = false;
     371        break;
     372      case PM_MODEL_LIMITS_LAX:
     373        paramsMinUse = paramsMinLax;
     374        paramsMaxUse = paramsMaxLax;
     375        limitsApply = true;
     376        break;
     377      case PM_MODEL_LIMITS_STRICT:
     378        paramsMinUse = paramsMinStrict;
     379        paramsMaxUse = paramsMaxStrict;
     380        limitsApply = true;
     381        break;
     382      default:
     383        psAbort("Unrecognised model limits type: %x", type);
     384    }
     385    return;
     386}
     387
    386388# undef PM_MODEL_FUNC
    387389# undef PM_MODEL_FLUX
    … …  
    392394# undef PM_MODEL_PARAMS_FROM_PSF
    393395# undef PM_MODEL_FIT_STATUS
     396# undef PM_MODEL_SET_LIMITS
Note: See TracChangeset for help on using the changeset viewer.