IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 5, 2012, 5:19:48 PM (14 years ago)
Author:
mhuber
Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:
branches/meh_branches/ppstack_test
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/meh_branches/ppstack_test

  • branches/meh_branches/ppstack_test/ppTranslate/src/ppMopsRead.c

    r30656 r33415  
    2020    psArray *detections = psArrayAlloc(num); // Array of detections, to return
    2121    for (int i = 0; i < num; i++) {
    22         psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file
    23 
     22        const char *name = inNames->data[i];
     23
     24        psFits *fits = psFitsOpen(name,  "r"); // FITS file
    2425        if (!fits) {
    2526            psError(PS_ERR_IO, false, "Unable to open input %d", i);
    2627            return false;
    2728        }
     29
    2830        psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
    2931        if (!header) {
    … …  
    7476        }
    7577        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     78        detections->data[i] = det;
     79        det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
     80        det->num = size;
     81        det->diffSkyfileId = diffSkyfileId;
    7682
    7783        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    … …  
    9096                                     psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    9197
    92         int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    93         int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     98        det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     99        det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    94100
    95101        psFree(header);
    96102
    97         psArray *table = psFitsReadTable(fits); // Table of interest
     103        psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
    98104        if (!table) {
    99105            psError(PS_ERR_IO, false, "Unable to read table %d", i);
    100             return false;
    101         }
     106            return NULL;
     107        }
     108        det->table = table;
    102109        psFitsClose(fits);
     110        if (args->version == 0) {
     111          if (skyChipPsfVersion < 2) {
     112            // XXX: TODO: Do we need to add dummy vectors for the missing columns?
     113           }
     114        }
     115
     116        psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
     117        psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     118
     119        det->raErr = psVectorAlloc(size, PS_TYPE_F64);
     120        det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     121        det->mask = psVectorAlloc(size, PS_TYPE_U8);
     122
     123        // convert ra and dec to radians for use in the purge duplicates function
     124        det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     125        det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     126        det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
     127        det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
     128        if (!det->ra || !det->dec || !det->x || !det->y) {
     129            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
     130            return NULL;
     131        }
     132
     133        // Add our new vectors to the table so that duplicates and masked items may be purged
     134        psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
     135        psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     136
     137        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
     138
     139        psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
     140        psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
     141        psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
     142        psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
     143        psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
     144        psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
     145        psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
    103146
    104147        double plateScale = 0.0;        // Plate scale
    105148        long numGood = 0;               // Number of good rows
    106         for (long j = 0; j < size; j++) {
    107             psMetadata *row = table->data[j]; // Row of interest
    108 
    109             psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
     149        for (long row = 0; row < size; row++) {
     150
     151            psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    110152            if (flags & SOURCE_MASK) {
    111                 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
     153                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     154                det->mask->data.U8[row] = 0xFF;
    112155                continue;
    113156            }
    114157
    115             det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
    116             det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
    117             det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
    118             det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
    119             det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
    120             det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
    121             det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");
    122             det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");
    123             det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");
    124             det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
    125             det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");
    126             det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");
    127             det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");
    128             det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");
    129             det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");
    130             det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");
    131             det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");
    132             det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");
    133             det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
    134             det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
    135             det->naxis1->data.S32[numGood] = naxis1;
    136             det->naxis2->data.S32[numGood] = naxis2;
    137             det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");
    138             det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");
    139             det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");
    140             det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");
    141             det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");
    142 
    143             //Additions of 2010-10-25
    144             if (args->version == 2) {
    145               //Values are set only if the version is 2
    146               if (skyChipPsfVersion == 2) {
    147                 det->psfInstFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX");
    148                 det->psfInstFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX_SIG");
    149                 det->apMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG");
    150                 det->apMagRaw->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RAW");
    151                 det->apMagRadius->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RADIUS");
    152                 det->apFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX");
    153                 det->apFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX_SIG");
    154                 det->peakFluxAsMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PEAK_FLUX_AS_MAG");
    155                 det->calPsfMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG");
    156                 det->calPsfMagSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG_SIG");
    157                 det->sky->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY");
    158                 det->skySig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY_SIGMA");
    159                 det->qualityPerfect->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF_PERFECT");
    160                 det->momentsR1->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_R1");
    161                 det->momentsRH->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_RH");
    162                 det->kronFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX");
    163                 det->kronFluxErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_ERR");
    164                 det->kronFluxInner->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_INNER");
    165                 det->kronFluxOuter->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_OUTER");
    166                 det->diffRP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_P");
    167                 det->diffSnP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_P");
    168                 det->diffRM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_M");
    169                 det->diffSnM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_M");
    170                 det->flags2->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS2");
    171                 det->ippIdet->data.U32[numGood] = psMetadataLookupU32(NULL, row, "IPP_IDET");
    172                 det->nFrames->data.U16[numGood] = psMetadataLookupU16(NULL, row, "N_FRAMES");
    173                 det->padding->data.S16[numGood] = psMetadataLookupS16(NULL, row, "PADDING");
    174               } else {
    175                 det->psfInstFlux->data.F32[numGood] = NAN;
    176                 det->psfInstFluxSig->data.F32[numGood] = NAN;
    177                 det->apMag->data.F32[numGood] = NAN;
    178                 det->apMagRaw->data.F32[numGood] = NAN;
    179                 det->apMagRadius->data.F32[numGood] = NAN;
    180                 det->apFlux->data.F32[numGood] = NAN;
    181                 det->apFluxSig->data.F32[numGood] = NAN;
    182                 det->peakFluxAsMag->data.F32[numGood] = NAN;
    183                 det->calPsfMag->data.F32[numGood] = NAN;
    184                 det->calPsfMagSig->data.F32[numGood] = NAN;
    185                 det->sky->data.F32[numGood] = NAN;
    186                 det->skySig->data.F32[numGood] = NAN;
    187                 det->qualityPerfect->data.F32[numGood] = NAN;
    188                 det->momentsR1->data.F32[numGood] = NAN;
    189                 det->momentsRH->data.F32[numGood] = NAN;
    190                 det->kronFlux->data.F32[numGood] = NAN;
    191                 det->kronFluxErr->data.F32[numGood] = NAN;
    192                 det->kronFluxInner->data.F32[numGood] = NAN;
    193                 det->kronFluxOuter->data.F32[numGood] = NAN;
    194                 det->diffRP->data.F32[numGood] = NAN;
    195                 det->diffSnP->data.F32[numGood] = NAN;
    196                 det->diffRM->data.F32[numGood] = NAN;
    197                 det->diffSnM->data.F32[numGood] = NAN;
    198                 det->flags2->data.U32[numGood] = 0;
    199                 det->ippIdet->data.U32[numGood] = 0;
    200                 det->nFrames->data.U16[numGood] = 0;
    201                 det->padding->data.S16[numGood] = 0;
    202               }
    203             }
    204 
    205158            // Calculate error in RA, Dec
    206             double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32
    207             double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32
    208             double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); //SC: Warning! Promotion of F32
    209             double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); //SC: Warning! Promotion of F32
    210 
    211             if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
    212                 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
    213                 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
     159
     160            double xErr = xErrV->data.F32[row];
     161            double yErr = yErrV->data.F32[row];
     162            double scale = scaleV->data.F32[row];
     163            double angle = angleV->data.F32[row];
     164
     165            if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
     166                !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
     167                !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
    214168                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    215169                psTrace("ppMops.read", 10,
    216170                        "Discarding row %ld from input %d because of non-finite values: "
    217171                        "%f %f %lf %lf %f %f %f %f %f %f",
    218                         j, i,
    219                         det->x->data.F32[numGood], det->y->data.F32[numGood],
    220                         det->ra->data.F64[numGood], det->dec->data.F64[numGood],
    221                         det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
     172                        row, i,
     173                        det->x->data.F32[row], det->y->data.F32[row],
     174                        det->ra->data.F64[row], det->dec->data.F64[row],
     175                        mag->data.F32[row], magErr->data.F32[row],
    222176                        xErr, yErr, scale, angle);
     177                det->mask->data.U8[row] = 0xFF;
    223178                continue;
    224179            }
    … …  
    229184            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    230185            double errScale = scale / 3600.0;
    231             det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    232             det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    233 
    234             det->mask->data.U8[numGood] = 0;
     186            det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     187            det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     188
     189            det->mask->data.U8[row] = 0;
    235190            plateScale += scale;
    236191            numGood++;
    … …  
    238193        det->seeing *= ((float) plateScale) / ((float) numGood);
    239194
    240         det->x->n = numGood;
    241         det->y->n = numGood;
    242         det->ra->n = numGood;
    243         det->dec->n = numGood;
    244         det->raErr->n = numGood;
    245         det->decErr->n = numGood;
    246         det->mag->n = numGood;
    247         det->magErr->n = numGood;
    248         det->chi2->n = numGood;
    249         det->dof->n = numGood;
    250         det->cr->n = numGood;
    251         det->extended->n = numGood;
    252         det->psfMajor->n = numGood;
    253         det->psfMinor->n = numGood;
    254         det->psfTheta->n = numGood;
    255         det->quality->n = numGood;
    256         det->numPix->n = numGood;
    257         det->xxMoment->n = numGood;
    258         det->xyMoment->n = numGood;
    259         det->yyMoment->n = numGood;
    260         det->flags->n = numGood;
    261         det->diffSkyfileId->n = numGood;
    262         det->naxis1->n = numGood;
    263         det->naxis2->n = numGood;
    264         det->mask->n = numGood;
    265         det->nPos->n = numGood;
    266         det->fPos->n = numGood;
    267         det->nRatioBad->n = numGood;
    268         det->nRatioMask->n = numGood;
    269         det->nRatioAll->n = numGood;
    270         det->psfInstFlux->n = numGood;
    271         det->psfInstFluxSig->n = numGood;
    272         det->apMag->n = numGood;
    273         det->apMagRaw->n = numGood;
    274         det->apMagRadius->n = numGood;
    275         det->apFlux->n = numGood;
    276         det->apFluxSig->n = numGood;
    277         det->peakFluxAsMag->n = numGood;
    278         det->calPsfMag->n = numGood;
    279         det->calPsfMagSig->n = numGood;
    280         det->sky->n = numGood;
    281         det->skySig->n = numGood;
    282         det->qualityPerfect->n = numGood;
    283         det->momentsR1->n = numGood;
    284         det->momentsRH->n = numGood;
    285         det->kronFlux->n = numGood;
    286         det->kronFluxErr->n = numGood;
    287         det->kronFluxInner->n = numGood;
    288         det->kronFluxOuter->n = numGood;
    289         det->diffRP->n = numGood;
    290         det->diffSnP->n = numGood;
    291         det->diffRM->n = numGood;
    292         det->diffSnM->n = numGood;
    293         det->flags2->n = numGood;
    294         det->ippIdet->n = numGood;
    295         det->nFrames->n = numGood;
    296         det->padding->n = numGood;
    297 
    298         det->num = numGood;
     195        // Are we using numGood for anything outside of this function?
     196        det->numGood = numGood;
    299197
    300198        if (isfinite(args->zp) && numGood > 0) {
    301             psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    302         }
    303 
    304         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    305 
    306         psFree(table);
    307         detections->data[i] = det;
     199            psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
     200        }
     201
     202        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
    308203    }
    309204
Note: See TracChangeset for help on using the changeset viewer.