Changeset 21363 for trunk/psModules/src/camera/pmFPARead.c

Timestamp:

Feb 5, 2009, 4:31:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:

• trunk/psModules/src/camera/pmFPARead.c (modified) (20 diffs)

trunk/psModules/src/camera/pmFPARead.c

@@ -9 +9 @@
 #include <pslib.h>
 
+#include "pmConfig.h"
+#include "pmConfigMask.h"
 #include "pmHDU.h"
 #include "pmFPA.h"
@@ -28 +30 @@
     FPA_READ_TYPE_IMAGE,                // Read image
     FPA_READ_TYPE_MASK,                 // Read mask
-    FPA_READ_TYPE_WEIGHT,               // Read weight map
+    FPA_READ_TYPE_VARIANCE,             // Read variance map
     FPA_READ_TYPE_HEADER                // Read header
 } fpaReadType;
@@ -54 +56 @@
       case FPA_READ_TYPE_MASK:
         return readout->thisMaskScan;
-      case FPA_READ_TYPE_WEIGHT:
-        return readout->thisWeightScan;
+      case FPA_READ_TYPE_VARIANCE:
+        return readout->thisVarianceScan;
       default:
         psAbort("Unknown read type: %x\n", type);
@@ -74 +76 @@
         readout->thisMaskScan = thisScan;
         return readout->lastMaskScan;
-      case FPA_READ_TYPE_WEIGHT:
-        readout->thisWeightScan = thisScan;
-        return readout->lastWeightScan;
+      case FPA_READ_TYPE_VARIANCE:
+        readout->thisVarianceScan = thisScan;
+        return readout->lastVarianceScan;
       default:
         psAbort("Unknown read type: %x\n", type);
@@ -93 +95 @@
       case FPA_READ_TYPE_MASK:
         return readout->lastMaskScan;
-      case FPA_READ_TYPE_WEIGHT:
-        return readout->lastWeightScan;
+      case FPA_READ_TYPE_VARIANCE:
+        return readout->lastVarianceScan;
       default:
         psAbort("Unknown read type: %x\n", type);
@@ -113 +115 @@
         readout->lastMaskScan = lastScan;
         return readout->lastMaskScan;
-      case FPA_READ_TYPE_WEIGHT:
-        readout->lastWeightScan = lastScan;
-        return readout->lastWeightScan;
+      case FPA_READ_TYPE_VARIANCE:
+        readout->lastVarianceScan = lastScan;
+        return readout->lastVarianceScan;
       default:
         psAbort("Unknown read type: %x\n", type);
@@ -132 +134 @@
       case FPA_READ_TYPE_MASK:
         return &readout->mask;
-      case FPA_READ_TYPE_WEIGHT:
-        return &readout->weight;
+      case FPA_READ_TYPE_VARIANCE:
+        return &readout->variance;
       default:
         psAbort("Unknown read type: %x\n", type);
@@ -193 +195 @@
 
     return naxis3;
+}
+
+// Determine whether a FITS file contains covariance matrices
+static bool hduCovariance(pmHDU *hdu,   // Header data unit
+                          psFits *fits  // FITS file
+    )
+{
+    if (hdu->extname && !psFitsMoveExtName(fits, hdu->extname)) {
+        psError(PS_ERR_IO, false, "Unable to move to extension %s", hdu->extname);
+        return false;
+    }
+    // Need to explicitly read the header, since the HDU may not contain the variance header
+    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
+    if (!header) {
+        psError(PS_ERR_IO, false, "Unable to read variance header.");
+        return false;
+    }
+    bool mdok;                          // Status of MD lookup
+    bool covar = psMetadataLookupBool(&mdok, header, PM_HDU_COVARIANCE_KEYWORD); // Got covariance?
+    psFree(header);
+    return covar;
 }
 
@@ -350 +373 @@
     *target = psImageSubset(image, region);
 
-    // Get the list of overscans: only for IMAGE types (no overscan for MASK and WEIGHT)
+    // Get the list of overscans: only for IMAGE types (no overscan for MASK and VARIANCE)
     if (type == FPA_READ_TYPE_IMAGE) {
         if (readout->bias->n != 0) {
@@ -515 +538 @@
     }
 
-    // XXX for IMAGE, we need the CELL.BAD value, but for MASK, we need the BAD mask value
-
-    float bad = 0;
-    if (type == FPA_READ_TYPE_MASK) {
-      bad = 1.0;
-    } else {
-      bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
+    // Need to set the invalid (unread) pixels appropriately, and to read the mask bits
+    float bad = 0;                      // Bad level
+    switch (type) {
+      case FPA_READ_TYPE_MASK: {
+          // Need to explicitly read the header, since what's in the pmHDU may not correspond to the mask
+          psMetadata *header = psFitsReadHeader(NULL, fits);
+          if (!header) {
+              psError(PS_ERR_IO, false, "Unable to read mask header.");
+              return false;
+          }
+          if (!pmConfigMaskReadHeader(config, header)) {
+              psError(PS_ERR_IO, false, "Unable to determine mask bits");
+              psFree(header);
+              return false;
+          }
+          psFree(header);
+          bad = pmConfigMaskGet("BAD", config);
+          break;
+      }
+      case FPA_READ_TYPE_IMAGE:
+      case FPA_READ_TYPE_VARIANCE:
+        bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD");
+        break;
+      default:
+        psAbort("Unrecognised type: %x", type);
     }
 
@@ -582 +623 @@
     *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]);
 
-    // Read overscans only for "image" type --- weights and masks shouldn't record overscans
+    // Read overscans only for "image" type --- variances and masks shouldn't record overscans
     if (type == FPA_READ_TYPE_IMAGE) {
         // Blow away existing data
@@ -608 +649 @@
 }
 
-// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadWeight
+// Read into an cell; this is the engine for pmCellRead, pmCellReadMask, pmCellReadVariance
 // Does most of the work for the reading --- reads the HDU, and portions the HDU into readouts.
 static bool cellRead(pmCell *cell,      // Cell into which to read
@@ -640 +681 @@
         dataPointer = hdu->masks;
         break;
-      case FPA_READ_TYPE_WEIGHT:
-        hduReadFunc = pmHDUReadWeight;
-        dataPointer = hdu->weights;
+      case FPA_READ_TYPE_VARIANCE:
+        hduReadFunc = pmHDUReadVariance;
+        dataPointer = hdu->variances;
         break;
       default:
@@ -680 +721 @@
         imageArray = hdu->masks;
         break;
-      case FPA_READ_TYPE_WEIGHT:
-        imageArray = hdu->weights;
+      case FPA_READ_TYPE_VARIANCE:
+        imageArray = hdu->variances;
         break;
       default:
@@ -729 +770 @@
 
 
-// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadWeight
+// Read into an chip; this is the engine for pmChipRead, pmChipReadMask, pmChipReadVariance
 // Iterates over component cells, reading each
 static bool chipRead(pmChip *chip,      // Chip into which to read
@@ -760 +801 @@
 
 
-// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadWeight
+// Read into an FPA; this is the engine for pmFPARead, pmFPAReadMask, pmFPAReadVariance
 // Iterates over component chips, reading each
 static bool fpaRead(pmFPA *fpa,         // FPA into which to read
@@ -852 +893 @@
     float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.BAD is not set --- assuming zero.\n");
+        psWarning("CELL.BAD is not set --- assuming zero.\n");
         bad = 0.0;
     }
@@ -1091 +1132 @@
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Reading the weight map
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-bool pmReadoutMoreWeight(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
+// Reading the variance map
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool pmReadoutMoreVariance(pmReadout *readout, psFits *fits, int z, int numScans, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_FITS_NON_NULL(fits, false);
 
-    return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_WEIGHT, config);
-}
-
-bool pmReadoutReadChunkWeight(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
+    return readoutMore(readout, fits, z, numScans, FPA_READ_TYPE_VARIANCE, config);
+}
+
+bool pmReadoutReadChunkVariance(pmReadout *readout, psFits *fits, int z, int numScans, int overlap,
                               pmConfig *config)
 {
@@ -1110 +1151 @@
     PS_ASSERT_INT_NONNEGATIVE(numScans, false);
 
-    return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_WEIGHT, config);
-}
-
-bool pmReadoutReadWeight(pmReadout *readout, psFits *fits, int z, pmConfig *config)
+    return readoutReadChunk(readout, fits, z, numScans, overlap, FPA_READ_TYPE_VARIANCE, config);
+}
+
+bool pmReadoutReadVariance(pmReadout *readout, psFits *fits, int z, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_FITS_NON_NULL(fits, false);
 
-    return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_WEIGHT, config);
-}
-
-bool pmCellReadWeight(pmCell *cell, psFits *fits, pmConfig *config)
+    return readoutReadChunk(readout, fits, z, 0, 0, FPA_READ_TYPE_VARIANCE, config);
+}
+
+bool pmCellReadVariance(pmCell *cell, psFits *fits, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(cell, false);
     PS_ASSERT_FITS_NON_NULL(fits, false);
 
-    return cellRead(cell, fits, config, FPA_READ_TYPE_WEIGHT);
-}
-
-bool pmChipReadWeight(pmChip *chip, psFits *fits, pmConfig *config)
+    if (!cellRead(cell, fits, config, FPA_READ_TYPE_VARIANCE)) {
+        return false;
+    }
+    pmHDU *hdu = pmHDUFromCell(cell);   // Header data unit
+    if (hduCovariance(hdu, fits)) {
+        return pmCellReadCovariance(cell, fits);
+    }
+    return true;
+}
+
+bool pmChipReadVariance(pmChip *chip, psFits *fits, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(chip, false);
     PS_ASSERT_FITS_NON_NULL(fits, false);
 
-    return chipRead(chip, fits, config, FPA_READ_TYPE_WEIGHT);
-}
-
-bool pmFPAReadWeight(pmFPA *fpa, psFits *fits, pmConfig *config)
+    if (!chipRead(chip, fits, config, FPA_READ_TYPE_VARIANCE)) {
+        return false;
+    }
+    pmHDU *hdu = pmHDUFromChip(chip);   // Header data unit
+    if (hduCovariance(hdu, fits)) {
+        return pmChipReadCovariance(chip, fits);
+    }
+    return true;
+}
+
+bool pmFPAReadVariance(pmFPA *fpa, psFits *fits, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(fpa, false);
     PS_ASSERT_FITS_NON_NULL(fits, false);
 
-    return fpaRead(fpa, fits, config, FPA_READ_TYPE_WEIGHT);
+    if (!fpaRead(fpa, fits, config, FPA_READ_TYPE_VARIANCE)) {
+        return false;
+    }
+    pmHDU *hdu = pmHDUFromFPA(fpa);     // Header data unit
+    if (hduCovariance(hdu, fits)) {
+        return pmFPAReadCovariance(fpa, fits);
+    }
+    return true;
 }
 
@@ -1268 +1330 @@
     return numRead;
 }
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Reading covariance matrices
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool pmCellReadCovariance(pmCell *cell, psFits *fits)
+{
+    PS_ASSERT_PTR_NON_NULL(cell, false);
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
+    psString extname = NULL;            // Extension name
+    psStringAppend(&extname, "COVAR_%s_%s", chipName, cellName);
+
+    if (!psFitsMoveExtName(fits, extname)) {
+        psError(PS_ERR_IO, false, "Unable to move to extension %s\n", extname);
+        psFree(extname);
+        return false;
+    }
+    psFree(extname);
+
+    psMetadata *header = psFitsReadHeader(NULL, fits); // The FITS header
+    if (!header) {
+        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", extname);
+        psFree(header);
+        return false;
+    }
+
+    bool mdok;                          // Status of MD lookup
+    int x0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.X"); // Centre of matrix in x
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre");
+        psFree(header);
+        return false;
+    }
+    int y0 = psMetadataLookupS32(&mdok, header, "COVARIANCE.CENTRE.Y"); // Centre of matrix in y
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to read covariance centre");
+        psFree(header);
+        return false;
+    }
+    psFree(header);
+
+    psArray *images = psFitsReadImageCube(fits, psRegionSet(0,0,0,0)); // Covariance matrices
+    if (!images) {
+        psError(PS_ERR_IO, false, "Unable to read covariance matrices for chip %s, cell %s",
+                chipName, cellName);
+        return false;
+    }
+
+    psArray *readouts = cell->readouts; // Readouts of cell
+    if (images->n != readouts->n) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Number of covariance matrices (%ld) doesn't match number of readouts (%ld)",
+                images->n, readouts->n);
+        psFree(images);
+        return false;
+    }
+
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *ro = readouts->data[i]; // Readout of interest
+        psImage *image = images->data[i]; // Image of interest
+        if (ro->covariance) {
+            psWarning("Clobbering extant covariance matrix in chip %s, cell %s, readout %d",
+                      chipName, cellName, i);
+            psFree(ro->covariance);
+        }
+        ro->covariance = psKernelAllocFromImage(image, x0, y0);
+    }
+    psFree(images);
+
+    return true;
+}
+
+
+bool pmChipReadCovariance(pmChip *chip, psFits *fits)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    psArray *cells = chip->cells;       // Array of cells
+    for (int i = 0; i < cells->n; i++) {
+        pmCell *cell = cells->data[i];  // Cell of interest
+        if (!pmCellReadCovariance(cell, fits)) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+bool pmFPAReadCovariance(pmFPA *fpa, psFits *fits)
+{
+    PS_ASSERT_PTR_NON_NULL(fpa, false);
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    psArray *chips = fpa->chips;        // Array of chips
+    for (int i = 0; i < chips->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!pmChipReadCovariance(chip, fits)) {
+            return false;
+        }
+    }
+
+    return true;
+}