IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 5:10:19 PM (17 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/20090820
Files:
6 deleted
37 edited
15 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090820

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

    r24875 r25766  
    5050        pmPSF_IO.c \
    5151        pmPSFtry.c \
     52        pmPSFtryModel.c \
     53        pmPSFtryFitEXT.c \
     54        pmPSFtryMakePSF.c \
     55        pmPSFtryFitPSF.c \
     56        pmPSFtryMetric.c \
    5257        pmTrend2D.c \
    5358        pmGrowthCurveGenerate.c \
    5459        pmGrowthCurve.c \
    55         pmSourceMatch.c
    56 
    57 EXTRA_DIST = \
     60        pmSourceMatch.c \
     61        pmDetEff.c \
    5862        models/pmModel_GAUSS.c \
    5963        models/pmModel_PGAUSS.c \
     64        models/pmModel_PS1_V1.c \
    6065        models/pmModel_QGAUSS.c \
    61         models/pmModel_SGAUSS.c \
    6266        models/pmModel_RGAUSS.c \
    6367        models/pmModel_SERSIC.c
    … …  
    9094        pmTrend2D.h \
    9195        pmGrowthCurve.h \
    92         pmSourceMatch.h
     96        pmSourceMatch.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
    93104
    94105CLEANFILES = *~

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

    r20001 r25766  
    11/******************************************************************************
    22 * this file defines the GAUSS source shape model.  Note that these model functions are loaded
    3  * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
     3 * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
    44 * models use a psVector to represent the set of parameters, with the sequence used to specify
    55 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
    6  * specifics of the model.  All models which are used a PSF representations share a few
     6 * specifics of the model.  All models which are used as a PSF representations share a few
    77 * parameters, for which # define names are listed in pmModel.h:
    88
    … …  
    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)
     56// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     57// values need to be pixel coords
    3158psF32 PM_MODEL_FUNC(psVector *deriv,
    3259                    const psVector *params,
    … …  
    6895bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    6996{
    70     float beta_lim = 0, params_min = 0, params_max = 0;
    71     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");
    72101
    73102    // we need to calculate the limits for SXY specially
     103    float q2 = NAN;
    74104    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);
     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;
    79109        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    80110        // angle and let f2,f1 fight it out
    81         q2  = 0.5*sqrt (q1);
     111        q2 = 0.5*sqrtf(q1);
    82112    }
    83113
    84114    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;
     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      }
    184159    default:
    185160        psAbort("invalid choice for limits");
    … …  
    190165
    191166// make an initial guess for parameters
     167// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    192168bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    193169{
    … …  
    205181    psEllipseShape shape = psEllipseAxesToShape (axes);
    206182
    207     PAR[PM_PAR_SKY]  = moments->Sky;
     183    PAR[PM_PAR_SKY]  = 0.0;
    208184    PAR[PM_PAR_I0]   = peak->flux;
    209185    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    257233    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    258234    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
     235    psAssert (isfinite(radius), "fix this code: radius should not be nan for %f", PAR[PM_PAR_I0]);
     236
    259237    return (radius);
    260238}
    … …  
    367345bool PM_MODEL_FIT_STATUS (pmModel *model)
    368346{
    369     psF32 dP;
    370347    bool  status;
    371348
    … …  
    373350    psF32 *dPAR = model->dparams->data.F32;
    374351
    375     dP = 0;
    376     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    377     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    378     dP = sqrt (dP);
    379 
    380352    status = true;
    381     status &= (dP < 0.5);
    382353    status &= (PAR[PM_PAR_I0] > 0);
    383354    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    384355
    385     if (status)
    386         return true;
    387     return false;
     356    return status;
     357}
     358
     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;
    388386}
    389387
    … …  
    396394# undef PM_MODEL_PARAMS_FROM_PSF
    397395# undef PM_MODEL_FIT_STATUS
     396# undef PM_MODEL_SET_LIMITS

  • branches/eam_branches/20090820/psModules/src/objects/models/pmModel_PGAUSS.c

    r20001 r25766  
    11/******************************************************************************
    22 * this file defines the PGAUSS source shape model.  Note that these model functions are loaded
    3  * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
     3 * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
    44 * models use a psVector to represent the set of parameters, with the sequence used to specify
    55 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
    6  * specifics of the model.  All models which are used a PSF representations share a few
     6 * specifics of the model.  All models which are used as a PSF representations share a few
    77 * parameters, for which # define names are listed in pmModel.h:
    88
    … …  
    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)
     56// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     57// values need to be pixel coords
    3158psF32 PM_MODEL_FUNC(psVector *deriv,
    3259                    const psVector *params,
    … …  
    6996bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    7097{
    71     float beta_lim = 0, params_min = 0, params_max = 0;
    72     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");
    73102
    74103    // we need to calculate the limits for SXY specially
     104    float q2 = NAN;
    75105    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);
     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;
    80110        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    81111        // angle and let f2,f1 fight it out
    82         q2  = 0.5*sqrt (q1);
     112        q2 = 0.5*sqrtf(q1);
    83113    }
    84114
    85115    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:
     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:
    186161        psAbort("invalid choice for limits");
    187162    }
    … …  
    190165}
    191166
     167
    192168// make an initial guess for parameters
     169// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    193170bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    194171{
    … …  
    205182    psEllipseShape shape = psEllipseAxesToShape (axes);
    206183
    207     PAR[PM_PAR_SKY]  = moments->Sky;
     184    PAR[PM_PAR_SKY]  = 0.0;
    208185    PAR[PM_PAR_I0]   = peak->flux;
    209186    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    284261    // choose a z value guaranteed to be beyond our limit
    285262    float z0 = pow((1.0 / limit), (1.0 / 3.0));
     263    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_I0]);
    286264    float z1 = (1.0 / limit);
     265    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_I0]);
    287266    z1 = PS_MAX (z0, z1);
    288267    z0 = 0.0;
    … …  
    415394bool PM_MODEL_FIT_STATUS (pmModel *model)
    416395{
    417     psF32 dP;
    418396    bool  status;
    419397
    … …  
    421399    psF32 *dPAR = model->dparams->data.F32;
    422400
    423     dP = 0;
    424     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    425     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    426     dP = sqrt (dP);
    427 
    428401    status = true;
    429     status &= (dP < 0.5);
    430402    status &= (PAR[PM_PAR_I0] > 0);
    431403    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    432404
    433405    return status;
     406}
     407
     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;
    434436}
    435437
    … …  
    442444# undef PM_MODEL_PARAMS_FROM_PSF
    443445# undef PM_MODEL_FIT_STATUS
     446# undef PM_MODEL_SET_LIMITS

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

    r23962 r25766  
    11/******************************************************************************
    2  * this file defines the PS1_V1 source shape model (XXX need a better name!).  Note that these
    3  * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
    4  * statements of their own.  The models use a psVector to represent the set of parameters, with
    5  * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
    6  * may thus vary depending on the specifics of the model.  All models which are used a PSF
    7  * representations share a few parameters, for which # define names are listed in pmModel.h:
     2 * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
     3 * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
     4 * models use a psVector to represent the set of parameters, with the sequence used to specify
     5 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
     6 * specifics of the model.  All models which are used as a PSF representations share a few
     7 * parameters, for which # define names are listed in pmModel.h:
    88
    99   power-law with fitted linear term
    … …  
    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// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
     45// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     46// 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?
    3363
    3464psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    84114bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    85115{
    86     float beta_lim = 0, params_min = 0, params_max = 0;
    87     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");
    88120
    89121    // we need to calculate the limits for SXY specially
     122    float q2 = NAN;
    90123    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);
     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);
    94127        q1 = (q1 < 0.0) ? 0.0 : q1;
    95128        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    96129        // angle and let f2,f1 fight it out
    97         q2  = 0.5*sqrt (q1);
     130        q2 = 0.5*sqrtf(q1);
    98131    }
    99132
    100133    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:
     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:
    210179        psAbort("invalid choice for limits");
    211180    }
    … …  
    216185
    217186// make an initial guess for parameters
     187// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    218188bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    219189{
    … …  
    240210    if (!isfinite(shape.sxy)) return false;
    241211
    242     // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
    243212    PAR[PM_PAR_SKY]  = 0.0;
    244213    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    299268    psF32 *PAR = params->data.F32;
    300269
    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);
     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    }
    307282
    308283    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    320295
    321296    // choose a z value guaranteed to be beyond our limit
    322     float z0 = pow((1.0 / limit), (1.0 / ALPHA));
    323     float z1 = (1.0 / limit) / PAR[PM_PAR_7];
    324     z1 = PS_MAX (z0, z1);
    325     z0 = 0.0;
     297    float z0 = 0.0;
     298    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
     299    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
     300    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
    326301
    327302    // perform a type of bisection to find the value
    … …  
    448423bool PM_MODEL_FIT_STATUS (pmModel *model)
    449424{
    450 
    451     psF32 dP;
    452425    bool  status;
    453426
    … …  
    455428    psF32 *dPAR = model->dparams->data.F32;
    456429
    457     dP = 0;
    458     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    459     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    460     dP = sqrt (dP);
    461 
    462430    status = true;
    463 //    status &= (dP < 0.5);
    464431    status &= (PAR[PM_PAR_I0] > 0);
    465432    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    … …  
    467434    return status;
    468435}
     436
     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
    469467
    470468# undef PM_MODEL_FUNC
    … …  
    476474# undef PM_MODEL_PARAMS_FROM_PSF
    477475# undef PM_MODEL_FIT_STATUS
     476# undef PM_MODEL_SET_LIMITS
    478477# undef ALPHA
    479478# undef ALPHA_M

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

    r20001 r25766  
    11/******************************************************************************
    22 * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
    3  * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
     3 * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
    44 * statements of their own.  The models use a psVector to represent the set of parameters, with
    55 * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
    … …  
    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// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
     42// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     43// 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?
    3059
    3160psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    79108# define AR_MAX 20.0
    80109# define AR_RATIO 0.99
     110
    81111bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    82112{
    83     float beta_lim = 0, params_min = 0, params_max = 0;
    84     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");
    85117
    86118    // we need to calculate the limits for SXY specially
     119    float q2 = NAN;
    87120    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);
     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);
    91124        q1 = (q1 < 0.0) ? 0.0 : q1;
    92125        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    93126        // angle and let f2,f1 fight it out
    94         q2  = 0.5*sqrt (q1);
     127        q2 = 0.5*sqrtf(q1);
    95128    }
    96129
    97130    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;
     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      }
    206175    default:
    207176        psAbort("invalid choice for limits");
    … …  
    213182
    214183// make an initial guess for parameters
     184// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    215185bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    216186{
    … …  
    237207    if (!isfinite(shape.sxy)) return false;
    238208
    239     // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
    240209    PAR[PM_PAR_SKY]  = 0.0;
    241210    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    318287    // choose a z value guaranteed to be beyond our limit
    319288    float z0 = pow((1.0 / limit), (1.0 / 2.25));
     289    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
    320290    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
     291    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
    321292    z1 = PS_MAX (z0, z1);
    322293    z0 = 0.0;
    … …  
    444415bool PM_MODEL_FIT_STATUS (pmModel *model)
    445416{
    446 
    447     psF32 dP;
    448417    bool  status;
    449418
    … …  
    451420    psF32 *dPAR = model->dparams->data.F32;
    452421
    453     dP = 0;
    454     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    455     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    456     dP = sqrt (dP);
    457 
    458422    status = true;
    459 //    status &= (dP < 0.5);
    460423    status &= (PAR[PM_PAR_I0] > 0);
    461424    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    462425
    463426    return status;
     427}
     428
     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;
    464457}
    465458
    … …  
    472465# undef PM_MODEL_PARAMS_FROM_PSF
    473466# undef PM_MODEL_FIT_STATUS
     467# undef PM_MODEL_SET_LIMITS

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

    r20001 r25766  
    11/******************************************************************************
    22 * this file defines the RGAUSS source shape model (XXX need a better name!).  Note that these
    3  * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
     3 * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
    44 * statements of their own.  The models use a psVector to represent the set of parameters, with
    55 * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
    6  * may thus vary depending on the specifics of the model.  All models which are used a PSF
     6 * may thus vary depending on the specifics of the model.  All models which are used as a PSF
    77 * representations share a few parameters, for which # define names are listed in pmModel.h:
    88
    … …  
    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// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
     42// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     43// 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?
    3059
    3160psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    6291        dPAR[PM_PAR_SXY] = -q*X*Y;
    6392
    64         // this model derivative is undefined at z = 0.0, but is actually 0.0
     93        // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0
    6594        dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -5.0*t*log(z)*p*z;
    6695    }
    … …  
    73102# define AR_MAX 20.0
    74103# define AR_RATIO 0.99
     104
    75105bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    76106{
    77     float beta_lim = 0, params_min = 0, params_max = 0;
    78     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");
    79111
    80112    // we need to calculate the limits for SXY specially
     113    float q2 = NAN;
    81114    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);
     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);
    85118        q1 = (q1 < 0.0) ? 0.0 : q1;
    86119        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    87120        // angle and let f2,f1 fight it out
    88         q2  = 0.5*sqrt (q1);
     121        q2 = 0.5*sqrtf(q1);
    89122    }
    90123
    91124    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:
     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:
    201170        psAbort("invalid choice for limits");
    202171    }
    … …  
    205174}
    206175
     176
    207177// make an initial guess for parameters
     178// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    208179bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    209180{
    … …  
    230201    if (!isfinite(shape.sxy)) return false;
    231202
    232     PAR[PM_PAR_SKY]  = moments->Sky;
     203    PAR[PM_PAR_SKY]  = 0.0;
    233204    PAR[PM_PAR_I0]   = peak->flux;
    234205    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    310281    // choose a z value guaranteed to be beyond our limit
    311282    float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7]));
     283    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
    312284    float z1 = (1.0 / limit);
     285    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
    313286    z1 = PS_MAX (z0, z1);
    314287    z0 = 0.0;
    … …  
    436409bool PM_MODEL_FIT_STATUS (pmModel *model)
    437410{
    438 
    439     psF32 dP;
    440411    bool  status;
    441412
    … …  
    443414    psF32 *dPAR = model->dparams->data.F32;
    444415
    445     dP = 0;
    446     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    447     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    448     dP = sqrt (dP);
    449 
    450416    status = true;
    451     status &= (dP < 0.5);
    452417    status &= (PAR[PM_PAR_I0] > 0);
    453418    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    454419
    455420    return status;
     421}
     422
     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;
    456451}
    457452
    … …  
    464459# undef PM_MODEL_PARAMS_FROM_PSF
    465460# undef PM_MODEL_FIT_STATUS
     461# undef PM_MODEL_SET_LIMITS

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

    r20001 r25766  
    11/******************************************************************************
    22 * this file defines the SERSIC source shape model.  Note that these model functions are loaded
    3  * by pmModelGroup.c using 'include', and thus need no 'include' statements of their own.  The
     3 * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
    44 * models use a psVector to represent the set of parameters, with the sequence used to specify
    55 * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
    6  * specifics of the model.  All models which are used a PSF representations share a few
     6 * specifics of the model.  All models which are used as a PSF representations share a few
    77 * parameters, for which # define names are listed in pmModel.h:
    88
    … …  
    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// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
     45// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
     46// 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?
    3362
    3463psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    91120bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    92121{
    93     float beta_lim = 0, params_min = 0, params_max = 0;
    94     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");
    95126
    96127    // we need to calculate the limits for SXY specially
     128    float q2 = NAN;
    97129    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);
     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);
    101133        q1 = (q1 < 0.0) ? 0.0 : q1;
    102134        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    103135        // angle and let f2,f1 fight it out
    104         q2  = 0.5*sqrt (q1);
     136        q2 = 0.5*sqrtf(q1);
    105137    }
    106138
    107139    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:
     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:
    217185        psAbort("invalid choice for limits");
    218186    }
    … …  
    221189}
    222190
    223 
    224191// make an initial guess for parameters
     192// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    225193bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    226194{
    … …  
    247215    if (!isfinite(shape.sxy)) return false;
    248216
    249     // XXX PAR[PM_PAR_SKY]  = moments->Sky;
    250217    PAR[PM_PAR_SKY]  = 0.0;
    251218    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    321288
    322289    psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7]));
     290    psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
    323291
    324292    psF64 radius = sigma * sqrt (2.0 * z);
     293    psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma);
    325294
    326295    if (isnan(radius))
    … …  
    429398bool PM_MODEL_FIT_STATUS (pmModel *model)
    430399{
    431 
    432     psF32 dP;
    433400    bool  status;
    434401
    … …  
    436403    psF32 *dPAR = model->dparams->data.F32;
    437404
    438     dP = 0;
    439     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    440     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    441     dP = sqrt (dP);
    442 
    443405    status = true;
    444 //    status &= (dP < 0.5);
    445406    status &= (PAR[PM_PAR_I0] > 0);
    446407    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    447408
    448     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]));
    450 
    451409    return status;
     410}
     411
     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;
    452440}
    453441
    … …  
    460448# undef PM_MODEL_PARAMS_FROM_PSF
    461449# undef PM_MODEL_FIT_STATUS
     450# undef PM_MODEL_SET_LIMITS

  • branches/eam_branches/20090820/psModules/src/objects/pmGrowthCurveGenerate.c

    r23989 r25766  
    6565
    6666            // use the center of the center pixel of the image
     67            // 0.5 PIX: is this offset needed? probably -- the psf model uses 0.5,0.5 as the center, double check
    6768            float xc = (int)(ix*readout->image->numCols + 0.5*readout->image->numCols) + readout->image->col0 + 0.5;
    6869            float yc = (int)(iy*readout->image->numRows + 0.5*readout->image->numRows) + readout->image->row0 + 0.5;
    … …  
    195196            return NULL;
    196197        }
    197         psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     198        psImageMaskPixels (mask, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    198199
    199200        // the 'ignore' mode is for testing

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

    r23187 r25766  
    5656
    5757    tmp->type = type;
    58     tmp->chisq = 0.0;
    59     tmp->chisqNorm = 0.0;
     58    tmp->mag = NAN;
     59    tmp->chisq = NAN;
     60    tmp->chisqNorm = NAN;
    6061    tmp->nDOF  = 0;
    6162    tmp->nPix  = 0;
    6263    tmp->nIter = 0;
    63     tmp->radiusFit = 0;
     64    tmp->fitRadius = 0;
    6465    tmp->flags = PM_MODEL_STATUS_NONE;
    6566    tmp->residuals = NULL;              // XXX should the model own this memory?
    … …  
    8687    tmp->modelParamsFromPSF = class->modelParamsFromPSF;
    8788    tmp->modelFitStatus     = class->modelFitStatus;
     89    tmp->modelSetLimits     = class->modelSetLimits;
    8890
    8991    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
    … …  
    108110    new->nIter     = model->nIter;
    109111    new->flags     = model->flags;
    110     new->radiusFit = model->radiusFit;
     112    new->fitRadius = model->fitRadius;
    111113
    112114    for (int i = 0; i < new->params->n; i++) {
    … …  
    189191    psVector *params = model->params;
    190192
    191     psS32 imageCol;
    192     psS32 imageRow;
    193     psF32 pixelValue;
     193    float imageCol;
     194    float imageRow;
     195    float pixelValue;
    194196
    195197    // save original values; restore before returning
    … …  
    232234    psF32 **Rx = NULL;
    233235    psF32 **Ry = NULL;
    234     psImageMaskType **Rm = NULL;
     236    pmResidMaskType **Rm = NULL;
    235237
    236238    if (model->residuals) {
    … …  
    240242        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
    241243        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;
     244        Rm = (model->residuals->mask) ? model->residuals->mask->data.PM_TYPE_RESID_MASK_DATA : NULL;
    243245        if (Ro) {
    244246            NX = model->residuals->Ro->numCols;
    … …  
    255257                continue;
    256258
    257             // XXX should we use using 0.5 pixel offset?
    258             // Convert to coordinate in parent image, with offset (dx,dy)
    259             imageCol = ix + image->col0 - dx;
    260             imageRow = iy + image->row0 - dy;
    261 
    262             x->data.F32[0] = (float) imageCol;
    263             x->data.F32[1] = (float) imageRow;
     259            // Convert to coordinate in parent image, with offset (dx,dy)
     260            // 0.5 PIX: the model take pixel coordinates so convert the pixel index here
     261            imageCol = ix + 0.5 + image->col0 - dx;
     262            imageRow = iy + 0.5 + image->row0 - dy;
     263
     264            x->data.F32[0] = imageCol;
     265            x->data.F32[1] = imageRow;
    264266
    265267            pixelValue = 0.0;
    … …  
    276278                float rx = xBin*ix + DX;
    277279
    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);
     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);
    327329            }
    328330
    329         skip:
     331        skip:
    330332            // add or subtract the value
    331333            if (add) {

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

    r21516 r25766  
    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
    … …  
    96107    int nIter;                          ///< number of iterations to reach min
    97108    pmModelStatus flags;                ///< model status flags
    98     float radiusFit;                    ///< fit radius actually used
     109    float fitRadius;                    ///< fit radius actually used
    99110    pmResiduals *residuals;             ///< normalized PSF residuals
    100111
    … …  
    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/20090820/psModules/src/objects/pmModelClass.c

    r20937 r25766  
    3636#include "pmErrorCodes.h"
    3737
    38 // XXX shouldn't these be defined for us in pslib.h ???
     38// XXX shouldn't these be defined for us in math.h ???
    3939double hypot(double x, double y);
    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/20090820/psModules/src/objects/pmModelClass.h

    r15697 r25766  
    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/20090820/psModules/src/objects/pmModelUtils.c

    r19999 r25766  
    4848        return NULL;
    4949    }
    50     // XXX note that model->residuals is just a reference
     50    // note that model->residuals is just a reference
    5151    modelPSF->residuals = psf->residuals;
    5252
    … …  
    6565    // set model parameters for this source based on PSF information
    6666    if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) {
    67         // XXX we do not want to raise an error here, just note that we failed
    68         // psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
    6967        psFree(modelPSF);
    7068        return NULL;
    7169    }
    7270
    73     // XXX note that model->residuals is just a reference
     71    // note that model->residuals is just a reference
    7472    modelPSF->residuals = psf->residuals;
    7573

  • branches/eam_branches/20090820/psModules/src/objects/pmPSF.h

    r24206 r25766  
    3838 *
    3939 */
    40 typedef struct
    41 {
     40typedef struct {
    4241    pmModelType type;                   ///< PSF Model in use
    4342    psArray *params;                    ///< Model parameters (psPolynomial2D)
    … …  
    6564    pmGrowthCurve *growth;              ///< apMag vs Radius
    6665    pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
    67 }
    68 pmPSF;
     66} pmPSF;
    6967
    7068typedef struct {
    … …  
    8179    bool          poissonErrorsPhotLin; ///< use poission errors for linear model fitting
    8280    bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
    83     float         radius;
     81    float         fitRadius;
     82    float         apRadius;
    8483    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
    8584} pmPSFOptions;

  • branches/eam_branches/20090820/psModules/src/objects/pmPSFtry.c

    r24206 r25766  
    3737#include "pmSourceVisual.h"
    3838
    39 bool printTrendMap (pmTrend2D *trend) {
    40 
    41     if (!trend->map) return false;
    42     if (!trend->map->map) return false;
    43 
    44     for (int j = 0; j < trend->map->map->numRows; j++) {
    45         for (int i = 0; i < trend->map->map->numCols; i++) {
    46             fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
    47         }
    48         fprintf (stderr, "\t\t\t");
    49         for (int i = 0; i < trend->map->map->numCols; i++) {
    50             fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
    51         }
    52         fprintf (stderr, "\n");
    53     }
    54     return true;
    55 }
    56 
    57 bool psImageMapCleanup (psImageMap *map) {
    58 
    59     if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;
    60 
    61     // find the weighted average of all pixels
    62     float Sum = 0.0;
    63     float Wt = 0.0;
    64     for (int j = 0; j < map->map->numRows; j++) {
    65         for (int i = 0; i < map->map->numCols; i++) {
    66             if (!isfinite(map->map->data.F32[j][i])) continue;
    67             Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
    68             Wt += map->error->data.F32[j][i];
    69         }
    70     }
    71 
    72     float Mean = Sum / Wt;
    73 
    74     // do any of the pixels in the map need to be repaired?
    75     // XXX for now, we are just replacing bad pixels with the Mean
    76     for (int j = 0; j < map->map->numRows; j++) {
    77         for (int i = 0; i < map->map->numCols; i++) {
    78             if (isfinite(map->map->data.F32[j][i]) &&
    79                 (map->error->data.F32[j][i] > 0.0)) continue;
    80             map->map->data.F32[j][i] = Mean;
    81         }
    82     }
    83     return true;
    84 }
    85 
    8639// ********  pmPSFtry functions  **************************************************
    8740// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
    … …  
    11063    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    11164
    112     test->psf       = pmPSFAlloc (options);
     65    test->psf       = NULL;
    11366    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
    11467    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    13689}
    13790
     91float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) {
    13892
    139 // build a pmPSFtry for the given model:
    140 // - fit each source with the free-floating model
    141 // - construct the pmPSF from the collection of models
    142 // - fit each source with the PSF-parameter models
    143 // - measure the pmPSF quality metric (dApResid)
     93    psAssert(residuals, "residuals cannot be NULL");
     94    psAssert(errors, "errors cannot be NULL");
     95    psAssert(residuals->n == errors->n, "residuals and errors must be the same length");
    14496
    145 // sources used in for pmPSFtry may be masked by the analysis
    146 // mask values indicate the reason the source was rejected:
     97    // given a vector of residuals and their formal errors, calculated the necessary systematic
     98    // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction
     99    // highest chi-square contributors (allowed outliers)
    147100
    148 // generate a pmPSFtry with a copy of the test PSF sources
    149 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
    150 {
    151     bool status;
    152     int Next = 0;
    153     int Npsf = 0;
     101    psVector *mask  = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK);
     102    psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32);
    154103
    155     // validate the requested model name
    156     options->type = pmModelClassGetType (modelName);
    157     if (options->type == -1) {
    158         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
    159         return NULL;
     104    // calculate the chisq vector:
     105    int Ngood = 0;
     106    for (int i = 0; i < residuals->n; i++) {
     107        chisq->data.F32[i] = PS_MAX_F32;
     108        if (!isfinite(residuals->data.F32[i])) continue;
     109        if (!isfinite(errors->data.F32[i])) continue;
     110        if (errors->data.F32[i] <= 0.0) continue;
     111        chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]);
     112        Ngood ++;
    160113    }
    161114
    162     pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
    163     if (psfTry == NULL) {
    164         psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    165         return NULL;
     115    psVector *index = psVectorSortIndex(NULL, chisq);
     116
     117    // toss out the clipFraction highest chisq values
     118    for (int i = 0; i < residuals->n; i++) {
     119        int n = index->data.S32[i];
     120        if (i < (1.0 - clipFraction)*Ngood) {
     121            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0;
     122        } else {
     123            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1;
     124        }
    166125    }
    167126
    168     // maskVal is used to test for rejected pixels, and must include markVal
    169     maskVal |= markVal;
     127    // Ndof ~= Ngood
     128    // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof
     129    // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0
     130   
     131    // use Newton-Raphson to solve for S2:
    170132
    171     // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    172     psTimerStart ("psf.fit");
    173     for (int i = 0; i < psfTry->sources->n; i++) {
     133    // use median sigma to calculate the initial guess for S2:
     134    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     135    psVectorStats (stats, errors, NULL, mask, 1);
     136    float errorMedian = stats->sampleMedian;
     137   
     138    float nPts = 0.0;
     139    float res2mean = 0.0;
     140    float ChiSq = 0.0;
     141    for (int i = 0; i < residuals->n; i++) {
     142        int n = index->data.S32[i];
     143        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
     144        res2mean += PS_SQR(residuals->data.F32[n]);
     145        ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]);
     146        nPts += 1.0;
     147    }
     148    res2mean /= nPts;
     149    ChiSq /= nPts;
     150   
     151    float S2guess = res2mean - PS_SQR(errorMedian);
    174152
    175         pmSource *source = psfTry->sources->data[i];
    176         if (!source->moments) {
    177             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    178             continue;
    179         }
    180         if (!source->moments->nPixels) {
    181             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    182             continue;
    183         }
     153    psLogMsg ("psModules", 10, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n",
     154              ChiSq, residuals->n, Ngood, nPts, S2guess);
    184155
    185         source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
    186         if (source->modelEXT == NULL) {
    187             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    188             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    189             continue;
    190         }
     156    for (int iter = 0; iter < 10; iter++) {
    191157
    192         // set object mask to define valid pixels
    193         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
     158        ChiSq = 0.0;
     159        float dRdS = 0.0;
     160        for (int i = 0; i < residuals->n; i++) {
     161            int n = index->data.S32[i];
     162            if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
     163            float error2 = PS_SQR(errors->data.F32[n]) + S2guess;
     164            ChiSq += PS_SQR(residuals->data.F32[n]) / error2;
     165            dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2);
     166        }
     167        ChiSq /= nPts;
     168        dRdS /= nPts;
    194169
    195         // fit model as EXT, not PSF
    196         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);
     170        // Note the sign on dS: dRdS above is -1 * dR/dS formally
     171        float dS = (ChiSq - 1.0) / dRdS;
     172        S2guess += dS;
     173        S2guess = PS_MAX(0.0, S2guess);
    197174
    198         // clear object mask to define valid pixels
    199         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    200 
    201         // exclude the poor fits
    202         if (!status) {
    203             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    204             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    205             continue;
    206         }
    207         Next ++;
    208     }
    209     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
    210     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
    211 
    212     if (Next == 0) {
    213         psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF.");
    214         psFree(psfTry);
    215         return NULL;
     175        psLogMsg ("psModules", 10, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);
    216176    }
    217177
    218     // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    219     if (!pmPSFFromPSFtry (psfTry)) {
    220         psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    221         psFree(psfTry);
    222         return NULL;
    223     }
     178    // free local allocations
     179    psFree (mask);
     180    psFree (chisq);
     181    psFree (stats);
     182    psFree (index);
    224183
    225     // stage 3: refit with fixed shape parameters
    226     psTimerStart ("psf.fit");
    227     for (int i = 0; i < psfTry->sources->n; i++) {
    228 
    229         pmSource *source = psfTry->sources->data[i];
    230 
    231         // masked for: bad model fit, outlier in parameters
    232         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
    233             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
    234             continue;
    235         }
    236 
    237         // set shape for this model based on PSF
    238         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
    239         if (source->modelPSF == NULL) {
    240             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
    241             abort();
    242             continue;
    243         }
    244         source->modelPSF->radiusFit = options->radius;
    245 
    246         // set object mask to define valid pixels
    247         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
    248 
    249         // fit the PSF model to the source
    250         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal);
    251 
    252         // skip poor fits
    253         if (!status) {
    254             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    255             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
    256             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    257             continue;
    258         }
    259 
    260         status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
    261         if (!status || isnan(source->apMag)) {
    262             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    263             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
    264             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    265             continue;
    266         }
    267 
    268         // clear object mask to define valid pixels
    269         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    270 
    271         psfTry->fitMag->data.F32[i] = source->psfMag;
    272         psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
    273         psfTry->metricErr->data.F32[i] = source->errMag;
    274 
    275         psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
    276         Npsf ++;
    277     }
    278     psfTry->psf->nPSFstars = Npsf;
    279 
    280     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
    281     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
    282 
    283     if (Npsf == 0) {
    284         psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
    285         psFree(psfTry);
    286         return NULL;
    287     }
    288 
    289     // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    290     // this should be linear for Poisson errors and quadratic for constant sky errors
    291     psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    292     psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    293     psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK);
    294 
    295     // generate the x and y vectors, and mask missing models
    296     for (int i = 0; i < psfTry->sources->n; i++) {
    297         pmSource *source = psfTry->sources->data[i];
    298         if (source->modelPSF == NULL) {
    299             flux->data.F32[i] = 0.0;
    300             chisq->data.F32[i] = 0.0;
    301             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    302         } else {
    303             flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
    304             chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
    305             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    306         }
    307     }
    308 
    309     // use 3hi/3lo sigma clipping on the chisq fit
    310     psStats *stats = options->stats;
    311 
    312     // linear clipped fit of chisq trend vs flux
    313     if (options->chiFluxTrend) {
    314         bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
    315                                                   mask, 0xff, chisq, NULL, flux);
    316         psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    317         psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    318 
    319         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
    320                   psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
    321 
    322         psFree(flux);
    323         psFree(mask);
    324         psFree(chisq);
    325 
    326         if (!result) {
    327             psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
    328             psFree(psfTry);
    329             return NULL;
    330         }
    331     }
    332 
    333     for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
    334         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
    335                   psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
    336                   psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
    337     }
    338 
    339     // XXX this function wants aperture radius for pmSourcePhotometry
    340     if (!pmPSFtryMetric (psfTry, options)) {
    341         psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
    342         psFree (psfTry);
    343         return NULL;
    344     }
    345 
    346     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
    347               modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    348 
    349     return (psfTry);
     184    return (sqrt(S2guess));
    350185}
    351186
    352 bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
    353 {
    354     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    355     PS_ASSERT_PTR_NON_NULL(options, false);
    356     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    357 
    358     float RADIUS = options->radius;
    359 
    360     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    361     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    362     //     dA = dAo + dsky/flux
    363     //   where flux is the flux of the star
    364     // we fit this trend to find the infinite flux aperture correction (dAo),
    365     //   the nominal sky bias (dsky), and the error on dAo
    366     // the values of dA are contaminated by stars with close neighbors in the aperture
    367     //   we use an outlier rejection to avoid this bias
    368 
    369     // r2rflux = radius^2 * ten(0.4*fitMag);
    370     psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    371 
    372     for (int i = 0; i < psfTry->sources->n; i++) {
    373         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
    374             continue;
    375         r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
    376     }
    377 
    378     // XXX test dump of aperture residual data
    379     if (psTraceGetLevel("psModules.objects") >= 5) {
    380         FILE *f = fopen ("apresid.dat", "w");
    381         for (int i = 0; i < psfTry->sources->n; i++) {
    382             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    383 
    384             pmSource *source = psfTry->sources->data[i];
    385             float x = source->peak->x;
    386             float y = source->peak->y;
    387 
    388             fprintf (f, "%d  %d, %f %f %f  %f %f %f \n",
    389                      i, keep, x, y,
    390                      psfTry->fitMag->data.F32[i],
    391                      r2rflux->data.F32[i],
    392                      psfTry->metric->data.F32[i],
    393                      psfTry->metricErr->data.F32[i]);
    394         }
    395         fclose (f);
    396     }
    397 
    398     // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
    399     // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
    400     // function of magnitude.  We re-measure the apResid statistics later in psphot using the
    401     // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
    402 
    403     // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    404     // linear clipped fit of ApResid to r2rflux
    405     psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    406     poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)
    407 
    408     // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    409     bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
    410                                               psfTry->metric, psfTry->metricErr, r2rflux);
    411     if (!result) {
    412         psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");
    413 
    414         psFree(poly);
    415         psFree(r2rflux);
    416 
    417         return false;
    418     }
    419     psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    420     psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
    421               psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
    422 
    423 
    424 
    425     // XXX test dump of fitted model (dump when tracing?)
    426     if (psTraceGetLevel("psModules.objects") >= 4) {
    427         FILE *f = fopen ("resid.dat", "w");
    428         psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux);
    429         for (int i = 0; i < psfTry->sources->n; i++) {
    430             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    431 
    432             pmSource *source = psfTry->sources->data[i];
    433             float x = source->peak->x;
    434             float y = source->peak->y;
    435 
    436             fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
    437                      i, keep, x, y,
    438                      psfTry->fitMag->data.F32[i],
    439                      r2rflux->data.F32[i],
    440                      psfTry->metric->data.F32[i],
    441                      psfTry->metricErr->data.F32[i],
    442                      apfit->data.F32[i]);
    443         }
    444         fclose (f);
    445         psFree (apfit);
    446     }
    447 
    448     // XXX drop the skyBias value, or include above??
    449     psfTry->psf->skyBias  = poly->coeff[1];
    450     psfTry->psf->ApResid  = poly->coeff[0];
    451     psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);
    452 
    453     psFree (r2rflux);
    454     psFree (poly);
    455 
    456     return true;
    457 }
    458 
    459 /*
    460   (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
    461   (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
    462   (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    463 */
    464 
    465 /*****************************************************************************
    466 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
    467 source->modelEXT entries.  The PSF ignores the first 4 (independent) model
    468 parameters and constructs a polynomial fit to the remaining as a function of
    469 image coordinate.
    470     input: psfTry with fitted source->modelEXT collection, pre-allocated psf
    471 Note: some of the array entries may be NULL (failed fits); ignore them.
    472  *****************************************************************************/
    473 bool pmPSFFromPSFtry (pmPSFtry *psfTry)
    474 {
    475     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    476     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    477 
    478     pmPSF *psf = psfTry->psf;
    479 
    480     // construct the fit vectors from the collection of objects
    481     psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    482     psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    483     psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    484 
    485     // construct the x,y terms
    486     for (int i = 0; i < psfTry->sources->n; i++) {
    487         pmSource *source = psfTry->sources->data[i];
    488         if (source->modelEXT == NULL)
    489             continue;
    490 
    491         x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    492         y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
    493     }
    494 
    495     if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
    496         psFree(x);
    497         psFree(y);
    498         psFree(z);
    499         return false;
    500     }
    501 
    502     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
    503     // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
    504     for (int i = 0; i < psf->params->n; i++) {
    505         switch (i) {
    506           case PM_PAR_SKY:
    507           case PM_PAR_I0:
    508           case PM_PAR_XPOS:
    509           case PM_PAR_YPOS:
    510           case PM_PAR_SXX:
    511           case PM_PAR_SYY:
    512           case PM_PAR_SXY:
    513             continue;
    514           default:
    515             break;
    516         }
    517 
    518         // select the per-object fitted data for this parameter
    519         for (int j = 0; j < psfTry->sources->n; j++) {
    520             pmSource *source = psfTry->sources->data[j];
    521             if (source->modelEXT == NULL) continue;
    522             z->data.F32[j] = source->modelEXT->params->data.F32[i];
    523         }
    524 
    525         psImageBinning *binning = psImageBinningAlloc();
    526         binning->nXruff = psf->trendNx;
    527         binning->nYruff = psf->trendNy;
    528         binning->nXfine = psf->fieldNx;
    529         binning->nYfine = psf->fieldNy;
    530 
    531         if (psf->psfTrendMode == PM_TREND_MAP) {
    532             psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    533             psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    534         }
    535 
    536         // free existing trend, re-alloc
    537         psFree (psf->params->data[i]);
    538         psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
    539         psFree (binning);
    540 
    541         // fit the collection of measured parameters to the PSF 2D model
    542         // the mask is carried from previous steps and updated with this operation
    543         // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
    544         if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
    545             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    546             psFree(x);
    547             psFree(y);
    548             psFree(z);
    549             return false;
    550         }
    551     }
    552 
    553     // test dump of star parameters vs position (compare with fitted values)
    554     if (psTraceGetLevel("psModules.objects") >= 4) {
    555         FILE *f = fopen ("params.dat", "w");
    556 
    557         for (int j = 0; j < psfTry->sources->n; j++) {
    558             pmSource *source = psfTry->sources->data[j];
    559             if (source == NULL) continue;
    560             if (source->modelEXT == NULL) continue;
    561 
    562             pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf);
    563 
    564             fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]);
    565 
    566             for (int i = 0; i < psf->params->n; i++) {
    567                 if (psf->params->data[i] == NULL) continue;
    568                 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
    569             }
    570             fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter);
    571 
    572             psFree(modelPSF);
    573         }
    574         fclose (f);
    575     }
    576 
    577     psFree (x);
    578     psFree (y);
    579     psFree (z);
    580     return true;
    581 }
    582 
    583 
    584 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
    585 
    586     // we are doing a robust fit.  after each pass, we drop points which are more deviant than
    587     // three sigma.  mask is currently updated for each pass.
    588 
    589     // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    590     // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
    591     // each source and fit this set of parameters.  These values are less tightly coupled, but
    592     // are still inter-related.  The fitted values do a good job of constraining the major axis
    593     // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
    594     // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
    595     // parameters, with the constraint that the minor axis must be greater than a minimum
    596     // threshold.
    597 
    598     // convert the measured source shape paramters to polarization terms
    599     psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
    600     psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
    601     psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
    602     psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);
    603 
    604     for (int i = 0; i < sources->n; i++) {
    605         // skip any masked sources (failed to fit one of the model steps or get a magnitude)
    606         if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    607 
    608         pmSource *source = sources->data[i];
    609         assert (source->modelEXT); // all unmasked sources should have modelEXT
    610 
    611         psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
    612 
    613         e0->data.F32[i] = pol.e0;
    614         e1->data.F32[i] = pol.e1;
    615         e2->data.F32[i] = pol.e2;
    616 
    617         float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
    618         mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
    619     }
    620 
    621     if (psf->psfTrendMode == PM_TREND_MAP) {
    622         float scatterTotal = 0.0;
    623         float scatterTotalMin = FLT_MAX;
    624         int entryMin = -1;
    625 
    626         psVector *dz = NULL;
    627         psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);
    628 
    629         // check the fit residuals and increase Nx,Ny until the error is minimized
    630         // pmPSFParamTrend increases the number along the longer of x or y
    631         for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    632 
    633             // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    634             for (int i = 0; i < mask->n; i++) {
    635                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    636             }
    637             if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    638                 break;
    639             }
    640 
    641             // store the resulting scatterTotal values and the scales, redo the best
    642             if (scatterTotal < scatterTotalMin) {
    643                 scatterTotalMin = scatterTotal;
    644                 entryMin = i;
    645             }
    646         }
    647         if (entryMin == -1) {
    648             psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params");
    649             return false;
    650         }
    651 
    652         // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    653         for (int i = 0; i < mask->n; i++) {
    654             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    655         }
    656         if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    657             psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    658         }
    659 
    660         pmTrend2D *trend = psf->params->data[PM_PAR_E0];
    661         psf->trendNx = trend->map->map->numCols;
    662         psf->trendNy = trend->map->map->numRows;
    663 
    664         // copy mask back to srcMask
    665         for (int i = 0; i < mask->n; i++) {
    666             srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    667         }
    668 
    669         psFree (mask);
    670         psFree (dz);
    671     } else {
    672 
    673         // XXX iterate Nx, Ny based on scatter?
    674         // XXX we force the x & y order to be the same
    675         // modify the order to correspond to the actual number of matched stars:
    676         int order = PS_MAX (psf->trendNx, psf->trendNy);
    677         if ((sources->n < 15) && (order >= 3)) order = 2;
    678         if ((sources->n < 11) && (order >= 2)) order = 1;
    679         if ((sources->n <  8) && (order >= 1)) order = 0;
    680         if ((sources->n <  3)) {
    681             psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    682             return false;
    683         }
    684         psf->trendNx = order;
    685         psf->trendNy = order;
    686 
    687         // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    688         // This way, the parameters masked by one of the fits will be applied to the others
    689         for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    690 
    691             psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    692             psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    693 
    694             pmTrend2D *trend = NULL;
    695             float mean, stdev;
    696 
    697             // XXX we are using the same stats structure on each pass: do we need to re-init it?
    698             bool status = true;
    699 
    700             trend = psf->params->data[PM_PAR_E0];
    701             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
    702             mean = psStatsGetValue (trend->stats, meanOption);
    703             stdev = psStatsGetValue (trend->stats, stdevOption);
    704             psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    705             pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    706 
    707             trend = psf->params->data[PM_PAR_E1];
    708             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
    709             mean = psStatsGetValue (trend->stats, meanOption);
    710             stdev = psStatsGetValue (trend->stats, stdevOption);
    711             psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    712             pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    713 
    714             trend = psf->params->data[PM_PAR_E2];
    715             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
    716             mean = psStatsGetValue (trend->stats, meanOption);
    717             stdev = psStatsGetValue (trend->stats, stdevOption);
    718             psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    719             pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    720 
    721             if (!status) {
    722                 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    723                 return false;
    724             }
    725         }
    726     }
    727 
    728     // test dump of the psf parameters
    729     if (psTraceGetLevel("psModules.objects") >= 4) {
    730         FILE *f = fopen ("pol.dat", "w");
    731         fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    732         for (int i = 0; i < e0->n; i++) {
    733             fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    734                      x->data.F32[i], y->data.F32[i],
    735                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
    736                      pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
    737                      pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
    738                      pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
    739                      srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    740         }
    741         fclose (f);
    742     }
    743 
    744     psFree (e0);
    745     psFree (e1);
    746     psFree (e2);
    747     psFree (mag);
    748     return true;
    749 }
    750 
    751 // fit the shape variations as a psImageMap for the given scale factor
    752 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
    753 
    754     int Nx, Ny;
    755 
    756     // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
    757     // that we have a single cell
    758     if (psf->fieldNx > psf->fieldNy) {
    759         Nx = scale;
    760         float AR = psf->fieldNy / (float) psf->fieldNx;
    761         Ny = (int) (Nx * AR + 0.5);
    762         Ny = PS_MAX (1, Ny);
    763     } else {
    764         Ny = scale;
    765         float AR = psf->fieldNx / (float) psf->fieldNy;
    766         Nx = (int) (Ny * AR + 0.5);
    767         Nx = PS_MAX (1, Nx);
    768     }
    769 
    770     // do we have enough sources for this fine of a grid?
    771     if (x->n < 10*Nx*Ny) {
    772         return false;
    773     }
    774 
    775     // XXX check this against the exising type
    776     pmTrend2DMode psfTrendMode = PM_TREND_MAP;
    777 
    778     psImageBinning *binning = psImageBinningAlloc();
    779     binning->nXruff = Nx;
    780     binning->nYruff = Ny;
    781     binning->nXfine = psf->fieldNx;
    782     binning->nYfine = psf->fieldNy;
    783     psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    784     psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    785 
    786     psFree (psf->params->data[PM_PAR_E0]);
    787     psFree (psf->params->data[PM_PAR_E1]);
    788     psFree (psf->params->data[PM_PAR_E2]);
    789 
    790     int nIter = psf->psfTrendStats->clipIter;
    791     psf->psfTrendStats->clipIter = 1;
    792     psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    793     psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    794     psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    795     psFree (binning);
    796 
    797     // if the map is 1x1 (a single value), we measure the resulting ensemble scatter
    798 
    799     // if the map is more finely sampled, divide the values into two sets: measure the fit from
    800     // one set and the scatter from the other set.
    801     psVector *x_fit = NULL;
    802     psVector *y_fit = NULL;
    803     psVector *x_tst = NULL;
    804     psVector *y_tst = NULL;
    805 
    806     psVector *e0obs_fit = NULL;
    807     psVector *e1obs_fit = NULL;
    808     psVector *e2obs_fit = NULL;
    809     psVector *e0obs_tst = NULL;
    810     psVector *e1obs_tst = NULL;
    811     psVector *e2obs_tst = NULL;
    812 
    813     if (scale == 1) {
    814         x_fit  = psMemIncrRefCounter (x);
    815         y_fit  = psMemIncrRefCounter (y);
    816         x_tst  = psMemIncrRefCounter (x);
    817         y_tst  = psMemIncrRefCounter (y);
    818         e0obs_fit = psMemIncrRefCounter (e0obs);
    819         e1obs_fit = psMemIncrRefCounter (e1obs);
    820         e2obs_fit = psMemIncrRefCounter (e2obs);
    821         e0obs_tst = psMemIncrRefCounter (e0obs);
    822         e1obs_tst = psMemIncrRefCounter (e1obs);
    823         e2obs_tst = psMemIncrRefCounter (e2obs);
    824     } else {
    825         x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    826         y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    827         x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    828         y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    829         e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    830         e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    831         e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    832         e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    833         e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    834         e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    835         for (int i = 0; i < e0obs_fit->n; i++) {
    836             // e0obs->n ==  8 or 9:
    837             // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
    838             // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
    839             x_fit->data.F32[i] = x->data.F32[2*i+0];
    840             x_tst->data.F32[i] = x->data.F32[2*i+1];
    841             y_fit->data.F32[i] = y->data.F32[2*i+0];
    842             y_tst->data.F32[i] = y->data.F32[2*i+1];
    843 
    844             e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];
    845             e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];
    846             e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
    847             e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
    848             e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
    849             e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
    850         }
    851     }
    852 
    853     // the mask marks the values not used to calculate the ApTrend
    854     psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK);
    855     // copy mask values to fitMask as a starting point
    856     for (int i = 0; i < fitMask->n; i++) {
    857         fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    858     }
    859 
    860     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    861     // This way, the parameters masked by one of the fits will be applied to the others
    862     for (int i = 0; i < nIter; i++) {
    863         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    864         psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    865         psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    866 
    867         pmTrend2D *trend = NULL;
    868         float mean, stdev;
    869 
    870         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    871         bool status = true;
    872 
    873         trend = psf->params->data[PM_PAR_E0];
    874         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    875         mean = psStatsGetValue (trend->stats, meanOption);
    876         stdev = psStatsGetValue (trend->stats, stdevOption);
    877         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    878         // printTrendMap (trend);
    879         psImageMapCleanup (trend->map);
    880         // printTrendMap (trend);
    881         pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    882 
    883         trend = psf->params->data[PM_PAR_E1];
    884         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    885         mean = psStatsGetValue (trend->stats, meanOption);
    886         stdev = psStatsGetValue (trend->stats, stdevOption);
    887         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    888         // printTrendMap (trend);
    889         psImageMapCleanup (trend->map);
    890         // printTrendMap (trend);
    891         pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    892 
    893         trend = psf->params->data[PM_PAR_E2];
    894         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    895         mean = psStatsGetValue (trend->stats, meanOption);
    896         stdev = psStatsGetValue (trend->stats, stdevOption);
    897         psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    898         // printTrendMap (trend);
    899         psImageMapCleanup (trend->map);
    900         // printTrendMap (trend);
    901         pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    902     }
    903     psf->psfTrendStats->clipIter = nIter; // restore default setting
    904 
    905     // construct the fitted values and the residuals
    906     psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
    907     psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
    908     psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);
    909 
    910     psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
    911     psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
    912     psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);
    913 
    914     // measure scatter for the unfitted points
    915     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
    916     // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
    917     pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
    918     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
    919     // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
    920 
    921     psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
    922     psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);
    923 
    924     psFree (x_fit);
    925     psFree (y_fit);
    926     psFree (x_tst);
    927     psFree (y_tst);
    928 
    929     psFree (e0obs_fit);
    930     psFree (e1obs_fit);
    931     psFree (e2obs_fit);
    932     psFree (e0obs_tst);
    933     psFree (e1obs_tst);
    934     psFree (e2obs_tst);
    935 
    936     psFree (e0fit);
    937     psFree (e1fit);
    938     psFree (e2fit);
    939 
    940     psFree (e0res);
    941     psFree (e1res);
    942     psFree (e2res);
    943 
    944     // XXX copy fitMask values back to mask
    945     for (int i = 0; i < fitMask->n; i++) {
    946         mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    947     }
    948     psFree (fitMask);
    949 
    950     return true;
    951 }
    952 
    953 // calculate the scatter of the parameters
    954 bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt)
    955 {
    956 
    957     // psStats *stats = psStatsAlloc(stdevOpt);
    958     psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);
    959 
    960     // calculate the root-mean-square of E0, E1, E2
    961     float dEsquare = 0.0;
    962     psStatsInit (stats);
    963     if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
    964         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    965         return false;
    966     }
    967     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    968 
    969     psStatsInit (stats);
    970     if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
    971         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    972         return false;
    973     }
    974     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    975 
    976     psStatsInit (stats);
    977     if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
    978         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    979         return false;
    980     }
    981     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    982 
    983     *scatterTotal = sqrtf(dEsquare);
    984 
    985     psFree(stats);
    986     return true;
    987 }
    988 
    989 // calculate the minimum scatter of the parameters
    990 bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
    991                             psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
    992 {
    993 
    994     psStats *statsS = psStatsAlloc(stdevOpt);
    995 
    996     // measure the trend in bins with 10 values each; if < 10 total, use them all
    997     int nBin = PS_MAX (mag->n / nGroup, 1);
    998 
    999     // use mag to group parameters in sequence
    1000     psVector *index = psVectorSortIndex (NULL, mag);
    1001 
    1002     // subset vectors for mag and dap values within the given range
    1003     psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1004     psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1005     psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1006     psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
    1007 
    1008     int n = 0;
    1009     float min = INFINITY;               // Minimum error
    1010     for (int i = 0; i < nBin; i++) {
    1011         int j;
    1012         int nValid = 0;
    1013         for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    1014             int N = index->data.U32[n];
    1015             dE0subset->data.F32[j] = e0res->data.F32[N];
    1016             dE1subset->data.F32[j] = e1res->data.F32[N];
    1017             dE2subset->data.F32[j] = e2res->data.F32[N];
    1018 
    1019             mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    1020             if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
    1021         }
    1022         if (nValid < 3) continue;
    1023 
    1024         dE0subset->n = j;
    1025         dE1subset->n = j;
    1026         dE2subset->n = j;
    1027         mkSubset->n = j;
    1028 
    1029         // calculate the root-mean-square of E0, E1, E2
    1030         float dEsquare = 0.0;
    1031         psStatsInit (statsS);
    1032         if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
    1033         }
    1034         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1035 
    1036         psStatsInit (statsS);
    1037         if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
    1038             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1039             return false;
    1040         }
    1041         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1042 
    1043         psStatsInit (statsS);
    1044         if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
    1045             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1046             return false;
    1047         }
    1048         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1049 
    1050         if (isfinite(dEsquare)) {
    1051             float err = sqrtf(dEsquare);
    1052             if (err < min) {
    1053                 min = err;
    1054             }
    1055         }
    1056     }
    1057     *errorFloor = min;
    1058 
    1059     psFree (dE0subset);
    1060     psFree (dE1subset);
    1061     psFree (dE2subset);
    1062     psFree (mkSubset);
    1063 
    1064     psFree(index);
    1065 
    1066     psFree(statsS);
    1067 
    1068     return true;
    1069 }

  • branches/eam_branches/20090820/psModules/src/objects/pmPSFtry.h

    r21183 r25766  
    8989 *
    9090 */
    91 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType mark);
     91pmPSFtry *pmPSFtryModel (
     92    const psArray *sources,             ///< PSF sources to use in the pmPSF model analysis
     93    const char *modelName,              ///< human-readable name of desired model
     94    pmPSFOptions *options,
     95    psImageMaskType maskVal,
     96    psImageMaskType mark
     97    );
     98
     99/** fit EXT models to all possible psf sources */
     100bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
     101
     102bool pmPSFtryMakePSF (pmPSFtry *psfTry);
     103
     104bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
    92105
    93106/** pmPSFtryMetric()
    … …  
    97110 *
    98111 */
    99 bool pmPSFtryMetric(
    100     pmPSFtry *psfTry,                  ///< Add comment.
    101     pmPSFOptions *options              ///< PSF fitting options
    102 );
     112bool pmPSFtryMetric(pmPSFtry *psfTry);
    103113
    104114/** pmPSFtryMetric_Alt()
    … …  
    112122    float RADIUS                        ///< Add comment.
    113123);
     124
     125bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
     126
     127float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
     128
     129/// @}
     130# endif
    114131
    115132/**
    … …  
    125142 *
    126143 */
    127 bool pmPSFFromPSFtry (pmPSFtry *psfTry);
    128144
    129 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
    130 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz);
    131 bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt);
    132 bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt);
    133 
    134 /// @}
    135 # endif

  • branches/eam_branches/20090820/psModules/src/objects/pmPeaks.c

    r24623 r25766  
    6060    // if min point is too deviant, use the peak value
    6161    // XXX need to calculate dx, dy correctly
     62    // 0.5 PIX: peaks are calculated using the pixel index and converted here to pixel coords
    6263    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
    63         peak->xf = min.x + ix + image->col0;
    64         peak->yf = min.y + iy + image->row0;
     64        peak->xf = min.x + ix + image->col0 + 0.5;
     65        peak->yf = min.y + iy + image->row0 + 0.5;
    6566
    6667        // These errors are fractional errors, and should be scaled by the
    … …  
    7374        peak->yf = PS_MAX (PS_MIN (peak->yf, image->numRows - 1), image->row0);
    7475    } else {
    75         peak->xf = ix;
    76         peak->yf = iy;
     76        peak->xf = ix + 0.5;
     77        peak->yf = iy + 0.5;
    7778        peak->dx = NAN;
    7879        peak->dy = NAN;

  • branches/eam_branches/20090820/psModules/src/objects/pmSource.c

    r24874 r25766  
    282282
    283283    psArray *peaks  = NULL;
    284     pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0};
    285     pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0};
     284    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0, 0, 0.0};
     285    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0, 0, 0.0};
    286286    pmPSFClump psfClump;
    287287
    … …  
    319319
    320320        // construct a sigma-plane image
    321         int numCols = SX_MAX / PSF_CLUMP_GRID_SCALE, numRows = SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
     321        int numCols = 1 + SX_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
     322        int numRows = 1 + SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
    322323        psTrace("psModules.objects", 10, "sigma-plane dimensions: %dx%d\n", numCols, numRows);
    323324        psImage *splane = psImageAlloc(numCols, numRows, PS_TYPE_F32); // sigma-plane image
    … …  
    332333            }
    333334
    334             int x = src->peak->x, y = src->peak->y; // Coordinates of peak
    335             if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
    336                 continue;
    337             }
     335            if (region) {
     336                int x = src->peak->x, y = src->peak->y; // Coordinates of peak
     337                if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
     338                    continue;
     339                }
     340            }
    338341
    339342            if (src->mode & PM_SOURCE_MODE_BLEND) {
    340343                continue;
    341344            }
     345
     346            if (!src->moments->nPixels) continue;
    342347
    343348            if (src->moments->SN < PSF_SN_LIM) {
    … …  
    384389
    385390        // find the peak in this image
    386         psStats *stats = psStatsAlloc (PS_STAT_MAX);
     391        psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV);
    387392        if (!psImageStats (stats, splane, NULL, 0)) {
    388393            psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
    … …  
    393398        peaks = pmPeaksInImage (splane, stats[0].max / 2);
    394399        psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
     400
     401        psfClump.nSigma = stats->sampleStdev;
    395402
    396403        const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
    … …  
    430437        psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);
    431438
     439        // XXX store the mean sigma?
     440        float meanSigma = psfClump.nSigma;
     441        psfClump.nStars = clump->value;
     442        psfClump.nSigma = clump->value / meanSigma;
     443
    432444        // define section window for clump
    433445        minSx = clump->x * PSF_CLUMP_GRID_SCALE - 2.0*PSF_CLUMP_GRID_SCALE;
    … …  
    452464                continue;
    453465
    454             if (tmpSrc->peak->x < region->x0) continue;
    455             if (tmpSrc->peak->x > region->x1) continue;
    456             if (tmpSrc->peak->y < region->y0) continue;
    457             if (tmpSrc->peak->y > region->y1) continue;
     466            if (region) {
     467                if (tmpSrc->peak->x < region->x0) continue;
     468                if (tmpSrc->peak->x > region->x1) continue;
     469                if (tmpSrc->peak->y < region->y0) continue;
     470                if (tmpSrc->peak->y > region->y1) continue;
     471            }
    458472
    459473            if (tmpSrc->moments->Mxx < minSx)
    … …  
    919933    if (model == NULL) return false;  // model must be defined
    920934
     935    bool addNoise = mode & PM_MODEL_OP_NOISE;
     936
    921937    if (source->modelFlux) {
    922938        // add in the pixels from the modelFlux image
    … …  
    940956
    941957        psF32 **target = source->pixels->data.F32;
    942         if (mode & PM_MODEL_OP_NOISE) {
    943             // XXX need to scale by the gain...
     958        if (addNoise) {
     959            // when adding noise, we assume the shape and Io have been modified
    944960            target = source->variance->data.F32;
    945961        }
    946962
    947         // XXX need to respect the source and model masks?
    948963        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    949964            int oy = iy + dY;
    … …  
    959974            }
    960975        }
     976        if (!addNoise) {
     977            if (add) {
     978                source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     979            } else {
     980                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     981            }
     982        }
    961983        return true;
    962984    }
    963985
    964986    psImage *target = source->pixels;
    965     if (mode & PM_MODEL_OP_NOISE) {
     987    if (addNoise) {
    966988        target = source->variance;
    967989    }
    968990
    969     if (add) {
    970         status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    971     } else {
    972         status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     991    if (!addNoise) {
     992        if (add) {
     993            status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     994            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     995        } else {
     996            status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     997            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     998        }
    973999    }
    9741000

  • branches/eam_branches/20090820/psModules/src/objects/pmSource.h

    r24579 r25766  
    3838
    3939typedef enum {
    40     PM_SOURCE_TMPF_MODEL_GUESS = 0x0001,
    41     PM_SOURCE_TMPF_SUBTRACTED  = 0x0002,
     40    PM_SOURCE_TMPF_MODEL_GUESS   = 0x0001,
     41    PM_SOURCE_TMPF_SUBTRACTED    = 0x0002,
     42    PM_SOURCE_TMPF_SIZE_MEASURED = 0x0004,
    4243} pmSourceTmpF;
    4344
    … …  
    8182    float extNsigma;                    ///< Nsigma deviation from PSF to EXT
    8283    float sky, skyErr;                  ///< The sky and its error at the center of the object
     84    float apRadius;
    8385    psRegion region;                    ///< area on image covered by selected pixels
    8486    pmSourceExtendedPars *extpars;      ///< extended source parameters
    … …  
    98100    float Y;
    99101    float dY;
     102    int nStars;
     103    float nSigma;
    100104}
    101105pmPSFClump;

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceExtendedPars.c

    r23487 r25766  
    1717#endif
    1818
    19 #include <stdio.h>
    20 #include <math.h>
    21 #include <string.h>
     19// #include <stdio.h>
     20// #include <math.h>
     21// #include <string.h>
    2222#include <pslib.h>
    23 #include "pmHDU.h"
    24 #include "pmFPA.h"
    25 #include "pmFPAMaskWeight.h"
    26 #include "pmSpan.h"
    27 #include "pmFootprint.h"
    28 #include "pmPeaks.h"
    29 #include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
    33 #include "pmPSF.h"
    34 #include "pmModel.h"
    35 #include "pmSource.h"
    36 
     23// #include "pmHDU.h"
     24// #include "pmFPA.h"
     25// #include "pmFPAMaskWeight.h"
     26// #include "pmSpan.h"
     27// #include "pmFootprint.h"
     28// #include "pmPeaks.h"
     29// #include "pmMoments.h"
     30// #include "pmResiduals.h"
     31// #include "pmGrowthCurve.h"
     32// #include "pmTrend2D.h"
     33// #include "pmPSF.h"
     34// #include "pmModel.h"
     35// #include "pmSource.h"
     36#include "pmSourceExtendedPars.h"
     37
     38// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
     39// intermediate data used to measure the profile
     40static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile)
     41{
     42    if (!profile) {
     43        return;
     44    }
     45    psFree(profile->radii);
     46    psFree(profile->fluxes);
     47    psFree(profile->theta);
     48    psFree(profile->isophotalRadii);
     49
     50    psFree(profile->radiusElliptical);
     51    psFree(profile->fluxElliptical);
     52
     53    psFree(profile->binSB);
     54    psFree(profile->binSBstdev);
     55    psFree(profile->binSBerror);
     56
     57    psFree(profile->radialBins);
     58    psFree(profile->area);
     59}
     60
     61pmSourceRadialProfile *pmSourceRadialProfileAlloc()
     62{
     63    pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile));
     64    psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree);
     65
     66    profile->radii = NULL;
     67    profile->fluxes = NULL;
     68    profile->theta = NULL;
     69    profile->isophotalRadii = NULL;
     70
     71    profile->radiusElliptical = NULL;
     72    profile->fluxElliptical = NULL;
     73
     74    profile->binSB = NULL;
     75    profile->binSBstdev = NULL;
     76    profile->binSBerror = NULL;
     77
     78    profile->radialBins = NULL;
     79    profile->area = NULL;
     80
     81    return profile;
     82}
     83
     84bool psMemCheckSourceRadialProfile(psPtr ptr)
     85{
     86    PS_ASSERT_PTR(ptr, false);
     87    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree);
     88}
     89
     90
     91// *** pmSourceRadialProfileFreeVectors frees the intermediate data values
     92bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) {
     93
     94    psFree(profile->radii);
     95    psFree(profile->fluxes);
     96    psFree(profile->theta);
     97    psFree(profile->isophotalRadii);
     98
     99    psFree(profile->radiusElliptical);
     100    psFree(profile->fluxElliptical);
     101
     102    // psFree(profile->binSB);
     103    // psFree(profile->binSBstdev);
     104    // psFree(profile->binSBerror);
     105   
     106    // psFree(profile->radialBins);
     107    psFree(profile->area);
     108
     109    profile->radii = NULL;
     110    profile->fluxes = NULL;
     111    profile->theta = NULL;
     112    profile->isophotalRadii = NULL;
     113
     114    profile->radiusElliptical = NULL;
     115    profile->fluxElliptical = NULL;
     116
     117    // profile->binSB = NULL;
     118    // profile->binSBstdev = NULL;
     119    // profile->binSBerror = NULL;
     120   
     121    // profile->radialBins = NULL;
     122    profile->area = NULL;
     123
     124    return true;
     125}
     126
     127// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
     128# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B])
     129# define SWAP_INDEX(TYPE,A,B) { \
     130  float tmp; \
     131  if (A != B) { \
     132    tmp = index->data.F32[A]; \
     133    index->data.F32[A] = index->data.F32[B]; \
     134    index->data.F32[B] = tmp; \
     135    tmp = extra->data.F32[A]; \
     136    extra->data.F32[A] = extra->data.F32[B]; \
     137    extra->data.F32[B] = tmp; \
     138  } \
     139}
     140
     141bool pmSourceRadialProfileSortPair (psVector *index, psVector *extra) {
     142
     143    // sort the vector set by the radius
     144    PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE);
     145    return true;
     146}
     147
     148// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source
    37149static void pmSourceExtendedParsFree (pmSourceExtendedPars *pars) {
    38150    if (!pars) return;
    39151
    40152    psFree(pars->profile);
    41     psFree(pars->annuli);
    42     psFree(pars->isophot);
    43     psFree(pars->petrosian);
    44     psFree(pars->kron);
     153    psFree(pars->petrosian_50);
     154    psFree(pars->petrosian_80);
    45155    return;
    46156}
    … …  
    51161
    52162    pars->profile = NULL;
    53     pars->annuli = NULL;
    54     pars->isophot = NULL;
    55     pars->petrosian = NULL;
    56     pars->kron = NULL;
     163    pars->petrosian_50 = NULL;
     164    pars->petrosian_80 = NULL;
    57165
    58166    return pars;
    … …  
    66174
    67175
    68 static void pmSourceRadialProfileFree (pmSourceRadialProfile *profile) {
    69     if (!profile) return;
    70 
    71     psFree(profile->radius);
    72     psFree(profile->flux);
    73     psFree(profile->variance);
     176// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind
     177static void pmSourceExtendedFluxFree (pmSourceExtendedFlux *flux) {
     178    if (!flux) return;
    74179    return;
    75180}
    76181
    77 pmSourceRadialProfile *pmSourceRadialProfileAlloc (void) {
    78 
    79     pmSourceRadialProfile *profile = (pmSourceRadialProfile *) psAlloc(sizeof(pmSourceRadialProfile));
    80     psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree);
    81 
    82     profile->radius = NULL;
    83     profile->flux = NULL;
    84     profile->variance = NULL;
    85 
    86     return profile;
    87 }
    88 
    89 bool psMemCheckSourceRadialProfile(psPtr ptr)
     182pmSourceExtendedFlux *pmSourceExtendedFluxAlloc (void) {
     183
     184    pmSourceExtendedFlux *flux = (pmSourceExtendedFlux *) psAlloc(sizeof(pmSourceExtendedFlux));
     185    psMemSetDeallocator(flux, (psFreeFunc) pmSourceExtendedFluxFree);
     186
     187    flux->flux = 0.0;
     188    flux->fluxErr = 0.0;
     189    flux->radius = 0.0;
     190    flux->radiusErr = 0.0;
     191
     192    return flux;
     193}
     194
     195
     196bool psMemCheckSourceExtendedFlux(psPtr ptr)
    90197{
    91198    PS_ASSERT_PTR(ptr, false);
    92     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree);
    93 }
    94 
    95 
    96 static void pmSourceIsophotalValuesFree (pmSourceIsophotalValues *isophot) {
    97     if (!isophot) return;
    98     return;
    99 }
    100 
    101 pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc (void) {
    102 
    103     pmSourceIsophotalValues *isophot = (pmSourceIsophotalValues *) psAlloc(sizeof(pmSourceIsophotalValues));
    104     psMemSetDeallocator(isophot, (psFreeFunc) pmSourceIsophotalValuesFree);
    105 
    106     isophot->mag = 0.0;
    107     isophot->magErr = 0.0;
    108     isophot->rad = 0.0;
    109     isophot->radErr = 0.0;
    110 
    111     return isophot;
    112 }
    113 
    114 
    115 bool psMemCheckSourceIsophotalValues(psPtr ptr)
    116 {
    117     PS_ASSERT_PTR(ptr, false);
    118     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceIsophotalValuesFree);
    119 }
    120 
    121 
    122 static void pmSourcePetrosianValuesFree (pmSourcePetrosianValues *petrosian) {
    123     if (!petrosian) return;
    124     return;
    125 }
    126 
    127 pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc (void) {
    128 
    129     pmSourcePetrosianValues *petrosian = (pmSourcePetrosianValues *) psAlloc(sizeof(pmSourcePetrosianValues));
    130     psMemSetDeallocator(petrosian, (psFreeFunc) pmSourcePetrosianValuesFree);
    131 
    132     petrosian->mag = 0.0;
    133     petrosian->magErr = 0.0;
    134     petrosian->rad = 0.0;
    135     petrosian->radErr = 0.0;
    136 
    137     return petrosian;
    138 }
    139 
    140 
    141 bool psMemCheckSourcePetrosianValues(psPtr ptr)
    142 {
    143     PS_ASSERT_PTR(ptr, false);
    144     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourcePetrosianValuesFree);
    145 }
    146 
    147 static void pmSourceKronValuesFree (pmSourceKronValues *kron) {
    148     if (!kron) return;
    149     return;
    150 }
    151 
    152 pmSourceKronValues *pmSourceKronValuesAlloc (void) {
    153 
    154     pmSourceKronValues *kron = (pmSourceKronValues *) psAlloc(sizeof(pmSourceKronValues));
    155     psMemSetDeallocator(kron, (psFreeFunc) pmSourceKronValuesFree);
    156 
    157     kron->mag = 0.0;
    158     kron->magErr = 0.0;
    159     kron->rad = 0.0;
    160     kron->radErr = 0.0;
    161 
    162     return kron;
    163 }
    164 
    165 
    166 bool psMemCheckSourceKronValues(psPtr ptr)
    167 {
    168     PS_ASSERT_PTR(ptr, false);
    169     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceKronValuesFree);
    170 }
    171 
    172 
    173 static void pmSourceAnnuliFree (pmSourceAnnuli *annuli) {
    174     if (!annuli) return;
    175 
    176     psFree (annuli->flux);
    177     psFree (annuli->fluxErr);
    178     psFree (annuli->fluxVar);
    179 
    180     return;
    181 }
    182 
    183 pmSourceAnnuli *pmSourceAnnuliAlloc (void) {
    184 
    185     pmSourceAnnuli *annuli = (pmSourceAnnuli *) psAlloc(sizeof(pmSourceAnnuli));
    186     psMemSetDeallocator(annuli, (psFreeFunc) pmSourceAnnuliFree);
    187 
    188     annuli->flux = NULL;
    189     annuli->fluxErr = NULL;
    190     annuli->fluxVar = NULL;
    191 
    192     return annuli;
    193 }
    194 
    195 
    196 bool psMemCheckSourceAnnuli(psPtr ptr)
    197 {
    198     PS_ASSERT_PTR(ptr, false);
    199     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceAnnuliFree);
    200 }
     199    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree);
     200}

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceExtendedPars.h

    r23487 r25766  
    1515
    1616typedef struct {
    17   psVector *radius;
    18   psVector *flux;
    19   psVector *variance;
     17    psArray  *radii;                    // radii for raw radial profiles at evenly-spaced angles
     18    psArray  *fluxes;                   // fluxes measured at above radii
     19    psVector *theta;                    // angles corresponding to above radial profiles
     20    psVector *isophotalRadii;           // isophotal radius for the above angles
     21
     22    psVector *radiusElliptical;         // normalized radial coordinates for all relevant pixels
     23    psVector *fluxElliptical;           // flux for the above radial coordinates
     24
     25    psVector *binSB;                    // mean surface brightness within radial bins
     26    psVector *binSBstdev;               // scatter of mean surface brightness within radial bins
     27    psVector *binSBerror;               // formal error on mean surface brightness within radial bins
     28
     29    psVector *radialBins;               // radii corresponding to above binnedBlux
     30    psVector *area;                     // differential area of the non-overlapping radial bins
     31
     32    psEllipseAxes axes;                 // shape of elliptical contour
    2033} pmSourceRadialProfile;
    2134
    2235typedef struct {
    23   psVector *flux;
    24   psVector *fluxErr;
    25   psVector *fluxVar;
    26 } pmSourceAnnuli;
    27 
    28 typedef struct {
    29   float mag;
    30   float magErr;
    31   float rad;
    32   float radErr;
    33 } pmSourceIsophotalValues;
    34 
    35 typedef struct {
    36   float mag;
    37   float magErr;
    38   float rad;
    39   float radErr;
    40 } pmSourcePetrosianValues;
    41 
    42 typedef struct {
    43   float mag;
    44   float magErr;
    45   float rad;
    46   float radErr;
    47 } pmSourceKronValues;
     36  float flux;
     37  float fluxErr;
     38  float radius;
     39  float radiusErr;
     40} pmSourceExtendedFlux;
    4841
    4942typedef struct {
    5043  pmSourceRadialProfile   *profile;
    51   pmSourceAnnuli          *annuli;
    52   pmSourceIsophotalValues *isophot;
    53   pmSourcePetrosianValues *petrosian;
    54   pmSourceKronValues      *kron;
     44  pmSourceExtendedFlux    *petrosian_50;
     45  pmSourceExtendedFlux    *petrosian_80;
    5546} pmSourceExtendedPars;
    5647
    57 pmSourceExtendedPars *pmSourceExtendedParsAlloc(void);
     48// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
     49// intermediate data used to measure the profile
     50pmSourceRadialProfile *pmSourceRadialProfileAlloc();
     51bool psMemCheckSourceRadialProfile(psPtr ptr);
     52
     53// *** pmSourceRadialProfileFreeVectors frees the intermediate data values
     54bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile);
     55
     56// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source
     57pmSourceExtendedPars *pmSourceExtendedParsAlloc (void);
    5858bool psMemCheckSourceExtendedPars(psPtr ptr);
    59 pmSourceRadialProfile *pmSourceRadialProfileAlloc(void);
    60 bool psMemCheckSourceRadialProfile(psPtr ptr);
    61 pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc(void);
    62 bool psMemCheckSourceIsophotalValues(psPtr ptr);
    63 pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc(void);
    64 bool psMemCheckSourcePetrosianValues(psPtr ptr);
    65 pmSourceKronValues *pmSourceKronValuesAlloc(void);
    66 bool psMemCheckSourceKronValues(psPtr ptr);
    67 pmSourceAnnuli *pmSourceAnnuliAlloc(void);
    68 bool psMemCheckSourceAnnuli(psPtr ptr);
     59
     60// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind
     61pmSourceExtendedFlux *pmSourceExtendedFluxAlloc(void);
     62bool psMemCheckSourceExtendedFlux(psPtr ptr);
     63
     64// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
     65bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra);
    6966
    7067/// @}

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceFitModel.c

    r23989 r25766  
    9696
    9797            // Convert i/j to image space:
    98             coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    99             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     98            // 0.5 PIX: the coordinate values must be in pixel coords, not index           
     99            coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0);
     100            coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
    100101            x->data[nPix] = (psPtr *) coord;
    101102            y->data.F32[nPix] = source->pixels->data.F32[i][j];
    … …  
    175176            continue;
    176177        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    177         if (psTraceGetLevel("psModules.objects") >= 4) {
    178             fprintf (stderr, "%f +/- %f\n", params->data.F32[i], dparams->data.F32[i]);
    179         }
     178        psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]);
    180179    }
    181180    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    … …  
    186185    model->nPix  = y->n;
    187186    model->nDOF  = y->n - nParams;
     187    model->chisqNorm = model->chisq / model->nDOF;
    188188    model->flags |= PM_MODEL_STATUS_FITTED;
    189189    if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceFitSet.c

    r23487 r25766  
    321321
    322322        for (int j = 0; j < model->params->n; j++, n++) {
    323             if (psTraceGetLevel("psModules.objects") >= 4) {
    324                 fprintf (stderr, "%f ", param->data.F32[n]);
    325             }
    326323            model->params->data.F32[j] = param->data.F32[n];
    327324            model->dparams->data.F32[j] = dparam->data.F32[n];
     325            psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]);
    328326        }
    329327        psTrace ("psModules.objects", 4, " src %d", i);
    … …  
    483481
    484482            // Convert i/j to image space:
    485             coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    486             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     483            // 0.5 PIX: the coordinate values must be in pixel coords, not index           
     484            coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0);
     485            coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
    487486            x->data[nPix] = (psPtr *) coord;
    488487            y->data.F32[nPix] = source->pixels->data.F32[i][j];

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO.c

    r24694 r25766  
    4242#include "pmSource.h"
    4343#include "pmModelClass.h"
     44#include "pmDetEff.h"
    4445#include "pmSourceIO.h"
    4546
    4647#define BLANK_HEADERS "BLANK.HEADERS"   // Name of metadata in camera configuration containing header names
    4748                                        // for putting values into a blank PHU
     49
     50// lookup the EXTNAME values used for table data and image header segments
     51static bool sourceExtensions(psString *headname, // Extension name for header
     52                             psString *dataname, // Extension name for data
     53                             psString *deteffname, // Extension name for detection efficiency
     54                             psString *xsrcname, // Extension name for extended sources
     55                             psString *xfitname, // Extension name for extended fits
     56                             const pmFPAfile *file, // File of interest
     57                             const pmFPAview *view // View to level of interest
     58                             )
     59{
     60    bool status;                        // Status of MD lookup
     61
     62    // Menu of EXTNAME rules
     63    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
     64    if (!menu) {
     65        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
     66        return false;
     67    }
     68
     69    // EXTNAME for image header
     70    if (headname) {
     71        const char *rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
     72        if (!rule) {
     73            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
     74            return false;
     75        }
     76        *headname = pmFPAfileNameFromRule(rule, file, view);
     77    }
     78
     79    // EXTNAME for table data
     80    if (dataname) {
     81        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
     82        if (!rule) {
     83            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
     84            return false;
     85        }
     86        *dataname = pmFPAfileNameFromRule(rule, file, view);
     87    }
     88
     89    // EXTNAME for detection efficiency
     90    if (deteffname) {
     91        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DETEFF");
     92        if (!rule) {
     93            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DETEFF in EXTNAME.RULES in camera.config");
     94            return false;
     95        }
     96        *deteffname = pmFPAfileNameFromRule(rule, file, view);
     97    }
     98
     99    // EXTNAME for extended source data table
     100    if (xsrcname) {
     101        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
     102        if (!rule) {
     103            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
     104            return false;
     105        }
     106        *xsrcname = pmFPAfileNameFromRule (rule, file, view);
     107    }
     108
     109    if (xfitname) {
     110        // EXTNAME for extended source data table
     111        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
     112        if (!rule) {
     113            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
     114            return false;
     115        }
     116        *xfitname = pmFPAfileNameFromRule (rule, file, view);
     117    }
     118
     119    return true;
     120}
     121
    48122
    49123// translations between psphot object types and dophot object types
    … …  
    271345
    272346    char *exttype  = NULL;
    273     char *dataname = NULL;
    274     char *xsrcname = NULL;
    275     char *xfitname = NULL;
    276     char *headname = NULL;
    277347
    278348    // if sources is NULL, write out an empty table
    … …  
    354424
    355425        // define the EXTNAME values for the different data segments:
    356         {
    357             // lookup the EXTNAME values used for table data and image header segments
    358             char *rule = NULL;
    359 
    360             // Menu of EXTNAME rules
    361             psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    362             if (!menu) {
    363                 psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
    364                 return false;
    365             }
    366 
    367             // EXTNAME for image header
    368             rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
    369             if (!rule) {
    370                 psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
    371                 return false;
    372             }
    373             headname = pmFPAfileNameFromRule (rule, file, view);
    374 
    375             // EXTNAME for table data
    376             rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
    377             if (!rule) {
    378                 psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
    379                 return false;
    380             }
    381             dataname = pmFPAfileNameFromRule (rule, file, view);
    382 
    383             if (XSRC_OUTPUT) {
    384               // EXTNAME for extended source data table
    385               rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
    386               if (!rule) {
    387                 psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
    388                 return false;
    389               }
    390               xsrcname = pmFPAfileNameFromRule (rule, file, view);
    391             }
    392             if (XFIT_OUTPUT) {
    393               // EXTNAME for extended source data table
    394               rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
    395               if (!rule) {
    396                 psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
    397                 return false;
    398               }
    399               xfitname = pmFPAfileNameFromRule (rule, file, view);
    400             }
     426        psString headname = NULL;
     427        psString dataname = NULL;
     428        psString deteffname = NULL;
     429        psString xsrcname = NULL;
     430        psString xfitname = NULL;
     431        if (!sourceExtensions(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL,
     432                              XFIT_OUTPUT ? &xfitname : NULL, file, view)) {
     433            return false;
    401434        }
    402435
    … …  
    480513
    481514            // XXX these are case-sensitive since the EXTYPE is case-sensitive
    482             status = false;
     515            status = true;
    483516            if (!strcmp (exttype, "SMPDATA")) {
    484                 status = pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
     517                status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
    485518            }
    486519            if (!strcmp (exttype, "PS1_DEV_0")) {
    487                 status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
     520                status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
    488521            }
    489522            if (!strcmp (exttype, "PS1_DEV_1")) {
    490                 status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
     523                status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
    491524            }
    492525            if (!strcmp (exttype, "PS1_CAL_0")) {
    493                 status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
     526                status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
    494527            }
    495528            if (!strcmp (exttype, "PS1_V1")) {
    496                 status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
     529                status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
    497530            }
    498531            if (!strcmp (exttype, "PS1_V2")) {
    499                 status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
    500             }
     532                status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
     533            }
     534
     535            if (deteffname) {
     536                status &= pmReadoutWriteDetEff(file->fits, readout, outhead, deteffname);
     537            }
     538
    501539            if (xsrcname) {
    502540              if (!strcmp (exttype, "PS1_DEV_1")) {
    503                   status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
     541                  status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
    504542              }
    505543              if (!strcmp (exttype, "PS1_CAL_0")) {
    506                   status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
     544                  status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    507545              }
    508546              if (!strcmp (exttype, "PS1_V1")) {
    509                   status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
     547                  status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
    510548              }
    511549              if (!strcmp (exttype, "PS1_V2")) {
    512                   status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
     550                  status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
    513551              }
    514552            }
    515553            if (xfitname) {
    516554              if (!strcmp (exttype, "PS1_DEV_1")) {
    517                   status = pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
     555                  status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
    518556              }
    519557              if (!strcmp (exttype, "PS1_CAL_0")) {
    520                   status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
     558                  status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
    521559              }
    522560              if (!strcmp (exttype, "PS1_V1")) {
    523                   status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
     561                  status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
    524562              }
    525563              if (!strcmp (exttype, "PS1_V2")) {
    526                   status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
     564                  status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
    527565              }
    528566            }
    … …  
    534572                psFree (xfitname);
    535573                psFree (outhead);
     574                psFree (deteffname);
    536575                return false;
    537576            }
    … …  
    545584        psFree (xfitname);
    546585        psFree (outhead);
     586        psFree (deteffname);
    547587        break;
    548588
    … …  
    572612    // not needed if only one chip
    573613    if (file->fpa->chips->n == 1) {
    574         pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
    575         return true;
     614        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
     615        return true;
    576616    }
    577617
    … …  
    885925        hdu = pmFPAviewThisHDU (view, file->fpa);
    886926
    887         // lookup the EXTNAME values used for table data and image header segments
    888         char *rule = NULL;
    889         // Menu of EXTNAME rules
    890         psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    891         if (!menu) {
    892             psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
    893             return false;
    894         }
    895         // EXTNAME for image header
    896         rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
    897         if (!rule) {
    898             psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
    899             return false;
    900         }
    901         char *headname = pmFPAfileNameFromRule (rule, file, view);
    902         // EXTNAME for table data
    903         rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
    904         if (!rule) {
    905             psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
    906             return false;
    907         }
    908         char *dataname = pmFPAfileNameFromRule (rule, file, view);
     927        // define the EXTNAME values for the different data segments:
     928        psString headname = NULL;
     929        psString dataname = NULL;
     930        psString deteffname = NULL;
     931        if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) {
     932            return false;
     933        }
    909934
    910935        // advance to the IMAGE HEADER extension
    … …  
    915940                psFree (headname);
    916941                psFree (dataname);
     942                psFree (deteffname);
    917943                return true;
    918944            }
    … …  
    958984                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
    959985            }
     986
     987            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
     988#if 0
     989                psError(PS_ERR_IO, false, "Unable to read detection efficiency");
     990                return false;
     991#else
     992                // No great loss
     993                psErrorClear();
     994#endif
     995            }
    960996        }
    961997
    … …  
    963999        psFree (headname);
    9641000        psFree (dataname);
     1001        psFree (deteffname);
    9651002        psFree (tableHeader);
    9661003        break;
    … …  
    10701107}
    10711108
    1072    
     1109

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r24401 r25766  
    128128            nDOF = model->nDOF;
    129129            nPix = model->nPix;
    130             apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
     130            apRadius = source->apRadius;
    131131            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    132132        } else {
    … …  
    308308        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    309309        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
     310        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    310311
    311312        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    313314        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
    314315        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
    315         model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    316316
    317317        source->moments = pmMomentsAlloc ();
    … …  
    353353
    354354    // which extended source analyses should we perform?
    355     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    356     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    357     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    358     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     355    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     356    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     357    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     358    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    359359
    360360    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    396396        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    397397
     398# if (0)
    398399        // Petrosian measurements
    399400        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    476477        }
    477478
     479# endif
    478480        psArrayAdd (table, 100, row);
    479481        psFree (row);

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r24403 r25766  
    128128            nDOF = model->nDOF;
    129129            nPix = model->nPix;
    130             apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
     130            apRadius = source->apRadius;
    131131            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    132132        } else {
    … …  
    308308        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    309309        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
     310        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    310311
    311312        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    313314        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
    314315        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
    315         model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    316316
    317317        source->moments = pmMomentsAlloc ();
    … …  
    353353
    354354    // which extended source analyses should we perform?
    355     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    356     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    357     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    358     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     355    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     356    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     357    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     358    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    359359
    360360    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    396396        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    397397
     398# if (0)
    398399        // Petrosian measurements
    399400        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    476477        }
    477478
     479# endif
    478480        psArrayAdd (table, 100, row);
    479481        psFree (row);

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

    r21516 r25766  
    349349
    350350    // which extended source analyses should we perform?
    351     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    352     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    353     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    354     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     351    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     352    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     353    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     354    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    355355
    356356    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    410410        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    411411
     412# if (0)
    412413        // Petrosian measurements
    413414        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    501502            }
    502503        }
     504# endif
    503505
    504506        psArrayAdd (table, 100, row);

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

    r21516 r25766  
    300300
    301301    // which extended source analyses should we perform?
    302     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    303     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    304     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    305     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     302    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     303    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     304    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     305    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    306306
    307307    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    343343        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    344344
     345// XXX disable these outputs until we clean up the names
     346# if (0)
    345347        // Petrosian measurements
    346348        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    422424            }
    423425        }
     426# endif
    424427
    425428        psArrayAdd (table, 100, row);

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceIO_RAW.c

    r24581 r25766  
    119119                 logChi, logChiNorm,
    120120                 source[0].peak->SN,
    121                  model[0].radiusFit,
     121                 source[0].apRadius,
    122122                 source[0].pixWeight,
    123123                 model[0].nDOF,
    … …  
    181181                 log10(model[0].chisqNorm/model[0].nDOF),
    182182                 source[0].peak->SN,
    183                  model[0].radiusFit,
     183                 source[0].apRadius,
    184184                 source[0].pixWeight,
    185185                 model[0].nDOF,

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceMatch.c

    r24262 r25766  
    221221        } else {
    222222            // Match with the master list
    223             psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources
     223            psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree
    224224            long numMatch = 0;          // Number of matches
    225225

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceMoments.c

    r24577 r25766  
    3636#include "pmSource.h"
    3737
    38 bool pmSourceMoments_Old(pmSource *source, psF32 radius);
    39 
    4038/******************************************************************************
    4139pmSourceMoments(source, radius): this function takes a subImage defined in the
    … …  
    5048    pmSource->mask (optional)
    5149
    52 XXX: The peak calculations are done in image coords, not subImage coords.
    53 
    54 this version clips input pixels on S/N
    55 XXX EAM : this version returns false for several reasons
     50The peak calculations are done in image coords, not subImage coords.
     51
     52this version optionally clips input pixels on S/N
    5653*****************************************************************************/
    5754# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    … …  
    6461    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    6562
    66     // XXX supply sky in a different way?
     63    // use sky from moments if defined, 0.0 otherwise
    6764    psF32 sky = 0.0;
    6865    if (source->moments == NULL) {
    … …  
    9188    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
    9289    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
    93     // minimal round-off error in the sums
    94 
    95     psF32 xOff  = source->peak->x;
    96     psF32 yOff  = source->peak->y;
    97     psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    98     psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
     90    // minimal round-off error in the sums.  since these values are subtracted just to minimize
     91    // the dynamic range and are added back below, the exact value does not matter. these are
     92    // (int) so they can be used in the image index below.
     93
     94    int xOff  = source->peak->x;
     95    int yOff  = source->peak->y;
     96    int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
     97    int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
     98
     99    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     100    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     101    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     102    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     103    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    99104
    100105    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    126131            psF32 wDiff = *vWgt;
    127132
    128             // skip pixels below specified significance level
     133            // skip pixels below specified significance level.  this is allowed, but should be
     134            // avoided -- the over-weights the wings of bright stars compared to those of faint
     135            // stars.
    129136            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    130137            if (pDiff < 0) continue;
    131138
     139            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     140            // weighting over weights the sky for faint sources
    132141            if (sigma > 0.0) {
    133               // apply a pseudo-gaussian weight
    134 
    135               // XXX a lot of extra flops; can we do pre-calculate?
    136               psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    137               assert (z >= 0.0);
    138               psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
    139               psF32 weight  = 1.0 / t;
    140 
    141               // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
    142 
    143               wDiff *= weight;
    144               pDiff *= weight;
     142                // XXX a lot of extra flops; can we pre-calculate?
     143                psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     144                assert (z >= 0.0);
     145                psF32 weight  = exp(-z);
     146
     147                wDiff *= weight;
     148                pDiff *= weight;
    145149            }
    146150
    … …  
    154158            Y1  += yWght;
    155159
    156             // fprintf (stderr, " : %6.1f %6.1f  %8.1f %8.1f\n", X1, Y1, Sum, Var);
    157 
    158             peakPixel = PS_MAX (*vPix, peakPixel);
    159             numPixels++;
    160         }
     160            peakPixel = PS_MAX (*vPix, peakPixel);
     161            numPixels++;
     162        }
    161163    }
    162164
    … …  
    164166    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    165167    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    166         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
    167         return (false);
     168        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
     169        return (false);
    168170    }
    169171
    … …  
    172174    float My = Y1/Sum;
    173175    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
    174         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
    175         return (false);
     176        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     177        return (false);
    176178    }
    177179
    … …  
    179181
    180182    // add back offset of peak in primary image
    181     source->moments->Mx = Mx + xOff;
    182     source->moments->My = My + yOff;
     183    // also offset from pixel index to pixel coordinate
     184    // (the calculation above uses pixel index instead of coordinate)
     185    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     186    source->moments->Mx = Mx + xOff + 0.5;
     187    source->moments->My = My + yOff + 0.5;
    183188
    184189    source->moments->Sum = Sum;
    … …  
    205210    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    206211
    207     // center of mass in subimage
     212    // center of mass in subimage.  Note: the calculation below uses pixel index, so we do not
     213    // correct xCM, yCM to pixel coordinate here.
    208214    psF32 xCM = Mx + xPeak; // coord of peak in subimage
    209215    psF32 yCM = My + yPeak; // coord of peak in subimage
    … …  
    211217    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    212218
    213         psF32 yDiff = row - yCM;
     219        psF32 yDiff = row - yCM;
    214220        if (fabs(yDiff) > radius) continue;
    215221
    216         psF32 *vPix = source->pixels->data.F32[row];
    217         psF32 *vWgt = source->variance->data.F32[row];
    218         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    219 
    220         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    221             if (vMsk) {
    222                 if (*vMsk) {
    223                     vMsk++;
    224                     continue;
    225                 }
    226                 vMsk++;
    227             }
    228             if (isnan(*vPix)) continue;
    229 
    230             psF32 xDiff = col - xCM;
     222        psF32 *vPix = source->pixels->data.F32[row];
     223        psF32 *vWgt = source->variance->data.F32[row];
     224        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     225
     226        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     227            if (vMsk) {
     228                if (*vMsk) {
     229                    vMsk++;
     230                    continue;
     231                }
     232                vMsk++;
     233            }
     234            if (isnan(*vPix)) continue;
     235
     236            psF32 xDiff = col - xCM;
    231237            if (fabs(xDiff) > radius) continue;
    232238
    233             // radius is just a function of (xDiff, yDiff)
    234             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    235             psF32 r  = sqrt(r2);
    236             if (r > radius) continue;
    237 
    238             psF32 pDiff = *vPix - sky;
    239             psF32 wDiff = *vWgt;
    240 
    241             // XXX EAM : should this limit be user-defined?
    242 
    243             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    244             if (pDiff < 0) continue;
    245 
     239            // radius is just a function of (xDiff, yDiff)
     240            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     241            psF32 r  = sqrt(r2);
     242            if (r > radius) continue;
     243
     244            psF32 pDiff = *vPix - sky;
     245            psF32 wDiff = *vWgt;
     246
     247            // skip pixels below specified significance level.  this is allowed, but should be
     248            // avoided -- the over-weights the wings of bright stars compared to those of faint
     249            // stars.
     250            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     251            if (pDiff < 0) continue;
     252
     253            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     254            // weighting over weights the sky for faint sources
    246255            if (sigma > 0.0) {
    247               // apply a pseudo-gaussian weight
    248 
    249               // XXX a lot of extra flops; can we do pre-calculate?
    250               psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    251               assert (z >= 0.0);
    252               psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
    253               psF32 weight  = 1.0 / t;
    254 
    255               // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
    256 
    257               wDiff *= weight;
    258               pDiff *= weight;
     256                // XXX a lot of extra flops; can we do pre-calculate?
     257                psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     258                assert (z >= 0.0);
     259                psF32 weight  = exp(-z);
     260
     261                wDiff *= weight;
     262                pDiff *= weight;
    259263            }
    260264
    261             Sum += pDiff;
    262 
    263             psF32 x = xDiff * pDiff;
    264             psF32 y = yDiff * pDiff;
    265 
    266             psF32 xx = xDiff * x;
    267             psF32 xy = xDiff * y;
    268             psF32 yy = yDiff * y;
    269 
    270             psF32 xxx = xDiff * xx / r;
    271             psF32 xxy = xDiff * xy / r;
    272             psF32 xyy = xDiff * yy / r;
    273             psF32 yyy = yDiff * yy / r;
    274 
    275             psF32 xxxx = xDiff * xxx / r2;
    276             psF32 xxxy = xDiff * xxy / r2;
    277             psF32 xxyy = xDiff * xyy / r2;
    278             psF32 xyyy = xDiff * yyy / r2;
    279             psF32 yyyy = yDiff * yyy / r2;
    280 
    281             XX  += xx;
    282             XY  += xy;
    283             YY  += yy;
    284 
    285             XXX  += xxx;
    286             XXY  += xxy;
    287             XYY  += xyy;
    288             YYY  += yyy;
    289 
    290             XXXX  += xxxx;
    291             XXXY  += xxxy;
    292             XXYY  += xxyy;
    293             XYYY  += xyyy;
    294             YYYY  += yyyy;
    295         }
     265            Sum += pDiff;
     266
     267            psF32 x = xDiff * pDiff;
     268            psF32 y = yDiff * pDiff;
     269
     270            psF32 xx = xDiff * x;
     271            psF32 xy = xDiff * y;
     272            psF32 yy = yDiff * y;
     273
     274            psF32 xxx = xDiff * xx / r;
     275            psF32 xxy = xDiff * xy / r;
     276            psF32 xyy = xDiff * yy / r;
     277            psF32 yyy = yDiff * yy / r;
     278
     279            psF32 xxxx = xDiff * xxx / r2;
     280            psF32 xxxy = xDiff * xxy / r2;
     281            psF32 xxyy = xDiff * xyy / r2;
     282            psF32 xyyy = xDiff * yyy / r2;
     283            psF32 yyyy = yDiff * yyy / r2;
     284
     285            XX  += xx;
     286            XY  += xy;
     287            YY  += yy;
     288
     289            XXX  += xxx;
     290            XXY  += xxy;
     291            XYY  += xyy;
     292            YYY  += yyy;
     293
     294            XXXX  += xxxx;
     295            XXXY  += xxxy;
     296            XXYY  += xxyy;
     297            XYYY  += xyyy;
     298            YYYY  += yyyy;
     299        }
    296300    }
    297301
    … …  
    312316
    313317    if (source->moments->Mxx < 0) {
    314       fprintf (stderr, "error: neg second moment??\n");
     318        fprintf (stderr, "error: neg second moment??\n");
    315319    }
    316320    if (source->moments->Myy < 0) {
    317       fprintf (stderr, "error: neg second moment??\n");
     321        fprintf (stderr, "error: neg second moment??\n");
    318322    }
    319323
    320324    psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
    321              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    322              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    323              source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     325             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     326             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     327             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
    324328
    325329    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
    326              source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
    327 
    328     if (source->moments->Mxx < 0) {
    329         fprintf (stderr, "error: neg second moment??\n");
    330     }
    331     if (source->moments->Myy < 0) {
    332         fprintf (stderr, "error: neg second moment??\n");
    333     }
    334 
    335     // XXX TEST:
    336     // pmSourceMoments_Old (source, radius);
     330             source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
     331
    337332    return(true);
    338333}
    339 
    340 
    341 bool pmSourceMoments_Old(pmSource *source, psF32 radius)
    342 {
    343     PS_ASSERT_PTR_NON_NULL(source, false);
    344     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    345     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    346     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    347 
    348     psF32 sky = 0.0;
    349     if (source->moments == NULL) {
    350         source->moments = pmMomentsAlloc();
    351     } else {
    352         sky = source->moments->Sky;
    353     }
    354 
    355     // Sum = SUM (z - sky)
    356     // X1  = SUM (x - xc)*(z - sky)
    357     // X2  = SUM (x - xc)^2 * (z - sky)
    358     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    359     psF32 peakPixel = -PS_MAX_F32;
    360     psS32 numPixels = 0;
    361     psF32 Sum = 0.0;
    362     psF32 Var = 0.0;
    363     psF32 X1 = 0.0;
    364     psF32 Y1 = 0.0;
    365     psF32 X2 = 0.0;
    366     psF32 Y2 = 0.0;
    367     psF32 XY = 0.0;
    368     psF32 x  = 0;
    369     psF32 y  = 0;
    370     psF32 R2 = PS_SQR(radius);
    371 
    372     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    373     // image coordinates.  the source->pixels image has an offset relative to its parent of
    374     // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
    375     // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
    376     // minimal round-off error in the sums
    377 
    378     psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    379     psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
    380 
    381     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    382 
    383         psF32 *vPix = source->pixels->data.F32[row];
    384         psF32 *vWgt = source->variance->data.F32[row];
    385         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    386 
    387         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    388             if (vMsk) {
    389                 if (*vMsk) {
    390                     vMsk++;
    391                     continue;
    392                 }
    393                 vMsk++;
    394             }
    395             if (isnan(*vPix)) continue;
    396 
    397             psF32 xDiff = col - xPeak;
    398             psF32 yDiff = row - yPeak;
    399 
    400             // radius is just a function of (xDiff, yDiff)
    401             if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
    402 
    403             psF32 pDiff = *vPix - sky;
    404             psF32 wDiff = *vWgt;
    405 
    406             // XXX EAM : should this limit be user-defined?
    407             if (PS_SQR(pDiff) < wDiff) continue;
    408 
    409             Var += wDiff;
    410             Sum += pDiff;
    411 
    412             psF32 xWght = xDiff * pDiff;
    413             psF32 yWght = yDiff * pDiff;
    414 
    415             X1  += xWght;
    416             Y1  += yWght;
    417 
    418             X2  += xDiff * xWght;
    419             XY  += xDiff * yWght;
    420             Y2  += yDiff * yWght;
    421 
    422             peakPixel = PS_MAX (*vPix, peakPixel);
    423             numPixels++;
    424         }
    425     }
    426 
    427     // if we have less than (1/4) of the possible pixels, force a retry
    428     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    429     if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    430         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
    431         return (false);
    432     }
    433 
    434     psTrace ("psModules.objects", 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    435              sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    436 
    437     //
    438     // first moment X  = X1/Sum + xc
    439     // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    440     // Sxy             = XY / Sum
    441     //
    442     x = X1/Sum;
    443     y = Y1/Sum;
    444     if ((fabs(x) > radius) || (fabs(y) > radius)) {
    445         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n",
    446                  source->peak->x, source->peak->y);
    447         return (false);
    448     }
    449 
    450 # if (PS_TRACE_ON)
    451     float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0);
    452     float Sxy = XY / Sum;
    453     float Syy = PS_MAX(Y2/Sum - PS_SQR(y), 0);
    454 
    455     psTrace ("psModules.objects", 3,
    456              "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
    457              sky, Sum, x, y, Sxx, Sxy, Syy);
    458 # endif
    459 
    460     return(true);
    461 }
    462 

  • branches/eam_branches/20090820/psModules/src/objects/pmSourcePhotometry.c

    r21511 r25766  
    8484
    8585    // we must have a valid model
     86    // XXX allow aperture magnitudes for sources without a model
    8687    model = pmSourceGetModel (&isPSF, source);
    8788    if (model == NULL) {
    … …  
    9091    }
    9192
     93    // XXX handle negative flux, low-significance
    9294    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
    9395        SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0];
    … …  
    101103
    102104    // measure PSF model photometry
    103     if (psf->FluxScale) {
     105    // XXX TEST: do not use flux scale
     106    if (0 && psf->FluxScale) {
    104107        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
    105108        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
    106         if (!isfinite(fluxScale)) {
    107             // XXX objects on the edge can be slightly outside -- if we get an
    108             // error, use the full fit.
    109             psErrorClear();
    110             status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    111         } else {
    112             if (isfinite(fluxScale) && (fluxScale > 0.0)) {
    113                 source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
    114             } else {
    115                 source->psfMag = NAN;
    116             }
    117         }
     109        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
     110        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
     111        source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
    118112    } else {
    119113        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
    120114    }
    121115
    122     // if we have a collection of model fits, one of them is a pointer to modelEXT?
    123     // XXX not guaranteed
     116    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
    124117    if (source->modelFits) {
    125       for (int i = 0; i < source->modelFits->n; i++) {
    126         pmModel *model = source->modelFits->data[i];
    127         status = pmSourcePhotometryModel (&model->mag, model);
    128       }
    129       if (source->modelEXT) {
    130         source->extMag = source->modelEXT->mag;
    131       }
     118        bool foundEXT = false;
     119        for (int i = 0; i < source->modelFits->n; i++) {
     120            pmModel *model = source->modelFits->data[i];
     121            status = pmSourcePhotometryModel (&model->mag, model);
     122            if (model == source->modelEXT) foundEXT = true;
     123        }
     124        if (foundEXT) {
     125            source->extMag = source->modelEXT->mag;
     126        } else {
     127            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     128        }
    132129    } else {
    133       if (source->modelEXT) {
    134         status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
    135       }
     130        if (source->modelEXT) {
     131            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     132        }
    136133    }
    137134
    … …  
    160157    }
    161158
    162     // replace source flux
    163     // XXX need to be certain which model and size of mask for prior subtraction
    164     // XXX full model or just analytical?
    165     // XXX use pmSourceAdd instead?
    166     if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    167         pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    168     }
     159    // if we measure aperture magnitudes, the source must not currently be subtracted!
     160    psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!");
    169161
    170162    // if we are measuring aperture photometry and applying the growth correction,
    … …  
    201193    if (isfinite (source->apMag) && isPSF && psf) {
    202194        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
    203             source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
     195            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
    204196        }
    205197        if (mode & PM_SOURCE_PHOT_APCORR) {
     198            // XXX this should be removed -- we no longer fit for the 'sky bias'
    206199            rflux   = pow (10.0, 0.4*source->psfMag);
    207             source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
     200            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    208201        }
    209202    }
    … …  
    211204        psFree(flux);
    212205        psFree(mask);
    213     }
    214 
    215     // if source was originally subtracted, re-subtract object, leave local sky
    216     // XXX replace with pmSourceSub...
    217     if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    218         pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    219206    }
    220207
    … …  
    763750    return flux;
    764751}
    765 
    766 // XXX this is test code to verify the shift is doing the right thing (seems to be)
    767 # if (0)
    768 // measure centroid of unshifted gaussian (should be 16.0,16.0)
    769         {
    770           psImage *image = source->pixels;
    771           float xo = 0.0;
    772           float yo = 0.0;
    773           float xo2 = 0.0;
    774           float yo2 = 0.0;
    775           float no = 0.0;
    776           for (int j = 0; j < image->numRows; j++)
    777           {
    778             for (int i = 0; i < image->numCols; i++) {
    779               xo += i*image->data.F32[j][i];
    780               yo += j*image->data.F32[j][i];
    781               xo2 += i*i*image->data.F32[j][i];
    782               yo2 += j*j*image->data.F32[j][i];
    783               no += image->data.F32[j][i];
    784             }
    785           }
    786           xo /= no;
    787           yo /= no;
    788           xo2 = sqrt (xo2/no - xo*xo);
    789           yo2 = sqrt (yo2/no - yo*yo);
    790           fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
    791         }
    792 
    793 // measure centroid of unshifted gaussian (should be 16.0,16.0)
    794         {
    795           psImage *image = flux;
    796           float xo = 0.0;
    797           float yo = 0.0;
    798           float xo2 = 0.0;
    799           float yo2 = 0.0;
    800           float no = 0.0;
    801           for (int j = 0; j < image->numRows; j++)
    802           {
    803             for (int i = 0; i < image->numCols; i++) {
    804               xo += i*image->data.F32[j][i];
    805               yo += j*image->data.F32[j][i];
    806               xo2 += i*i*image->data.F32[j][i];
    807               yo2 += j*j*image->data.F32[j][i];
    808               no += image->data.F32[j][i];
    809             }
    810           }
    811           xo /= no;
    812           yo /= no;
    813           xo2 = sqrt (xo2/no - xo*xo);
    814           yo2 = sqrt (yo2/no - yo*yo);
    815           fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
    816         }
    817 # endif

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceVisual.c

    r23242 r25766  
    55#include <pslib.h>
    66#include "pmTrend2D.h"
     7#include "pmPSF.h"
     8#include "pmPSFtry.h"
     9#include "pmSource.h"
    710#include "pmSourceVisual.h"
    811
    … …  
    1518
    1619static int kapa1 = -1;
     20static int kapa2 = -1;
    1721static bool plotPSF = true;
    18 // static int kapa2 = -1;
    1922// static int kapa3 = -1;
    2023
    … …  
    2730bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi);
    2831
    29 
    30 bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
    31 
    32     KapaSection section;  // put the positive profile in one and the residuals in another?
     32bool pmSourceVisualPlotPSFMetric (pmPSFtry *psfTry) {
    3333
    3434    Graphdata graphdata;
    3535
    36     if (!pmVisualIsVisual() || !plotPSF) return true;
     36    if (!pmVisualIsVisual()) return true;
    3737
    3838    if (kapa1 == -1) {
    … …  
    4545    }
    4646
     47    KapaClearSections (kapa1);
     48    KapaInitGraph (&graphdata);
     49
     50    psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     51    psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     52    psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32);
     53
     54    graphdata.xmin = +32.0;
     55    graphdata.xmax = -32.0;
     56    graphdata.ymin = +32.0;
     57    graphdata.ymax = -32.0;
     58
     59    // construct the plot vectors
     60    int n = 0;
     61    for (int i = 0; i < psfTry->sources->n; i++) {
     62        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     63        x->data.F32[n] = psfTry->fitMag->data.F32[i];
     64        y->data.F32[n] = psfTry->metric->data.F32[i];
     65        dy->data.F32[n] = psfTry->metricErr->data.F32[i];
     66        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     67        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     68        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     69        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     70        n++;
     71    }
     72    x->n = y->n = dy->n = n;
     73
     74    float range;
     75    range = graphdata.xmax - graphdata.xmin;
     76    graphdata.xmax += 0.05*range;
     77    graphdata.xmin -= 0.05*range;
     78    range = graphdata.ymax - graphdata.ymin;
     79    graphdata.ymax += 0.05*range;
     80    graphdata.ymin -= 0.05*range;
     81
     82    // better choice for range?
     83    // graphdata.xmin = -17.0;
     84    // graphdata.xmax =  -9.0;
     85    graphdata.ymin = -0.51;
     86    graphdata.ymax = +0.51;
     87
     88    KapaSetLimits (kapa1, &graphdata);
     89
     90    KapaSetFont (kapa1, "helvetica", 14);
     91    KapaBox (kapa1, &graphdata);
     92    KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM);
     93    KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     94
     95    graphdata.color = KapaColorByName ("black");
     96    graphdata.ptype = 2;
     97    graphdata.size = 0.5;
     98    graphdata.style = 2;
     99    graphdata.etype |= 0x01;
     100
     101    KapaPrepPlot (kapa1, n, &graphdata);
     102    KapaPlotVector (kapa1, n, x->data.F32, "x");
     103    KapaPlotVector (kapa1, n, y->data.F32, "y");
     104    KapaPlotVector (kapa1, n, dy->data.F32, "dym");
     105    KapaPlotVector (kapa1, n, dy->data.F32, "dyp");
     106
     107    psFree (x);
     108    psFree (y);
     109    psFree (dy);
     110
     111    pmVisualAskUser(NULL);
     112    return true;
     113}
     114
     115// to see the structure of the psf model, place the sources in a fake image 1/10th the size
     116// at their appropriate relative location. later sources stomp on earlier sources
     117bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) {
     118
     119    if (!pmVisualIsVisual()) return true;
     120
     121    if (kapa2 == -1) {
     122        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
     123        if (kapa2 == -1) {
     124            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     125            pmVisualSetVisual(false);
     126            return false;
     127        }
     128    }
     129
     130    // create images 1/10 scale:
     131    psImage *image = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     132    psImage *model = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     133    psImage *resid = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     134    psImageInit (image, 0.0);
     135    psImageInit (model, 0.0);
     136    psImageInit (resid, 0.0);
     137
     138    for (int i = sources->n - 1; i >= 0; i--) {
     139        pmSource *source = sources->data[i];
     140        if (!source) continue;
     141        if (!source->pixels) continue;
     142
     143        pmSourceCacheModel (source, maskVal);
     144        if (!source->modelFlux) continue;
     145
     146        pmModel *srcModel = pmSourceGetModel (NULL, source);
     147        if (!model) continue;
     148
     149        float norm = srcModel->params->data.F32[PM_PAR_I0];
     150
     151        int Xo = 0.1*srcModel->params->data.F32[PM_PAR_XPOS];
     152        int Yo = 0.1*srcModel->params->data.F32[PM_PAR_YPOS];
     153
     154        // insert source pixels in the image at 1/10th offset
     155        for (int iy = 0; iy < source->pixels->numRows; iy++) {
     156            int jy = iy + Yo;
     157            if (jy >= image->numRows) continue;
     158            for (int ix = 0; ix < source->pixels->numCols; ix++) {
     159                int jx = ix + Xo;
     160                if (jx >= image->numCols) continue;
     161                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue;
     162                if (source->modelFlux->data.F32[iy][ix] < 0.001) continue;
     163                image->data.F32[jy][jx] = source->pixels->data.F32[iy][ix];
     164                model->data.F32[jy][jx] = source->modelFlux->data.F32[iy][ix];
     165                resid->data.F32[jy][jx] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     166            }
     167        }
     168    }
     169
     170    // KapaClearSections (kapa2);
     171    pmVisualScaleImage (kapa2, image, "image", 0, true);
     172    pmVisualScaleImage (kapa2, model, "model", 1, true);
     173    pmVisualScaleImage (kapa2, resid, "resid", 2, true);
     174
     175# ifdef DEBUG
     176    {
     177        psFits *fits = psFitsOpen ("image.fits", "w");
     178        psFitsWriteImage (fits, NULL, image, 0, NULL);
     179        psFitsClose (fits);
     180        fits = psFitsOpen ("model.fits", "w");
     181        psFitsWriteImage (fits, NULL, model, 0, NULL);
     182        psFitsClose (fits);
     183        fits = psFitsOpen ("resid.fits", "w");
     184        psFitsWriteImage (fits, NULL, resid, 0, NULL);
     185        psFitsClose (fits);
     186    }
     187# endif
     188
     189    psFree (image);
     190    psFree (model);
     191    psFree (resid);
     192
     193    pmVisualAskUser(NULL);
     194    return true;
     195}
     196
     197bool pmSourceVisualShowModelFit (pmSource *source) {
     198
     199    if (!pmVisualIsVisual()) return true;
     200    if (!source->pixels) return false;
     201    if (!source->modelFlux) return false;
     202
     203    if (kapa2 == -1) {
     204        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
     205        if (kapa2 == -1) {
     206            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     207            pmVisualSetVisual(false);
     208            return false;
     209        }
     210    }
     211
     212    // KapaClearSections (kapa2);
     213    pmVisualScaleImage (kapa2, source->pixels, "source", 0, false);
     214    pmVisualScaleImage (kapa2, source->modelFlux, "model", 1, false);
     215
     216    pmModel *model = pmSourceGetModel (NULL, source);
     217    float norm = model->params->data.F32[PM_PAR_I0];
     218
     219    psImage *resid = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     220    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     221        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     222            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     223                resid->data.F32[iy][ix] = NAN;
     224                continue;
     225            }
     226            resid->data.F32[iy][ix] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     227        }
     228    }
     229    pmVisualScaleImage (kapa2, resid, "resid", 2, false);
     230
     231    psFree (resid);
     232
     233    pmVisualAskUser(NULL);
     234    return true;
     235}
     236
     237bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
     238
     239    KapaSection section;  // put the positive profile in one and the residuals in another?
     240
     241    Graphdata graphdata;
     242
     243    if (!pmVisualIsVisual() || !plotPSF) return true;
     244
     245    if (kapa1 == -1) {
     246        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
     247        if (kapa1 == -1) {
     248            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     249            pmVisualSetVisual(false);
     250            return false;
     251        }
     252    }
     253
    47254    KapaClearPlots (kapa1);
    48255    KapaInitGraph (&graphdata);
    49256
    50     float min = +1e32;
    51     float max = -1e32;
    52     float Min = +1e32;
    53     float Max = -1e32;
     257    float Xmin = +1e32;
     258    float Xmax = -1e32;
     259    float Ymin = +1e32;
     260    float Ymax = -1e32;
     261    float Fmin = +1e32;
     262    float Fmax = -1e32;
    54263
    55264    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32);
    56265    psVector *model = psVectorAlloc (x->n, PS_TYPE_F32);
    57266
     267    psVector *xm = psVectorAlloc (x->n, PS_TYPE_F32);
     268    psVector *ym = psVectorAlloc (x->n, PS_TYPE_F32);
     269    psVector *Fm = psVectorAlloc (x->n, PS_TYPE_F32);
     270
     271    int n = 0;
    58272    for (int i = 0; i < x->n; i++) {
    59273        model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
    60274        resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
    61275        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    62         min = PS_MIN (min, resid->data.F32[i]);
    63         max = PS_MAX (max, resid->data.F32[i]);
    64         Min = PS_MIN (min, param->data.F32[i]);
    65         Max = PS_MAX (max, param->data.F32[i]);
    66     }
    67 
    68     psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32);
    69     psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32);
    70     psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32);
    71     psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32);
    72     psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
    73     for (int i = 0; i < x->n; i++) {
    74         xn->data.F32[i] = x->data.F32[i] / 5000.0;
    75         yn->data.F32[i] = y->data.F32[i] / 5000.0;
    76         zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
    77         Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
    78         Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
    79     }
     276        Xmin = PS_MIN (Xmin, x->data.F32[i]);
     277        Xmax = PS_MAX (Xmax, x->data.F32[i]);
     278        Ymin = PS_MIN (Ymin, y->data.F32[i]);
     279        Ymax = PS_MAX (Ymax, y->data.F32[i]);
     280        Fmin = PS_MIN (Fmin, param->data.F32[i]);
     281        Fmax = PS_MAX (Fmax, param->data.F32[i]);
     282        xm->data.F32[n] = x->data.F32[i];
     283        ym->data.F32[n] = y->data.F32[i];
     284        Fm->data.F32[n] = param->data.F32[i];
     285        n++;
     286    }
     287    xm->n = ym->n = Fm->n = n;
    80288
    81289    // view 1 on resid
    82     section.dx = 0.5;
    83     section.dy = 0.33;
     290    section.dx = 1.0;
     291    section.dy = 0.5;
    84292    section.x = 0.0;
    85293    section.y = 0.0;
    … …  
    88296    KapaSetSection (kapa1, &section);
    89297    psFree (section.name);
    90     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     298
     299    graphdata.color = KapaColorByName ("black");
     300    graphdata.xmin = Xmin;
     301    graphdata.xmax = Xmax;
     302    graphdata.ymin = Fmin;
     303    graphdata.ymax = Fmax;
     304
     305    {
     306        float range;
     307        range = graphdata.xmax - graphdata.xmin;
     308        graphdata.xmax += 0.05*range;
     309        graphdata.xmin -= 0.05*range;
     310        range = graphdata.ymax - graphdata.ymin;
     311        graphdata.ymax += 0.05*range;
     312        graphdata.ymin -= 0.05*range;
     313    }
     314
     315    KapaSetLimits (kapa1, &graphdata);
     316    KapaSetFont (kapa1, "helvetica", 14);
     317    KapaBox (kapa1, &graphdata);
     318    KapaSendLabel (kapa1, "X (pixels)", KAPA_LABEL_XM);
     319    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     320
     321    graphdata.ptype = 2;
     322    graphdata.size = 1.0;
     323    graphdata.style = 2;
     324    KapaPrepPlot (kapa1,   x->n, &graphdata);
     325    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     326    KapaPlotVector (kapa1, x->n, param->data.F32, "y");
     327
     328    graphdata.color = KapaColorByName ("red");
     329    graphdata.ptype = 1;
     330    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     331    KapaPlotVector (kapa1, xm->n, xm->data.F32, "x");
     332    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     333
     334    graphdata.color = KapaColorByName ("blue");
     335    graphdata.ptype = 1;
     336    KapaPrepPlot (kapa1,   x->n, &graphdata);
     337    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     338    KapaPlotVector (kapa1, x->n, model->data.F32, "y");
    91339
    92340    // view 2 on resid
    93     section.dx = 0.5;
    94     section.dy = 0.33;
    95     section.x = 0.5;
    96     section.y = 0.0;
     341    section.dx = 1.0;
     342    section.dy = 0.5;
     343    section.x = 0.0;
     344    section.y = 0.5;
    97345    section.name = NULL;
    98346    psStringAppend (&section.name, "a2");
    99347    KapaSetSection (kapa1, &section);
    100348    psFree (section.name);
    101     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    102 
    103     // view 3 on resid
    104     section.dx = 0.5;
    105     section.dy = 0.33;
    106     section.x = 0.0;
    107     section.y = 0.33;
    108     section.name = NULL;
    109     psStringAppend (&section.name, "a3");
    110     KapaSetSection (kapa1, &section);
    111     psFree (section.name);
    112     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    113 
    114     // view 4 on resid
    115     section.dx = 0.5;
    116     section.dy = 0.33;
    117     section.x = 0.5;
    118     section.y = 0.33;
    119     section.name = NULL;
    120     psStringAppend (&section.name, "a4");
    121     KapaSetSection (kapa1, &section);
    122     psFree (section.name);
    123     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    124 
    125     // view 5 on resid
    126     section.dx = 0.5;
    127     section.dy = 0.33;
    128     section.x = 0.0;
    129     section.y = 0.66;
    130     section.name = NULL;
    131     psStringAppend (&section.name, "a5");
    132     KapaSetSection (kapa1, &section);
    133     psFree (section.name);
    134     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    135 
    136     // view 6 on resid
    137     section.dx = 0.5;
    138     section.dy = 0.33;
    139     section.x = 0.5;
    140     section.y = 0.66;
    141     section.name = NULL;
    142     psStringAppend (&section.name, "a6");
    143     KapaSetSection (kapa1, &section);
    144     psFree (section.name);
    145     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     349
     350    graphdata.color = KapaColorByName ("black");
     351    graphdata.xmin = Ymin;
     352    graphdata.xmax = Ymax;
     353    graphdata.ymin = Fmin;
     354    graphdata.ymax = Fmax;
     355    {
     356        float range;
     357        range = graphdata.xmax - graphdata.xmin;
     358        graphdata.xmax += 0.05*range;
     359        graphdata.xmin -= 0.05*range;
     360        range = graphdata.ymax - graphdata.ymin;
     361        graphdata.ymax += 0.05*range;
     362        graphdata.ymin -= 0.05*range;
     363    }
     364
     365    KapaSetLimits (kapa1, &graphdata);
     366    KapaSetFont (kapa1, "helvetica", 14);
     367    KapaBox (kapa1, &graphdata);
     368    KapaSendLabel (kapa1, "Y (pixels)", KAPA_LABEL_XM);
     369    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     370
     371    graphdata.ptype = 2;
     372    graphdata.size = 1.0;
     373    graphdata.style = 2;
     374    KapaPrepPlot (kapa1,   y->n, &graphdata);
     375    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     376    KapaPlotVector (kapa1, y->n, param->data.F32, "y");
     377
     378    graphdata.color = KapaColorByName ("red");
     379    graphdata.ptype = 1;
     380    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     381    KapaPlotVector (kapa1, xm->n, ym->data.F32, "x");
     382    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     383
     384    graphdata.color = KapaColorByName ("blue");
     385    graphdata.ptype = 1;
     386    KapaPrepPlot (kapa1,   y->n, &graphdata);
     387    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     388    KapaPlotVector (kapa1, y->n, model->data.F32, "y");
    146389
    147390    psFree (resid);
    148 
    149     psFree (xn);
    150     psFree (yn);
    151     psFree (zn);
    152     psFree (Zn);
     391    psFree (model);
    153392
    154393    // pause and wait for user input:
    … …  
    159398}
    160399
    161 // send in normalized points
     400// Somewhat broken 3D plotting function (was used by pmSourceVisualPSFModelResid, but not anymore)
    162401bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi) {
     402
     403    return true;
    163404
    164405    psVector *xv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32);
    … …  
    192433    KapaSetLimits (myKapa, graphdata);
    193434
    194     // KapaSetFont (myKapa, "helvetica", 14);
    195     // KapaBox (myKapa, graphdata);
    196     // KapaSendLabel (myKapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    197     // KapaSendLabel (myKapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    198 
    199435    graphdata->color = KapaColorByName ("black");
    200436    graphdata->ptype = 100;

  • branches/eam_branches/20090820/psModules/src/objects/pmSourceVisual.h

    r23242 r25766  
    1818
    1919bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);
     20bool pmSourceVisualPlotPSFMetric (pmPSFtry *try);
     21bool pmSourceVisualShowModelFit (pmSource *source);
     22bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal);
    2023
    2124/// @}

  • branches/eam_branches/20090820/psModules/src/objects/pmTrend2D.c

    r21183 r25766  
    298298    return PM_TREND_NONE;
    299299}
     300
     301bool pmTrend2DPrintMap (pmTrend2D *trend) {
     302
     303    if (!trend->map) return false;
     304    if (!trend->map->map) return false;
     305
     306    for (int j = 0; j < trend->map->map->numRows; j++) {
     307        for (int i = 0; i < trend->map->map->numCols; i++) {
     308            fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
     309        }
     310        fprintf (stderr, "\t\t\t");
     311        for (int i = 0; i < trend->map->map->numCols; i++) {
     312            fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
     313        }
     314        fprintf (stderr, "\n");
     315    }
     316    return true;
     317}
     318

  • branches/eam_branches/20090820/psModules/src/objects/pmTrend2D.h

    r21183 r25766  
    9797pmTrend2DMode pmTrend2DModeFromString(psString name);
    9898
     99bool pmTrend2DPrintMap (pmTrend2D *trend);
     100
    99101/// @}
    100102# endif
Note: See TracChangeset for help on using the changeset viewer.