IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 2, 2012, 10:32:09 AM (14 years ago)
Author:
Serge CHASTEL
Message:

MOPS trails fitting integration

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsWrite.c

    r32415 r34279  
    99#include "ppTranslateVersion.h"
    1010
    11 static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32);
     11static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, int extension, char *outColName, char *inColName, bool convertTo32);
    1212static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName);
    1313
    1414bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
    1515{
    16 
    17     psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    18     if (!fits) {
    19         psError(PS_ERR_IO, false, "Unable to open output file.");
    20         return false;
    21     }
    22 
    23     psMetadata *header = psMetadataAlloc(); // Header to write
    24     psString source = ppTranslateSource(), version = ppTranslateVersion();
    25     psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
    26     psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
    27     ppTranslateVersionHeader(header);
    28     psFree(source);
    29     psFree(version);
    30 
    31     psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
    32     psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
    33     psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
    34     psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
    35     psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
    36     psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
    37     psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
    38     psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
    39 
    40     // Get these header words from the first input non null input
    41     ppMopsDetections *det = NULL;
    42     for (int d = 0; d < detections->n && det == NULL; d++) {
    43         det = detections->data[d];
    44     }
    45     if (det != NULL) {
    46         psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
    47         psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
    48         psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
    49         psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
    50         psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
    51         psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
    52         psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
    53         psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
    54         psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
    55         psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
    56     } else {
    57         psWarning("no inputs with surviving detections. output header will be incomplete");
    58     }
    59     psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    60     psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
    61     psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
    62     psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
    63 
    64     //field in header that tells about the CMF version
    65     char cmfVersion[8];
    66     sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
    67     psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
    68 
    69     // Find the total number of detections
    70 
    71     long total = 0;
     16  psFits *fits = psFitsOpen(args->output, "w"); // FITS file
     17  if (!fits) {
     18    psError(PS_ERR_IO, false, "Unable to open output file.");
     19    return false;
     20  }
     21
     22  psMetadata *header = psMetadataAlloc(); // Header to write
     23  psString source = ppTranslateSource(), version = ppTranslateVersion();
     24  psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
     25  psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
     26  ppTranslateVersionHeader(header);
     27  psFree(source);
     28  psFree(version);
     29
     30  psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
     31  psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
     32  psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
     33  psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
     34  psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
     35  psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
     36  psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
     37  psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
     38
     39  // Get these header words from the first input non null input
     40  ppMopsDetections *det = NULL;
     41  for (int d = 0; d < detections->n && det == NULL; d++) {
     42    det = detections->data[d];
     43  }
     44  if (det != NULL) {
     45    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
     46    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
     47    psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
     48    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
     49    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
     50    psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
     51    psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
     52    psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
     53    psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     54    psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
     55  } else {
     56    psWarning("no inputs with surviving detections. output header will be incomplete");
     57  }
     58  psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
     59  psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
     60  psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
     61  psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
     62
     63  //field in header that tells about the CMF version
     64  char cmfVersion[8];
     65  sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
     66  psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
     67
     68  // Find the total number of detections
     69
     70  long total = 0;
     71  for (long i=0; i<detections->n; i++) {
     72    ppMopsDetections *det = detections->data[i];
     73    if (!det) {
     74      continue;
     75    }
     76    total += det->num;
     77  }
     78
     79  psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
     80
     81  if (total == 0) {
     82    // Write dummy table
     83    psMetadata *row = psMetadataAlloc(); // Output row
     84    psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
     85    psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
     86    psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
     87    psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
     88    psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
     89    psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
     90    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
     91    psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
     92    psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
     93    psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
     94    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
     95    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
     96    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
     97    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
     98    psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
     99    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
     100    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
     101    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
     102    psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
     103    psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
     104
     105    psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
     106    psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
     107    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
     108    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
     109    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
     110
     111    if (args->version == 2) {
     112      // Write data of version 2 (see ICD)
     113      psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
     114                     NAN);
     115      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
     116                     NAN);
     117      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
     118                     NAN);
     119      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
     120                     NAN);
     121      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
     122                     NAN);
     123      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
     124                     NAN);
     125      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
     126                     NAN);
     127      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
     128                     NAN);
     129      psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
     130                     NAN);
     131      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
     132                     NAN);
     133      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
     134                     NAN);
     135      psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
     136                     NAN);
     137      psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
     138                     NAN);
     139      psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
     140                     NAN);
     141      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
     142                     NAN);
     143      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
     144                     NAN);
     145      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
     146                     NAN);
     147      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
     148                     NAN);
     149      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
     150                     NAN);
     151      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
     152                     NAN);
     153      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
     154                     NAN);
     155      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
     156                     NAN);
     157      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
     158                     NAN);
     159      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
     160                     NAN);
     161      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
     162                     0);
     163      psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
     164                     0);
     165      psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
     166                     0);
     167      if (args->version == 3) {
     168        // TODO: Write data of version 3 (see ICD)
     169      }
     170    }
     171
     172    if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
     173      psErrorStackPrint(stderr, "Unable to write empty table.");
     174      psFree(header);
     175      psFree(row);
     176      return false;
     177    }
     178    psFree(row);
     179  } else {
     180
     181#define addColumn(_outName, _inName, _convertTo32)                      \
     182    if (!addOutputColumn(table, detections, total, 0, _outName, _inName, _convertTo32)) { \
     183      psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
     184      return false;                                                     \
     185    }
     186
     187    // Allocate the output table
     188    psMetadata *table = psMetadataAlloc();
     189    addColumn("RA", "RA_PSF", 0);
     190    addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     191    addColumn("DEC", "DEC_PSF", 0);
     192    addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     193    addColumn("MAG", "PSF_INST_MAG", 0);
     194    addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
     195    addColumn("PSF_CHI2", "PSF_CHISQ", 0);
     196    addColumn("PSF_DOF", "PSF_NDOF", 1);
     197    addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
     198    addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
     199    addColumn("PSF_MAJOR", NULL, 0);
     200    addColumn("PSF_MINOR", NULL, 0);     
     201    addColumn("PSF_THETA", NULL, 0);   
     202    addColumn("PSF_QUALITY", "PSF_QF", 0);
     203    addColumn("PSF_NPIX", NULL, 1);       
     204    addColumn("MOMENTS_XX", NULL, 0);
     205    addColumn("MOMENTS_XY", NULL, 0);
     206    addColumn("MOMENTS_YY", NULL, 0);
     207    addColumn("N_POS", "DIFF_NPOS", 1);
     208    addColumn("F_POS", "DIFF_FRATIO", 0);
     209    addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
     210    addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
     211    addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
     212    addColumn("FLAGS", "FLAGS", 1);
     213    addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
     214    if (args->version >= 2) {
     215      addColumn("IPP_IDET", NULL, 1);
     216      addColumn("PSF_INST_FLUX", NULL, 0);
     217      addColumn("PSF_INST_FLUX_SIG", NULL, 0);
     218      addColumn("AP_MAG", NULL, 0);
     219      addColumn("AP_MAG_RAW", NULL, 0);
     220      addColumn("AP_MAG_RADIUS", NULL, 0);
     221      addColumn("AP_FLUX", NULL, 0);
     222      addColumn("AP_FLUX_SIG", NULL, 0);
     223      addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
     224      addColumn("CAL_PSF_MAG", NULL, 0);
     225      addColumn("CAL_PSF_MAG_SIG", NULL, 0);
     226      addColumn("SKY", NULL, 0);
     227      addColumn("SKY_SIGMA", NULL, 0);
     228      addColumn("PSF_QF_PERFECT", NULL, 0);
     229      addColumn("MOMENTS_R1", NULL, 0);
     230      addColumn("MOMENTS_RH", NULL, 0);
     231      addColumn("KRON_FLUX", NULL, 0);
     232      addColumn("KRON_FLUX_ERR", NULL, 0);
     233      addColumn("KRON_FLUX_INNER", NULL, 0);
     234      addColumn("KRON_FLUX_OUTER", NULL, 0);
     235      addColumn("DIFF_R_P", NULL, 0);
     236      addColumn("DIFF_SN_P", NULL, 0);
     237      addColumn("DIFF_R_M", NULL, 0);
     238      addColumn("DIFF_SN_M", NULL, 0);
     239      addColumn("FLAGS2", NULL, 1);
     240      addColumn("IPP_IDET", NULL, 0);
     241      addColumn("N_FRAMES", NULL, 0);
     242      addColumn("PADDING", NULL, 0);
     243      if (det->fittedTrails != NULL) {
     244#define addFittedTrailColumn(_outName, _inName, _convertTo32)                   \
     245        if (!addOutputColumn(table, detections, total, 1, _outName, _inName, _convertTo32)) { \
     246          psTrace("ppMops.write", 1, "Failed to add column %s -> Replaced with NAN", _outName); \
     247        }
     248        addFittedTrailColumn("X_EXT", NULL, 0);
     249        addFittedTrailColumn("Y_EXT", NULL, 0);
     250        addFittedTrailColumn("X_EXT_SIG", NULL, 0);
     251        addFittedTrailColumn("Y_EXT_SIG", NULL, 0);
     252        addFittedTrailColumn("EXT_INST_MAG", NULL, 0);
     253        addFittedTrailColumn("EXT_INST_MAG_SIG", NULL, 0);
     254        addFittedTrailColumn("NPARAMS", NULL, 0);
     255        addFittedTrailColumn("EXT_WIDTH_MAJ", NULL, 0);
     256        addFittedTrailColumn("EXT_WIDTH_MIN", NULL, 0);
     257        addFittedTrailColumn("EXT_THETA", NULL, 0);
     258        addFittedTrailColumn("EXT_WIDTH_MAJ_ERR", NULL, 0);
     259        addFittedTrailColumn("EXT_WIDTH_MIN_ERR", NULL, 0);
     260        addFittedTrailColumn("EXT_THETA_ERR", NULL, 0);
     261      }
     262    }
     263    if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
     264      psError(psErrorCodeLast(), false, "Unable to write table");
     265      return false;
     266    }
     267    psFree(table);
     268  }
     269
     270  psFree(header);
     271  psFitsClose(fits);
     272
     273  psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
     274
     275  return true;
     276}
     277
     278// extension parameter values:
     279// 0: SkyChip.psf
     280// 1: SkyChip.xfit
     281// Any other value is ignored
     282static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, int extension, char *outColumnName, char *inColumnName, bool convertTo32)
     283{
     284  if (inColumnName == NULL) {
     285    inColumnName = outColumnName;
     286  }
     287  psVector *out = NULL;
     288  if (convertTo32) {
     289    // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
     290    // We want to write 32 bits to the output.
     291    int next = 0;
    72292    for (long i=0; i<detections->n; i++) {
    73         ppMopsDetections *det = detections->data[i];
    74         if (!det) {
    75             continue;
     293      ppMopsDetections *det = detections->data[i];
     294      if (!det || det->num == 0) {
     295        // no detections survived for this input
     296        continue;
     297      }
     298      psVector *in = NULL;
     299      if (extension == 0) {
     300        in = psMetadataLookupVector(NULL, det->table, inColumnName);
     301      } else if (extension == 1) {
     302        in = psMetadataLookupVector(NULL, det->fittedTrails, inColumnName);
     303      } else {
     304        //Error?
     305      }
     306      if (!in) {
     307        psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is true)", inColumnName);
     308        return false;
     309      }
     310      if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
     311        psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
     312                inColumnName, in->type.type);
     313        return false;
     314      }
     315      if (out == NULL) {
     316        // First time through set up the output vector and the copy parameters
     317        if (in->type.type == PS_TYPE_S64) {
     318          out = psVectorAlloc(outputSize, PS_TYPE_S32);
     319        } else {
     320          out = psVectorAlloc(outputSize, PS_TYPE_U32);
    76321        }
    77         total += det->num;
    78     }
    79 
    80     psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
    81 
    82     if (total == 0) {
    83         // Write dummy table
    84         psMetadata *row = psMetadataAlloc(); // Output row
    85         psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
    86         psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
    87         psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
    88         psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
    89         psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
    90         psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
    91         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
    92         psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
    93         psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
    94         psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
    95         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
    96         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
    97         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
    98         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
    99         psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
    100         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
    101         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
    102         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
    103         psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
    104         psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
    105 
    106         psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
    107         psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
    108         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
    109         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
    110         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
    111 
    112         if (args->version == 2) {
    113               // Write data of version 2 (see ICD)
    114               psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
    115                              NAN);
    116               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
    117                              NAN);
    118               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
    119                              NAN);
    120               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
    121                              NAN);
    122               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
    123                              NAN);
    124               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
    125                              NAN);
    126               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
    127                              NAN);
    128               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
    129                              NAN);
    130               psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
    131                              NAN);
    132               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
    133                              NAN);
    134               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
    135                              NAN);
    136               psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
    137                              NAN);
    138               psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
    139                              NAN);
    140               psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
    141                              NAN);
    142               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
    143                              NAN);
    144               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
    145                              NAN);
    146               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
    147                              NAN);
    148               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
    149                              NAN);
    150               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
    151                              NAN);
    152               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
    153                              NAN);
    154               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
    155                              NAN);
    156               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
    157                              NAN);
    158               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
    159                              NAN);
    160               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
    161                              NAN);
    162               psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
    163                              0);
    164               psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
    165                              0);
    166               psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
    167                              0);
     322      }
     323      for (long d=0; d < det->num; d++) {
     324        if (in->type.type == PS_TYPE_S64) {
     325          out->data.S32[next++] = in->data.S64[d];
     326        } else {
     327          out->data.U32[next++] = in->data.U64[d];
    168328        }
    169 
    170         if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
    171             psErrorStackPrint(stderr, "Unable to write empty table.");
    172             psFree(header);
    173             psFree(row);
    174             return false;
    175         }
    176         psFree(row);
    177     } else {
    178 
    179 #define addColumn(_outName, _inName, _convertTo32) \
    180         if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
    181             psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
    182             return false; \
    183         }
    184 
    185         // Allocate the output table
    186         psMetadata *table = psMetadataAlloc();
    187         addColumn("RA", "RA_PSF", 0);
    188         addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
    189         addColumn("DEC", "DEC_PSF", 0);
    190         addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
    191         addColumn("MAG", "PSF_INST_MAG", 0);
    192         addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
    193         addColumn("PSF_CHI2", "PSF_CHISQ", 0);
    194         addColumn("PSF_DOF", "PSF_NDOF", 1);
    195         addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
    196         addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
    197         addColumn("PSF_MAJOR", NULL, 0);
    198         addColumn("PSF_MINOR", NULL, 0);     
    199         addColumn("PSF_THETA", NULL, 0);   
    200         addColumn("PSF_QUALITY", "PSF_QF", 0);
    201         addColumn("PSF_NPIX", NULL, 1);       
    202         addColumn("MOMENTS_XX", NULL, 0);
    203         addColumn("MOMENTS_XY", NULL, 0);
    204         addColumn("MOMENTS_YY", NULL, 0);
    205         addColumn("N_POS", "DIFF_NPOS", 1);
    206         addColumn("F_POS", "DIFF_FRATIO", 0);
    207         addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
    208         addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
    209         addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
    210         addColumn("FLAGS", "FLAGS", 1);
    211         addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
    212         if (args->version == 2) {
    213             addColumn("IPP_IDET", NULL, 1);
    214             addColumn("PSF_INST_FLUX", NULL, 0);
    215             addColumn("PSF_INST_FLUX_SIG", NULL, 0);
    216             addColumn("AP_MAG", NULL, 0);
    217             addColumn("AP_MAG_RAW", NULL, 0);
    218             addColumn("AP_MAG_RADIUS", NULL, 0);
    219             addColumn("AP_FLUX", NULL, 0);
    220             addColumn("AP_FLUX_SIG", NULL, 0);
    221             addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
    222             addColumn("CAL_PSF_MAG", NULL, 0);
    223             addColumn("CAL_PSF_MAG_SIG", NULL, 0);
    224             addColumn("SKY", NULL, 0);
    225             addColumn("SKY_SIGMA", NULL, 0);
    226             addColumn("PSF_QF_PERFECT", NULL, 0);
    227             addColumn("MOMENTS_R1", NULL, 0);
    228             addColumn("MOMENTS_RH", NULL, 0);
    229             addColumn("KRON_FLUX", NULL, 0);
    230             addColumn("KRON_FLUX_ERR", NULL, 0);
    231             addColumn("KRON_FLUX_INNER", NULL, 0);
    232             addColumn("KRON_FLUX_OUTER", NULL, 0);
    233             addColumn("DIFF_R_P", NULL, 0);
    234             addColumn("DIFF_SN_P", NULL, 0);
    235             addColumn("DIFF_R_M", NULL, 0);
    236             addColumn("DIFF_SN_M", NULL, 0);
    237             addColumn("FLAGS2", NULL, 1);
    238             addColumn("IPP_IDET", NULL, 0);
    239             addColumn("N_FRAMES", NULL, 0);
    240             addColumn("PADDING", NULL, 0);
    241         }
    242         if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
    243             psError(psErrorCodeLast(), false, "Unable to write table");
    244             return false;
    245         }
    246         psFree(table);
    247     }
    248 
    249     psFree(header);
    250     psFitsClose(fits);
    251 
    252     psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
    253 
    254     return true;
     329      }
     330    }
     331  } else {
     332    void *next = NULL;
     333    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
     334    for (long i=0; i<detections->n; i++) {
     335      ppMopsDetections *det = detections->data[i];
     336      if (!det || det->num == 0) {
     337        // no detections survived for this input
     338        continue;
     339      }
     340      psVector *in = NULL;
     341      if (extension == 0) {
     342        psTrace("ppMops.write", 1, "Using det->table for %s", inColumnName);
     343        in = psMetadataLookupVector(NULL, det->table, inColumnName);
     344      } else if (extension == 1) {
     345        psTrace("ppMops.write", 1, "Using det->fittedTrails for %s", inColumnName);
     346        in = psMetadataLookupVector(NULL, det->fittedTrails, inColumnName);
     347      } else {
     348        //Error?
     349      }
     350      if (!in) {
     351        psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is false)", inColumnName);
     352        out = psVectorAlloc(outputSize, PS_TYPE_F32);
     353        psVectorInit(out, NAN);
     354        psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     355        return false;
     356      }
     357      if (out == NULL) {
     358        // First time through set up the output vector and the copy parameters
     359        out = psVectorAlloc(outputSize, in->type.type);
     360        next = (void *) out->data.U8;
     361        switch (in->type.type) {
     362        case PS_TYPE_S8:
     363          elementSize = sizeof(psS8);
     364          break;
     365        case PS_TYPE_U8:
     366          elementSize = sizeof(psU8);
     367          break;
     368        case PS_TYPE_S16:
     369          elementSize = sizeof(psS16);
     370          break;
     371        case PS_TYPE_U16:
     372          elementSize = sizeof(psU16);
     373          break;
     374        case PS_TYPE_S32:
     375          elementSize = sizeof(psS32);
     376          break;
     377        case PS_TYPE_U32:
     378          elementSize = sizeof(psU32);
     379          break;
     380        case PS_TYPE_S64:
     381          elementSize = sizeof(psS64);
     382          break;
     383        case PS_TYPE_U64:
     384          elementSize = sizeof(psU64);
     385          break;
     386        case PS_TYPE_F32:
     387          elementSize = sizeof(psF32);
     388          break;
     389        case PS_TYPE_F64:
     390          elementSize = sizeof(psF64);
     391          break;
     392        default:
     393          psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
     394          return false;
     395        }
     396      }
     397      // We are doing nasty things here so we can use memcpy.
     398      // It would be safer to do a proper loop over the elements.
     399      long toCopy = det->num * elementSize;
     400      memcpy(next, in->data.U8, toCopy);
     401      next += toCopy;
     402    }
     403  }
     404  // Finally add the new column to the output table
     405  psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     406  psFree(out);    // drop reference
     407  return true;
    255408}
    256409
    257 static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
    258 {
    259     if (inColumnName == NULL) {
    260         inColumnName = outColumnName;
    261     }
    262 
    263     psVector *out = NULL;
    264     if (convertTo32) {
    265         // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
    266         // We want to write 32 bits to the output.
    267         int next = 0;
    268         for (long i=0; i<detections->n; i++) {
    269             ppMopsDetections *det = detections->data[i];
    270             if (!det || det->num == 0) {
    271                 // no detections survived for this input
    272                 continue;
    273             }
    274             psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
    275             if (!in) {
    276                 psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
    277                 return false;
    278             }
    279             if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
    280                 psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
    281                     inColumnName, in->type.type);
    282                 return false;
    283             }
    284             if (out == NULL) {
    285                 // First time through set up the output vector and the copy parameters
    286                 if (in->type.type == PS_TYPE_S64) {
    287                     out = psVectorAlloc(outputSize, PS_TYPE_S32);
    288                 } else {
    289                     out = psVectorAlloc(outputSize, PS_TYPE_U32);
    290                 }
    291             }
    292             for (long d=0; d < det->num; d++) {
    293                 if (in->type.type == PS_TYPE_S64) {
    294                     out->data.S32[next++] = in->data.S64[d];
    295                 } else {
    296                     out->data.U32[next++] = in->data.U64[d];
    297                 }
    298             }
    299         }
    300     } else {
    301         void *next = NULL;
    302         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
    303         for (long i=0; i<detections->n; i++) {
    304             ppMopsDetections *det = detections->data[i];
    305             if (!det || det->num == 0) {
    306                 // no detections survived for this input
    307                 continue;
    308             }
    309             psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
    310             if (!in) {
    311                 psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
    312                 return false;
    313             }
    314             if (out == NULL) {
    315                 // First time through set up the output vector and the copy parameters
    316                 out = psVectorAlloc(outputSize, in->type.type);
    317                 next = (void *) out->data.U8;
    318                 switch (in->type.type) {
    319                     case PS_TYPE_S8:
    320                     elementSize = sizeof(psS8);
    321                     break;
    322                 case PS_TYPE_U8:
    323                     elementSize = sizeof(psU8);
    324                     break;
    325                 case PS_TYPE_S16:
    326                     elementSize = sizeof(psS16);
    327                     break;
    328                 case PS_TYPE_U16:
    329                     elementSize = sizeof(psU16);
    330                     break;
    331                 case PS_TYPE_S32:
    332                     elementSize = sizeof(psS32);
    333                     break;
    334                 case PS_TYPE_U32:
    335                     elementSize = sizeof(psU32);
    336                     break;
    337                 case PS_TYPE_S64:
    338                     elementSize = sizeof(psS64);
    339                     break;
    340                 case PS_TYPE_U64:
    341                     elementSize = sizeof(psU64);
    342                     break;
    343                 case PS_TYPE_F32:
    344                     elementSize = sizeof(psF32);
    345                     break;
    346                 case PS_TYPE_F64:
    347                     elementSize = sizeof(psF64);
    348                     break;
    349                 default:
    350                     psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
    351                     return false;
    352 
    353                 }
    354             }
    355             // We are doing nasty things here so we can use memcpy.
    356             // It would be safer to do a proper loop over the elements.
    357             long toCopy = det->num * elementSize;
    358             memcpy(next, in->data.U8, toCopy);
    359             next += toCopy;
    360         }
    361     }
    362 
    363     // Finally add the new column to the output table
    364     psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
    365     psFree(out);    // drop reference
    366 
    367     return true;
    368 }
    369410static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
    370411{
    371     psVector *out = psVectorAlloc(total, PS_TYPE_S64);
    372     long next = 0;
    373     for (long i = 0; i<detections->n; i++) {
    374         ppMopsDetections *det = detections->data[i];
    375         if (!det) {
    376             continue;
    377         }
    378         psS64 diffSkyfileId = det->diffSkyfileId;
    379         for (long j = 0; j < det->num; j++) {
    380             out->data.S64[next++] = diffSkyfileId;
    381         }
    382     }
    383     psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
    384     psFree(out);
    385     return true;
     412  psVector *out = psVectorAlloc(total, PS_TYPE_S64);
     413  long next = 0;
     414  for (long i = 0; i<detections->n; i++) {
     415    ppMopsDetections *det = detections->data[i];
     416    if (!det) {
     417      continue;
     418    }
     419    psS64 diffSkyfileId = det->diffSkyfileId;
     420    for (long j = 0; j < det->num; j++) {
     421      out->data.S64[next++] = diffSkyfileId;
     422    }
     423  }
     424  psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
     425  psFree(out);
     426  return true;
    386427}
Note: See TracChangeset for help on using the changeset viewer.