Changeset 32971 for trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

Timestamp:

Dec 19, 2011, 1:22:04 PM (14 years ago)

Author:

bills

Message:

Add and call functions to read the extended source extensions for the PS1_SV1
cmf format

File:

• trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (modified) (3 diffs)

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -539 +539 @@
 }
 
+bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n");
+        return false;
+    }
+
+    // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point
+    float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP");
+    float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME");
+    float magOffset = zeropt + 2.5*log10(exptime);
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        pmSource *source = NULL;
+#ifndef ASSUME_SORTED
+        long j = seq < sources->n ? seq : sources->n - 1;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+#else
+        long j = sourceIndex[seq];
+        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
+        source = sources->data[j];
+#endif
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+
+        if (!source->extpars) {
+            source->extpars = pmSourceExtendedParsAlloc ();
+        }
+        pmSourceExtendedPars *extpars = source->extpars;
+
+        // Assume that X_EXT Y_EXT and sigmas match the psf src so skip
+
+        // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by
+        // psf size or something?
+        extpars->axes.major = 1.0;
+        extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO");
+        extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
+
+        float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
+        float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
+        if (isfinite(mag)) {
+            extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
+            if (isfinite(magErr)) {
+                extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
+            }
+        }
+
+        extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
+        extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR");
+        extpars->petrosianR50      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50");
+        extpars->petrosianR50Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR");
+        extpars->petrosianR90      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90");
+        extpars->petrosianR90Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR");
+        extpars->petrosianFill     = psMetadataLookupF32(&status, row, "PETRO_FILL");
+
+        psVector *radSB   = psMetadataLookupVector(&status, row, "PROF_SB");
+        psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX");
+        psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL");
+
+        if (radSB && radSB->n > 0) {
+            extpars->radProfile = pmSourceRadialProfileAlloc();
+            extpars->radProfile->binSB   = psMemIncrRefCounter(radSB);
+            extpars->radProfile->binSum   = psMemIncrRefCounter(radFlux);
+            extpars->radProfile->binFill = psMemIncrRefCounter(radFill);
+        }
+
+        psFree(row);
+    }
+
+    return true;
+}
+
 // XXX this layout is still the same as PS1_DEV_1
 bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
@@ -684 +778 @@
     psFree (outhead);
     psFree (table);
+    return true;
+}
+
+bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n");
+        return false;
+    }
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq.
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        long j = seq < sources->n ? seq : sources->n - 1;
+        pmSource *source = NULL;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+        if (!source->modelFits) {
+            // XXX: where to find the number of models to expect?
+            source->modelFits = psArrayAllocEmpty(5);
+        }
+        psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE");
+        if (!modelName) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i);
+            psFree(row);
+            return false;
+        }
+        pmModelType modelType = pmModelClassGetType(modelName);
+        if (modelType < 0) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName);
+            psFree(row);
+            return false;
+        }
+        pmModel *model = pmModelAlloc(modelType);
+
+        psF32 *PAR = model->params->data.F32;
+        psF32 *dPAR = model->dparams->data.F32;
+
+        PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT");
+        PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT");
+        dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG");
+        dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG");
+
+        model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG");
+        model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
+
+        psEllipseAxes axes;
+        axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
+        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
+        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
+        if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
+            // Do we need to fail here or can this happen?
+            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
+            psFree(model);
+            psFree(row);
+            return false;
+        }
+        // XXX: clean this up
+        if (model->params->n > 7) {
+            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
+        }
+        // read the covariance matrix
+        int nparams = model->params->n;
+        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
+        for (int y = 0; y < nparams; y++) {
+            for (int x = 0; x < nparams; x++) {
+                char name[64];
+                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
+                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+            }
+        }
+        model->covar = covar;
+
+        psArrayAdd(source->modelFits, 1, model);
+        psFree(model);
+
+        psFree(row);
+    }
+
     return true;
 }
@@ -831 +1025 @@
     return true;
 }
+
+bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n");
+        return false;
+    }
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq.
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        long j = seq < sources->n ? seq : sources->n - 1;
+        pmSource *source = NULL;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+        if (!source->radialAper) {
+            // XXX: where to find the number of models to expect?
+            source->radialAper = psArrayAllocEmpty(5);
+        }
+        pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc();
+
+        radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX"));
+        radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV"));
+        radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR"));
+        radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL"));
+
+        psArrayAdd(source->radialAper, 1, radialAper);
+
+        psFree(radialAper);
+        psFree(row);
+    }
+
+    return true;
+}