Index: trunk/ppTranslate/src/ppMopsWrite.c =================================================================== --- trunk/ppTranslate/src/ppMopsWrite.c (revision 29567) +++ 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; +}