IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25754


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

check in changes from genes development branch : extensive changes to moments calculation, psf model generation, aperture residuals

Location:
trunk/psModules
Files:
6 deleted
43 edited
7 copied

Legend:

Unmodified
Added
Removed

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

    r24792 r25754  
    387387                psFree(regionString);
    388388                psFree(readout);
     389                psFree(iter);
    389390                return false;
    390391            }

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

    r24769 r25754  
    611611            if (cells->n == 0) {
    612612                // Nothing to do
     613                psFree (cells);
    613614                return true;
    614615            }
    … …  
    660661            if (cells->n == 0) {
    661662                // Nothing to do
     663                psFree (cells);
    662664                return true;
    663665            }
    … …  
    707709            if (cells->n == 0) {
    708710                // Nothing to do
     711                psFree (cells);
    709712                return true;
    710713            }
    711 
    712             return generateHDU(hdu, cells);
     714           
     715            bool status = generateHDU(hdu, cells);
     716            psFree (cells);
     717            return status;
    713718        }
    714719    default:

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

    r25738 r25754  
    9797                continue;
    9898            }
     99            // check that all params are valid:
     100            bool validParams = true;
     101            for (int n = 0; validParams && (n < normModel->params->n); n++) {
     102                if (n == PM_PAR_SKY) continue;
     103                if (n == PM_PAR_I0) continue;
     104                if (n == PM_PAR_XPOS) continue;
     105                if (n == PM_PAR_YPOS) continue;
     106                if (!isfinite(normModel->params->data.F32[n])) validParams = false;
     107            }
     108            if (!validParams) {
     109                psFree(normModel);
     110                continue;
     111            }           
    99112            if (circularise && !circulariseModel(normModel)) {
    100113                psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
    … …  
    112125            continue;
    113126        }
     127        // check that all params are valid:
     128        bool validParams = true;
     129        for (int n = 0; validParams && (n < fakeModel->params->n); n++) {
     130            if (n == PM_PAR_SKY) continue;
     131            if (n == PM_PAR_I0) continue;
     132            if (n == PM_PAR_XPOS) continue;
     133            if (n == PM_PAR_YPOS) continue;
     134            if (!isfinite(fakeModel->params->data.F32[n])) validParams = false;
     135        }
     136        if (!validParams) {
     137            psFree(fakeModel);
     138            continue;
     139        }               
    114140        if (circularise && !circulariseModel(fakeModel)) {
    115141            psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");

  • trunk/psModules/src/extras/pmVisual.c

    r23989 r25754  
    2222#include "pmSubtractionStamps.h"
    2323#include "pmTrend2D.h"
     24#include "pmPSF.h"
     25#include "pmPSFtry.h"
     26#include "pmSource.h"
    2427#include "pmFPAExtent.h"
    2528
    … …  
    8689{
    8790    char key[10];
    88     fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
     91    if (plotFlag) {
     92        fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) ");
     93    } else {
     94        fprintf (stdout, "[c]ontinue? [a]bort all visual plots? (c) ");
     95    }
    8996    if (!fgets(key, 8, stdin)) {
    9097        psWarning("Unable to read option");
    9198    }
    92     if (key[0] == 's') {
     99    if (plotFlag && (key[0] == 's')) {
    93100        *plotFlag = false;
    94101    }

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

    r25738 r25754  
    381381    options->poissonErrorsParams = true;
    382382    options->stats = psStatsAlloc(PSF_STATS);
    383     options->radius = maxRadius;
     383    options->fitRadius = maxRadius;
     384    options->apRadius = maxRadius; // XXX need to decide if aperture mags need a different radius
    384385    options->psfTrendMode = PM_TREND_MAP;
    385386    options->psfTrendNx = xOrder;

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

    r25738 r25754  
    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 \

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

    r25738 r25754  
    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
    … …  
    5454
    5555// 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
    5658psF32 PM_MODEL_FUNC(psVector *deriv,
    5759                    const psVector *params,
    … …  
    163165
    164166// make an initial guess for parameters
     167// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    165168bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    166169{
    … …  
    178181    psEllipseShape shape = psEllipseAxesToShape (axes);
    179182
    180     PAR[PM_PAR_SKY]  = moments->Sky;
     183    PAR[PM_PAR_SKY]  = 0.0;
    181184    PAR[PM_PAR_I0]   = peak->flux;
    182185    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    230233    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    231234    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
    232237    return (radius);
    233238}
    … …  
    340345bool PM_MODEL_FIT_STATUS (pmModel *model)
    341346{
    342     psF32 dP;
    343347    bool  status;
    344348
    … …  
    346350    psF32 *dPAR = model->dparams->data.F32;
    347351
    348     dP = 0;
    349     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    350     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    351     dP = sqrt (dP);
    352 
    353352    status = true;
    354     status &= (dP < 0.5);
    355353    status &= (PAR[PM_PAR_I0] > 0);
    356354    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    357355
    358     if (status)
    359         return true;
    360     return false;
     356    return status;
    361357}
    362358

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

    r25738 r25754  
    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
    … …  
    5454
    5555// 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
    5658psF32 PM_MODEL_FUNC(psVector *deriv,
    5759                    const psVector *params,
    … …  
    165167
    166168// make an initial guess for parameters
     169// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    167170bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    168171{
    … …  
    179182    psEllipseShape shape = psEllipseAxesToShape (axes);
    180183
    181     PAR[PM_PAR_SKY]  = moments->Sky;
     184    PAR[PM_PAR_SKY]  = 0.0;
    182185    PAR[PM_PAR_I0]   = peak->flux;
    183186    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    258261    // choose a z value guaranteed to be beyond our limit
    259262    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]);
    260264    float z1 = (1.0 / limit);
     265    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_I0]);
    261266    z1 = PS_MAX (z0, z1);
    262267    z0 = 0.0;
    … …  
    389394bool PM_MODEL_FIT_STATUS (pmModel *model)
    390395{
    391     psF32 dP;
    392396    bool  status;
    393397
    … …  
    395399    psF32 *dPAR = model->dparams->data.F32;
    396400
    397     dP = 0;
    398     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    399     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    400     dP = sqrt (dP);
    401 
    402401    status = true;
    403     status &= (dP < 0.5);
    404402    status &= (PAR[PM_PAR_I0] > 0);
    405403    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);

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

    r25738 r25754  
    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
    … …  
    4242# define ALPHA_M 0.666
    4343
     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
    4448// Lax parameter limits
    4549static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
    … …  
    5761
    5862static bool limitsApply = true;         // Apply limits?
    59 
    6063
    6164psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    182185
    183186// make an initial guess for parameters
     187// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    184188bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    185189{
    … …  
    206210    if (!isfinite(shape.sxy)) return false;
    207211
    208     // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
    209212    PAR[PM_PAR_SKY]  = 0.0;
    210213    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    292295
    293296    // choose a z value guaranteed to be beyond our limit
    294     float z0 = pow((1.0 / limit), (1.0 / ALPHA));
    295     float z1 = (1.0 / limit) / PAR[PM_PAR_7];
    296     z1 = PS_MAX (z0, z1);
    297     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;
    298301
    299302    // perform a type of bisection to find the value
    … …  
    420423bool PM_MODEL_FIT_STATUS (pmModel *model)
    421424{
    422 
    423     psF32 dP;
    424425    bool  status;
    425426
    … …  
    427428    psF32 *dPAR = model->dparams->data.F32;
    428429
    429     dP = 0;
    430     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    431     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    432     dP = sqrt (dP);
    433 
    434430    status = true;
    435 //    status &= (dP < 0.5);
    436431    status &= (PAR[PM_PAR_I0] > 0);
    437432    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);

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

    r25738 r25754  
    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
    … …  
    3838# define PM_MODEL_FIT_STATUS      pmModelFitStatus_QGAUSS
    3939# 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
    4044
    4145// Lax parameter limits
    … …  
    178182
    179183// make an initial guess for parameters
     184// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    180185bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    181186{
    … …  
    202207    if (!isfinite(shape.sxy)) return false;
    203208
    204     // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
    205209    PAR[PM_PAR_SKY]  = 0.0;
    206210    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    283287    // choose a z value guaranteed to be beyond our limit
    284288    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]);
    285290    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]);
    286292    z1 = PS_MAX (z0, z1);
    287293    z0 = 0.0;
    … …  
    409415bool PM_MODEL_FIT_STATUS (pmModel *model)
    410416{
    411 
    412     psF32 dP;
    413417    bool  status;
    414418
    … …  
    416420    psF32 *dPAR = model->dparams->data.F32;
    417421
    418     dP = 0;
    419     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    420     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    421     dP = sqrt (dP);
    422 
    423422    status = true;
    424 //    status &= (dP < 0.5);
    425423    status &= (PAR[PM_PAR_I0] > 0);
    426424    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);

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

    r25738 r25754  
    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
    … …  
    3939# define PM_MODEL_SET_LIMITS      pmModelSetLimits_RGAUSS
    4040
     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
    4145// Lax parameter limits
    4246static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
    … …  
    8791        dPAR[PM_PAR_SXY] = -q*X*Y;
    8892
    89         // 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
    9094        dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -5.0*t*log(z)*p*z;
    9195    }
    … …  
    172176
    173177// make an initial guess for parameters
     178// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    174179bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    175180{
    … …  
    196201    if (!isfinite(shape.sxy)) return false;
    197202
    198     PAR[PM_PAR_SKY]  = moments->Sky;
     203    PAR[PM_PAR_SKY]  = 0.0;
    199204    PAR[PM_PAR_I0]   = peak->flux;
    200205    PAR[PM_PAR_XPOS] = peak->xf;
    … …  
    276281    // choose a z value guaranteed to be beyond our limit
    277282    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]);
    278284    float z1 = (1.0 / limit);
     285    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
    279286    z1 = PS_MAX (z0, z1);
    280287    z0 = 0.0;
    … …  
    402409bool PM_MODEL_FIT_STATUS (pmModel *model)
    403410{
    404 
    405     psF32 dP;
    406411    bool  status;
    407412
    … …  
    409414    psF32 *dPAR = model->dparams->data.F32;
    410415
    411     dP = 0;
    412     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    413     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    414     dP = sqrt (dP);
    415 
    416416    status = true;
    417     status &= (dP < 0.5);
    418417    status &= (PAR[PM_PAR_I0] > 0);
    419418    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);

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

    r25738 r25754  
    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
    … …  
    4141# define PM_MODEL_FIT_STATUS      pmModelFitStatus_SERSIC
    4242# 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
    4347
    4448// Lax parameter limits
    … …  
    186190
    187191// make an initial guess for parameters
     192// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    188193bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
    189194{
    … …  
    210215    if (!isfinite(shape.sxy)) return false;
    211216
    212     // XXX PAR[PM_PAR_SKY]  = moments->Sky;
    213217    PAR[PM_PAR_SKY]  = 0.0;
    214218    PAR[PM_PAR_I0]   = peak->flux;
    … …  
    284288
    285289    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]);
    286291
    287292    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);
    288294
    289295    if (isnan(radius))
    … …  
    392398bool PM_MODEL_FIT_STATUS (pmModel *model)
    393399{
    394 
    395     psF32 dP;
    396400    bool  status;
    397401
    … …  
    399403    psF32 *dPAR = model->dparams->data.F32;
    400404
    401     dP = 0;
    402     dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    403     dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    404     dP = sqrt (dP);
    405 
    406405    status = true;
    407 //    status &= (dP < 0.5);
    408406    status &= (PAR[PM_PAR_I0] > 0);
    409407    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    410408
     409<<<<<<< .working
    411410    fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n",
    412411             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    413412
     413=======
     414>>>>>>> .merge-right.r25750
    414415    return status;
    415416}

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

    r23989 r25754  
    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

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

    r25738 r25754  
    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?
    … …  
    109110    new->nIter     = model->nIter;
    110111    new->flags     = model->flags;
    111     new->radiusFit = model->radiusFit;
     112    new->fitRadius = model->fitRadius;
    112113
    113114    for (int i = 0; i < new->params->n; i++) {
    … …  
    190191    psVector *params = model->params;
    191192
    192     psS32 imageCol;
    193     psS32 imageRow;
    194     psF32 pixelValue;
     193    float imageCol;
     194    float imageRow;
     195    float pixelValue;
    195196
    196197    // save original values; restore before returning
    … …  
    233234    psF32 **Rx = NULL;
    234235    psF32 **Ry = NULL;
    235     psImageMaskType **Rm = NULL;
     236    pmResidMaskType **Rm = NULL;
    236237
    237238    if (model->residuals) {
    238         DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
    239         DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
    240         Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
    241         Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
    242         Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
    243         Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
    244         if (Ro) {
    245             NX = model->residuals->Ro->numCols;
    246             NY = model->residuals->Ro->numRows;
    247         }
     239        DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
     240        DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
     241        Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
     242        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
     243        Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
     244        Rm = (model->residuals->mask) ? model->residuals->mask->data.PM_TYPE_RESID_MASK_DATA : NULL;
     245        if (Ro) {
     246            NX = model->residuals->Ro->numCols;
     247            NY = model->residuals->Ro->numRows;
     248        }           
    248249    }
    249250
    … …  
    256257                continue;
    257258
    258             // XXX should we use using 0.5 pixel offset?
    259259            // Convert to coordinate in parent image, with offset (dx,dy)
    260             imageCol = ix + image->col0 - dx;
    261             imageRow = iy + image->row0 - dy;
    262 
    263             x->data.F32[0] = (float) imageCol;
    264             x->data.F32[1] = (float) imageRow;
     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;
    265266
    266267            pixelValue = 0.0;

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

    r25738 r25754  
    107107    int nIter;                          ///< number of iterations to reach min
    108108    pmModelStatus flags;                ///< model status flags
    109     float radiusFit;                    ///< fit radius actually used
     109    float fitRadius;                    ///< fit radius actually used
    110110    pmResiduals *residuals;             ///< normalized PSF residuals
    111111

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

    r25738 r25754  
    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);

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

    r19999 r25754  
    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

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

    r24206 r25754  
    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;

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

    r24206 r25754  
    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 }

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

    r21183 r25754  
    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

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

    r24623 r25754  
    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;

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

    r24874 r25754  
    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

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

    r24579 r25754  
    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;

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

    r23487 r25754  
    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}

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

    r23487 r25754  
    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/// @}

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

    r23989 r25754  
    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;

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

    r23487 r25754  
    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];

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

    r25415 r25754  
    572572                psFree (xfitname);
    573573                psFree (outhead);
     574                psFree (deteffname);
    574575                return false;
    575576            }
    … …  
    583584        psFree (xfitname);
    584585        psFree (outhead);
     586        psFree (deteffname);
    585587        break;
    586588
    … …  
    938940                psFree (headname);
    939941                psFree (dataname);
     942                psFree (deteffname);
    940943                return true;
    941944            }
    … …  
    996999        psFree (headname);
    9971000        psFree (dataname);
     1001        psFree (deteffname);
    9981002        psFree (tableHeader);
    9991003        break;

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

    r24401 r25754  
    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);

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

    r24403 r25754  
    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);

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

    r21516 r25754  
    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);

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

    r21516 r25754  
    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);

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

    r24581 r25754  
    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,

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

    r24577 r25754  
    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 

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

    r21511 r25754  
    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

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

    r23242 r25754  
    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;

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

    r23242 r25754  
    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/// @}

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

    r21183 r25754  
    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

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

    r21183 r25754  
    9797pmTrend2DMode pmTrend2DModeFromString(psString name);
    9898
     99bool pmTrend2DPrintMap (pmTrend2D *trend);
     100
    99101/// @}
    100102# endif

  • trunk/psModules/test/objects/tap_pmGrowthCurve.c

    r24851 r25754  
    131131        source->mode = PM_SOURCE_MODE_PSFSTAR;
    132132
    133         source->modelPSF->radiusFit = 15.0;
     133        source->apRadius = 15.0;
    134134
    135135        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    136136        double refMag = source->apMag;
    137137
    138         source->modelPSF->radiusFit = 10.0;
    139         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    140         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    141 
    142         source->modelPSF->radiusFit = 8.0;
    143         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    144         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    145 
    146         source->modelPSF->radiusFit = 6.0;
     138        source->apRadius = 10.0;
     139        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     140        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     141
     142        source->apRadius = 8.0;
     143        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     144        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     145
     146        source->apRadius = 6.0;
    147147        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    148148        ok_float_tol(refMag - source->apMag, +0.0003, 0.0001, "growth offset is is %f", refMag - source->apMag);
    149149
    150         source->modelPSF->radiusFit = 4.0;
     150        source->apRadius = 4.0;
    151151        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    152152        ok_float_tol(refMag - source->apMag, +0.0020, 0.0001, "growth offset is is %f", refMag - source->apMag);
    153153
    154         source->modelPSF->radiusFit = 3.0;
     154        source->apRadius = 3.0;
    155155        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    156156        ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag);
    157157
    158         source->modelPSF->radiusFit = 2.0;
     158        source->apRadius = 2.0;
    159159        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    160160        ok_float_tol(refMag - source->apMag, -0.0075, 0.0001, "growth offset is is %f", refMag - source->apMag);
    … …  
    234234        source->mode = PM_SOURCE_MODE_PSFSTAR;
    235235
    236         source->modelPSF->radiusFit = 15.0;
     236        source->apRadius = 15.0;
    237237        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    238238        double refMag = source->apMag;
    239239
    240         source->modelPSF->radiusFit = 10.0;
    241         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    242         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    243 
    244         source->modelPSF->radiusFit = 8.0;
    245         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    246         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    247 
    248         source->modelPSF->radiusFit = 6.0;
     240        source->apRadius = 10.0;
     241        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     242        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     243
     244        source->apRadius = 8.0;
     245        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     246        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     247
     248        source->apRadius = 6.0;
    249249        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    250250        ok_float_tol(refMag - source->apMag, +0.0004, 0.0001, "growth offset is is %f", refMag - source->apMag);
    251251
    252         source->modelPSF->radiusFit = 4.0;
     252        source->apRadius = 4.0;
    253253        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    254254        ok_float_tol(refMag - source->apMag, +0.0026, 0.0001, "growth offset is is %f", refMag - source->apMag);
    255255
    256         source->modelPSF->radiusFit = 3.0;
     256        source->apRadius = 3.0;
    257257        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    258258        ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag);
    259259
    260         source->modelPSF->radiusFit = 2.0;
     260        source->apRadius = 2.0;
    261261        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    262262        ok_float_tol(refMag - source->apMag, -0.0103, 0.0001, "growth offset is is %f", refMag - source->apMag);
    … …  
    336336        source->mode = PM_SOURCE_MODE_PSFSTAR;
    337337
    338         source->modelPSF->radiusFit = 15.0;
     338        source->apRadius = 15.0;
    339339        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    340340        double refMag = source->apMag;
    341341
    342         source->modelPSF->radiusFit = 10.0;
    343         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    344         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    345 
    346         source->modelPSF->radiusFit = 8.0;
    347         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    348         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    349 
    350         source->modelPSF->radiusFit = 6.0;
     342        source->apRadius = 10.0;
     343        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     344        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     345
     346        source->apRadius = 8.0;
     347        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     348        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     349
     350        source->apRadius = 6.0;
    351351        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    352352        ok_float_tol(refMag - source->apMag, +0.0006, 0.0001, "growth offset is is %f", refMag - source->apMag);
    353353
    354         source->modelPSF->radiusFit = 4.0;
     354        source->apRadius = 4.0;
    355355        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    356356        ok_float_tol(refMag - source->apMag, +0.0038, 0.0001, "growth offset is is %f", refMag - source->apMag);
    357357
    358         source->modelPSF->radiusFit = 3.0;
    359         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    360         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    361 
    362         source->modelPSF->radiusFit = 2.0;
     358        source->apRadius = 3.0;
     359        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     360        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     361
     362        source->apRadius = 2.0;
    363363        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    364364        ok_float_tol(refMag - source->apMag, -0.0164, 0.0001, "growth offset is is %f", refMag - source->apMag);

  • trunk/psModules/test/objects/tap_pmModel.c

    r21471 r25754  
    7777        ok(model->nDOF == 0, "pmModelAlloc() set pmModel->nDOF correctly");
    7878        ok(model->nIter == 0, "pmModelAlloc() set pmModel->nIter correctly");
    79         ok(model->radiusFit == 0, "pmModelAlloc() set pmModel->radiusFit correctly");
     79        ok(model->fitRadius == 0, "pmModelAlloc() set pmModel->fitRadius correctly");
    8080        ok(model->flags == PM_MODEL_STATUS_NONE, "pmModelAlloc() set pmModel->flags correctly");
    8181        ok(model->residuals == NULL, "pmModelAlloc() set pmModel->residuals correctly");
    … …  
    132132        modelSrc->nIter = 3;
    133133        modelSrc->flags = PM_MODEL_STATUS_NONE;
    134         modelSrc->radiusFit = 4;
     134        modelSrc->fitRadius = 4;
    135135        pmModel *modelDst = pmModelCopy(modelSrc);
    136136        ok(modelDst != NULL && psMemCheckModel(modelDst), "pmModelCopy() returned a non-NULL pmModel");
    … …  
    139139        ok(modelDst->nIter == 3, "pmModelCopy() set the pmModel->nIter member correctly");
    140140        ok(modelDst->flags == PM_MODEL_STATUS_NONE, "pmModelCopy() set the pmModel->flags member correctly");
    141         ok(modelDst->radiusFit == 4, "pmModelCopy() set the pmModel->radiusFit member correctly");
     141        ok(modelDst->fitRadius == 4, "pmModelCopy() set the pmModel->fitRadius member correctly");
    142142
    143143        psFree(modelSrc);

  • trunk/psModules/test/objects/tap_pmModelUtils.c

    r15985 r25754  
    8181        ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly");
    8282        ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly");
    83         ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly");
     83        ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly");
    8484        ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly");
    8585        ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly");
    … …  
    140140        ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly");
    141141        ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly");
    142         ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly");
     142        ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly");
    143143        ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly");
    144144        ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly");

  • trunk/psModules/test/objects/tap_pmSourcePhotometry.c

    r9922 r25754  
    9696    source->modelPSF = pmModelFromPSF (modelRef, psf);
    9797    source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1;
    98     source->modelPSF->radiusFit = radius;
     98    source->apRadius = radius;
    9999
    100100    // measure photometry for centered source (fractional pix : 0.5,0.5)
Note: See TracChangeset for help on using the changeset viewer.