IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 19, 2012, 5:24:19 PM (14 years ago)
Author:
mhuber
Message:

merging latest r34040 trunk changes to branch

Location:
branches/meh_branches/ppstack_test
Files:
32 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/meh_branches/ppstack_test

  • branches/meh_branches/ppstack_test/psModules

  • branches/meh_branches/ppstack_test/psModules/src/camera/Makefile.am

    r24836 r34041  
    3030        pmReadoutStack.c \
    3131        pmReadoutFake.c \
    32         pmFPABin.c
     32        pmFPABin.c \
     33        pmFPAExpNumIO.c
    3334
    3435pkginclude_HEADERS = \
    … …  
    5960        pmReadoutStack.h \
    6061        pmReadoutFake.h \
    61         pmFPABin.h
     62        pmFPABin.h \
     63        pmFPAExpNumIO.h
    6264
    6365CLEANFILES = *~

  • branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.c

    r31671 r34041  
    382382        // Number of the file in list
    383383        psString num = NULL;            // Number to use
    384         psStringAppend(&num, "%03d", file->fileID);
     384        psStringAppend(&num, "%" PRId64, file->fileID);
    385385        psStringSubstitute(&newRule, num, "{FILE.ID}");
    386386        psFree(num);
    … …  
    543543    if (!strcasecmp(type, "PATTERN")) {
    544544        return PM_FPA_FILE_PATTERN;
     545    }
     546    if (!strcasecmp(type, "EXPNUM")) {
     547        return PM_FPA_FILE_EXPNUM;
    545548    }
    546549
    … …  
    589592      case PM_FPA_FILE_PATTERN:
    590593        return "PATTERN";
     594      case PM_FPA_FILE_EXPNUM:
     595        return "EXPNUM";
    591596      default:
    592597        return ("NONE");

  • branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.h

    r29833 r34041  
    5151    PM_FPA_FILE_PATTERN,
    5252    PM_FPA_FILE_LINEARITY,
     53    PM_FPA_FILE_EXPNUM,
    5354} pmFPAfileType;
    5455
    … …  
    113114
    114115    int fileIndex;                      // Index of file
    115     int fileID;                         // internal sequence number
     116    psS64 fileID;                       // internal sequence number
    116117
    117118    psS64 imageId, sourceId;            // Image and source identifiers

  • branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfileIO.c

    r29935 r34041  
    5353#include "pmSubtractionIO.h"
    5454#include "pmPatternIO.h"
     55#include "pmFPAExpNumIO.h"
    5556#include "pmConcepts.h"
    5657#include "pmConfigRun.h"
    … …  
    231232        status = pmAstromModelReadForView (view, file, config);
    232233        break;
     234      case PM_FPA_FILE_EXPNUM:
     235        status = pmExpNumRead(view, file, config);
     236        break;
    233237      case PM_FPA_FILE_ASTROM_REFSTARS:
    234238      case PM_FPA_FILE_JPEG:
    … …  
    291295      case PM_FPA_FILE_DARK:
    292296      case PM_FPA_FILE_PATTERN:
     297      case PM_FPA_FILE_EXPNUM:
    293298        {
    294299            // create FPA structure component based on view
    … …  
    500505        status = pmFPAviewWriteSourcePlot (view, file, config);
    501506        break;
     507
     508      case PM_FPA_FILE_EXPNUM:
     509        // when ppStack output's EXPNUM file it uses a file rule where the file type is MASK
     510        psError(PS_ERR_IO, true, "cannot write type EXPNUM (%s)", file->name);
     511        return false;
    502512
    503513      case PM_FPA_FILE_WCS:
    … …  
    562572      case PM_FPA_FILE_ASTROM_REFSTARS:
    563573      case PM_FPA_FILE_LINEARITY:
     574      case PM_FPA_FILE_EXPNUM:
    564575        psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
    565576        status = psFitsClose (file->fits);
    … …  
    636647      case PM_FPA_FILE_ASTROM_MODEL:
    637648      case PM_FPA_FILE_ASTROM_REFSTARS:
     649      case PM_FPA_FILE_EXPNUM:
    638650        psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
    639651        return true;
    … …  
    797809      case PM_FPA_FILE_ASTROM_REFSTARS:
    798810      case PM_FPA_FILE_LINEARITY:
     811      case PM_FPA_FILE_EXPNUM:
    799812        psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n",
    800813                 file->filename, file->name, view->chip, view->cell, view->readout);
    … …  
    9961009        status = pmAstromRefstarsWritePHU (view, file, config);
    9971010        break;
     1011      case PM_FPA_FILE_EXPNUM:
    9981012      case PM_FPA_FILE_ASTROM_MODEL:
    9991013      case PM_FPA_FILE_SX:

  • branches/meh_branches/ppstack_test/psModules/src/camera/pmReadoutFake.c

    r33415 r34041  
    162162            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
    163163        }
     164        //MEH - add some outputs to test --
     165        psTrace("psModules.camera", 10, "TESTING:: source at %f,%f: radius %f fakeradius %f minflux %f psf params later \n",
     166                fakeModel->params->data.F32[PM_PAR_XPOS],fakeModel->params->data.F32[PM_PAR_YPOS],radius,fakeRadius,minFlux);
     167        //MEH -- seems radius is passed as float 0.0.. hardcode test -
     168
    164169        if (radius > 0) {
    165170            fakeRadius = PS_MAX(fakeRadius, radius);
    … …  
    210215    const pmPSF *psf = args->data[7];         // PSF
    211216    float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux
    212     float radius = PS_SCALAR_VALUE(args->data[9], F32);  // Minimum radius
     217    float radius = PS_SCALAR_VALUE(args->data[9], S32);  // Minimum radius - typecast to float from S32 outside of PS_SCALAR_VALUE otherwise sets 0.0
    213218    bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF?
    214219    bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak?
    215220    int groupIndex = PS_SCALAR_VALUE(args->data[12], S32); // Group index
    216221    int cellIndex = PS_SCALAR_VALUE(args->data[13], S32);  // Cell index
     222
     223    //MEH
     224    psTrace("psModules.camera", 10, "TESTING:thread: radius %f minflux %f  \n",radius,minFlux);
    217225
    218226    return readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
    … …  
    242250    return old;
    243251}
    244 
    245252
    246253bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows,
    … …  
    287294    }
    288295
     296    psTrace("psModules.camera", 10, "TESTING:prethread: radius %d minflux %f threaded %d \n",radius,minFlux,threaded);
    289297    if (threaded) {
    290298        for (int i = 0; i < groups->groups->n; i++) {
    … …  
    303311                psArrayAdd(args, 1, (pmPSF*)psf);
    304312                PS_ARRAY_ADD_SCALAR(args, minFlux, PS_TYPE_F32);
     313                //MEH - S32 is proper type here
    305314                PS_ARRAY_ADD_SCALAR(args, radius, PS_TYPE_S32);
    306315                PS_ARRAY_ADD_SCALAR(args, circularise, PS_TYPE_U8);
    … …  
    349358}
    350359
    351 
    352360bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
    353361                              pmSourceMode sourceMask, const psVector *xOffset, const psVector *yOffset,
    … …  
    392400    mag->n = numGood;
    393401
     402    //MEH -- convert int->float for calling downstream since nothing calls Vectors directly? psphot does.. argh!
     403    psTrace("psModules.camera", 10, "TESTING:: radius %d minflux %f  \n",radius,minFlux);
     404
    394405    bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf,
    395406                                           minFlux, radius, circularise, normalisePeak);

  • branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.c

    r31067 r34041  
    107107        DETREND_STRING_CASE(ASTROM);
    108108        DETREND_STRING_CASE(NOISEMAP);
     109        DETREND_STRING_CASE(VIDEOMASK);
     110        DETREND_STRING_CASE(VIDEODARK);
    109111        DETREND_STRING_CASE(LINEARITY);
    110112    default:

  • branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.h

    r29833 r34041  
    3636    PM_DETREND_TYPE_ASTROM,
    3737    PM_DETREND_TYPE_NOISEMAP,
     38    PM_DETREND_TYPE_VIDEOMASK,
     39    PM_DETREND_TYPE_VIDEODARK,
    3840    PM_DETREND_TYPE_LINEARITY,
    3941} pmDetrendType;

  • branches/meh_branches/ppstack_test/psModules/src/extras/psVectorBracket.c

    r10610 r34041  
    5151        }
    5252    }
    53     // at this point, index[Nhi] >= key > index[Nlo]
    54     N = Nhi;
     53    N = (Nhi >= index->n) ? Nhi - 1 : Nhi;
     54    // at this point, index[N] >= key > index[Nlo]
    5555    while ((index->data.F32[N] >= key) && (N > Nlo)) {
    5656        N--;

  • branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtraction.c

    r33467 r34041  
    835835//#ifdef TESTING
    836836    //MEH fix index issue - index is supposed to be an image index but from where?
     837    //MEH - index conflict or changed in past?
    837838    for (int j = 0; j < kernels->num; j++) {
    838839        if (stamp->convolutions1) {

  • branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtractionStamps.c

    r33596 r34041  
    12321232        // XXX this is somewhat arbitrary...
    12331233        if (source->psfMagErr > 0.05) continue;
    1234         if (fabs(source->psfMag - source->apMag) > 0.5) continue;
     1234        if (isfinite(source->apMag)) {
     1235            if (fabs(source->psfMag - source->apMag) > 0.5) continue;
     1236        } else if (isfinite(source->apMagRaw)) {
     1237            if (fabs(source->psfMag - source->apMagRaw) > 0.5) continue;
     1238        } else {
     1239            // XXX: Should we carry on or drop this source?
     1240            // drop it for now
     1241            continue;
     1242        }
    12351243
    12361244        if (source->modelPSF) {

  • branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_QGAUSS.c

    r33415 r34041  
    402402    assert (psf->params->n > PM_PAR_YPOS);
    403403    assert (psf->params->n > PM_PAR_XPOS);
     404
     405    if (! isfinite(Io)) {
     406        fprintf(stderr, "non-finite Io passed to PM_MODEL_PARAMS_FROM_PSF\n");
     407        return false;
     408    }
    404409
    405410    PAR[PM_PAR_SKY]  = 0.0;

  • branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_SERSIC.c

    r33415 r34041  
    192192    psF32 z0 = PAR[PM_PAR_I0]*f1;
    193193    psF32 f0 = PAR[PM_PAR_SKY] + z0;
     194
     195    if (!isfinite(z0)) {
     196        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
     197        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
     198            PAR[5], PAR[6], PAR[7]);
     199    }
    194200
    195201    assert (isfinite(f2));

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmFootprintCullPeaks.c

    r33415 r34041  
    178178            psArray *myFP = pmFootprintsFind(subImg, threshold, 5);
    179179            if (!myFP) {
    180                 psWarning ("missing footprint?");
     180                psWarning ("missing footprint? threshold: %.f", threshold);
    181181                continue;
    182182            }
    183183            if (!myFP->n) {
    184                 psWarning ("empty footprint?");
     184                psWarning ("empty footprint? threshold: %.f", threshold);
    185185                psFree (myFP);
    186186                continue;

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmModelFuncs.h

    r29004 r34041  
    2929# define PM_MODEL_FUNCS_H
    3030
     31# define HAVE_MODEL_VAR 1
     32
    3133/// @addtogroup Objects Object Detection / Analysis Functions
    3234/// @{
    … …  
    4648
    4749typedef enum {
    48     PM_MODEL_OP_NONE    = 0x00,
    49     PM_MODEL_OP_FUNC    = 0x01,
    50     PM_MODEL_OP_RES0    = 0x02,
    51     PM_MODEL_OP_RES1    = 0x04,
    52     PM_MODEL_OP_FULL    = 0x07,
    53     PM_MODEL_OP_SKY     = 0x08,
    54     PM_MODEL_OP_CENTER  = 0x10,
    55     PM_MODEL_OP_NORM    = 0x20,
    56     PM_MODEL_OP_NOISE   = 0x40,
     50    PM_MODEL_OP_NONE     = 0x00,
     51    PM_MODEL_OP_FUNC     = 0x01,
     52    PM_MODEL_OP_RES0     = 0x02,
     53    PM_MODEL_OP_RES1     = 0x04,
     54    PM_MODEL_OP_FULL     = 0x07,
     55    PM_MODEL_OP_SKY      = 0x08,
     56    PM_MODEL_OP_CENTER   = 0x10,
     57    PM_MODEL_OP_NORM     = 0x20,
     58    PM_MODEL_OP_NOISE    = 0x40,
     59# if (HAVE_MODEL_VAR)
     60    PM_MODEL_OP_MODELVAR = 0x80,
     61# endif
    5762} pmModelOpMode;
    5863

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmModelUtils.c

    r33415 r34041  
    144144
    145145    *Io = source->peak->rawFlux;
     146
     147#ifndef ALLOW_NONFINITE_PEAK
     148    // Gene says fail of peak !finite
     149    if (!isfinite(*Io)) return false;
     150#else
     151    // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline
     152    // causing assertion failures
    146153    if (!isfinite(*Io) && !source->moments) return false;
    147154
    148155    *Io = source->moments->Peak;
    149156    if (!isfinite(*Io)) return false;
     157#endif
    150158
    151159    return true;

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c

    r33415 r34041  
    136136            sum += value;
    137137        }
     138    }
     139
     140    if (!(sum > 0.0)) {
     141        // Crazy PSF image print out some debugging information ...
     142        fprintf(stderr, "invalid kernel sum %f found by pmPCMkernelFromPSF\n", sum);    for (int j = psf->yMin; j <= psf->yMax; j++) {
     143            fprintf(stderr, "Row %d\n", j);
     144            for (int i = psf->xMin; i <= psf->xMax; i++) {
     145                double value = source->psfImage->data.F32[y0 + j][x0 + i];
     146                fprintf(stderr, "  %d %f\n", i, value);
     147            }
     148        }
     149        fflush(stderr);
     150        // ... but avoid the asssertion two lines down by escaping
     151        goto escape;
    138152    }
    139153    assert (sum > 0.0);

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmPSF.c

    r33415 r34041  
    455455        return NAN;
    456456    }
    457 
     457    //MEH -- multiple ways of doing this in past? add param details to log useful for ppStack?
     458   
    458459    // get the model full-width at half-max
    459460    float fwhmMajor = 2*model->modelRadius (model->params, 0.5);

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmPSFtryFitEXT.c

    r30621 r34041  
    7373            continue;
    7474        }
     75        // If mask object does not exist, mark the source as bad.
     76        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
     77        // psphot to exit with a fault.
     78        if (source->maskObj == NULL) {
     79            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
     80            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     81            continue;
     82        }
    7583
    7684        source->modelEXT = pmSourceModelGuess (source, options->type);
    … …  
    8997
    9098        // clear object mask to define valid pixels
    91         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     99        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    92100
    93101        // exclude the poor fits

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.c

    r33415 r34041  
    3939#include "pmSourceExtendedPars.h"
    4040#include "pmSourceDiffStats.h"
     41#include "pmSourcePhotometry.h"
    4142#include "pmSource.h"
    4243
    … …  
    5051    psFree(tmp->pixels);
    5152    psFree(tmp->variance);
     53# if (HAVE_MODEL_VAR)
     54    psFree(tmp->modelVar);
     55# endif
    5256    psFree(tmp->maskObj);
    5357    psFree(tmp->maskView);
    … …  
    7680    psFree (source->pixels);
    7781    psFree (source->variance);
     82# if (HAVE_MODEL_VAR)
     83    psFree (source->modelVar);
     84# endif
    7885    psFree (source->maskObj);
    7986    psFree (source->maskView);
    … …  
    8390    source->pixels = NULL;
    8491    source->variance = NULL;
     92# if (HAVE_MODEL_VAR)
     93    source->modelVar = NULL;
     94# endif
    8595    source->maskObj = NULL;
    8696    source->maskView = NULL;
    … …  
    112122    source->pixels = NULL;
    113123    source->variance = NULL;
     124# if (HAVE_MODEL_VAR)
     125    source->modelVar = NULL;
     126# endif
    114127    source->maskObj = NULL;
    115128    source->maskView = NULL;
    … …  
    159172    source->radialAper = NULL;
    160173    source->parent = NULL;
     174    source->tmpPtr = NULL;
    161175    source->imageID = -1;
     176    source->nFrames = 0;
    162177
    163178    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
    … …  
    197212    source->pixels   = in->pixels   ? psImageCopyView(NULL, in->pixels)   : NULL;
    198213    source->variance = in->variance ? psImageCopyView(NULL, in->variance) : NULL;
     214# if (HAVE_MODEL_VAR)
     215    source->modelVar = NULL;
     216# endif
    199217    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
    200218
    … …  
    10261044    bool addNoise = mode & PM_MODEL_OP_NOISE;
    10271045
     1046# if (HAVE_MODEL_VAR)
     1047    bool addModelVar = mode & PM_MODEL_OP_MODELVAR;
     1048    if (addModelVar) psAssert (source->modelVar, "programming error");
     1049# endif
     1050
    10281051    // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
    10291052    if (!addNoise && source->modelFlux) {
    … …  
    10471070        }
    10481071
     1072# if (HAVE_MODEL_VAR)
     1073        psF32 **target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
     1074# else
    10491075        psF32 **target = source->pixels->data.F32;
     1076# endif
    10501077
    10511078        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    … …  
    10621089            }
    10631090        }
     1091# if (HAVE_MODEL_VAR)
     1092        if (!addModelVar) {
     1093            if (add) {
     1094                source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1095            } else {
     1096                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     1097            }
     1098        }
     1099# else
    10641100        if (add) {
    10651101            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    … …  
    10671103            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    10681104        }
     1105# endif
    10691106        return true;
    10701107    }
    … …  
    10741111        target = source->variance;
    10751112    }
     1113# if (HAVE_MODEL_VAR)
     1114    if (addModelVar) {
     1115        target = source->modelVar;
     1116    }
     1117# endif
    10761118
    10771119    if (add) {
    10781120        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     1121# if (HAVE_MODEL_VAR)
     1122        if (!addNoise && !addModelVar) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1123# else
    10791124        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1125# endif
    10801126    } else {
    10811127        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     1128# if (HAVE_MODEL_VAR)
     1129        if (!addNoise && !addModelVar) source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     1130# else
    10821131        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     1132# endif
    10831133    }
    10841134
    … …  
    11561206    PAR[PM_PAR_SYY] = oldshape.sy;
    11571207    PAR[PM_PAR_SXY] = oldshape.sxy;
     1208
     1209    return true;
     1210}
     1211
     1212bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy)
     1213{
     1214//    assert (mode & PM_MODEL_OP_NOISE);
     1215    PS_ASSERT_PTR_NON_NULL(source, false);
     1216    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1217    PS_ASSERT_PTR_NON_NULL(target, false);
     1218
     1219    if (add) {
     1220        psTrace ("psphot", 3, "adding smoothed object at %f,%f\n", source->peak->xf, source->peak->yf);
     1221    } else {
     1222        psTrace ("psphot", 3, "removing smooted object at %f,%f\n", source->peak->xf, source->peak->yf);
     1223    }
     1224
     1225    pmModel *model = pmSourceGetModel(NULL, source);
     1226    if (!model) {
     1227        return false;
     1228    }
     1229    pmSourceSmoothOpModel (model, source, mode, target, sigma, add, maskVal, dx, dy);
     1230
     1231    return true;
     1232}
     1233
     1234bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy)
     1235{
     1236    bool status;
     1237    psEllipseShape oldshape;
     1238    psEllipseShape newshape;
     1239    psEllipseAxes axes;
     1240
     1241    if (add) {
     1242        psTrace ("psphot", 4, "adding smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1243    } else {
     1244        psTrace ("psphot", 4, "removing smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1245    }
     1246
     1247    psF32 *PAR = model->params->data.F32;
     1248
     1249    // Isn't this hanging around somewhere?
     1250    float oldFlux = NAN;
     1251    if (!pmSourcePhotometryModel (NULL, &oldFlux, model)) return false;
     1252
     1253    // save original values
     1254    float oldI0  = PAR[PM_PAR_I0];
     1255    float oldsxx = PAR[PM_PAR_SXX];
     1256    float oldsyy = PAR[PM_PAR_SYY];
     1257    float oldsxy = PAR[PM_PAR_SXY];
     1258
     1259    // Since we are going to scale the flux correctly we need to get our
     1260    // factors of sqrt(2) right
     1261    oldshape.sx  = oldsxx / M_SQRT2;
     1262    oldshape.sy  = oldsyy / M_SQRT2;
     1263    oldshape.sxy = oldsxy;
     1264
     1265    // XXX can this be done more intelligently?
     1266    if (oldI0 == 0.0) return false;
     1267    if (!isfinite(oldI0)) return false;
     1268
     1269    // increase size and height of source
     1270    axes = psEllipseShapeToAxes (oldshape, 20.0);
     1271    axes.major = sqrt(PS_SQR(axes.major) + PS_SQR(sigma));
     1272    axes.minor = sqrt(PS_SQR(axes.minor) + PS_SQR(sigma));
     1273    newshape = psEllipseAxesToShape (axes);
     1274    PAR[PM_PAR_SXX] = newshape.sx * M_SQRT2;
     1275    PAR[PM_PAR_SYY] = newshape.sy * M_SQRT2;
     1276    PAR[PM_PAR_SXY] = newshape.sxy;
     1277
     1278    PAR[PM_PAR_I0]  = 1.0;
     1279
     1280    float newFlux;
     1281    if (!pmSourcePhotometryModel (NULL, &newFlux, model)) return false;
     1282
     1283    PAR[PM_PAR_I0]  = oldFlux / newFlux;
     1284
     1285    if (add) {
     1286        status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1287    } else {
     1288        status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1289    }
     1290
     1291    // restore original values
     1292    PAR[PM_PAR_I0]  = oldI0;
     1293    PAR[PM_PAR_SXX] = oldsxx;
     1294    PAR[PM_PAR_SYY] = oldsyy;
     1295    PAR[PM_PAR_SXY] = oldsxy;
    11581296
    11591297    return true;

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.h

    r33415 r34041  
    7272    psImage *pixels;                    ///< Rectangular region including object pixels.
    7373    psImage *variance;                  ///< Image variance.
     74# if (HAVE_MODEL_VAR)
     75    psImage *modelVar;                  ///< variance based on current models
     76# endif
    7477    psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
    7578    psImage *maskView;                  ///< view into global image mask for this object region
    … …  
    117120    psArray *radialAper;                ///< radial flux in circular apertures
    118121    pmSource *parent;                   ///< reference to the master source from which this is derived
     122    psPtr *tmpPtr;                      ///< pointer that may be used to store data in a particular module. e.g. psphotKronIterate.
    119123    int imageID;
     124    psU16 nFrames;
    120125};
    121126
    … …  
    301306bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
    302307
     308bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
     309bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
     310
    303311bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
    304312bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO.c

    r33415 r34041  
    10691069                for (long i = sources->n -1; i >= 0; i--) {
    10701070                    pmSource *source = sources->data[i];
     1071                    if (source->seq < 0) {
     1072                        // This can happen cmf files that have been corrupted
     1073                        psError(PS_ERR_IO, true, "seq < 0 for source %ld: Suspect %s is corrupt", i, file->origname);
     1074                        return false;
     1075                    }
    10711076                    if (source->seq > seq_max) {
    10721077                        seq_max = source->seq;

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF.c.in

    r33415 r34041  
    8080    table = psArrayAllocEmpty (sources->n);
    8181
    82     short nImageOverlap;
    8382    float magOffset;
    8483    float zeroptErr;
    8584    float fwhmMajor;
    8685    float fwhmMinor;
    87     pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     86    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    8887
    8988    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    179178
    180179        // XXX not sure how to get this : need to load Nimages with weight?
    181         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
     180        @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    182181        @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    183182

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

    r33415 r34041  
    8181    }
    8282
    83     short nImageOverlap;
    8483    float magOffset;
    8584    float zeroptErr;
    8685    float fwhmMajor;
    8786    float fwhmMinor;
    88     pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     87    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    8988
    9089    table = psArrayAllocEmpty (sources->n);
    … …  
    161160        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
    162161
    163         // XXX not sure how to get this : need to load Nimages with weight?
    164         psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
     162        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    165163        psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    166164

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

    r33415 r34041  
    8282    table = psArrayAllocEmpty (sources->n);
    8383
    84     short nImageOverlap;
    8584    float magOffset;
    8685    float zeroptErr;
    8786    float fwhmMajor;
    8887    float fwhmMinor;
    89     pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     88    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    9089
    9190    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    180179
    181180        // XXX not sure how to get this : need to load Nimages with weight?
    182         psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
     181        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    183182        psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    184183

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r33415 r34041  
    8484    table = psArrayAllocEmpty (sources->n);
    8585
    86     short nImageOverlap;
    8786    float magOffset;
    8887    float zeroptErr;
    8988    float fwhmMajor;
    9089    float fwhmMinor;
    91     pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     90    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    9291
    9392    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    170169        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    171170
    172         // XXX not sure how to get this : need to load Nimages with weight?
    173         psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
     171        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    174172        psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    175173

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceMoments.c

    r33415 r34041  
    319319    }
    320320
    321     source->moments->Mrf = RF/RS;
    322321    source->moments->Mrh = RH/RS;
    323322
    324     // if Mrf (first radial moment) is very small, we are getting into low-significance
     323    // if Mrf = RF/RS (first radial moment) is very small, we are getting into low-significance
    325324    // territory.  saturate at minKronRadius.  conversely, if Mrf is >> radius for faint
    326325    // sources, we are clearly making an error.  saturate at radius.
    327     float kronRefRadius = MAX(minKronRadius, source->moments->Mrf);
     326    float kronRefRadius = MAX(minKronRadius, RF/RS);
    328327    if (source->moments->SN < 10) {
    329328        kronRefRadius = MIN(radius, kronRefRadius);
    330329    }
     330    source->moments->Mrf = kronRefRadius;
    331331
    332332    // *** now calculate the kron flux values using the 1st radial moment

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.c

    r33415 r34041  
    4040#include "pmSourceOutputs.h"
    4141
    42 bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) {
     42bool pmSourceOutputsCommonValues (float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) {
    4343
    4444    pmFPA  *fpa  = readout->parent->parent->parent;
    … …  
    7272    }
    7373
    74     *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS");
    7574    return true;
    7675

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.h

    r33415 r34041  
    5656} pmSourceOutputsMoments;
    5757
    58 bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header);
     58bool pmSourceOutputsCommonValues (float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header);
    5959
    6060bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset);

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.c

    r33415 r34041  
    899899}
    900900
     901# if (HAVE_MODEL_VAR)
     902double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
     903# else
    901904double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
     905# endif
    902906{
    903907    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
     908# if (HAVE_MODEL_VAR)
     909    double flux = 0, wt = 1.0, factor = 0;
     910# else
    904911    double flux = 0, wt = 0, factor = 0;
     912# endif
    905913
    906914    const psImage *Pi = Mi->modelFlux;
    907915    assert (Pi != NULL);
     916
     917# if (HAVE_MODEL_VAR)
     918    const psImage *Wi = NULL;
     919    switch (fitVarMode) {
     920      case PM_SOURCE_PHOTFIT_CONST:
     921        break;
     922      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     923        Wi = Mi->variance;
     924        psAssert (Wi, "programming error");
     925        break;
     926      case PM_SOURCE_PHOTFIT_MODEL_VAR:
     927        Wi = Mi->modelVar;
     928        psAssert (Wi, "programming error");
     929        break;
     930      case PM_SOURCE_PHOTFIT_NONE:
     931        psAbort("programming error");
     932    }   
     933# else
    908934    const psImage *Wi = Mi->variance;
    909935    if (!unweighted_sum) {
    910936        assert (Wi != NULL);
    911937    }
     938# endif
    912939    const psImage *Ti = Mi->maskObj;
    913940    assert (Ti != NULL);
    … …  
    917944            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    918945                continue;
     946# if (HAVE_MODEL_VAR)
     947            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
     948                wt = covarFactor * Wi->data.F32[yi][xi];
     949                if (wt == 0) continue;
     950            }
     951# else
    919952            if (!unweighted_sum) {
    920953                wt = covarFactor * Wi->data.F32[yi][xi];
    … …  
    922955                    continue;
    923956            }
     957# endif
    924958
    925959            switch (term) {
    … …  
    937971            }
    938972
     973# if (HAVE_MODEL_VAR)
     974            // wt is 1.0 for CONST
     975            flux += (factor * Pi->data.F32[yi][xi]) / wt;
     976# else
    939977            if (unweighted_sum) {
    940978                flux += (factor * Pi->data.F32[yi][xi]);
    … …  
    942980                flux += (factor * Pi->data.F32[yi][xi]) / wt;
    943981            }
     982# endif
    944983        }
    945984    }
    … …  
    947986}
    948987
     988# if (HAVE_MODEL_VAR)
     989double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
     990# else
    949991double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
     992# endif
    950993{
    951994    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    955998    int xIs, xJs, yIs, yJs;
    956999    int xIe, yIe;
     1000# if (HAVE_MODEL_VAR)
     1001    double flux;
     1002    double wt = 1.0;
     1003# else
    9571004    double flux, wt;
     1005# endif
    9581006
    9591007    const psImage *Pi = Mi->modelFlux;
    … …  
    9621010    assert (Pj != NULL);
    9631011
     1012# if (HAVE_MODEL_VAR)
     1013    const psImage *Wi = NULL;
     1014    switch (fitVarMode) {
     1015      case PM_SOURCE_PHOTFIT_CONST:
     1016        break;
     1017      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1018        Wi = Mi->variance;
     1019        psAssert (Wi, "programming error");
     1020        break;
     1021      case PM_SOURCE_PHOTFIT_MODEL_VAR:
     1022        Wi = Mi->modelVar;
     1023        psAssert (Wi, "programming error");
     1024        break;
     1025      case PM_SOURCE_PHOTFIT_NONE:
     1026        psAbort("programming error");
     1027    }   
     1028# else
    9641029    const psImage *Wi = Mi->variance;
    9651030    if (!unweighted_sum) {
    9661031        assert (Wi != NULL);
    9671032    }
     1033# endif
    9681034
    9691035    const psImage *Ti = Mi->maskObj;
    … …  
    9951061                continue;
    9961062
     1063# if (HAVE_MODEL_VAR)
     1064            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
     1065            switch (fitVarMode) {
     1066              case PM_SOURCE_PHOTFIT_CONST:
     1067                wt = 1.0;
     1068                break;
     1069              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1070              case PM_SOURCE_PHOTFIT_MODEL_VAR:
     1071                wt = covarFactor * Wi->data.F32[yi][xi];
     1072                break;
     1073              case PM_SOURCE_PHOTFIT_NONE:
     1074                psAbort("programming error");
     1075            }
     1076            // skip pixels with nonsense weight values
     1077            if (wt <= 0) continue;
     1078
     1079            flux += value / wt;
     1080# else
    9971081            // XXX skip the nonsense weight pixels?
    9981082            if (unweighted_sum) {
    … …  
    10041088                }
    10051089            }
     1090# endif
    10061091        }
    10071092    }
    … …  
    10091094}
    10101095
     1096# if (HAVE_MODEL_VAR)
     1097double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
     1098# else
    10111099double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
     1100# endif
    10121101{
    10131102    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    … …  
    10171106    int xIs, xJs, yIs, yJs;
    10181107    int xIe, yIe;
     1108# if (HAVE_MODEL_VAR)
     1109    double flux;
     1110    double wt = 1.0;
     1111# else
    10191112    double flux, wt;
     1113# endif
    10201114
    10211115    const psImage *Pi = Mi->pixels;
    … …  
    10241118    assert (Pj != NULL);
    10251119
     1120# if (HAVE_MODEL_VAR)
     1121    const psImage *Wi = NULL;
     1122    switch (fitVarMode) {
     1123      case PM_SOURCE_PHOTFIT_CONST:
     1124        break;
     1125      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1126        Wi = Mi->variance;
     1127        psAssert (Wi, "programming error");
     1128        break;
     1129      case PM_SOURCE_PHOTFIT_MODEL_VAR:
     1130        Wi = Mi->modelVar;
     1131        psAssert (Wi, "programming error");
     1132        break;
     1133      case PM_SOURCE_PHOTFIT_NONE:
     1134        psAbort("programming error");
     1135    }   
     1136# else
    10261137    const psImage *Wi = Mi->variance;
    10271138    if (!unweighted_sum) {
    10281139        assert (Wi != NULL);
    10291140    }
     1141# endif
    10301142
    10311143    const psImage *Ti = Mi->maskObj;
    … …  
    10571169                continue;
    10581170
     1171# if (HAVE_MODEL_VAR)
     1172            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
     1173            switch (fitVarMode) {
     1174              case PM_SOURCE_PHOTFIT_CONST:
     1175                wt = 1.0;
     1176                break;
     1177              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1178              case PM_SOURCE_PHOTFIT_MODEL_VAR:
     1179                wt = covarFactor * Wi->data.F32[yi][xi];
     1180                break;
     1181              case PM_SOURCE_PHOTFIT_NONE:
     1182                psAbort("programming error");
     1183            }
     1184            // skip pixels with nonsense weight values
     1185            if (wt <= 0) continue;
     1186
     1187            flux += value / wt;
     1188
     1189# else
    10591190            // XXX skip the nonsense weight pixels?
    10601191            if (unweighted_sum) {
    … …  
    10661197                }
    10671198            }
     1199# endif
    10681200        }
    10691201    }

  • branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.h

    r31670 r34041  
    3838} pmSourcePhotometryMode;
    3939
     40# if (HAVE_MODEL_VAR)
     41typedef enum {
     42    PM_SOURCE_PHOTFIT_NONE       = 0,
     43    PM_SOURCE_PHOTFIT_CONST      = 1,
     44    PM_SOURCE_PHOTFIT_IMAGE_VAR  = 2,
     45    PM_SOURCE_PHOTFIT_MODEL_VAR  = 3,
     46} pmSourceFitVarMode;
     47# endif
     48
    4049bool pmSourcePhotometryModel(
    4150    float *fitMag,                      ///< integrated fit magnitude
    … …  
    7584bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    7685
     86# if (HAVE_MODEL_VAR)
     87double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
     88double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
     89double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
     90# else
    7791double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    7892double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    7993double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
     94# endif
    8095
    8196bool pmSourceNeighborFlags (pmSource *source);
Note: See TracChangeset for help on using the changeset viewer.