Changeset 21264

Timestamp:

Feb 2, 2009, 8:05:38 AM (17 years ago)

Author:

Paul Price

Message:

Adding work on integrating covariances into the read/write system. Not sure I like this framework (putting covariances for an HDU into a common table); might rework to write a covariance image for each cell. I think I started doing that and ran into problems, but I can't remember what they were and how they stopped me.

Location:

branches/pap_branch_20090128/psModules/src/camera

Files:

• pmFPARead.c (modified) (4 diffs)
• pmFPAWrite.c (modified) (4 diffs)
• pmFPAfileIO.c (modified) (2 diffs)
• pmHDU.c (modified) (3 diffs)
• pmHDU.h (modified) (2 diffs)
• pmHDUCovariances.c (added)
• pmHDUCovariances.h (added)
• pmHDUGenerate.c (modified) (11 diffs)

branches/pap_branch_20090128/psModules/src/camera/pmFPARead.c

@@ -28 +28 @@
     FPA_READ_TYPE_IMAGE,                // Read image
     FPA_READ_TYPE_MASK,                 // Read mask
-    FPA_READ_TYPE_VARIANCE,               // Read variance map
+    FPA_READ_TYPE_VARIANCE,             // Read variance map
     FPA_READ_TYPE_HEADER                // Read header
 } fpaReadType;
@@ -515 +515 @@
     }
 
-    // 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: {
+          pmHDU *phu = pmHDUGetHighest(cell->parent->parent, cell->parent, cell); // Primary header
+          if (!pmConfigMaskReadHeader(config, phu->header)) {
+              psError(PS_ERR_IO, false, "Unable to determine mask bits");
+              return false;
+          }
+          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);
     }
 
@@ -721 +732 @@
             return NULL;
         }
+
+        if (type == FPA_READ_TYPE_VARIANCE && hdu->covariances) {
+            psArray *covariances = hdu->covariances; // Covariances in HDU
+            for (int i = 0; i < covariances->n; i++) {
+                pmHDUCovariances *covar = covariances->data[i]; // Covariance information
+
         psFree(readout);                // Drop reference
     }
@@ -852 +869 @@
     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;
     }

branches/pap_branch_20090128/psModules/src/camera/pmFPAWrite.c

@@ -29 +29 @@
     FPA_WRITE_TYPE_IMAGE,               // Write image
     FPA_WRITE_TYPE_MASK,                // Write mask
-    FPA_WRITE_TYPE_VARIANCE               // Write variance map
+    FPA_WRITE_TYPE_VARIANCE             // Write variance map
 } fpaWriteType;
 
@@ -74 +74 @@
 }
 
+
+
+
+// Some type-specific additions to the header
+static bool writeUpdateHeader(pmFPA *fpa, // FPA of interest
+                              pmChip *chip, // Chip of interest, or NULL
+                              pmCell *cell, // Cell of interest, or NULL
+                              fpaWriteType type // Type to write
+    )
+{
+    switch (type) {
+      case FPA_WRITE_TYPE_MASK: {
+          pmHDU *phu = pmHDUGetHighest(fpa, chip, cell); // Primary header
+          if (!pmConfigMaskWriteHeader(config, phu->header)) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to set the mask names in the PHU header");
+              return false;
+          }
+          break;
+      }
+      case FPA_WRITE_TYPE_VARIANCE: {
+          pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // Header being written
+
+
+
 // Write a cell image/mask/variance
 static bool cellWrite(pmCell *cell,     // Cell to write
@@ -97 +121 @@
     // generate the HDU, but only copies the structure.
     if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.\n");
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.");
         return false;
     }
@@ -106 +130 @@
 
     if (writeBlank || writeImage) {
-
         pmFPAUpdateNames(cell->parent->parent, cell->parent, cell);
         pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS |
                                  PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE;
         if (!pmConceptsWriteCell(cell, source, true, config)) {
-            psError(PS_ERR_IO, false, "Unable to write concepts for cell.\n");
+            psError(PS_ERR_IO, false, "Unable to write concepts for cell.");
             return false;
         }
+
+
         if (!appropriateWriteFunc(hdu, fits, config, type)) {
             psError(PS_ERR_IO, false, "Unable to write HDU for cell.\n");

branches/pap_branch_20090128/psModules/src/camera/pmFPAfileIO.c

@@ -771 +771 @@
           }
 
-          // XXX if we have a mask file, then we need to read the mask bit names
-          // defined for this file
-          if (file->type == PM_FPA_FILE_MASK) {
-            if (!pmConfigMaskReadHeader (config, phu)) {
-                psError(PS_ERR_IO, false, "error in mask bits");
-                return false;
-            }
-          }
-
           // determine the current format from the header
           // determine camera if not specified already
@@ -786 +777 @@
           psString formatName = NULL;
           psString cameraName = NULL;
-          file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true);
+          file->format = pmConfigCameraFormatFromHeader(&camera, &cameraName, &formatName, config, phu, true);
           if (!file->format) {
             psError(PS_ERR_IO, false, "Failed to read CCD format from %s\n", file->filename);

branches/pap_branch_20090128/psModules/src/camera/pmHDU.c

@@ -12 +12 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+
+#define COVARIANCE_INDICATOR "PS_COVAR" // FITS keyword to indicate if covariance matrices are present
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -160 +162 @@
     PS_ASSERT_PTR_NON_NULL(fits, false);
 
-    return hduRead(hdu, &hdu->variances, fits);
+    if (!hduRead(hdu, &hdu->variances, fits)) {
+        return false;
+    }
+
+    return true;
 }
 
@@ -240 +246 @@
     PS_ASSERT_PTR_NON_NULL(fits, false);
 
+    pmHDUCovariancesWriteHeader(hdu);
+
     psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask
-    return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits);
-}
+    if (!hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits)) {
+        return false;
+    }
+    pmHDUCovarianceClearHeader(hdu);
+
+    if (!pmHDUCovarianceWrite(hdu, fits)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to write covariance");
+        return false;
+    }
+
+    return true;
+}

branches/pap_branch_20090128/psModules/src/camera/pmHDU.h

@@ -4 +4 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.9.22.1 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-29 00:33:51 $
+ * @version $Revision: 1.9.22.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-02 18:05:38 $
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -22 +22 @@
 /// Of course, it is not an exact replica of a FITS HDU --- they have no mask and variance data, but these are
 /// stored here for convenience --- it keeps all the relevant data about the image in one place.
-typedef struct
-{
+typedef struct {
     psString extname;                   ///< The extension name
     bool blankPHU;                      ///< Is this a blank FITS Primary Header Unit, i.e., no data?
     psMetadata *format;                 ///< The camera format
     psMetadata *header;                 ///< The FITS header, or NULL if primary for FITS; or section info
-    psArray *images;                    ///< The pixel data
-    psArray *variances;                   ///< The pixel data
-    psArray *masks;                     ///< The pixel data
-}
-pmHDU;
+    psArray *images;                    ///< Pixel data
+    psArray *variances;                 ///< Variance in the pixel data, or NULL
+    psArray *masks;                     ///< Mask for the pixel data, or NULL
+    psArray *covariances;               ///< Covariance matrices (pmHDUCovariance), or NULL
+} pmHDU;
 
 

branches/pap_branch_20090128/psModules/src/camera/pmHDUGenerate.c

@@ -355 +355 @@
     psElemType imageType = 0;           // Type of readout images
     psElemType maskType = 0;            // Type of readout masks
-    psElemType varianceType = 0;          // Type of readout variances
+    psElemType varianceType = 0;        // Type of readout variances
+    psElemType covarianceType = 0;      // Type of readout covariances
     {
         psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells
@@ -384 +385 @@
                     varianceType = checkTypes(varianceType, readout->variance->type.type);
                 }
+                if (!hdu->covariances && readout->covariance) {
+                    covarianceType = checkTypes(covarianceType, PS_TYPE_F32);
+                }
             }
         }
         psFree(iter);
     }
-    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) {
+    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0 && covarianceType == 0)) {
         // Nothing from which to create an HDU
         psFree(cells);
@@ -428 +432 @@
             hdu->variances->data[i] = variance;
         }
+    }
+    if (covarianceType) {
+        hdu->covariances = psHashAlloc(psListLength(cells));
     }
 
@@ -448 +455 @@
 
             psArray *readouts = cell->readouts; // Array of readouts
+
+            // Check size of covariances
+            int xMinCovar = INT_MAX, xMaxCovar = INT_MIN, yMinCovar = INT_MAX, yMaxCovar = INT_MIN; // Size
+            if (covarianceType) {
+                for (int i = 0; i < readouts->n; i++) {
+                    pmReadout *readout = readouts->data[i]; // The readout of interest
+                    if (!readout) {
+                        continue;
+                    }
+                    psKernel *covar = readout->covariance; // Covariance matrix
+                    xMinCovar = PS_MIN(xMinCovar, covar->xMin);
+                    xMaxCovar = PS_MAX(xMaxCovar, covar->xMax);
+                    yMinCovar = PS_MIN(yMinCovar, covar->yMin);
+                    yMaxCovar = PS_MAX(yMaxCovar, covar->yMax);
+                }
+            }
+
             psArray *hduImages = hdu->images; // Array of images in the HDU
             psArray *hduMasks = hdu->masks; // Array of masks in the HDU
             psArray *hduVariances = hdu->variances; // Array of variances in the HDU
+            psArray *covariances = (covarianceType ? psArrayAlloc(readouts->n) : NULL; // Covariance images
             for (int i = 0; i < readouts->n; i++) {
                 pmReadout *readout = readouts->data[i]; // The readout of interest
@@ -471 +496 @@
                     psFree(readout->variance);
                     readout->variance = new;
+                }
+                if (readout->covariance) {
+                    psKernel *covar = readout->covariance; // Covariance matrix
+                    int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax;// Size
+                    if (xMin == xMinCovar && xMax == xMaxCovar && yMin == yMinCovar && yMax == yMaxCovar) {
+                        variances->data[i] = psMemIncrRefCounter(readout->covariance);
+                    } else {
+                        psKernel *new = psKernelAlloc(xMinCovar, xMaxCovar, yMinCovar, yMaxCovar);// New covar
+                        psImageInit(covar->image, 0);
+                        psImageOverlaySection(new->image, covar->image,
+                                              xMinCovar - xMin, yMinCovar - yMin, "=");
+                        variances->data[i] = new;
+                    }
                 }
 
@@ -495 +533 @@
             }
             psFree(biassecsIter);
+
+            if (covariances) {
+                const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
+                pmChip *chip = cell->parent; // Parent chip
+                const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
+                if (!cellName || !chipName) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to find cell or chip name.");
+                    return false;
+                }
+                psString name = NULL; // Name for covariance extension
+                psStringAppend(&name, "COVAR_%s_%s", chipName, cellName);
+                if (psHashLookup(hdu->covariances, name)) {
+                    psHashRemove(hash, name);
+                }
+                psHashAdd(hdu->covariances, name, covariance);
+                psFree(name);
+            }
         } // Iterating over cells within the HDU
         psFree(iter);
@@ -504 +559 @@
 
 // Return the level that an extension applies to
-static pmFPALevel extensionLevel(pmHDU *hdu // HDU to check
+static pmFPALevel extensionLevel(const pmHDU *hdu // HDU to check
                                 )
 {
@@ -512 +567 @@
     }
     bool mdok = true;                   // Status of MD lookup
-    psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File information for camera format
+    psMetadata *file = psMetadataLookupMetadata(&mdok, hdu->format, "FILE"); // File info for camera format
     if (!mdok || !file) {
         psError(PS_ERR_UNEXPECTED_NULL, true, "Can't file FILE information for camera format "
@@ -581 +636 @@
         return false;
     }
-    if (hdu->images && hdu->masks && hdu->variances) {
+    if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) {
         // It's already here!
         return true;
@@ -631 +686 @@
         return generateForCells(chip);
     }
-    if (hdu->images && hdu->masks && hdu->variances) {
+    if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) {
         // It's already here!
         return true;
@@ -680 +735 @@
         return generateForChips(fpa);
     }
-    if (hdu->images && hdu->masks && hdu->variances) {
+    if (hdu->images && hdu->masks && hdu->variances && hdu->covariances) {
         // It's already here!
         return true;