Index: /trunk/ppTranslate/src/ppMops.h =================================================================== --- /trunk/ppTranslate/src/ppMops.h (revision 28208) +++ /trunk/ppTranslate/src/ppMops.h (revision 28209) @@ -7,7 +7,4 @@ #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 @@ -27,35 +24,9 @@ } 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 { + psMetadata *header; // FITS header psString raBoresight, decBoresight; // RA,Dec of telescope boresight psString filter; // Filter for exposure @@ -64,45 +35,24 @@ double posangle; // Position angle double alt, az; // Telescope altitude and azimuth - double mjd; // Modified Julian Date + double mjd; // Modified Julian Date of exposure mid-point float seeing; // Seeing of exposure + int naxis1, naxis2; // Size of image long num; // Number of detections + psMetadata *table; // Columns of data psVector *x, *y; // Image coordinates psVector *ra, *dec; // Sky coordinates - 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; - -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); - +// Allocator +ppMopsDetections *ppMopsDetectionsAlloc(void); /// Read detections psArray *ppMopsRead(const ppMopsArguments *args); -/// Merge detections -ppMopsDetections *ppMopsMerge(const psArray *detections); +/// Purge duplicate detections +ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections); /// Write detections -bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args); #endif Index: /trunk/ppTranslate/src/ppMopsDetections.c =================================================================== --- /trunk/ppTranslate/src/ppMopsDetections.c (revision 28208) +++ /trunk/ppTranslate/src/ppMopsDetections.c (revision 28209) @@ -13,39 +13,15 @@ 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->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(long num) +ppMopsDetections *ppMopsDetectionsAlloc(void) { ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return @@ -62,235 +38,11 @@ det->mjd = NAN; det->seeing = NAN; - det->num = 0; - 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); + det->header = NULL; + det->table = NULL; + det->x = NULL; + det->y = NULL; + det->ra = NULL; + det->dec = NULL; 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: /trunk/ppTranslate/src/ppMopsMerge.c =================================================================== --- /trunk/ppTranslate/src/ppMopsMerge.c (revision 28208) +++ /trunk/ppTranslate/src/ppMopsMerge.c (revision 28209) @@ -28,9 +28,20 @@ -ppMopsDetections *ppMopsMerge(const psArray *detections) +ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections) { PS_ASSERT_ARRAY_NON_NULL(detections, NULL); - psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); + int numInputs = detections->n; // Number of inputs + psTrace("ppMops.merge", 1, "Checking detections from %ld inputs\n", numInputs); + + psArray *dupes = psArrayAlloc(numInputs); // Vector of duplicate bits for each input + for (int i = 0; i < numInputs; i++) { + ppMopsDetections *det = detections->data[i]; // Detections from + + + + + + dupes->data[i] = psVector ppMopsDetections *merged = NULL; // Merged list Index: /trunk/ppTranslate/src/ppMopsRead.c =================================================================== --- /trunk/ppTranslate/src/ppMopsRead.c (revision 28208) +++ /trunk/ppTranslate/src/ppMopsRead.c (revision 28209) @@ -27,10 +27,4 @@ } - 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; - } - if (!psFitsMoveExtName(fits, "SkyChip.psf")) { psError(PS_ERR_IO, false, "Unable to move to HDU with detections"); @@ -48,6 +42,5 @@ } ppMopsDetections *det = ppMopsDetectionsAlloc(size); - - psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); + det->header = header; det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); @@ -64,131 +57,25 @@ psMetadataLookupF32(NULL, header, "FWHM_MIN")); - int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns - int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows + det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); + det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); - psFree(header); - - psArray *table = psFitsReadTable(fits); // Table of interest + det->table = psFitsReadTableAllColumns(fits); // Table of interest if (!table) { - psError(PS_ERR_IO, false, "Unable to read table %d", i); + 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 + psTrace("ppMops.read", 2, "Read %ld rows from %s\n", numGood, (const char*)inNames->data[i]); - 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)); + det->ra = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "RA_PSF")); + det->dec = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "DEC_PSF")); + 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 good rows from %s\n", numGood, (const char*)inNames->data[i]); - - psFree(table); detections->data[i] = det; } Index: /trunk/ppTranslate/src/ppMopsWrite.c =================================================================== --- /trunk/ppTranslate/src/ppMopsWrite.c (revision 28208) +++ /trunk/ppTranslate/src/ppMopsWrite.c (revision 28209) @@ -9,7 +9,7 @@ #include "ppTranslateVersion.h" -bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args) { - psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output); + psTrace("ppMops.write", 1, "Writing %ld extensions to %s", detections->n, args->output); psFits *fits = psFitsOpen(args->output, "w"); // FITS file @@ -19,134 +19,39 @@ } + // 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); - 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); - - 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); + if (!psFitsWriteBlank(fits, header)) { + psError(psErrorCodeLast(), false, "Unable to write primary header"); return false; } - 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]); + psFree(header); + } - table->data[i] = row; - } - if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { - psErrorStackPrint(stderr, "Unable to write table."); - psFree(header); - 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); + if (!psFitsWriteTableAllColumns(fits, det->header, det->table, NULL)) { + psError(psErrorCodeLast(), false, "Unable to write extension %d", i); return false; } - psFree(table); } - - psFree(header); psFitsClose(fits); - psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output); + psTrace("ppMops.write", 1, "Done writing to %s", args->output); return true;