Index: /branches/pap_mops/ppMops/src/ppMops.h =================================================================== --- /branches/pap_mops/ppMops/src/ppMops.h (revision 25180) +++ /branches/pap_mops/ppMops/src/ppMops.h (revision 25181) @@ -49,7 +49,16 @@ psVector *flags; // psphot flags psVector *diffSkyfileId; // Identifier for source image + psVector *naxis1, *naxis2; // Size of image + psVector *mask; // Mask for detections } 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); + /// Read detections Index: /branches/pap_mops/ppMops/src/ppMopsDetections.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsDetections.c (revision 25180) +++ /branches/pap_mops/ppMops/src/ppMopsDetections.c (revision 25181) @@ -28,4 +28,7 @@ psFree(det->flags); psFree(det->diffSkyfileId); + psFree(det->naxis1); + psFree(det->naxis2); + psFree(det->mask); return; } @@ -56,10 +59,13 @@ det->magErr = psVectorAlloc(num, PS_TYPE_F32); det->extended = psVectorAlloc(num, PS_TYPE_F32); - det->angle = psVectorAlloc(num, PS_TYPE_F64); - det->angleErr = psVectorAlloc(num, PS_TYPE_F64); + det->angle = psVectorAlloc(num, PS_TYPE_F32); + det->angleErr = psVectorAlloc(num, PS_TYPE_F32); det->length = psVectorAlloc(num, PS_TYPE_F32); det->lengthErr = psVectorAlloc(num, PS_TYPE_F32); det->flags = psVectorAlloc(num, PS_TYPE_U32); det->diffSkyfileId = psVectorAlloc(num, PS_TYPE_S64); + det->naxis1 = psVectorAlloc(num, PS_TYPE_S32); + det->naxis2 = psVectorAlloc(num, PS_TYPE_S32); + det->mask = psVectorAlloc(num, PS_TYPE_U8); return det; @@ -67,2 +73,131 @@ +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->extended = psVectorRealloc(det->extended, num); + det->angle = psVectorRealloc(det->angle, num); + det->angleErr = psVectorRealloc(det->angleErr, num); + det->length = psVectorRealloc(det->length, num); + det->lengthErr = psVectorRealloc(det->lengthErr, 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); + + return det; +} + + +bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec, + double raErr, double decErr, float mag, float magErr, float extended, + float angle, float angleErr, float length, float lengthErr, + psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2) +{ + 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->extended, extended); + psVectorAppend(det->angle, angle); + psVectorAppend(det->angleErr, angleErr); + psVectorAppend(det->length, length); + psVectorAppend(det->lengthErr, lengthErr); + psVectorAppend(det->flags, flags); + psVectorAppend(det->diffSkyfileId, diffSkyfileId); + psVectorAppend(det->naxis1, naxis1); + psVectorAppend(det->naxis2, naxis2); + psVectorAppend(det->mask, 0); + 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->extended, source->extended->data.F32[index]); + psVectorAppend(target->angle, source->angle->data.F32[index]); + psVectorAppend(target->angleErr, source->angleErr->data.F32[index]); + psVectorAppend(target->length, source->length->data.F32[index]); + psVectorAppend(target->lengthErr, source->lengthErr->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); + 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->extended->data.F32[num] = det->extended->data.F32[i]; + det->angle->data.F32[num] = det->angle->data.F32[i]; + det->angleErr->data.F32[num] = det->angleErr->data.F32[i]; + det->length->data.F32[num] = det->length->data.F32[i]; + det->lengthErr->data.F32[num] = det->lengthErr->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; + 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->extended->n = num; + det->angle->n = num; + det->angleErr->n = num; + det->length->n = num; + det->lengthErr->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; + return true; +} + Index: /branches/pap_mops/ppMops/src/ppMopsMerge.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsMerge.c (revision 25180) +++ /branches/pap_mops/ppMops/src/ppMopsMerge.c (revision 25181) @@ -5,10 +5,141 @@ #include #include +#include #include "ppMops.h" +#define LEAF_SIZE 4 // Size of leaf +#define MATCH_RADIUS 1.0/3600 // Matching radius + +// Get distance from detection to centre of image +static float mergeDistance(const ppMopsDetections *detections, // Detections of interest + long index // Index for source of interest + ) +{ + float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0; + return PS_SQR(dx) + PS_SQR(dy); +} + + ppMopsDetections *ppMopsMerge(const psArray *detections) { + PS_ASSERT_ARRAY_NON_NULL(detections, NULL); - return NULL; + ppMopsDetections *merged = psMemIncrRefCounter(detections->data[0]); // Merged list + int num = 1; // Number of merged files + for (int i = 1; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections of interest + if (!det) { + continue; + } + num++; + + // 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 (merged->airmass != det->airmass) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", + merged->airmass, det->airmass); + return NULL; + } + if (merged->exptime != det->exptime) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", + merged->exptime, det->exptime); + return NULL; + } + if (merged->posangle != det->posangle) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", + merged->posangle, det->posangle); + return NULL; + } + if (merged->alt != det->alt) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", + merged->alt, det->alt); + return NULL; + } + if (merged->az != det->az) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", + merged->az, det->az); + return NULL; + } + if (merged->mjd != det->mjd) { + 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; + } + + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest + for (int j = 0; j < det->num; j++) { + coords->data.F64[0] = det->ra->data.F64[j]; + coords->data.F64[1] = det->dec->data.F64[j]; + psVector *indices = psTreeAllWithin(tree, coords, MATCH_RADIUS); // Indices for matching sources + if (!indices) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to search for matches"); + psFree(coords); + psFree(tree); + psFree(merged); + return NULL; + } + if (indices->n == 0) { + psFree(indices); + continue; + } + + // Which one do we keep? + float bestDistance = INFINITY; // Best distance to centre + 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 + if (distance < bestDistance) { + bestDistance = distance; + bestIndex = index; + } + } + + float distance = mergeDistance(det, j); // Distance to centre of image + if (distance < bestDistance) { + // Blow away existing sources + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + merged->mask->data.U8[index] = 0xFF; + } + + ppMopsDetectionsCopySingle(merged, det, j); + } + psFree(indices); + } + + psFree(tree); + ppMopsDetectionsPurge(merged); + } + + merged->seeing /= num; + + return merged; } + Index: /branches/pap_mops/ppMops/src/ppMopsRead.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsRead.c (revision 25180) +++ /branches/pap_mops/ppMops/src/ppMopsRead.c (revision 25181) @@ -60,4 +60,7 @@ psMetadataLookupF32(NULL, header, "FWHM_MIN")); + int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows + psFree(header); @@ -85,4 +88,7 @@ det->flags->data.F32[j] = psMetadataLookupU32(NULL, row, "FLAGS"); det->diffSkyfileId->data.F32[j] = diffSkyfileId; + det->naxis1->data.S32[j] = naxis1; + det->naxis2->data.S32[j] = naxis2; + det->mask->data.U8[j] = 0; // Calculate error in RA, Dec Index: /branches/pap_mops/ppMops/src/ppMopsWrite.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsWrite.c (revision 25180) +++ /branches/pap_mops/ppMops/src/ppMopsWrite.c (revision 25181) @@ -8,7 +8,104 @@ #include "ppMops.h" -bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args) +bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) { + psFits *fits = psFitsOpen(args->output, "w"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open output file."); + return false; + } + psMetadata *header = psMetadataAlloc(); // Header to write + psString source = ppMopsSource(), version = ppMopsVersion(); + psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); + psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); + ppMopsVersionHeader(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); + psMetadataAddF64(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", NAN); + psMetadataAddS32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0); + if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write empty table."); + psFree(header); + psFree(row); + 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)", 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)", det->dec->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", + det->decErr->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", det->extended->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", + det->angle->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", + det->angleErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", + det->length->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", + det->lengthErr->data.F32[i]); + psMetadataAddS32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.S32[i]); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", + det->diffSkyfileId->data.S64[i]); + table->data[i] = row; + } + if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(header); + psFree(table); + return false; + } + psFree(table); + } + + psFree(header); + psFitsClose(fits); return true; }