IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 15, 2011, 4:50:36 PM (15 years ago)
Author:
bills
Message:

Modify ppMops to use psFitsReadTableAllColumns and psFitsWriteTableAllColumns. These
functions store the columns in psVectors. This avoids the memory explosion that
results by storing the table as an array of metadatas representing rows.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsWrite.c

    r29567 r32406  
    99#include "ppTranslateVersion.h"
    1010
    11 bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args)
     11static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32);
     12static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName);
     13
     14bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
    1215{
    13     psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output);
    1416
    1517    psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    … …  
    3638    psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
    3739
     40    // Get these header words from the first input
     41    ppMopsDetections *det = detections->data[0];
    3842    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
    3943    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
    … …  
    5559    psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
    5660
    57     if (det->num == 0) {
     61    // Find the total number of detections
     62
     63    long total = 0;
     64    for (long i=0; i<detections->n; i++) {
     65        ppMopsDetections *det = detections->data[i];
     66        total += det->num;
     67    }
     68
     69    psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
     70
     71    if (total == 0) {
    5872        // Write dummy table
    5973        psMetadata *row = psMetadataAlloc(); // Output row
    … …  
    151165        psFree(row);
    152166    } else {
    153         psArray *table = psArrayAlloc(det->num); // Table to write
    154         for (long i = 0; i < det->num; i++) {
    155             psMetadata *row = psMetadataAlloc(); // Output row
    156             psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)",
    157                              RAD_TO_DEG(det->ra->data.F64[i]));
    158             psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)",
    159                              det->raErr->data.F64[i]);
    160             psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)",
    161                              RAD_TO_DEG(det->dec->data.F64[i]));
    162             psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)",
    163                              det->decErr->data.F64[i]);
    164             psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]);
    165             psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]);
    166             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", det->chi2->data.F32[i]);
    167             psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit",
    168                              det->dof->data.S32[i]);
    169             psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR",
    170                              det->cr->data.F32[i]);
    171             psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness",
    172                              det->extended->data.F32[i]);
    173             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", det->psfMajor->data.F32[i]);
    174             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", det->psfMinor->data.F32[i]);
    175             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)",
    176                              det->psfTheta->data.F32[i]);
    177             psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor",
    178                              det->quality->data.F32[i]);
    179             psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF",
    180                              det->numPix->data.S32[i]);
    181             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", det->xxMoment->data.F32[i]);
    182             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", det->xyMoment->data.F32[i]);
    183             psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", det->yyMoment->data.F32[i]);
    184             psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels",
    185                              det->nPos->data.S32[i]);
    186             psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels",
    187                              det->fPos->data.F32[i]);
    188             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative",
    189                              det->nRatioBad->data.F32[i]);
    190             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked",
    191                              det->nRatioMask->data.F32[i]);
    192             psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all",
    193                              det->nRatioAll->data.F32[i]);
    194             psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]);
    195             psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile",
    196                              det->diffSkyfileId->data.S64[i]);
    197 
    198             if (args->version == 2) {
    199               // Write data of version 2 (see ICD)
    200               psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
    201                              det->ippIdet->data.U32[i]);
    202               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
    203                              det->psfInstFlux->data.F32[i]);
    204               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
    205                              det->psfInstFluxSig->data.F32[i]);
    206               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
    207                              det->apMag->data.F32[i]);
    208               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
    209                              det->apMagRaw->data.F32[i]);
    210               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
    211                              det->apMagRadius->data.F32[i]);
    212               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
    213                              det->apFlux->data.F32[i]);
    214               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
    215                              det->apFluxSig->data.F32[i]);
    216               psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
    217                              det->peakFluxAsMag->data.F32[i]);
    218               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
    219                              det->calPsfMag->data.F32[i]);
    220               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
    221                              det->calPsfMagSig->data.F32[i]);
    222               psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
    223                              det->sky->data.F32[i]);
    224               psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
    225                              det->skySig->data.F32[i]);
    226               psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
    227                              det->qualityPerfect->data.F32[i]);
    228               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
    229                              det->momentsR1->data.F32[i]);
    230               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
    231                              det->momentsRH->data.F32[i]);
    232               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
    233                              det->kronFlux->data.F32[i]);
    234               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
    235                              det->kronFluxErr->data.F32[i]);
    236               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
    237                              det->kronFluxInner->data.F32[i]);
    238               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
    239                              det->kronFluxOuter->data.F32[i]);
    240               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
    241                              det->diffRP->data.F32[i]);
    242               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
    243                              det->diffSnP->data.F32[i]);
    244               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
    245                              det->diffRM->data.F32[i]);
    246               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
    247                              det->diffSnM->data.F32[i]);
    248               psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
    249                              det->flags2->data.U32[i]);
    250               psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
    251                              det->nFrames->data.U16[i]);
    252               psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
    253                              det->padding->data.S16[i]);
    254             }
    255 
    256             //Update with the table with the current row
    257             table->data[i] = row;
    258         }
    259         if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) {
    260             psErrorStackPrint(stderr, "Unable to write table.");
    261             psFree(header);
    262             psFree(table);
     167
     168#define addColumn(_outName, _inName, _convertTo32) \
     169        if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
     170            psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
     171            return false; \
     172        }
     173
     174        // Allocate the output table
     175        psMetadata *table = psMetadataAlloc();
     176        addColumn("RA", "RA_PSF", 0);
     177        addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     178        addColumn("DEC", "DEC_PSF", 0);
     179        addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     180        addColumn("MAG", "PSF_INST_MAG", 0);
     181        addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
     182        addColumn("PSF_CHI2", "PSF_CHISQ", 0);
     183        addColumn("PSF_DOF", "PSF_NDOF", 1);
     184        addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
     185        addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
     186        addColumn("PSF_MAJOR", NULL, 0);
     187        addColumn("PSF_MINOR", NULL, 0);     
     188        addColumn("PSF_THETA", NULL, 0);   
     189        addColumn("PSF_QUALITY", "PSF_QF", 0);
     190        addColumn("PSF_NPIX", NULL, 1);       
     191        addColumn("MOMENTS_XX", NULL, 0);
     192        addColumn("MOMENTS_XY", NULL, 0);
     193        addColumn("MOMENTS_YY", NULL, 0);
     194        addColumn("N_POS", "DIFF_NPOS", 1);
     195        addColumn("F_POS", "DIFF_FRATIO", 0);
     196        addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
     197        addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
     198        addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
     199        addColumn("FLAGS", "FLAGS", 1);
     200        addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
     201        if (args->version == 2) {
     202            addColumn("IPP_IDET", NULL, 1);
     203            addColumn("PSF_INST_FLUX", NULL, 0);
     204            addColumn("PSF_INST_FLUX_SIG", NULL, 0);
     205            addColumn("AP_MAG", NULL, 0);
     206            addColumn("AP_MAG_RAW", NULL, 0);
     207            addColumn("AP_MAG_RADIUS", NULL, 0);
     208            addColumn("AP_FLUX", NULL, 0);
     209            addColumn("AP_FLUX_SIG", NULL, 0);
     210            addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
     211            addColumn("CAL_PSF_MAG", NULL, 0);
     212            addColumn("CAL_PSF_MAG_SIG", NULL, 0);
     213            addColumn("SKY", NULL, 0);
     214            addColumn("SKY_SIGMA", NULL, 0);
     215            addColumn("PSF_QF_PERFECT", NULL, 0);
     216            addColumn("MOMENTS_R1", NULL, 0);
     217            addColumn("MOMENTS_RH", NULL, 0);
     218            addColumn("KRON_FLUX", NULL, 0);
     219            addColumn("KRON_FLUX_ERR", NULL, 0);
     220            addColumn("KRON_FLUX_INNER", NULL, 0);
     221            addColumn("KRON_FLUX_OUTER", NULL, 0);
     222            addColumn("DIFF_R_P", NULL, 0);
     223            addColumn("DIFF_SN_P", NULL, 0);
     224            addColumn("DIFF_R_M", NULL, 0);
     225            addColumn("DIFF_SN_M", NULL, 0);
     226            addColumn("FLAGS2", NULL, 1);
     227            addColumn("IPP_IDET", NULL, 0);
     228            addColumn("N_FRAMES", NULL, 0);
     229            addColumn("PADDING", NULL, 0);
     230        }
     231        if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
     232            psError(psErrorCodeLast(), false, "Unable to write table");
    263233            return false;
    264234        }
    … …  
    273243    return true;
    274244}
     245
     246static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
     247{
     248    if (inColumnName == NULL) {
     249        inColumnName = outColumnName;
     250    }
     251
     252    psVector *out = NULL;
     253    if (convertTo32) {
     254        // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
     255        // We want to write 32 bits to the output.
     256        int next = 0;
     257        for (long i=0; i<detections->n; i++) {
     258            ppMopsDetections *det = detections->data[i];
     259            if (det->num == 0) {
     260                // no detections survived for this input
     261                continue;
     262            }
     263            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     264            if (!in) {
     265                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     266                return false;
     267            }
     268            if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
     269                psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
     270                    inColumnName, in->type.type);
     271                return false;
     272            }
     273            if (out == NULL) {
     274                // First time through set up the output vector and the copy parameters
     275                if (in->type.type == PS_TYPE_S64) {
     276                    out = psVectorAlloc(outputSize, PS_TYPE_S32);
     277                } else {
     278                    out = psVectorAlloc(outputSize, PS_TYPE_U32);
     279                }
     280            }
     281            for (long d=0; d < det->num; d++) {
     282                if (in->type.type == PS_TYPE_S64) {
     283                    out->data.S32[next++] = in->data.S64[d];
     284                } else {
     285                    out->data.U32[next++] = in->data.U64[d];
     286                }
     287            }
     288        }
     289    } else {
     290        void *next = NULL;
     291        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
     292        for (long i=0; i<detections->n; i++) {
     293            ppMopsDetections *det = detections->data[i];
     294            if (det->num == 0) {
     295                // no detections survived for this input
     296                continue;
     297            }
     298            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
     299            if (!in) {
     300                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
     301                return false;
     302            }
     303            if (out == NULL) {
     304                // First time through set up the output vector and the copy parameters
     305                out = psVectorAlloc(outputSize, in->type.type);
     306                next = (void *) out->data.U8;
     307                switch (in->type.type) {
     308                    case PS_TYPE_S8:
     309                    elementSize = sizeof(psS8);
     310                    break;
     311                case PS_TYPE_U8:
     312                    elementSize = sizeof(psU8);
     313                    break;
     314                case PS_TYPE_S16:
     315                    elementSize = sizeof(psS16);
     316                    break;
     317                case PS_TYPE_U16:
     318                    elementSize = sizeof(psU16);
     319                    break;
     320                case PS_TYPE_S32:
     321                    elementSize = sizeof(psS32);
     322                    break;
     323                case PS_TYPE_U32:
     324                    elementSize = sizeof(psU32);
     325                    break;
     326                case PS_TYPE_S64:
     327                    elementSize = sizeof(psS64);
     328                    break;
     329                case PS_TYPE_U64:
     330                    elementSize = sizeof(psU64);
     331                    break;
     332                case PS_TYPE_F32:
     333                    elementSize = sizeof(psF32);
     334                    break;
     335                case PS_TYPE_F64:
     336                    elementSize = sizeof(psF64);
     337                    break;
     338                default:
     339                    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
     340                    return false;
     341
     342                }
     343            }
     344            // We are doing nasty things here so we can use memcpy.
     345            // It would be safer to do a proper loop over the elements.
     346            long toCopy = det->num * elementSize;
     347            memcpy(next, in->data.U8, toCopy);
     348            next += toCopy;
     349        }
     350    }
     351
     352    // Finally add the new column to the output table
     353    psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     354    psFree(out);    // drop reference
     355
     356    return true;
     357}
     358static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
     359{
     360    psVector *out = psVectorAlloc(total, PS_TYPE_S64);
     361    long next = 0;
     362    for (long i = 0; i<detections->n; i++) {
     363        ppMopsDetections *det = detections->data[i];
     364        psS64 diffSkyfileId = det->diffSkyfileId;
     365        for (long j = 0; j < det->num; j++) {
     366            out->data.S64[next++] = diffSkyfileId;
     367        }
     368    }
     369    psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
     370    psFree(out);
     371    return true;
     372}
Note: See TracChangeset for help on using the changeset viewer.