Changeset 7875 for trunk/psModules/src/detrend/pmFringeStats.c

Timestamp:

Jul 11, 2006, 6:44:01 PM (20 years ago)

Author:

Paul Price

Message:

• Allow user to specify a random number generator for pmFringeRegionsCreatePoints.
• Split pmFringeStatsWriteFits into pmFringeRegionsWriteFits and pmFringeStatsWriteFits.
• Split pmFringeStatsReadFits into pmFringeRegionsReadFits and pmFringeStatsReadFits.

File:

• trunk/psModules/src/detrend/pmFringeStats.c (modified) (12 diffs)

trunk/psModules/src/detrend/pmFringeStats.c

@@ -3 +3 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-07-07 03:26:22 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-07-12 04:44:01 $
  *
  *  Copyright 2004 IfA
@@ -56 +56 @@
 }
 
-bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image)
+bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, const psImage *image, psRandom *random)
 {
     PS_ASSERT_PTR_NON_NULL(fringe, false);
@@ -65 +65 @@
     // create fringe->nRequested
 
-    psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    psRandom *rng;
+    if (random) {
+        rng = psMemIncrRefCounter(random);
+    } else {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
 
     fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32);
@@ -84 +89 @@
     // generate random points located within image bounds
     for (int i = 0; i < fringe->nRequested; i++) {
-        frnd = psRandomUniform(rnd);
+        frnd = psRandomUniform(rng);
         xPt[i] = (nX - 2*dX)* frnd + dX;
-        frnd = psRandomUniform(rnd);
+        frnd = psRandomUniform(rng);
         yPt[i] = (nY - 2*dY)* frnd + dY;
     }
+
+    psFree(rng);
+
     return true;
 }
 
+bool pmFringeRegionsWriteFits(psFits *fits, psMetadata *header,
+                              const pmFringeRegions *regions, const char *extname)
+{
+    // Make sure the input is well-behaved
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(regions, false);
+    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, numRows, false);
+    PS_ASSERT_VECTOR_SIZE(y, 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, numRows, false);
+    }
+
+    // We need to write:
+    // Scalars: dX, dY, nX, nY
+    // Vectors: x, y
+
+    psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header
+    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
+                     regions->dX);
+    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height",
+                     regions->dY);
+    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x",
+                     regions->nX);
+    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y",
+                     regions->nY);
+
+    psArray *table = psArrayAlloc(numRows); // The table
+    table->n = numRows;
+    // 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]);
+        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;
+}
+
+pmFringeRegions *pmFringeRegionsReadFits(psMetadata *header, const psFits *fits, const char *extname)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, NULL);
+
+    if (extname && strlen(extname) > 0) {
+        if (!psFitsMoveExtName(fits, extname)) {
+            psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
+            return NULL;
+        }
+    } else if (!psFitsMoveExtNum(fits, 0, false)) {
+        psError(PS_ERR_IO, false, "Unable to move to PHU\n");
+        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(header);
+        return NULL;
+    }
+
+    // Read the scalars from the header
+    #define READ_SCALAR(SCALAR, NAME) \
+    int SCALAR = psMetadataLookupS32(&mdok, header, NAME); \
+    if (!mdok || SCALAR <= 0) { \
+        psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s.\n", extname); \
+        psFree(header); \
+        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(header);
+
+    // Now the vectors: x, y
+    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
+    x->n = y->n = numRows;
+    regions->x = x;
+    regions->y = y;
+
+    #define READ_REGIONS_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"); \
+        psFree(table); \
+        psFree(regions); \
+        return NULL; \
+    }
+
+    // Translate the table into vectors
+    for (long i = 0; i < numRows; i++) {
+        psMetadata *row = table->data[i]; // Table row
+        READ_REGIONS_ROW(x, F32, "x", "x position");
+        READ_REGIONS_ROW(y, F32, "y", "y position");
+    }
+    psFree(table);
+
+    return regions;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -129 +264 @@
     if (!fringe->x || !fringe->y) {
         // create the fringe vectors for this image
-        pmFringeRegionsCreatePoints(fringe, readout->image);
+        pmFringeRegionsCreatePoints(fringe, readout->image, NULL);
     }
 
@@ -167 +302 @@
             psImageStats(median, subImage, subMask, maskVal);
             sky->data.F32[i][j] = median->sampleMedian;
+            psFree(subImage);
+            psFree(subMask);
         }
     }
@@ -193 +330 @@
                 (int)region.y0, (int)region.y1, fPt[i], dfPt[i]);
     }
+    psFree(sky);
+    psFree(median);
+    psFree(medianSd);
 
     return measurements;
 }
 
-bool pmFringeStatsWriteFits(psFits *fits, const pmFringeStats *fringe, const char *extname)
+bool pmFringeStatsWriteFits(psFits *fits,
+                            psMetadata *header,
+                            const pmFringeStats *fringe,
+                            const char *extname
+                           )
 {
     // Make sure the input is well-behaved
@@ -204 +348 @@
     pmFringeRegions *regions = fringe->regions; // The fringe regions
     PS_ASSERT_PTR_NON_NULL(regions, false);
-    psVector *x = regions->x;           // The x positions
-    psVector *y = regions->y;           // The y positions
-    psVector *mask = regions->mask;     // The region mask
-    int numRows = x->n;                 // Number of rows in the table
-    PS_ASSERT_INT_POSITIVE(regions->nRequested, 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(y, 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, numRows, false);
-    }
+    int numRows = regions->nRequested;  // Number of rows in the table
+    PS_ASSERT_INT_POSITIVE(numRows, false);
     psVector *f = fringe->f;            // The fringe measurements
     psVector *df = fringe->df;      // The fringe standard deviatiations
@@ -229 +360 @@
 
     // We need to write:
-    // Scalars: dX, dY, nX, nY
-    // Vectors: x, y, mask, f, df
-
-    psMetadata *scalars = psMetadataAlloc(); // Metadata to hold the scalars; will be the header
-    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDX", PS_META_REPLACE, "Median box half-width",
-                     regions->dX);
-    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGDY", PS_META_REPLACE, "Median box half-height",
-                     regions->dY);
-    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNX", PS_META_REPLACE, "Large-scale smoothing in x",
-                     regions->nX);
-    psMetadataAddF32(scalars, PS_LIST_TAIL, "PSFRNGNY", PS_META_REPLACE, "Large-scale smoothing in y",
-                     regions->nY);
-
+    // Vectors: f, df
     psArray *table = psArrayAlloc(numRows); // The table
     table->n = numRows;
@@ -247 +366 @@
     for (long i = 0; i < numRows; i++) {
         psMetadata *row = psMetadataAlloc();
-        psArraySet(table, i, row);
-        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]);
-        if (mask) {
-            psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", mask->data.F32[i]);
-        } else {
-            psMetadataAddU8(row, PS_LIST_TAIL, "mask", PS_META_REPLACE, "Mask", 0);
-        }
         psMetadataAddF32(row, PS_LIST_TAIL, "f", PS_META_REPLACE, "Fringe measurement", f->data.F32[i]);
         psMetadataAddF32(row, PS_LIST_TAIL, "df", PS_META_REPLACE, "Fringe stdev", df->data.F32[i]);
-        psFree(row);                    // Drop reference
-    }
-
-    bool success;                       // Success of operation
-    if (!(success = psFitsWriteTable(fits, scalars, table, extname))) {
+        table->data[i] = row;
+    }
+
+    if (!psFitsWriteTable(fits, header, table, extname)) {
         psError(PS_ERR_IO, false, "Unable to write fringe data to extension %s\n", extname);
-    }
-    psFree(scalars);
+        psFree(table);
+        return false;
+    }
+
     psFree(table);
-
-    return success;
-}
-
-pmFringeStats *pmFringeStatsReadFits(const psFits *fits, const char *extname)
+    return true;
+}
+
+pmFringeStats *pmFringeStatsReadFits(psMetadata *header, const psFits *fits, const char *extname,
+                                     pmFringeRegions *regions)
 {
     PS_ASSERT_PTR_NON_NULL(fits, NULL);
+    PS_ASSERT_PTR_NON_NULL(regions, NULL);
+    PS_ASSERT_INT_POSITIVE(regions->nRequested, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(regions->x, regions->y, NULL);
+    PS_ASSERT_VECTOR_SIZE(regions->x, regions->nRequested, NULL);
 
     if (extname && strlen(extname) > 0) {
@@ -284 +400 @@
     }
 
-    psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header
-    if (!header) {
+    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_SCALAR(SCALAR, NAME) \
-    int SCALAR = psMetadataLookupS32(&mdok, header, NAME); \
-    if (!mdok || SCALAR <= 0) { \
-        psError(PS_ERR_IO, true, "Unable to find " NAME " in header of extension %s.\n", extname); \
-        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, "PSFRNGDX");
-    READ_SCALAR(nX, "PSFRNGDX");
-    READ_SCALAR(nY, "PSFRNGDX");
-    psFree(header);
-
-    // Now the vectors: x, y, mask, f, df
+    psFree(headerCopy);
+
+    // Now the vectors: 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
-    x->n = y->n = mask->n = numRows;
-    regions->x = x;
-    regions->y = y;
-    regions->mask = mask;
     pmFringeStats *fringes = pmFringeStatsAlloc(regions); // The fringe measurements
-    psFree(regions);                    // Drop reference
     psVector *f = fringes->f;           // fringe measurement
-    psVector *df = fringes->df;     // fringe stdev
-
-    #define READ_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \
+    psVector *df = fringes->df;         // fringe stdev
+
+    #define READ_STATS_ROW(VECTOR, TYPE, NAME, DESCRIPTION) \
     VECTOR->data.TYPE[i] = psMetadataLookup##TYPE(&mdok, row, NAME); \
     if (!mdok) { \
@@ -333 +428 @@
 
     // Translate the table into vectors
+    bool mdok;                          // Status of MD lookup
     for (long i = 0; i < numRows; i++) {
         psMetadata *row = table->data[i]; // Table row
-        READ_ROW(x, F32, "x", "x position");
-        READ_ROW(y, F32, "y", "y position");
-        READ_ROW(mask, U8, "mask", "mask value");
-        READ_ROW(f, F32, "f", "fringe measurement");
-        READ_ROW(df, F32, "df", "fringe standard deviation");
+        READ_STATS_ROW(f, F32, "f", "fringe measurement");
+        READ_STATS_ROW(df, F32, "df", "fringe standard deviation");
     }
     psFree(table);