IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25745


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

merge changes from trunk

Location:
branches/eam_branches/20090715/psModules/src
Files:
17 edited
6 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/astrom/pmAstrometryVisual.c

    r25022 r25745  
    12091209    KapaClearPlots (kapa2);
    12101210
    1211     graphdata.color = KapaColorByName ("black");
    1212     graphdata.ptype = 2;
     1211    KapaSendLabel (kapa2, "X", KAPA_LABEL_XM);
     1212    KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM);
     1213    KapaSendLabel (kapa2, "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars", KAPA_LABEL_XP);
     1214
     1215    // X vs Y by mag (ref)
     1216    graphdata.color = KapaColorByName ("red");
     1217    graphdata.ptype = 7;
    12131218    graphdata.style = 2;
    12141219
    … …  
    12171222    psFree (zVec);
    12181223
    1219     xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
    1220     yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
    1221     zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
    1222 
    1223     // X vs Y by mag (raw)
    1224     n = 0;
    1225     for (int i = 0; i < rawstars->n; i++) {
    1226         pmAstromObj *raw = rawstars->data[i];
    1227         if (!isfinite(raw->Mag)) continue;
    1228         if (raw->Mag < iMagMin) continue;
    1229         if (raw->Mag > iMagMax) continue;
    1230 
    1231         xVec->data.F32[n] = raw->chip->x;
    1232         yVec->data.F32[n] = raw->chip->y;
    1233         zVec->data.F32[n] = raw->Mag;
    1234         n++;
    1235     }
    1236     xVec->n = yVec->n = zVec->n = n;
    1237 
    1238     KapaSendLabel (kapa2, "X", KAPA_LABEL_XM);
    1239     KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM);
    1240     KapaSendLabel (kapa2,
    1241                    "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars"
    1242                    , KAPA_LABEL_XP);
    1243     pmVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false);
    1244 
    1245     // X vs Y by mag (ref)
    1246     psFree (xVec);
    1247     psFree (yVec);
    1248     psFree (zVec);
    1249 
    12501224    xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
    12511225    yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
    12521226    zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
    1253 
    1254     graphdata.color = KapaColorByName ("red");
    1255     graphdata.ptype = 7;
    1256     graphdata.style = 2;
    12571227
    12581228    n = 0;
    … …  
    12691239    }
    12701240    xVec->n = yVec->n = zVec->n = n;
    1271     pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
     1241    pmVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false);
    12721242
    12731243    //rescale the graph to include all points
    … …  
    12821252    graphdata.ymax = PS_MAX(ymax, graphdata.ymax);
    12831253    KapaSetLimits (kapa2, &graphdata);
     1254
     1255    bool plotTweak;
     1256    pmVisualAskUser(&plotTweak);
     1257
     1258    // X vs Y by mag (raw)
     1259    graphdata.color = KapaColorByName ("black");
     1260    graphdata.ptype = 2;
     1261    graphdata.style = 2;
     1262
     1263    psFree (xVec);
     1264    psFree (yVec);
     1265    psFree (zVec);
     1266
     1267    xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1268    yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1269    zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1270
     1271    n = 0;
     1272    for (int i = 0; i < rawstars->n; i++) {
     1273        pmAstromObj *raw = rawstars->data[i];
     1274        if (!isfinite(raw->Mag)) continue;
     1275        if (raw->Mag < iMagMin) continue;
     1276        if (raw->Mag > iMagMax) continue;
     1277
     1278        xVec->data.F32[n] = raw->chip->x;
     1279        yVec->data.F32[n] = raw->chip->y;
     1280        zVec->data.F32[n] = raw->Mag;
     1281        n++;
     1282    }
     1283    xVec->n = yVec->n = zVec->n = n;
     1284    pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
    12841285
    12851286    //overplot matched stars in blue

  • branches/eam_branches/20090715/psModules/src/camera/pmReadoutFake.c

  • branches/eam_branches/20090715/psModules/src/config/pmConfigMask.c

    r23498 r25745  
    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

  • branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c

    r25626 r25745  
    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);

  • branches/eam_branches/20090715/psModules/src/objects/Makefile.am

    r25476 r25745  
    5959        pmGrowthCurve.c \
    6060        pmSourceMatch.c \
    61         pmDetEff.c
    62 
    63 EXTRA_DIST = \
     61        pmDetEff.c \
    6462        models/pmModel_GAUSS.c \
    6563        models/pmModel_PGAUSS.c \
     64        models/pmModel_PS1_V1.c \
    6665        models/pmModel_QGAUSS.c \
    67         models/pmModel_SGAUSS.c \
    6866        models/pmModel_RGAUSS.c \
    6967        models/pmModel_SERSIC.c
    … …  
    9795        pmGrowthCurve.h \
    9896        pmSourceMatch.h \
    99         pmDetEff.h
     97        pmDetEff.h \
     98        models/pmModel_GAUSS.h \
     99        models/pmModel_PGAUSS.h \
     100        models/pmModel_PS1_V1.h \
     101        models/pmModel_QGAUSS.h \
     102        models/pmModel_RGAUSS.h \
     103        models/pmModel_SERSIC.h
    100104
    101105CLEANFILES = *~

  • 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

  • branches/eam_branches/20090715/psModules/src/objects/models/pmModel_PGAUSS.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_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)
    … …  
    7196bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    7297{
    73     float beta_lim = 0, params_min = 0, params_max = 0;
    74     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     98    if (!limitsApply) {
     99        return true;
     100    }
     101    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    75102
    76103    // we need to calculate the limits for SXY specially
     104    float q2 = NAN;
    77105    if (nParam == PM_PAR_SXY) {
    78         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    79         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    80         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    81         q1 = PS_MAX (0.0, q1);
     106        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     107        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     108        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     109        q1 = (q1 < 0.0) ? 0.0 : q1;
    82110        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    83111        // angle and let f2,f1 fight it out
    84         q2  = 0.5*sqrt (q1);
     112        q2 = 0.5*sqrtf(q1);
    85113    }
    86114
    87115    switch (mode) {
    88     case PS_MINIMIZE_BETA_LIMIT:
    89         switch (nParam) {
    90         case PM_PAR_SKY:
    91             beta_lim = 1000;
    92             break;
    93         case PM_PAR_I0:
    94             beta_lim = 3e6;
    95             break;
    96         case PM_PAR_XPOS:
    97             beta_lim = 5;
    98             break;
    99         case PM_PAR_YPOS:
    100             beta_lim = 5;
    101             break;
    102         case PM_PAR_SXX:
    103             beta_lim = 2.0;
    104             break;
    105         case PM_PAR_SYY:
    106             beta_lim = 2.0;
    107             break;
    108         case PM_PAR_SXY:
    109             beta_lim =  0.5*q2;
    110             break;
    111         default:
    112             psAbort("invalid parameter %d for beta test", nParam);
    113         }
    114         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    115             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    116             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    117                      nParam, beta[nParam], beta_lim);
    118             return false;
    119         }
    120         return true;
    121     case PS_MINIMIZE_PARAM_MIN:
    122         switch (nParam) {
    123         case PM_PAR_SKY:
    124             params_min = -1000;
    125             break;
    126         case PM_PAR_I0:
    127             params_min =  0.01;
    128             break;
    129         case PM_PAR_XPOS:
    130             params_min =  -100;
    131             break;
    132         case PM_PAR_YPOS:
    133             params_min =  -100;
    134             break;
    135         case PM_PAR_SXX:
    136             params_min =   0.5;
    137             break;
    138         case PM_PAR_SYY:
    139             params_min =   0.5;
    140             break;
    141         case PM_PAR_SXY:
    142             params_min =   -q2;
    143             break;
    144         default:
    145             psAbort("invalid parameter %d for param min test", nParam);
    146         }
    147         if (params[nParam] < params_min) {
    148             params[nParam] = params_min;
    149             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    150                      nParam, params[nParam], params_min);
    151             return false;
    152         }
    153         return true;
    154     case PS_MINIMIZE_PARAM_MAX:
    155         switch (nParam) {
    156         case PM_PAR_SKY:
    157             params_max =   1e5;
    158             break;
    159         case PM_PAR_I0:
    160             params_max =   1e8;
    161             break;
    162         case PM_PAR_XPOS:
    163             params_max =   1e4;
    164             break;
    165         case PM_PAR_YPOS:
    166             params_max =   1e4;
    167             break;
    168         case PM_PAR_SXX:
    169             params_max =   100;
    170             break;
    171         case PM_PAR_SYY:
    172             params_max =   100;
    173             break;
    174         case PM_PAR_SXY:
    175             params_max =   +q2;
    176             break;
    177         default:
    178             psAbort("invalid parameter %d for param max test", nParam);
    179         }
    180         if (params[nParam] > params_max) {
    181             params[nParam] = params_max;
    182             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    183                      nParam, params[nParam], params_max);
    184             return false;
    185         }
    186         return true;
    187     default:
     116      case PS_MINIMIZE_BETA_LIMIT: {
     117          psAssert(beta, "Require beta to limit beta");
     118          float limit = betaUse[nParam];
     119          if (nParam == PM_PAR_SXY) {
     120              limit *= q2;
     121          }
     122          if (fabs(beta[nParam]) > fabs(limit)) {
     123              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     124              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     125                      nParam, beta[nParam], limit);
     126              return false;
     127          }
     128          return true;
     129      }
     130      case PS_MINIMIZE_PARAM_MIN: {
     131          psAssert(params, "Require parameters to limit parameters");
     132          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     133          float limit = paramsMinUse[nParam];
     134          if (nParam == PM_PAR_SXY) {
     135              limit *= q2;
     136          }
     137          if (params[nParam] < limit) {
     138              params[nParam] = limit;
     139              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     140                      nParam, params[nParam], limit);
     141              return false;
     142          }
     143          return true;
     144      }
     145      case PS_MINIMIZE_PARAM_MAX: {
     146          psAssert(params, "Require parameters to limit parameters");
     147          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     148          float limit = paramsMaxUse[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_max; %g v. %g",
     155                      nParam, params[nParam], limit);
     156              return false;
     157          }
     158          return true;
     159      }
     160      default:
    188161        psAbort("invalid choice for limits");
    189162    }
    … …  
    191164    return false;
    192165}
     166
    193167
    194168// make an initial guess for parameters
    … …  
    432406}
    433407
     408
     409void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     410{
     411    switch (type) {
     412      case PM_MODEL_LIMITS_NONE:
     413        paramsMinUse = NULL;
     414        paramsMaxUse = NULL;
     415        limitsApply = true;
     416        break;
     417      case PM_MODEL_LIMITS_IGNORE:
     418        paramsMinUse = NULL;
     419        paramsMaxUse = NULL;
     420        limitsApply = false;
     421        break;
     422      case PM_MODEL_LIMITS_LAX:
     423        paramsMinUse = paramsMinLax;
     424        paramsMaxUse = paramsMaxLax;
     425        limitsApply = true;
     426        break;
     427      case PM_MODEL_LIMITS_STRICT:
     428        paramsMinUse = paramsMinStrict;
     429        paramsMaxUse = paramsMaxStrict;
     430        limitsApply = true;
     431        break;
     432      default:
     433        psAbort("Unrecognised model limits type: %x", type);
     434    }
     435    return;
     436}
     437
    434438# undef PM_MODEL_FUNC
    435439# undef PM_MODEL_FLUX
    … …  
    440444# undef PM_MODEL_PARAMS_FROM_PSF
    441445# undef PM_MODEL_FIT_STATUS
     446# undef PM_MODEL_SET_LIMITS

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

    r25683 r25745  
    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
    … …  
    3545// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    3646// values need to be pixel coords
     47
     48// Lax parameter limits
     49static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
     50static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     51
     52// Strict parameter limits
     53// k = PAR_7 < 0 is very undesirable (big divot in the middle)
     54static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
     55static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     56
     57// Parameter limits to use
     58static float *paramsMinUse = paramsMinLax;
     59static float *paramsMaxUse = paramsMaxLax;
     60static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     61
     62static bool limitsApply = true;         // Apply limits?
     63
    3764psF32 PM_MODEL_FUNC (psVector *deriv,
    3865                     const psVector *params,
    … …  
    87114bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    88115{
    89     float beta_lim = 0, params_min = 0, params_max = 0;
    90     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     116    if (!limitsApply) {
     117        return true;
     118    }
     119    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    91120
    92121    // we need to calculate the limits for SXY specially
     122    float q2 = NAN;
    93123    if (nParam == PM_PAR_SXY) {
    94         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    95         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    96         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     124        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     125        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     126        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    97127        q1 = (q1 < 0.0) ? 0.0 : q1;
    98128        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    99129        // angle and let f2,f1 fight it out
    100         q2  = 0.5*sqrt (q1);
     130        q2 = 0.5*sqrtf(q1);
    101131    }
    102132
    103133    switch (mode) {
    104     case PS_MINIMIZE_BETA_LIMIT:
    105         switch (nParam) {
    106         case PM_PAR_SKY:
    107             beta_lim = 1000;
    108             break;
    109         case PM_PAR_I0:
    110             beta_lim = 3e6;
    111             break;
    112         case PM_PAR_XPOS:
    113             beta_lim = 5;
    114             break;
    115         case PM_PAR_YPOS:
    116             beta_lim = 5;
    117             break;
    118         case PM_PAR_SXX:
    119             beta_lim = 1.0;
    120             break;
    121         case PM_PAR_SYY:
    122             beta_lim = 1.0;
    123             break;
    124         case PM_PAR_SXY:
    125             beta_lim =  0.5*q2;
    126             break;
    127         case PM_PAR_7:
    128             beta_lim = 2.0;
    129             break;
    130         default:
    131             psAbort("invalid parameter %d for beta test", nParam);
    132         }
    133         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    134             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    135             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    136                      nParam, beta[nParam], beta_lim);
    137             return false;
    138         }
    139         return true;
    140     case PS_MINIMIZE_PARAM_MIN:
    141         switch (nParam) {
    142         case PM_PAR_SKY:
    143             params_min = -1000;
    144             break;
    145         case PM_PAR_I0:
    146             params_min =   0.01;
    147             break;
    148         case PM_PAR_XPOS:
    149             params_min =  -100;
    150             break;
    151         case PM_PAR_YPOS:
    152             params_min =  -100;
    153             break;
    154         case PM_PAR_SXX:
    155             params_min =   0.5;
    156             break;
    157         case PM_PAR_SYY:
    158             params_min =   0.5;
    159             break;
    160         case PM_PAR_SXY:
    161             params_min =  -q2;
    162             break;
    163         case PM_PAR_7:
    164             params_min =  -1.0;
    165             break;
    166         default:
    167             psAbort("invalid parameter %d for param min test", nParam);
    168         }
    169         if (params[nParam] < params_min) {
    170             params[nParam] = params_min;
    171             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    172                      nParam, params[nParam], params_min);
    173             return false;
    174         }
    175         return true;
    176     case PS_MINIMIZE_PARAM_MAX:
    177         switch (nParam) {
    178         case PM_PAR_SKY:
    179             params_max =   1e5;
    180             break;
    181         case PM_PAR_I0:
    182             params_max =   1e8;
    183             break;
    184         case PM_PAR_XPOS:
    185             params_max =   1e4;
    186             break;
    187         case PM_PAR_YPOS:
    188             params_max =   1e4;
    189             break;
    190         case PM_PAR_SXX:
    191             params_max =   100;
    192             break;
    193         case PM_PAR_SYY:
    194             params_max =   100;
    195             break;
    196         case PM_PAR_SXY:
    197             params_max =  +q2;
    198             break;
    199         case PM_PAR_7:
    200             params_max =  20.0;
    201             break;
    202         default:
    203             psAbort("invalid parameter %d for param max test", nParam);
    204         }
    205         if (params[nParam] > params_max) {
    206             params[nParam] = params_max;
    207             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    208                      nParam, params[nParam], params_max);
    209             return false;
    210         }
    211         return true;
    212     default:
     134      case PS_MINIMIZE_BETA_LIMIT: {
     135          psAssert(beta, "Require beta to limit beta");
     136          float limit = betaUse[nParam];
     137          if (nParam == PM_PAR_SXY) {
     138              limit *= q2;
     139          }
     140          if (fabs(beta[nParam]) > fabs(limit)) {
     141              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     142              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     143                      nParam, beta[nParam], limit);
     144              return false;
     145          }
     146          return true;
     147      }
     148      case PS_MINIMIZE_PARAM_MIN: {
     149          psAssert(params, "Require parameters to limit parameters");
     150          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     151          float limit = paramsMinUse[nParam];
     152          if (nParam == PM_PAR_SXY) {
     153              limit *= q2;
     154          }
     155          if (params[nParam] < limit) {
     156              params[nParam] = limit;
     157              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     158                      nParam, params[nParam], limit);
     159              return false;
     160          }
     161          return true;
     162      }
     163      case PS_MINIMIZE_PARAM_MAX: {
     164          psAssert(params, "Require parameters to limit parameters");
     165          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     166          float limit = paramsMaxUse[nParam];
     167          if (nParam == PM_PAR_SXY) {
     168              limit *= q2;
     169          }
     170          if (params[nParam] > limit) {
     171              params[nParam] = limit;
     172              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     173                      nParam, params[nParam], limit);
     174              return false;
     175          }
     176          return true;
     177      }
     178      default:
    213179        psAbort("invalid choice for limits");
    214180    }
    … …  
    302268    psF32 *PAR = params->data.F32;
    303269
    304     if (flux <= 0)
    305         return (1.0);
    306     if (PAR[PM_PAR_I0] <= 0)
    307         return (1.0);
    308     if (flux >= PAR[PM_PAR_I0])
    309         return (1.0);
     270    if (flux <= 0) {
     271        return 1.0;
     272    }
     273    if (PAR[PM_PAR_I0] <= 0) {
     274        return 1.0;
     275    }
     276    if (flux >= PAR[PM_PAR_I0]) {
     277        return 1.0;
     278    }
     279    if (PAR[PM_PAR_7] == 0.0) {
     280        return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
     281    }
    310282
    311283    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    463435}
    464436
     437
     438void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     439{
     440    switch (type) {
     441      case PM_MODEL_LIMITS_NONE:
     442        paramsMinUse = NULL;
     443        paramsMaxUse = NULL;
     444        limitsApply = true;
     445        break;
     446      case PM_MODEL_LIMITS_IGNORE:
     447        paramsMinUse = NULL;
     448        paramsMaxUse = NULL;
     449        limitsApply = false;
     450        break;
     451      case PM_MODEL_LIMITS_LAX:
     452        paramsMinUse = paramsMinLax;
     453        paramsMaxUse = paramsMaxLax;
     454        limitsApply = true;
     455        break;
     456      case PM_MODEL_LIMITS_STRICT:
     457        paramsMinUse = paramsMinStrict;
     458        paramsMaxUse = paramsMaxStrict;
     459        limitsApply = true;
     460        break;
     461      default:
     462        psAbort("Unrecognised model limits type: %x", type);
     463    }
     464    return;
     465}
     466
     467
    465468# undef PM_MODEL_FUNC
    466469# undef PM_MODEL_FLUX
    … …  
    471474# undef PM_MODEL_PARAMS_FROM_PSF
    472475# undef PM_MODEL_FIT_STATUS
     476# undef PM_MODEL_SET_LIMITS
    473477# undef ALPHA
    474478# undef ALPHA_M

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

    r25683 r25745  
    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
    3040
    3141// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    3242// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    3343// values need to be pixel coords
     44
     45// Lax parameter limits
     46static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
     47static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
     48
     49// Strict parameter limits
     50static float *paramsMinStrict = paramsMinLax;
     51static float *paramsMaxStrict = paramsMaxLax;
     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 };
     57
     58static bool limitsApply = true;         // Apply limits?
     59
    3460psF32 PM_MODEL_FUNC (psVector *deriv,
    3561                     const psVector *params,
    … …  
    82108# define AR_MAX 20.0
    83109# define AR_RATIO 0.99
     110
    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 =   0.1;
    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;
     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      }
    209175    default:
    210176        psAbort("invalid choice for limits");
    … …  
    461427}
    462428
     429
     430void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     431{
     432    switch (type) {
     433      case PM_MODEL_LIMITS_NONE:
     434        paramsMinUse = NULL;
     435        paramsMaxUse = NULL;
     436        limitsApply = true;
     437        break;
     438      case PM_MODEL_LIMITS_IGNORE:
     439        paramsMinUse = NULL;
     440        paramsMaxUse = NULL;
     441        limitsApply = false;
     442        break;
     443      case PM_MODEL_LIMITS_LAX:
     444        paramsMinUse = paramsMinLax;
     445        paramsMaxUse = paramsMaxLax;
     446        limitsApply = true;
     447        break;
     448      case PM_MODEL_LIMITS_STRICT:
     449        paramsMinUse = paramsMinStrict;
     450        paramsMaxUse = paramsMaxStrict;
     451        limitsApply = true;
     452        break;
     453      default:
     454        psAbort("Unrecognised model limits type: %x", type);
     455    }
     456    return;
     457}
     458
    463459# undef PM_MODEL_FUNC
    464460# undef PM_MODEL_FLUX
    … …  
    469465# undef PM_MODEL_PARAMS_FROM_PSF
    470466# undef PM_MODEL_FIT_STATUS
     467# undef PM_MODEL_SET_LIMITS

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

    r25683 r25745  
    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
    3040
    3141// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    3242// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    3343// values need to be pixel coords
     44
     45// Lax parameter limits
     46static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
     47static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     48
     49// Strict parameter limits
     50static float *paramsMinStrict = paramsMinLax;
     51static float *paramsMaxStrict = paramsMaxLax;
     52
     53// Parameter limits to use
     54static float *paramsMinUse = paramsMinLax;
     55static float *paramsMaxUse = paramsMaxLax;
     56static float betaUse[] = { 1000, 3e6, 5, 5, 0.5, 0.5, 0.5, 0.5 };
     57
     58static bool limitsApply = true;         // Apply limits?
     59
    3460psF32 PM_MODEL_FUNC (psVector *deriv,
    3561                     const psVector *params,
    … …  
    76102# define AR_MAX 20.0
    77103# define AR_RATIO 0.99
     104
    78105bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    79106{
    80     float beta_lim = 0, params_min = 0, params_max = 0;
    81     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     107    if (!limitsApply) {
     108        return true;
     109    }
     110    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    82111
    83112    // we need to calculate the limits for SXY specially
     113    float q2 = NAN;
    84114    if (nParam == PM_PAR_SXY) {
    85         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    86         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    87         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     115        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     116        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     117        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    88118        q1 = (q1 < 0.0) ? 0.0 : q1;
    89119        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    90120        // angle and let f2,f1 fight it out
    91         q2  = 0.5*sqrt (q1);
     121        q2 = 0.5*sqrtf(q1);
    92122    }
    93123
    94124    switch (mode) {
    95     case PS_MINIMIZE_BETA_LIMIT:
    96         switch (nParam) {
    97         case PM_PAR_SKY:
    98             beta_lim = 1000;
    99             break;
    100         case PM_PAR_I0:
    101             beta_lim = 3e6;
    102             break;
    103         case PM_PAR_XPOS:
    104             beta_lim = 5;
    105             break;
    106         case PM_PAR_YPOS:
    107             beta_lim = 5;
    108             break;
    109         case PM_PAR_SXX:
    110             beta_lim = 0.5;
    111             break;
    112         case PM_PAR_SYY:
    113             beta_lim = 0.5;
    114             break;
    115         case PM_PAR_SXY:
    116             beta_lim =  0.5*q2;
    117             break;
    118         case PM_PAR_7:
    119             beta_lim = 0.5;
    120             break;
    121         default:
    122             psAbort("invalid parameter %d for beta test", nParam);
    123         }
    124         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    125             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    126             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    127                      nParam, beta[nParam], beta_lim);
    128             return false;
    129         }
    130         return true;
    131     case PS_MINIMIZE_PARAM_MIN:
    132         switch (nParam) {
    133         case PM_PAR_SKY:
    134             params_min = -1000;
    135             break;
    136         case PM_PAR_I0:
    137             params_min =   0.01;
    138             break;
    139         case PM_PAR_XPOS:
    140             params_min =  -100;
    141             break;
    142         case PM_PAR_YPOS:
    143             params_min =  -100;
    144             break;
    145         case PM_PAR_SXX:
    146             params_min =   0.5;
    147             break;
    148         case PM_PAR_SYY:
    149             params_min =   0.5;
    150             break;
    151         case PM_PAR_SXY:
    152             params_min =  -q2;
    153             break;
    154         case PM_PAR_7:
    155             params_min =   1.25;
    156             break;
    157         default:
    158             psAbort("invalid parameter %d for param min test", nParam);
    159         }
    160         if (params[nParam] < params_min) {
    161             params[nParam] = params_min;
    162             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    163                      nParam, params[nParam], params_min);
    164             return false;
    165         }
    166         return true;
    167     case PS_MINIMIZE_PARAM_MAX:
    168         switch (nParam) {
    169         case PM_PAR_SKY:
    170             params_max =   1e5;
    171             break;
    172         case PM_PAR_I0:
    173             params_max =   1e8;
    174             break;
    175         case PM_PAR_XPOS:
    176             params_max =   1e4;
    177             break;
    178         case PM_PAR_YPOS:
    179             params_max =   1e4;
    180             break;
    181         case PM_PAR_SXX:
    182             params_max =   100;
    183             break;
    184         case PM_PAR_SYY:
    185             params_max =   100;
    186             break;
    187         case PM_PAR_SXY:
    188             params_max =  +q2;
    189             break;
    190         case PM_PAR_7:
    191             params_max =  4.0;
    192             break;
    193         default:
    194             psAbort("invalid parameter %d for param max test", nParam);
    195         }
    196         if (params[nParam] > params_max) {
    197             params[nParam] = params_max;
    198             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    199                      nParam, params[nParam], params_max);
    200             return false;
    201         }
    202         return true;
    203     default:
     125      case PS_MINIMIZE_BETA_LIMIT: {
     126          psAssert(beta, "Require beta to limit beta");
     127          float limit = betaUse[nParam];
     128          if (nParam == PM_PAR_SXY) {
     129              limit *= q2;
     130          }
     131          if (fabs(beta[nParam]) > fabs(limit)) {
     132              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     133              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     134                      nParam, beta[nParam], limit);
     135              return false;
     136          }
     137          return true;
     138      }
     139      case PS_MINIMIZE_PARAM_MIN: {
     140          psAssert(params, "Require parameters to limit parameters");
     141          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     142          float limit = paramsMinUse[nParam];
     143          if (nParam == PM_PAR_SXY) {
     144              limit *= q2;
     145          }
     146          if (params[nParam] < limit) {
     147              params[nParam] = limit;
     148              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     149                      nParam, params[nParam], limit);
     150              return false;
     151          }
     152          return true;
     153      }
     154      case PS_MINIMIZE_PARAM_MAX: {
     155          psAssert(params, "Require parameters to limit parameters");
     156          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     157          float limit = paramsMaxUse[nParam];
     158          if (nParam == PM_PAR_SXY) {
     159              limit *= q2;
     160          }
     161          if (params[nParam] > limit) {
     162              params[nParam] = limit;
     163              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     164                      nParam, params[nParam], limit);
     165              return false;
     166          }
     167          return true;
     168      }
     169      default:
    204170        psAbort("invalid choice for limits");
    205171    }
    … …  
    207173    return false;
    208174}
     175
    209176
    210177// make an initial guess for parameters
    … …  
    454421}
    455422
     423
     424void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     425{
     426    switch (type) {
     427      case PM_MODEL_LIMITS_NONE:
     428        paramsMinUse = NULL;
     429        paramsMaxUse = NULL;
     430        limitsApply = true;
     431        break;
     432      case PM_MODEL_LIMITS_IGNORE:
     433        paramsMinUse = NULL;
     434        paramsMaxUse = NULL;
     435        limitsApply = false;
     436        break;
     437      case PM_MODEL_LIMITS_LAX:
     438        paramsMinUse = paramsMinLax;
     439        paramsMaxUse = paramsMaxLax;
     440        limitsApply = true;
     441        break;
     442      case PM_MODEL_LIMITS_STRICT:
     443        paramsMinUse = paramsMinStrict;
     444        paramsMaxUse = paramsMaxStrict;
     445        limitsApply = true;
     446        break;
     447      default:
     448        psAbort("Unrecognised model limits type: %x", type);
     449    }
     450    return;
     451}
     452
    456453# undef PM_MODEL_FUNC
    457454# undef PM_MODEL_FLUX
    … …  
    462459# undef PM_MODEL_PARAMS_FROM_PSF
    463460# undef PM_MODEL_FIT_STATUS
     461# undef PM_MODEL_SET_LIMITS

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

    r25683 r25745  
    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
    3343
    3444// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    3545// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    3646// values need to be pixel coords
     47
     48// Lax parameter limits
     49static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
     50static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     51
     52// Strict parameter limits
     53static float *paramsMinStrict = paramsMinLax;
     54static float *paramsMaxStrict = paramsMaxLax;
     55
     56// Parameter limits to use
     57static float *paramsMinUse = paramsMinLax;
     58static float *paramsMaxUse = paramsMaxLax;
     59static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     60
     61static bool limitsApply = true;         // Apply limits?
     62
    3763psF32 PM_MODEL_FUNC (psVector *deriv,
    3864                     const psVector *params,
    … …  
    94120bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    95121{
    96     float beta_lim = 0, params_min = 0, params_max = 0;
    97     float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
     122    if (!limitsApply) {
     123        return true;
     124    }
     125    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
    98126
    99127    // we need to calculate the limits for SXY specially
     128    float q2 = NAN;
    100129    if (nParam == PM_PAR_SXY) {
    101         f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    102         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
    103         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
     130        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
     131        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     132        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    104133        q1 = (q1 < 0.0) ? 0.0 : q1;
    105134        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    106135        // angle and let f2,f1 fight it out
    107         q2  = 0.5*sqrt (q1);
     136        q2 = 0.5*sqrtf(q1);
    108137    }
    109138
    110139    switch (mode) {
    111     case PS_MINIMIZE_BETA_LIMIT:
    112         switch (nParam) {
    113         case PM_PAR_SKY:
    114             beta_lim = 1000;
    115             break;
    116         case PM_PAR_I0:
    117             beta_lim = 3e6;
    118             break;
    119         case PM_PAR_XPOS:
    120             beta_lim = 5;
    121             break;
    122         case PM_PAR_YPOS:
    123             beta_lim = 5;
    124             break;
    125         case PM_PAR_SXX:
    126             beta_lim = 1.0;
    127             break;
    128         case PM_PAR_SYY:
    129             beta_lim = 1.0;
    130             break;
    131         case PM_PAR_SXY:
    132             beta_lim =  0.5*q2;
    133             break;
    134         case PM_PAR_7:
    135             beta_lim = 2.0;
    136             break;
    137         default:
    138             psAbort("invalid parameter %d for beta test", nParam);
    139         }
    140         if (fabs(beta[nParam]) > fabs(beta_lim)) {
    141             beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
    142             psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
    143                      nParam, beta[nParam], beta_lim);
    144             return false;
    145         }
    146         return true;
    147     case PS_MINIMIZE_PARAM_MIN:
    148         switch (nParam) {
    149         case PM_PAR_SKY:
    150             params_min = -1000;
    151             break;
    152         case PM_PAR_I0:
    153             params_min =     0.01;
    154             break;
    155         case PM_PAR_XPOS:
    156             params_min =  -100;
    157             break;
    158         case PM_PAR_YPOS:
    159             params_min =  -100;
    160             break;
    161         case PM_PAR_SXX:
    162             params_min =   0.05;
    163             break;
    164         case PM_PAR_SYY:
    165             params_min =   0.05;
    166             break;
    167         case PM_PAR_SXY:
    168             params_min =  -q2;
    169             break;
    170         case PM_PAR_7:
    171             params_min =   0.05;
    172             break;
    173         default:
    174             psAbort("invalid parameter %d for param min test", nParam);
    175         }
    176         if (params[nParam] < params_min) {
    177             params[nParam] = params_min;
    178             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
    179                      nParam, params[nParam], params_min);
    180             return false;
    181         }
    182         return true;
    183     case PS_MINIMIZE_PARAM_MAX:
    184         switch (nParam) {
    185         case PM_PAR_SKY:
    186             params_max =   1e5;
    187             break;
    188         case PM_PAR_I0:
    189             params_max =   1e8;
    190             break;
    191         case PM_PAR_XPOS:
    192             params_max =   1e4;
    193             break;
    194         case PM_PAR_YPOS:
    195             params_max =   1e4;
    196             break;
    197         case PM_PAR_SXX:
    198             params_max =   100;
    199             break;
    200         case PM_PAR_SYY:
    201             params_max =   100;
    202             break;
    203         case PM_PAR_SXY:
    204             params_max =  +q2;
    205             break;
    206         case PM_PAR_7:
    207             params_max =   4.0;
    208             break;
    209         default:
    210             psAbort("invalid parameter %d for param max test", nParam);
    211         }
    212         if (params[nParam] > params_max) {
    213             params[nParam] = params_max;
    214             psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
    215                      nParam, params[nParam], params_max);
    216             return false;
    217         }
    218         return true;
    219     default:
     140      case PS_MINIMIZE_BETA_LIMIT: {
     141          psAssert(beta, "Require beta to limit beta");
     142          float limit = betaUse[nParam];
     143          if (nParam == PM_PAR_SXY) {
     144              limit *= q2;
     145          }
     146          if (fabs(beta[nParam]) > fabs(limit)) {
     147              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
     148              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
     149                      nParam, beta[nParam], limit);
     150              return false;
     151          }
     152          return true;
     153      }
     154      case PS_MINIMIZE_PARAM_MIN: {
     155          psAssert(params, "Require parameters to limit parameters");
     156          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
     157          float limit = paramsMinUse[nParam];
     158          if (nParam == PM_PAR_SXY) {
     159              limit *= q2;
     160          }
     161          if (params[nParam] < limit) {
     162              params[nParam] = limit;
     163              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
     164                      nParam, params[nParam], limit);
     165              return false;
     166          }
     167          return true;
     168      }
     169      case PS_MINIMIZE_PARAM_MAX: {
     170          psAssert(params, "Require parameters to limit parameters");
     171          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
     172          float limit = paramsMaxUse[nParam];
     173          if (nParam == PM_PAR_SXY) {
     174              limit *= q2;
     175          }
     176          if (params[nParam] > limit) {
     177              params[nParam] = limit;
     178              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
     179                      nParam, params[nParam], limit);
     180              return false;
     181          }
     182          return true;
     183      }
     184      default:
    220185        psAbort("invalid choice for limits");
    221186    }
    … …  
    223188    return false;
    224189}
    225 
    226190
    227191// make an initial guess for parameters
    … …  
    446410}
    447411
     412
     413void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
     414{
     415    switch (type) {
     416      case PM_MODEL_LIMITS_NONE:
     417        paramsMinUse = NULL;
     418        paramsMaxUse = NULL;
     419        limitsApply = true;
     420        break;
     421      case PM_MODEL_LIMITS_IGNORE:
     422        paramsMinUse = NULL;
     423        paramsMaxUse = NULL;
     424        limitsApply = false;
     425        break;
     426      case PM_MODEL_LIMITS_LAX:
     427        paramsMinUse = paramsMinLax;
     428        paramsMaxUse = paramsMaxLax;
     429        limitsApply = true;
     430        break;
     431      case PM_MODEL_LIMITS_STRICT:
     432        paramsMinUse = paramsMinStrict;
     433        paramsMaxUse = paramsMaxStrict;
     434        limitsApply = true;
     435        break;
     436      default:
     437        psAbort("Unrecognised model limits type: %x", type);
     438    }
     439    return;
     440}
     441
    448442# undef PM_MODEL_FUNC
    449443# undef PM_MODEL_FLUX
    … …  
    454448# undef PM_MODEL_PARAMS_FROM_PSF
    455449# undef PM_MODEL_FIT_STATUS
     450# undef PM_MODEL_SET_LIMITS

  • branches/eam_branches/20090715/psModules/src/objects/pmModel.c

    r25503 r25745  
    8787    tmp->modelParamsFromPSF = class->modelParamsFromPSF;
    8888    tmp->modelFitStatus     = class->modelFitStatus;
     89    tmp->modelSetLimits     = class->modelSetLimits;
    8990
    9091    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
    … …  
    256257                continue;
    257258
    258             // Convert to coordinate in parent image, with offset (dx,dy)
     259            // Convert to coordinate in parent image, with offset (dx,dy)
    259260            // 0.5 PIX: the model take pixel coordinates so convert the pixel index here
    260261            imageCol = ix + 0.5 + image->col0 - dx;
    … …  
    277278                float rx = xBin*ix + DX;
    278279
    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);
     280                int rx0 = rx - 0.5;
     281                int rx1 = rx + 0.5;
     282                int ry0 = ry - 0.5;
     283                int ry1 = ry + 0.5;
     284
     285                if (rx0 < 0) goto skip;
     286                if (ry0 < 0) goto skip;
     287                if (rx1 >= NX) goto skip;
     288                if (ry1 >= NY) goto skip;
     289
     290                // these go from 0.0 to 1.0 between the centers of the pixels
     291                float fx = rx - 0.5 - rx0;
     292                float Fx = 1.0 - fx;
     293                float fy = ry - 0.5 - ry0;
     294                float Fy = 1.0 - fy;
     295
     296                // check the residual image mask (if set). give up if any of the 4 pixels are masked.
     297                if (Rm) {
     298                    if (Rm[ry0][rx0]) goto skip;
     299                    if (Rm[ry0][rx1]) goto skip;
     300                    if (Rm[ry1][rx0]) goto skip;
     301                    if (Rm[ry1][rx1]) goto skip;
     302                }
     303
     304                // a possible further optimization if we re-use these values
     305                // XXX allow for masked pixels, and add pixel weights
     306                float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
     307                float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
     308                float Vo = V0*Fy + V1*fy;
     309                if (!isfinite(Vo)) goto skip;
     310
     311                float Vx = 0.0;
     312                float Vy = 0.0;
     313
     314                // skip Rx,Ry if Ro is masked
     315                if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
     316                    V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
     317                    V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
     318                    Vx = V0*Fy + V1*fy;
     319
     320                    V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
     321                    V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
     322                    Vy = V0*Fy + V1*fy;
     323                }
     324                if (!isfinite(Vx)) goto skip;
     325                if (!isfinite(Vy)) goto skip;
     326
     327                // 2D residual variations are set for the true source position
     328                pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
    328329            }
    329330
    330         skip:
     331        skip:
    331332            // add or subtract the value
    332333            if (add) {

  • branches/eam_branches/20090715/psModules/src/objects/pmModel.h

    r25496 r25745  
    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 */

  • branches/eam_branches/20090715/psModules/src/objects/pmModelClass.c

    r25349 r25745  
    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

  • branches/eam_branches/20090715/psModules/src/objects/pmModelClass.h

    r15697 r25745  
    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 */

  • branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.c

  • branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.h

Note: See TracChangeset for help on using the changeset viewer.