IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 8, 2011, 2:44:12 PM (15 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20110906
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110906

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMops.c

    r32176 r32630  
    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();

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMops.h

    r29567 r32630  
    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:

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMopsDetections.c

    r30660 r32630  
    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 }

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMopsMerge.c

    r32175 r32630  
    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

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMopsRead.c

    r30656 r32630  
    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         psArray *table = psFitsReadTable(fits); // Table of interest
     103        psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
    98104        if (!table) {
    99105            psError(PS_ERR_IO, false, "Unable to read table %d", i);
    100             return false;
    101         }
     106            return NULL;
     107        }
     108        det->table = table;
    102109        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");
    103146
    104147        double plateScale = 0.0;        // Plate scale
    105148        long numGood = 0;               // Number of good rows
    106         for (long j = 0; j < size; j++) {
    107             psMetadata *row = table->data[j]; // Row of interest
    108 
    109             psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
     149        for (long row = 0; row < size; row++) {
     150
     151            psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    110152            if (flags & SOURCE_MASK) {
    111                 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
     153                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     154                det->mask->data.U8[row] = 0xFF;
    112155                continue;
    113156            }
    114157
    115             det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
    116             det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
    117             det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
    118             det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
    119             det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
    120             det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
    121             det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");
    122             det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");
    123             det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");
    124             det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
    125             det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");
    126             det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");
    127             det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");
    128             det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");
    129             det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");
    130             det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");
    131             det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");
    132             det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");
    133             det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
    134             det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
    135             det->naxis1->data.S32[numGood] = naxis1;
    136             det->naxis2->data.S32[numGood] = naxis2;
    137             det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");
    138             det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");
    139             det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");
    140             det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");
    141             det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");
    142 
    143             //Additions of 2010-10-25
    144             if (args->version == 2) {
    145               //Values are set only if the version is 2
    146               if (skyChipPsfVersion == 2) {
    147                 det->psfInstFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX");
    148                 det->psfInstFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX_SIG");
    149                 det->apMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG");
    150                 det->apMagRaw->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RAW");
    151                 det->apMagRadius->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RADIUS");
    152                 det->apFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX");
    153                 det->apFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX_SIG");
    154                 det->peakFluxAsMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PEAK_FLUX_AS_MAG");
    155                 det->calPsfMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG");
    156                 det->calPsfMagSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG_SIG");
    157                 det->sky->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY");
    158                 det->skySig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY_SIGMA");
    159                 det->qualityPerfect->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF_PERFECT");
    160                 det->momentsR1->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_R1");
    161                 det->momentsRH->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_RH");
    162                 det->kronFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX");
    163                 det->kronFluxErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_ERR");
    164                 det->kronFluxInner->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_INNER");
    165                 det->kronFluxOuter->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_OUTER");
    166                 det->diffRP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_P");
    167                 det->diffSnP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_P");
    168                 det->diffRM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_M");
    169                 det->diffSnM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_M");
    170                 det->flags2->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS2");
    171                 det->ippIdet->data.U32[numGood] = psMetadataLookupU32(NULL, row, "IPP_IDET");
    172                 det->nFrames->data.U16[numGood] = psMetadataLookupU16(NULL, row, "N_FRAMES");
    173                 det->padding->data.S16[numGood] = psMetadataLookupS16(NULL, row, "PADDING");
    174               } else {
    175                 det->psfInstFlux->data.F32[numGood] = NAN;
    176                 det->psfInstFluxSig->data.F32[numGood] = NAN;
    177                 det->apMag->data.F32[numGood] = NAN;
    178                 det->apMagRaw->data.F32[numGood] = NAN;
    179                 det->apMagRadius->data.F32[numGood] = NAN;
    180                 det->apFlux->data.F32[numGood] = NAN;
    181                 det->apFluxSig->data.F32[numGood] = NAN;
    182                 det->peakFluxAsMag->data.F32[numGood] = NAN;
    183                 det->calPsfMag->data.F32[numGood] = NAN;
    184                 det->calPsfMagSig->data.F32[numGood] = NAN;
    185                 det->sky->data.F32[numGood] = NAN;
    186                 det->skySig->data.F32[numGood] = NAN;
    187                 det->qualityPerfect->data.F32[numGood] = NAN;
    188                 det->momentsR1->data.F32[numGood] = NAN;
    189                 det->momentsRH->data.F32[numGood] = NAN;
    190                 det->kronFlux->data.F32[numGood] = NAN;
    191                 det->kronFluxErr->data.F32[numGood] = NAN;
    192                 det->kronFluxInner->data.F32[numGood] = NAN;
    193                 det->kronFluxOuter->data.F32[numGood] = NAN;
    194                 det->diffRP->data.F32[numGood] = NAN;
    195                 det->diffSnP->data.F32[numGood] = NAN;
    196                 det->diffRM->data.F32[numGood] = NAN;
    197                 det->diffSnM->data.F32[numGood] = NAN;
    198                 det->flags2->data.U32[numGood] = 0;
    199                 det->ippIdet->data.U32[numGood] = 0;
    200                 det->nFrames->data.U16[numGood] = 0;
    201                 det->padding->data.S16[numGood] = 0;
    202               }
    203             }
    204 
    205158            // Calculate error in RA, Dec
    206             double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32
    207             double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32
    208             double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); //SC: Warning! Promotion of F32
    209             double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); //SC: Warning! Promotion of F32
    210 
    211             if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
    212                 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
    213                 !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]) ||
    214168                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    215169                psTrace("ppMops.read", 10,
    216170                        "Discarding row %ld from input %d because of non-finite values: "
    217171                        "%f %f %lf %lf %f %f %f %f %f %f",
    218                         j, i,
    219                         det->x->data.F32[numGood], det->y->data.F32[numGood],
    220                         det->ra->data.F64[numGood], det->dec->data.F64[numGood],
    221                         det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
     172                        row, i,
     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],
    222176                        xErr, yErr, scale, angle);
     177                det->mask->data.U8[row] = 0xFF;
    223178                continue;
    224179            }
    … …  
    229184            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    230185            double errScale = scale / 3600.0;
    231             det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    232             det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    233 
    234             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;
    235190            plateScale += scale;
    236191            numGood++;
    … …  
    238193        det->seeing *= ((float) plateScale) / ((float) numGood);
    239194
    240         det->x->n = numGood;
    241         det->y->n = numGood;
    242         det->ra->n = numGood;
    243         det->dec->n = numGood;
    244         det->raErr->n = numGood;
    245         det->decErr->n = numGood;
    246         det->mag->n = numGood;
    247         det->magErr->n = numGood;
    248         det->chi2->n = numGood;
    249         det->dof->n = numGood;
    250         det->cr->n = numGood;
    251         det->extended->n = numGood;
    252         det->psfMajor->n = numGood;
    253         det->psfMinor->n = numGood;
    254         det->psfTheta->n = numGood;
    255         det->quality->n = numGood;
    256         det->numPix->n = numGood;
    257         det->xxMoment->n = numGood;
    258         det->xyMoment->n = numGood;
    259         det->yyMoment->n = numGood;
    260         det->flags->n = numGood;
    261         det->diffSkyfileId->n = numGood;
    262         det->naxis1->n = numGood;
    263         det->naxis2->n = numGood;
    264         det->mask->n = numGood;
    265         det->nPos->n = numGood;
    266         det->fPos->n = numGood;
    267         det->nRatioBad->n = numGood;
    268         det->nRatioMask->n = numGood;
    269         det->nRatioAll->n = numGood;
    270         det->psfInstFlux->n = numGood;
    271         det->psfInstFluxSig->n = numGood;
    272         det->apMag->n = numGood;
    273         det->apMagRaw->n = numGood;
    274         det->apMagRadius->n = numGood;
    275         det->apFlux->n = numGood;
    276         det->apFluxSig->n = numGood;
    277         det->peakFluxAsMag->n = numGood;
    278         det->calPsfMag->n = numGood;
    279         det->calPsfMagSig->n = numGood;
    280         det->sky->n = numGood;
    281         det->skySig->n = numGood;
    282         det->qualityPerfect->n = numGood;
    283         det->momentsR1->n = numGood;
    284         det->momentsRH->n = numGood;
    285         det->kronFlux->n = numGood;
    286         det->kronFluxErr->n = numGood;
    287         det->kronFluxInner->n = numGood;
    288         det->kronFluxOuter->n = numGood;
    289         det->diffRP->n = numGood;
    290         det->diffSnP->n = numGood;
    291         det->diffRM->n = numGood;
    292         det->diffSnM->n = numGood;
    293         det->flags2->n = numGood;
    294         det->ippIdet->n = numGood;
    295         det->nFrames->n = numGood;
    296         det->padding->n = numGood;
    297 
    298         det->num = numGood;
     195        // Are we using numGood for anything outside of this function?
     196        det->numGood = numGood;
    299197
    300198        if (isfinite(args->zp) && numGood > 0) {
    301             psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    302         }
    303 
    304         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    305 
    306         psFree(table);
    307         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);
    308203    }
    309204

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMopsWrite.c

    r29567 r32630  
    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
    38     psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
    39     psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
    40     psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
    41     psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
    42     psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
    43     psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
    44     psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
    45     psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
    46     psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     40    // Get these header words from the first input non null input
     41    ppMopsDetections *det = NULL;
     42    for (int d = 0; d < detections->n && det == NULL; d++) {
     43        det = detections->data[d];
     44    }
     45    if (det != NULL) {
     46        psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
     47        psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
     48        psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
     49        psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
     50        psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
     51        psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
     52        psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
     53        psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
     54        psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     55        psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
     56    } else {
     57        psWarning("no inputs with surviving detections. output header will be incomplete");
     58    }
    4759    psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    48     psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
    4960    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
    5061    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
    5162    psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
    52     //New field in header that tells about the CMF version
     63
     64    //field in header that tells about the CMF version
    5365    char cmfVersion[8];
    5466    sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
    5567    psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
    5668
    57     if (det->num == 0) {
     69    // Find the total number of detections
     70
     71    long total = 0;
     72    for (long i=0; i<detections->n; i++) {
     73        ppMopsDetections *det = detections->data[i];
     74        if (!det) {
     75            continue;
     76        }
     77        total += det->num;
     78    }
     79
     80    psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
     81
     82    if (total == 0) {
    5883        // Write dummy table
    5984        psMetadata *row = psMetadataAlloc(); // Output row
    … …  
    85110        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
    86111
    87             if (args->version == 2) {
     112        if (args->version == 2) {
    88113              // Write data of version 2 (see ICD)
    89114              psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
    … …  
    141166              psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
    142167                             0);
    143             }
     168        }
    144169
    145170        if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
    … …  
    151176        psFree(row);
    152177    } 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);
     178
     179#define addColumn(_outName, _inName, _convertTo32) \
     180        if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
     181            psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
     182            return false; \
     183        }
     184
     185        // Allocate the output table
     186        psMetadata *table = psMetadataAlloc();
     187        addColumn("RA", "RA_PSF", 0);
     188        addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     189        addColumn("DEC", "DEC_PSF", 0);
     190        addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     191        addColumn("MAG", "PSF_INST_MAG", 0);
     192        addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
     193        addColumn("PSF_CHI2", "PSF_CHISQ", 0);
     194        addColumn("PSF_DOF", "PSF_NDOF", 1);
     195        addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
     196        addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
     197        addColumn("PSF_MAJOR", NULL, 0);
     198        addColumn("PSF_MINOR", NULL, 0);     
     199        addColumn("PSF_THETA", NULL, 0);   
     200        addColumn("PSF_QUALITY", "PSF_QF", 0);
     201        addColumn("PSF_NPIX", NULL, 1);       
     202        addColumn("MOMENTS_XX", NULL, 0);
     203        addColumn("MOMENTS_XY", NULL, 0);
     204        addColumn("MOMENTS_YY", NULL, 0);
     205        addColumn("N_POS", "DIFF_NPOS", 1);
     206        addColumn("F_POS", "DIFF_FRATIO", 0);
     207        addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
     208        addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
     209        addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
     210        addColumn("FLAGS", "FLAGS", 1);
     211        addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
     212        if (args->version == 2) {
     213            addColumn("IPP_IDET", NULL, 1);
     214            addColumn("PSF_INST_FLUX", NULL, 0);
     215            addColumn("PSF_INST_FLUX_SIG", NULL, 0);
     216            addColumn("AP_MAG", NULL, 0);
     217            addColumn("AP_MAG_RAW", NULL, 0);
     218            addColumn("AP_MAG_RADIUS", NULL, 0);
     219            addColumn("AP_FLUX", NULL, 0);
     220            addColumn("AP_FLUX_SIG", NULL, 0);
     221            addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
     222            addColumn("CAL_PSF_MAG", NULL, 0);
     223            addColumn("CAL_PSF_MAG_SIG", NULL, 0);
     224            addColumn("SKY", NULL, 0);
     225            addColumn("SKY_SIGMA", NULL, 0);
     226            addColumn("PSF_QF_PERFECT", NULL, 0);
     227            addColumn("MOMENTS_R1", NULL, 0);
     228            addColumn("MOMENTS_RH", NULL, 0);
     229            addColumn("KRON_FLUX", NULL, 0);
     230            addColumn("KRON_FLUX_ERR", NULL, 0);
     231            addColumn("KRON_FLUX_INNER", NULL, 0);
     232            addColumn("KRON_FLUX_OUTER", NULL, 0);
     233            addColumn("DIFF_R_P", NULL, 0);
     234            addColumn("DIFF_SN_P", NULL, 0);
     235            addColumn("DIFF_R_M", NULL, 0);
     236            addColumn("DIFF_SN_M", NULL, 0);
     237            addColumn("FLAGS2", NULL, 1);
     238            addColumn("IPP_IDET", NULL, 0);
     239            addColumn("N_FRAMES", NULL, 0);
     240            addColumn("PADDING", NULL, 0);
     241        }
     242        if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
     243            psError(psErrorCodeLast(), false, "Unable to write table");
    263244            return false;
    264245        }
    … …  
    273254    return true;
    274255}
     256
     257static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
     258{
     259    if (inColumnName == NULL) {
     260        inColumnName = outColumnName;
     261    }
     262
     263    psVector *out = NULL;
     264    if (convertTo32) {
     265        // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
     266        // We want to write 32 bits to the output.
     267        int next = 0;
     268        for (long i=0; i<detections->n; i++) {
     269            ppMopsDetections *det = detections->data[i];
     270            if (!det || det->num == 0) {
     271                // no detections survived for this input
     272                continue;
     273            }
     274            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     275            if (!in) {
     276                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     277                return false;
     278            }
     279            if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
     280                psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
     281                    inColumnName, in->type.type);
     282                return false;
     283            }
     284            if (out == NULL) {
     285                // First time through set up the output vector and the copy parameters
     286                if (in->type.type == PS_TYPE_S64) {
     287                    out = psVectorAlloc(outputSize, PS_TYPE_S32);
     288                } else {
     289                    out = psVectorAlloc(outputSize, PS_TYPE_U32);
     290                }
     291            }
     292            for (long d=0; d < det->num; d++) {
     293                if (in->type.type == PS_TYPE_S64) {
     294                    out->data.S32[next++] = in->data.S64[d];
     295                } else {
     296                    out->data.U32[next++] = in->data.U64[d];
     297                }
     298            }
     299        }
     300    } else {
     301        void *next = NULL;
     302        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
     303        for (long i=0; i<detections->n; i++) {
     304            ppMopsDetections *det = detections->data[i];
     305            if (!det || det->num == 0) {
     306                // no detections survived for this input
     307                continue;
     308            }
     309            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     310            if (!in) {
     311                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     312                return false;
     313            }
     314            if (out == NULL) {
     315                // First time through set up the output vector and the copy parameters
     316                out = psVectorAlloc(outputSize, in->type.type);
     317                next = (void *) out->data.U8;
     318                switch (in->type.type) {
     319                    case PS_TYPE_S8:
     320                    elementSize = sizeof(psS8);
     321                    break;
     322                case PS_TYPE_U8:
     323                    elementSize = sizeof(psU8);
     324                    break;
     325                case PS_TYPE_S16:
     326                    elementSize = sizeof(psS16);
     327                    break;
     328                case PS_TYPE_U16:
     329                    elementSize = sizeof(psU16);
     330                    break;
     331                case PS_TYPE_S32:
     332                    elementSize = sizeof(psS32);
     333                    break;
     334                case PS_TYPE_U32:
     335                    elementSize = sizeof(psU32);
     336                    break;
     337                case PS_TYPE_S64:
     338                    elementSize = sizeof(psS64);
     339                    break;
     340                case PS_TYPE_U64:
     341                    elementSize = sizeof(psU64);
     342                    break;
     343                case PS_TYPE_F32:
     344                    elementSize = sizeof(psF32);
     345                    break;
     346                case PS_TYPE_F64:
     347                    elementSize = sizeof(psF64);
     348                    break;
     349                default:
     350                    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
     351                    return false;
     352
     353                }
     354            }
     355            // We are doing nasty things here so we can use memcpy.
     356            // It would be safer to do a proper loop over the elements.
     357            long toCopy = det->num * elementSize;
     358            memcpy(next, in->data.U8, toCopy);
     359            next += toCopy;
     360        }
     361    }
     362
     363    // Finally add the new column to the output table
     364    psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     365    psFree(out);    // drop reference
     366
     367    return true;
     368}
     369static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
     370{
     371    psVector *out = psVectorAlloc(total, PS_TYPE_S64);
     372    long next = 0;
     373    for (long i = 0; i<detections->n; i++) {
     374        ppMopsDetections *det = detections->data[i];
     375        if (!det) {
     376            continue;
     377        }
     378        psS64 diffSkyfileId = det->diffSkyfileId;
     379        for (long j = 0; j < det->num; j++) {
     380            out->data.S64[next++] = diffSkyfileId;
     381        }
     382    }
     383    psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
     384    psFree(out);
     385    return true;
     386}
Note: See TracChangeset for help on using the changeset viewer.