Index: trunk/ppTranslate/src/ppMops.c =================================================================== --- trunk/ppTranslate/src/ppMops.c (revision 32396) +++ trunk/ppTranslate/src/ppMops.c (revision 32406) @@ -3,4 +3,24 @@ #include "ppMops.h" + +void test() +{ + psMetadata *md = psMetadataAlloc(); + + psVector *vec = psVectorAlloc(42, PS_TYPE_S32); + + psMetadataAddVector(md, PS_LIST_TAIL, "TEST", 0, NULL, vec); + + psFree(vec); + psFree(md); + + psLibFinalize(); + + fprintf (stderr, "found %d leaks at %s\n", + psMemCheckLeaks2 (0, + NULL, stdout, false, 500), "ppMops"); + + exit(0); +} /* @@ -39,4 +59,6 @@ psLibInit(NULL); + // test(); + ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments if (!args) { @@ -45,4 +67,5 @@ } + psArray *detections = ppMopsRead(args); // Detections from each input if (!detections) { @@ -51,19 +74,19 @@ } - ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections - psFree(detections); - if (!merged) { + if (!ppMopsPurgeDuplicates(detections)) { psErrorStackPrint(stderr, "Unable to merge detections"); exit(PS_EXIT_SYS_ERROR); } - if (!ppMopsWrite(merged, args)) { + if (!ppMopsWrite(detections, args)) { psErrorStackPrint(stderr, "Unable to write detections"); exit(PS_EXIT_SYS_ERROR); } - psFree(merged); + for (int i = 0; i < detections->n; i++) { + psFree(detections->data[i]); + } + psFree(detections); psFree(args); - psFree(detections); psLibFinalize(); Index: trunk/ppTranslate/src/ppMops.h =================================================================== --- trunk/ppTranslate/src/ppMops.h (revision 32396) +++ trunk/ppTranslate/src/ppMops.h (revision 32406) @@ -62,4 +62,5 @@ typedef struct { + psString component; // skycell_id for these detections psString raBoresight, decBoresight; // RA,Dec of telescope boresight psString filter; // Filter for exposure @@ -70,53 +71,16 @@ double mjd; // Modified Julian Date float seeing; // Seeing of exposure + int naxis1, naxis2; // size of the image long num; // Number of detections + long numGood; // Number of "good" detections + psS64 diffSkyfileId; // unique id for input skyfile + psMetadata *table; // Columns from the input file 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 - psVector *psfInstFlux; // PSF fit instrumental magnitude - psVector *psfInstFluxSig; // Sigma of PSF instrumental magnitude - psVector *apMag; // Magnitude in standard aperture - psVector *apMagRadius; // Radius used for aperture mags - psVector *apMagRaw; // Magnitude in real aperture - psVector *apFlux; // Instrumental flux in standard aperture - psVector *apFluxSig; // Aperture flux error - psVector *peakFluxAsMag; // Peak flux expressed as magnitude - psVector *calPsfMag; // PSF Magnitude using supplied calibration - psVector *calPsfMagSig; // Measured scatter of zero point calibration - psVector *sky; // Sky level - psVector *skySig; // Sigma of sky level - psVector *qualityPerfect; // PSF coverage/quality factor (poor) - psVector *momentsR1; // First radial moment - psVector *momentsRH; // Half radial moment - psVector *kronFlux; // Kron Flux (in 2.5 R1) - psVector *kronFluxErr; // Kron Flux Error - psVector *kronFluxInner; // Kron Flux (in 1.0 R1) - psVector *kronFluxOuter; // Kron Flux (in 4.0 R1) - psVector *diffRP; // Distance to positive match source - psVector *diffSnP; // Signal-to-noise of pos match src - psVector *diffRM; // Distance to negative match source - psVector *diffSnM; // Signal-to-noise of neg match src - psVector *flags2; // psphot flags (group 2) - psVector *ippIdet; // IPP detection identifier index - psVector *nFrames; // Number of frames overlapping source center - psVector *padding; // Padding } ppMopsDetections; -ppMopsDetections *ppMopsDetectionsAlloc(long num); +ppMopsDetections *ppMopsDetectionsAlloc(); /// Copy a detection @@ -130,8 +94,9 @@ /// Merge detections -ppMopsDetections *ppMopsMerge(const psArray *detections); +// ppMopsDetections *ppMopsMerge(const psArray *detections); +bool ppMopsPurgeDuplicates(const psArray *detections); /// Write detections -bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args); /// Get the version contained in EXTTYPE of the "SkyChip.psf" table: Index: trunk/ppTranslate/src/ppMopsDetections.c =================================================================== --- trunk/ppTranslate/src/ppMopsDetections.c (revision 32396) +++ trunk/ppTranslate/src/ppMopsDetections.c (revision 32406) @@ -8,7 +8,9 @@ static void mopsDetectionsFree(ppMopsDetections *det) { + psFree(det->component); psFree(det->raBoresight); psFree(det->decBoresight); psFree(det->filter); + psFree(det->table); psFree(det->x); psFree(det->y); @@ -17,62 +19,13 @@ 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); - psFree(det->psfInstFlux); - psFree(det->psfInstFluxSig); - psFree(det->apMag); - psFree(det->apMagRadius); - psFree(det->apMagRaw); - psFree(det->apFlux); - psFree(det->apFluxSig); - psFree(det->peakFluxAsMag); - psFree(det->calPsfMag); - psFree(det->calPsfMagSig); - psFree(det->sky); - psFree(det->skySig); - psFree(det->qualityPerfect); - psFree(det->momentsR1); - psFree(det->momentsRH); - psFree(det->kronFlux); - psFree(det->kronFluxErr); - psFree(det->kronFluxInner); - psFree(det->kronFluxOuter); - psFree(det->diffRP); - psFree(det->diffSnP); - psFree(det->diffRM); - psFree(det->diffSnM); - psFree(det->flags2); - psFree(det->ippIdet); - psFree(det->nFrames); - psFree(det->padding); return; } -ppMopsDetections *ppMopsDetectionsAlloc(long num) +ppMopsDetections *ppMopsDetectionsAlloc() { ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree); + det->component = NULL; det->raBoresight = NULL; det->decBoresight = NULL; @@ -86,396 +39,13 @@ 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->psfInstFlux = psVectorAllocEmpty(num, PS_TYPE_F32); - det->psfInstFluxSig = psVectorAllocEmpty(num, PS_TYPE_F32); - det->apMag = psVectorAllocEmpty(num, PS_TYPE_F32); - det->apMagRadius = psVectorAllocEmpty(num, PS_TYPE_F32); - det->apMagRaw = psVectorAllocEmpty(num, PS_TYPE_F32); - det->apFlux = psVectorAllocEmpty(num, PS_TYPE_F32); - det->apFluxSig = psVectorAllocEmpty(num, PS_TYPE_F32); - det->peakFluxAsMag = psVectorAllocEmpty(num, PS_TYPE_F32); - det->calPsfMag = psVectorAllocEmpty(num, PS_TYPE_F32); - det->calPsfMagSig = psVectorAllocEmpty(num, PS_TYPE_F32); - det->sky = psVectorAllocEmpty(num, PS_TYPE_F32); - det->skySig = psVectorAllocEmpty(num, PS_TYPE_F32); - det->qualityPerfect = psVectorAllocEmpty(num, PS_TYPE_F32); - det->momentsR1 = psVectorAllocEmpty(num, PS_TYPE_F32); - det->momentsRH = psVectorAllocEmpty(num, PS_TYPE_F32); - det->kronFlux = psVectorAllocEmpty(num, PS_TYPE_F32); - det->kronFluxErr = psVectorAllocEmpty(num, PS_TYPE_F32); - det->kronFluxInner = psVectorAllocEmpty(num, PS_TYPE_F32); - det->kronFluxOuter = psVectorAllocEmpty(num, PS_TYPE_F32); - det->diffRP = psVectorAllocEmpty(num, PS_TYPE_F32); - det->diffSnP = psVectorAllocEmpty(num, PS_TYPE_F32); - det->diffRM = psVectorAllocEmpty(num, PS_TYPE_F32); - det->diffSnM = psVectorAllocEmpty(num, PS_TYPE_F32); - det->flags2 = psVectorAllocEmpty(num, PS_TYPE_U32); - det->ippIdet = psVectorAllocEmpty(num, PS_TYPE_U32); - det->nFrames = psVectorAllocEmpty(num, PS_TYPE_U16); - det->padding = psVectorAllocEmpty(num, PS_TYPE_S16); + det->table = NULL; + det->x = NULL; + det->y = NULL; + det->ra = NULL; + det->dec = NULL; + det->raErr = NULL; + det->decErr = NULL; + det->mask = NULL; + det->diffSkyfileId = 0; 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); - det->psfInstFlux = psVectorRealloc(det->psfInstFlux, num); - det->psfInstFluxSig = psVectorRealloc(det->psfInstFluxSig, num); - det->apMag = psVectorRealloc(det->apMag, num); - det->apMagRadius = psVectorRealloc(det->apMagRadius, num); - det->apMagRaw = psVectorRealloc(det->apMagRadius, num); - det->apFlux = psVectorRealloc(det->apFlux, num); - det->apFluxSig = psVectorRealloc(det->apFluxSig, num); - det->peakFluxAsMag = psVectorRealloc(det->peakFluxAsMag, num); - det->calPsfMag = psVectorRealloc(det->calPsfMag, num); - det->calPsfMagSig = psVectorRealloc(det->calPsfMagSig, num); - det->sky = psVectorRealloc(det->sky, num); - det->skySig = psVectorRealloc(det->skySig, num); - det->qualityPerfect = psVectorRealloc(det->qualityPerfect, num); - det->momentsR1 = psVectorRealloc(det->momentsR1, num); - det->momentsRH = psVectorRealloc(det->momentsRH, num); - det->kronFlux = psVectorRealloc(det->kronFlux, num); - det->kronFluxErr = psVectorRealloc(det->kronFluxErr, num); - det->kronFluxInner = psVectorRealloc(det->kronFluxInner, num); - det->kronFluxOuter = psVectorRealloc(det->kronFluxOuter, num); - det->diffRP = psVectorRealloc(det->diffRP, num); - det->diffSnP = psVectorRealloc(det->diffSnP, num); - det->diffRM = psVectorRealloc(det->diffRM, num); - det->diffSnM = psVectorRealloc(det->diffSnM, num); - det->flags2 = psVectorRealloc(det->flags2, num); - det->ippIdet = psVectorRealloc(det->ippIdet, num); - det->nFrames = psVectorRealloc(det->nFrames, num); - det->padding = psVectorRealloc(det->padding, 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, - float psfInstFlux, float psfInstFluxSig, - float apMag, float apMagRadius, float apMagRaw, float apFlux, float apFluxSig, - float peakFluxAsMag, float calPsfMag, float calPsfMagSig, - float sky, float skySig, float qualityPerfect, - float momentsR1, float momentsRH, - float kronFlux, float kronFluxErr, float kronFluxInner, float kronFluxOuter, - float diffRP, float diffSnP, float diffRM, float diffSnM, - psU32 flags2, psU32 ippIdet, psU16 nFrames, psS16 padding) -{ - 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); - psVectorAppend(det->psfInstFlux, psfInstFlux); - psVectorAppend(det->psfInstFluxSig, psfInstFluxSig); - psVectorAppend(det->apMag, apMag); - psVectorAppend(det->apMagRadius, apMagRadius); - psVectorAppend(det->apMagRaw, apMagRaw); - psVectorAppend(det->apFlux, apFlux); - psVectorAppend(det->apFluxSig, apFluxSig); - psVectorAppend(det->peakFluxAsMag, peakFluxAsMag); - psVectorAppend(det->calPsfMag, calPsfMag); - psVectorAppend(det->calPsfMagSig, calPsfMagSig); - psVectorAppend(det->sky, sky); - psVectorAppend(det->skySig, skySig); - psVectorAppend(det->qualityPerfect, qualityPerfect); - psVectorAppend(det->momentsR1, momentsR1); - psVectorAppend(det->momentsRH, momentsRH); - psVectorAppend(det->kronFlux, kronFlux); - psVectorAppend(det->kronFluxErr, kronFluxErr); - psVectorAppend(det->kronFluxInner, kronFluxInner); - psVectorAppend(det->kronFluxOuter, kronFluxOuter); - psVectorAppend(det->diffRP, diffRP); - psVectorAppend(det->diffSnP, diffSnP); - psVectorAppend(det->diffRM, diffRM); - psVectorAppend(det->diffSnM, diffSnM); - psVectorAppend(det->flags2, flags2); - psVectorAppend(det->ippIdet, ippIdet); - psVectorAppend(det->nFrames, nFrames); - psVectorAppend(det->padding, padding); - 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]); - psVectorAppend(target->psfInstFlux, source->psfInstFlux->data.F32[index]); - psVectorAppend(target->psfInstFluxSig, source->psfInstFluxSig->data.F32[index]); - psVectorAppend(target->apMag, source->apMag->data.F32[index]); - psVectorAppend(target->apMagRadius, source->apMagRadius->data.F32[index]); - psVectorAppend(target->apMagRaw, source->apMagRaw->data.F32[index]); - psVectorAppend(target->apFlux, source->apFlux->data.F32[index]); - psVectorAppend(target->apFluxSig, source->apFluxSig->data.F32[index]); - psVectorAppend(target->peakFluxAsMag, source->peakFluxAsMag->data.F32[index]); - psVectorAppend(target->calPsfMag, source->calPsfMag->data.F32[index]); - psVectorAppend(target->calPsfMagSig, source->calPsfMagSig->data.F32[index]); - psVectorAppend(target->sky, source->sky->data.F32[index]); - psVectorAppend(target->skySig, source->skySig->data.F32[index]); - psVectorAppend(target->qualityPerfect, source->qualityPerfect->data.F32[index]); - psVectorAppend(target->momentsR1, source->momentsR1->data.F32[index]); - psVectorAppend(target->momentsRH, source->momentsRH->data.F32[index]); - psVectorAppend(target->kronFlux, source->kronFlux->data.F32[index]); - psVectorAppend(target->kronFluxErr, source->kronFluxErr->data.F32[index]); - psVectorAppend(target->kronFluxInner, source->kronFluxInner->data.F32[index]); - psVectorAppend(target->kronFluxOuter, source->kronFluxOuter->data.F32[index]); - psVectorAppend(target->diffRP, source->diffRP->data.F32[index]); - psVectorAppend(target->diffSnP, source->diffSnP->data.F32[index]); - psVectorAppend(target->diffRM, source->diffRM->data.F32[index]); - psVectorAppend(target->diffSnM, source->diffSnM->data.F32[index]); - psVectorAppend(target->flags2, source->flags2->data.U32[index]); - psVectorAppend(target->ippIdet, source->ippIdet->data.U32[index]); - psVectorAppend(target->nFrames, source->nFrames->data.U16[index]); - psVectorAppend(target->padding, source->padding->data.S16[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]; - det->psfInstFlux->data.F32[num] = det->psfInstFlux->data.F32[i]; - det->psfInstFluxSig->data.F32[num] = det->psfInstFluxSig->data.F32[i]; - det->apMag->data.F32[num] = det->apMag->data.F32[i]; - det->apMagRadius->data.F32[num] = det->apMagRadius->data.F32[i]; - det->apMagRaw->data.F32[num] = det->apMagRaw->data.F32[i]; - det->apFlux->data.F32[num] = det->apFlux->data.F32[i]; - det->apFluxSig->data.F32[num] = det->apFluxSig->data.F32[i]; - det->peakFluxAsMag->data.F32[num] = det->peakFluxAsMag->data.F32[i]; - det->calPsfMag->data.F32[num] = det->calPsfMag->data.F32[i]; - det->calPsfMagSig->data.F32[num] = det->calPsfMagSig->data.F32[i]; - det->sky->data.F32[num] = det->sky->data.F32[i]; - det->skySig->data.F32[num] = det->skySig->data.F32[i]; - det->qualityPerfect->data.F32[num] = det->qualityPerfect->data.F32[i]; - det->momentsR1->data.F32[num] = det->momentsR1->data.F32[i]; - det->momentsRH->data.F32[num] = det->momentsRH->data.F32[i]; - det->kronFlux->data.F32[num] = det->kronFlux->data.F32[i]; - det->kronFluxErr->data.F32[num] = det->kronFluxErr->data.F32[i]; - det->kronFluxInner->data.F32[num] = det->kronFluxInner->data.F32[i]; - det->kronFluxOuter->data.F32[num] = det->kronFluxOuter->data.F32[i]; - det->diffRP->data.F32[num] = det->diffRP->data.F32[i]; - det->diffSnP->data.F32[num] = det->diffSnP->data.F32[i]; - det->diffRM->data.F32[num] = det->diffRM->data.F32[i]; - det->diffSnM->data.F32[num] = det->diffSnM->data.F32[i]; - det->flags2->data.U32[num] = det->flags2->data.U32[i]; - det->ippIdet->data.U32[num] = det->ippIdet->data.U32[i]; - det->nFrames->data.U16[num] = det->nFrames->data.U16[i]; - det->padding->data.S16[num] = det->padding->data.S16[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; - det->psfInstFlux->n = num; - det->psfInstFluxSig->n = num; - det->apMag->n = num; - det->apMagRadius->n = num; - det->apMagRaw->n = num; - det->apFlux->n = num; - det->apFluxSig->n = num; - det->peakFluxAsMag->n = num; - det->calPsfMag->n = num; - det->calPsfMagSig->n = num; - det->sky->n = num; - det->skySig->n = num; - det->qualityPerfect->n = num; - det->momentsR1->n = num; - det->momentsRH->n = num; - det->kronFlux->n = num; - det->kronFluxErr->n = num; - det->kronFluxInner->n = num; - det->kronFluxOuter->n = num; - det->diffRP->n = num; - det->diffSnP->n = num; - det->diffRM->n = num; - det->diffSnM->n = num; - det->flags2->n = num; - det->ippIdet->n = num; - det->nFrames->n = num; - det->padding->n = num; - return true; -} Index: trunk/ppTranslate/src/ppMopsMerge.c =================================================================== --- trunk/ppTranslate/src/ppMopsMerge.c (revision 32396) +++ trunk/ppTranslate/src/ppMopsMerge.c (revision 32406) @@ -17,4 +17,12 @@ #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 @@ -22,94 +30,148 @@ ) { - float dx = (float) (detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0); - float dy = (float) (detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0); - return PS_SQR(dx) + PS_SQR(dy); + float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0; + return PS_SQR(dx) + PS_SQR(dy); } -ppMopsDetections *ppMopsMerge(const psArray *detections) +bool ppMopsPurgeDuplicates(const psArray *detections) { PS_ASSERT_ARRAY_NON_NULL(detections, NULL); - psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); - - ppMopsDetections *merged = NULL; // Merged list - int num = 0; // Number of merged files - psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest + int numInputs = detections->n; // Number of inputs + psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs); + + 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 *det = detections->data[i]; // Detections of interest + if (!det) { + psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i); + continue; + } else if (det->num == 0) { + psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i); + continue; + } + + // 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; + } + } + + 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; + } + num += det->num; + } + + 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; + } + for (int i = 0; i < detections->n; i++) { ppMopsDetections *det = detections->data[i]; // Detections of interest if (!det) { - psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i); continue; } else if (det->num == 0) { - psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i); - 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); - - // 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; - } - - 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; - } - 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; - } - - merged->seeing += det->seeing; // Taking average - - 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; - } - + continue; + } + 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++) { + // XXX: added by Bill + if (det->mask->data.U8[j]) { + // we marked this source marked as bad when we read it, mark it as a duplicate + dupes->data.U8[j] = 0xFF; + continue; + } + coords->data.F64[0] = det->ra->data.F64[j]; coords->data.F64[1] = det->dec->data.F64[j]; @@ -119,14 +181,14 @@ psFree(coords); psFree(tree); - psFree(merged); - return NULL; - } - if (indices->n == 0) { - psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i); + return false; + } + 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); - ppMopsDetectionsCopySingle(merged, det, j); continue; } - psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); // Which one do we keep? @@ -134,6 +196,16 @@ long bestIndex = -1; // Index with best distance for (int k = 0; k < indices->n; k++) { - long index = indices->data.S64[k]; // Index of point - float distance = mergeDistance(merged, index); // Distance to centre of image + 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 if (distance < bestDistance) { bestDistance = distance; @@ -142,36 +214,100 @@ } - 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; + 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]++; } - ppMopsDetectionsCopySingle(merged, det, j); - } else { - psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); } psFree(indices); } - - psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num); - - psFree(tree); - ppMopsDetectionsPurge(merged); - } - psFree(coords); - - if (num == 0) { - //All detections were NULL?! - psTrace("ppMops.merge", 3, "All %ld detections were NULL\n", detections->n); - return NULL; - } - psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); - - merged->seeing /= (float) num; - - return merged; + 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); + +#define VECTOR_PURGE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + long j = 0; \ + for (long i = 0; i < vector->n; i++) { \ + if (!dupes->data.U8[i]) { \ + if (i == j) { \ + j++; \ + continue; \ + } \ + vector->data.TYPE[j++] = vector->data.TYPE[i]; \ + } \ + } \ + vector->n = j; \ + } \ + break + + psVector *dupes = duplicates->data[i]; // Duplicates + psArray *table = psListToArray(det->table->list); // Table of data + long newLength = -1; + for (int t = 0; t < table->n; t++) { + psMetadataItem *item = table->data[t]; // 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); + } + if (newLength == -1) { + newLength = vector->n; + } else if (vector->n != newLength) { + psAbort("Unexpected new length found : %ld expected: %ld", + vector->n, newLength); + } + + } + // XXX IS this safe? Perhaps save in numGood? + det->num = newLength; + psFree(table); + } + psFree(dupNum); + psFree(duplicates); + + return true; } Index: trunk/ppTranslate/src/ppMopsRead.c =================================================================== --- trunk/ppTranslate/src/ppMopsRead.c (revision 32396) +++ trunk/ppTranslate/src/ppMopsRead.c (revision 32406) @@ -20,10 +20,12 @@ psArray *detections = psArrayAlloc(num); // Array of detections, to return for (int i = 0; i < num; i++) { - psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file - + const char *name = inNames->data[i]; + + psFits *fits = psFitsOpen(name, "r"); // FITS file if (!fits) { psError(PS_ERR_IO, false, "Unable to open input %d", i); return false; } + psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header if (!header) { @@ -74,4 +76,8 @@ } ppMopsDetections *det = ppMopsDetectionsAlloc(size); + detections->data[i] = det; + det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension + det->num = size; + det->diffSkyfileId = diffSkyfileId; psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); @@ -90,18 +96,52 @@ 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"); // Number of columns + det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows psFree(header); -#ifndef USE_OLD_READ_TABLE - // use new fits table functions - - psFitsTable *table = psFitsReadTableNew(fits); // Table of interest + psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest if (!table) { psError(PS_ERR_IO, false, "Unable to read table %d", i); - return false; - } + return NULL; + } + det->table = table; psFitsClose(fits); + if (args->version == 0) { + if (skyChipPsfVersion < 2) { + // XXX: TODO: Do we need to add dummy vectors for the missing columns? + } + } + + psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); + psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); + + det->raErr = psVectorAlloc(size, PS_TYPE_F64); + det->decErr = psVectorAlloc(size, PS_TYPE_F64); + det->mask = psVectorAlloc(size, PS_TYPE_U8); + + // convert ra and dec to radians for use in the purge duplicates function + 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 NULL; + } + + // Add our new vectors to the table so that duplicates and masked items may be purged + psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr); + psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr); + + psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); + + psVector *mag = psMetadataLookupVector(NULL, table, "PSF_INST_MAG"); + psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG"); + psVector *xErrV = psMetadataLookupVector(NULL, table, "X_PSF_SIG"); + psVector *yErrV = psMetadataLookupVector(NULL, table, "Y_PSF_SIG"); + psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE"); + psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE"); + psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS"); double plateScale = 0.0; // Plate scale @@ -109,109 +149,21 @@ for (long row = 0; row < size; row++) { - psU32 flags = psFitsTableGetU32(NULL, table, row, "FLAGS"); + psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS"); if (flags & SOURCE_MASK) { psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags); + det->mask->data.U8[row] = 0xFF; continue; } - det->x->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "X_PSF"); - det->y->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "Y_PSF"); - det->ra->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "RA_PSF")); - det->dec->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "DEC_PSF")); - det->mag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG"); - det->magErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG_SIG"); - det->chi2->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_CHISQ"); - det->dof->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NDOF"); - det->cr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CR_NSIGMA"); - det->extended->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "EXT_NSIGMA"); - det->psfMajor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MAJOR"); - det->psfMinor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MINOR"); - det->psfTheta->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_THETA"); - det->quality->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF"); - det->numPix->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NPIX"); - det->xxMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XX"); - det->xyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XY"); - det->yyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_YY"); - det->flags->data.U32[numGood] = psFitsTableGetU32(NULL, table, 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] = psFitsTableGetS32(NULL, table, row, "DIFF_NPOS"); - det->fPos->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_FRATIO"); - det->nRatioBad->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_BAD"); - det->nRatioMask->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_MASK"); - det->nRatioAll->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_ALL"); - - //Additions of 2010-10-25 - if (args->version == 2) { - //Values are set only if the version is 2 - if (skyChipPsfVersion == 2) { - det->psfInstFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX"); - det->psfInstFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX_SIG"); - det->apMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG"); - det->apMagRaw->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RAW"); - det->apMagRadius->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RADIUS"); - det->apFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX"); - det->apFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX_SIG"); - det->peakFluxAsMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PEAK_FLUX_AS_MAG"); - det->calPsfMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG"); - det->calPsfMagSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG_SIG"); - det->sky->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY"); - det->skySig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY_SIGMA"); - det->qualityPerfect->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF_PERFECT"); - det->momentsR1->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_R1"); - det->momentsRH->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_RH"); - det->kronFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX"); - det->kronFluxErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_ERR"); - det->kronFluxInner->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_INNER"); - det->kronFluxOuter->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_OUTER"); - det->diffRP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_P"); - det->diffSnP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_P"); - det->diffRM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_M"); - det->diffSnM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_M"); - det->flags2->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "FLAGS2"); - det->ippIdet->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "IPP_IDET"); - det->nFrames->data.U16[numGood] = psFitsTableGetU16(NULL, table, row, "N_FRAMES"); - det->padding->data.S16[numGood] = psFitsTableGetS16(NULL, table, row, "PADDING"); - } else { - det->psfInstFlux->data.F32[numGood] = NAN; - det->psfInstFluxSig->data.F32[numGood] = NAN; - det->apMag->data.F32[numGood] = NAN; - det->apMagRaw->data.F32[numGood] = NAN; - det->apMagRadius->data.F32[numGood] = NAN; - det->apFlux->data.F32[numGood] = NAN; - det->apFluxSig->data.F32[numGood] = NAN; - det->peakFluxAsMag->data.F32[numGood] = NAN; - det->calPsfMag->data.F32[numGood] = NAN; - det->calPsfMagSig->data.F32[numGood] = NAN; - det->sky->data.F32[numGood] = NAN; - det->skySig->data.F32[numGood] = NAN; - det->qualityPerfect->data.F32[numGood] = NAN; - det->momentsR1->data.F32[numGood] = NAN; - det->momentsRH->data.F32[numGood] = NAN; - det->kronFlux->data.F32[numGood] = NAN; - det->kronFluxErr->data.F32[numGood] = NAN; - det->kronFluxInner->data.F32[numGood] = NAN; - det->kronFluxOuter->data.F32[numGood] = NAN; - det->diffRP->data.F32[numGood] = NAN; - det->diffSnP->data.F32[numGood] = NAN; - det->diffRM->data.F32[numGood] = NAN; - det->diffSnM->data.F32[numGood] = NAN; - det->flags2->data.U32[numGood] = 0; - det->ippIdet->data.U32[numGood] = 0; - det->nFrames->data.U16[numGood] = 0; - det->padding->data.S16[numGood] = 0; - } - } - // Calculate error in RA, Dec - double xErr = psFitsTableGetF64(NULL, table, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32 - double yErr = psFitsTableGetF64(NULL, table, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32 - double scale = psFitsTableGetF64(NULL, table, row, "PLTSCALE"); //SC: Warning! Promotion of F32 - double angle = psFitsTableGetF64(NULL, table, row, "POSANGLE"); //SC: Warning! Promotion of F32 - - 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]) || + + double xErr = xErrV->data.F32[row]; + double yErr = yErrV->data.F32[row]; + double scale = scaleV->data.F32[row]; + double angle = angleV->data.F32[row]; + + if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) || + !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) || + !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) || !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { psTrace("ppMops.read", 10, @@ -219,8 +171,9 @@ "%f %f %lf %lf %f %f %f %f %f %f", row, 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], + det->x->data.F32[row], det->y->data.F32[row], + det->ra->data.F64[row], det->dec->data.F64[row], + mag->data.F32[row], magErr->data.F32[row], xErr, yErr, scale, angle); + det->mask->data.U8[row] = 0xFF; continue; } @@ -231,226 +184,21 @@ 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; + det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[row] = 0; plateScale += scale; numGood++; } -#else - // OLD way of reading fits tables read it into an array of metadata objects, one for each row. Uses lots and lots - // and lots of memory. - psArray *table = psFitsReadTable(fits); // Table of interest - if (!table) { - psError(PS_ERR_IO, 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"); - - //Additions of 2010-10-25 - if (args->version == 2) { - //Values are set only if the version is 2 - if (skyChipPsfVersion == 2) { - det->psfInstFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX"); - det->psfInstFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX_SIG"); - det->apMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG"); - det->apMagRaw->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RAW"); - det->apMagRadius->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RADIUS"); - det->apFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX"); - det->apFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX_SIG"); - det->peakFluxAsMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PEAK_FLUX_AS_MAG"); - det->calPsfMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG"); - det->calPsfMagSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG_SIG"); - det->sky->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY"); - det->skySig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY_SIGMA"); - det->qualityPerfect->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF_PERFECT"); - det->momentsR1->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_R1"); - det->momentsRH->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_RH"); - det->kronFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX"); - det->kronFluxErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_ERR"); - det->kronFluxInner->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_INNER"); - det->kronFluxOuter->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_OUTER"); - det->diffRP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_P"); - det->diffSnP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_P"); - det->diffRM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_M"); - det->diffSnM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_M"); - det->flags2->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS2"); - det->ippIdet->data.U32[numGood] = psMetadataLookupU32(NULL, row, "IPP_IDET"); - det->nFrames->data.U16[numGood] = psMetadataLookupU16(NULL, row, "N_FRAMES"); - det->padding->data.S16[numGood] = psMetadataLookupS16(NULL, row, "PADDING"); - } else { - det->psfInstFlux->data.F32[numGood] = NAN; - det->psfInstFluxSig->data.F32[numGood] = NAN; - det->apMag->data.F32[numGood] = NAN; - det->apMagRaw->data.F32[numGood] = NAN; - det->apMagRadius->data.F32[numGood] = NAN; - det->apFlux->data.F32[numGood] = NAN; - det->apFluxSig->data.F32[numGood] = NAN; - det->peakFluxAsMag->data.F32[numGood] = NAN; - det->calPsfMag->data.F32[numGood] = NAN; - det->calPsfMagSig->data.F32[numGood] = NAN; - det->sky->data.F32[numGood] = NAN; - det->skySig->data.F32[numGood] = NAN; - det->qualityPerfect->data.F32[numGood] = NAN; - det->momentsR1->data.F32[numGood] = NAN; - det->momentsRH->data.F32[numGood] = NAN; - det->kronFlux->data.F32[numGood] = NAN; - det->kronFluxErr->data.F32[numGood] = NAN; - det->kronFluxInner->data.F32[numGood] = NAN; - det->kronFluxOuter->data.F32[numGood] = NAN; - det->diffRP->data.F32[numGood] = NAN; - det->diffSnP->data.F32[numGood] = NAN; - det->diffRM->data.F32[numGood] = NAN; - det->diffSnM->data.F32[numGood] = NAN; - det->flags2->data.U32[numGood] = 0; - det->ippIdet->data.U32[numGood] = 0; - det->nFrames->data.U16[numGood] = 0; - det->padding->data.S16[numGood] = 0; - } - } - - // Calculate error in RA, Dec - double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32 - double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32 - double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); //SC: Warning! Promotion of F32 - double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); //SC: Warning! Promotion of F32 - - 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++; - } -#endif det->seeing *= ((float) plateScale) / ((float) 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->psfInstFlux->n = numGood; - det->psfInstFluxSig->n = numGood; - det->apMag->n = numGood; - det->apMagRaw->n = numGood; - det->apMagRadius->n = numGood; - det->apFlux->n = numGood; - det->apFluxSig->n = numGood; - det->peakFluxAsMag->n = numGood; - det->calPsfMag->n = numGood; - det->calPsfMagSig->n = numGood; - det->sky->n = numGood; - det->skySig->n = numGood; - det->qualityPerfect->n = numGood; - det->momentsR1->n = numGood; - det->momentsRH->n = numGood; - det->kronFlux->n = numGood; - det->kronFluxErr->n = numGood; - det->kronFluxInner->n = numGood; - det->kronFluxOuter->n = numGood; - det->diffRP->n = numGood; - det->diffSnP->n = numGood; - det->diffRM->n = numGood; - det->diffSnM->n = numGood; - det->flags2->n = numGood; - det->ippIdet->n = numGood; - det->nFrames->n = numGood; - det->padding->n = numGood; - - det->num = numGood; + // Are we using numGood for anything outside of this function? + det->numGood = numGood; if (isfinite(args->zp) && numGood > 0) { - psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); - } - - psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); - - psFree(table); - detections->data[i] = det; + psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); + } + + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name); } Index: trunk/ppTranslate/src/ppMopsWrite.c =================================================================== --- trunk/ppTranslate/src/ppMopsWrite.c (revision 32396) +++ trunk/ppTranslate/src/ppMopsWrite.c (revision 32406) @@ -9,7 +9,9 @@ #include "ppTranslateVersion.h" -bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) +static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32); +static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName); + +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args) { - psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output); psFits *fits = psFitsOpen(args->output, "w"); // FITS file @@ -36,4 +38,6 @@ psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); + // Get these header words from the first input + ppMopsDetections *det = detections->data[0]; 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); @@ -55,5 +59,15 @@ psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion); - if (det->num == 0) { + // Find the total number of detections + + long total = 0; + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + total += det->num; + } + + psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output); + + if (total == 0) { // Write dummy table psMetadata *row = psMetadataAlloc(); // Output row @@ -151,114 +165,70 @@ 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]); - - if (args->version == 2) { - // Write data of version 2 (see ICD) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", - det->ippIdet->data.U32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", - det->psfInstFlux->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", - det->psfInstFluxSig->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", - det->apMag->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", - det->apMagRaw->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", - det->apMagRadius->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", - det->apFlux->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", - det->apFluxSig->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", - det->peakFluxAsMag->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", - det->calPsfMag->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", - det->calPsfMagSig->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", - det->sky->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", - det->skySig->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", - det->qualityPerfect->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", - det->momentsR1->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", - det->momentsRH->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", - det->kronFlux->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", - det->kronFluxErr->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", - det->kronFluxInner->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", - det->kronFluxOuter->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", - det->diffRP->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", - det->diffSnP->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", - det->diffRM->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", - det->diffSnM->data.F32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags (group 2)", - det->flags2->data.U32[i]); - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", - det->nFrames->data.U16[i]); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", - det->padding->data.S16[i]); - } - - //Update with the table with the current row - table->data[i] = row; - } - if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { - psErrorStackPrint(stderr, "Unable to write table."); - psFree(header); - psFree(table); + +#define addColumn(_outName, _inName, _convertTo32) \ + if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \ + psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \ + return false; \ + } + + // Allocate the output table + psMetadata *table = psMetadataAlloc(); + addColumn("RA", "RA_PSF", 0); + addColumn("RA_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG + addColumn("DEC", "DEC_PSF", 0); + addColumn("DEC_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG + addColumn("MAG", "PSF_INST_MAG", 0); + addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0); + addColumn("PSF_CHI2", "PSF_CHISQ", 0); + addColumn("PSF_DOF", "PSF_NDOF", 1); + addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0); + addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0); + addColumn("PSF_MAJOR", NULL, 0); + addColumn("PSF_MINOR", NULL, 0); + addColumn("PSF_THETA", NULL, 0); + addColumn("PSF_QUALITY", "PSF_QF", 0); + addColumn("PSF_NPIX", NULL, 1); + addColumn("MOMENTS_XX", NULL, 0); + addColumn("MOMENTS_XY", NULL, 0); + addColumn("MOMENTS_YY", NULL, 0); + addColumn("N_POS", "DIFF_NPOS", 1); + addColumn("F_POS", "DIFF_FRATIO", 0); + addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0); + addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0); + addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0); + addColumn("FLAGS", "FLAGS", 1); + addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID"); + if (args->version == 2) { + addColumn("IPP_IDET", NULL, 1); + addColumn("PSF_INST_FLUX", NULL, 0); + addColumn("PSF_INST_FLUX_SIG", NULL, 0); + addColumn("AP_MAG", NULL, 0); + addColumn("AP_MAG_RAW", NULL, 0); + addColumn("AP_MAG_RADIUS", NULL, 0); + addColumn("AP_FLUX", NULL, 0); + addColumn("AP_FLUX_SIG", NULL, 0); + addColumn("PEAK_FLUX_AS_MAG", NULL, 0); + addColumn("CAL_PSF_MAG", NULL, 0); + addColumn("CAL_PSF_MAG_SIG", NULL, 0); + addColumn("SKY", NULL, 0); + addColumn("SKY_SIGMA", NULL, 0); + addColumn("PSF_QF_PERFECT", NULL, 0); + addColumn("MOMENTS_R1", NULL, 0); + addColumn("MOMENTS_RH", NULL, 0); + addColumn("KRON_FLUX", NULL, 0); + addColumn("KRON_FLUX_ERR", NULL, 0); + addColumn("KRON_FLUX_INNER", NULL, 0); + addColumn("KRON_FLUX_OUTER", NULL, 0); + addColumn("DIFF_R_P", NULL, 0); + addColumn("DIFF_SN_P", NULL, 0); + addColumn("DIFF_R_M", NULL, 0); + addColumn("DIFF_SN_M", NULL, 0); + addColumn("FLAGS2", NULL, 1); + addColumn("IPP_IDET", NULL, 0); + addColumn("N_FRAMES", NULL, 0); + addColumn("PADDING", NULL, 0); + } + if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) { + psError(psErrorCodeLast(), false, "Unable to write table"); return false; } @@ -273,2 +243,130 @@ return true; } + +static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32) +{ + if (inColumnName == NULL) { + inColumnName = outColumnName; + } + + psVector *out = NULL; + if (convertTo32) { + // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers + // We want to write 32 bits to the output. + int next = 0; + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (det->num == 0) { + // no detections survived for this input + continue; + } + psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName); + if (!in) { + psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName); + return false; + } + if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) { + psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d", + inColumnName, in->type.type); + return false; + } + if (out == NULL) { + // First time through set up the output vector and the copy parameters + if (in->type.type == PS_TYPE_S64) { + out = psVectorAlloc(outputSize, PS_TYPE_S32); + } else { + out = psVectorAlloc(outputSize, PS_TYPE_U32); + } + } + for (long d=0; d < det->num; d++) { + if (in->type.type == PS_TYPE_S64) { + out->data.S32[next++] = in->data.S64[d]; + } else { + out->data.U32[next++] = in->data.U64[d]; + } + } + } + } else { + void *next = NULL; + int elementSize = 0; // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (det->num == 0) { + // no detections survived for this input + continue; + } + psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName); + if (!in) { + psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName); + return false; + } + if (out == NULL) { + // First time through set up the output vector and the copy parameters + out = psVectorAlloc(outputSize, in->type.type); + next = (void *) out->data.U8; + switch (in->type.type) { + case PS_TYPE_S8: + elementSize = sizeof(psS8); + break; + case PS_TYPE_U8: + elementSize = sizeof(psU8); + break; + case PS_TYPE_S16: + elementSize = sizeof(psS16); + break; + case PS_TYPE_U16: + elementSize = sizeof(psU16); + break; + case PS_TYPE_S32: + elementSize = sizeof(psS32); + break; + case PS_TYPE_U32: + elementSize = sizeof(psU32); + break; + case PS_TYPE_S64: + elementSize = sizeof(psS64); + break; + case PS_TYPE_U64: + elementSize = sizeof(psU64); + break; + case PS_TYPE_F32: + elementSize = sizeof(psF32); + break; + case PS_TYPE_F64: + elementSize = sizeof(psF64); + break; + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type); + return false; + + } + } + // We are doing nasty things here so we can use memcpy. + // It would be safer to do a proper loop over the elements. + long toCopy = det->num * elementSize; + memcpy(next, in->data.U8, toCopy); + next += toCopy; + } + } + + // Finally add the new column to the output table + psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out); + psFree(out); // drop reference + + return true; +} +static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName) +{ + psVector *out = psVectorAlloc(total, PS_TYPE_S64); + long next = 0; + for (long i = 0; in; i++) { + ppMopsDetections *det = detections->data[i]; + psS64 diffSkyfileId = det->diffSkyfileId; + for (long j = 0; j < det->num; j++) { + out->data.S64[next++] = diffSkyfileId; + } + } + psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out); + psFree(out); + return true; +}