Index: /tags/ipp-20100623/ppTranslate/src/ppMops.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMops.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMops.c (revision 28431) @@ -20,15 +20,17 @@ } - if (!ppMopsPurgeDuplicates(detections)) { + ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections + psFree(detections); + if (!merged) { psErrorStackPrint(stderr, "Unable to merge detections"); exit(PS_EXIT_SYS_ERROR); } - if (!ppMopsWrite(detections, args)) { + if (!ppMopsWrite(merged, args)) { psErrorStackPrint(stderr, "Unable to write detections"); exit(PS_EXIT_SYS_ERROR); } - psFree(detections); + psFree(merged); psFree(args); Index: /tags/ipp-20100623/ppTranslate/src/ppMops.h =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMops.h (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMops.h (revision 28431) @@ -7,4 +7,7 @@ #define IN_EXTNAME "SkyChip.psf" // Extension name for data in input #define OBSERVATORY_CODE "F51" // IAU Observatory Code +#define OUT_EXTNAME "MOPS_TRANSIENT_DETECTIONS" // Extension name for data in output +#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_SATURATED | \ + PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SKY_FAILURE) // Flags to exclude // Configuration data @@ -24,10 +27,35 @@ } ppMopsArguments; +#if 0 +TTYPE19 = 'PSF_CHISQ' / label for field 19 +TFORM19 = '1E ' / data format of field: 4-byte REAL +TTYPE20 = 'CR_NSIGMA' / label for field 20 +TFORM20 = '1E ' / data format of field: 4-byte REAL +TTYPE21 = 'EXT_NSIGMA' / label for field 21 +TFORM21 = '1E ' / data format of field: 4-byte REAL +TTYPE22 = 'PSF_MAJOR' / label for field 22 +TFORM22 = '1E ' / data format of field: 4-byte REAL +TTYPE23 = 'PSF_MINOR' / label for field 23 +TFORM23 = '1E ' / data format of field: 4-byte REAL +TTYPE24 = 'PSF_THETA' / label for field 24 +TFORM24 = '1E ' / data format of field: 4-byte REAL +TTYPE25 = 'PSF_QF ' / label for field 25 +TFORM25 = '1E ' / data format of field: 4-byte REAL +TTYPE26 = 'PSF_NDOF' / label for field 26 +TFORM26 = '1J ' / data format of field: 4-byte INTEGER +TTYPE27 = 'PSF_NPIX' / label for field 27 +TFORM27 = '1J ' / data format of field: 4-byte INTEGER +TTYPE28 = 'MOMENTS_XX' / label for field 28 +TFORM28 = '1E ' / data format of field: 4-byte REAL +TTYPE29 = 'MOMENTS_XY' / label for field 29 +TFORM29 = '1E ' / data format of field: 4-byte REAL +TTYPE30 = 'MOMENTS_YY' / label for field 30 +TFORM30 = '1E ' / data format of field: 4-byte REAL +#endif + /// Parse arguments ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]); typedef struct { - psString component; // Component name - psMetadata *header; // FITS header (extension names *.hdr) psString raBoresight, decBoresight; // RA,Dec of telescope boresight psString filter; // Filter for exposure @@ -36,27 +64,45 @@ double posangle; // Position angle double alt, az; // Telescope altitude and azimuth - double mjd; // Modified Julian Date of exposure mid-point + double mjd; // Modified Julian Date float seeing; // Seeing of exposure - int naxis1, naxis2; // Size of image long num; // Number of detections - psMetadata *psfHeader; // FITS header (extension names *.psf) - psMetadata *table; // Columns of data psVector *x, *y; // Image coordinates psVector *ra, *dec; // Sky coordinates - psMetadata *deteffHeader; // Detection efficiency header (extension names *.deteff) - psMetadata *deteffTable; // Detection efficiency table + psVector *raErr, *decErr; // Error in sky coordinates + psVector *mag, *magErr; // Magnitude and associated error + psVector *chi2, *dof; // Chi^2 from fitting, with associated degrees of freedom + psVector *cr, *extended; // Measures of CR-ness and extendedness + psVector *psfMajor, *psfMinor, *psfTheta; // PSF major and minor axes, and position angle + psVector *quality, *numPix; // PSF quality factor and number of pixels + psVector *xxMoment, *xyMoment, *yyMoment; // Moments + psVector *flags; // psphot flags + psVector *diffSkyfileId; // Identifier for source image + psVector *naxis1, *naxis2; // Size of image + psVector *mask; // Mask for detections + psVector *nPos; // Number of positive pixels + psVector *fPos; // Fraction of positive flux + psVector *nRatioBad; // Fraction of positive pixels to negative + psVector *nRatioMask; // Fraction of positive pixels to masked + psVector *nRatioAll; // Fraction of positive pixels to all } ppMopsDetections; -// Allocator -ppMopsDetections *ppMopsDetectionsAlloc(void); + +ppMopsDetections *ppMopsDetectionsAlloc(long num); + +/// Copy a detection +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index); + +/// Purge the detections list of masked detections +bool ppMopsDetectionsPurge(ppMopsDetections *detections); + /// Read detections psArray *ppMopsRead(const ppMopsArguments *args); -/// Purge duplicate detections -bool ppMopsPurgeDuplicates(const psArray *detections); +/// Merge detections +ppMopsDetections *ppMopsMerge(const psArray *detections); /// Write detections -bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args); +bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); #endif Index: /tags/ipp-20100623/ppTranslate/src/ppMopsDetections.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMopsDetections.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMopsDetections.c (revision 28431) @@ -10,26 +10,46 @@ static void mopsDetectionsFree(ppMopsDetections *det) { - psFree(det->component); psFree(det->raBoresight); psFree(det->decBoresight); psFree(det->filter); - psFree(det->header); - psFree(det->table); psFree(det->x); psFree(det->y); psFree(det->ra); psFree(det->dec); - psFree(det->deteffHeader); - psFree(det->deteffTable); + psFree(det->raErr); + psFree(det->decErr); + psFree(det->mag); + psFree(det->magErr); + psFree(det->chi2); + psFree(det->dof); + psFree(det->cr); + psFree(det->extended); + psFree(det->psfMajor); + psFree(det->psfMinor); + psFree(det->psfTheta); + psFree(det->quality); + psFree(det->numPix); + psFree(det->xxMoment); + psFree(det->xyMoment); + psFree(det->yyMoment); + psFree(det->flags); + psFree(det->diffSkyfileId); + psFree(det->naxis1); + psFree(det->naxis2); + psFree(det->mask); + psFree(det->nPos); + psFree(det->fPos); + psFree(det->nRatioBad); + psFree(det->nRatioMask); + psFree(det->nRatioAll); return; } -ppMopsDetections *ppMopsDetectionsAlloc(void) +ppMopsDetections *ppMopsDetectionsAlloc(long num) { ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree); - det->component = NULL; det->raBoresight = NULL; det->decBoresight = NULL; @@ -43,13 +63,234 @@ det->seeing = NAN; det->num = 0; - det->header = NULL; - det->table = NULL; - det->x = NULL; - det->y = NULL; - det->ra = NULL; - det->dec = NULL; - det->deteffHeader = NULL; - det->deteffTable = NULL; + det->x = psVectorAllocEmpty(num, PS_TYPE_F32); + det->y = psVectorAllocEmpty(num, PS_TYPE_F32); + det->ra = psVectorAllocEmpty(num, PS_TYPE_F64); + det->dec = psVectorAllocEmpty(num, PS_TYPE_F64); + det->raErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->decErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->mag = psVectorAllocEmpty(num, PS_TYPE_F32); + det->magErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->chi2 = psVectorAllocEmpty(num, PS_TYPE_F32); + det->dof = psVectorAllocEmpty(num, PS_TYPE_S32); + det->cr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->extended = psVectorAllocEmpty(num, PS_TYPE_F32); + det->psfMajor = psVectorAllocEmpty(num, PS_TYPE_F32); + det->psfMinor = psVectorAllocEmpty(num, PS_TYPE_F32); + det->psfTheta = psVectorAllocEmpty(num, PS_TYPE_F32); + det->quality = psVectorAllocEmpty(num, PS_TYPE_F32); + det->numPix = psVectorAllocEmpty(num, PS_TYPE_S32); + det->xxMoment = psVectorAllocEmpty(num, PS_TYPE_F32); + det->xyMoment = psVectorAllocEmpty(num, PS_TYPE_F32); + det->yyMoment = psVectorAllocEmpty(num, PS_TYPE_F32); + det->flags = psVectorAllocEmpty(num, PS_TYPE_U32); + det->diffSkyfileId = psVectorAllocEmpty(num, PS_TYPE_S64); + det->naxis1 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->naxis2 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->mask = psVectorAllocEmpty(num, PS_TYPE_U8); + det->nPos = psVectorAllocEmpty(num, PS_TYPE_S32); + det->fPos = psVectorAllocEmpty(num, PS_TYPE_F32); + det->nRatioBad = psVectorAllocEmpty(num, PS_TYPE_F32); + det->nRatioMask = psVectorAllocEmpty(num, PS_TYPE_F32); + det->nRatioAll = psVectorAllocEmpty(num, PS_TYPE_F32); return det; } + + +ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num) +{ + det->x = psVectorRealloc(det->x, num); + det->y = psVectorRealloc(det->y, num); + det->ra = psVectorRealloc(det->ra, num); + det->dec = psVectorRealloc(det->dec, num); + det->raErr = psVectorRealloc(det->raErr, num); + det->decErr = psVectorRealloc(det->decErr, num); + det->mag = psVectorRealloc(det->mag, num); + det->magErr = psVectorRealloc(det->magErr, num); + det->chi2 = psVectorRealloc(det->chi2, num); + det->dof = psVectorRealloc(det->dof, num); + det->cr = psVectorRealloc(det->cr, num); + det->extended = psVectorRealloc(det->extended, num); + det->psfMajor = psVectorRealloc(det->psfMajor, num); + det->psfMinor = psVectorRealloc(det->psfMinor, num); + det->psfTheta = psVectorRealloc(det->psfTheta, num); + det->quality = psVectorRealloc(det->quality, num); + det->numPix = psVectorRealloc(det->numPix, num); + det->xxMoment = psVectorRealloc(det->xxMoment, num); + det->xyMoment = psVectorRealloc(det->xyMoment, num); + det->yyMoment = psVectorRealloc(det->yyMoment, num); + det->flags = psVectorRealloc(det->flags, num); + det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num); + det->naxis1 = psVectorRealloc(det->naxis1, num); + det->naxis2 = psVectorRealloc(det->naxis2, num); + det->mask = psVectorRealloc(det->mask, num); + det->nPos = psVectorRealloc(det->nPos, num); + det->fPos = psVectorRealloc(det->fPos, num); + det->nRatioBad = psVectorRealloc(det->nRatioBad, num); + det->nRatioMask = psVectorRealloc(det->nRatioMask, num); + det->nRatioAll = psVectorRealloc(det->nRatioAll, num); + + return det; +} + + +bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec, + double raErr, double decErr, float mag, float magErr, + float chi2, int dof, float cr, float extended, float psfMajor, + float psfMinor, float psfTheta, float quality, int numPix, + float xxMoment, float xyMoment, float yyMoment, + psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2, + int nPos, float fPos, float nRatioBad, float nRatioMask, float nRatioAll) +{ + psVectorAppend(det->x, x); + psVectorAppend(det->y, y); + psVectorAppend(det->ra, ra); + psVectorAppend(det->dec, dec); + psVectorAppend(det->raErr, raErr); + psVectorAppend(det->decErr, decErr); + psVectorAppend(det->mag, mag); + psVectorAppend(det->magErr, magErr); + psVectorAppend(det->chi2, chi2); + psVectorAppend(det->dof, dof); + psVectorAppend(det->cr, cr); + psVectorAppend(det->extended, extended); + psVectorAppend(det->psfMajor, psfMajor); + psVectorAppend(det->psfMinor, psfMinor); + psVectorAppend(det->psfTheta, psfTheta); + psVectorAppend(det->quality, quality); + psVectorAppend(det->numPix, numPix); + psVectorAppend(det->xxMoment, xxMoment); + psVectorAppend(det->xyMoment, xyMoment); + psVectorAppend(det->yyMoment, yyMoment); + psVectorAppend(det->flags, flags); + psVectorAppend(det->diffSkyfileId, diffSkyfileId); + psVectorAppend(det->naxis1, naxis1); + psVectorAppend(det->naxis2, naxis2); + psVectorAppend(det->mask, 0); + psVectorAppend(det->nPos, nPos); + psVectorAppend(det->fPos, fPos); + psVectorAppend(det->nRatioBad, nRatioBad); + psVectorAppend(det->nRatioMask, nRatioMask); + psVectorAppend(det->nRatioAll, nRatioAll); + + return true; +} + + +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index) +{ + psVectorAppend(target->x, source->x->data.F32[index]); + psVectorAppend(target->y, source->y->data.F32[index]); + psVectorAppend(target->ra, source->ra->data.F64[index]); + psVectorAppend(target->dec, source->dec->data.F64[index]); + psVectorAppend(target->raErr, source->raErr->data.F64[index]); + psVectorAppend(target->decErr, source->decErr->data.F64[index]); + psVectorAppend(target->mag, source->mag->data.F32[index]); + psVectorAppend(target->magErr, source->magErr->data.F32[index]); + psVectorAppend(target->chi2, source->chi2->data.F32[index]); + psVectorAppend(target->dof, source->dof->data.S32[index]); + psVectorAppend(target->cr, source->cr->data.F32[index]); + psVectorAppend(target->extended, source->extended->data.F32[index]); + psVectorAppend(target->psfMajor, source->psfMajor->data.F32[index]); + psVectorAppend(target->psfMinor, source->psfMinor->data.F32[index]); + psVectorAppend(target->psfTheta, source->psfTheta->data.F32[index]); + psVectorAppend(target->quality, source->quality->data.F32[index]); + psVectorAppend(target->numPix, source->numPix->data.S32[index]); + psVectorAppend(target->xxMoment, source->xxMoment->data.F32[index]); + psVectorAppend(target->xyMoment, source->xyMoment->data.F32[index]); + psVectorAppend(target->yyMoment, source->yyMoment->data.F32[index]); + psVectorAppend(target->flags, source->flags->data.U32[index]); + psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]); + psVectorAppend(target->naxis1, source->naxis1->data.S32[index]); + psVectorAppend(target->naxis2, source->naxis2->data.S32[index]); + psVectorAppend(target->mask, 0); + psVectorAppend(target->nPos, source->nPos->data.S32[index]); + psVectorAppend(target->fPos, source->fPos->data.F32[index]); + psVectorAppend(target->nRatioBad, source->nRatioBad->data.F32[index]); + psVectorAppend(target->nRatioMask, source->nRatioMask->data.F32[index]); + psVectorAppend(target->nRatioAll, source->nRatioAll->data.F32[index]); + + target->num++; + + return true; +} + + +bool ppMopsDetectionsPurge(ppMopsDetections *det) +{ + long num = 0; + for (long i = 0; i < det->num; i++) { + if (!det->mask->data.U8[i]) { + if (i == num) { + // No need to copy + num++; + continue; + } + det->x->data.F32[num] = det->x->data.F32[i]; + det->y->data.F32[num] = det->y->data.F32[i]; + det->ra->data.F64[num] = det->ra->data.F64[i]; + det->dec->data.F64[num] = det->dec->data.F64[i]; + det->raErr->data.F64[num] = det->raErr->data.F64[i]; + det->decErr->data.F64[num] = det->decErr->data.F64[i]; + det->mag->data.F32[num] = det->mag->data.F32[i]; + det->magErr->data.F32[num] = det->magErr->data.F32[i]; + det->chi2->data.F32[num] = det->chi2->data.F32[i]; + det->dof->data.S32[num] = det->dof->data.S32[i]; + det->cr->data.F32[num] = det->cr->data.F32[i]; + det->extended->data.F32[num] = det->extended->data.F32[i]; + det->psfMajor->data.F32[num] = det->psfMajor->data.F32[i]; + det->psfMinor->data.F32[num] = det->psfMinor->data.F32[i]; + det->psfTheta->data.F32[num] = det->psfTheta->data.F32[i]; + det->quality->data.F32[num] = det->quality->data.F32[i]; + det->numPix->data.S32[num] = det->numPix->data.S32[i]; + det->xxMoment->data.F32[num] = det->xxMoment->data.F32[i]; + det->xyMoment->data.F32[num] = det->xyMoment->data.F32[i]; + det->yyMoment->data.F32[num] = det->yyMoment->data.F32[i]; + det->flags->data.U32[num] = det->flags->data.U32[i]; + det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i]; + det->naxis1->data.S32[num] = det->naxis1->data.S32[i]; + det->naxis2->data.S32[num] = det->naxis2->data.S32[i]; + det->mask->data.U8[num] = 0; + det->nPos->data.S32[num] = det->nPos->data.S32[i]; + det->fPos->data.F32[num] = det->fPos->data.F32[i]; + det->nRatioBad->data.F32[num] = det->nRatioBad->data.F32[i]; + det->nRatioMask->data.F32[num] = det->nRatioMask->data.F32[i]; + det->nRatioAll->data.F32[num] = det->nRatioAll->data.F32[i]; + num++; + } + } + det->x->n = num; + det->y->n = num; + det->ra->n = num; + det->dec->n = num; + det->raErr->n = num; + det->decErr->n = num; + det->mag->n = num; + det->magErr->n = num; + det->chi2->n = num; + det->dof->n = num; + det->cr->n = num; + det->extended->n = num; + det->psfMajor->n = num; + det->psfMinor->n = num; + det->psfTheta->n = num; + det->quality->n = num; + det->numPix->n = num; + det->xxMoment->n = num; + det->xyMoment->n = num; + det->yyMoment->n = num; + det->flags->n = num; + det->diffSkyfileId->n = num; + det->naxis1->n = num; + det->naxis2->n = num; + det->mask->n = num; + det->num = num; + det->nPos->n = num; + det->fPos->n = num; + det->nRatioBad->n = num; + det->nRatioMask->n = num; + det->nRatioAll->n = num; + + return true; +} + Index: /tags/ipp-20100623/ppTranslate/src/ppMopsMerge.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMopsMerge.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMopsMerge.c (revision 28431) @@ -9,5 +9,5 @@ #include "ppMops.h" -#define LEAF_SIZE 2 // Size of leaf +#define LEAF_SIZE 4 // Size of leaf #define MATCH_RADIUS SEC_TO_RAD(1.0) // Matching radius #define MJD_TOL 1.0/3600.0/24.0 // Tolerance for MJD matching @@ -17,12 +17,4 @@ #define AIRMASS_TOL 1.0e-3 // Tolerance for airmass matching -#if 0 -#undef psTrace -#define psTrace(facil, level, ...) \ - if (level <= 5) { \ - fprintf(stderr, __VA_ARGS__); \ - } -#endif - // Get distance from detection to centre of image static float mergeDistance(const ppMopsDetections *detections, // Detections of interest @@ -30,45 +22,19 @@ ) { - float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0; - float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0; + float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0; return PS_SQR(dx) + PS_SQR(dy); } -bool ppMopsPurgeDuplicates(const psArray *detections) +ppMopsDetections *ppMopsMerge(const psArray *detections) { PS_ASSERT_ARRAY_NON_NULL(detections, NULL); - int numInputs = detections->n; // Number of inputs - psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs); + psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); - long total = 0; // Total number of sources - psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input - psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input - psVectorInit(dupNum, 0); - for (int i = 0; i < numInputs; i++) { - ppMopsDetections *det = detections->data[i]; // Detections from - if (!det || det->num == 0) { - continue; - } - psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8); - psVectorInit(dupes, 0); - total += det->num; - } - psTrace("ppMops.merge", 2, "Total detections: %ld\n", total); - - psVector *raMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged RAs - psVector *decMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged Decs - psVector *sourceMerged = psVectorAlloc(total, PS_TYPE_U16); // Source image of merged sources - psVector *indexMerged = psVectorAlloc(total, PS_TYPE_U32); // Index of merged sources - long num = 0; // Number of merged sources - - const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates - const char *filter = NULL; // Filter name - float airmass = NAN, exptime = NAN; // Exposure details - double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing - double mjd = NAN; // Time of exposure - - for (int i = 0; i < numInputs; i++) { + ppMopsDetections *merged = NULL; // Merged list + int num = 1; // Number of merged files + for (int i = 0; i < detections->n; i++) { ppMopsDetections *det = detections->data[i]; // Detections of interest if (!det) { @@ -79,90 +45,69 @@ continue; } + num++; + if (!merged) { + psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); + merged = psMemIncrRefCounter(det); + continue; + } + psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i); - // Check exposure characteristics - if (num == 0) { - raBoresight = det->raBoresight; - decBoresight = det->decBoresight; - filter = det->filter; - airmass = det->airmass; - exptime = det->exptime; - posangle = det->posangle; - alt = det->alt; - az = det->az; - } else { - if (strcmp(raBoresight, det->raBoresight) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s", - raBoresight, det->raBoresight); - return false; - } - if (strcmp(decBoresight, det->decBoresight) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s", - decBoresight, det->decBoresight); - return false; - } - if (strcmp(filter, det->filter) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s", - filter, det->filter); - return false; - } - if (fabsf(airmass - det->airmass) > AIRMASS_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", - airmass, det->airmass); - return false; - } - if (fabsf(exptime - det->exptime) > EXPTIME_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", - exptime, det->exptime); - return false; - } - if (fabs(posangle - det->posangle) > POSANGLE_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", - posangle, det->posangle); - return false; - } - if (fabs(alt - det->alt) > BORESIGHT_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", - alt, det->alt); - return false; - } - if (fabs(az - det->az) > BORESIGHT_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", - az, det->az); - return false; - } - if (fabs(mjd - det->mjd) > MJD_TOL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf", - mjd, det->mjd); - return false; - } + // XXX compare exposure properties + if (strcmp(merged->raBoresight, det->raBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s", + merged->raBoresight, det->raBoresight); + return NULL; + } + if (strcmp(merged->decBoresight, det->decBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s", + merged->decBoresight, det->decBoresight); + return NULL; + } + if (strcmp(merged->filter, det->filter) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s", + merged->filter, det->filter); + return NULL; } - psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); - memcpy(&raMerged->data.F64[num], det->ra->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - memcpy(&decMerged->data.F64[num], det->dec->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - for (long j = 0, k = num; j < det->num; j++, k++) { - sourceMerged->data.U16[k] = i; - indexMerged->data.U32[k] = j; + if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", + merged->airmass, det->airmass); + return NULL; } - num += det->num; - } + if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", + merged->exptime, det->exptime); + return NULL; + } + if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", + merged->posangle, det->posangle); + return NULL; + } + if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", + merged->alt, det->alt); + return NULL; + } + if (fabs(merged->az - det->az) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", + merged->az, det->az); + return NULL; + } + if (fabs(merged->mjd - det->mjd) > MJD_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf", + merged->mjd, det->mjd); + return NULL; + } - psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num); - psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree - if (!tree) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree"); - return false; - } + merged->seeing += det->seeing; // Taking average - for (int i = 0; i < detections->n; i++) { - ppMopsDetections *det = detections->data[i]; // Detections of interest - if (!det) { - continue; - } else if (det->num == 0) { - continue; + psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree + if (!tree) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree"); + psFree(merged); + return NULL; } - psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i); - psVector *dupes = duplicates->data[i]; // Duplicates list psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest for (int j = 0; j < det->num; j++) { @@ -174,14 +119,14 @@ psFree(coords); psFree(tree); - return false; + psFree(merged); + return NULL; + } + if (indices->n == 0) { + psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i); + psFree(indices); + ppMopsDetectionsCopySingle(merged, det, j); + continue; } psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); - psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i); - - if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) { - // It's myself - psFree(indices); - continue; - } // Which one do we keep? @@ -189,16 +134,6 @@ long bestIndex = -1; // Index with best distance for (int k = 0; k < indices->n; k++) { - long mergeIndex = indices->data.S64[k]; // Index of point in merged list - int source = sourceMerged->data.U16[mergeIndex]; // Source image - if (source == i) { - continue; - } - long index = indexMerged->data.U32[mergeIndex]; // Index in source - psVector *dupes = duplicates->data[source]; // Duplicates list - if (dupes->data.U8[index]) { - continue; - } - - float distance = mergeDistance(detections->data[source], index); // Distance to centre of image + long index = indices->data.S64[k]; // Index of point + float distance = mergeDistance(merged, index); // Distance to centre of image if (distance < bestDistance) { bestDistance = distance; @@ -207,84 +142,30 @@ } - if (bestIndex >= 0) { - float distance = mergeDistance(det, j); // Distance to centre of image - if (bestIndex >= 0 && distance < bestDistance) { - psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); - for (int k = 0; k < indices->n; k++) { - long mergeIndex = indices->data.S64[k]; // Index of point - int source = sourceMerged->data.U16[mergeIndex]; // Source image - if (source == i) { - continue; - } - long index = indexMerged->data.U32[mergeIndex]; // Index in source - psVector *dupes = duplicates->data[source]; // Duplicates list - if (!dupes->data.U8[index]) { - dupes->data.U8[index] = 0xFF; - dupNum->data.U32[source]++; - } - } - } else { - psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); - dupes->data.U8[j] = 0xFF; - dupNum->data.U32[i]++; + float distance = mergeDistance(det, j); // Distance to centre of image + if (distance < bestDistance) { + psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); + // Blow away existing sources + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + merged->mask->data.U8[index] = 0xFF; } + ppMopsDetectionsCopySingle(merged, det, j); + } else { + psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); } psFree(indices); } - psFree(coords); - } - psFree(tree); - psFree(raMerged); - psFree(decMerged); - psFree(sourceMerged); - psFree(indexMerged); - // Remove duplicates - for (int i = 0; i < detections->n; i++) { - ppMopsDetections *det = detections->data[i]; // Detections of interest - if (!det) { - continue; - } else if (det->num == 0) { - continue; - } - psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i); + psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num); -#define VECTOR_PURGE_CASE(TYPE) \ - case PS_TYPE_##TYPE: \ - for (int i = 0, j = 0; i < vector->n; i++) { \ - if (!dupes->data.U8[i]) { \ - if (i == j) { \ - j++; \ - continue; \ - } \ - vector->data.TYPE[j] = vector->data.TYPE[i]; \ - } \ - } \ - break; - - psVector *dupes = duplicates->data[i]; // Duplicates - psArray *table = psListToArray(det->table->list); // Table of data - for (int j = 0; j < table->n; j++) { - psMetadataItem *item = table->data[j]; // Table item - psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type); - psVector *vector = item->data.V; // Vector to purge - switch (vector->type.type) { - VECTOR_PURGE_CASE(U8); - VECTOR_PURGE_CASE(U16); - VECTOR_PURGE_CASE(U32); - VECTOR_PURGE_CASE(U64); - VECTOR_PURGE_CASE(S8); - VECTOR_PURGE_CASE(S16); - VECTOR_PURGE_CASE(S32); - VECTOR_PURGE_CASE(S64); - VECTOR_PURGE_CASE(F32); - VECTOR_PURGE_CASE(F64); - default: - psAbort("Unrecognised vector type: %x", vector->type.type); - } - } + psFree(tree); + ppMopsDetectionsPurge(merged); } - return true; + psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); + + merged->seeing /= num; + + return merged; } Index: /tags/ipp-20100623/ppTranslate/src/ppMopsRead.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMopsRead.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMopsRead.c (revision 28431) @@ -4,5 +4,4 @@ #include -#include #include @@ -17,6 +16,5 @@ psArray *detections = psArrayAlloc(num); // Array of detections, to return for (int i = 0; i < num; i++) { - const char *name = inNames->data[i]; // File name - psFits *fits = psFitsOpen(name, "r"); // FITS file + psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file if (!fits) { psError(PS_ERR_IO, false, "Unable to open input %d", i); @@ -26,4 +24,10 @@ if (!header) { psError(PS_ERR_IO, false, "Unable to read header %d", i); + return false; + } + + psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image + if (diffSkyfileId == 0) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i); return false; } @@ -43,8 +47,7 @@ continue; } - ppMopsDetections *det = detections->data[i] = ppMopsDetectionsAlloc(); - det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension - det->header = header; - det->num = size; + ppMopsDetections *det = ppMopsDetectionsAlloc(size); + + psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); @@ -61,55 +64,132 @@ psMetadataLookupF32(NULL, header, "FWHM_MIN")); - det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); - det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); + int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows - det->psfHeader = psFitsReadHeader(NULL, fits); - if (!det->psfHeader) { - psError(psErrorCodeLast(), false, "Unable to read header %d", i); + psFree(header); + + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); return false; } - psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest - if (!table) { - psError(psErrorCodeLast(), false, "Unable to read table %d", i); - return false; + psFitsClose(fits); + + double plateScale = 0.0; // Plate scale + long numGood = 0; // Number of good rows + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + + psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); + if (flags & SOURCE_MASK) { + psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags); + continue; + } + + det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF"); + det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF")); + det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF")); + det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); + det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); + det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ"); + det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF"); + det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA"); + det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); + det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR"); + det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR"); + det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA"); + det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF"); + det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX"); + det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX"); + det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY"); + det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY"); + det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS"); + det->diffSkyfileId->data.S64[numGood] = diffSkyfileId; + det->naxis1->data.S32[numGood] = naxis1; + det->naxis2->data.S32[numGood] = naxis2; + + det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS"); + det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO"); + det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD"); + det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK"); + det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL"); + + // Calculate error in RA, Dec + double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); + double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); + double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); + double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); + + if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || + !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || + !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || + !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { + psTrace("ppMops.read", 10, + "Discarding row %ld from input %d because of non-finite values: " + "%f %f %lf %lf %f %f %f %f %f %f", + j, i, + det->x->data.F32[numGood], det->y->data.F32[numGood], + det->ra->data.F64[numGood], det->dec->data.F64[numGood], + det->mag->data.F32[numGood], det->magErr->data.F32[numGood], + xErr, yErr, scale, angle); + continue; + } + + // XXX Not at all sure I've got the angles around the right way here... + double cosAngle = cos(angle), sinAngle = sin(angle); + double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); + double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); + double errScale = scale / 3600.0; + det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[numGood] = 0; + plateScale += scale; + numGood++; + } + det->seeing *= plateScale / numGood; + + det->x->n = numGood; + det->y->n = numGood; + det->ra->n = numGood; + det->dec->n = numGood; + det->raErr->n = numGood; + det->decErr->n = numGood; + det->mag->n = numGood; + det->magErr->n = numGood; + det->chi2->n = numGood; + det->dof->n = numGood; + det->cr->n = numGood; + det->extended->n = numGood; + det->psfMajor->n = numGood; + det->psfMinor->n = numGood; + det->psfTheta->n = numGood; + det->quality->n = numGood; + det->numPix->n = numGood; + det->xxMoment->n = numGood; + det->xyMoment->n = numGood; + det->yyMoment->n = numGood; + det->flags->n = numGood; + det->diffSkyfileId->n = numGood; + det->naxis1->n = numGood; + det->naxis2->n = numGood; + det->mask->n = numGood; + det->nPos->n = numGood; + det->fPos->n = numGood; + det->nRatioBad->n = numGood; + det->nRatioMask->n = numGood; + det->nRatioAll->n = numGood; + + det->num = numGood; + + if (isfinite(args->zp) && numGood > 0) { + psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); } - psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); - psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); - det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); - det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); - if (!det->ra || !det->dec || !det->x || !det->y) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); - return false; - } - - psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); - - if (!psFitsMoveExtName(fits, "SkyChip.deteff")) { - psWarning("No detection efficiencies included in %s", det->component); - psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component); - psErrorClear(); - continue; - } - - det->deteffHeader = psFitsReadHeader(NULL, fits); - if (!det->deteffHeader) { - psWarning("Unable to read detection efficiency header in %s", det->component); - psErrorClear(); - continue; - } - det->deteffTable = psFitsReadTableAllColumns(fits); - if (!det->deteffTable) { - psFree(det->deteffHeader); - det->deteffHeader = NULL; - psWarning("Unable to read detection efficiency table in %s", det->component); - psErrorClear(); - continue; - } - psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component); - psFitsClose(fits); + psFree(table); + detections->data[i] = det; } Index: /tags/ipp-20100623/ppTranslate/src/ppMopsWrite.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppMopsWrite.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppMopsWrite.c (revision 28431) @@ -9,7 +9,7 @@ #include "ppTranslateVersion.h" -bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args) +bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) { - psTrace("ppMops.write", 1, "Writing %ld extensions to %s", detections->n, args->output); + psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output); psFits *fits = psFitsOpen(args->output, "w"); // FITS file @@ -19,67 +19,134 @@ } - // Primary header - { - psMetadata *header = psMetadataAlloc(); // Header to write - ppTranslateVersionHeader(header); - psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name); - psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id); - psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id); - psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id); - psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id); - psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id); - psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id); - psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); - psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp); - psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr); - psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom); - psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); - ppMopsDetections *det = detections->data[0]; // Representative set of detections - psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight); - psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight); - psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); - psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass at boresight", det->airmass); - psMetadataAddF32(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); - psMetadataAddF64(header, PS_LIST_TAIL, "POSANG", 0, "Position angle", det->posangle); - psMetadataAddF64(header, PS_LIST_TAIL, "ALT", 0, "Altitude of boresight", det->alt); - psMetadataAddF64(header, PS_LIST_TAIL, "AZ", 0, "Azimuth of boresight", det->az); - psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of observation", det->mjd); + psMetadata *header = psMetadataAlloc(); // Header to write + psString source = ppTranslateSource(), version = ppTranslateVersion(); + psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); + psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); + ppTranslateVersionHeader(header); + psFree(source); + psFree(version); - if (!psFitsWriteBlank(fits, header, NULL)) { - psError(psErrorCodeLast(), false, "Unable to write primary header"); + psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name); + psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id); + psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id); + psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id); + psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); + + psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd); + psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight); + psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az); + psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); + psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle); + psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); + psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass); + psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr); + psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom); + + if (det->num == 0) { + // Write dummy table + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0); + + psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN); + if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write empty table."); + psFree(header); + psFree(row); return false; } - psFree(header); + psFree(row); + } else { + psArray *table = psArrayAlloc(det->num); // Table to write + for (long i = 0; i < det->num; i++) { + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", + RAD_TO_DEG(det->ra->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", + det->raErr->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", + RAD_TO_DEG(det->dec->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", + det->decErr->data.F64[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", det->chi2->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", + det->dof->data.S32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", + det->cr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", + det->extended->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", det->psfMajor->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", det->psfMinor->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", + det->psfTheta->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", + det->quality->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", + det->numPix->data.S32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", det->xxMoment->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", det->xyMoment->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", det->yyMoment->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", + det->nPos->data.S32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", + det->fPos->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", + det->nRatioBad->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", + det->nRatioMask->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", + det->nRatioAll->data.F32[i]); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", + det->diffSkyfileId->data.S64[i]); + + table->data[i] = row; + } + if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(header); + psFree(table); + return false; + } + psFree(table); } - for (int i = 0; i < detections->n; i++) { - ppMopsDetections *det = detections->data[i]; // Detections for extension - psTrace("ppMops.write", 1, "Writing extension %d to %s", i, args->output); - psString hdrName = NULL, psfName = NULL, deteffName = NULL; - psStringAppend(&hdrName, "%s.hdr", det->component); - psStringAppend(&psfName, "%s.psf", det->component); - psStringAppend(&deteffName, "%s.deteff", det->component); - - if (!psFitsWriteBlank(fits, det->header, hdrName)) { - psError(psErrorCodeLast(), false, "Unable to write header %d", i); - return false; - } - if (!psFitsWriteTableAllColumns(fits, det->psfHeader, det->table, psfName)) { - psError(psErrorCodeLast(), false, "Unable to write table %d", i); - return false; - } - if (det->deteffHeader && det->deteffTable && - !psFitsWriteTableAllColumns(fits, det->deteffHeader, det->deteffTable, deteffName)) { - psError(psErrorCodeLast(), false, "Unable to write detection efficiency %d", i); - return false; - } - psFree(hdrName); - psFree(psfName); - psFree(deteffName); - } + psFree(header); psFitsClose(fits); - psTrace("ppMops.write", 1, "Done writing to %s", args->output); + psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output); return true; Index: /tags/ipp-20100623/ppTranslate/src/ppTranslateVersion.c =================================================================== --- /tags/ipp-20100623/ppTranslate/src/ppTranslateVersion.c (revision 28430) +++ /tags/ipp-20100623/ppTranslate/src/ppTranslateVersion.c (revision 28431) @@ -67,6 +67,4 @@ psString source = ppTranslateSource(); // Software source - psMetadataAddStr(header, PS_LIST_TAIL, "TRANSL_V", PS_META_REPLACE, NULL, PPTRANSLATE_VERSION); - psStringPrepend(&version, "ppTranslate version: "); psStringPrepend(&source, "ppTranslate source: ");