IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28623


Ignore:
Timestamp:
Jul 7, 2010, 10:44:29 AM (16 years ago)
Author:
Paul Price
Message:

Reverse merging to the 'new old version' of ppMops (r28043), since this is what MOPS wants right now. The 'new version' is at /branches/ppTranslate/

Location:
trunk/ppTranslate/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMops.c

    r28212 r28623  
    2020    }
    2121
    22     if (!ppMopsPurgeDuplicates(detections)) {
     22    ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections
     23    psFree(detections);
     24    if (!merged) {
    2325        psErrorStackPrint(stderr, "Unable to merge detections");
    2426        exit(PS_EXIT_SYS_ERROR);
    2527    }
    2628
    27     if (!ppMopsWrite(detections, args)) {
     29    if (!ppMopsWrite(merged, args)) {
    2830        psErrorStackPrint(stderr, "Unable to write detections");
    2931        exit(PS_EXIT_SYS_ERROR);
    3032    }
    3133
    32     psFree(detections);
     34    psFree(merged);
    3335    psFree(args);
    3436

  • trunk/ppTranslate/src/ppMops.h

    r28243 r28623  
    77#define IN_EXTNAME "SkyChip.psf"        // Extension name for data in input
    88#define OBSERVATORY_CODE "F51"          // IAU Observatory Code
     9#define OUT_EXTNAME "MOPS_TRANSIENT_DETECTIONS" // Extension name for data in output
     10#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_SATURATED | \
     11                     PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SKY_FAILURE) // Flags to exclude
    912
    1013// Configuration data
    … …  
    2427} ppMopsArguments;
    2528
     29#if 0
     30TTYPE19 = 'PSF_CHISQ'          / label for field  19
     31TFORM19 = '1E      '           / data format of field: 4-byte REAL
     32TTYPE20 = 'CR_NSIGMA'          / label for field  20
     33TFORM20 = '1E      '           / data format of field: 4-byte REAL
     34TTYPE21 = 'EXT_NSIGMA'         / label for field  21
     35TFORM21 = '1E      '           / data format of field: 4-byte REAL
     36TTYPE22 = 'PSF_MAJOR'          / label for field  22
     37TFORM22 = '1E      '           / data format of field: 4-byte REAL
     38TTYPE23 = 'PSF_MINOR'          / label for field  23
     39TFORM23 = '1E      '           / data format of field: 4-byte REAL
     40TTYPE24 = 'PSF_THETA'          / label for field  24
     41TFORM24 = '1E      '           / data format of field: 4-byte REAL
     42TTYPE25 = 'PSF_QF  '           / label for field  25
     43TFORM25 = '1E      '           / data format of field: 4-byte REAL
     44TTYPE26 = 'PSF_NDOF'           / label for field  26
     45TFORM26 = '1J      '           / data format of field: 4-byte INTEGER
     46TTYPE27 = 'PSF_NPIX'           / label for field  27
     47TFORM27 = '1J      '           / data format of field: 4-byte INTEGER
     48TTYPE28 = 'MOMENTS_XX'         / label for field  28
     49TFORM28 = '1E      '           / data format of field: 4-byte REAL
     50TTYPE29 = 'MOMENTS_XY'         / label for field  29
     51TFORM29 = '1E      '           / data format of field: 4-byte REAL
     52TTYPE30 = 'MOMENTS_YY'         / label for field  30
     53TFORM30 = '1E      '           / data format of field: 4-byte REAL
     54#endif
     55
    2656/// Parse arguments
    2757ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]);
    2858
    2959typedef struct {
    30     psString component;                 // Component name
    31     psMetadata *header;                 // FITS header (extension names *.hdr)
    3260    psString raBoresight, decBoresight; // RA,Dec of telescope boresight
    3361    psString filter;                    // Filter for exposure
    … …  
    3664    double posangle;                    // Position angle
    3765    double alt, az;                     // Telescope altitude and azimuth
    38     double mjd;                         // Modified Julian Date of exposure mid-point
     66    double mjd;                         // Modified Julian Date
    3967    float seeing;                       // Seeing of exposure
    40     int naxis1, naxis2;                 // Size of image
    4168    long num;                           // Number of detections
    42     psMetadata *psfHeader;              // FITS header (extension names *.psf)
    43     psMetadata *table;                  // Columns of data
    4469    psVector *x, *y;                    // Image coordinates
    4570    psVector *ra, *dec;                 // Sky coordinates
    46     psMetadata *deteffHeader;           // Detection efficiency header (extension names *.deteff)
    47     psMetadata *deteffTable;            // Detection efficiency table
     71    psVector *raErr, *decErr;           // Error in sky coordinates
     72    psVector *mag, *magErr;             // Magnitude and associated error
     73    psVector *chi2, *dof;               // Chi^2 from fitting, with associated degrees of freedom
     74    psVector *cr, *extended;            // Measures of CR-ness and extendedness
     75    psVector *psfMajor, *psfMinor, *psfTheta; // PSF major and minor axes, and position angle
     76    psVector *quality, *numPix;               // PSF quality factor and number of pixels
     77    psVector *xxMoment, *xyMoment, *yyMoment; // Moments
     78    psVector *flags;                    // psphot flags
     79    psVector *diffSkyfileId;            // Identifier for source image
     80    psVector *naxis1, *naxis2;          // Size of image
     81    psVector *mask;                     // Mask for detections
     82    psVector *nPos;                     // Number of positive pixels
     83    psVector *fPos;                     // Fraction of positive flux
     84    psVector *nRatioBad;                // Fraction of positive pixels to negative
     85    psVector *nRatioMask;               // Fraction of positive pixels to masked
     86    psVector *nRatioAll;                // Fraction of positive pixels to all
    4887} ppMopsDetections;
    4988
    50 // Allocator
    51 ppMopsDetections *ppMopsDetectionsAlloc(void);
     89
     90ppMopsDetections *ppMopsDetectionsAlloc(long num);
     91
     92/// Copy a detection
     93bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index);
     94
     95/// Purge the detections list of masked detections
     96bool ppMopsDetectionsPurge(ppMopsDetections *detections);
     97
    5298
    5399/// Read detections
    54100psArray *ppMopsRead(const ppMopsArguments *args);
    55101
    56 /// Purge duplicate detections
    57 bool ppMopsPurgeDuplicates(const psArray *detections);
     102/// Merge detections
     103ppMopsDetections *ppMopsMerge(const psArray *detections);
    58104
    59105/// Write detections
    60 bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args);
     106bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args);
    61107
    62108#endif

  • trunk/ppTranslate/src/ppMopsDetections.c

    r28243 r28623  
    1010static void mopsDetectionsFree(ppMopsDetections *det)
    1111{
    12     psFree(det->component);
    1312    psFree(det->raBoresight);
    1413    psFree(det->decBoresight);
    1514    psFree(det->filter);
    16     psFree(det->header);
    17     psFree(det->table);
    1815    psFree(det->x);
    1916    psFree(det->y);
    2017    psFree(det->ra);
    2118    psFree(det->dec);
    22     psFree(det->deteffHeader);
    23     psFree(det->deteffTable);
     19    psFree(det->raErr);
     20    psFree(det->decErr);
     21    psFree(det->mag);
     22    psFree(det->magErr);
     23    psFree(det->chi2);
     24    psFree(det->dof);
     25    psFree(det->cr);
     26    psFree(det->extended);
     27    psFree(det->psfMajor);
     28    psFree(det->psfMinor);
     29    psFree(det->psfTheta);
     30    psFree(det->quality);
     31    psFree(det->numPix);
     32    psFree(det->xxMoment);
     33    psFree(det->xyMoment);
     34    psFree(det->yyMoment);
     35    psFree(det->flags);
     36    psFree(det->diffSkyfileId);
     37    psFree(det->naxis1);
     38    psFree(det->naxis2);
     39    psFree(det->mask);
     40    psFree(det->nPos);
     41    psFree(det->fPos);
     42    psFree(det->nRatioBad);
     43    psFree(det->nRatioMask);
     44    psFree(det->nRatioAll);
    2445
    2546    return;
    2647}
    2748
    28 ppMopsDetections *ppMopsDetectionsAlloc(void)
     49ppMopsDetections *ppMopsDetectionsAlloc(long num)
    2950{
    3051    ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return
    3152    psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree);
    3253
    33     det->component = NULL;
    3454    det->raBoresight = NULL;
    3555    det->decBoresight = NULL;
    … …  
    4363    det->seeing = NAN;
    4464    det->num = 0;
    45     det->header = NULL;
    46     det->table = NULL;
    47     det->x = NULL;
    48     det->y = NULL;
    49     det->ra = NULL;
    50     det->dec = NULL;
    51     det->deteffHeader = NULL;
    52     det->deteffTable = NULL;
     65    det->x = psVectorAllocEmpty(num, PS_TYPE_F32);
     66    det->y = psVectorAllocEmpty(num, PS_TYPE_F32);
     67    det->ra = psVectorAllocEmpty(num, PS_TYPE_F64);
     68    det->dec = psVectorAllocEmpty(num, PS_TYPE_F64);
     69    det->raErr = psVectorAllocEmpty(num, PS_TYPE_F64);
     70    det->decErr = psVectorAllocEmpty(num, PS_TYPE_F64);
     71    det->mag = psVectorAllocEmpty(num, PS_TYPE_F32);
     72    det->magErr = psVectorAllocEmpty(num, PS_TYPE_F32);
     73    det->chi2 = psVectorAllocEmpty(num, PS_TYPE_F32);
     74    det->dof = psVectorAllocEmpty(num, PS_TYPE_S32);
     75    det->cr = psVectorAllocEmpty(num, PS_TYPE_F32);
     76    det->extended = psVectorAllocEmpty(num, PS_TYPE_F32);
     77    det->psfMajor = psVectorAllocEmpty(num, PS_TYPE_F32);
     78    det->psfMinor = psVectorAllocEmpty(num, PS_TYPE_F32);
     79    det->psfTheta = psVectorAllocEmpty(num, PS_TYPE_F32);
     80    det->quality = psVectorAllocEmpty(num, PS_TYPE_F32);
     81    det->numPix = psVectorAllocEmpty(num, PS_TYPE_S32);
     82    det->xxMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
     83    det->xyMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
     84    det->yyMoment = psVectorAllocEmpty(num, PS_TYPE_F32);
     85    det->flags = psVectorAllocEmpty(num, PS_TYPE_U32);
     86    det->diffSkyfileId = psVectorAllocEmpty(num, PS_TYPE_S64);
     87    det->naxis1 = psVectorAllocEmpty(num, PS_TYPE_S32);
     88    det->naxis2 = psVectorAllocEmpty(num, PS_TYPE_S32);
     89    det->mask = psVectorAllocEmpty(num, PS_TYPE_U8);
     90    det->nPos = psVectorAllocEmpty(num, PS_TYPE_S32);
     91    det->fPos = psVectorAllocEmpty(num, PS_TYPE_F32);
     92    det->nRatioBad = psVectorAllocEmpty(num, PS_TYPE_F32);
     93    det->nRatioMask = psVectorAllocEmpty(num, PS_TYPE_F32);
     94    det->nRatioAll = psVectorAllocEmpty(num, PS_TYPE_F32);
    5395
    5496    return det;
    5597}
     98
     99
     100ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num)
     101{
     102    det->x = psVectorRealloc(det->x, num);
     103    det->y = psVectorRealloc(det->y, num);
     104    det->ra = psVectorRealloc(det->ra, num);
     105    det->dec = psVectorRealloc(det->dec, num);
     106    det->raErr = psVectorRealloc(det->raErr, num);
     107    det->decErr = psVectorRealloc(det->decErr, num);
     108    det->mag = psVectorRealloc(det->mag, num);
     109    det->magErr = psVectorRealloc(det->magErr, num);
     110    det->chi2 = psVectorRealloc(det->chi2, num);
     111    det->dof = psVectorRealloc(det->dof, num);
     112    det->cr = psVectorRealloc(det->cr, num);
     113    det->extended = psVectorRealloc(det->extended, num);
     114    det->psfMajor = psVectorRealloc(det->psfMajor, num);
     115    det->psfMinor = psVectorRealloc(det->psfMinor, num);
     116    det->psfTheta = psVectorRealloc(det->psfTheta, num);
     117    det->quality = psVectorRealloc(det->quality, num);
     118    det->numPix = psVectorRealloc(det->numPix, num);
     119    det->xxMoment = psVectorRealloc(det->xxMoment, num);
     120    det->xyMoment = psVectorRealloc(det->xyMoment, num);
     121    det->yyMoment = psVectorRealloc(det->yyMoment, num);
     122    det->flags = psVectorRealloc(det->flags, num);
     123    det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num);
     124    det->naxis1 = psVectorRealloc(det->naxis1, num);
     125    det->naxis2 = psVectorRealloc(det->naxis2, num);
     126    det->mask = psVectorRealloc(det->mask, num);
     127    det->nPos = psVectorRealloc(det->nPos, num);
     128    det->fPos = psVectorRealloc(det->fPos, num);
     129    det->nRatioBad = psVectorRealloc(det->nRatioBad, num);
     130    det->nRatioMask = psVectorRealloc(det->nRatioMask, num);
     131    det->nRatioAll = psVectorRealloc(det->nRatioAll, num);
     132
     133    return det;
     134}
     135
     136
     137bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec,
     138                         double raErr, double decErr, float mag, float magErr,
     139                         float chi2, int dof, float cr, float extended, float psfMajor,
     140                         float psfMinor, float psfTheta, float quality, int numPix,
     141                         float xxMoment, float xyMoment, float yyMoment,
     142                         psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2,
     143                         int nPos, float fPos, float nRatioBad, float nRatioMask, float nRatioAll)
     144{
     145    psVectorAppend(det->x, x);
     146    psVectorAppend(det->y, y);
     147    psVectorAppend(det->ra, ra);
     148    psVectorAppend(det->dec, dec);
     149    psVectorAppend(det->raErr, raErr);
     150    psVectorAppend(det->decErr, decErr);
     151    psVectorAppend(det->mag, mag);
     152    psVectorAppend(det->magErr, magErr);
     153    psVectorAppend(det->chi2, chi2);
     154    psVectorAppend(det->dof, dof);
     155    psVectorAppend(det->cr, cr);
     156    psVectorAppend(det->extended, extended);
     157    psVectorAppend(det->psfMajor, psfMajor);
     158    psVectorAppend(det->psfMinor, psfMinor);
     159    psVectorAppend(det->psfTheta, psfTheta);
     160    psVectorAppend(det->quality, quality);
     161    psVectorAppend(det->numPix, numPix);
     162    psVectorAppend(det->xxMoment, xxMoment);
     163    psVectorAppend(det->xyMoment, xyMoment);
     164    psVectorAppend(det->yyMoment, yyMoment);
     165    psVectorAppend(det->flags, flags);
     166    psVectorAppend(det->diffSkyfileId, diffSkyfileId);
     167    psVectorAppend(det->naxis1, naxis1);
     168    psVectorAppend(det->naxis2, naxis2);
     169    psVectorAppend(det->mask, 0);
     170    psVectorAppend(det->nPos, nPos);
     171    psVectorAppend(det->fPos, fPos);
     172    psVectorAppend(det->nRatioBad, nRatioBad);
     173    psVectorAppend(det->nRatioMask, nRatioMask);
     174    psVectorAppend(det->nRatioAll, nRatioAll);
     175
     176    return true;
     177}
     178
     179
     180bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index)
     181{
     182    psVectorAppend(target->x, source->x->data.F32[index]);
     183    psVectorAppend(target->y, source->y->data.F32[index]);
     184    psVectorAppend(target->ra, source->ra->data.F64[index]);
     185    psVectorAppend(target->dec, source->dec->data.F64[index]);
     186    psVectorAppend(target->raErr, source->raErr->data.F64[index]);
     187    psVectorAppend(target->decErr, source->decErr->data.F64[index]);
     188    psVectorAppend(target->mag, source->mag->data.F32[index]);
     189    psVectorAppend(target->magErr, source->magErr->data.F32[index]);
     190    psVectorAppend(target->chi2, source->chi2->data.F32[index]);
     191    psVectorAppend(target->dof, source->dof->data.S32[index]);
     192    psVectorAppend(target->cr, source->cr->data.F32[index]);
     193    psVectorAppend(target->extended, source->extended->data.F32[index]);
     194    psVectorAppend(target->psfMajor, source->psfMajor->data.F32[index]);
     195    psVectorAppend(target->psfMinor, source->psfMinor->data.F32[index]);
     196    psVectorAppend(target->psfTheta, source->psfTheta->data.F32[index]);
     197    psVectorAppend(target->quality, source->quality->data.F32[index]);
     198    psVectorAppend(target->numPix, source->numPix->data.S32[index]);
     199    psVectorAppend(target->xxMoment, source->xxMoment->data.F32[index]);
     200    psVectorAppend(target->xyMoment, source->xyMoment->data.F32[index]);
     201    psVectorAppend(target->yyMoment, source->yyMoment->data.F32[index]);
     202    psVectorAppend(target->flags, source->flags->data.U32[index]);
     203    psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]);
     204    psVectorAppend(target->naxis1, source->naxis1->data.S32[index]);
     205    psVectorAppend(target->naxis2, source->naxis2->data.S32[index]);
     206    psVectorAppend(target->mask, 0);
     207    psVectorAppend(target->nPos, source->nPos->data.S32[index]);
     208    psVectorAppend(target->fPos, source->fPos->data.F32[index]);
     209    psVectorAppend(target->nRatioBad, source->nRatioBad->data.F32[index]);
     210    psVectorAppend(target->nRatioMask, source->nRatioMask->data.F32[index]);
     211    psVectorAppend(target->nRatioAll, source->nRatioAll->data.F32[index]);
     212
     213    target->num++;
     214
     215    return true;
     216}
     217
     218
     219bool ppMopsDetectionsPurge(ppMopsDetections *det)
     220{
     221    long num = 0;
     222    for (long i = 0; i < det->num; i++) {
     223        if (!det->mask->data.U8[i]) {
     224            if (i == num) {
     225                // No need to copy
     226                num++;
     227                continue;
     228            }
     229            det->x->data.F32[num] = det->x->data.F32[i];
     230            det->y->data.F32[num] = det->y->data.F32[i];
     231            det->ra->data.F64[num] = det->ra->data.F64[i];
     232            det->dec->data.F64[num] = det->dec->data.F64[i];
     233            det->raErr->data.F64[num] = det->raErr->data.F64[i];
     234            det->decErr->data.F64[num] = det->decErr->data.F64[i];
     235            det->mag->data.F32[num] = det->mag->data.F32[i];
     236            det->magErr->data.F32[num] = det->magErr->data.F32[i];
     237            det->chi2->data.F32[num] = det->chi2->data.F32[i];
     238            det->dof->data.S32[num] = det->dof->data.S32[i];
     239            det->cr->data.F32[num] = det->cr->data.F32[i];
     240            det->extended->data.F32[num] = det->extended->data.F32[i];
     241            det->psfMajor->data.F32[num] = det->psfMajor->data.F32[i];
     242            det->psfMinor->data.F32[num] = det->psfMinor->data.F32[i];
     243            det->psfTheta->data.F32[num] = det->psfTheta->data.F32[i];
     244            det->quality->data.F32[num] = det->quality->data.F32[i];
     245            det->numPix->data.S32[num] = det->numPix->data.S32[i];
     246            det->xxMoment->data.F32[num] = det->xxMoment->data.F32[i];
     247            det->xyMoment->data.F32[num] = det->xyMoment->data.F32[i];
     248            det->yyMoment->data.F32[num] = det->yyMoment->data.F32[i];
     249            det->flags->data.U32[num] = det->flags->data.U32[i];
     250            det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i];
     251            det->naxis1->data.S32[num] = det->naxis1->data.S32[i];
     252            det->naxis2->data.S32[num] = det->naxis2->data.S32[i];
     253            det->mask->data.U8[num] = 0;
     254            det->nPos->data.S32[num] = det->nPos->data.S32[i];
     255            det->fPos->data.F32[num] = det->fPos->data.F32[i];
     256            det->nRatioBad->data.F32[num] = det->nRatioBad->data.F32[i];
     257            det->nRatioMask->data.F32[num] = det->nRatioMask->data.F32[i];
     258            det->nRatioAll->data.F32[num] = det->nRatioAll->data.F32[i];
     259            num++;
     260        }
     261    }
     262    det->x->n = num;
     263    det->y->n = num;
     264    det->ra->n = num;
     265    det->dec->n = num;
     266    det->raErr->n = num;
     267    det->decErr->n = num;
     268    det->mag->n = num;
     269    det->magErr->n = num;
     270    det->chi2->n = num;
     271    det->dof->n = num;
     272    det->cr->n = num;
     273    det->extended->n = num;
     274    det->psfMajor->n = num;
     275    det->psfMinor->n = num;
     276    det->psfTheta->n = num;
     277    det->quality->n = num;
     278    det->numPix->n = num;
     279    det->xxMoment->n = num;
     280    det->xyMoment->n = num;
     281    det->yyMoment->n = num;
     282    det->flags->n = num;
     283    det->diffSkyfileId->n = num;
     284    det->naxis1->n = num;
     285    det->naxis2->n = num;
     286    det->mask->n = num;
     287    det->num = num;
     288    det->nPos->n = num;
     289    det->fPos->n = num;
     290    det->nRatioBad->n = num;
     291    det->nRatioMask->n = num;
     292    det->nRatioAll->n = num;
     293
     294    return true;
     295}
     296

  • trunk/ppTranslate/src/ppMopsMerge.c

    r28220 r28623  
    99#include "ppMops.h"
    1010
    11 #define LEAF_SIZE 2                     // Size of leaf
     11#define LEAF_SIZE 4                     // Size of leaf
    1212#define MATCH_RADIUS SEC_TO_RAD(1.0)    // Matching radius
    1313#define MJD_TOL 1.0/3600.0/24.0         // Tolerance for MJD matching
    … …  
    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 
    2719// Get distance from detection to centre of image
    2820static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
    … …  
    3022    )
    3123{
    32     float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
    33     float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
     24    float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0;
     25    float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0;
    3426    return PS_SQR(dx) + PS_SQR(dy);
    3527}
    3628
    3729
    38 bool ppMopsPurgeDuplicates(const psArray *detections)
     30ppMopsDetections *ppMopsMerge(const psArray *detections)
    3931{
    4032    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    4133
    42     int numInputs = detections->n;                // Number of inputs
    43     psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
     34    psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
    4435
    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++) {
     36    ppMopsDetections *merged = NULL;    // Merged list
     37    int num = 1;                                                         // Number of merged files
     38    for (int i = 0; i < detections->n; i++) {
    7339        ppMopsDetections *det = detections->data[i]; // Detections of interest
    7440        if (!det) {
    … …  
    7945            continue;
    8046        }
     47        num++;
     48        if (!merged) {
     49            psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
     50            merged = psMemIncrRefCounter(det);
     51            continue;
     52        }
     53        psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
    8154
    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             }
     55        // XXX compare exposure properties
     56        if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
     57            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     58                    merged->raBoresight, det->raBoresight);
     59            return NULL;
     60        }
     61        if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
     62            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     63                    merged->decBoresight, det->decBoresight);
     64            return NULL;
     65        }
     66        if (strcmp(merged->filter, det->filter) != 0) {
     67            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     68                    merged->filter, det->filter);
     69            return NULL;
    13870        }
    13971
    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;
     72        if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
     73            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     74                    merged->airmass, det->airmass);
     75            return NULL;
    14676        }
    147         num += det->num;
    148     }
     77        if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
     78            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     79                    merged->exptime, det->exptime);
     80            return NULL;
     81        }
     82        if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
     83            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     84                    merged->posangle, det->posangle);
     85            return NULL;
     86        }
     87        if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
     88            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     89                    merged->alt, det->alt);
     90            return NULL;
     91        }
     92        if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
     93            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     94                    merged->az, det->az);
     95            return NULL;
     96        }
     97        if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
     98            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     99                    merged->mjd, det->mjd);
     100            return NULL;
     101        }
    149102
    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     }
     103        merged->seeing += det->seeing;  // Taking average
    156104
    157     for (int i = 0; i < detections->n; i++) {
    158         ppMopsDetections *det = detections->data[i]; // Detections of interest
    159         if (!det) {
    160             continue;
    161         } else if (det->num == 0) {
    162             continue;
     105        psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
     106        if (!tree) {
     107            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     108            psFree(merged);
     109            return NULL;
    163110        }
    164         psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
    165111
    166         psVector *dupes = duplicates->data[i];            // Duplicates list
    167112        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
    168113        for (int j = 0; j < det->num; j++) {
    … …  
    174119                psFree(coords);
    175120                psFree(tree);
    176                 return false;
     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);
     126                psFree(indices);
     127                ppMopsDetectionsCopySingle(merged, det, j);
     128                continue;
    177129            }
    178130            psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
    179             psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
    180 
    181             if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) {
    182                 // It's myself
    183                 psFree(indices);
    184                 continue;
    185             }
    186131
    187132            // Which one do we keep?
    … …  
    189134            long bestIndex = -1;           // Index with best distance
    190135            for (int k = 0; k < indices->n; k++) {
    191                 long mergeIndex = indices->data.S64[k]; // Index of point in merged list
    192                 int source = sourceMerged->data.U16[mergeIndex]; // Source image
    193                 if (source == i) {
    194                     continue;
    195                 }
    196                 long index = indexMerged->data.U32[mergeIndex];  // Index in source
    197                 psVector *dupes = duplicates->data[source];     // Duplicates list
    198                 if (dupes->data.U8[index]) {
    199                     continue;
    200                 }
    201 
    202                 float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
     136                long index = indices->data.S64[k]; // Index of point
     137                float distance = mergeDistance(merged, index); // Distance to centre of image
    203138                if (distance < bestDistance) {
    204139                    bestDistance = distance;
    … …  
    207142            }
    208143
    209             if (bestIndex >= 0) {
    210                 float distance = mergeDistance(det, j); // Distance to centre of image
    211                 if (bestIndex >= 0 && distance < bestDistance) {
    212                     psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
    213                     for (int k = 0; k < indices->n; k++) {
    214                         long mergeIndex = indices->data.S64[k]; // Index of point
    215                         int source = sourceMerged->data.U16[mergeIndex]; // Source image
    216                         if (source == i) {
    217                             continue;
    218                         }
    219                         long index = indexMerged->data.U32[mergeIndex];  // Index in source
    220                         psVector *dupes = duplicates->data[source];      // Duplicates list
    221                         if (!dupes->data.U8[index]) {
    222                             dupes->data.U8[index] = 0xFF;
    223                             dupNum->data.U32[source]++;
    224                         }
    225                     }
    226                 } else {
    227                     psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    228                     dupes->data.U8[j] = 0xFF;
    229                     dupNum->data.U32[i]++;
     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;
    230151                }
     152                ppMopsDetectionsCopySingle(merged, det, j);
     153            } else {
     154                psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    231155            }
    232156            psFree(indices);
    233157        }
    234         psFree(coords);
    235     }
    236     psFree(tree);
    237     psFree(raMerged);
    238     psFree(decMerged);
    239     psFree(sourceMerged);
    240     psFree(indexMerged);
    241158
    242     // Remove duplicates
    243     for (int i = 0; i < detections->n; i++) {
    244         ppMopsDetections *det = detections->data[i]; // Detections of interest
    245         if (!det) {
    246             continue;
    247         } else if (det->num == 0) {
    248             continue;
    249         }
    250         psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
     159        psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
    251160
    252 #define VECTOR_PURGE_CASE(TYPE)                                      \
    253         case PS_TYPE_##TYPE:                                         \
    254           for (int i = 0, j = 0; i < vector->n; i++) {               \
    255               if (!dupes->data.U8[i]) {                              \
    256                   if (i == j) {                                      \
    257                       j++;                                           \
    258                       continue;                                      \
    259                   }                                                  \
    260                   vector->data.TYPE[j] = vector->data.TYPE[i];       \
    261               }                                                      \
    262           }                                                          \
    263           break;
    264 
    265         psVector *dupes = duplicates->data[i]; // Duplicates
    266         psArray *table = psListToArray(det->table->list); // Table of data
    267         for (int j = 0; j < table->n; j++) {
    268             psMetadataItem *item = table->data[j]; // Table item
    269             psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type);
    270             psVector *vector = item->data.V; // Vector to purge
    271             switch (vector->type.type) {
    272                 VECTOR_PURGE_CASE(U8);
    273                 VECTOR_PURGE_CASE(U16);
    274                 VECTOR_PURGE_CASE(U32);
    275                 VECTOR_PURGE_CASE(U64);
    276                 VECTOR_PURGE_CASE(S8);
    277                 VECTOR_PURGE_CASE(S16);
    278                 VECTOR_PURGE_CASE(S32);
    279                 VECTOR_PURGE_CASE(S64);
    280                 VECTOR_PURGE_CASE(F32);
    281                 VECTOR_PURGE_CASE(F64);
    282               default:
    283                 psAbort("Unrecognised vector type: %x", vector->type.type);
    284             }
    285         }
     161        psFree(tree);
     162        ppMopsDetectionsPurge(merged);
    286163    }
    287164
    288     return true;
     165    psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
     166
     167    merged->seeing /= num;
     168
     169    return merged;
    289170}
    290171

  • trunk/ppTranslate/src/ppMopsRead.c

    r28243 r28623  
    44
    55#include <stdio.h>
    6 #include <string.h>
    76#include <pslib.h>
    87
    … …  
    1716    psArray *detections = psArrayAlloc(num); // Array of detections, to return
    1817    for (int i = 0; i < num; i++) {
    19         const char *name = inNames->data[i]; // File name
    20         psFits *fits = psFitsOpen(name, "r"); // FITS file
     18        psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file
    2119        if (!fits) {
    2220            psError(PS_ERR_IO, false, "Unable to open input %d", i);
    … …  
    2624        if (!header) {
    2725            psError(PS_ERR_IO, false, "Unable to read header %d", i);
     26            return false;
     27        }
     28
     29        psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
     30        if (diffSkyfileId == 0) {
     31            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
    2832            return false;
    2933        }
    … …  
    4347            continue;
    4448        }
    45         ppMopsDetections *det = detections->data[i] = ppMopsDetectionsAlloc();
    46         det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
    47         det->header = header;
    48         det->num = size;
     49        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     50
     51        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    4952
    5053        det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
    … …  
    6164                             psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    6265
    63         det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1");
    64         det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2");
     66        int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     67        int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    6568
    66         det->psfHeader = psFitsReadHeader(NULL, fits);
    67         if (!det->psfHeader) {
    68             psError(psErrorCodeLast(), false, "Unable to read header %d", i);
     69        psFree(header);
     70
     71        psArray *table = psFitsReadTable(fits); // Table of interest
     72        if (!table) {
     73            psError(PS_ERR_IO, false, "Unable to read table %d", i);
    6974            return false;
    7075        }
    71         psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest
    72         if (!table) {
    73             psError(psErrorCodeLast(), false, "Unable to read table %d", i);
    74             return false;
     76        psFitsClose(fits);
     77
     78        double plateScale = 0.0;        // Plate scale
     79        long numGood = 0;               // Number of good rows
     80        for (long j = 0; j < size; j++) {
     81            psMetadata *row = table->data[j]; // Row of interest
     82
     83            psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
     84            if (flags & SOURCE_MASK) {
     85                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
     86                continue;
     87            }
     88
     89            det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
     90            det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
     91            det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
     92            det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
     93            det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
     94            det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
     95            det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");
     96            det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");
     97            det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");
     98            det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
     99            det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");
     100            det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");
     101            det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");
     102            det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");
     103            det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");
     104            det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");
     105            det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");
     106            det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");
     107            det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
     108            det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
     109            det->naxis1->data.S32[numGood] = naxis1;
     110            det->naxis2->data.S32[numGood] = naxis2;
     111
     112            det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");
     113            det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");
     114            det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");
     115            det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");
     116            det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");
     117
     118            // Calculate error in RA, Dec
     119            double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG");
     120            double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG");
     121            double scale = psMetadataLookupF64(NULL, row, "PLTSCALE");
     122            double angle = psMetadataLookupF64(NULL, row, "POSANGLE");
     123
     124            if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
     125                !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
     126                !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
     127                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
     128                psTrace("ppMops.read", 10,
     129                        "Discarding row %ld from input %d because of non-finite values: "
     130                        "%f %f %lf %lf %f %f %f %f %f %f",
     131                        j, i,
     132                        det->x->data.F32[numGood], det->y->data.F32[numGood],
     133                        det->ra->data.F64[numGood], det->dec->data.F64[numGood],
     134                        det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
     135                        xErr, yErr, scale, angle);
     136                continue;
     137            }
     138
     139            // XXX Not at all sure I've got the angles around the right way here...
     140            double cosAngle = cos(angle), sinAngle = sin(angle);
     141            double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
     142            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
     143            double errScale = scale / 3600.0;
     144            det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     145            det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     146
     147            det->mask->data.U8[numGood] = 0;
     148            plateScale += scale;
     149            numGood++;
     150        }
     151        det->seeing *= plateScale / numGood;
     152
     153        det->x->n = numGood;
     154        det->y->n = numGood;
     155        det->ra->n = numGood;
     156        det->dec->n = numGood;
     157        det->raErr->n = numGood;
     158        det->decErr->n = numGood;
     159        det->mag->n = numGood;
     160        det->magErr->n = numGood;
     161        det->chi2->n = numGood;
     162        det->dof->n = numGood;
     163        det->cr->n = numGood;
     164        det->extended->n = numGood;
     165        det->psfMajor->n = numGood;
     166        det->psfMinor->n = numGood;
     167        det->psfTheta->n = numGood;
     168        det->quality->n = numGood;
     169        det->numPix->n = numGood;
     170        det->xxMoment->n = numGood;
     171        det->xyMoment->n = numGood;
     172        det->yyMoment->n = numGood;
     173        det->flags->n = numGood;
     174        det->diffSkyfileId->n = numGood;
     175        det->naxis1->n = numGood;
     176        det->naxis2->n = numGood;
     177        det->mask->n = numGood;
     178        det->nPos->n = numGood;
     179        det->fPos->n = numGood;
     180        det->nRatioBad->n = numGood;
     181        det->nRatioMask->n = numGood;
     182        det->nRatioAll->n = numGood;
     183
     184        det->num = numGood;
     185
     186        if (isfinite(args->zp) && numGood > 0) {
     187            psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    75188        }
    76189
    77         psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    78         psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     190        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    79191
    80         det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    81         det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    82         det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
    83         det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
    84         if (!det->ra || !det->dec || !det->x || !det->y) {
    85             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
    86             return false;
    87         }
    88 
    89         psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    90 
    91         if (!psFitsMoveExtName(fits, "SkyChip.deteff")) {
    92             psWarning("No detection efficiencies included in %s", det->component);
    93             psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component);
    94             psErrorClear();
    95             continue;
    96         }
    97 
    98         det->deteffHeader = psFitsReadHeader(NULL, fits);
    99         if (!det->deteffHeader) {
    100             psWarning("Unable to read detection efficiency header in %s", det->component);
    101             psErrorClear();
    102             continue;
    103         }
    104         det->deteffTable = psFitsReadTableAllColumns(fits);
    105         if (!det->deteffTable) {
    106             psFree(det->deteffHeader);
    107             det->deteffHeader = NULL;
    108             psWarning("Unable to read detection efficiency table in %s", det->component);
    109             psErrorClear();
    110             continue;
    111         }
    112         psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component);
    113         psFitsClose(fits);
     192        psFree(table);
     193        detections->data[i] = det;
    114194    }
    115195

  • trunk/ppTranslate/src/ppMopsWrite.c

    r28264 r28623  
    99#include "ppTranslateVersion.h"
    1010
    11 bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
     11bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args)
    1212{
    13     psTrace("ppMops.write", 1, "Writing %ld extensions to %s", detections->n, args->output);
     13    psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output);
    1414
    1515    psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    … …  
    1919    }
    2020
    21     // Primary header
    22     {
    23         psMetadata *header = psMetadataAlloc(); // Header to write
    24         ppTranslateVersionHeader(header);
    25         psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
    26         psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
    27         psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
    28         psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
    29         psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
    30         psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
    31         psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
    32         psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
    33         psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
    34         psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
    35         psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
    36         psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    3721
    38         ppMopsDetections *det = detections->data[0]; // Representative set of detections
    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         psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
    42         psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass at boresight", det->airmass);
    43         psMetadataAddF32(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
    44         psMetadataAddF64(header, PS_LIST_TAIL, "POSANG", 0, "Position angle", det->posangle);
    45         psMetadataAddF64(header, PS_LIST_TAIL, "ALT", 0, "Altitude of boresight", det->alt);
    46         psMetadataAddF64(header, PS_LIST_TAIL, "AZ", 0, "Azimuth of boresight", det->az);
    47         psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of observation", det->mjd);
     22    psMetadata *header = psMetadataAlloc(); // Header to write
     23    psString source = ppTranslateSource(), version = ppTranslateVersion();
     24    psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
     25    psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
     26    ppTranslateVersionHeader(header);
     27    psFree(source);
     28    psFree(version);
    4829
    49         if (!psFitsWriteBlank(fits, header, NULL)) {
    50             psError(psErrorCodeLast(), false, "Unable to write primary header");
     30    psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
     31    psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
     32    psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
     33    psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
     34    psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
     35    psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
     36    psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
     37    psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
     38
     39    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
     40    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
     41    psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
     42    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
     43    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
     44    psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
     45    psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
     46    psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
     47    psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     48    psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
     49    psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
     50    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
     51    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
     52    psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
     53
     54    if (det->num == 0) {
     55        // Write dummy table
     56        psMetadata *row = psMetadataAlloc(); // Output row
     57        psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
     58        psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
     59        psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
     60        psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
     61        psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
     62        psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
     63        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
     64        psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
     65        psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
     66        psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
     67        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
     68        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
     69        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
     70        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
     71        psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
     72        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
     73        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
     74        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
     75        psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
     76        psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
     77
     78        psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
     79        psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
     80        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
     81        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
     82        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
     83        if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
     84            psErrorStackPrint(stderr, "Unable to write empty table.");
     85            psFree(header);
     86            psFree(row);
    5187            return false;
    5288        }
    53         psFree(header);
     89        psFree(row);
     90    } else {
     91        psArray *table = psArrayAlloc(det->num); // Table to write
     92        for (long i = 0; i < det->num; i++) {
     93            psMetadata *row = psMetadataAlloc(); // Output row
     94            psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)",
     95                             RAD_TO_DEG(det->ra->data.F64[i]));
     96            psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)",
     97                             det->raErr->data.F64[i]);
     98            psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)",
     99                             RAD_TO_DEG(det->dec->data.F64[i]));
     100            psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)",
     101                             det->decErr->data.F64[i]);
     102            psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]);
     103            psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]);
     104            psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", det->chi2->data.F32[i]);
     105            psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit",
     106                             det->dof->data.S32[i]);
     107            psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR",
     108                             det->cr->data.F32[i]);
     109            psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness",
     110                             det->extended->data.F32[i]);
     111            psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", det->psfMajor->data.F32[i]);
     112            psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", det->psfMinor->data.F32[i]);
     113            psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)",
     114                             det->psfTheta->data.F32[i]);
     115            psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor",
     116                             det->quality->data.F32[i]);
     117            psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF",
     118                             det->numPix->data.S32[i]);
     119            psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", det->xxMoment->data.F32[i]);
     120            psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", det->xyMoment->data.F32[i]);
     121            psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", det->yyMoment->data.F32[i]);
     122            psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels",
     123                             det->nPos->data.S32[i]);
     124            psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels",
     125                             det->fPos->data.F32[i]);
     126            psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative",
     127                             det->nRatioBad->data.F32[i]);
     128            psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked",
     129                             det->nRatioMask->data.F32[i]);
     130            psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all",
     131                             det->nRatioAll->data.F32[i]);
     132            psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]);
     133            psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile",
     134                             det->diffSkyfileId->data.S64[i]);
     135
     136            table->data[i] = row;
     137        }
     138        if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) {
     139            psErrorStackPrint(stderr, "Unable to write table.");
     140            psFree(header);
     141            psFree(table);
     142            return false;
     143        }
     144        psFree(table);
    54145    }
    55146
    56     for (int i = 0; i < detections->n; i++) {
    57         ppMopsDetections *det = detections->data[i]; // Detections for extension
    58         psTrace("ppMops.write", 1, "Writing extension %d to %s", i, args->output);
    59         psString hdrName = NULL, psfName = NULL, deteffName = NULL;
    60         psStringAppend(&hdrName, "%s.hdr", det->component);
    61         psStringAppend(&psfName, "%s.psf", det->component);
    62         psStringAppend(&deteffName, "%s.deteff", det->component);
    63 
    64         if (!psFitsWriteBlank(fits, det->header, hdrName)) {
    65             psError(psErrorCodeLast(), false, "Unable to write header %d", i);
    66             return false;
    67         }
    68         if (!psFitsWriteTableAllColumns(fits, det->psfHeader, det->table, psfName)) {
    69             psError(psErrorCodeLast(), false, "Unable to write table %d", i);
    70             return false;
    71         }
    72         if (det->deteffHeader && det->deteffTable &&
    73             !psFitsWriteTableAllColumns(fits, det->deteffHeader, det->deteffTable, deteffName)) {
    74             psError(psErrorCodeLast(), false, "Unable to write detection efficiency %d", i);
    75             return false;
    76         }
    77         psFree(hdrName);
    78         psFree(psfName);
    79         psFree(deteffName);
    80     }
     147    psFree(header);
    81148    psFitsClose(fits);
    82149
    83     psTrace("ppMops.write", 1, "Done writing to %s", args->output);
     150    psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
    84151
    85152    return true;
Note: See TracChangeset for help on using the changeset viewer.