IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28209


Ignore:
Timestamp:
Jun 3, 2010, 5:38:55 PM (16 years ago)
Author:
Paul Price
Message:

In the process of reworking ppMops to conform to new interface format, which is SMF-like.

Location:
trunk/ppTranslate/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMops.h

    r28027 r28209  
    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
    129
    1310// Configuration data
    … …  
    2724} ppMopsArguments;
    2825
    29 #if 0
    30 TTYPE19 = 'PSF_CHISQ'          / label for field  19
    31 TFORM19 = '1E      '           / data format of field: 4-byte REAL
    32 TTYPE20 = 'CR_NSIGMA'          / label for field  20
    33 TFORM20 = '1E      '           / data format of field: 4-byte REAL
    34 TTYPE21 = 'EXT_NSIGMA'         / label for field  21
    35 TFORM21 = '1E      '           / data format of field: 4-byte REAL
    36 TTYPE22 = 'PSF_MAJOR'          / label for field  22
    37 TFORM22 = '1E      '           / data format of field: 4-byte REAL
    38 TTYPE23 = 'PSF_MINOR'          / label for field  23
    39 TFORM23 = '1E      '           / data format of field: 4-byte REAL
    40 TTYPE24 = 'PSF_THETA'          / label for field  24
    41 TFORM24 = '1E      '           / data format of field: 4-byte REAL
    42 TTYPE25 = 'PSF_QF  '           / label for field  25
    43 TFORM25 = '1E      '           / data format of field: 4-byte REAL
    44 TTYPE26 = 'PSF_NDOF'           / label for field  26
    45 TFORM26 = '1J      '           / data format of field: 4-byte INTEGER
    46 TTYPE27 = 'PSF_NPIX'           / label for field  27
    47 TFORM27 = '1J      '           / data format of field: 4-byte INTEGER
    48 TTYPE28 = 'MOMENTS_XX'         / label for field  28
    49 TFORM28 = '1E      '           / data format of field: 4-byte REAL
    50 TTYPE29 = 'MOMENTS_XY'         / label for field  29
    51 TFORM29 = '1E      '           / data format of field: 4-byte REAL
    52 TTYPE30 = 'MOMENTS_YY'         / label for field  30
    53 TFORM30 = '1E      '           / data format of field: 4-byte REAL
    54 #endif
    55 
    5626/// Parse arguments
    5727ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]);
    5828
    5929typedef struct {
     30    psMetadata *header;                 // FITS header
    6031    psString raBoresight, decBoresight; // RA,Dec of telescope boresight
    6132    psString filter;                    // Filter for exposure
    … …  
    6435    double posangle;                    // Position angle
    6536    double alt, az;                     // Telescope altitude and azimuth
    66     double mjd;                         // Modified Julian Date
     37    double mjd;                         // Modified Julian Date of exposure mid-point
    6738    float seeing;                       // Seeing of exposure
     39    int naxis1, naxis2;                 // Size of image
    6840    long num;                           // Number of detections
     41    psMetadata *table;                  // Columns of data
    6942    psVector *x, *y;                    // Image coordinates
    7043    psVector *ra, *dec;                 // Sky coordinates
    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
    8744} ppMopsDetections;
    8845
    89 
    90 ppMopsDetections *ppMopsDetectionsAlloc(long num);
    91 
    92 /// Copy a detection
    93 bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index);
    94 
    95 /// Purge the detections list of masked detections
    96 bool ppMopsDetectionsPurge(ppMopsDetections *detections);
    97 
     46// Allocator
     47ppMopsDetections *ppMopsDetectionsAlloc(void);
    9848
    9949/// Read detections
    10050psArray *ppMopsRead(const ppMopsArguments *args);
    10151
    102 /// Merge detections
    103 ppMopsDetections *ppMopsMerge(const psArray *detections);
     52/// Purge duplicate detections
     53ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections);
    10454
    10555/// Write detections
    106 bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args);
     56bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args);
    10757
    10858#endif

  • trunk/ppTranslate/src/ppMopsDetections.c

    r28027 r28209  
    1313    psFree(det->decBoresight);
    1414    psFree(det->filter);
     15    psFree(det->header);
     16    psFree(det->table);
    1517    psFree(det->x);
    1618    psFree(det->y);
    1719    psFree(det->ra);
    1820    psFree(det->dec);
    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);
    4521
    4622    return;
    4723}
    4824
    49 ppMopsDetections *ppMopsDetectionsAlloc(long num)
     25ppMopsDetections *ppMopsDetectionsAlloc(void)
    5026{
    5127    ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return
    … …  
    6238    det->mjd = NAN;
    6339    det->seeing = NAN;
    64     det->num = 0;
    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);
     40    det->header = NULL;
     41    det->table = NULL;
     42    det->x = NULL;
     43    det->y = NULL;
     44    det->ra = NULL;
     45    det->dec = NULL;
    9546
    9647    return det;
    9748}
    98 
    99 
    100 ppMopsDetections *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 
    137 bool 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 
    180 bool 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 
    219 bool 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

    r25256 r28209  
    2828
    2929
    30 ppMopsDetections *ppMopsMerge(const psArray *detections)
     30ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections)
    3131{
    3232    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    3333
    34     psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
     34    int numInputs = detections->n;                // Number of inputs
     35    psTrace("ppMops.merge", 1, "Checking detections from %ld inputs\n", numInputs);
     36
     37    psArray *dupes = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
     38    for (int i = 0; i < numInputs; i++) {
     39        ppMopsDetections *det = detections->data[i]; // Detections from
     40
     41
     42
     43
     44
     45        dupes->data[i] = psVector
    3546
    3647    ppMopsDetections *merged = NULL;    // Merged list

  • trunk/ppTranslate/src/ppMopsRead.c

    r28027 r28209  
    2727        }
    2828
    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);
    32             return false;
    33         }
    34 
    3529        if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
    3630            psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
    … …  
    4842        }
    4943        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
    50 
    51         psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
     44        det->header = header;
    5245
    5346        det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
    … …  
    6457                             psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    6558
    66         int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    67         int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     59        det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1");
     60        det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2");
    6861
    69         psFree(header);
    70 
    71         psArray *table = psFitsReadTable(fits); // Table of interest
     62        det->table = psFitsReadTableAllColumns(fits); // Table of interest
    7263        if (!table) {
    73             psError(PS_ERR_IO, false, "Unable to read table %d", i);
     64            psError(psErrorCodeLast(), false, "Unable to read table %d", i);
    7465            return false;
    7566        }
    7667        psFitsClose(fits);
    7768
    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
     69        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", numGood, (const char*)inNames->data[i]);
    8270
    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));
     71        det->ra = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "RA_PSF"));
     72        det->dec = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "DEC_PSF"));
     73        det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
     74        det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
     75        if (!det->ra || !det->dec || !det->x || !det->y) {
     76            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
     77            return false;
    18878        }
    18979
    190         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    191 
    192         psFree(table);
    19380        detections->data[i] = det;
    19481    }

  • trunk/ppTranslate/src/ppMopsWrite.c

    r28027 r28209  
    99#include "ppTranslateVersion.h"
    1010
    11 bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args)
     11bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
    1212{
    13     psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output);
     13    psTrace("ppMops.write", 1, "Writing %ld extensions to %s", detections->n, 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);
    2137
    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);
    29 
    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);
     38        if (!psFitsWriteBlank(fits, header)) {
     39            psError(psErrorCodeLast(), false, "Unable to write primary header");
    8740            return false;
    8841        }
    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]);
     42        psFree(header);
     43    }
    13544
    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);
     45    for (int i = 0; i < detections->n; i++) {
     46        ppMopsDetections *det = detections->data[i]; // Detections for extension
     47        psTrace("ppMops.write", 1, "Writing extension %d to %s", i, args->output);
     48        if (!psFitsWriteTableAllColumns(fits, det->header, det->table, NULL)) {
     49            psError(psErrorCodeLast(), false, "Unable to write extension %d", i);
    14250            return false;
    14351        }
    144         psFree(table);
    14552    }
    146 
    147     psFree(header);
    14853    psFitsClose(fits);
    14954
    150     psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
     55    psTrace("ppMops.write", 1, "Done writing to %s", args->output);
    15156
    15257    return true;
Note: See TracChangeset for help on using the changeset viewer.