IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25738


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:
35 edited
6 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/archive/ducttape

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/catalyst

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/catalyst/2007.0-specs

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/catalyst/2008.0-specs

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/catalyst/2008.0-specs/hardened

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/doc/misc/docgen.pl

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/hardware

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/ippScripts/scripts/magic_destreak_revert.pl

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/ippTests/tap

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/ippTools/share/disttool_revertcomponent.sql

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/ippTools/share/disttool_revertrun.sql

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/ppStack/src/ppStackLoop.c

    r23576 r25738  
    8080        return false;
    8181    }
    82     psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
    83              psTimerClear("PPSTACK_STEPS"));
    8482    ppStackMemDump("convolve");
    85     psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
     83    psLogMsg("ppStack", PS_LOG_INFO, "Stage 3: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
    8684
    8785
    … …  
    9492        return false;
    9593    }
    96     psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
     94    psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
    9795    ppStackMemDump("reject");
    9896
    … …  
    106104        return false;
    107105    }
    108     psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
     106    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
    109107    ppStackMemDump("final");
    110108
    … …  
    118116        return false;
    119117    }
    120     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     118    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
    121119    ppStackMemDump("cleanup");
    122120
    … …  
    131129        return false;
    132130    }
    133     psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
     131    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
    134132    ppStackMemDump("photometry");
    135133
    … …  
    142140        return false;
    143141    }
    144     psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
     142    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
    145143    ppStackMemDump("finish");
    146144

  • trunk/psModules/src/camera/pmReadoutFake.c

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/psModules/src/config/pmConfigMask.c

    r23498 r25738  
    2222    // Non-astronomical structures
    2323    { "CR",        NULL,       0x08, true  }, // Pixel contains a cosmic ray
    24     { "SPIKE",     NULL,       0x08, true  }, // Pixel contains a diffraction spike
    25     { "GHOST",     NULL,       0x08, true  }, // Pixel contains an optical ghost
    26     { "STREAK",    NULL,       0x08, true  }, // Pixel contains streak data
    27     { "STARCORE",  NULL,       0x08, true  }, // Pixel contains a bright star core
     24    { "SPIKE",     NULL,       0x08, false  }, // Pixel contains a diffraction spike
     25    { "GHOST",     NULL,       0x08, false  }, // Pixel contains an optical ghost
     26    { "STREAK",    NULL,       0x08, false  }, // Pixel contains streak data
     27    { "STARCORE",  NULL,       0x08, false  }, // Pixel contains a bright star core
    2828    // Effects of convolution and interpolation
    2929    { "CONV.BAD",  NULL,       0x02, true  }, // Pixel is bad after convolution with a bad pixel

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r25491 r25738  
    149149        }
    150150
     151        // Test PSF
     152        {
     153            bool goodPSF = true;                                                                // Good PSF?
     154            pmModelClassSetLimits(PM_MODEL_LIMITS_IGNORE);
     155            pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, PEAK_FLUX); // Test model
     156            model->modelSetLimits(PM_MODEL_LIMITS_STRICT);
     157            for (int j = 0; j < model->params->n && goodPSF; j++) {
     158                if (!model->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) ||
     159                    !model->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) {
     160                    goodPSF = false;
     161                }
     162            }
     163            psFree(model);
     164            if (!goodPSF) {
     165                psWarning("PSF %d is bad --- not including in envelope calculation.", i);
     166                continue;
     167            }
     168        }
     169
    151170        pmResiduals *resid = psf->residuals;// PSF residuals
    152171        psf->residuals = NULL;
    153172        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
    154                                       NAN, radius, true, true)) {
     173                                      NAN, radius, true, false)) {
    155174            psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
    156175            psFree(envelope);

  • trunk/psModules/src/objects/Makefile.am

    r25383 r25738  
    5454        pmGrowthCurve.c \
    5555        pmSourceMatch.c \
    56         pmDetEff.c
    57 
    58 EXTRA_DIST = \
     56        pmDetEff.c \
    5957        models/pmModel_GAUSS.c \
    6058        models/pmModel_PGAUSS.c \
     59        models/pmModel_PS1_V1.c \
    6160        models/pmModel_QGAUSS.c \
    62         models/pmModel_SGAUSS.c \
    6361        models/pmModel_RGAUSS.c \
    6462        models/pmModel_SERSIC.c
    … …  
    9290        pmGrowthCurve.h \
    9391        pmSourceMatch.h \
    94         pmDetEff.h
     92        pmDetEff.h \
     93        models/pmModel_GAUSS.h \
     94        models/pmModel_PGAUSS.h \
     95        models/pmModel_PS1_V1.h \
     96        models/pmModel_QGAUSS.h \
     97        models/pmModel_RGAUSS.h \
     98        models/pmModel_SERSIC.h
    9599
    96100CLEANFILES = *~

  • trunk/psModules/src/objects/models/pmModel_GAUSS.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_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)
    … …  
    6893bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    6994{
    70     float beta_lim = 0, params_min = 0, params_max = 0;
    71     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     95    if (!limitsApply) {
     96        return true;
     97    }
     98    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    7299
    73100    // we need to calculate the limits for SXY specially
     101    float q2 = NAN;
    74102    if (nParam == PM_PAR_SXY) {
    75         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    76         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    77         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    78         q1 = PS_MAX (0.0, q1);
     103        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     104        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     105        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     106        q1 = (q1 < 0.0) ? 0.0 : q1;
    79107        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    80108        // angle and let f2,f1 fight it out
    81         q2  = 0.5*sqrt (q1);
     109        q2 = 0.5*sqrtf(q1);
    82110    }
    83111
    84112    switch (mode) {
    85     case PS_MINIMIZE_BETA_LIMIT:
    86         switch (nParam) {
    87         case PM_PAR_SKY:
    88             beta_lim = 1000;
    89             break;
    90         case PM_PAR_I0:
    91             beta_lim = 3e6;
    92             break;
    93         case PM_PAR_XPOS:
    94             beta_lim = 5;
    95             break;
    96         case PM_PAR_YPOS:
    97             beta_lim = 5;
    98             break;
    99         case PM_PAR_SXX:
    100             beta_lim = 2.0;
    101             break;
    102         case PM_PAR_SYY:
    103             beta_lim = 2.0;
    104             break;
    105         case PM_PAR_SXY:
    106             beta_lim =  0.5*q2;
    107             break;
    108         default:
    109             psAbort("invalid parameter %d for beta test", nParam);
    110         }
    111         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    112             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    113             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    114                      nParam, beta[nParam], beta_lim);
    115             return false;
    116         }
    117         return true;
    118     case PS_MINIMIZE_PARAM_MIN:
    119         switch (nParam) {
    120         case PM_PAR_SKY:
    121             params_min = -1000;
    122             break;
    123         case PM_PAR_I0:
    124             params_min =   0.01;
    125             break;
    126         case PM_PAR_XPOS:
    127             params_min =  -100;
    128             break;
    129         case PM_PAR_YPOS:
    130             params_min =  -100;
    131             break;
    132         case PM_PAR_SXX:
    133             params_min =   0.5;
    134             break;
    135         case PM_PAR_SYY:
    136             params_min =   0.5;
    137             break;
    138         case PM_PAR_SXY:
    139             params_min =   -q2;
    140             break;
    141         default:
    142             psAbort("invalid parameter %d for param min test", nParam);
    143         }
    144         if (params[nParam] < params_min) {
    145             params[nParam] = params_min;
    146             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    147                      nParam, params[nParam], params_min);
    148             return false;
    149         }
    150         return true;
    151     case PS_MINIMIZE_PARAM_MAX:
    152         switch (nParam) {
    153         case PM_PAR_SKY:
    154             params_max =   1e5;
    155             break;
    156         case PM_PAR_I0:
    157             params_max =   1e8;
    158             break;
    159         case PM_PAR_XPOS:
    160             params_max =   1e4;
    161             break;
    162         case PM_PAR_YPOS:
    163             params_max =   1e4;
    164             break;
    165         case PM_PAR_SXX:
    166             params_max =   100;
    167             break;
    168         case PM_PAR_SYY:
    169             params_max =   100;
    170             break;
    171         case PM_PAR_SXY:
    172             params_max =   +q2;
    173             break;
    174         default:
    175             psAbort("invalid parameter %d for param max test", nParam);
    176         }
    177         if (params[nParam] > params_max) {
    178             params[nParam] = params_max;
    179             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    180                      nParam, params[nParam], params_max);
    181             return false;
    182         }
    183         return true;
     113      case PS_MINIMIZE_BETA_LIMIT: {
     114          psAssert(beta, "Require beta to limit beta");
     115          float limit = betaUse[nParam];
     116          if (nParam == PM_PAR_SXY) {
     117              limit *= q2;
     118          }
     119          if (fabs(beta[nParam]) > fabs(limit)) {
     120              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     121              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     122                      nParam, beta[nParam], limit);
     123              return false;
     124          }
     125          return true;
     126      }
     127      case PS_MINIMIZE_PARAM_MIN: {
     128          psAssert(params, "Require parameters to limit parameters");
     129          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     130          float limit = paramsMinUse[nParam];
     131          if (nParam == PM_PAR_SXY) {
     132              limit *= q2;
     133          }
     134          if (params[nParam] < limit) {
     135              params[nParam] = limit;
     136              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     137                      nParam, params[nParam], limit);
     138              return false;
     139          }
     140          return true;
     141      }
     142      case PS_MINIMIZE_PARAM_MAX: {
     143          psAssert(params, "Require parameters to limit parameters");
     144          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     145          float limit = paramsMaxUse[nParam];
     146          if (nParam == PM_PAR_SXY) {
     147              limit *= q2;
     148          }
     149          if (params[nParam] > limit) {
     150              params[nParam] = limit;
     151              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     152                      nParam, params[nParam], limit);
     153              return false;
     154          }
     155          return true;
     156      }
    184157    default:
    185158        psAbort("invalid choice for limits");
    … …  
    388361}
    389362
     363void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     364{
     365    switch (type) {
     366      case PM_MODEL_LIMITS_NONE:
     367        paramsMinUse = NULL;
     368        paramsMaxUse = NULL;
     369        limitsApply = true;
     370        break;
     371      case PM_MODEL_LIMITS_IGNORE:
     372        paramsMinUse = NULL;
     373        paramsMaxUse = NULL;
     374        limitsApply = false;
     375        break;
     376      case PM_MODEL_LIMITS_LAX:
     377        paramsMinUse = paramsMinLax;
     378        paramsMaxUse = paramsMaxLax;
     379        limitsApply = true;
     380        break;
     381      case PM_MODEL_LIMITS_STRICT:
     382        paramsMinUse = paramsMinStrict;
     383        paramsMaxUse = paramsMaxStrict;
     384        limitsApply = true;
     385        break;
     386      default:
     387        psAbort("Unrecognised model limits type: %x", type);
     388    }
     389    return;
     390}
     391
    390392# undef PM_MODEL_FUNC
    391393# undef PM_MODEL_FLUX
    … …  
    396398# undef PM_MODEL_PARAMS_FROM_PSF
    397399# undef PM_MODEL_FIT_STATUS
     400# undef PM_MODEL_SET_LIMITS

  • 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

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

    r23962 r25738  
    2020   *****************************************************************************/
    2121
     22#include <stdio.h>
     23#include <pslib.h>
     24
     25#include "pmMoments.h"
     26#include "pmPeaks.h"
     27#include "pmSource.h"
     28#include "pmModel.h"
     29#include "pmModel_PS1_V1.h"
     30
    2231# define PM_MODEL_FUNC            pmModelFunc_PS1_V1
    2332# define PM_MODEL_FLUX            pmModelFlux_PS1_V1
    … …  
    2837# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PS1_V1
    2938# define PM_MODEL_FIT_STATUS      pmModelFitStatus_PS1_V1
     39# define PM_MODEL_SET_LIMITS      pmModelSetLimits_PS1_V1
    3040
    3141# define ALPHA   1.666
    3242# define ALPHA_M 0.666
     43
     44// Lax parameter limits
     45static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
     46static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     47
     48// Strict parameter limits
     49// k = PAR_7 < 0 is very undesirable (big divot in the middle)
     50static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
     51static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     52
     53// Parameter limits to use
     54static float *paramsMinUse = paramsMinLax;
     55static float *paramsMaxUse = paramsMaxLax;
     56static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     57
     58static bool limitsApply = true;         // Apply limits?
     59
    3360
    3461psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    84111bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    85112{
    86     float beta_lim = 0, params_min = 0, params_max = 0;
    87     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     113    if (!limitsApply) {
     114        return true;
     115    }
     116    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    88117
    89118    // we need to calculate the limits for SXY specially
     119    float q2 = NAN;
    90120    if (nParam == PM_PAR_SXY) {
    91         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    92         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    93         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     121        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     122        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     123        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    94124        q1 = (q1 < 0.0) ? 0.0 : q1;
    95125        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    96126        // angle and let f2,f1 fight it out
    97         q2  = 0.5*sqrt (q1);
     127        q2 = 0.5*sqrtf(q1);
    98128    }
    99129
    100130    switch (mode) {
    101     case PS_MINIMIZE_BETA_LIMIT:
    102         switch (nParam) {
    103         case PM_PAR_SKY:
    104             beta_lim = 1000;
    105             break;
    106         case PM_PAR_I0:
    107             beta_lim = 3e6;
    108             break;
    109         case PM_PAR_XPOS:
    110             beta_lim = 5;
    111             break;
    112         case PM_PAR_YPOS:
    113             beta_lim = 5;
    114             break;
    115         case PM_PAR_SXX:
    116             beta_lim = 1.0;
    117             break;
    118         case PM_PAR_SYY:
    119             beta_lim = 1.0;
    120             break;
    121         case PM_PAR_SXY:
    122             beta_lim =  0.5*q2;
    123             break;
    124         case PM_PAR_7:
    125             beta_lim = 2.0;
    126             break;
    127         default:
    128             psAbort("invalid parameter %d for beta test", nParam);
    129         }
    130         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    131             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    132             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    133                      nParam, beta[nParam], beta_lim);
    134             return false;
    135         }
    136         return true;
    137     case PS_MINIMIZE_PARAM_MIN:
    138         switch (nParam) {
    139         case PM_PAR_SKY:
    140             params_min = -1000;
    141             break;
    142         case PM_PAR_I0:
    143             params_min =   0.01;
    144             break;
    145         case PM_PAR_XPOS:
    146             params_min =  -100;
    147             break;
    148         case PM_PAR_YPOS:
    149             params_min =  -100;
    150             break;
    151         case PM_PAR_SXX:
    152             params_min =   0.5;
    153             break;
    154         case PM_PAR_SYY:
    155             params_min =   0.5;
    156             break;
    157         case PM_PAR_SXY:
    158             params_min =  -q2;
    159             break;
    160         case PM_PAR_7:
    161             params_min =  -1.0;
    162             break;
    163         default:
    164             psAbort("invalid parameter %d for param min test", nParam);
    165         }
    166         if (params[nParam] < params_min) {
    167             params[nParam] = params_min;
    168             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    169                      nParam, params[nParam], params_min);
    170             return false;
    171         }
    172         return true;
    173     case PS_MINIMIZE_PARAM_MAX:
    174         switch (nParam) {
    175         case PM_PAR_SKY:
    176             params_max =   1e5;
    177             break;
    178         case PM_PAR_I0:
    179             params_max =   1e8;
    180             break;
    181         case PM_PAR_XPOS:
    182             params_max =   1e4;
    183             break;
    184         case PM_PAR_YPOS:
    185             params_max =   1e4;
    186             break;
    187         case PM_PAR_SXX:
    188             params_max =   100;
    189             break;
    190         case PM_PAR_SYY:
    191             params_max =   100;
    192             break;
    193         case PM_PAR_SXY:
    194             params_max =  +q2;
    195             break;
    196         case PM_PAR_7:
    197             params_max =  20.0;
    198             break;
    199         default:
    200             psAbort("invalid parameter %d for param max test", nParam);
    201         }
    202         if (params[nParam] > params_max) {
    203             params[nParam] = params_max;
    204             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    205                      nParam, params[nParam], params_max);
    206             return false;
    207         }
    208         return true;
    209     default:
     131      case PS_MINIMIZE_BETA_LIMIT: {
     132          psAssert(beta, "Require beta to limit beta");
     133          float limit = betaUse[nParam];
     134          if (nParam == PM_PAR_SXY) {
     135              limit *= q2;
     136          }
     137          if (fabs(beta[nParam]) > fabs(limit)) {
     138              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     139              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     140                      nParam, beta[nParam], limit);
     141              return false;
     142          }
     143          return true;
     144      }
     145      case PS_MINIMIZE_PARAM_MIN: {
     146          psAssert(params, "Require parameters to limit parameters");
     147          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     148          float limit = paramsMinUse[nParam];
     149          if (nParam == PM_PAR_SXY) {
     150              limit *= q2;
     151          }
     152          if (params[nParam] < limit) {
     153              params[nParam] = limit;
     154              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     155                      nParam, params[nParam], limit);
     156              return false;
     157          }
     158          return true;
     159      }
     160      case PS_MINIMIZE_PARAM_MAX: {
     161          psAssert(params, "Require parameters to limit parameters");
     162          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     163          float limit = paramsMaxUse[nParam];
     164          if (nParam == PM_PAR_SXY) {
     165              limit *= q2;
     166          }
     167          if (params[nParam] > limit) {
     168              params[nParam] = limit;
     169              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     170                      nParam, params[nParam], limit);
     171              return false;
     172          }
     173          return true;
     174      }
     175      default:
    210176        psAbort("invalid choice for limits");
    211177    }
    … …  
    299265    psF32 *PAR = params->data.F32;
    300266
    301     if (flux <= 0)
    302         return (1.0);
    303     if (PAR[PM_PAR_I0] <= 0)
    304         return (1.0);
    305     if (flux >= PAR[PM_PAR_I0])
    306         return (1.0);
     267    if (flux <= 0) {
     268        return 1.0;
     269    }
     270    if (PAR[PM_PAR_I0] <= 0) {
     271        return 1.0;
     272    }
     273    if (flux >= PAR[PM_PAR_I0]) {
     274        return 1.0;
     275    }
     276    if (PAR[PM_PAR_7] == 0.0) {
     277        return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
     278    }
    307279
    308280    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    468440}
    469441
     442
     443void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     444{
     445    switch (type) {
     446      case PM_MODEL_LIMITS_NONE:
     447        paramsMinUse = NULL;
     448        paramsMaxUse = NULL;
     449        limitsApply = true;
     450        break;
     451      case PM_MODEL_LIMITS_IGNORE:
     452        paramsMinUse = NULL;
     453        paramsMaxUse = NULL;
     454        limitsApply = false;
     455        break;
     456      case PM_MODEL_LIMITS_LAX:
     457        paramsMinUse = paramsMinLax;
     458        paramsMaxUse = paramsMaxLax;
     459        limitsApply = true;
     460        break;
     461      case PM_MODEL_LIMITS_STRICT:
     462        paramsMinUse = paramsMinStrict;
     463        paramsMaxUse = paramsMaxStrict;
     464        limitsApply = true;
     465        break;
     466      default:
     467        psAbort("Unrecognised model limits type: %x", type);
     468    }
     469    return;
     470}
     471
     472
    470473# undef PM_MODEL_FUNC
    471474# undef PM_MODEL_FLUX
    … …  
    476479# undef PM_MODEL_PARAMS_FROM_PSF
    477480# undef PM_MODEL_FIT_STATUS
     481# undef PM_MODEL_SET_LIMITS
    478482# undef ALPHA
    479483# undef ALPHA_M

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

    r20001 r25738  
    2020   *****************************************************************************/
    2121
     22#include <stdio.h>
     23#include <pslib.h>
     24
     25#include "pmMoments.h"
     26#include "pmPeaks.h"
     27#include "pmSource.h"
     28#include "pmModel.h"
     29#include "pmModel_QGAUSS.h"
     30
    2231# define PM_MODEL_FUNC            pmModelFunc_QGAUSS
    2332# define PM_MODEL_FLUX            pmModelFlux_QGAUSS
    … …  
    2837# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS
    2938# define PM_MODEL_FIT_STATUS      pmModelFitStatus_QGAUSS
     39# define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
     40
     41// Lax parameter limits
     42static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
     43static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     44
     45// Strict parameter limits
     46static float *paramsMinStrict = paramsMinLax;
     47static float *paramsMaxStrict = paramsMaxLax;
     48
     49// Parameter limits to use
     50static float *paramsMinUse = paramsMinLax;
     51static float *paramsMaxUse = paramsMaxLax;
     52static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     53
     54static bool limitsApply = true;         // Apply limits?
    3055
    3156psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    79104# define AR_MAX 20.0
    80105# define AR_RATIO 0.99
     106
    81107bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    82108{
    83     float beta_lim = 0, params_min = 0, params_max = 0;
    84     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     109    if (!limitsApply) {
     110        return true;
     111    }
     112    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    85113
    86114    // we need to calculate the limits for SXY specially
     115    float q2 = NAN;
    87116    if (nParam == PM_PAR_SXY) {
    88         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    89         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    90         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     117        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     118        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     119        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    91120        q1 = (q1 < 0.0) ? 0.0 : q1;
    92121        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    93122        // angle and let f2,f1 fight it out
    94         q2  = 0.5*sqrt (q1);
     123        q2 = 0.5*sqrtf(q1);
    95124    }
    96125
    97126    switch (mode) {
    98     case PS_MINIMIZE_BETA_LIMIT:
    99         switch (nParam) {
    100         case PM_PAR_SKY:
    101             beta_lim = 1000;
    102             break;
    103         case PM_PAR_I0:
    104             beta_lim = 3e6;
    105             break;
    106         case PM_PAR_XPOS:
    107             beta_lim = 5;
    108             break;
    109         case PM_PAR_YPOS:
    110             beta_lim = 5;
    111             break;
    112         case PM_PAR_SXX:
    113             beta_lim = 1.0;
    114             break;
    115         case PM_PAR_SYY:
    116             beta_lim = 1.0;
    117             break;
    118         case PM_PAR_SXY:
    119             beta_lim =  0.5*q2;
    120             break;
    121         case PM_PAR_7:
    122             beta_lim = 2.0;
    123             break;
    124         default:
    125             psAbort("invalid parameter %d for beta test", nParam);
    126         }
    127         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    128             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    129             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    130                      nParam, beta[nParam], beta_lim);
    131             return false;
    132         }
    133         return true;
    134     case PS_MINIMIZE_PARAM_MIN:
    135         switch (nParam) {
    136         case PM_PAR_SKY:
    137             params_min = -1000;
    138             break;
    139         case PM_PAR_I0:
    140             params_min =   0.01;
    141             break;
    142         case PM_PAR_XPOS:
    143             params_min =  -100;
    144             break;
    145         case PM_PAR_YPOS:
    146             params_min =  -100;
    147             break;
    148         case PM_PAR_SXX:
    149             params_min =   0.5;
    150             break;
    151         case PM_PAR_SYY:
    152             params_min =   0.5;
    153             break;
    154         case PM_PAR_SXY:
    155             params_min =  -q2;
    156             break;
    157         case PM_PAR_7:
    158             params_min =   0.1;
    159             break;
    160         default:
    161             psAbort("invalid parameter %d for param min test", nParam);
    162         }
    163         if (params[nParam] < params_min) {
    164             params[nParam] = params_min;
    165             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    166                      nParam, params[nParam], params_min);
    167             return false;
    168         }
    169         return true;
    170     case PS_MINIMIZE_PARAM_MAX:
    171         switch (nParam) {
    172         case PM_PAR_SKY:
    173             params_max =   1e5;
    174             break;
    175         case PM_PAR_I0:
    176             params_max =   1e8;
    177             break;
    178         case PM_PAR_XPOS:
    179             params_max =   1e4;
    180             break;
    181         case PM_PAR_YPOS:
    182             params_max =   1e4;
    183             break;
    184         case PM_PAR_SXX:
    185             params_max =   100;
    186             break;
    187         case PM_PAR_SYY:
    188             params_max =   100;
    189             break;
    190         case PM_PAR_SXY:
    191             params_max =  +q2;
    192             break;
    193         case PM_PAR_7:
    194             params_max =  20.0;
    195             break;
    196         default:
    197             psAbort("invalid parameter %d for param max test", nParam);
    198         }
    199         if (params[nParam] > params_max) {
    200             params[nParam] = params_max;
    201             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    202                      nParam, params[nParam], params_max);
    203             return false;
    204         }
    205         return true;
     127      case PS_MINIMIZE_BETA_LIMIT: {
     128          psAssert(beta, "Require beta to limit beta");
     129          float limit = betaUse[nParam];
     130          if (nParam == PM_PAR_SXY) {
     131              limit *= q2;
     132          }
     133          if (fabs(beta[nParam]) > fabs(limit)) {
     134              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     135              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     136                      nParam, beta[nParam], limit);
     137              return false;
     138          }
     139          return true;
     140      }
     141      case PS_MINIMIZE_PARAM_MIN: {
     142          psAssert(params, "Require parameters to limit parameters");
     143          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     144          float limit = paramsMinUse[nParam];
     145          if (nParam == PM_PAR_SXY) {
     146              limit *= q2;
     147          }
     148          if (params[nParam] < limit) {
     149              params[nParam] = limit;
     150              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     151                      nParam, params[nParam], limit);
     152              return false;
     153          }
     154          return true;
     155      }
     156      case PS_MINIMIZE_PARAM_MAX: {
     157          psAssert(params, "Require parameters to limit parameters");
     158          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     159          float limit = paramsMaxUse[nParam];
     160          if (nParam == PM_PAR_SXY) {
     161              limit *= q2;
     162          }
     163          if (params[nParam] > limit) {
     164              params[nParam] = limit;
     165              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     166                      nParam, params[nParam], limit);
     167              return false;
     168          }
     169          return true;
     170      }
    206171    default:
    207172        psAbort("invalid choice for limits");
    … …  
    464429}
    465430
     431
     432void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     433{
     434    switch (type) {
     435      case PM_MODEL_LIMITS_NONE:
     436        paramsMinUse = NULL;
     437        paramsMaxUse = NULL;
     438        limitsApply = true;
     439        break;
     440      case PM_MODEL_LIMITS_IGNORE:
     441        paramsMinUse = NULL;
     442        paramsMaxUse = NULL;
     443        limitsApply = false;
     444        break;
     445      case PM_MODEL_LIMITS_LAX:
     446        paramsMinUse = paramsMinLax;
     447        paramsMaxUse = paramsMaxLax;
     448        limitsApply = true;
     449        break;
     450      case PM_MODEL_LIMITS_STRICT:
     451        paramsMinUse = paramsMinStrict;
     452        paramsMaxUse = paramsMaxStrict;
     453        limitsApply = true;
     454        break;
     455      default:
     456        psAbort("Unrecognised model limits type: %x", type);
     457    }
     458    return;
     459}
     460
    466461# undef PM_MODEL_FUNC
    467462# undef PM_MODEL_FLUX
    … …  
    472467# undef PM_MODEL_PARAMS_FROM_PSF
    473468# undef PM_MODEL_FIT_STATUS
     469# undef PM_MODEL_SET_LIMITS

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

    r20001 r25738  
    2020 *****************************************************************************/
    2121
     22#include <stdio.h>
     23#include <pslib.h>
     24
     25#include "pmMoments.h"
     26#include "pmPeaks.h"
     27#include "pmSource.h"
     28#include "pmModel.h"
     29#include "pmModel_RGAUSS.h"
     30
    2231# define PM_MODEL_FUNC            pmModelFunc_RGAUSS
    2332# define PM_MODEL_FLUX            pmModelFlux_RGAUSS
    … …  
    2837# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS
    2938# define PM_MODEL_FIT_STATUS      pmModelFitStatus_RGAUSS
     39# define PM_MODEL_SET_LIMITS      pmModelSetLimits_RGAUSS
     40
     41// Lax parameter limits
     42static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
     43static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     44
     45// Strict parameter limits
     46static float *paramsMinStrict = paramsMinLax;
     47static float *paramsMaxStrict = paramsMaxLax;
     48
     49// Parameter limits to use
     50static float *paramsMinUse = paramsMinLax;
     51static float *paramsMaxUse = paramsMaxLax;
     52static float betaUse[] = { 1000, 3e6, 5, 5, 0.5, 0.5, 0.5, 0.5 };
     53
     54static bool limitsApply = true;         // Apply limits?
    3055
    3156psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    7398# define AR_MAX 20.0
    7499# define AR_RATIO 0.99
     100
    75101bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    76102{
    77     float beta_lim = 0, params_min = 0, params_max = 0;
    78     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     103    if (!limitsApply) {
     104        return true;
     105    }
     106    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    79107
    80108    // we need to calculate the limits for SXY specially
     109    float q2 = NAN;
    81110    if (nParam == PM_PAR_SXY) {
    82         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    83         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    84         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     111        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     112        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     113        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    85114        q1 = (q1 < 0.0) ? 0.0 : q1;
    86115        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    87116        // angle and let f2,f1 fight it out
    88         q2  = 0.5*sqrt (q1);
     117        q2 = 0.5*sqrtf(q1);
    89118    }
    90119
    91120    switch (mode) {
    92     case PS_MINIMIZE_BETA_LIMIT:
    93         switch (nParam) {
    94         case PM_PAR_SKY:
    95             beta_lim = 1000;
    96             break;
    97         case PM_PAR_I0:
    98             beta_lim = 3e6;
    99             break;
    100         case PM_PAR_XPOS:
    101             beta_lim = 5;
    102             break;
    103         case PM_PAR_YPOS:
    104             beta_lim = 5;
    105             break;
    106         case PM_PAR_SXX:
    107             beta_lim = 0.5;
    108             break;
    109         case PM_PAR_SYY:
    110             beta_lim = 0.5;
    111             break;
    112         case PM_PAR_SXY:
    113             beta_lim =  0.5*q2;
    114             break;
    115         case PM_PAR_7:
    116             beta_lim = 0.5;
    117             break;
    118         default:
    119             psAbort("invalid parameter %d for beta test", nParam);
    120         }
    121         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    122             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    123             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    124                      nParam, beta[nParam], beta_lim);
    125             return false;
    126         }
    127         return true;
    128     case PS_MINIMIZE_PARAM_MIN:
    129         switch (nParam) {
    130         case PM_PAR_SKY:
    131             params_min = -1000;
    132             break;
    133         case PM_PAR_I0:
    134             params_min =   0.01;
    135             break;
    136         case PM_PAR_XPOS:
    137             params_min =  -100;
    138             break;
    139         case PM_PAR_YPOS:
    140             params_min =  -100;
    141             break;
    142         case PM_PAR_SXX:
    143             params_min =   0.5;
    144             break;
    145         case PM_PAR_SYY:
    146             params_min =   0.5;
    147             break;
    148         case PM_PAR_SXY:
    149             params_min =  -q2;
    150             break;
    151         case PM_PAR_7:
    152             params_min =   1.25;
    153             break;
    154         default:
    155             psAbort("invalid parameter %d for param min test", nParam);
    156         }
    157         if (params[nParam] < params_min) {
    158             params[nParam] = params_min;
    159             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    160                      nParam, params[nParam], params_min);
    161             return false;
    162         }
    163         return true;
    164     case PS_MINIMIZE_PARAM_MAX:
    165         switch (nParam) {
    166         case PM_PAR_SKY:
    167             params_max =   1e5;
    168             break;
    169         case PM_PAR_I0:
    170             params_max =   1e8;
    171             break;
    172         case PM_PAR_XPOS:
    173             params_max =   1e4;
    174             break;
    175         case PM_PAR_YPOS:
    176             params_max =   1e4;
    177             break;
    178         case PM_PAR_SXX:
    179             params_max =   100;
    180             break;
    181         case PM_PAR_SYY:
    182             params_max =   100;
    183             break;
    184         case PM_PAR_SXY:
    185             params_max =  +q2;
    186             break;
    187         case PM_PAR_7:
    188             params_max =  4.0;
    189             break;
    190         default:
    191             psAbort("invalid parameter %d for param max test", nParam);
    192         }
    193         if (params[nParam] > params_max) {
    194             params[nParam] = params_max;
    195             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    196                      nParam, params[nParam], params_max);
    197             return false;
    198         }
    199         return true;
    200     default:
     121      case PS_MINIMIZE_BETA_LIMIT: {
     122          psAssert(beta, "Require beta to limit beta");
     123          float limit = betaUse[nParam];
     124          if (nParam == PM_PAR_SXY) {
     125              limit *= q2;
     126          }
     127          if (fabs(beta[nParam]) > fabs(limit)) {
     128              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     129              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     130                      nParam, beta[nParam], limit);
     131              return false;
     132          }
     133          return true;
     134      }
     135      case PS_MINIMIZE_PARAM_MIN: {
     136          psAssert(params, "Require parameters to limit parameters");
     137          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     138          float limit = paramsMinUse[nParam];
     139          if (nParam == PM_PAR_SXY) {
     140              limit *= q2;
     141          }
     142          if (params[nParam] < limit) {
     143              params[nParam] = limit;
     144              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     145                      nParam, params[nParam], limit);
     146              return false;
     147          }
     148          return true;
     149      }
     150      case PS_MINIMIZE_PARAM_MAX: {
     151          psAssert(params, "Require parameters to limit parameters");
     152          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     153          float limit = paramsMaxUse[nParam];
     154          if (nParam == PM_PAR_SXY) {
     155              limit *= q2;
     156          }
     157          if (params[nParam] > limit) {
     158              params[nParam] = limit;
     159              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     160                      nParam, params[nParam], limit);
     161              return false;
     162          }
     163          return true;
     164      }
     165      default:
    201166        psAbort("invalid choice for limits");
    202167    }
    … …  
    204169    return false;
    205170}
     171
    206172
    207173// make an initial guess for parameters
    … …  
    456422}
    457423
     424
     425void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     426{
     427    switch (type) {
     428      case PM_MODEL_LIMITS_NONE:
     429        paramsMinUse = NULL;
     430        paramsMaxUse = NULL;
     431        limitsApply = true;
     432        break;
     433      case PM_MODEL_LIMITS_IGNORE:
     434        paramsMinUse = NULL;
     435        paramsMaxUse = NULL;
     436        limitsApply = false;
     437        break;
     438      case PM_MODEL_LIMITS_LAX:
     439        paramsMinUse = paramsMinLax;
     440        paramsMaxUse = paramsMaxLax;
     441        limitsApply = true;
     442        break;
     443      case PM_MODEL_LIMITS_STRICT:
     444        paramsMinUse = paramsMinStrict;
     445        paramsMaxUse = paramsMaxStrict;
     446        limitsApply = true;
     447        break;
     448      default:
     449        psAbort("Unrecognised model limits type: %x", type);
     450    }
     451    return;
     452}
     453
    458454# undef PM_MODEL_FUNC
    459455# undef PM_MODEL_FLUX
    … …  
    464460# undef PM_MODEL_PARAMS_FROM_PSF
    465461# undef PM_MODEL_FIT_STATUS
     462# undef PM_MODEL_SET_LIMITS

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

    r20001 r25738  
    2323   *****************************************************************************/
    2424
     25#include <stdio.h>
     26#include <pslib.h>
     27
     28#include "pmMoments.h"
     29#include "pmPeaks.h"
     30#include "pmSource.h"
     31#include "pmModel.h"
     32#include "pmModel_SERSIC.h"
     33
    2534# define PM_MODEL_FUNC            pmModelFunc_SERSIC
    2635# define PM_MODEL_FLUX            pmModelFlux_SERSIC
    … …  
    3140# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC
    3241# define PM_MODEL_FIT_STATUS      pmModelFitStatus_SERSIC
     42# define PM_MODEL_SET_LIMITS      pmModelSetLimits_SERSIC
     43
     44// Lax parameter limits
     45static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
     46static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     47
     48// Strict parameter limits
     49static float *paramsMinStrict = paramsMinLax;
     50static float *paramsMaxStrict = paramsMaxLax;
     51
     52// Parameter limits to use
     53static float *paramsMinUse = paramsMinLax;
     54static float *paramsMaxUse = paramsMaxLax;
     55static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     56
     57static bool limitsApply = true;         // Apply limits?
    3358
    3459psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    91116bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    92117{
    93     float beta_lim = 0, params_min = 0, params_max = 0;
    94     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     118    if (!limitsApply) {
     119        return true;
     120    }
     121    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    95122
    96123    // we need to calculate the limits for SXY specially
     124    float q2 = NAN;
    97125    if (nParam == PM_PAR_SXY) {
    98         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    99         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    100         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     126        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     127        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     128        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    101129        q1 = (q1 < 0.0) ? 0.0 : q1;
    102130        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    103131        // angle and let f2,f1 fight it out
    104         q2  = 0.5*sqrt (q1);
     132        q2 = 0.5*sqrtf(q1);
    105133    }
    106134
    107135    switch (mode) {
    108     case PS_MINIMIZE_BETA_LIMIT:
    109         switch (nParam) {
    110         case PM_PAR_SKY:
    111             beta_lim = 1000;
    112             break;
    113         case PM_PAR_I0:
    114             beta_lim = 3e6;
    115             break;
    116         case PM_PAR_XPOS:
    117             beta_lim = 5;
    118             break;
    119         case PM_PAR_YPOS:
    120             beta_lim = 5;
    121             break;
    122         case PM_PAR_SXX:
    123             beta_lim = 1.0;
    124             break;
    125         case PM_PAR_SYY:
    126             beta_lim = 1.0;
    127             break;
    128         case PM_PAR_SXY:
    129             beta_lim =  0.5*q2;
    130             break;
    131         case PM_PAR_7:
    132             beta_lim = 2.0;
    133             break;
    134         default:
    135             psAbort("invalid parameter %d for beta test", nParam);
    136         }
    137         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    138             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    139             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    140                      nParam, beta[nParam], beta_lim);
    141             return false;
    142         }
    143         return true;
    144     case PS_MINIMIZE_PARAM_MIN:
    145         switch (nParam) {
    146         case PM_PAR_SKY:
    147             params_min = -1000;
    148             break;
    149         case PM_PAR_I0:
    150             params_min =     0.01;
    151             break;
    152         case PM_PAR_XPOS:
    153             params_min =  -100;
    154             break;
    155         case PM_PAR_YPOS:
    156             params_min =  -100;
    157             break;
    158         case PM_PAR_SXX:
    159             params_min =   0.05;
    160             break;
    161         case PM_PAR_SYY:
    162             params_min =   0.05;
    163             break;
    164         case PM_PAR_SXY:
    165             params_min =  -q2;
    166             break;
    167         case PM_PAR_7:
    168             params_min =   0.05;
    169             break;
    170         default:
    171             psAbort("invalid parameter %d for param min test", nParam);
    172         }
    173         if (params[nParam] < params_min) {
    174             params[nParam] = params_min;
    175             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    176                      nParam, params[nParam], params_min);
    177             return false;
    178         }
    179         return true;
    180     case PS_MINIMIZE_PARAM_MAX:
    181         switch (nParam) {
    182         case PM_PAR_SKY:
    183             params_max =   1e5;
    184             break;
    185         case PM_PAR_I0:
    186             params_max =   1e8;
    187             break;
    188         case PM_PAR_XPOS:
    189             params_max =   1e4;
    190             break;
    191         case PM_PAR_YPOS:
    192             params_max =   1e4;
    193             break;
    194         case PM_PAR_SXX:
    195             params_max =   100;
    196             break;
    197         case PM_PAR_SYY:
    198             params_max =   100;
    199             break;
    200         case PM_PAR_SXY:
    201             params_max =  +q2;
    202             break;
    203         case PM_PAR_7:
    204             params_max =   4.0;
    205             break;
    206         default:
    207             psAbort("invalid parameter %d for param max test", nParam);
    208         }
    209         if (params[nParam] > params_max) {
    210             params[nParam] = params_max;
    211             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    212                      nParam, params[nParam], params_max);
    213             return false;
    214         }
    215         return true;
    216     default:
     136      case PS_MINIMIZE_BETA_LIMIT: {
     137          psAssert(beta, "Require beta to limit beta");
     138          float limit = betaUse[nParam];
     139          if (nParam == PM_PAR_SXY) {
     140              limit *= q2;
     141          }
     142          if (fabs(beta[nParam]) > fabs(limit)) {
     143              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     144              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     145                      nParam, beta[nParam], limit);
     146              return false;
     147          }
     148          return true;
     149      }
     150      case PS_MINIMIZE_PARAM_MIN: {
     151          psAssert(params, "Require parameters to limit parameters");
     152          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     153          float limit = paramsMinUse[nParam];
     154          if (nParam == PM_PAR_SXY) {
     155              limit *= q2;
     156          }
     157          if (params[nParam] < limit) {
     158              params[nParam] = limit;
     159              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     160                      nParam, params[nParam], limit);
     161              return false;
     162          }
     163          return true;
     164      }
     165      case PS_MINIMIZE_PARAM_MAX: {
     166          psAssert(params, "Require parameters to limit parameters");
     167          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     168          float limit = paramsMaxUse[nParam];
     169          if (nParam == PM_PAR_SXY) {
     170              limit *= q2;
     171          }
     172          if (params[nParam] > limit) {
     173              params[nParam] = limit;
     174              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     175                      nParam, params[nParam], limit);
     176              return false;
     177          }
     178          return true;
     179      }
     180      default:
    217181        psAbort("invalid choice for limits");
    218182    }
    … …  
    220184    return false;
    221185}
    222 
    223186
    224187// make an initial guess for parameters
    … …  
    447410
    448411    fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n",
    449              dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
     412             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    450413
    451414    return status;
     415}
     416
     417
     418void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     419{
     420    switch (type) {
     421      case PM_MODEL_LIMITS_NONE:
     422        paramsMinUse = NULL;
     423        paramsMaxUse = NULL;
     424        limitsApply = true;
     425        break;
     426      case PM_MODEL_LIMITS_IGNORE:
     427        paramsMinUse = NULL;
     428        paramsMaxUse = NULL;
     429        limitsApply = false;
     430        break;
     431      case PM_MODEL_LIMITS_LAX:
     432        paramsMinUse = paramsMinLax;
     433        paramsMaxUse = paramsMaxLax;
     434        limitsApply = true;
     435        break;
     436      case PM_MODEL_LIMITS_STRICT:
     437        paramsMinUse = paramsMinStrict;
     438        paramsMaxUse = paramsMaxStrict;
     439        limitsApply = true;
     440        break;
     441      default:
     442        psAbort("Unrecognised model limits type: %x", type);
     443    }
     444    return;
    452445}
    453446
    … …  
    460453# undef PM_MODEL_PARAMS_FROM_PSF
    461454# undef PM_MODEL_FIT_STATUS
     455# undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/pmModel.c

    r23187 r25738  
    8686    tmp->modelParamsFromPSF = class->modelParamsFromPSF;
    8787    tmp->modelFitStatus     = class->modelFitStatus;
     88    tmp->modelSetLimits     = class->modelSetLimits;
    8889
    8990    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
    … …  
    235236
    236237    if (model->residuals) {
    237         DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
    238         DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
    239         Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
    240         Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
    241         Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
    242         Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
    243         if (Ro) {
    244             NX = model->residuals->Ro->numCols;
    245             NY = model->residuals->Ro->numRows;
    246         }           
     238        DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
     239        DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
     240        Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
     241        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
     242        Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
     243        Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
     244        if (Ro) {
     245            NX = model->residuals->Ro->numCols;
     246            NY = model->residuals->Ro->numRows;
     247        }
    247248    }
    248249
    … …  
    256257
    257258            // XXX should we use using 0.5 pixel offset?
    258             // Convert to coordinate in parent image, with offset (dx,dy)
     259            // Convert to coordinate in parent image, with offset (dx,dy)
    259260            imageCol = ix + image->col0 - dx;
    260261            imageRow = iy + image->row0 - dy;
    … …  
    276277                float rx = xBin*ix + DX;
    277278
    278                 int rx0 = rx - 0.5;
    279                 int rx1 = rx + 0.5;
    280                 int ry0 = ry - 0.5;
    281                 int ry1 = ry + 0.5;
    282 
    283                 if (rx0 < 0) goto skip;
    284                 if (ry0 < 0) goto skip;
    285                 if (rx1 >= NX) goto skip;
    286                 if (ry1 >= NY) goto skip;
    287 
    288                 // these go from 0.0 to 1.0 between the centers of the pixels
    289                 float fx = rx - 0.5 - rx0;
    290                 float Fx = 1.0 - fx;
    291                 float fy = ry - 0.5 - ry0;
    292                 float Fy = 1.0 - fy;
    293 
    294                 // check the residual image mask (if set). give up if any of the 4 pixels are masked.
    295                 if (Rm) {
    296                     if (Rm[ry0][rx0]) goto skip;
    297                     if (Rm[ry0][rx1]) goto skip;
    298                     if (Rm[ry1][rx0]) goto skip;
    299                     if (Rm[ry1][rx1]) goto skip;
    300                 }
    301 
    302                 // a possible further optimization if we re-use these values
    303                 // XXX allow for masked pixels, and add pixel weights
    304                 float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
    305                 float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
    306                 float Vo = V0*Fy + V1*fy;
    307                 if (!isfinite(Vo)) goto skip;
    308 
    309                 float Vx = 0.0;
    310                 float Vy = 0.0;
    311 
    312                 // skip Rx,Ry if Ro is masked
    313                 if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
    314                     V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
    315                     V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
    316                     Vx = V0*Fy + V1*fy;
    317 
    318                     V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
    319                     V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
    320                     Vy = V0*Fy + V1*fy;
    321                 }
    322                 if (!isfinite(Vx)) goto skip;
    323                 if (!isfinite(Vy)) goto skip;
    324 
    325                 // 2D residual variations are set for the true source position
    326                 pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
     279                int rx0 = rx - 0.5;
     280                int rx1 = rx + 0.5;
     281                int ry0 = ry - 0.5;
     282                int ry1 = ry + 0.5;
     283
     284                if (rx0 < 0) goto skip;
     285                if (ry0 < 0) goto skip;
     286                if (rx1 >= NX) goto skip;
     287                if (ry1 >= NY) goto skip;
     288
     289                // these go from 0.0 to 1.0 between the centers of the pixels
     290                float fx = rx - 0.5 - rx0;
     291                float Fx = 1.0 - fx;
     292                float fy = ry - 0.5 - ry0;
     293                float Fy = 1.0 - fy;
     294
     295                // check the residual image mask (if set). give up if any of the 4 pixels are masked.
     296                if (Rm) {
     297                    if (Rm[ry0][rx0]) goto skip;
     298                    if (Rm[ry0][rx1]) goto skip;
     299                    if (Rm[ry1][rx0]) goto skip;
     300                    if (Rm[ry1][rx1]) goto skip;
     301                }
     302
     303                // a possible further optimization if we re-use these values
     304                // XXX allow for masked pixels, and add pixel weights
     305                float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
     306                float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
     307                float Vo = V0*Fy + V1*fy;
     308                if (!isfinite(Vo)) goto skip;
     309
     310                float Vx = 0.0;
     311                float Vy = 0.0;
     312
     313                // skip Rx,Ry if Ro is masked
     314                if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
     315                    V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
     316                    V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
     317                    Vx = V0*Fy + V1*fy;
     318
     319                    V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
     320                    V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
     321                    Vy = V0*Fy + V1*fy;
     322                }
     323                if (!isfinite(Vx)) goto skip;
     324                if (!isfinite(Vy)) goto skip;
     325
     326                // 2D residual variations are set for the true source position
     327                pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
    327328            }
    328329
    329         skip:
     330        skip:
    330331            // add or subtract the value
    331332            if (add) {

  • trunk/psModules/src/objects/pmModel.h

    r21516 r25738  
    4444} pmModelOpMode;
    4545
     46/// Parameter limit types
     47typedef enum {
     48    PM_MODEL_LIMITS_NONE,               ///< Apply no limits: suitable for debugging
     49    PM_MODEL_LIMITS_IGNORE,             ///< Ignore all limits: fit can go to town
     50    PM_MODEL_LIMITS_LAX,                ///< Lax limits: attempting to reproduce mildly bad data
     51    PM_MODEL_LIMITS_STRICT,             ///< Strict limits: good quality data
     52} pmModelLimitsType;
     53
    4654typedef struct pmModel pmModel;
    4755typedef struct pmSource pmSource;
    … …  
    7482//  This function returns the success / failure status of the given model fit
    7583typedef bool (*pmModelFitStatusFunc)(pmModel *model);
     84
     85//  This function sets the parameter limits for the given model
     86typedef bool (*pmModelSetLimitsFunc)(pmModelLimitsType type);
    7687
    7788/** pmModel data structure
    … …  
    108119    pmModelParamsFromPSF modelParamsFromPSF;
    109120    pmModelFitStatusFunc modelFitStatus;
     121    pmModelSetLimitsFunc modelSetLimits;
    110122};
    111123
    … …  
    151163    pmModel *model,                     ///< The input pmModel
    152164    pmModelOpMode mode,                 ///< mode to control how the model is added into the image
    153     psImageMaskType maskVal             ///< Value to mask
     165    psImageMaskType maskVal             ///< Value to mask
    154166);
    155167
    … …  
    169181    pmModel *model,                     ///< The input pmModel
    170182    pmModelOpMode mode,                 ///< mode to control how the model is added into the image
    171     psImageMaskType maskVal             ///< Value to mask
     183    psImageMaskType maskVal             ///< Value to mask
    172184);
    173185
    … …  
    202214);
    203215
     216
     217/// Set the model parameter limits for the given model
     218///
     219/// Wraps the model-specific pmModelSetLimitsFunc function.
     220bool pmModelSetLimits(
     221    const pmModel *model,               ///< Model of interest
     222    pmModelLimits type                  ///< Type of limits
     223    );
     224
     225
    204226/// @}
    205227# endif /* PM_MODEL_H */

  • trunk/psModules/src/objects/pmModelClass.c

    r20937 r25738  
    4040double sqrt (double x);
    4141
    42 # include "models/pmModel_GAUSS.c"
    43 # include "models/pmModel_PGAUSS.c"
    44 # include "models/pmModel_QGAUSS.c"
    45 # include "models/pmModel_PS1_V1.c"
    46 # include "models/pmModel_RGAUSS.c"
    47 # include "models/pmModel_SERSIC.c"
     42# include "models/pmModel_GAUSS.h"
     43# include "models/pmModel_PGAUSS.h"
     44# include "models/pmModel_QGAUSS.h"
     45# include "models/pmModel_PS1_V1.h"
     46# include "models/pmModel_RGAUSS.h"
     47# include "models/pmModel_SERSIC.h"
    4848
    4949static pmModelClass defaultModels[] = {
    50     {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelParamsFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
    51     {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelParamsFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
    52     {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelParamsFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
    53     {"PS_MODEL_PS1_V1",       8, pmModelFunc_PS1_V1,  pmModelFlux_PS1_V1,  pmModelRadius_PS1_V1,  pmModelLimits_PS1_V1,  pmModelGuess_PS1_V1, pmModelFromPSF_PS1_V1, pmModelParamsFromPSF_PS1_V1, pmModelFitStatus_PS1_V1},
    54     {"PS_MODEL_RGAUSS",       8, pmModelFunc_RGAUSS,  pmModelFlux_RGAUSS,  pmModelRadius_RGAUSS,  pmModelLimits_RGAUSS,  pmModelGuess_RGAUSS, pmModelFromPSF_RGAUSS, pmModelParamsFromPSF_RGAUSS, pmModelFitStatus_RGAUSS},
    55     {"PS_MODEL_SERSIC",       8, pmModelFunc_SERSIC,  pmModelFlux_SERSIC,  pmModelRadius_SERSIC,  pmModelLimits_SERSIC,  pmModelGuess_SERSIC, pmModelFromPSF_SERSIC, pmModelParamsFromPSF_SERSIC, pmModelFitStatus_SERSIC}
     50    {"PS_MODEL_GAUSS",        7, (pmModelFunc)pmModelFunc_GAUSS,   (pmModelFlux)pmModelFlux_GAUSS,   (pmModelRadius)pmModelRadius_GAUSS,   (pmModelLimits)pmModelLimits_GAUSS,   (pmModelGuessFunc)pmModelGuess_GAUSS,  (pmModelFromPSFFunc)pmModelFromPSF_GAUSS,  (pmModelParamsFromPSF)pmModelParamsFromPSF_GAUSS,  (pmModelFitStatusFunc)pmModelFitStatus_GAUSS,  (pmModelSetLimitsFunc)pmModelSetLimits_GAUSS  },
     51    {"PS_MODEL_PGAUSS",       7, (pmModelFunc)pmModelFunc_PGAUSS,  (pmModelFlux)pmModelFlux_PGAUSS,  (pmModelRadius)pmModelRadius_PGAUSS,  (pmModelLimits)pmModelLimits_PGAUSS,  (pmModelGuessFunc)pmModelGuess_PGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_PGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_PGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_PGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_PGAUSS },
     52    {"PS_MODEL_QGAUSS",       8, (pmModelFunc)pmModelFunc_QGAUSS,  (pmModelFlux)pmModelFlux_QGAUSS,  (pmModelRadius)pmModelRadius_QGAUSS,  (pmModelLimits)pmModelLimits_QGAUSS,  (pmModelGuessFunc)pmModelGuess_QGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_QGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_QGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_QGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_QGAUSS },
     53    {"PS_MODEL_PS1_V1",       8, (pmModelFunc)pmModelFunc_PS1_V1,  (pmModelFlux)pmModelFlux_PS1_V1,  (pmModelRadius)pmModelRadius_PS1_V1,  (pmModelLimits)pmModelLimits_PS1_V1,  (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 },
     54    {"PS_MODEL_RGAUSS",       8, (pmModelFunc)pmModelFunc_RGAUSS,  (pmModelFlux)pmModelFlux_RGAUSS,  (pmModelRadius)pmModelRadius_RGAUSS,  (pmModelLimits)pmModelLimits_RGAUSS,  (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS },
     55    {"PS_MODEL_SERSIC",       8, (pmModelFunc)pmModelFunc_SERSIC,  (pmModelFlux)pmModelFlux_SERSIC,  (pmModelRadius)pmModelRadius_SERSIC,  (pmModelLimits)pmModelLimits_SERSIC,  (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }
    5656};
    5757
    … …  
    168168    return (models[type].name);
    169169}
     170
     171
     172void pmModelClassSetLimits(pmModelLimitsType type)
     173{
     174    if (!models) {
     175        pmModelClassInit();
     176    }
     177
     178    for (int i = 0; i < Nmodels; i++) {
     179        if (models[i].modelSetLimits) {
     180            models[i].modelSetLimits(type);
     181        }
     182    }
     183
     184}
     185

  • trunk/psModules/src/objects/pmModelClass.h

    r15697 r25738  
    2929# define PM_MODEL_CLASS_H
    3030
     31#include <pmModel.h>
     32
    3133/// @addtogroup Objects Object Detection / Analysis Functions
    3234/// @{
    3335
    34 typedef struct
    35 {
     36typedef struct {
    3637    char *name;
    3738    int nParams;
    … …  
    4445    pmModelParamsFromPSF modelParamsFromPSF;
    4546    pmModelFitStatusFunc modelFitStatus;
     47    pmModelSetLimitsFunc modelSetLimits;
    4648} pmModelClass;
    4749
    … …  
    7375pmModelType pmModelClassGetType (const char *name);
    7476
     77/// Set parameter limits for all models
     78void pmModelClassSetLimits(pmModelLimitsType type);
     79
     80
    7581/// @}
    7682# endif /* PM_MODEL_CLASS_H */

  • trunk/psphot/doc/efficiency.txt

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/psphot/src/psphotEfficiency.c

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/psphot/src/psphotFake.c

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/psphot/src/psphotModelGroupInit.c

    r14655 r25738  
    11# include "psphotInternal.h"
    22
    3 // Add locally-defined models here.  As these mature, they can be moved to 
     3// Add locally-defined models here.  As these mature, they can be moved to
    44// psModule/src/objects/models
    55
    … …  
    88
    99static pmModelClass userModels[] = {
    10     {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1,  pmModelFlux_TEST1,  pmModelRadius_TEST1,  pmModelLimits_TEST1,  pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1},
    11     {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL,  pmModelFlux_STRAIL,  pmModelRadius_STRAIL,  pmModelLimits_STRAIL,  pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL},
     10    {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1,  pmModelFlux_TEST1,  pmModelRadius_TEST1,  pmModelLimits_TEST1,  pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1, NULL},
     11    {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL,  pmModelFlux_STRAIL,  pmModelRadius_STRAIL,  pmModelLimits_STRAIL,  pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL, NULL},
    1212};
    1313
    14 void psphotModelClassInit (void) 
    15 { 
     14void psphotModelClassInit (void)
     15{
    1616
    1717    // if pmModelClassInit returns false, we have already init'ed
    … …  
    2020    int Nmodels = sizeof (userModels) / sizeof (pmModelClass);
    2121    for (int i = 0; i < Nmodels; i++) {
    22         pmModelClassAdd (&userModels[i]);
     22        pmModelClassAdd (&userModels[i]);
    2323    }
    2424    return;

  • trunk/psphot/src/psphotReadout.c

    r25383 r25738  
    1717    // jointly by the multiple threads, not the total time used by all threads.
    1818    psTimerStart ("psphotReadout");
     19
     20    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
    1921
    2022    // select the current recipe

  • trunk/psphot/src/psphotReadoutMinimal.c

    r25383 r25738  
    1414    // jointly by the multiple threads, not the total time used by all threads.
    1515    psTimerStart ("psphotReadout");
     16
     17    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
    1618
    1719    // select the current recipe

  • trunk/psphot/src/psphotVersion.c

    r23805 r25738  
    3030psString psphotVersionLong(void)
    3131{
    32     psString version = psLibVersion();  // Version, to return
    33     psString source = psLibSource();    // Source
     32    psString version = psphotVersion();  // Version, to return
     33    psString source = psphotSource();    // Source
    3434
    3535    psStringPrepend(&version, "psphot ");

  • trunk/pswarp/src/pswarpLoop.c

    r25518 r25738  
    397397        }
    398398
     399        pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
     400
    399401        // measure the PSF using these sources
    400402        if (!psphotReadoutFindPSF(config, view, sources)) {
Note: See TracChangeset for help on using the changeset viewer.