IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30621


Ignore:
Timestamp:
Feb 13, 2011, 11:59:42 AM (15 years ago)
Author:
eugene
Message:

do not cull peaks with only one entry; skip the section that culls duplicate peaks -- this was dropping many valid peaks and seems unneeded; save radial apertures in their own extension; fixed logic on culling of faint peaks

Location:
trunk/psModules/src/objects
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects

    • Property svn:mergeinfo deleted

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

    r29004 r30621  
    255255}
    256256
     257// create a new footprint with the same span set as the input footprint
     258pmFootprint *pmFootprintCopyData(pmFootprint *inFoot, psImage *image) {
     259
     260    pmFootprint *outFoot = pmFootprintAlloc(inFoot->nspans, image);
     261    for (int i = 0; i < inFoot->nspans; i++) {
     262        pmSpan *span = inFoot->spans->data[i];
     263        pmFootprintAddSpan(outFoot, span->y, span->x0, span->x1);
     264    }
     265    return outFoot;
     266}
     267
    257268/************************************************************************************************************/

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

    r29004 r30621  
    2727
    2828bool pmFootprintAllocEmptySpans (pmFootprint *footprint, int nSpans);
     29
     30pmFootprint *pmFootprintCopyData(pmFootprint *inFoot, psImage *image);
    2931
    3032pmFootprint *pmFootprintNormalize(pmFootprint *fp);

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

    r29004 r30621  
    2424/*
    2525 * Given a psArray of pmFootprints and another of pmPeaks, assign the peaks to the
    26  * footprints in which that fall; if they _don't_ fall in a footprint, add a suitable
     26 * footprints in which they fall; if they _don't_ fall in a footprint, add a suitable
    2727 * one to the list.
    2828 */
    … …  
    8383
    8484        // XXX are we allowed to have peak-less footprints??
    85         if (!fp->peaks->n) continue;
     85        if (fp->peaks->n == 0) continue;
     86        if (fp->peaks->n == 1) continue;
    8687
    8788        fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
    8889
     90        // XXX check for an assert on duplicates (I don't think they can happen, but
     91        // let's double check for now)
     92
     93        for (int j = 1; j < fp->peaks->n; j++) {
     94            psAssert (fp->peaks->data[j] != fp->peaks->data[j-1], "duplicate peak!");
     95        }
     96        continue;
     97
     98        // XXX WHY am I culling duplicates?  how can there be duplicates?
    8999        // XXX EAM : the algorithm below should be much faster than using psArrayRemove if
    90100        // the number of peaks in the footprint is large, or if there are no duplicates.
    … …  
    93103
    94104        // track the number of good peaks in the footprint
    95         int nGood = 1;
     105        int lastGood = 0;
    96106
    97107        // check for duplicates
    … …  
    100110        // (if sorted list has A, B, A, B ...)
    101111        for (int j = 1; j < fp->peaks->n; j++) {
    102             if (fp->peaks->data[j] == fp->peaks->data[nGood]) {
     112            if (fp->peaks->data[j] == fp->peaks->data[lastGood]) {
    103113                // everything on the array has its own mem reference; free and drop this one
    104114                psFree (fp->peaks->data[j]);
    105115                fp->peaks->data[j] = NULL;
    106116            } else {
    107                 nGood ++;
     117                lastGood ++;
    108118            }
    109119        }
     120        int nGood = lastGood + 1;
    110121
    111122        // no deleted peaks, go to next footprint

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

    r29004 r30621  
    2828  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    2929  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    30   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
     30  * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
    3131  * starting point, discard the peak.
    3232  */
    … …  
    9999        float threshold = flux - nsigma_delta*stdev_pad;
    100100
    101         if (isnan(threshold) || threshold < min_threshold) {
     101        if (isnan(threshold)) {
    102102            // min_threshold is assumed to be below the detection threshold,
    103103            // so all the peaks are pmFootprint, and this isn't the brightest
    … …  
    109109            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    110110        }
     111
     112        if (threshold < min_threshold) {
     113            threshold = min_threshold;
     114        }
    111115
    112116        // init peakFootprint here?

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

    r29004 r30621  
    7070    tmp->fitRadius = 0;
    7171    tmp->flags = PM_MODEL_STATUS_NONE;
    72     tmp->residuals = NULL;              // XXX should the model own this memory?
     72    tmp->residuals = NULL;              // do not free: the model does not own this memory
     73    tmp->isPCM = false;
    7374
    7475    psS32 Nparams = pmModelClassParameterCount(type);
    … …  
    108109pmModel *pmModelCopy (pmModel *model)
    109110{
    110     PS_ASSERT_PTR_NON_NULL(model, NULL);
    111 
     111    if (model == NULL) {
     112        return NULL;
     113    }
    112114    pmModel *new = pmModelAlloc (model->type);
    113115

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

    r29004 r30621  
    4343    float fitRadius;                    ///< fit radius actually used
    4444    pmResiduals *residuals;             ///< normalized PSF residuals
     45    bool isPCM;                         ///< is this model fitted with PSF-convolution?
    4546
    4647    // functions for this model which depend on the model class

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

    r29004 r30621  
    8383    pmPCMdata *pcm);
    8484
     85psKernel *pmPCMkernelFromPSF (pmSource *source, int nPix);
     86
    8587bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    8688

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

    r29004 r30621  
    8484
    8585    for (int i = 0; i < sources->n; i++) {
    86         test->sources->data[i] = pmSourceCopy (sources->data[i]);
     86        pmSource *sourceOld = sources->data[i];
     87        pmSource *sourceNew = pmSourceCopy (sourceOld);
     88
     89        // save a reference so we can get back to the original
     90        // this is specifically used in psphotChooosePSF to unflag the candidate PSF sources
     91        // which were not actually used to generate a PSF model
     92        sourceNew->parent = sourceOld;
     93
     94        test->sources->data[i] = sourceNew;
    8795    }
    8896

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

    r29004 r30621  
    6565        if (!source->moments) {
    6666            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     67            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
    6768            continue;
    6869        }
    6970        if (!source->moments->nPixels) {
     71            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
    7072            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    7173            continue;
    … …  
    7476        source->modelEXT = pmSourceModelGuess (source, options->type);
    7577        if (source->modelEXT == NULL) {
     78            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    7679            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    77             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    7880            continue;
    7981        }
    … …  
    9193        // exclude the poor fits
    9294        if (!status) {
     95            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    9396            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    94             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    9597            continue;
    9698        }

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

    r29004 r30621  
    8181        if (source->modelPSF == NULL) {
    8282            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
     83            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
    8384            continue;
    8485        }

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

    r30031 r30621  
    5252// mask values indicate the reason the source was rejected:
    5353
     54// XXX some test code (delete eventually)
     55// bool countMaskedSources(pmPSFtry *try, char *msg) {
     56//
     57//     int N1 = 0;
     58//     int N2 = 0;
     59//     for (int i = 0; i < try->sources->n; i++) {
     60//         pmSource *source = try->sources->data[i];
     61//      // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR));
     62//      if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
     63//          N1 ++;
     64//         }
     65//      if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
     66//          N2 ++;
     67//         }
     68//     }
     69//     fprintf (stderr, "%s : masked: %d or %d of %ld\n", msg, N1, N2, try->sources->n);
     70//     return true;
     71// }
     72
    5473// generate a pmPSFtry with a copy of the test PSF sources
    5574pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
    … …  
    203222            chisq->data.F32[i] = 0.0;
    204223            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
     224            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : no PSF model\n", i, source->peak->x, source->peak->y);
    205225        } else {
    206226            flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];

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

    r29935 r30621  
    6161    psFree(tmp->moments);
    6262    psFree(tmp->diffStats);
    63     psFree(tmp->radial);
     63    psFree(tmp->radialAper);
    6464    psFree(tmp->blends);
    6565    psTrace("psModules.objects", 10, "---- end ----\n");
    … …  
    126126    source->extpars = NULL;
    127127    source->diffStats = NULL;
    128     source->radial = NULL;
     128    source->radialAper = NULL;
     129    source->parent = NULL;
    129130
    130131    source->region = psRegionSet(NAN, NAN, NAN, NAN);
    … …  
    159160affecting the original.  This Copy can be used to allow multiple fit attempts on the same
    160161object.  The pixels, variance, and mask arrays all point to the same original subarrays.  The
    161 peak and moments point at the original values.
     162peak and moments point at the original values.  The models, blends, and XXX are NOT copied
    162163*****************************************************************************/
    163164pmSource *pmSourceCopy(pmSource *in)
    … …  
    167168    }
    168169    pmSource *source = pmSourceAlloc ();
    169 
    170     // keep the original ID so we can find map back to the original
    171     P_PM_SOURCE_SET_ID(source, in->id);
    172170
    173171    // peak has the same values as the original
    … …  
    203201}
    204202
    205 /******************************************************************************
    206 pmSourceCopyData(): this creates a new, duplicate source with the same parameters as the
    207 original (but is actually a new source at the same location)
    208 *****************************************************************************/
    209 pmSource *pmSourceCopyData(pmSource *in)
    210 {
    211     if (in == NULL) {
    212         return(NULL);
    213     }
    214     // this copy is used to allow multiple fit attempts on the
    215     // same object.  the pixels, variance, and mask arrays all point to
    216     // the same original subarrays.  the peak and moments point at
    217     // the original values.
    218     pmSource *source = pmSourceAlloc ();
    219 
    220     // this is actually the same peak as the original, is this the correct way to handle this?
    221     if (in->peak != NULL) {
    222         source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
    223         source->peak->xf = in->peak->xf;
    224         source->peak->yf = in->peak->yf;
    225         source->peak->flux = in->peak->flux;
    226         source->peak->SN = in->peak->SN;
    227     }
    228 
    229     // copy the values in the moments structure
    230     if (in->moments != NULL) {
    231         source->moments  =  pmMomentsAlloc();
    232         *source->moments = *in->moments;
    233     }
    234 
    235     // These images are all views to the parent.
    236     // We want a new view, but pointing at the same pixels.
    237     source->pixels   = psImageCopyView(NULL, in->pixels);
    238     source->variance   = psImageCopyView(NULL, in->variance);
    239     source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
    240 
    241     // the maskObj is a unique mask array; create a new mask image
    242     source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL;
    243 
    244     source->type = in->type;
    245     source->mode = in->mode;
    246     source->imageID = in->imageID;
    247 
    248     return(source);
    249 }
    250 
    251203// x,y are defined in the parent image coords of readout->image
    252204bool pmSourceDefinePixels(pmSource *mySource,
    … …  
    300252    newRegion = psRegionForImage (readout->image, newRegion);
    301253
     254    // re-define if required by region or absence of pixels
    302255    extend = false;
    303256    extend |= (int)(newRegion.x0) < (int)(mySource->region.x0);
    … …  
    335288    }
    336289    return extend;
     290}
     291
     292bool pmSourceRedefinePixelsByRegion(pmSource *mySource,
     293                                    const pmReadout *readout,
     294                                    psRegion newRegion)
     295{
     296    PS_ASSERT_PTR_NON_NULL(mySource, false);
     297    PS_ASSERT_PTR_NON_NULL(readout, false);
     298    PS_ASSERT_PTR_NON_NULL(readout->image, false);
     299
     300    // re-create the subimage
     301    psFree (mySource->pixels);
     302    psFree (mySource->variance);
     303    psFree (mySource->maskView);
     304       
     305    mySource->pixels   = psImageSubset(readout->image,    newRegion);
     306    mySource->variance = psImageSubset(readout->variance, newRegion);
     307    mySource->maskView = psImageSubset(readout->mask,     newRegion);
     308    mySource->region   = newRegion;
     309
     310    // re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
     311    // pixels (eg, with psImageKeepCircle)
     312    mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_IMAGE_MASK);
     313
     314    // drop the old modelFlux pixels and force the user to re-create
     315    psFree (mySource->modelFlux);
     316    mySource->modelFlux = NULL;
     317
     318    // drop the old psfImage pixels and force the user to re-create
     319    psFree (mySource->psfImage);
     320    mySource->psfImage = NULL;
     321
     322    return true;
    337323}
    338324
    … …  
    675661            if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
    676662                source->type = PM_SOURCE_TYPE_STAR;
    677                 source->mode |= PM_SOURCE_MODE_PSFSTAR;
     663                source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
    678664                Npsf ++;
    679665                continue;
    … …  
    934920}
    935921
    936 // construct a realization of the source model
    937 // XXX this function should optionally save an existing psf image from modelFlux
     922// construct a realization of the PSF model at the location of this source
    938923bool pmSourceCachePSF (pmSource *source, psImageMaskType maskVal) {
    939924    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    10201005    }
    10211006
    1022     if (!addNoise) {
    1023         if (add) {
    1024             status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    1025             source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1026         } else {
    1027             status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    1028             source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1029         }
     1007    if (add) {
     1008        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     1009        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1010    } else {
     1011        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     1012        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    10301013    }
    10311014

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

    r29546 r30621  
    3535    PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008,
    3636    PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
     37    PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020,
    3738} pmSourceTmpF;
    3839
    … …  
    4243 *  simplest measurement of a source is the location and flux of the peak pixel
    4344 *  associated with the source:
     45 *
     46 * a pmSource is the information about a (possible) blob of flux in a specific image.  A source
     47 * may represent an insignificant or undetected source.  There may be multiple representations
     48 * of an image (eg, alternate smoothed copies); sources on alternate images may have a pointer
     49 * to the version on the primary image (source->parent).  A set of sources on different, but
     50 * related images (eg, multiple exposures or different filters) which (may) represent the same
     51 * astronomical object are grouped together with the pmPhotObj type (set pmPhotObj.h).
     52 *
     53 * A single source may be fitted with multiple models (not at the same time!).  The PSF model
     54 * fit is a fit of the position (optionally) and the flux to the PSF model at the location of
     55 * the source.  Alternate model fits are extended source models. The best model fit is used to
     56 * subtract the object from the image.
    4457 *
    4558 *  XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits??
    … …  
    5467    pmPeak  *peak;                      ///< Description of peak pixel.
    5568    psImage *pixels;                    ///< Rectangular region including object pixels.
    56     psImage *variance;                    ///< Image variance.
     69    psImage *variance;                  ///< Image variance.
    5770    psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
    5871    psImage *maskView;                  ///< view into global image mask for this object region
    5972    psImage *modelFlux;                 ///< cached copy of the best model for this source
    60     psImage *psfImage;                   ///< cached copy of the psf model for this source
     73    psImage *psfImage;                  ///< cached copy of the psf model for this source
    6174    pmMoments *moments;                 ///< Basic moments measured for the object.
    6275    pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
    … …  
    89102    pmSourceExtendedPars *extpars;      ///< extended source parameters
    90103    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
    91     pmSourceRadialApertures *radial;    ///< radial flux in circular apertures
     104    psArray *radialAper;                ///< radial flux in circular apertures
     105    pmSource *parent;                   ///< reference to the master source from which this is derived
    92106    int imageID;
    93107};
    … …  
    162176);
    163177
     178bool pmSourceRedefinePixelsByRegion (
     179    pmSource *mySource,   ///< source to be re-defined
     180    const pmReadout *readout,   ///< base the source on this readout
     181    psRegion newRegion ///< region for source pixel definition
     182);
     183
    164184/** pmSourcePSFClump()
    165185 *

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

    r28013 r30621  
    1717#endif
    1818
    19 // #include <stdio.h>
    20 // #include <math.h>
    21 // #include <string.h>
    2219#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"
    3620#include "pmSourceExtendedPars.h"
    3721

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

    r29004 r30621  
    8989    if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
    9090
     91    // once we have fitted a model, we need to record that this model is a PCM model:
     92    pcm->modelConv->isPCM = true;
     93
    9194    // models can go insane: reject these
    9295    bool onPic = true;

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

    r29546 r30621  
    5959
    6060// lookup the EXTNAME values used for table data and image header segments
    61 static bool sourceExtensions(psString *headname, // Extension name for header
    62                              psString *dataname, // Extension name for data
    63                              psString *deteffname, // Extension name for detection efficiency
    64                              psString *xsrcname, // Extension name for extended sources
    65                              psString *xfitname, // Extension name for extended fits
    66                              const pmFPAfile *file, // File of interest
    67                              const pmFPAview *view // View to level of interest
    68                              )
     61bool pmSourceIOextnames(psString *headname,    // Extension name for image header
     62                        psString *dataname,    // Extension name for PSF table data
     63                        psString *deteffname,  // Extension name for detection efficiency
     64                        psString *xsrcname,    // Extension name for extended non-parametric measurements
     65                        psString *xfitname,    // Extension name for extended fitted measurements
     66                        psString *xradname,    // Extension name for radial apertures
     67                        const pmFPAfile *file, // File of interest
     68                        const pmFPAview *view  // View to level of interest
     69    )
    6970{
    7071    bool status;                        // Status of MD lookup
    … …  
    8788    }
    8889
    89     // EXTNAME for table data
     90    // EXTNAME for PSF table data
    9091    if (dataname) {
    9192        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
    … …  
    107108    }
    108109
    109     // EXTNAME for extended source data table
     110    // EXTNAME for extended source non-parametric measurements
    110111    if (xsrcname) {
    111112        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
    … …  
    117118    }
    118119
     120    // EXTNAME for extended source fitted measurements
    119121    if (xfitname) {
    120         // EXTNAME for extended source data table
    121122        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
    122123        if (!rule) {
    … …  
    125126        }
    126127        *xfitname = pmFPAfileNameFromRule (rule, file, view);
     128    }
     129
     130    // EXTNAME for radial apertures
     131    if (xradname) {
     132        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XRAD");
     133        if (!rule) {
     134            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XRAD in EXTNAME.RULES in camera.config");
     135            return false;
     136        }
     137        *xradname = pmFPAfileNameFromRule (rule, file, view);
    127138    }
    128139
    … …  
    344355            status &= pmSourcesWrite_##TYPE##_XFIT (file->fits, readout, sources, file->header, xfitname); \
    345356        }                                                               \
     357        if (xradname) {                                                 \
     358            status &= pmSourcesWrite_##TYPE##_XRAD (file->fits, readout, sources, file->header, xradname, recipe); \
     359        }                                                               \
    346360    }
    347361
    … …  
    449463        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
    450464        bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
     465        bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
    451466
    452467        // define the EXTNAME values for the different data segments:
    … …  
    456471        psString xsrcname = NULL;
    457472        psString xfitname = NULL;
    458         if (!sourceExtensions(&headname, &dataname, &deteffname,
    459                               XSRC_OUTPUT ? &xsrcname : NULL,
    460                               XFIT_OUTPUT ? &xfitname : NULL,
    461                               file, view)) {
     473        psString xradname = NULL;
     474        if (!pmSourceIOextnames(&headname, &dataname, &deteffname,
     475                                XSRC_OUTPUT ? &xsrcname : NULL,
     476                                XFIT_OUTPUT ? &xfitname : NULL,
     477                                XRAD_OUTPUT ? &xradname : NULL,
     478                                file, view)) {
    462479            return false;
    463480        }
    … …  
    539556                psMetadataAddStr (outhead, PS_LIST_TAIL, "XFITNAME", PS_META_REPLACE, "name of XFIT table extension", xfitname);
    540557            }
    541    
     558            if (xradname) {
     559                psMetadataAddStr (outhead, PS_LIST_TAIL, "XRADNAME", PS_META_REPLACE, "name of XRAD table extension", xradname);
     560            }
    542561
    543562            // these are case-sensitive since the EXTYPE is case-sensitive
    … …  
    563582        }
    564583
    565 
    566584        // write out the detection efficiency TABLE segments
    567585        if (deteffname) {
    … …  
    583601        psFree (xsrcname);
    584602        psFree (xfitname);
     603        psFree (xradname);
    585604        psFree (deteffname);
    586605
    … …  
    593612        psFree (xsrcname);
    594613        psFree (xfitname);
     614        psFree (xradname);
    595615        psFree (deteffname);
    596616        return false;
    … …  
    939959        psString dataname = NULL;
    940960        psString deteffname = NULL;
    941         if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) {
     961        if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) {
    942962            return false;
    943963        }

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

    r29546 r30621  
    3030bool pmSourcesWrite_SMPDATA_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    3131bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     32bool pmSourcesWrite_SMPDATA_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    3233
    3334bool pmSourcesWrite_PS1_DEV_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    3435bool pmSourcesWrite_PS1_DEV_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    3536bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     37bool pmSourcesWrite_PS1_DEV_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    3638
    3739bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    3840bool pmSourcesWrite_PS1_DEV_1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    3941bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     42bool pmSourcesWrite_PS1_DEV_1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    4043
    4144bool pmSourcesWrite_PS1_CAL_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    4245bool pmSourcesWrite_PS1_CAL_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    4346bool pmSourcesWrite_PS1_CAL_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     47bool pmSourcesWrite_PS1_CAL_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    4448
    4549bool pmSourcesWrite_CMF_PS1_V1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    4650bool pmSourcesWrite_CMF_PS1_V1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    4751bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     52bool pmSourcesWrite_CMF_PS1_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    4853
    4954bool pmSourcesWrite_CMF_PS1_V2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    5055bool pmSourcesWrite_CMF_PS1_V2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    5156bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     57bool pmSourcesWrite_CMF_PS1_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    5258
    5359bool pmSourcesWrite_CMF_PS1_V3(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    5460bool pmSourcesWrite_CMF_PS1_V3_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    5561bool pmSourcesWrite_CMF_PS1_V3_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     62bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    5663
    5764bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    5865bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    5966bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     67bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    6068
    6169bool pmSourcesWrite_CMF_PS1_DV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    6270bool pmSourcesWrite_CMF_PS1_DV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    6371bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     72bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    6473
    6574bool pmSourcesWrite_CMF_PS1_DV2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    6675bool pmSourcesWrite_CMF_PS1_DV2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    6776bool pmSourcesWrite_CMF_PS1_DV2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     77bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    6878
    6979psArray *pmSourcesReadCMP (char *filename, psMetadata *header);

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

    r29004 r30621  
    676676    return true;
    677677}
     678
     679bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     680{
     681    return true;
     682}

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

    r29620 r30621  
    737737    return true;
    738738}
     739
     740bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     741{
     742    return true;
     743}

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

    r29019 r30621  
    2020
    2121#include "pmConfig.h"
     22#include "pmErrorCodes.h"
    2223#include "pmDetrendDB.h"
    2324
    … …  
    7172    psF32 errMag, chisq, apRadius;
    7273    psS32 nPix, nDOF;
    73     char keyword1[80], keyword2[80];
    7474
    7575    pmChip *chip = readout->parent->parent;
    … …  
    103103    table = psArrayAllocEmpty (sources->n);
    104104
     105# if (0)
    105106    // we use this just to define the output vectors (which must be present for all objects)
    106107    bool status = false;
    … …  
    118119      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
    119120    }
    120 
     121# endif
    121122
    122123    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    258259        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    259260
     261# if (0)
     262        // XXX if we have raw radial apertures, write them out here
    260263        psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    261264        psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
    … …  
    284287        psFree (radFluxErr);
    285288        psFree (radFill);
     289# endif
    286290
    287291        // XXX not sure how to get this : need to load Nimages with weight?
    … …  
    781785    return true;
    782786}
     787
     788// **** write out the radial flux values for the sources for a given matched-PSF image
     789// **** how do we distinguish the matched-PSF images from the non-matched version
     790bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     791{
     792
     793    psArray *table;
     794    psMetadata *row;
     795    psF32 xPos, yPos;
     796    char keyword1[80], keyword2[80];
     797
     798    // create a header to hold the output data
     799    psMetadata *outhead = psMetadataAlloc ();
     800
     801    // write the links to the image header
     802    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname);
     803
     804    // we use this just to define the output vectors (which must be present for all objects)
     805    bool status = false;
     806    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     807    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     808    psAssert (radMax, "this must have been defined and tested earlier!");
     809    psAssert (radMax->n, "this must have been defined and tested earlier!");
     810    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     811
     812    // write the radial profile apertures to header
     813    for (int i = 0; i < radMax->n; i++) {
     814      sprintf (keyword1, "RMIN_%02d", i);
     815      sprintf (keyword2, "RMAX_%02d", i);
     816      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     817      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     818    }
     819
     820    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     821    if (!fwhmValues) {
     822        psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
     823        return false;
     824    }
     825
     826    // let's write these out in S/N order
     827    sources = psArraySort (sources, pmSourceSortBySN);
     828
     829    table = psArrayAllocEmpty (sources->n);
     830
     831    // we write out all sources, regardless of quality.  the source flags tell us the state
     832    for (int i = 0; i < sources->n; i++) {
     833
     834        pmSource *source = sources->data[i];
     835
     836        // skip sources without radial aper measurements (or insufficient)
     837        if (source->radialAper == NULL) continue;
     838        psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
     839
     840        for (int entry = 0; entry < fwhmValues->n; entry++) {
     841
     842            // choose the convolved EXT model, if available, otherwise the simple one
     843            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
     844            assert (radialAper);
     845
     846            bool useMoments = true;
     847            useMoments = (useMoments && source->moments);          // can't if there are no moments
     848            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
     849            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
     850
     851            if (useMoments) {
     852                xPos = source->moments->Mx;
     853                yPos = source->moments->My;
     854            } else {
     855                xPos = source->peak->xf;
     856                yPos = source->peak->yf;
     857            }
     858
     859            row = psMetadataAlloc ();
     860
     861            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     862            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     863            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
     864            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
     865            psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
     866
     867            // XXX if we have raw radial apertures, write them out here
     868            psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
     869            psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
     870            psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
     871            psVectorInit (radFlux,    NAN);
     872            psVectorInit (radFluxErr, NAN);
     873            psVectorInit (radFill,    NAN);
     874            if (!radialAper->flux) goto write_annuli;
     875            if (!radialAper->fill) goto write_annuli;
     876            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
     877            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
     878
     879            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     880            for (int j = 0; j < radialAper->flux->n; j++) {
     881                radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
     882                radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
     883                radFill->data.F32[j]    = radialAper->fill->data.F32[j];
     884            }
     885
     886        write_annuli:
     887            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
     888            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
     889            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
     890            psFree (radFlux);
     891            psFree (radFluxErr);
     892            psFree (radFill);
     893
     894            psArrayAdd (table, 100, row);
     895            psFree (row);
     896        }
     897    }
     898
     899    if (table->n == 0) {
     900        if (!psFitsWriteBlank (fits, outhead, extname)) {
     901            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     902            psFree(outhead);
     903            psFree(table);
     904            return false;
     905        }
     906        psFree (outhead);
     907        psFree (table);
     908        return true;
     909    }
     910
     911    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     912    if (!psFitsWriteTable (fits, outhead, table, extname)) {
     913        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     914        psFree (outhead);
     915        psFree(table);
     916        return false;
     917    }
     918    psFree (outhead);
     919    psFree (table);
     920    return true;
     921}

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

    r29004 r30621  
    721721    return false;
    722722}
     723
     724bool pmSourcesWrite_CMF_PS1_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     725{
     726    return true;
     727}

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

    r29004 r30621  
    724724    return true;
    725725}
     726
     727bool pmSourcesWrite_CMF_PS1_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     728{
     729    return true;
     730}

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

    r29620 r30621  
    762762    return true;
    763763}
     764
     765bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     766{
     767    return true;
     768}

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

    r29004 r30621  
    702702    return true;
    703703}
     704
     705bool pmSourcesWrite_PS1_CAL_0_XRAD (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     706{
     707    return true;
     708}

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

    r29004 r30621  
    248248    return true;
    249249}
     250
     251bool pmSourcesWrite_PS1_DEV_0_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     252{
     253    return true;
     254}

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

    r29004 r30621  
    588588    return true;
    589589}
     590
     591bool pmSourcesWrite_PS1_DEV_1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     592{
     593    return true;
     594}

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

    r29004 r30621  
    222222    return true;
    223223}
     224
     225bool pmSourcesWrite_SMPDATA_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     226{
     227    return true;
     228}
Note: See TracChangeset for help on using the changeset viewer.