Changeset 9822

Timestamp:

Nov 1, 2006, 12:39:53 PM (20 years ago)

Author:

Paul Price

Message:

Adding functions to read/write an array of fringe statistics to a single table. This allows us to use multiple readouts as the different fringe components, and then add a table containing the fringe measurements.

Location:

trunk/psModules/src/detrend

Files:

• pmFringeStats.c (modified) (2 diffs)
• pmFringeStats.h (modified) (2 diffs)

trunk/psModules/src/detrend/pmFringeStats.c

@@ -122 +122 @@
     // Vectors: x, y, mask
 
-    psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header
+    psMetadata *scalars = psMemIncrRefCounter(header); // Metadata to hold the scalars; will be the header
+    if (!scalars) {
+        scalars = psMetadataAlloc();
+    }
     psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
                      regions->dX);
@@ -503 +506 @@
 }
 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// pmFringeIO
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool pmFringesWriteFits(psFits *fits, psMetadata *header, const psArray *fringes, const char *extname)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_ARRAY_NON_NULL(fringes, false);
+
+    // Check the regions are all identical
+    pmFringeRegions *regions = ((pmFringeStats*)fringes->data[0])->regions; // First region
+    for (int i = i; i < fringes->n; i++) {
+        pmFringeStats *stats = fringes->data[i];
+        if (stats->regions != regions) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Regions for fringe statistics are not identical.\n");
+            return false;
+        }
+    }
+
+    // Ensure the region is legit
+    psVector *x = regions->x;           // The x positions
+    psVector *y = regions->y;           // The y positions
+    psVector *mask = regions->mask;     // The region mask
+    int numRows = regions->nRequested;  // Number of rows in the table
+    PS_ASSERT_INT_POSITIVE(numRows, false);
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_SIZE(x, (long)numRows, false);
+    PS_ASSERT_VECTOR_SIZE(y, (long)numRows, false);
+    if (mask) {
+        PS_ASSERT_VECTOR_NON_NULL(mask, false);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
+        PS_ASSERT_VECTOR_SIZE(mask, (long)numRows, false);
+    }
+
+    // We need to write:
+    // Scalars: dX, dY, nX, nY
+    // Vectors: x, y, mask, f, df
+
+    psMetadata *scalars = psMemIncrRefCounter(header); // Metadata to hold the scalars; will be the header
+    if (!scalars) {
+        scalars = psMetadataAlloc();
+    }
+    psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
+                     regions->dX);
+    psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height",
+                     regions->dY);
+    psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x",
+                     regions->nX);
+    psMetadataAddS32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y",
+                     regions->nY);
+
+    psArray *table = psArrayAlloc(numRows); // The table
+    // Translate the vectors into the required format for psFitsWriteTable()
+    for (long i = 0; i < numRows; i++) {
+        psMetadata *row = psMetadataAlloc();
+        psMetadataAddF32(row, PS_LIST_TAIL, "x", PS_META_REPLACE, "Fringe position in x", x->data.F32[i]);
+        psMetadataAddF32(row, PS_LIST_TAIL, "y", PS_META_REPLACE, "Fringe position in y", y->data.F32[i]);
+        psU8 maskValue = 0;             // Mask value
+        if (mask && mask->data.U8[i]) {
+            maskValue = 0xff;
+        }
+
+        psVector *f = psVectorAlloc(fringes->n, PS_TYPE_F32); // Measurements for each fringe component
+        psVector *df = psVectorAlloc(fringes->n, PS_TYPE_F32); // Errors in measurements
+        for (long j = 0; j < fringes->n; j++) {
+            pmFringeStats *stats = fringes->data[j]; // Fringe statistics of interest
+            f->data.F32[j] = stats->f->data.F32[i];
+            df->data.F32[j] = stats->df->data.F32[i];
+            if (!isfinite(f->data.F32[j]) || !isfinite(df->data.F32[j])) {
+                maskValue = 0xff;
+            }
+        }
+        psMetadataAdd(row, PS_LIST_TAIL, "f", PS_DATA_VECTOR | PS_META_REPLACE, "Fringe measurements", f);
+        psMetadataAdd(row, PS_LIST_TAIL, "df", PS_DATA_VECTOR | PS_META_REPLACE, "Fringe errors", df);
+        // Drop references
+        psFree(f);
+        psFree(df);
+
+        psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", maskValue);
+        table->data[i] = row;
+    }
+
+    bool success;                       // Success of operation
+    if (!(success = psFitsWriteTable(fits, scalars, table, extname))) {
+        psError(PS_ERR_IO, false, "Unable to write fringe data to extension %s\n", extname);
+    }
+    psFree(scalars);
+    psFree(table);
+
+    return success;
+}
+
+
+psArray *pmFringesReadFits(psMetadata *header, psFits *fits, const char *extname)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, NULL);
+
+    if (extname && strlen(extname) > 0 && !psFitsMoveExtName(fits, extname)) {
+        psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
+        return NULL;
+    }
+
+    psMetadata *headerCopy = psMemIncrRefCounter(header); // Copy of the header, or NULL
+    headerCopy = psFitsReadHeader(headerCopy, fits); // The FITS header
+    if (!headerCopy) {
+        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname);
+        psFree(headerCopy);
+        return NULL;
+    }
+
+    // Read the scalars from the header
+    #define READ_FRINGES_SCALAR(SCALAR, NAME) \
+    int SCALAR = psMetadataLookupS32(&mdok, headerCopy, NAME); \
+    if (!mdok || SCALAR <= 0) { \
+        psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s\n", extname); \
+        psFree(headerCopy); \
+        return NULL; \
+    }
+
+    // Need to retrieve the scalars: dX, dY, nX, nY
+    bool mdok = true;                   // Status of MD lookup
+    READ_SCALAR(dX, "PSFRNGDX");
+    READ_SCALAR(dY, "PSFRNGDY");
+    READ_SCALAR(nX, "PSFRNGNX");
+    READ_SCALAR(nY, "PSFRNGNY");
+    psFree(headerCopy);
+
+    // Now the vectors: x, y, mask, f, df
+    psArray *table = psFitsReadTable(fits); // The table
+    long numRows = table->n;            // Number of rows
+
+    pmFringeRegions *regions = pmFringeRegionsAlloc(numRows, dX, dY, nX, nY); // The fringe regions
+    psVector *x = psVectorAlloc(numRows, PS_TYPE_F32); // x position
+    psVector *y = psVectorAlloc(numRows, PS_TYPE_F32); // y position
+    psVector *mask = psVectorAlloc(numRows, PS_TYPE_U8); // mask
+    regions->x = x;
+    regions->y = y;
+    regions->mask = mask;
+    psArray *f = psArrayAlloc(numRows); // Array of fringe measurements
+    psArray *df = psArrayAlloc(numRows);// Array of errors
+    psArray *fringes = NULL; // Array of fringes, to return
+
+    #define READ_FRINGES_VECTOR_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \
+    VECTOR->data.TYPE[i] = psMetadataLookup##TYPE(&mdok, row, NAME); \
+    if (!mdok) { \
+        psError(PS_ERR_IO, true, "Unable to find " #DESCRIPTION " .\n"); \
+        goto READ_FRINGES_DONE; \
+    }
+
+    #define READ_FRINGES_ARRAY_ROW(ARRAY, NAME, DESCRIPTION) \
+    ARRAY->data[i] = psMetadataLookupPtr(&mdok, row, NAME); \
+    if (!ARRAY->data[i] || !mdok) { \
+        psError(PS_ERR_IO, true, "Unable to find " #DESCRIPTION " .\n"); \
+        goto READ_FRINGES_DONE; \
+    }
+
+    // Translate the table into vectors
+    for (long i = 0; i < numRows; i++) {
+        psMetadata *row = table->data[i]; // Table row
+        READ_FRINGES_VECTOR_ROW(x, F32, "x", "x position");
+        READ_FRINGES_VECTOR_ROW(y, F32, "y", "y position");
+        READ_FRINGES_VECTOR_ROW(mask, U8, "mask", "mask");
+        READ_FRINGES_ARRAY_ROW(f, "f", "fringe measurement");
+        READ_FRINGES_ARRAY_ROW(df, "df", "fringe error");
+    }
+
+    // Get f,df into pmFringeStats
+    long numFringes = ((psVector*)(f->data[0]))->n; // Number of fringe components
+    fringes = psArrayAlloc(numFringes);
+    for (int j = 0; j < numFringes; j++) {
+        fringes->data[j] = pmFringeStatsAlloc(regions);
+    }
+
+    for (long i = 0; i < numRows; i++) {
+        psVector *measurements = f->data[i]; // Vector of measurements
+        psVector *errors = df->data[i]; // Vector of errors
+        for (int j = 0; j < numFringes; j++) {
+            pmFringeStats *fringe = fringes->data[j];
+            fringe->f->data.F32[i] = measurements->data.F32[j];
+            fringe->df->data.F32[i] = errors->data.F32[j];
+        }
+    }
+
+READ_FRINGES_DONE:
+    psFree(table);
+    psFree(regions);
+    psFree(x);
+    psFree(y);
+    psFree(mask);
+    psFree(f);
+    psFree(df);
+
+    return fringes;
+}
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////

trunk/psModules/src/detrend/pmFringeStats.h

@@ -8 +8 @@
 /// @author Paul Price, IfA
 ///
-/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2006-10-17 22:03:32 $
+/// @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2006-11-01 22:39:53 $
 ///
 /// Copyright 2004-2006 Institute for Astronomy, University of Hawaii
@@ -142 +142 @@
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Input/output for multiple pmFringeStats
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/// Write an array of fringes measurements to a FITS table.
+///
+/// Writes an array of fringe measurements for the same FPA component as a FITS table.  The array of fringe
+/// statistics must all use the same fringe regions (or there is no point in storing them all together).  The
+/// header is supplemented with scalar values dX, dY, nX and nY (as PSFRNGDX, PSFRNGDY, PSFRNGNX, PSFRNGNY)
+/// from the fringe regions, while the fringe coordinates and mask are written as a FITS table (as x, y,
+/// mask, f, df; f and df are vectors).
+bool pmFringesWriteFits(psFits *fits,    ///< FITS file to which to write
+                        psMetadata *header, ///< Header, or NULL
+                        const psArray *fringes, ///< Array of pmFringeStats, all for the same pmFringeRegion
+                        const char *extname ///< Extension name for table
+                       );
+
+/// Reads an array of fringes measurements from a FITS table.
+///
+/// The fringes are read from the FITS file, at the given extension name.  The table provides the region and
+/// the (possibly multiple) fringe statistics for that region.  The current extension is used if the extension
+/// name is not provided.
+psArray *pmFringesReadFits(psMetadata *header, ///< Header to read, or NULL
+                           psFits *fits,///< FITS file from which to read
+                           const char *extname ///< Extension name to read, or NULL
+                          );
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // pmFringeScale
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////