IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32407


Ignore:
Timestamp:
Sep 15, 2011, 4:54:20 PM (15 years ago)
Author:
bills
Message:

use psFitsTableWriteAllCollums to avoid memory explosions

Location:
tags/ipp-20110622/ppTranslate
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • tags/ipp-20110622/ppTranslate

  • tags/ipp-20110622/ppTranslate/src

    • Property svn:mergeinfo deleted

  • tags/ipp-20110622/ppTranslate/src/ppMops.c

    r32189 r32407  
    33
    44#include "ppMops.h"
     5
     6void test()
     7{
     8    psMetadata *md = psMetadataAlloc();
     9
     10    psVector *vec = psVectorAlloc(42, PS_TYPE_S32);
     11
     12    psMetadataAddVector(md, PS_LIST_TAIL, "TEST", 0, NULL, vec);
     13
     14    psFree(vec);
     15    psFree(md);
     16
     17    psLibFinalize();
     18
     19    fprintf (stderr, "found %d leaks at %s\n",
     20        psMemCheckLeaks2 (0,
     21                NULL, stdout, false, 500), "ppMops");
     22
     23    exit(0);
     24}
    525
    626/*
    … …  
    3959    psLibInit(NULL);
    4060
     61    // test();
     62
    4163    ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments
    4264    if (!args) {
    … …  
    4567    }
    4668
     69
    4770    psArray *detections = ppMopsRead(args); // Detections from each input
    4871    if (!detections) {
    … …  
    5174    }
    5275
    53     ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections
    54     psFree(detections);
    55     if (!merged) {
     76    if (!ppMopsPurgeDuplicates(detections)) {
    5677        psErrorStackPrint(stderr, "Unable to merge detections");
    5778        exit(PS_EXIT_SYS_ERROR);
    5879    }
    5980
    60     if (!ppMopsWrite(merged, args)) {
     81    if (!ppMopsWrite(detections, args)) {
    6182        psErrorStackPrint(stderr, "Unable to write detections");
    6283        exit(PS_EXIT_SYS_ERROR);
    6384    }
    6485
    65     psFree(merged);
     86    for (int i = 0; i < detections->n; i++) {
     87        psFree(detections->data[i]);
     88    }
     89    psFree(detections);
    6690    psFree(args);
    67     psFree(detections);
    6891
    6992    psLibFinalize();

  • tags/ipp-20110622/ppTranslate/src/ppMops.h

    r29567 r32407  
    6262
    6363typedef struct {
     64  psString component;                 // skycell_id for these detections
    6465  psString raBoresight, decBoresight; // RA,Dec of telescope boresight
    6566  psString filter;                    // Filter for exposure
    … …  
    7071  double mjd;                         // Modified Julian Date
    7172  float seeing;                       // Seeing of exposure
     73  int   naxis1, naxis2;               // size of the image
    7274  long num;                           // Number of detections
     75  long numGood;                       // Number of "good" detections
     76  psS64 diffSkyfileId;                // unique id for input skyfile
     77  psMetadata *table;                  // Columns from the input file
    7378  psVector *x, *y;                    // Image coordinates
    7479  psVector *ra, *dec;                 // Sky coordinates
    7580  psVector *raErr, *decErr;           // Error in sky coordinates
    76   psVector *mag, *magErr;             // Magnitude and associated error
    77   psVector *chi2, *dof;               // Chi^2 from fitting, with associated degrees of freedom
    78   psVector *cr, *extended;            // Measures of CR-ness and extendedness
    79   psVector *psfMajor, *psfMinor, *psfTheta; // PSF major and minor axes, and position angle
    80   psVector *quality, *numPix;               // PSF quality factor and number of pixels
    81   psVector *xxMoment, *xyMoment, *yyMoment; // Moments
    82   psVector *flags;                    // psphot flags
    83   psVector *diffSkyfileId;            // Identifier for source image
    84   psVector *naxis1, *naxis2;          // Size of image
    8581  psVector *mask;                     // Mask for detections
    86   psVector *nPos;                     // Number of positive pixels
    87   psVector *fPos;                     // Fraction of positive flux
    88   psVector *nRatioBad;                // Fraction of positive pixels to negative
    89   psVector *nRatioMask;               // Fraction of positive pixels to masked
    90   psVector *nRatioAll;                // Fraction of positive pixels to all
    91   psVector *psfInstFlux;              // PSF fit instrumental magnitude
    92   psVector *psfInstFluxSig;           // Sigma of PSF instrumental magnitude
    93   psVector *apMag;                    // Magnitude in standard aperture
    94   psVector *apMagRadius;              // Radius used for aperture mags
    95   psVector *apMagRaw;                 // Magnitude in real aperture
    96   psVector *apFlux;                   // Instrumental flux in standard aperture
    97   psVector *apFluxSig;                // Aperture flux error
    98   psVector *peakFluxAsMag;            // Peak flux expressed as magnitude
    99   psVector *calPsfMag;                // PSF Magnitude using supplied calibration
    100   psVector *calPsfMagSig;             // Measured scatter of zero point calibration
    101   psVector *sky;                      // Sky level
    102   psVector *skySig;                   // Sigma of sky level
    103   psVector *qualityPerfect;           // PSF coverage/quality factor (poor)
    104   psVector *momentsR1;                // First radial moment
    105   psVector *momentsRH;                // Half radial moment
    106   psVector *kronFlux;                 // Kron Flux (in 2.5 R1)
    107   psVector *kronFluxErr;              // Kron Flux Error
    108   psVector *kronFluxInner;            // Kron Flux (in 1.0 R1)
    109   psVector *kronFluxOuter;            // Kron Flux (in 4.0 R1)
    110   psVector *diffRP;                   // Distance to positive match source
    111   psVector *diffSnP;                  // Signal-to-noise of pos match src
    112   psVector *diffRM;                   // Distance to negative match source
    113   psVector *diffSnM;                  // Signal-to-noise of neg match src
    114   psVector *flags2;                   // psphot flags (group 2)
    115   psVector *ippIdet;                  // IPP detection identifier index
    116   psVector *nFrames;                  // Number of frames overlapping source center
    117   psVector *padding;                  // Padding
    11882} ppMopsDetections;
    11983
    120 ppMopsDetections *ppMopsDetectionsAlloc(long num);
     84ppMopsDetections *ppMopsDetectionsAlloc();
    12185
    12286/// Copy a detection
    … …  
    13094
    13195/// Merge detections
    132 ppMopsDetections *ppMopsMerge(const psArray *detections);
     96// ppMopsDetections *ppMopsMerge(const psArray *detections);
     97bool ppMopsPurgeDuplicates(const psArray *detections);
    13398
    13499/// Write detections
    135 bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args);
     100bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args);
    136101
    137102/// Get the version contained in EXTTYPE of the "SkyChip.psf" table:

  • tags/ipp-20110622/ppTranslate/src/ppMopsDetections.c

    r30660 r32407  
    88static void mopsDetectionsFree(ppMopsDetections *det)
    99{
     10    psFree(det->component);
    1011    psFree(det->raBoresight);
    1112    psFree(det->decBoresight);
    1213    psFree(det->filter);
     14    psFree(det->table);
    1315    psFree(det->x);
    1416    psFree(det->y);
    … …  
    1719    psFree(det->raErr);
    1820    psFree(det->decErr);
    19     psFree(det->mag);
    20     psFree(det->magErr);
    21     psFree(det->chi2);
    22     psFree(det->dof);
    23     psFree(det->cr);
    24     psFree(det->extended);
    25     psFree(det->psfMajor);
    26     psFree(det->psfMinor);
    27     psFree(det->psfTheta);
    28     psFree(det->quality);
    29     psFree(det->numPix);
    30     psFree(det->xxMoment);
    31     psFree(det->xyMoment);
    32     psFree(det->yyMoment);
    33     psFree(det->flags);
    34     psFree(det->diffSkyfileId);
    35     psFree(det->naxis1);
    36     psFree(det->naxis2);
    3721    psFree(det->mask);
    38     psFree(det->nPos);
    39     psFree(det->fPos);
    40     psFree(det->nRatioBad);
    41     psFree(det->nRatioMask);
    42     psFree(det->nRatioAll);
    43     psFree(det->psfInstFlux);
    44     psFree(det->psfInstFluxSig);
    45     psFree(det->apMag);
    46     psFree(det->apMagRadius);
    47     psFree(det->apMagRaw);
    48     psFree(det->apFlux);
    49     psFree(det->apFluxSig);
    50     psFree(det->peakFluxAsMag);
    51     psFree(det->calPsfMag);
    52     psFree(det->calPsfMagSig);
    53     psFree(det->sky);
    54     psFree(det->skySig);
    55     psFree(det->qualityPerfect);
    56     psFree(det->momentsR1);
    57     psFree(det->momentsRH);
    58     psFree(det->kronFlux);
    59     psFree(det->kronFluxErr);
    60     psFree(det->kronFluxInner);
    61     psFree(det->kronFluxOuter);
    62     psFree(det->diffRP);
    63     psFree(det->diffSnP);
    64     psFree(det->diffRM);
    65     psFree(det->diffSnM);
    66     psFree(det->flags2);
    67     psFree(det->ippIdet);
    68     psFree(det->nFrames);
    69     psFree(det->padding);
    7022    return;
    7123}
    7224
    73 ppMopsDetections *ppMopsDetectionsAlloc(long num)
     25ppMopsDetections *ppMopsDetectionsAlloc()
    7426{
    7527    ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return
    7628    psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree);
     29    det->component = NULL;
    7730    det->raBoresight = NULL;
    7831    det->decBoresight = NULL;
    … …  
    8639    det->seeing = NAN;
    8740    det->num = 0;
    88     det->x = psVectorAllocEmpty(num, PS_TYPE_F32);
    89     det->y = psVectorAllocEmpty(num, PS_TYPE_F32);
    90     det->ra = psVectorAllocEmpty(num, PS_TYPE_F64);
    91     det->dec = psVectorAllocEmpty(num, PS_TYPE_F64);
    92     det->raErr = psVectorAllocEmpty(num, PS_TYPE_F64);
    93     det->decErr = psVectorAllocEmpty(num, PS_TYPE_F64);
    94     det->mag = psVectorAllocEmpty(num, PS_TYPE_F32);
    95     det->magErr = psVectorAllocEmpty(num, PS_TYPE_F32);
    96     det->chi2 = psVectorAllocEmpty(num, PS_TYPE_F32);
    97     det->dof = psVectorAllocEmpty(num, PS_TYPE_S32);
    98     det->cr = psVectorAllocEmpty(num, PS_TYPE_F32);
    99     det->extended = psVectorAllocEmpty(num, PS_TYPE_F32);
    100     det->psfMajor = psVectorAllocEmpty(num, PS_TYPE_F32);
    101     det->psfMinor = psVectorAllocEmpty(num, PS_TYPE_F32);
    102     det->psfTheta = psVectorAllocEmpty(num, PS_TYPE_F32);
    103     det->quality = psVectorAllocEmpty(num, PS_TYPE_F32);
    104     det->numPix = psVectorAllocEmpty(num, PS_TYPE_S32);
    105     det->xxMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
    106     det->xyMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
    107     det->yyMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
    108     det->flags = psVectorAllocEmpty(num, PS_TYPE_U32);
    109     det->diffSkyfileId = psVectorAllocEmpty(num, PS_TYPE_S64);
    110     det->naxis1 = psVectorAllocEmpty(num, PS_TYPE_S32);
    111     det->naxis2 = psVectorAllocEmpty(num, PS_TYPE_S32);
    112     det->mask = psVectorAllocEmpty(num, PS_TYPE_U8);
    113     det->nPos = psVectorAllocEmpty(num, PS_TYPE_S32);
    114     det->fPos = psVectorAllocEmpty(num, PS_TYPE_F32);
    115     det->nRatioBad = psVectorAllocEmpty(num, PS_TYPE_F32);
    116     det->nRatioMask = psVectorAllocEmpty(num, PS_TYPE_F32);
    117     det->nRatioAll = psVectorAllocEmpty(num, PS_TYPE_F32);
    118     det->psfInstFlux = psVectorAllocEmpty(num, PS_TYPE_F32);
    119     det->psfInstFluxSig = psVectorAllocEmpty(num, PS_TYPE_F32);
    120     det->apMag = psVectorAllocEmpty(num, PS_TYPE_F32);
    121     det->apMagRadius  = psVectorAllocEmpty(num, PS_TYPE_F32);
    122     det->apMagRaw = psVectorAllocEmpty(num, PS_TYPE_F32);
    123     det->apFlux = psVectorAllocEmpty(num, PS_TYPE_F32);
    124     det->apFluxSig = psVectorAllocEmpty(num, PS_TYPE_F32);
    125     det->peakFluxAsMag = psVectorAllocEmpty(num, PS_TYPE_F32);
    126     det->calPsfMag = psVectorAllocEmpty(num, PS_TYPE_F32);
    127     det->calPsfMagSig = psVectorAllocEmpty(num, PS_TYPE_F32);
    128     det->sky = psVectorAllocEmpty(num, PS_TYPE_F32);
    129     det->skySig = psVectorAllocEmpty(num, PS_TYPE_F32);
    130     det->qualityPerfect = psVectorAllocEmpty(num, PS_TYPE_F32);
    131     det->momentsR1 = psVectorAllocEmpty(num, PS_TYPE_F32);
    132     det->momentsRH = psVectorAllocEmpty(num, PS_TYPE_F32);
    133     det->kronFlux = psVectorAllocEmpty(num, PS_TYPE_F32);
    134     det->kronFluxErr = psVectorAllocEmpty(num, PS_TYPE_F32);
    135     det->kronFluxInner = psVectorAllocEmpty(num, PS_TYPE_F32);
    136     det->kronFluxOuter = psVectorAllocEmpty(num, PS_TYPE_F32);
    137     det->diffRP = psVectorAllocEmpty(num, PS_TYPE_F32);
    138     det->diffSnP = psVectorAllocEmpty(num, PS_TYPE_F32);
    139     det->diffRM = psVectorAllocEmpty(num, PS_TYPE_F32);
    140     det->diffSnM = psVectorAllocEmpty(num, PS_TYPE_F32);
    141     det->flags2 = psVectorAllocEmpty(num, PS_TYPE_U32);
    142     det->ippIdet = psVectorAllocEmpty(num, PS_TYPE_U32);
    143     det->nFrames = psVectorAllocEmpty(num, PS_TYPE_U16);
    144     det->padding = psVectorAllocEmpty(num, PS_TYPE_S16);
     41    det->table = NULL;
     42    det->x = NULL;
     43    det->y = NULL;
     44    det->ra = NULL;
     45    det->dec = NULL;
     46    det->raErr = NULL;
     47    det->decErr = NULL;
     48    det->mask = NULL;
     49    det->diffSkyfileId = 0;
    14550    return det;
    14651}
    147 
    148 ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num)
    149 {
    150     det->x = psVectorRealloc(det->x, num);
    151     det->y = psVectorRealloc(det->y, num);
    152     det->ra = psVectorRealloc(det->ra, num);
    153     det->dec = psVectorRealloc(det->dec, num);
    154     det->raErr = psVectorRealloc(det->raErr, num);
    155     det->decErr = psVectorRealloc(det->decErr, num);
    156     det->mag = psVectorRealloc(det->mag, num);
    157     det->magErr = psVectorRealloc(det->magErr, num);
    158     det->chi2 = psVectorRealloc(det->chi2, num);
    159     det->dof = psVectorRealloc(det->dof, num);
    160     det->cr = psVectorRealloc(det->cr, num);
    161     det->extended = psVectorRealloc(det->extended, num);
    162     det->psfMajor = psVectorRealloc(det->psfMajor, num);
    163     det->psfMinor = psVectorRealloc(det->psfMinor, num);
    164     det->psfTheta = psVectorRealloc(det->psfTheta, num);
    165     det->quality = psVectorRealloc(det->quality, num);
    166     det->numPix = psVectorRealloc(det->numPix, num);
    167     det->xxMoment = psVectorRealloc(det->xxMoment, num);
    168     det->xyMoment = psVectorRealloc(det->xyMoment, num);
    169     det->yyMoment = psVectorRealloc(det->yyMoment, num);
    170     det->flags = psVectorRealloc(det->flags, num);
    171     det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num);
    172     det->naxis1 = psVectorRealloc(det->naxis1, num);
    173     det->naxis2 = psVectorRealloc(det->naxis2, num);
    174     det->mask = psVectorRealloc(det->mask, num);
    175     det->nPos = psVectorRealloc(det->nPos, num);
    176     det->fPos = psVectorRealloc(det->fPos, num);
    177     det->nRatioBad = psVectorRealloc(det->nRatioBad, num);
    178     det->nRatioMask = psVectorRealloc(det->nRatioMask, num);
    179     det->nRatioAll = psVectorRealloc(det->nRatioAll, num);
    180     det->psfInstFlux = psVectorRealloc(det->psfInstFlux, num);
    181     det->psfInstFluxSig = psVectorRealloc(det->psfInstFluxSig, num);
    182     det->apMag = psVectorRealloc(det->apMag, num);
    183     det->apMagRadius = psVectorRealloc(det->apMagRadius, num);
    184     det->apMagRaw = psVectorRealloc(det->apMagRadius, num);
    185     det->apFlux = psVectorRealloc(det->apFlux, num);
    186     det->apFluxSig = psVectorRealloc(det->apFluxSig, num);
    187     det->peakFluxAsMag = psVectorRealloc(det->peakFluxAsMag, num);
    188     det->calPsfMag = psVectorRealloc(det->calPsfMag, num);
    189     det->calPsfMagSig = psVectorRealloc(det->calPsfMagSig, num);
    190     det->sky = psVectorRealloc(det->sky, num);
    191     det->skySig = psVectorRealloc(det->skySig, num);
    192     det->qualityPerfect = psVectorRealloc(det->qualityPerfect, num);
    193     det->momentsR1 = psVectorRealloc(det->momentsR1, num);
    194     det->momentsRH = psVectorRealloc(det->momentsRH, num);
    195     det->kronFlux = psVectorRealloc(det->kronFlux, num);
    196     det->kronFluxErr = psVectorRealloc(det->kronFluxErr, num);
    197     det->kronFluxInner = psVectorRealloc(det->kronFluxInner, num);
    198     det->kronFluxOuter = psVectorRealloc(det->kronFluxOuter, num);
    199     det->diffRP = psVectorRealloc(det->diffRP, num);
    200     det->diffSnP = psVectorRealloc(det->diffSnP, num);
    201     det->diffRM = psVectorRealloc(det->diffRM, num);
    202     det->diffSnM = psVectorRealloc(det->diffSnM, num);
    203     det->flags2 = psVectorRealloc(det->flags2, num);
    204     det->ippIdet = psVectorRealloc(det->ippIdet, num);
    205     det->nFrames = psVectorRealloc(det->nFrames, num);
    206     det->padding = psVectorRealloc(det->padding, num);
    207     return det;
    208 }
    209 
    210 bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec,
    211                          double raErr, double decErr, float mag, float magErr,
    212                          float chi2, int dof, float cr, float extended, float psfMajor,
    213                          float psfMinor, float psfTheta, float quality, int numPix,
    214                          float xxMoment, float xyMoment, float yyMoment,
    215                          psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2,
    216                          int nPos, float fPos, float nRatioBad, float nRatioMask, float nRatioAll,
    217                          float psfInstFlux, float psfInstFluxSig,
    218                          float apMag, float apMagRadius, float apMagRaw, float apFlux, float apFluxSig,
    219                          float peakFluxAsMag, float calPsfMag, float calPsfMagSig,
    220                          float sky, float skySig, float qualityPerfect,
    221                          float momentsR1, float momentsRH,
    222                          float kronFlux, float kronFluxErr, float kronFluxInner, float kronFluxOuter,
    223                          float diffRP, float diffSnP, float diffRM, float diffSnM,
    224                          psU32 flags2, psU32 ippIdet, psU16 nFrames, psS16 padding)
    225 {
    226     psVectorAppend(det->x, x);
    227     psVectorAppend(det->y, y);
    228     psVectorAppend(det->ra, ra);
    229     psVectorAppend(det->dec, dec);
    230     psVectorAppend(det->raErr, raErr);
    231     psVectorAppend(det->decErr, decErr);
    232     psVectorAppend(det->mag, mag);
    233     psVectorAppend(det->magErr, magErr);
    234     psVectorAppend(det->chi2, chi2);
    235     psVectorAppend(det->dof, dof);
    236     psVectorAppend(det->cr, cr);
    237     psVectorAppend(det->extended, extended);
    238     psVectorAppend(det->psfMajor, psfMajor);
    239     psVectorAppend(det->psfMinor, psfMinor);
    240     psVectorAppend(det->psfTheta, psfTheta);
    241     psVectorAppend(det->quality, quality);
    242     psVectorAppend(det->numPix, numPix);
    243     psVectorAppend(det->xxMoment, xxMoment);
    244     psVectorAppend(det->xyMoment, xyMoment);
    245     psVectorAppend(det->yyMoment, yyMoment);
    246     psVectorAppend(det->flags, flags);
    247     psVectorAppend(det->diffSkyfileId, diffSkyfileId);
    248     psVectorAppend(det->naxis1, naxis1);
    249     psVectorAppend(det->naxis2, naxis2);
    250     psVectorAppend(det->mask, 0);
    251     psVectorAppend(det->nPos, nPos);
    252     psVectorAppend(det->fPos, fPos);
    253     psVectorAppend(det->nRatioBad, nRatioBad);
    254     psVectorAppend(det->nRatioMask, nRatioMask);
    255     psVectorAppend(det->nRatioAll, nRatioAll);
    256     psVectorAppend(det->psfInstFlux, psfInstFlux);
    257     psVectorAppend(det->psfInstFluxSig, psfInstFluxSig);
    258     psVectorAppend(det->apMag, apMag);
    259     psVectorAppend(det->apMagRadius, apMagRadius);
    260     psVectorAppend(det->apMagRaw, apMagRaw);
    261     psVectorAppend(det->apFlux, apFlux);
    262     psVectorAppend(det->apFluxSig, apFluxSig);
    263     psVectorAppend(det->peakFluxAsMag, peakFluxAsMag);
    264     psVectorAppend(det->calPsfMag, calPsfMag);
    265     psVectorAppend(det->calPsfMagSig, calPsfMagSig);
    266     psVectorAppend(det->sky, sky);
    267     psVectorAppend(det->skySig, skySig);
    268     psVectorAppend(det->qualityPerfect, qualityPerfect);
    269     psVectorAppend(det->momentsR1, momentsR1);
    270     psVectorAppend(det->momentsRH, momentsRH);
    271     psVectorAppend(det->kronFlux, kronFlux);
    272     psVectorAppend(det->kronFluxErr, kronFluxErr);
    273     psVectorAppend(det->kronFluxInner, kronFluxInner);
    274     psVectorAppend(det->kronFluxOuter, kronFluxOuter);
    275     psVectorAppend(det->diffRP, diffRP);
    276     psVectorAppend(det->diffSnP, diffSnP);
    277     psVectorAppend(det->diffRM, diffRM);
    278     psVectorAppend(det->diffSnM, diffSnM);
    279     psVectorAppend(det->flags2, flags2);
    280     psVectorAppend(det->ippIdet, ippIdet);
    281     psVectorAppend(det->nFrames, nFrames);
    282     psVectorAppend(det->padding, padding);
    283     return true;
    284 }
    285 
    286 bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index)
    287 {
    288     psVectorAppend(target->x, source->x->data.F32[index]);
    289     psVectorAppend(target->y, source->y->data.F32[index]);
    290     psVectorAppend(target->ra, source->ra->data.F64[index]);
    291     psVectorAppend(target->dec, source->dec->data.F64[index]);
    292     psVectorAppend(target->raErr, source->raErr->data.F64[index]);
    293     psVectorAppend(target->decErr, source->decErr->data.F64[index]);
    294     psVectorAppend(target->mag, source->mag->data.F32[index]);
    295     psVectorAppend(target->magErr, source->magErr->data.F32[index]);
    296     psVectorAppend(target->chi2, source->chi2->data.F32[index]);
    297     psVectorAppend(target->dof, source->dof->data.S32[index]);
    298     psVectorAppend(target->cr, source->cr->data.F32[index]);
    299     psVectorAppend(target->extended, source->extended->data.F32[index]);
    300     psVectorAppend(target->psfMajor, source->psfMajor->data.F32[index]);
    301     psVectorAppend(target->psfMinor, source->psfMinor->data.F32[index]);
    302     psVectorAppend(target->psfTheta, source->psfTheta->data.F32[index]);
    303     psVectorAppend(target->quality, source->quality->data.F32[index]);
    304     psVectorAppend(target->numPix, source->numPix->data.S32[index]);
    305     psVectorAppend(target->xxMoment, source->xxMoment->data.F32[index]);
    306     psVectorAppend(target->xyMoment, source->xyMoment->data.F32[index]);
    307     psVectorAppend(target->yyMoment, source->yyMoment->data.F32[index]);
    308     psVectorAppend(target->flags, source->flags->data.U32[index]);
    309     psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]);
    310     psVectorAppend(target->naxis1, source->naxis1->data.S32[index]);
    311     psVectorAppend(target->naxis2, source->naxis2->data.S32[index]);
    312     psVectorAppend(target->mask, 0);
    313     psVectorAppend(target->nPos, source->nPos->data.S32[index]);
    314     psVectorAppend(target->fPos, source->fPos->data.F32[index]);
    315     psVectorAppend(target->nRatioBad, source->nRatioBad->data.F32[index]);
    316     psVectorAppend(target->nRatioMask, source->nRatioMask->data.F32[index]);
    317     psVectorAppend(target->nRatioAll, source->nRatioAll->data.F32[index]);
    318     psVectorAppend(target->psfInstFlux, source->psfInstFlux->data.F32[index]);
    319     psVectorAppend(target->psfInstFluxSig, source->psfInstFluxSig->data.F32[index]);
    320     psVectorAppend(target->apMag, source->apMag->data.F32[index]);
    321     psVectorAppend(target->apMagRadius, source->apMagRadius->data.F32[index]);
    322     psVectorAppend(target->apMagRaw, source->apMagRaw->data.F32[index]);
    323     psVectorAppend(target->apFlux, source->apFlux->data.F32[index]);
    324     psVectorAppend(target->apFluxSig, source->apFluxSig->data.F32[index]);
    325     psVectorAppend(target->peakFluxAsMag, source->peakFluxAsMag->data.F32[index]);
    326     psVectorAppend(target->calPsfMag, source->calPsfMag->data.F32[index]);
    327     psVectorAppend(target->calPsfMagSig, source->calPsfMagSig->data.F32[index]);
    328     psVectorAppend(target->sky, source->sky->data.F32[index]);
    329     psVectorAppend(target->skySig, source->skySig->data.F32[index]);
    330     psVectorAppend(target->qualityPerfect, source->qualityPerfect->data.F32[index]);
    331     psVectorAppend(target->momentsR1, source->momentsR1->data.F32[index]);
    332     psVectorAppend(target->momentsRH, source->momentsRH->data.F32[index]);
    333     psVectorAppend(target->kronFlux, source->kronFlux->data.F32[index]);
    334     psVectorAppend(target->kronFluxErr, source->kronFluxErr->data.F32[index]);
    335     psVectorAppend(target->kronFluxInner, source->kronFluxInner->data.F32[index]);
    336     psVectorAppend(target->kronFluxOuter, source->kronFluxOuter->data.F32[index]);
    337     psVectorAppend(target->diffRP, source->diffRP->data.F32[index]);
    338     psVectorAppend(target->diffSnP, source->diffSnP->data.F32[index]);
    339     psVectorAppend(target->diffRM, source->diffRM->data.F32[index]);
    340     psVectorAppend(target->diffSnM, source->diffSnM->data.F32[index]);
    341     psVectorAppend(target->flags2, source->flags2->data.U32[index]);
    342     psVectorAppend(target->ippIdet, source->ippIdet->data.U32[index]);
    343     psVectorAppend(target->nFrames, source->nFrames->data.U16[index]);
    344     psVectorAppend(target->padding, source->padding->data.S16[index]);
    345 
    346     target->num++;
    347 
    348     return true;
    349 }
    350 
    351 
    352 bool ppMopsDetectionsPurge(ppMopsDetections *det)
    353 {
    354     long num = 0;
    355     for (long i = 0; i < det->num; i++) {
    356         if (!det->mask->data.U8[i]) {
    357             if (i == num) {
    358                 // No need to copy
    359                 num++;
    360                 continue;
    361             }
    362             det->x->data.F32[num] = det->x->data.F32[i];
    363             det->y->data.F32[num] = det->y->data.F32[i];
    364             det->ra->data.F64[num] = det->ra->data.F64[i];
    365             det->dec->data.F64[num] = det->dec->data.F64[i];
    366             det->raErr->data.F64[num] = det->raErr->data.F64[i];
    367             det->decErr->data.F64[num] = det->decErr->data.F64[i];
    368             det->mag->data.F32[num] = det->mag->data.F32[i];
    369             det->magErr->data.F32[num] = det->magErr->data.F32[i];
    370             det->chi2->data.F32[num] = det->chi2->data.F32[i];
    371             det->dof->data.S32[num] = det->dof->data.S32[i];
    372             det->cr->data.F32[num] = det->cr->data.F32[i];
    373             det->extended->data.F32[num] = det->extended->data.F32[i];
    374             det->psfMajor->data.F32[num] = det->psfMajor->data.F32[i];
    375             det->psfMinor->data.F32[num] = det->psfMinor->data.F32[i];
    376             det->psfTheta->data.F32[num] = det->psfTheta->data.F32[i];
    377             det->quality->data.F32[num] = det->quality->data.F32[i];
    378             det->numPix->data.S32[num] = det->numPix->data.S32[i];
    379             det->xxMoment->data.F32[num] = det->xxMoment->data.F32[i];
    380             det->xyMoment->data.F32[num] = det->xyMoment->data.F32[i];
    381             det->yyMoment->data.F32[num] = det->yyMoment->data.F32[i];
    382             det->flags->data.U32[num] = det->flags->data.U32[i];
    383             det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i];
    384             det->naxis1->data.S32[num] = det->naxis1->data.S32[i];
    385             det->naxis2->data.S32[num] = det->naxis2->data.S32[i];
    386             det->mask->data.U8[num] = 0;
    387             det->nPos->data.S32[num] = det->nPos->data.S32[i];
    388             det->fPos->data.F32[num] = det->fPos->data.F32[i];
    389             det->nRatioBad->data.F32[num] = det->nRatioBad->data.F32[i];
    390             det->nRatioMask->data.F32[num] = det->nRatioMask->data.F32[i];
    391             det->nRatioAll->data.F32[num] = det->nRatioAll->data.F32[i];
    392             det->psfInstFlux->data.F32[num] = det->psfInstFlux->data.F32[i];
    393             det->psfInstFluxSig->data.F32[num] = det->psfInstFluxSig->data.F32[i];
    394             det->apMag->data.F32[num] = det->apMag->data.F32[i];
    395             det->apMagRadius->data.F32[num] = det->apMagRadius->data.F32[i];
    396             det->apMagRaw->data.F32[num] = det->apMagRaw->data.F32[i];
    397             det->apFlux->data.F32[num] = det->apFlux->data.F32[i];
    398             det->apFluxSig->data.F32[num] = det->apFluxSig->data.F32[i];
    399             det->peakFluxAsMag->data.F32[num] = det->peakFluxAsMag->data.F32[i];
    400             det->calPsfMag->data.F32[num] = det->calPsfMag->data.F32[i];
    401             det->calPsfMagSig->data.F32[num] = det->calPsfMagSig->data.F32[i];
    402             det->sky->data.F32[num] = det->sky->data.F32[i];
    403             det->skySig->data.F32[num] = det->skySig->data.F32[i];
    404             det->qualityPerfect->data.F32[num] = det->qualityPerfect->data.F32[i];
    405             det->momentsR1->data.F32[num] = det->momentsR1->data.F32[i];
    406             det->momentsRH->data.F32[num] = det->momentsRH->data.F32[i];
    407             det->kronFlux->data.F32[num] = det->kronFlux->data.F32[i];
    408             det->kronFluxErr->data.F32[num] = det->kronFluxErr->data.F32[i];
    409             det->kronFluxInner->data.F32[num] = det->kronFluxInner->data.F32[i];
    410             det->kronFluxOuter->data.F32[num] = det->kronFluxOuter->data.F32[i];
    411             det->diffRP->data.F32[num] = det->diffRP->data.F32[i];
    412             det->diffSnP->data.F32[num] = det->diffSnP->data.F32[i];
    413             det->diffRM->data.F32[num] = det->diffRM->data.F32[i];
    414             det->diffSnM->data.F32[num] = det->diffSnM->data.F32[i];
    415             det->flags2->data.U32[num] = det->flags2->data.U32[i];
    416             det->ippIdet->data.U32[num] = det->ippIdet->data.U32[i];
    417             det->nFrames->data.U16[num] = det->nFrames->data.U16[i];
    418             det->padding->data.S16[num] = det->padding->data.S16[i];
    419             num++;
    420         }
    421     }
    422     det->x->n = num;
    423     det->y->n = num;
    424     det->ra->n = num;
    425     det->dec->n = num;
    426     det->raErr->n = num;
    427     det->decErr->n = num;
    428     det->mag->n = num;
    429     det->magErr->n = num;
    430     det->chi2->n = num;
    431     det->dof->n = num;
    432     det->cr->n = num;
    433     det->extended->n = num;
    434     det->psfMajor->n = num;
    435     det->psfMinor->n = num;
    436     det->psfTheta->n = num;
    437     det->quality->n = num;
    438     det->numPix->n = num;
    439     det->xxMoment->n = num;
    440     det->xyMoment->n = num;
    441     det->yyMoment->n = num;
    442     det->flags->n = num;
    443     det->diffSkyfileId->n = num;
    444     det->naxis1->n = num;
    445     det->naxis2->n = num;
    446     det->mask->n = num;
    447     det->num = num;
    448     det->nPos->n = num;
    449     det->fPos->n = num;
    450     det->nRatioBad->n = num;
    451     det->nRatioMask->n = num;
    452     det->nRatioAll->n = num;
    453     det->psfInstFlux->n = num;
    454     det->psfInstFluxSig->n = num;
    455     det->apMag->n = num;
    456     det->apMagRadius->n = num;
    457     det->apMagRaw->n = num;
    458     det->apFlux->n = num;
    459     det->apFluxSig->n = num;
    460     det->peakFluxAsMag->n = num;
    461     det->calPsfMag->n = num;
    462     det->calPsfMagSig->n = num;
    463     det->sky->n = num;
    464     det->skySig->n = num;
    465     det->qualityPerfect->n = num;
    466     det->momentsR1->n = num;
    467     det->momentsRH->n = num;
    468     det->kronFlux->n = num;
    469     det->kronFluxErr->n = num;
    470     det->kronFluxInner->n = num;
    471     det->kronFluxOuter->n = num;
    472     det->diffRP->n = num;
    473     det->diffSnP->n = num;
    474     det->diffRM->n = num;
    475     det->diffSnM->n = num;
    476     det->flags2->n = num;
    477     det->ippIdet->n = num;
    478     det->nFrames->n = num;
    479     det->padding->n = num;
    480     return true;
    481 }

  • tags/ipp-20110622/ppTranslate/src/ppMopsMerge.c

    r32189 r32407  
    1717#define AIRMASS_TOL 1.0e-3              // Tolerance for airmass matching
    1818
     19#if 0
     20#undef psTrace
     21#define psTrace(facil, level, ...)              \
     22    if (level <= 5) {                           \
     23        fprintf(stderr, __VA_ARGS__);           \
     24    }
     25#endif
     26
    1927// Get distance from detection to centre of image
    2028static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
    … …  
    2230    )
    2331{
    24   float dx = (float) (detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0);
    25   float dy = (float) (detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0);
    26   return PS_SQR(dx) + PS_SQR(dy);
     32    float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
     33    float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
     34    return PS_SQR(dx) + PS_SQR(dy);
    2735}
    2836
    2937
    30 ppMopsDetections *ppMopsMerge(const psArray *detections)
     38bool ppMopsPurgeDuplicates(const psArray *detections)
    3139{
    3240    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    3341
    34     psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
    35 
    36     ppMopsDetections *merged = NULL;    // Merged list
    37     int num = 0;                                                         // Number of merged files
    38     psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
     42    int numInputs = detections->n;                // Number of inputs
     43    psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
     44
     45    long total = 0;                                // Total number of sources
     46    psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
     47    psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input
     48    psVectorInit(dupNum, 0);
     49    for (int i = 0; i < numInputs; i++) {
     50        ppMopsDetections *det = detections->data[i]; // Detections from
     51        if (!det || det->num == 0) {
     52            continue;
     53        }
     54        psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8);
     55        psVectorInit(dupes, 0);
     56        total += det->num;
     57    }
     58    psTrace("ppMops.merge", 2, "Total detections: %ld\n", total);
     59
     60    psVector *raMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged RAs
     61    psVector *decMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged Decs
     62    psVector *sourceMerged = psVectorAlloc(total, PS_TYPE_U16); // Source image of merged sources
     63    psVector *indexMerged = psVectorAlloc(total, PS_TYPE_U32); // Index of merged sources
     64    long num = 0;                                                   // Number of merged sources
     65
     66    const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates
     67    const char *filter = NULL;                    // Filter name
     68    float airmass = NAN, exptime = NAN;           // Exposure details
     69    double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing
     70    double mjd = NAN;                           // Time of exposure
     71
     72    for (int i = 0; i < numInputs; i++) {
     73        ppMopsDetections *det = detections->data[i]; // Detections of interest
     74        if (!det) {
     75            psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i);
     76            continue;
     77        } else if (det->num == 0) {
     78            psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i);
     79            continue;
     80        }
     81
     82        // Check exposure characteristics
     83        if (num == 0) {
     84            raBoresight = det->raBoresight;
     85            decBoresight = det->decBoresight;
     86            filter = det->filter;
     87            airmass = det->airmass;
     88            exptime = det->exptime;
     89            posangle = det->posangle;
     90            alt = det->alt;
     91            az = det->az;
     92        } else {
     93            if (strcmp(raBoresight, det->raBoresight) != 0) {
     94                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     95                        raBoresight, det->raBoresight);
     96                return false;
     97            }
     98            if (strcmp(decBoresight, det->decBoresight) != 0) {
     99                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     100                        decBoresight, det->decBoresight);
     101                return false;
     102            }
     103            if (strcmp(filter, det->filter) != 0) {
     104                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     105                        filter, det->filter);
     106                return false;
     107            }
     108            if (fabsf(airmass - det->airmass) > AIRMASS_TOL) {
     109                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     110                        airmass, det->airmass);
     111                return false;
     112            }
     113            if (fabsf(exptime - det->exptime) > EXPTIME_TOL) {
     114                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     115                        exptime, det->exptime);
     116                return false;
     117            }
     118            if (fabs(posangle - det->posangle) > POSANGLE_TOL) {
     119                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     120                        posangle, det->posangle);
     121                return false;
     122            }
     123            if (fabs(alt - det->alt) > BORESIGHT_TOL) {
     124                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     125                        alt, det->alt);
     126                return false;
     127            }
     128            if (fabs(az - det->az) > BORESIGHT_TOL) {
     129                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     130                        az, det->az);
     131                return false;
     132            }
     133            if (fabs(mjd - det->mjd) > MJD_TOL) {
     134                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     135                        mjd, det->mjd);
     136                return false;
     137            }
     138        }
     139
     140        psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
     141        memcpy(&raMerged->data.F64[num], det->ra->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     142        memcpy(&decMerged->data.F64[num], det->dec->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     143        for (long j = 0, k = num; j < det->num; j++, k++) {
     144            sourceMerged->data.U16[k] = i;
     145            indexMerged->data.U32[k] = j;
     146        }
     147        num += det->num;
     148    }
     149
     150    psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num);
     151    psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree
     152    if (!tree) {
     153        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     154        return false;
     155    }
     156
    39157    for (int i = 0; i < detections->n; i++) {
    40158        ppMopsDetections *det = detections->data[i]; // Detections of interest
    41159        if (!det) {
    42             psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i);
    43160            continue;
    44161        } else if (det->num == 0) {
    45             psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i);
    46             continue;
    47         }
    48         num++;
    49         if (!merged) {
    50             psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
    51             merged = psMemIncrRefCounter(det);
    52             continue;
    53         }
    54         psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
    55 
    56         // XXX compare exposure properties
    57         if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
    58             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
    59                     merged->raBoresight, det->raBoresight);
    60             return NULL;
    61         }
    62         if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
    63             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
    64                     merged->decBoresight, det->decBoresight);
    65             return NULL;
    66         }
    67         if (strcmp(merged->filter, det->filter) != 0) {
    68             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
    69                     merged->filter, det->filter);
    70             return NULL;
    71         }
    72 
    73         if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
    74             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
    75                     merged->airmass, det->airmass);
    76             return NULL;
    77         }
    78         if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
    79             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
    80                     merged->exptime, det->exptime);
    81             return NULL;
    82         }
    83         if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
    84             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
    85                     merged->posangle, det->posangle);
    86             return NULL;
    87         }
    88         if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
    89             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
    90                     merged->alt, det->alt);
    91             return NULL;
    92         }
    93         if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
    94             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
    95                     merged->az, det->az);
    96             return NULL;
    97         }
    98         if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
    99             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
    100                     merged->mjd, det->mjd);
    101             return NULL;
    102         }
    103 
    104         merged->seeing += det->seeing;  // Taking average
    105 
    106         psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
    107         if (!tree) {
    108             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
    109             psFree(merged);
    110             return NULL;
    111         }
    112 
     162            continue;
     163        }
     164        psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
     165
     166        psVector *dupes = duplicates->data[i];            // Duplicates list
     167        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
    113168        for (int j = 0; j < det->num; j++) {
     169            // XXX: added by Bill
     170            if (det->mask->data.U8[j]) {
     171                // we marked this source marked as bad when we read it, mark it as a duplicate
     172                dupes->data.U8[j] = 0xFF;
     173                continue;
     174            }
     175
    114176            coords->data.F64[0] = det->ra->data.F64[j];
    115177            coords->data.F64[1] = det->dec->data.F64[j];
    … …  
    119181                psFree(coords);
    120182                psFree(tree);
    121                 psFree(merged);
    122                 return NULL;
    123             }
    124             if (indices->n == 0) {
    125                 psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i);
     183                return false;
     184            }
     185            psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
     186            psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
     187
     188            if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) {
     189                // It's myself
    126190                psFree(indices);
    127                 ppMopsDetectionsCopySingle(merged, det, j);
    128191                continue;
    129192            }
    130             psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
    131193
    132194            // Which one do we keep?
    … …  
    134196            long bestIndex = -1;           // Index with best distance
    135197            for (int k = 0; k < indices->n; k++) {
    136                 long index = indices->data.S64[k]; // Index of point
    137                 float distance = mergeDistance(merged, index); // Distance to centre of image
     198                long mergeIndex = indices->data.S64[k]; // Index of point in merged list
     199                int source = sourceMerged->data.U16[mergeIndex]; // Source image
     200                if (source == i) {
     201                    continue;
     202                }
     203                long index = indexMerged->data.U32[mergeIndex];  // Index in source
     204                psVector *dupes = duplicates->data[source];     // Duplicates list
     205                if (dupes->data.U8[index]) {
     206                    continue;
     207                }
     208
     209                float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
    138210                if (distance < bestDistance) {
    139211                    bestDistance = distance;
    … …  
    142214            }
    143215
    144             float distance = mergeDistance(det, j); // Distance to centre of image
    145             if (distance < bestDistance) {
    146                 psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
    147                 // Blow away existing sources
    148                 for (int k = 0; k < indices->n; k++) {
    149                     long index = indices->data.S64[k]; // Index of point
    150                     merged->mask->data.U8[index] = 0xFF;
     216            if (bestIndex >= 0) {
     217                float distance = mergeDistance(det, j); // Distance to centre of image
     218                if (bestIndex >= 0 && distance < bestDistance) {
     219                    psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
     220                    for (int k = 0; k < indices->n; k++) {
     221                        long mergeIndex = indices->data.S64[k]; // Index of point
     222                        int source = sourceMerged->data.U16[mergeIndex]; // Source image
     223                        if (source == i) {
     224                            continue;
     225                        }
     226                        long index = indexMerged->data.U32[mergeIndex];  // Index in source
     227                        psVector *dupes = duplicates->data[source];      // Duplicates list
     228                        if (!dupes->data.U8[index]) {
     229                            dupes->data.U8[index] = 0xFF;
     230                            dupNum->data.U32[source]++;
     231                        }
     232                    }
     233                } else {
     234                    psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
     235                    dupes->data.U8[j] = 0xFF;
     236                    dupNum->data.U32[i]++;
    151237                }
    152                 ppMopsDetectionsCopySingle(merged, det, j);
    153             } else {
    154                 psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    155238            }
    156239            psFree(indices);
    157240        }
    158 
    159         psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
    160 
    161         psFree(tree);
    162         ppMopsDetectionsPurge(merged);
    163     }
    164     psFree(coords);
    165 
    166     if (num == 0) {
    167         //All detections were NULL?!
    168         psTrace("ppMops.merge", 3, "All %ld detections were NULL\n", detections->n);
    169         return NULL;
    170     }
    171     psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
    172 
    173     merged->seeing /= (float) num;
    174 
    175     return merged;
     241        psFree(coords);
     242    }
     243    psFree(tree);
     244    psFree(raMerged);
     245    psFree(decMerged);
     246    psFree(sourceMerged);
     247    psFree(indexMerged);
     248
     249    // Remove duplicates
     250    for (int i = 0; i < detections->n; i++) {
     251        ppMopsDetections *det = detections->data[i]; // Detections of interest
     252        if (!det) {
     253            continue;
     254        } else if (det->num == 0) {
     255            continue;
     256        }
     257        psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
     258
     259#define VECTOR_PURGE_CASE(TYPE)                                     \
     260        case PS_TYPE_##TYPE:    {                                   \
     261          long j = 0;                                               \
     262          for (long i = 0; i < vector->n; i++) {                    \
     263              if (!dupes->data.U8[i]) {                             \
     264                  if (i == j) {                                     \
     265                      j++;                                          \
     266                      continue;                                     \
     267                  }                                                 \
     268                  vector->data.TYPE[j++] = vector->data.TYPE[i];    \
     269              }                                                     \
     270          }                                                         \
     271          vector->n = j;                                            \
     272        }                                                           \
     273        break
     274
     275        psVector *dupes = duplicates->data[i]; // Duplicates
     276        psArray *table = psListToArray(det->table->list); // Table of data
     277        long newLength = -1;
     278        for (int t = 0; t < table->n; t++) {
     279            psMetadataItem *item = table->data[t]; // Table item
     280            psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type);
     281            psVector *vector = item->data.V; // Vector to purge
     282            switch (vector->type.type) {
     283                VECTOR_PURGE_CASE(U8);
     284                VECTOR_PURGE_CASE(U16);
     285                VECTOR_PURGE_CASE(U32);
     286                VECTOR_PURGE_CASE(U64);
     287                VECTOR_PURGE_CASE(S8);
     288                VECTOR_PURGE_CASE(S16);
     289                VECTOR_PURGE_CASE(S32);
     290                VECTOR_PURGE_CASE(S64);
     291                VECTOR_PURGE_CASE(F32);
     292                VECTOR_PURGE_CASE(F64);
     293              default:
     294                psAbort("Unrecognised vector type: %x", vector->type.type);
     295            }
     296            if (newLength == -1) {
     297                newLength = vector->n;
     298            } else if (vector->n != newLength) {
     299                psAbort("Unexpected new length found : %ld expected: %ld",
     300                                                vector->n, newLength);
     301            }
     302               
     303        }
     304        // XXX IS this safe? Perhaps save in numGood?
     305        det->num = newLength;
     306        psFree(table);
     307    }
     308    psFree(dupNum);
     309    psFree(duplicates);
     310
     311    return true;
    176312}
    177313

  • tags/ipp-20110622/ppTranslate/src/ppMopsRead.c

    r32397 r32407  
    2020    psArray *detections = psArrayAlloc(num); // Array of detections, to return
    2121    for (int i = 0; i < num; i++) {
    22         psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file
    23 
     22        const char *name = inNames->data[i];
     23
     24        psFits *fits = psFitsOpen(name,  "r"); // FITS file
    2425        if (!fits) {
    2526            psError(PS_ERR_IO, false, "Unable to open input %d", i);
    2627            return false;
    2728        }
     29
    2830        psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
    2931        if (!header) {
    … …  
    7476        }
    7577        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     78        detections->data[i] = det;
     79        det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
     80        det->num = size;
     81        det->diffSkyfileId = diffSkyfileId;
    7682
    7783        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    … …  
    9096                                     psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    9197
    92         int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    93         int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     98        det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     99        det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    94100
    95101        psFree(header);
    96102
    97 #ifndef USE_OLD_READ_TABLE
    98         // use new fits table functions
    99 
    100         psFitsTable *table = psFitsReadTableNew(fits); // Table of interest
     103        psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
    101104        if (!table) {
    102105            psError(PS_ERR_IO, false, "Unable to read table %d", i);
    103             return false;
    104         }
     106            return NULL;
     107        }
     108        det->table = table;
    105109        psFitsClose(fits);
     110        if (args->version == 0) {
     111          if (skyChipPsfVersion < 2) {
     112            // XXX: TODO: Do we need to add dummy vectors for the missing columns?
     113           }
     114        }
     115
     116        psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
     117        psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     118
     119        det->raErr = psVectorAlloc(size, PS_TYPE_F64);
     120        det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     121        det->mask = psVectorAlloc(size, PS_TYPE_U8);
     122
     123        // convert ra and dec to radians for use in the purge duplicates function
     124        det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     125        det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     126        det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
     127        det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
     128        if (!det->ra || !det->dec || !det->x || !det->y) {
     129            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
     130            return NULL;
     131        }
     132
     133        // Add our new vectors to the table so that duplicates and masked items may be purged
     134        psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
     135        psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     136
     137        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
     138
     139        psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
     140        psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
     141        psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
     142        psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
     143        psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
     144        psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
     145        psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
    106146
    107147        double plateScale = 0.0;        // Plate scale
    … …  
    109149        for (long row = 0; row < size; row++) {
    110150
    111             psU32 flags = psFitsTableGetU32(NULL, table, row, "FLAGS");
     151            psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    112152            if (flags & SOURCE_MASK) {
    113153                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     154                det->mask->data.U8[row] = 0xFF;
    114155                continue;
    115156            }
    116157
    117             det->x->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "X_PSF");
    118             det->y->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "Y_PSF");
    119             det->ra->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "RA_PSF"));
    120             det->dec->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "DEC_PSF"));
    121             det->mag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG");
    122             det->magErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG_SIG");
    123             det->chi2->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_CHISQ");
    124             det->dof->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NDOF");
    125             det->cr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CR_NSIGMA");
    126             det->extended->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "EXT_NSIGMA");
    127             det->psfMajor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MAJOR");
    128             det->psfMinor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MINOR");
    129             det->psfTheta->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_THETA");
    130             det->quality->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF");
    131             det->numPix->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NPIX");
    132             det->xxMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XX");
    133             det->xyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XY");
    134             det->yyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_YY");
    135             det->flags->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "FLAGS");
    136             det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
    137             det->naxis1->data.S32[numGood] = naxis1;
    138             det->naxis2->data.S32[numGood] = naxis2;
    139             det->nPos->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "DIFF_NPOS");
    140             det->fPos->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_FRATIO");
    141             det->nRatioBad->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_BAD");
    142             det->nRatioMask->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_MASK");
    143             det->nRatioAll->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_ALL");
    144 
    145             //Additions of 2010-10-25
    146             if (args->version == 2) {
    147               //Values are set only if the version is 2
    148               if (skyChipPsfVersion == 2) {
    149                 det->psfInstFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX");
    150                 det->psfInstFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX_SIG");
    151                 det->apMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG");
    152                 det->apMagRaw->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RAW");
    153                 det->apMagRadius->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RADIUS");
    154                 det->apFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX");
    155                 det->apFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX_SIG");
    156                 det->peakFluxAsMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PEAK_FLUX_AS_MAG");
    157                 det->calPsfMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG");
    158                 det->calPsfMagSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG_SIG");
    159                 det->sky->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY");
    160                 det->skySig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY_SIGMA");
    161                 det->qualityPerfect->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF_PERFECT");
    162                 det->momentsR1->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_R1");
    163                 det->momentsRH->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_RH");
    164                 det->kronFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX");
    165                 det->kronFluxErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_ERR");
    166                 det->kronFluxInner->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_INNER");
    167                 det->kronFluxOuter->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_OUTER");
    168                 det->diffRP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_P");
    169                 det->diffSnP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_P");
    170                 det->diffRM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_M");
    171                 det->diffSnM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_M");
    172                 det->flags2->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "FLAGS2");
    173                 det->ippIdet->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "IPP_IDET");
    174                 det->nFrames->data.U16[numGood] = psFitsTableGetU16(NULL, table, row, "N_FRAMES");
    175                 det->padding->data.S16[numGood] = psFitsTableGetS16(NULL, table, row, "PADDING");
    176               } else {
    177                 det->psfInstFlux->data.F32[numGood] = NAN;
    178                 det->psfInstFluxSig->data.F32[numGood] = NAN;
    179                 det->apMag->data.F32[numGood] = NAN;
    180                 det->apMagRaw->data.F32[numGood] = NAN;
    181                 det->apMagRadius->data.F32[numGood] = NAN;
    182                 det->apFlux->data.F32[numGood] = NAN;
    183                 det->apFluxSig->data.F32[numGood] = NAN;
    184                 det->peakFluxAsMag->data.F32[numGood] = NAN;
    185                 det->calPsfMag->data.F32[numGood] = NAN;
    186                 det->calPsfMagSig->data.F32[numGood] = NAN;
    187                 det->sky->data.F32[numGood] = NAN;
    188                 det->skySig->data.F32[numGood] = NAN;
    189                 det->qualityPerfect->data.F32[numGood] = NAN;
    190                 det->momentsR1->data.F32[numGood] = NAN;
    191                 det->momentsRH->data.F32[numGood] = NAN;
    192                 det->kronFlux->data.F32[numGood] = NAN;
    193                 det->kronFluxErr->data.F32[numGood] = NAN;
    194                 det->kronFluxInner->data.F32[numGood] = NAN;
    195                 det->kronFluxOuter->data.F32[numGood] = NAN;
    196                 det->diffRP->data.F32[numGood] = NAN;
    197                 det->diffSnP->data.F32[numGood] = NAN;
    198                 det->diffRM->data.F32[numGood] = NAN;
    199                 det->diffSnM->data.F32[numGood] = NAN;
    200                 det->flags2->data.U32[numGood] = 0;
    201                 det->ippIdet->data.U32[numGood] = 0;
    202                 det->nFrames->data.U16[numGood] = 0;
    203                 det->padding->data.S16[numGood] = 0;
    204               }
    205             }
    206 
    207158            // Calculate error in RA, Dec
    208             double xErr = psFitsTableGetF64(NULL, table, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32
    209             double yErr = psFitsTableGetF64(NULL, table, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32
    210             double scale = psFitsTableGetF64(NULL, table, row, "PLTSCALE"); //SC: Warning! Promotion of F32
    211             double angle = psFitsTableGetF64(NULL, table, row, "POSANGLE"); //SC: Warning! Promotion of F32
    212 
    213             if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
    214                 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
    215                 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
     159
     160            double xErr = xErrV->data.F32[row];
     161            double yErr = yErrV->data.F32[row];
     162            double scale = scaleV->data.F32[row];
     163            double angle = angleV->data.F32[row];
     164
     165            if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
     166                !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
     167                !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
    216168                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    217169                psTrace("ppMops.read", 10,
    … …  
    219171                        "%f %f %lf %lf %f %f %f %f %f %f",
    220172                        row, i,
    221                         det->x->data.F32[numGood], det->y->data.F32[numGood],
    222                         det->ra->data.F64[numGood], det->dec->data.F64[numGood],
    223                         det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
     173                        det->x->data.F32[row], det->y->data.F32[row],
     174                        det->ra->data.F64[row], det->dec->data.F64[row],
     175                        mag->data.F32[row], magErr->data.F32[row],
    224176                        xErr, yErr, scale, angle);
     177                det->mask->data.U8[row] = 0xFF;
    225178                continue;
    226179            }
    … …  
    231184            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    232185            double errScale = scale / 3600.0;
    233             det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    234             det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    235 
    236             det->mask->data.U8[numGood] = 0;
     186            det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     187            det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     188
     189            det->mask->data.U8[row] = 0;
    237190            plateScale += scale;
    238191            numGood++;
    239192        }
    240 #else
    241     // OLD way of reading fits tables read it into an array of metadata objects, one for each row. Uses lots and lots
    242     // and lots of memory.
    243         psArray *table = psFitsReadTable(fits); // Table of interest
    244         if (!table) {
    245             psError(PS_ERR_IO, false, "Unable to read table %d", i);
    246             return false;
    247         }
    248         psFitsClose(fits);
    249 
    250         double plateScale = 0.0;        // Plate scale
    251         long numGood = 0;               // Number of good rows
    252         for (long j = 0; j < size; j++) {
    253             psMetadata *row = table->data[j]; // Row of interest
    254 
    255             psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
    256             if (flags & SOURCE_MASK) {
    257                 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
    258                 continue;
    259             }
    260 
    261             det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
    262             det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
    263             det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
    264             det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
    265             det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
    266             det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
    267             det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");
    268             det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");
    269             det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");
    270             det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
    271             det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");
    272             det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");
    273             det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");
    274             det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");
    275             det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");
    276             det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");
    277             det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");
    278             det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");
    279             det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
    280             det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
    281             det->naxis1->data.S32[numGood] = naxis1;
    282             det->naxis2->data.S32[numGood] = naxis2;
    283             det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");
    284             det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");
    285             det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");
    286             det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");
    287             det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");
    288 
    289             //Additions of 2010-10-25
    290             if (args->version == 2) {
    291               //Values are set only if the version is 2
    292               if (skyChipPsfVersion == 2) {
    293                 det->psfInstFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX");
    294                 det->psfInstFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX_SIG");
    295                 det->apMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG");
    296                 det->apMagRaw->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RAW");
    297                 det->apMagRadius->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RADIUS");
    298                 det->apFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX");
    299                 det->apFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX_SIG");
    300                 det->peakFluxAsMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PEAK_FLUX_AS_MAG");
    301                 det->calPsfMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG");
    302                 det->calPsfMagSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG_SIG");
    303                 det->sky->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY");
    304                 det->skySig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY_SIGMA");
    305                 det->qualityPerfect->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF_PERFECT");
    306                 det->momentsR1->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_R1");
    307                 det->momentsRH->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_RH");
    308                 det->kronFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX");
    309                 det->kronFluxErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_ERR");
    310                 det->kronFluxInner->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_INNER");
    311                 det->kronFluxOuter->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_OUTER");
    312                 det->diffRP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_P");
    313                 det->diffSnP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_P");
    314                 det->diffRM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_M");
    315                 det->diffSnM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_M");
    316                 det->flags2->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS2");
    317                 det->ippIdet->data.U32[numGood] = psMetadataLookupU32(NULL, row, "IPP_IDET");
    318                 det->nFrames->data.U16[numGood] = psMetadataLookupU16(NULL, row, "N_FRAMES");
    319                 det->padding->data.S16[numGood] = psMetadataLookupS16(NULL, row, "PADDING");
    320               } else {
    321                 det->psfInstFlux->data.F32[numGood] = NAN;
    322                 det->psfInstFluxSig->data.F32[numGood] = NAN;
    323                 det->apMag->data.F32[numGood] = NAN;
    324                 det->apMagRaw->data.F32[numGood] = NAN;
    325                 det->apMagRadius->data.F32[numGood] = NAN;
    326                 det->apFlux->data.F32[numGood] = NAN;
    327                 det->apFluxSig->data.F32[numGood] = NAN;
    328                 det->peakFluxAsMag->data.F32[numGood] = NAN;
    329                 det->calPsfMag->data.F32[numGood] = NAN;
    330                 det->calPsfMagSig->data.F32[numGood] = NAN;
    331                 det->sky->data.F32[numGood] = NAN;
    332                 det->skySig->data.F32[numGood] = NAN;
    333                 det->qualityPerfect->data.F32[numGood] = NAN;
    334                 det->momentsR1->data.F32[numGood] = NAN;
    335                 det->momentsRH->data.F32[numGood] = NAN;
    336                 det->kronFlux->data.F32[numGood] = NAN;
    337                 det->kronFluxErr->data.F32[numGood] = NAN;
    338                 det->kronFluxInner->data.F32[numGood] = NAN;
    339                 det->kronFluxOuter->data.F32[numGood] = NAN;
    340                 det->diffRP->data.F32[numGood] = NAN;
    341                 det->diffSnP->data.F32[numGood] = NAN;
    342                 det->diffRM->data.F32[numGood] = NAN;
    343                 det->diffSnM->data.F32[numGood] = NAN;
    344                 det->flags2->data.U32[numGood] = 0;
    345                 det->ippIdet->data.U32[numGood] = 0;
    346                 det->nFrames->data.U16[numGood] = 0;
    347                 det->padding->data.S16[numGood] = 0;
    348               }
    349             }
    350 
    351             // Calculate error in RA, Dec
    352             double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32
    353             double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32
    354             double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); //SC: Warning! Promotion of F32
    355             double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); //SC: Warning! Promotion of F32
    356 
    357             if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
    358                 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
    359                 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
    360                 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    361                 psTrace("ppMops.read", 10,
    362                         "Discarding row %ld from input %d because of non-finite values: "
    363                         "%f %f %lf %lf %f %f %f %f %f %f",
    364                         j, i,
    365                         det->x->data.F32[numGood], det->y->data.F32[numGood],
    366                         det->ra->data.F64[numGood], det->dec->data.F64[numGood],
    367                         det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
    368                         xErr, yErr, scale, angle);
    369                 continue;
    370             }
    371 
    372             // XXX Not at all sure I've got the angles around the right way here...
    373             double cosAngle = cos(angle), sinAngle = sin(angle);
    374             double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
    375             double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    376             double errScale = scale / 3600.0;
    377             det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    378             det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    379 
    380             det->mask->data.U8[numGood] = 0;
    381             plateScale += scale;
    382             numGood++;
    383         }
    384 #endif
    385193        det->seeing *= ((float) plateScale) / ((float) numGood);
    386194
    387         det->x->n = numGood;
    388         det->y->n = numGood;
    389         det->ra->n = numGood;
    390         det->dec->n = numGood;
    391         det->raErr->n = numGood;
    392         det->decErr->n = numGood;
    393         det->mag->n = numGood;
    394         det->magErr->n = numGood;
    395         det->chi2->n = numGood;
    396         det->dof->n = numGood;
    397         det->cr->n = numGood;
    398         det->extended->n = numGood;
    399         det->psfMajor->n = numGood;
    400         det->psfMinor->n = numGood;
    401         det->psfTheta->n = numGood;
    402         det->quality->n = numGood;
    403         det->numPix->n = numGood;
    404         det->xxMoment->n = numGood;
    405         det->xyMoment->n = numGood;
    406         det->yyMoment->n = numGood;
    407         det->flags->n = numGood;
    408         det->diffSkyfileId->n = numGood;
    409         det->naxis1->n = numGood;
    410         det->naxis2->n = numGood;
    411         det->mask->n = numGood;
    412         det->nPos->n = numGood;
    413         det->fPos->n = numGood;
    414         det->nRatioBad->n = numGood;
    415         det->nRatioMask->n = numGood;
    416         det->nRatioAll->n = numGood;
    417         det->psfInstFlux->n = numGood;
    418         det->psfInstFluxSig->n = numGood;
    419         det->apMag->n = numGood;
    420         det->apMagRaw->n = numGood;
    421         det->apMagRadius->n = numGood;
    422         det->apFlux->n = numGood;
    423         det->apFluxSig->n = numGood;
    424         det->peakFluxAsMag->n = numGood;
    425         det->calPsfMag->n = numGood;
    426         det->calPsfMagSig->n = numGood;
    427         det->sky->n = numGood;
    428         det->skySig->n = numGood;
    429         det->qualityPerfect->n = numGood;
    430         det->momentsR1->n = numGood;
    431         det->momentsRH->n = numGood;
    432         det->kronFlux->n = numGood;
    433         det->kronFluxErr->n = numGood;
    434         det->kronFluxInner->n = numGood;
    435         det->kronFluxOuter->n = numGood;
    436         det->diffRP->n = numGood;
    437         det->diffSnP->n = numGood;
    438         det->diffRM->n = numGood;
    439         det->diffSnM->n = numGood;
    440         det->flags2->n = numGood;
    441         det->ippIdet->n = numGood;
    442         det->nFrames->n = numGood;
    443         det->padding->n = numGood;
    444 
    445         det->num = numGood;
     195        // Are we using numGood for anything outside of this function?
     196        det->numGood = numGood;
    446197
    447198        if (isfinite(args->zp) && numGood > 0) {
    448             psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    449         }
    450 
    451         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    452 
    453         psFree(table);
    454         detections->data[i] = det;
     199            psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
     200        }
     201
     202        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
    455203    }
    456204

  • tags/ipp-20110622/ppTranslate/src/ppMopsWrite.c

    r29567 r32407  
    99#include "ppTranslateVersion.h"
    1010
    11 bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args)
     11static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32);
     12static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName);
     13
     14bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
    1215{
    13     psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output);
    1416
    1517    psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    … …  
    3638    psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
    3739
     40    // Get these header words from the first input
     41    ppMopsDetections *det = detections->data[0];
    3842    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
    3943    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
    … …  
    5559    psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
    5660
    57     if (det->num == 0) {
     61    // Find the total number of detections
     62
     63    long total = 0;
     64    for (long i=0; i<detections->n; i++) {
     65        ppMopsDetections *det = detections->data[i];
     66        total += det->num;
     67    }
     68
     69    psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
     70
     71    if (total == 0) {
    5872        // Write dummy table
    5973        psMetadata *row = psMetadataAlloc(); // Output row
    … …  
    151165        psFree(row);
    152166    } else {
    153         psArray *table = psArrayAlloc(det->num); // Table to write
    154         for (long i = 0; i < det->num; i++) {
    155             psMetadata *row = psMetadataAlloc(); // Output row
    156             psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)",
    157                              RAD_TO_DEG(det->ra->data.F64[i]));
    158             psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)",
    159                              det->raErr->data.F64[i]);
    160             psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)",
    161                              RAD_TO_DEG(det->dec->data.F64[i]));
    162             psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)",
    163                              det->decErr->data.F64[i]);
    164             psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]);
    165             psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]);
    166             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", det->chi2->data.F32[i]);
    167             psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit",
    168                              det->dof->data.S32[i]);
    169             psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR",
    170                              det->cr->data.F32[i]);
    171             psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness",
    172                              det->extended->data.F32[i]);
    173             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", det->psfMajor->data.F32[i]);
    174             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", det->psfMinor->data.F32[i]);
    175             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)",
    176                              det->psfTheta->data.F32[i]);
    177             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor",
    178                              det->quality->data.F32[i]);
    179             psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF",
    180                              det->numPix->data.S32[i]);
    181             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", det->xxMoment->data.F32[i]);
    182             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", det->xyMoment->data.F32[i]);
    183             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", det->yyMoment->data.F32[i]);
    184             psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels",
    185                              det->nPos->data.S32[i]);
    186             psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels",
    187                              det->fPos->data.F32[i]);
    188             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative",
    189                              det->nRatioBad->data.F32[i]);
    190             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked",
    191                              det->nRatioMask->data.F32[i]);
    192             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all",
    193                              det->nRatioAll->data.F32[i]);
    194             psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]);
    195             psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile",
    196                              det->diffSkyfileId->data.S64[i]);
    197 
    198             if (args->version == 2) {
    199               // Write data of version 2 (see ICD)
    200               psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
    201                              det->ippIdet->data.U32[i]);
    202               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
    203                              det->psfInstFlux->data.F32[i]);
    204               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
    205                              det->psfInstFluxSig->data.F32[i]);
    206               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
    207                              det->apMag->data.F32[i]);
    208               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
    209                              det->apMagRaw->data.F32[i]);
    210               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
    211                              det->apMagRadius->data.F32[i]);
    212               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
    213                              det->apFlux->data.F32[i]);
    214               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
    215                              det->apFluxSig->data.F32[i]);
    216               psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
    217                              det->peakFluxAsMag->data.F32[i]);
    218               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
    219                              det->calPsfMag->data.F32[i]);
    220               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
    221                              det->calPsfMagSig->data.F32[i]);
    222               psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
    223                              det->sky->data.F32[i]);
    224               psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
    225                              det->skySig->data.F32[i]);
    226               psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
    227                              det->qualityPerfect->data.F32[i]);
    228               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
    229                              det->momentsR1->data.F32[i]);
    230               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
    231                              det->momentsRH->data.F32[i]);
    232               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
    233                              det->kronFlux->data.F32[i]);
    234               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
    235                              det->kronFluxErr->data.F32[i]);
    236               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
    237                              det->kronFluxInner->data.F32[i]);
    238               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
    239                              det->kronFluxOuter->data.F32[i]);
    240               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
    241                              det->diffRP->data.F32[i]);
    242               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
    243                              det->diffSnP->data.F32[i]);
    244               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
    245                              det->diffRM->data.F32[i]);
    246               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
    247                              det->diffSnM->data.F32[i]);
    248               psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
    249                              det->flags2->data.U32[i]);
    250               psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
    251                              det->nFrames->data.U16[i]);
    252               psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
    253                              det->padding->data.S16[i]);
    254             }
    255 
    256             //Update with the table with the current row
    257             table->data[i] = row;
    258         }
    259         if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) {
    260             psErrorStackPrint(stderr, "Unable to write table.");
    261             psFree(header);
    262             psFree(table);
     167
     168#define addColumn(_outName, _inName, _convertTo32) \
     169        if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
     170            psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
     171            return false; \
     172        }
     173
     174        // Allocate the output table
     175        psMetadata *table = psMetadataAlloc();
     176        addColumn("RA", "RA_PSF", 0);
     177        addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     178        addColumn("DEC", "DEC_PSF", 0);
     179        addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     180        addColumn("MAG", "PSF_INST_MAG", 0);
     181        addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
     182        addColumn("PSF_CHI2", "PSF_CHISQ", 0);
     183        addColumn("PSF_DOF", "PSF_NDOF", 1);
     184        addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
     185        addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
     186        addColumn("PSF_MAJOR", NULL, 0);
     187        addColumn("PSF_MINOR", NULL, 0);     
     188        addColumn("PSF_THETA", NULL, 0);   
     189        addColumn("PSF_QUALITY", "PSF_QF", 0);
     190        addColumn("PSF_NPIX", NULL, 1);       
     191        addColumn("MOMENTS_XX", NULL, 0);
     192        addColumn("MOMENTS_XY", NULL, 0);
     193        addColumn("MOMENTS_YY", NULL, 0);
     194        addColumn("N_POS", "DIFF_NPOS", 1);
     195        addColumn("F_POS", "DIFF_FRATIO", 0);
     196        addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
     197        addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
     198        addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
     199        addColumn("FLAGS", "FLAGS", 1);
     200        addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
     201        if (args->version == 2) {
     202            addColumn("IPP_IDET", NULL, 1);
     203            addColumn("PSF_INST_FLUX", NULL, 0);
     204            addColumn("PSF_INST_FLUX_SIG", NULL, 0);
     205            addColumn("AP_MAG", NULL, 0);
     206            addColumn("AP_MAG_RAW", NULL, 0);
     207            addColumn("AP_MAG_RADIUS", NULL, 0);
     208            addColumn("AP_FLUX", NULL, 0);
     209            addColumn("AP_FLUX_SIG", NULL, 0);
     210            addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
     211            addColumn("CAL_PSF_MAG", NULL, 0);
     212            addColumn("CAL_PSF_MAG_SIG", NULL, 0);
     213            addColumn("SKY", NULL, 0);
     214            addColumn("SKY_SIGMA", NULL, 0);
     215            addColumn("PSF_QF_PERFECT", NULL, 0);
     216            addColumn("MOMENTS_R1", NULL, 0);
     217            addColumn("MOMENTS_RH", NULL, 0);
     218            addColumn("KRON_FLUX", NULL, 0);
     219            addColumn("KRON_FLUX_ERR", NULL, 0);
     220            addColumn("KRON_FLUX_INNER", NULL, 0);
     221            addColumn("KRON_FLUX_OUTER", NULL, 0);
     222            addColumn("DIFF_R_P", NULL, 0);
     223            addColumn("DIFF_SN_P", NULL, 0);
     224            addColumn("DIFF_R_M", NULL, 0);
     225            addColumn("DIFF_SN_M", NULL, 0);
     226            addColumn("FLAGS2", NULL, 1);
     227            addColumn("IPP_IDET", NULL, 0);
     228            addColumn("N_FRAMES", NULL, 0);
     229            addColumn("PADDING", NULL, 0);
     230        }
     231        if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
     232            psError(psErrorCodeLast(), false, "Unable to write table");
    263233            return false;
    264234        }
    … …  
    273243    return true;
    274244}
     245
     246static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
     247{
     248    if (inColumnName == NULL) {
     249        inColumnName = outColumnName;
     250    }
     251
     252    psVector *out = NULL;
     253    if (convertTo32) {
     254        // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
     255        // We want to write 32 bits to the output.
     256        int next = 0;
     257        for (long i=0; i<detections->n; i++) {
     258            ppMopsDetections *det = detections->data[i];
     259            if (det->num == 0) {
     260                // no detections survived for this input
     261                continue;
     262            }
     263            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     264            if (!in) {
     265                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     266                return false;
     267            }
     268            if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
     269                psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
     270                    inColumnName, in->type.type);
     271                return false;
     272            }
     273            if (out == NULL) {
     274                // First time through set up the output vector and the copy parameters
     275                if (in->type.type == PS_TYPE_S64) {
     276                    out = psVectorAlloc(outputSize, PS_TYPE_S32);
     277                } else {
     278                    out = psVectorAlloc(outputSize, PS_TYPE_U32);
     279                }
     280            }
     281            for (long d=0; d < det->num; d++) {
     282                if (in->type.type == PS_TYPE_S64) {
     283                    out->data.S32[next++] = in->data.S64[d];
     284                } else {
     285                    out->data.U32[next++] = in->data.U64[d];
     286                }
     287            }
     288        }
     289    } else {
     290        void *next = NULL;
     291        int elementSize = 0;    // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy
     292        for (long i=0; i<detections->n; i++) {
     293            ppMopsDetections *det = detections->data[i];
     294            if (det->num == 0) {
     295                // no detections survived for this input
     296                continue;
     297            }
     298            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     299            if (!in) {
     300                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     301                return false;
     302            }
     303            if (out == NULL) {
     304                // First time through set up the output vector and the copy parameters
     305                out = psVectorAlloc(outputSize, in->type.type);
     306                next = (void *) out->data.U8;
     307                switch (in->type.type) {
     308                    case PS_TYPE_S8:
     309                    elementSize = sizeof(psS8);
     310                    break;
     311                case PS_TYPE_U8:
     312                    elementSize = sizeof(psU8);
     313                    break;
     314                case PS_TYPE_S16:
     315                    elementSize = sizeof(psS16);
     316                    break;
     317                case PS_TYPE_U16:
     318                    elementSize = sizeof(psU16);
     319                    break;
     320                case PS_TYPE_S32:
     321                    elementSize = sizeof(psS32);
     322                    break;
     323                case PS_TYPE_U32:
     324                    elementSize = sizeof(psU32);
     325                    break;
     326                case PS_TYPE_S64:
     327                    elementSize = sizeof(psS64);
     328                    break;
     329                case PS_TYPE_U64:
     330                    elementSize = sizeof(psU64);
     331                    break;
     332                case PS_TYPE_F32:
     333                    elementSize = sizeof(psF32);
     334                    break;
     335                case PS_TYPE_F64:
     336                    elementSize = sizeof(psF64);
     337                    break;
     338                default:
     339                    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
     340                    return false;
     341
     342                }
     343            }
     344            // We are doing nasty things here so we can use memcpy.
     345            // It would be safer to do a proper loop over the elements.
     346            long toCopy = det->num * elementSize;
     347            memcpy(next, in->data.U8, toCopy);
     348            next += toCopy;
     349        }
     350    }
     351
     352    // Finally add the new column to the output table
     353    psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     354    psFree(out);    // drop reference
     355
     356    return true;
     357}
     358static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
     359{
     360    psVector *out = psVectorAlloc(total, PS_TYPE_S64);
     361    long next = 0;
     362    for (long i = 0; i<detections->n; i++) {
     363        ppMopsDetections *det = detections->data[i];
     364        psS64 diffSkyfileId = det->diffSkyfileId;
     365        for (long j = 0; j < det->num; j++) {
     366            out->data.S64[next++] = diffSkyfileId;
     367        }
     368    }
     369    psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
     370    psFree(out);
     371    return true;
     372}
Note: See TracChangeset for help on using the changeset viewer.