Changeset 21363 for trunk/psModules/src/camera

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.

Location:

trunk/psModules/src/camera

Files:

• pmFPA.c (modified) (9 diffs)
• pmFPA.h (modified) (6 diffs)
• pmFPAConstruct.c (modified) (3 diffs)
• pmFPACopy.c (modified) (3 diffs)
• pmFPAExtent.c (modified) (4 diffs)
• pmFPAMaskWeight.c (modified) (37 diffs)
• pmFPAMaskWeight.h (modified) (7 diffs)
• pmFPAMosaic.c (modified) (24 diffs)
• pmFPARead.c (modified) (20 diffs)
• pmFPARead.h (modified) (11 diffs)
• pmFPAWrite.c (modified) (13 diffs)
• pmFPAWrite.h (modified) (9 diffs)
• pmFPAfile.c (modified) (2 diffs)
• pmFPAfile.h (modified) (2 diffs)
• pmFPAfileFitsIO.c (modified) (4 diffs)
• pmFPAfileFitsIO.h (modified) (3 diffs)
• pmFPAfileIO.c (modified) (10 diffs)
• pmHDU.c (modified) (5 diffs)
• pmHDU.h (modified) (4 diffs)
• pmHDUGenerate.c (modified) (13 diffs)
• pmHDUUtils.c (modified) (2 diffs)
• pmReadoutStack.c (modified) (7 diffs)

trunk/psModules/src/camera/pmFPA.c

@@ -32 +32 @@
     psFree(readout->image);
     psFree(readout->mask);
-    psFree(readout->weight);
+    psFree(readout->variance);
+    psFree(readout->covariance);
     psFree(readout->analysis);
     psFree(readout->bias);
@@ -163 +164 @@
     psFree(readout->image);
     psFree(readout->mask);
-    psFree(readout->weight);
+    psFree(readout->variance);
+    psFree(readout->covariance);
     psFree(readout->bias);
 
@@ -169 +171 @@
 
     readout->image = NULL;
-    readout->weight = NULL;
+    readout->variance = NULL;
+    readout->covariance = NULL;
     readout->mask = NULL;
 
@@ -179 +182 @@
     readout->thisImageScan = 0;
     readout->thisMaskScan = 0;
-    readout->thisWeightScan = 0;
+    readout->thisVarianceScan = 0;
 
     readout->lastImageScan = 0;
     readout->lastMaskScan = 0;
-    readout->lastWeightScan = 0;
+    readout->lastVarianceScan = 0;
 }
 
@@ -197 +200 @@
     if (cell->hdu) {
         psFree(cell->hdu->images);
-        psFree(cell->hdu->weights);
+        psFree(cell->hdu->variances);
         psFree(cell->hdu->masks);
         // psFree(cell->hdu->header);
 
         cell->hdu->images = NULL;
-        cell->hdu->weights = NULL;
+        cell->hdu->variances = NULL;
         cell->hdu->masks = NULL;
         // cell->hdu->header = NULL;
@@ -219 +222 @@
     if (chip->hdu) {
         psFree(chip->hdu->images);
-        psFree(chip->hdu->weights);
+        psFree(chip->hdu->variances);
         psFree(chip->hdu->masks);
         // psFree(chip->hdu->header);
 
         chip->hdu->images = NULL;
-        chip->hdu->weights = NULL;
+        chip->hdu->variances = NULL;
         chip->hdu->masks = NULL;
         // chip->hdu->header = NULL;
@@ -241 +244 @@
     if (fpa->hdu) {
         psFree(fpa->hdu->images);
-        psFree(fpa->hdu->weights);
+        psFree(fpa->hdu->variances);
         psFree(fpa->hdu->masks);
         // psFree(fpa->hdu->header);
 
         fpa->hdu->images = NULL;
-        fpa->hdu->weights = NULL;
+        fpa->hdu->variances = NULL;
         fpa->hdu->masks = NULL;
         // fpa->hdu->header = NULL;
@@ -259 +262 @@
     tmpReadout->image = NULL;
     tmpReadout->mask = NULL;
-    tmpReadout->weight = NULL;
+    tmpReadout->variance = NULL;
+    tmpReadout->covariance = NULL;
     tmpReadout->bias = psListAlloc(NULL);
     tmpReadout->analysis = psMetadataAlloc();
@@ -276 +280 @@
     tmpReadout->thisImageScan = 0;
     tmpReadout->thisMaskScan = 0;
-    tmpReadout->thisWeightScan = 0;
+    tmpReadout->thisVarianceScan = 0;
 
     tmpReadout->lastImageScan = 0;
     tmpReadout->lastMaskScan = 0;
-    tmpReadout->lastWeightScan = 0;
+    tmpReadout->lastVarianceScan = 0;
 
     tmpReadout->forceScan = false;

trunk/psModules/src/camera/pmFPA.h

@@ -6 +6 @@
  * @author Eugene Magnier, IfA
  *
- * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-08-12 02:51:20 $
+ * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -113 +113 @@
 ///
 /// A readout corresponds to an individual read of a cell (e.g., a single image as part of a video sequence,
-/// or one of multiple coadds).  It contains the actual pixels used in analysis (along with mask and weight
+/// or one of multiple coadds).  It contains the actual pixels used in analysis (along with mask and variance
 /// maps).  When reading from a FITS file, the images are subimages (from CELL.TRIMSEC) of the pixels read
 /// from the appropriate HDU (at the FPA, chip or cell level).  The readout also contains a list of bias
@@ -123 +123 @@
     psImage *image;                     ///< Imaging area of readout (corresponds to CELL.TRIMSEC region)
     psImage *mask;                      ///< Mask of input image (corresponds to CELL.TRIMSEC region)
-    psImage *weight;                    ///< Weight of input image (corresponds to CELL.TRIMSEC region)
+    psImage *variance;                  ///< Variance of input image (corresponds to CELL.TRIMSEC region)
+    psKernel *covariance;               ///< Covariance pseudo-matrix (covariance factors for single pixel)
     psList *bias;                       ///< List of bias (prescan/overscan) images
     psMetadata *analysis;               ///< Readout-level analysis metadata
@@ -130 +131 @@
     bool file_exists;                   ///< Does the file for this readout exist (read case only)?
     bool data_exists;                   ///< Does the data for this readout exist (read case only)?
-    int thisImageScan;                  ///< start scan for next/current read
-    int lastImageScan;                  ///< start scan of the last read
-    int thisMaskScan;                   ///< start scan for next/current read
-    int lastMaskScan;                   ///< start scan of the last read
-    int thisWeightScan;                 ///< start scan for next/current read
-    int lastWeightScan;                 ///< start scan of the last read
+    int thisImageScan;                  ///< start scan for next/current read of image
+    int lastImageScan;                  ///< start scan of the last read of image
+    int thisMaskScan;                   ///< start scan for next/current read of mask
+    int lastMaskScan;                   ///< start scan of the last read of mask
+    int thisVarianceScan;               ///< start scan for next/current read of variance
+    int lastVarianceScan;               ///< start scan of the last read of variance
     bool forceScan;                     ///< Force pmFPARead to obey the above commanded this and last scans.
 } pmReadout;
@@ -217 +218 @@
         } \
     } \
-    psImage *weight = (READOUT)->weight; /* Weight map pixels */ \
-    if (weight) { \
-        PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL); \
-        if ((numCols != 0 || numRows != 0) && (weight->numCols != numCols || weight->numRows != numRows)) { \
-            psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Weight in readout %s has wrong size (%dx%d vs %dx%d)", \
-                    #READOUT, weight->numCols, weight->numRows, numCols, numRows); \
+    psImage *variance = (READOUT)->variance; /* Variance map pixels */ \
+    if (variance) { \
+        PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL); \
+        if ((numCols != 0 || numRows != 0) && \
+            (variance->numCols != numCols || variance->numRows != numRows)) { \
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true, \
+                    "Variance in readout %s has wrong size (%dx%d vs %dx%d)", \
+                    #READOUT, variance->numCols, variance->numRows, numCols, numRows); \
             return RETVAL; \
         } \
@@ -236 +239 @@
     PS_ASSERT_IMAGE_NON_NULL((READOUT)->mask, RETVAL);
 
-/// Assert that a readout contains a weight map
-#define PM_ASSERT_READOUT_WEIGHT(READOUT, RETVAL) \
-    PS_ASSERT_IMAGE_NON_NULL((READOUT)->weight, RETVAL);
+/// Assert that a readout contains a variance map
+#define PM_ASSERT_READOUT_VARIANCE(READOUT, RETVAL) \
+    PS_ASSERT_IMAGE_NON_NULL((READOUT)->variance, RETVAL);
+
+/// Assert that a readout contains a covariance matrix
+#define PM_ASSERT_READOUT_COVARIANCE(READOUT, RETVAL) \
+    PS_ASSERT_KERNEL_NON_NULL((READOUT)->covariance, RETVAL);
 
 /// @}

trunk/psModules/src/camera/pmFPAConstruct.c

@@ -1499 +1499 @@
                 psImage *image = readout->image; // The image
                 psImage *mask = readout->mask; // The mask
-                psImage *weight = readout->weight; // The weight
+                psImage *variance = readout->variance; // The variance
                 psList *bias = readout->bias; // The list of bias images
                 if (image) {
                     INDENT(fd, 7);
-                    fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n", image->col0, image->col0 +
-                            image->numCols, image->row0, image->row0 + image->numRows, image->numCols,
-                            image->numRows);
+                    fprintf(fd, "Image: [%d:%d,%d:%d] (%dx%d)\n",
+                            image->col0, image->col0 + image->numCols,
+                            image->row0, image->row0 + image->numRows,
+                            image->numCols, image->numRows);
                 }
                 if (bias) {
@@ -1512 +1513 @@
                     while ((biasImage = psListGetAndIncrement(biasIter))) {
                         INDENT(fd, 7);
-                        fprintf(fd, "Bias:  [%d:%d,%d:%d] (%dx%d)\n", biasImage->col0,
-                                biasImage->col0 + biasImage->numCols, biasImage->row0,
-                                biasImage->row0 + biasImage->numRows, biasImage->numCols, biasImage->numRows);
+                        fprintf(fd, "Bias:  [%d:%d,%d:%d] (%dx%d)\n",
+                                biasImage->col0, biasImage->col0 + biasImage->numCols,
+                                biasImage->row0, biasImage->row0 + biasImage->numRows,
+                                biasImage->numCols, biasImage->numRows);
                     }
                     psFree(biasIter);
@@ -1520 +1522 @@
                 if (mask) {
                     INDENT(fd, 7);
-                    fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n", mask->col0, mask->col0 +
-                            mask->numCols, mask->row0, mask->row0 + mask->numRows, mask->numCols,
-                            mask->numRows);
+                    fprintf(fd, "Mask: [%d:%d,%d:%d] (%dx%d)\n",
+                            mask->col0, mask->col0 + mask->numCols,
+                            mask->row0, mask->row0 + mask->numRows,
+                            mask->numCols, mask->numRows);
                 }
-                if (weight) {
+                if (variance) {
                     INDENT(fd, 7);
-                    fprintf(fd, "Weight: [%d:%d,%d:%d] (%dx%d)\n", weight->col0, weight->col0 +
-                            weight->numCols, weight->row0, weight->row0 + weight->numRows, weight->numCols,
-                            weight->numRows);
+                    fprintf(fd, "Variance: [%d:%d,%d:%d] (%dx%d)\n",
+                            variance->col0, variance->col0 + variance->numCols,
+                            variance->row0, variance->row0 + variance->numRows,
+                            variance->numCols, variance->numRows);
                 }
             } // Iterating over cell

trunk/psModules/src/camera/pmFPACopy.c

@@ -224 +224 @@
         targetReadout->data_exists = sourceReadout->data_exists;
 
-        // Copy all three image components (image, mask, weight)
-        readoutCopyComponent(&targetReadout->image,  sourceReadout->image,  binning, xFlip, yFlip, pixels);
-        readoutCopyComponent(&targetReadout->mask,   sourceReadout->mask,   binning, xFlip, yFlip, pixels);
-        readoutCopyComponent(&targetReadout->weight, sourceReadout->weight, binning, xFlip, yFlip, pixels);
+        // Copy all three image components (image, mask, variance)
+        readoutCopyComponent(&targetReadout->image, sourceReadout->image, binning, xFlip, yFlip, pixels);
+        readoutCopyComponent(&targetReadout->mask, sourceReadout->mask, binning, xFlip, yFlip, pixels);
+        readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip,
+                             pixels);
 
         // Copy bias
@@ -328 +329 @@
         int xSize = psMetadataLookupS32(NULL, source->concepts, "CELL.XSIZE"); // CELL.XSIZE of source
 
-        psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source");
+        psAssert (abs(xParity) == 1, "CELL.XPARITY not set for source");
 
         if (sourceBin == 0) {
@@ -478 +479 @@
         pmCell *targetCell = pmCellAlloc (targetChip, cellName);
         int xParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.XPARITY"); // Target x parity
-        psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target");
+        psAssert (abs(xParityTarget) == 1, "CELL.XPARITY not set for target");
         psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "", xParityTarget);
         int yParityTarget = psMetadataLookupS32(&status, sourceCell->concepts, "CELL.YPARITY"); // Target y parity
-        psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target");
+        psAssert (abs(xParityTarget) == 1, "CELL.YPARITY not set for target");
         psMetadataAddS32 (targetCell->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "", yParityTarget);
         if (!cellCopy(targetCell, sourceCell, true, 1, 1)) {

trunk/psModules/src/camera/pmFPAExtent.c

@@ -16 +16 @@
     }
     if (!image) {
-        image = readout->weight;
+        image = readout->variance;
     }
 
@@ -23 +23 @@
 
     if (!image) {
-        pmHDU *hdu = pmHDUFromReadout (readout);
-        if (hdu && hdu->header) {
-            bool status;
-            xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1");
-            if (!status) {
-                xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1");
-                if (!status) return NULL;
-            }           
-            ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2");
-            if (!status) {
-                ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2");
-                if (!status) return NULL;
-            }           
-        } else {
+        pmHDU *hdu = pmHDUFromReadout (readout);
+        if (hdu && hdu->header) {
+            bool status;
+            xSize = psMetadataLookupS32(&status, hdu->header, "NAXIS1");
+            if (!status) {
+                xSize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS1");
+                if (!status) return NULL;
+            }
+            ySize = psMetadataLookupS32(&status, hdu->header, "NAXIS2");
+            if (!status) {
+                ySize = psMetadataLookupS32(&status, hdu->header, "IMNAXIS2");
+                if (!status) return NULL;
+            }
+        } else {
         // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent)
-            xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE");
-            ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE");
-        }
+            xSize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.XSIZE");
+            ySize = psMetadataLookupS32(NULL, readout->parent->concepts, "CELL.YSIZE");
+        }
         return psRegionAlloc(0, xSize, 0, ySize);
     }
@@ -70 +70 @@
     // Don't have anything to base the true extent on, so have to give the hardwired value (largest possible extent)
     if (readouts->n == 0) {
-        int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
-        int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");
-        cellExtent->x0 = 0;
-        cellExtent->x1 = xSize;
-        cellExtent->y0 = 0;
-        cellExtent->y1 = ySize;
+        int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
+        int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");
+        cellExtent->x0 = 0;
+        cellExtent->x1 = xSize;
+        cellExtent->y0 = 0;
+        cellExtent->y1 = ySize;
     }
 
@@ -97 +97 @@
     // CELL.X0,Y0 are the coordinate of the amp on the chip, subtract size if parity flipped
     if (xParityCell > 0) {
-        cellExtent->x0 += x0;
-        cellExtent->x1 += x0;
+        cellExtent->x0 += x0;
+        cellExtent->x1 += x0;
     } else {
-        float x0Cell = x0 - cellExtent->x1;
-        float x1Cell = x0 - cellExtent->x0;
-        cellExtent->x0 = x0Cell;
-        cellExtent->x1 = x1Cell;
+        float x0Cell = x0 - cellExtent->x1;
+        float x1Cell = x0 - cellExtent->x0;
+        cellExtent->x0 = x0Cell;
+        cellExtent->x1 = x1Cell;
     }
     if (yParityCell > 0) {
-        cellExtent->y0 += y0;
-        cellExtent->y1 += y0;
+        cellExtent->y0 += y0;
+        cellExtent->y1 += y0;
     } else {
-        float y0Cell = y0 - cellExtent->y1;
-        float y1Cell = y0 - cellExtent->y0;
-        cellExtent->y0 = y0Cell;
-        cellExtent->y1 = y1Cell;
+        float y0Cell = y0 - cellExtent->y1;
+        float y1Cell = y0 - cellExtent->y0;
+        cellExtent->y0 = y0Cell;
+        cellExtent->y1 = y1Cell;
     }
 

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -49 +49 @@
 }
 
-// Create the parent weight images that reside in the HDU
-static void createParentWeights(pmHDU *hdu // The HDU for which to create
+// Create the parent variance images that reside in the HDU
+static void createParentVariances(pmHDU *hdu // The HDU for which to create
                                )
 {
@@ -58 +58 @@
     // Generate the parent mask images
     psArray *images = hdu->images;      // Array of images
-    psArray *weights = hdu->weights;    // Array of weight images
-    if (!weights) {
-        weights = psArrayAlloc(images->n);
-        hdu->weights = weights;
+    psArray *variances = hdu->variances;    // Array of variance images
+    if (!variances) {
+        variances = psArrayAlloc(images->n);
+        hdu->variances = variances;
     }
 
     for (long i = 0; i < images->n; i++) {
         psImage *image = images->data[i]; // The image for this readout
-        if (!image || weights->data[i]) {
+        if (!image || variances->data[i]) {
             continue;
         }
-        weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+        variances->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     }
 
@@ -199 +199 @@
 }
 
-bool pmReadoutSetWeight(pmReadout *readout, bool poisson)
+bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -209 +209 @@
     float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
     if (!mdok || isnan(gain)) {
-        psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n");
         return false;
     }
     float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
     if (!mdok || isnan(readnoise)) {
-        psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n");
         return false;
     }
@@ -223 +223 @@
 
     if (poisson) {
-        // Set weight image to the variance in ADU = f/g + rn^2
+        // Set variance image to the variance in ADU = f/g + rn^2
         psImage *image = readout->image;    // The image pixels
-        readout->weight = (psImage*)psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32));
-
-        // a negative weight is non-sensical. if the image value drops below 1, the weight must be 1.
-        readout->weight = (psImage*)psUnaryOp(readout->weight, readout->weight, "abs");
-        readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "max",
+        readout->variance = (psImage*)psBinaryOp(readout->variance, image, "/", psScalarAlloc(gain, PS_TYPE_F32));
+
+        // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
+        readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
                                                psScalarAlloc(1, PS_TYPE_F32));
     } else {
         // Just use the read noise
-        if (!readout->weight) {
-            readout->weight = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
-        }
-        psImageInit(readout->weight, 0.0);
-    }
-
-    readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "+",
+        if (!readout->variance) {
+            readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+        }
+        psImageInit(readout->variance, 0.0);
+    }
+
+    readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+",
                                            psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
 
@@ -245 +245 @@
 }
 
-// this function creates the weight pixels, or uses the existing weight pixels.  it will set
-// the noise pixel values only if the weight image is not supplied
-bool pmReadoutGenerateWeight(pmReadout *readout, bool poisson)
+// this function creates the variance pixels, or uses the existing variance pixels.  it will set
+// the noise pixel values only if the variance image is not supplied
+bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -254 +254 @@
     bool mdok = true;                   // Status of MD lookup
 
-    // Create the weight image if required
-    if (readout->weight)
+    // Create the variance image if required
+    if (readout->variance)
         return true;
 
     psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
     if (!mdok || psRegionIsNaN(*trimsec)) {
-        psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n");
         return false;
     }
@@ -271 +271 @@
     }
 
-    createParentWeights(hdu);
+    createParentVariances(hdu);
 
     // Need to identify which readout we're working with....
@@ -280 +280 @@
     }
 
-    psImage *weight = psImageSubset(hdu->weights->data[index], *trimsec); // The weight pixels
-    if (!weight) {
+    psImage *variance = psImageSubset(hdu->variances->data[index], *trimsec); // The variance pixels
+    if (!variance) {
         psString trimsecString = psRegionToString(*trimsec);
-        psError(PS_ERR_UNKNOWN, false, "Unable to set weight from HDU with trimsec: %s.\n",
+        psError(PS_ERR_UNKNOWN, false, "Unable to set variance from HDU with trimsec: %s.\n",
                 trimsecString);
         psFree(trimsecString);
         return false;
     }
-    psImageInit(weight, 0);
-    readout->weight = weight;
-
-    return pmReadoutSetWeight(readout, poisson);
-}
-
-bool pmReadoutGenerateMaskWeight(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
+    psImageInit(variance, 0);
+    readout->variance = variance;
+
+    return pmReadoutSetVariance(readout, poisson);
+}
+
+bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -301 +301 @@
 
     success &= pmReadoutGenerateMask(readout, satMask, badMask);
-    success &= pmReadoutGenerateWeight(readout, poisson);
+    success &= pmReadoutGenerateVariance(readout, poisson);
 
     return success;
 }
 
-bool pmCellGenerateMaskWeight(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
+bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(cell, false);
@@ -314 +314 @@
     for (int i = 0; i < readouts->n; i++) {
         pmReadout *readout = readouts->data[i]; // The readout
-        success &= pmReadoutGenerateMaskWeight(readout, poisson, satMask, badMask);
+        success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask);
     }
 
@@ -321 +321 @@
 
 
-bool pmReadoutWeightRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
+bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
                                  psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
     PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_WEIGHT(readout, false);
-
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout parts
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
+
+    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
 
     if (!psMemIncrRefCounter(rng)) {
@@ -334 +334 @@
     }
 
-    psStats *weightStats = psStatsAlloc(meanStat);// Statistics for mean
-    if (!psImageBackground(weightStats, NULL, weight, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean weight for image");
-        psFree(weightStats);
+    psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean
+    if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image");
+        psFree(varianceStats);
         psFree(rng);
         return false;
     }
-    float meanVariance = weightStats->robustMedian; // Mean variance
-    psFree(weightStats);
+    float meanVariance = varianceStats->robustMedian; // Mean variance
+    psFree(varianceStats);
 
     psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
@@ -356 +356 @@
 
     float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", correction);
-    psBinaryOp(weight, weight, "*", psScalarAlloc(correction, PS_TYPE_F32));
+    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction);
+    psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32));
 
     return true;
@@ -363 +363 @@
 
 
-bool pmReadoutWeightRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
+bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
                                psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
     PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_WEIGHT(readout, false);
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
 
     if (!psMemIncrRefCounter(rng)) {
@@ -374 +374 @@
     }
 
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
+    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
 
     // Measure background
@@ -424 +424 @@
         float sumNoise = 0.0;       // Sum for noise measurement
         float sumSource = 0.0;      // Sum for source measurement
-        float sumWeight = 0.0;      // Sum for weight measurement
+        float sumVariance = 0.0;      // Sum for variance measurement
         float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels
         for (int v = 0, y = yPix - size; v < fullSize; v++, y++) {
             float xSumNoise = 0.0;  // Sum for noise measurement in x
             float xSumSource = 0.0; // Sum for source measurement in x
-            float xSumWeight = 0.0; // Sum for weight measurement in x
+            float xSumVariance = 0.0; // Sum for variance measurement in x
             float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x
             float yGauss = gauss->data.F32[v]; // Value of Gaussian in y
@@ -441 +441 @@
                 xSumNoise += value * xGauss;
                 xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss;
-                xSumWeight += weight->data.F32[y][x] * xGauss2;
+                xSumVariance += variance->data.F32[y][x] * xGauss2;
                 xSumGauss += xGauss;
                 xSumGauss2 += xGauss2;
@@ -448 +448 @@
             sumNoise += xSumNoise * yGauss;
             sumSource += xSumSource * yGauss;
-            sumWeight += xSumWeight * yGauss2;
+            sumVariance += xSumVariance * yGauss2;
             sumGauss += xSumGauss * yGauss;
             sumGauss2 += xSumGauss2 * yGauss2;
@@ -454 +454 @@
 
         photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
-                                                isfinite(sumWeight) && sumGauss > 0 && sumGauss2 > 0) ?
+                                                isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ?
                                                0 : 0xFF);
 
         float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise
         float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source
-        float smoothWeight = sumWeight / sumGauss2; // Value of smoothed weight pixel
+        float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel
 
         noise->data.F32[i] = smoothImageNoise;
         source->data.F32[i] = smoothImageSource;
-        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothWeight : 0.0;
+        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0;
         psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
-                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothWeight);
+                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
     }
     psFree(gauss);
@@ -516 +516 @@
     psFree(photMask);
 
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", meanRatio);
-    psBinaryOp(weight, weight, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
+    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio);
+    psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
 
     return true;
@@ -523 +523 @@
 
 
-bool pmReadoutWeightRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
+bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
                            psStatsOptions stdevStat, int width, psRandom *rng)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
     PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_WEIGHT(readout, false);
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
     PS_ASSERT_INT_POSITIVE(width, false);
 
@@ -535 +535 @@
     }
 
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
+    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
     int numCols = image->numCols, numRows = image->numRows; // Size of images
     int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions
@@ -554 +554 @@
             psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
             psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels
-            psImage *subWeight = psImageSubset(weight, region); // Sub image of the weight pixels
+            psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels
             psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels
 
             if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) ||
-                !psImageBackground(meanStats, &buffer, subWeight, subMask, maskVal, rng)) {
+                !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) {
                 // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it.
                 psString regionStr = psRegionToString(region); // String with region
@@ -564 +564 @@
                 psFree(regionStr);
                 psErrorClear();
-                psImageInit(subWeight, NAN);
+                psImageInit(subVariance, NAN);
                 if (subMask) {
                     psImageInit(subMask, maskVal);
@@ -574 +574 @@
                         "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n",
                         xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar);
-                psBinaryOp(subWeight, subWeight, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
+                psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
             }
 
             psFree(subImage);
-            psFree(subWeight);
+            psFree(subVariance);
             psFree(subMask);
         }
@@ -597 +597 @@
 
     psImage *image = readout->image;    // Readout's image
-    psImage *weight = readout->weight;  // Readout's weight
+    psImage *variance = readout->variance;  // Readout's variance
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
@@ -607 +607 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (!isfinite(image->data.F32[y][x]) || (weight && !isfinite(weight->data.F32[y][x]))) {
+            if (!isfinite(image->data.F32[y][x]) || (variance && !isfinite(variance->data.F32[y][x]))) {
                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
             }
@@ -627 +627 @@
     psImageMaskType **maskData = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Dereference mask
     psF32 **imageData = readout->image->data.F32;// Dereference image
-    psF32 **weightData = readout->weight ? readout->weight->data.F32 : NULL; // Dereference weight map
+    psF32 **varianceData = readout->variance ? readout->variance->data.F32 : NULL; // Dereference variance map
 
     for (int y = 0; y < numRows; y++) {
@@ -633 +633 @@
             if (maskData[y][x] & maskVal) {
                 imageData[y][x] = NAN;
-                if (weightData) {
-                    weightData[y][x] = NAN;
+                if (varianceData) {
+                    varianceData[y][x] = NAN;
                 }
             }
@@ -656 +656 @@
     psImage *image = readout->image;    // Image
     psImage *mask = readout->mask;      // Mask
-    psImage *weight = readout->weight;  // Weight map
-
-    psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal,
+    psImage *variance = readout->variance;  // Variance map
+
+    psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal,
                                                              NAN, NAN, maskBad, maskPoor, poorFrac, 0);
     interp->shifting = false;           // Turn off "exact shifts" so we get proper interpolation
@@ -666 +666 @@
     psPixels *pixels = psPixelsAllocEmpty(PIXELS_BUFFER); // Pixels that have been interpolated
     psVector *imagePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for image
-    psVector *weightPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for weight
+    psVector *variancePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for variance
     psVector *maskPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_IMAGE_MASK); // Corresponding values for mask
     // NOTE: maskPix carries the actual image mask values -- do NOT use
@@ -675 +675 @@
         for (int x = 0; x < numCols; x++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
-                double imageValue, weightValue; // Image and weight value from interpolation
+                double imageValue, varianceValue; // Image and variance value from interpolation
                 psImageMaskType maskValue = 0; // Mask value from interpolation
-                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &weightValue, &maskValue, x, y, interp);
+                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp);
                 if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y);
@@ -683 +683 @@
                     psFree(pixels);
                     psFree(imagePix);
-                    psFree(weightPix);
+                    psFree(variancePix);
                     psFree(maskPix);
                     return false;
@@ -694 +694 @@
                 pixels = psPixelsAdd(pixels, PIXELS_BUFFER, x, y);
                 imagePix = psVectorExtend(imagePix, PIXELS_BUFFER, 1);
-                weightPix = psVectorExtend(weightPix, PIXELS_BUFFER, 1);
+                variancePix = psVectorExtend(variancePix, PIXELS_BUFFER, 1);
                 maskPix = psVectorExtend(maskPix, PIXELS_BUFFER, 1);
                 imagePix->data.F32[numBad] = imageValue;
-                weightPix->data.F32[numBad] = weightValue;
+                variancePix->data.F32[numBad] = varianceValue;
                 maskPix->data.PS_TYPE_IMAGE_MASK_DATA[numBad] = maskValue;
                 numBad++;
@@ -709 +709 @@
         int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of pixel
         image->data.F32[y][x] = imagePix->data.F32[i];
-        weight->data.F32[y][x] = weightPix->data.F32[i];
+        variance->data.F32[y][x] = variancePix->data.F32[i];
         mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskPix->data.PS_TYPE_IMAGE_MASK_DATA[i];
     }
@@ -715 +715 @@
     psFree(pixels);
     psFree(imagePix);
-    psFree(weightPix);
+    psFree(variancePix);
     psFree(maskPix);
 

trunk/psModules/src/camera/pmFPAMaskWeight.h

@@ -5 +5 @@
  * @author Eugene Magnier, IfA
  *
- * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-27 06:39:38 $
+ * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -47 +47 @@
 
 
-/// Set a temporary readout weight map using CELL.GAIN and CELL.READNOISE
+/// Set a temporary readout variance map using CELL.GAIN and CELL.READNOISE
 ///
-/// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain
-/// (CELL.GAIN) and read noise (CELL.READNOISE).  The weight map that is produced within the readout is
-/// temporary --- it is not added to the HDU.  This is intended for when the user is iterating using
-/// pmReadoutReadNext, in which case the HDU can't be generated.
-bool pmReadoutSetWeight(pmReadout *readout, ///< Readout for which to set weight
-                        bool poisson    ///< Use poisson weights (in addition to read noise)?
-                       );
+/// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise
+/// (CELL.READNOISE).  The weight map that is produced within the readout is temporary --- it is not added to
+/// the HDU.  This is intended for when the user is iterating using pmReadoutReadNext, in which case the HDU
+/// can't be generated.
+bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance
+                          bool poisson    ///< Include poisson variance (in addition to read noise)?
+    );
 
 /// Generate a readout mask (suitable for output) using CELL.SATURATION and CELL.BAD
@@ -66 +66 @@
     );
 
-/// Generate a weight map (suitable for output) using CELL.GAIN and CELL.READNOISE
+/// Generate a variance map (suitable for output) using CELL.GAIN and CELL.READNOISE
 ///
-/// Calculates weights (actually variances) for each pixel using photon statistics and the cell gain
-/// (CELL.GAIN) and read noise (CELL.READNOISE).  The weight map that is produced within the readout is
-/// suitable for output (complete with HDU entry).  This is intended for most operations.
-bool pmReadoutGenerateWeight(pmReadout *readout, ///< Readout for which to generate weight
-                             bool poisson    ///< Use poisson weights (in addition to read noise)?
-                            );
+/// Calculates variances for each pixel using photon statistics and the cell gain (CELL.GAIN) and read noise
+/// (CELL.READNOISE).  The variance map that is produced within the readout is suitable for output (complete
+/// with HDU entry).  This is intended for most operations.
+bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance
+                               bool poisson    ///< Include poisson variance (in addition to read noise)?
+    );
 
-/// Generate mask and weight map for a readout
+/// Generate mask and variance map for a readout
 ///
-/// Calls pmReadoutGenerateMask and pmReadoutGenerateWeight for the readout
-bool pmReadoutGenerateMaskWeight(pmReadout *readout, ///< Readout for which to generate mask and weights
-                                 psImageMaskType sat, ///< Mask value to give saturated pixels
-                                 psImageMaskType bad, ///< Mask value to give bad (low) pixels
-                                 bool poisson ///< Use poisson weights (in addition to read noise)?
-                                );
+/// Calls pmReadoutGenerateMask and pmReadoutGenerateVariance for the readout
+bool pmReadoutGenerateMaskVariance(pmReadout *readout, ///< Readout for which to generate mask and variance
+                                   psImageMaskType sat, ///< Mask value to give saturated pixels
+                                   psImageMaskType bad, ///< Mask value to give bad (low) pixels
+                                   bool poisson ///< Include poisson variance (in addition to read noise)?
+    );
 
-/// Generate mask and weight maps for all readouts within a cell
+/// Generate mask and variance maps for all readouts within a cell
 ///
-/// Calls pmReadoutGenerateMaskWeight for each readout within the cell.
-bool pmCellGenerateMaskWeight(pmCell *cell, ///< Cell for which to generate mask and weights
-                              psImageMaskType sat, ///< Mask value to give saturated pixels
-                              psImageMaskType bad, ///< Mask value to give bad (low) pixels
-                              bool poisson ///< Use poisson weights (in addition to read noise)?
-                             );
+/// Calls pmReadoutGenerateMaskVariance for each readout within the cell.
+bool pmCellGenerateMaskVariance(pmCell *cell, ///< Cell for which to generate mask and variance
+                                psImageMaskType sat, ///< Mask value to give saturated pixels
+                                psImageMaskType bad, ///< Mask value to give bad (low) pixels
+                                bool poisson ///< Include poisson variance (in addition to read noise)?
+    );
 
-/// Renormalise the weight map to match the actual pixel variance
+/// Renormalise the variance map to match the actual pixel variance
 ///
-/// The weight (variance) map is adjusted so that the mean matches the actual pixel variance in the image
-bool pmReadoutWeightRenormPixels(
+/// The variance map is adjusted so that the mean matches the actual pixel variance in the image
+bool pmReadoutVarianceRenormPixels(
     const pmReadout *readout,           ///< Readout to normalise
     psImageMaskType maskVal,                 ///< Value to mask
@@ -104 +104 @@
     );
 
-/// Renormalise the weight map to match the actual photometry variance
+/// Renormalise the variance map to match the actual photometry variance
 ///
-/// The weight (variance) map is adjusted so that the actual significance of fake sources matches the
+/// The variance map is adjusted so that the actual significance of fake sources matches the
 /// guestimated significance
-bool pmReadoutWeightRenormPhot(
+bool pmReadoutVarianceRenormPhot(
     const pmReadout *readout,           ///< Readout to normalise
-    psImageMaskType maskVal,                 ///< Value to mask
+    psImageMaskType maskVal,            ///< Value to mask
     int num,                            ///< Number of instances to measure over the image
     float width,                        ///< Photometry width
@@ -118 +118 @@
     );
 
-/// Renormalise the weight map to match the actual variance
+/// Renormalise the variance map to match the actual variance
 ///
 /// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for
 /// that patch corresponds.
-bool pmReadoutWeightRenorm(const pmReadout *readout, // Readout to normalise
-                           psImageMaskType maskVal, // Value to mask
-                           psStatsOptions meanStat, // Statistic to measure the mean (of the variance map)
-                           psStatsOptions stdevStat, // Statistic to measure the stdev (of the image)
-                           int width,   // Width of patch (pixels)
-                           psRandom *rng // Random number generator (for sub-sampling images)
+bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise
+                             psImageMaskType maskVal, // Value to mask
+                             psStatsOptions meanStat, // Statistic to measure the mean (of the variance map)
+                             psStatsOptions stdevStat, // Statistic to measure the stdev (of the image)
+                             int width,   // Width of patch (pixels)
+                             psRandom *rng // Random number generator (for sub-sampling images)
     );
 
@@ -133 +133 @@
 ///
 /// Since unmasked non-finite pixels can occur (e.g., by out-of-range in quantisation), it is sometimes
-/// necessary to mask them explicitly.  Non-finite pixels in the image or weight have their mask OR-ed with
+/// necessary to mask them explicitly.  Non-finite pixels in the image or variance have their mask OR-ed with
 /// the provided value.
 bool pmReadoutMaskNonfinite(pmReadout *readout, ///< Readout to mask
@@ -139 +139 @@
     );
 
-/// Apply a mask to the image and weight map
+/// Apply a mask to the image and variance map
 ///
 /// Unfortunately, image subtraction may result in a bi-modal image in masked areas, which can upset image
 /// statistics (very important for quantising images so that a product can be written out!).  This function
-/// sets masked areas to NAN in the image and weight.
+/// sets masked areas to NAN in the image and variance.
 bool pmReadoutMaskApply(pmReadout *readout, ///< Readout to mask
                         psImageMaskType maskVal ///< Mask value for which to apply mask

trunk/psModules/src/camera/pmFPAMosaic.c

@@ -554 +554 @@
         // We have to do all of the hard work ourselves
         switch (type) {
-            FILL_IN(U8);
-            FILL_IN(U16);
-            FILL_IN(U32);
-            FILL_IN(U64);
-            FILL_IN(S8);
-            FILL_IN(S16);
-            FILL_IN(S32);
-            FILL_IN(S64);
-            FILL_IN(F32);
-            FILL_IN(F64);
+            FILL_IN(U8);
+            FILL_IN(U16);
+            FILL_IN(U32);
+            FILL_IN(U64);
+            FILL_IN(S8);
+            FILL_IN(S16);
+            FILL_IN(S32);
+            FILL_IN(S64);
+            FILL_IN(F32);
+            FILL_IN(F64);
           default:
             psAbort("Should never get here.\n");
@@ -575 +575 @@
 static bool addCell(psArray *images,    // Array of images
                     psArray *masks,     // Array of masks
-                    psArray *weights,   // Array of weights
+                    psArray *variances,   // Array of variances
                     psVector *x0,       // Array of X0
                     psVector *y0,       // Array of Y0
@@ -604 +604 @@
         images  = psArrayRealloc(images,  index + CELL_LIST_BUFFER);
         masks   = psArrayRealloc(masks,   index + CELL_LIST_BUFFER);
-        weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER);
+        variances = psArrayRealloc(variances, index + CELL_LIST_BUFFER);
         x0    = psVectorRealloc(x0,    index+ CELL_LIST_BUFFER);
         y0    = psVectorRealloc(y0,    index+ CELL_LIST_BUFFER);
@@ -615 +615 @@
     images->n = index + 1;
     masks->n = index + 1;
-    weights->n = index + 1;
+    variances->n = index + 1;
     x0->n = index + 1;
     y0->n = index + 1;
@@ -726 +726 @@
     // The images to put into the mosaic
     images->data[index]  = psMemIncrRefCounter(readout->image);
-    weights->data[index] = psMemIncrRefCounter(readout->weight);
+    variances->data[index] = psMemIncrRefCounter(readout->variance);
     masks->data[index]   = psMemIncrRefCounter(readout->mask);
 
@@ -740 +740 @@
 static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                        psImage **mosaicMask, // The mosaic mask, to be returned
-                       psImage **mosaicWeights, // The mosaic weights, to be returned
+                       psImage **mosaicVariances, // The mosaic variances, to be returned
                        int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
                        const pmChip *chip, // Chip to mosaic
@@ -749 +749 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicWeights);
+    assert(mosaicVariances);
     assert(xBinChip);
     assert(yBinChip);
@@ -756 +756 @@
 
     psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
-    psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
+    psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked
     psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
     psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
@@ -802 +802 @@
             continue;
         }
-        allGood &= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
+        allGood &= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip,
                            cell, xBinChip, yBinChip, false, x0Target, y0Target,
                            xParityTarget, yParityTarget);
@@ -826 +826 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
-        *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
+        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank);
     }
@@ -832 +832 @@
     // Clean up
     psFree(images);
-    psFree(weights);
+    psFree(variances);
     psFree(masks);
     psFree(xFlip);
@@ -847 +847 @@
 static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                       psImage **mosaicMask, // The mosaic mask, to be returned
-                      psImage **mosaicWeights, // The mosaic weights, to be returned
+                      psImage **mosaicVariances, // The mosaic variances, to be returned
                       int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
                       const pmFPA *fpa,  // FPA to mosaic
@@ -857 +857 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicWeights);
+    assert(mosaicVariances);
     assert(xBinFPA);
     assert(yBinFPA);
@@ -865 +865 @@
 
     psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
-    psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
+    psArray *variances = psArrayAlloc(0); // Array of variance images to be mosaicked
     psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
     psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
@@ -941 +941 @@
                 continue;
             }
-            allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
+            allGood |= addCell(images, masks, variances, x0, y0, xBin, yBin, xFlip, yFlip,
                                cell, xBinFPA, yBinFPA, true, x0Target, y0Target,
                                xParityTarget, yParityTarget);
@@ -960 +960 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
-        *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
+        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank);
     }
@@ -966 +966 @@
     // Clean up
     psFree(images);
-    psFree(weights);
+    psFree(variances);
     psFree(masks);
     psFree(xFlip);
@@ -1025 +1025 @@
     psImage *mosaicImage   = NULL;      // The mosaic image
     psImage *mosaicMask    = NULL;      // The mosaic mask
-    psImage *mosaicWeights = NULL;      // The mosaic weights
+    psImage *mosaicVariances = NULL;      // The mosaic variances
 
     // Find the HDU
@@ -1051 +1051 @@
             }
         }
-        if (hdu->weights) {
-            mosaicWeights = psImageSubset(hdu->weights->data[0], bounds);
-            if (!mosaicWeights) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to select weight pixels.\n");
+        if (hdu->variances) {
+            mosaicVariances = psImageSubset(hdu->variances->data[0], bounds);
+            if (!mosaicVariances) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n");
                 return false;
             }
@@ -1061 +1061 @@
         // Case 2 --- we need to mosaic by cut and paste
         psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source, targetCell, blank)) {
+        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
             return false;
@@ -1094 +1094 @@
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
-    newReadout->weight = mosaicWeights;
+    newReadout->variance = mosaicVariances;
     psFree(newReadout);                 // Drop reference
 
@@ -1262 +1262 @@
     psImage *mosaicImage   = NULL;      // The mosaic image
     psImage *mosaicMask    = NULL;      // The mosaic mask
-    psImage *mosaicWeights = NULL;      // The mosaic weights
+    psImage *mosaicVariances = NULL;      // The mosaic variances
 
     // Find the HDU
@@ -1275 +1275 @@
             mosaicMask = psImageSubset(hdu->masks->data[0], *fpaRegion);
         }
-        if (hdu->weights) {
-            mosaicWeights = psImageSubset(hdu->weights->data[0], *fpaRegion);
+        if (hdu->variances) {
+            mosaicVariances = psImageSubset(hdu->variances->data[0], *fpaRegion);
         }
     } else {
         // Case 2 --- we need to mosaic by cut and paste
         psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source,
+        if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source,
                        targetChip, targetCell, blank)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n");
@@ -1338 +1338 @@
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
-    newReadout->weight = mosaicWeights;
+    newReadout->variance = mosaicVariances;
     psFree(newReadout);                 // Drop reference
 

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;
+}

trunk/psModules/src/camera/pmFPARead.h

@@ -4 +4 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-06-17 22:16:38 $
+ * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -60 +60 @@
                        int numRows,      // The number of rows to read
                        pmConfig *config
-                      );
+    );
 
 /// Return the number of readouts within a cell
 ///
-/// This function is type-independent (doesn't matter if you are interested in the image/mask/weight).
+/// This function is type-independent (doesn't matter if you are interested in the image/mask/variance).
 int pmCellNumReadouts(pmCell *cell, psFits *fits, pmConfig *config);
 
@@ -74 +74 @@
                 psFits *fits,           // FITS file from which to read
                 pmConfig *config        // Configuration
-               );
+    );
 
 /// Read a chip
@@ -82 +82 @@
                 psFits *fits,           // FITS file from which to read
                 pmConfig *config        // Configuration
-               );
+    );
 
 /// Read an FPA
@@ -90 +90 @@
                psFits *fits,            // FITS file from which to read
                pmConfig *config         // Configuration
-              );
+    );
 
 // Mask functions follow
@@ -126 +126 @@
                     psFits *fits,       // FITS file from which to read
                     pmConfig *config    // Configuration
-                   );
+    );
 
 /// Read an entire chip into the mask
@@ -134 +134 @@
                     psFits *fits,       // FITS file from which to read
                     pmConfig *config    // Configuration
-                   );
+    );
 
 /// Read an entire FPA into the mask
@@ -142 +142 @@
                    psFits *fits,        // FITS file from which to read
                    pmConfig *config     // Configuration
-                  );
-
-// Weight functions follow
-
-/// Check to see if there is more to read when reading a readout incrementally into the weight
-bool pmReadoutMoreWeight(pmReadout *readout, ///< Readout of interest
-                         psFits *fits,  ///< FITS file from which to read
-                         int z,         ///< Readout number/plane; zero-offset indexing
-                         int numScans,  ///< Number of scans (rows/cols) to read
-                         pmConfig *config ///< Configuration
-    );
-
-/// Read a chunk of a readout into the weight
+    );
+
+// Variance functions follow
+
+/// Check to see if there is more to read when reading a readout incrementally into the variance
+bool pmReadoutMoreVariance(pmReadout *readout, ///< Readout of interest
+                           psFits *fits, ///< FITS file from which to read
+                           int z,       ///< Readout number/plane; zero-offset indexing
+                           int numScans, ///< Number of scans (rows/cols) to read
+                           pmConfig *config ///< Configuration
+    );
+
+/// Read a chunk of a readout into the variance
 ///
 /// Allows reading the readout incrementally
-bool pmReadoutReadChunkWeight(pmReadout *readout, ///< Readout of interest
-                              psFits *fits, ///< FITS file from which to read
-                              int z,    ///< Readout number/plane; zero-offset indexing
-                              int numScans, ///< Number of scans (rows/cols) to read
-                              int overlap, ///< Overlap between consecutive reads
-                              pmConfig *config ///< Configuration
-    );
-
-/// Read the entire readout into the weight
-bool pmReadoutReadWeight(pmReadout *readout, ///< Readout of interest
-                         psFits *fits,  ///< FITS file from which to read
-                         int z,         ///< Readout number/plane; zero-offset indexing
-                         pmConfig *config ///< Configuration
-    );
-
-/// Read an entire cell into the weight
-///
-/// Same as pmCellRead, but reads into the weight element of the readouts.
-bool pmCellReadWeight(pmCell *cell,     // Cell to read into
-                      psFits *fits,     // FITS file from which to read
-                      pmConfig *config  // Configuration
-                     );
-
-/// Read an entire chip into the weight
-///
-/// Same as pmChipRead, but reads into the weight element of the readouts.
-bool pmChipReadWeight(pmChip *chip,     // Chip to read into
-                      psFits *fits,     // FITS file from which to read
-                      pmConfig *config  // Configuration
-                     );
-
-/// Read an entire FPA into the weight
-///
-/// Same as pmFPARead, but reads into the weight element of the readouts.
-bool pmFPAReadWeight(pmFPA *fpa,        // FPA to read into
-                     psFits *fits,      // FITS file from which to read
-                     pmConfig *config   // Configuration
-                    );
+bool pmReadoutReadChunkVariance(pmReadout *readout, ///< Readout of interest
+                                psFits *fits, ///< FITS file from which to read
+                                int z,    ///< Readout number/plane; zero-offset indexing
+                                int numScans, ///< Number of scans (rows/cols) to read
+                                int overlap, ///< Overlap between consecutive reads
+                                pmConfig *config ///< Configuration
+    );
+
+/// Read the entire readout into the variance
+bool pmReadoutReadVariance(pmReadout *readout, ///< Readout of interest
+                           psFits *fits,  ///< FITS file from which to read
+                           int z,         ///< Readout number/plane; zero-offset indexing
+                           pmConfig *config ///< Configuration
+    );
+
+/// Read an entire cell into the variance
+///
+/// Same as pmCellRead, but reads into the variance element of the readouts.
+bool pmCellReadVariance(pmCell *cell,     // Cell to read into
+                        psFits *fits,     // FITS file from which to read
+                        pmConfig *config  // Configuration
+    );
+
+/// Read an entire chip into the variance
+///
+/// Same as pmChipRead, but reads into the variance element of the readouts.
+bool pmChipReadVariance(pmChip *chip,     // Chip to read into
+                        psFits *fits,     // FITS file from which to read
+                        pmConfig *config  // Configuration
+    );
+
+/// Read an entire FPA into the variance
+///
+/// Same as pmFPARead, but reads into the variance element of the readouts.
+bool pmFPAReadVariance(pmFPA *fpa,        // FPA to read into
+                       psFits *fits,      // FITS file from which to read
+                       pmConfig *config   // Configuration
+    );
 
 /// Read cell headers
@@ -226 +226 @@
 /// also read and included in the cell analysis metadata under "name.HEADER".
 int pmCellReadTable(pmCell *cell,       ///< Cell for which to read table
-                    psFits *fits,       ///< FITS file from which the table
+                    psFits *fits,       ///< FITS file from which to read the table
                     const char *name    ///< Specifies the extension name, and target in the analysis metadata
                    );
@@ -233 +233 @@
 ///
 /// Iterates over component cells, calling pmCellReadTable.
-int pmChipReadTable(pmChip *chip,       ///< Cell for which to read table
-                    psFits *fits,       ///< FITS file from which the table
+int pmChipReadTable(pmChip *chip,       ///< Chip for which to read table
+                    psFits *fits,       ///< FITS file from which to read the table
                     const char *name    ///< Specifies the extension name, and target in the analysis metadata
                    );
@@ -241 +241 @@
 ///
 /// Iterates over component chips, calling pmChipReadTable.
-int pmFPAReadTable(pmFPA *fpa,          ///< Cell for which to read table
-                   psFits *fits,        ///< FITS file from which the table
+int pmFPAReadTable(pmFPA *fpa,          ///< FPA for which to read table
+                   psFits *fits,        ///< FITS file from which to read the table
                    const char *name     ///< Specifies the extension name, and target in the analysis metadata
                   );
+
+/// Read covariance matrices for a cell
+bool pmCellReadCovariance(pmCell *cell, ///< Cell for which to read covariance matrices
+                          psFits *fits  ///< FITS file from which to read
+    );
+
+/// Read covariance matrices for a chip
+bool pmChipReadCovariance(pmChip *chip, ///< Chip for which to read covariance matrices
+                          psFits *fits  ///< FITS file from which to read
+    );
+
+/// Read covariance matrices for a cell
+bool pmFPAReadCovariance(pmFPA *fpa,    ///< FPA for which to read covariance matrices
+                         psFits *fits   ///< FITS file from which to read
+    );
+
+
+
 /// @}
 #endif

trunk/psModules/src/camera/pmFPAWrite.c

@@ -9 +9 @@
 
 #include "pmConfig.h"
+#include "pmConfigMask.h"
 #include "pmHDU.h"
 #include "pmFPA.h"
@@ -29 +30 @@
     FPA_WRITE_TYPE_IMAGE,               // Write image
     FPA_WRITE_TYPE_MASK,                // Write mask
-    FPA_WRITE_TYPE_WEIGHT               // Write weight map
+    FPA_WRITE_TYPE_VARIANCE             // Write variance map
 } fpaWriteType;
 
@@ -46 +47 @@
     case FPA_WRITE_TYPE_MASK:
         return &hdu->masks;
-    case FPA_WRITE_TYPE_WEIGHT:
-        return &hdu->weights;
+    case FPA_WRITE_TYPE_VARIANCE:
+        return &hdu->variances;
     default:
         psAbort("Unknown write type: %x\n", type);
@@ -66 +67 @@
     case FPA_WRITE_TYPE_MASK:
         return pmHDUWriteMask(hdu, fits, config);
-    case FPA_WRITE_TYPE_WEIGHT:
-        return pmHDUWriteWeight(hdu, fits, config);
+    case FPA_WRITE_TYPE_VARIANCE:
+        return pmHDUWriteVariance(hdu, fits, config);
     default:
         psAbort("Unknown write type: %x\n", type);
@@ -74 +75 @@
 }
 
-// Write a cell image/mask/weight
+// Indicate whether a covariance matrix is defined
+static bool readoutSearchCovariances(pmReadout *ro)
+{
+    return ro->covariance ? true : false;
+}
+
+// Search for a covariance matrix
+#define SEARCH_COVARIANCES(NAME, PARENT, CHILD, CHILDREN, TESTFUNC) \
+static bool NAME(PARENT *parent) \
+{ \
+    if (!parent || !parent->CHILDREN) { \
+        return false; \
+    } \
+    psArray *children = parent->CHILDREN; /* Array of children */ \
+    for (int i = 0; i < children->n; i++) { \
+        CHILD *child = children->data[i]; /* Child of interest */ \
+        if (child && TESTFUNC(child)) { \
+            return true; \
+        } \
+    } \
+    return false; \
+}
+
+SEARCH_COVARIANCES(cellSearchCovariances, pmCell, pmReadout, readouts, readoutSearchCovariances);
+SEARCH_COVARIANCES(chipSearchCovariances, pmChip, pmCell,    cells,    cellSearchCovariances);
+SEARCH_COVARIANCES(fpaSearchCovariances,  pmFPA,  pmChip,    chips,    chipSearchCovariances);
+
+// 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
+                              pmConfig *config // Configuration
+                              )
+{
+    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: {
+          bool covar = false;           // Are covariances present?
+          if ((cell && cellSearchCovariances(cell)) ||
+              (!cell && ((chip && chipSearchCovariances(chip)) ||
+                         (!chip && fpa && fpaSearchCovariances(fpa))))) {
+              covar = true;
+          }
+
+          pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // Header being written
+          psMetadataAddBool(hdu->header, PS_LIST_TAIL, PM_HDU_COVARIANCE_KEYWORD, PS_META_REPLACE,
+                            "Is a covariance matrix present?", covar);
+          break;
+      }
+      default:
+        break;
+    }
+
+    return true;
+}
+
+
+// Write a cell image/mask/variance
 static bool cellWrite(pmCell *cell,     // Cell to write
                       psFits *fits,     // FITS file to which to write
@@ -97 +163 @@
     // 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 +172 @@
 
     if (writeBlank || writeImage) {
-
         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 (!writeUpdateHeader(NULL, NULL, cell, type, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
             return false;
         }
@@ -123 +192 @@
 }
 
-// Write a chip image/mask/weight
+// Write a chip image/mask/variance
 static bool chipWrite(pmChip *chip,     // Chip to write
                       psFits *fits,     // FITS file to which to write
@@ -162 +231 @@
                 return false;
             }
+
+            if (!writeUpdateHeader(NULL, chip, NULL, type, config)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
+                return false;
+            }
+
             if (!appropriateWriteFunc(hdu, fits, config, type)) {
                 psError(PS_ERR_IO, false, "Unable to write HDU for chip.\n");
@@ -186 +261 @@
 
 
-// Write an FPA image/mask/weight
+// Write an FPA image/mask/variance
 static bool fpaWrite(pmFPA *fpa,        // FPA to write
                      psFits *fits,      // FITS file to which to write
@@ -223 +298 @@
             if (!pmConceptsWriteFPA(fpa, source, true, config)) {
                 psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n");
+                return false;
+            }
+            if (!writeUpdateHeader(fpa, NULL, NULL, type, config)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to update header for writing");
                 return false;
             }
@@ -433 +512 @@
 
 
-bool pmCellWriteWeight(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
+bool pmCellWriteVariance(pmCell *cell, psFits *fits, pmConfig *config, bool blank)
 {
     PS_ASSERT_PTR_NON_NULL(cell, false);
     PS_ASSERT_PTR_NON_NULL(fits, false);
-    return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_WEIGHT);
-}
-
-bool pmChipWriteWeight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
+    if (!cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_VARIANCE)) {
+        return false;
+    }
+    if (!pmCellWriteCovariance(fits, cell)) {
+        return false;
+    }
+    return true;
+}
+
+bool pmChipWriteVariance(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse)
 {
     PS_ASSERT_PTR_NON_NULL(chip, false);
     PS_ASSERT_PTR_NON_NULL(fits, false);
-    return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT);
-}
-
-bool pmFPAWriteWeight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
+    if (!chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) {
+        return false;
+    }
+    if (!pmChipWriteCovariance(fits, chip)) {
+        return false;
+    }
+    return true;
+}
+
+bool pmFPAWriteVariance(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse)
 {
     PS_ASSERT_PTR_NON_NULL(fpa, false);
     PS_ASSERT_PTR_NON_NULL(fits, false);
-    return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT);
+    if (!fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_VARIANCE)) {
+        return false;
+    }
+    if (!pmFPAWriteCovariance(fits, fpa)) {
+        return false;
+    }
+    return true;
 }
 
@@ -526 +623 @@
     return numWrite;
 }
+
+bool pmCellWriteCovariance(psFits *fits, const pmCell *cell)
+{
+    PS_ASSERT_PTR_NON_NULL(cell, false);
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+
+    int numCovar = 0;
+    psArray *readouts = cell->readouts; // Array of readouts
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+        if (readout && readout->covariance) {
+            numCovar++;
+        }
+    }
+    if (numCovar == 0) {
+        return true;
+    }
+    if (numCovar != readouts->n) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Number of covariances (%d) doesn't match number of readouts (%ld)",
+                numCovar, readouts->n);
+        return false;
+    }
+
+    // Check size of covariances
+    int xMinCovar = INT_MAX, xMaxCovar = INT_MIN, yMinCovar = INT_MAX, yMaxCovar = INT_MIN; // Size
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+        psAssert(readout, "Should be defined.");
+        psKernel *covar = readout->covariance; // Covariance matrix
+        psAssert(covar, "Should be defined.");
+        xMinCovar = PS_MIN(xMinCovar, covar->xMin);
+        xMaxCovar = PS_MAX(xMaxCovar, covar->xMax);
+        yMinCovar = PS_MIN(yMinCovar, covar->yMin);
+        yMaxCovar = PS_MAX(yMaxCovar, covar->yMax);
+    }
+
+    // Correct covariances to common size
+    psArray *images = psArrayAlloc(numCovar); // Array of images
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+        psAssert(readout, "Should be defined.");
+        psKernel *covar = readout->covariance; // Covariance matrix
+        psAssert(covar, "Should be defined.");
+        int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax;// Size
+        if (xMin == xMinCovar && xMax == xMaxCovar && yMin == yMinCovar && yMax == yMaxCovar) {
+            images->data[i] = psMemIncrRefCounter(covar->image);
+        } else {
+            psImage *new = psImageAlloc(xMaxCovar - xMinCovar + 1, yMaxCovar - yMinCovar + 1, PS_TYPE_F32);
+            psImageInit(new, 0);
+            psImageOverlaySection(new, covar->image, xMinCovar - xMin, yMinCovar - yMin, "=");
+            images->data[i] = new;
+        }
+    }
+
+    // Determine extension name
+    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);
+
+    // Generate header
+    pmHDU *hdu = pmHDUFromCell(cell);   // HDU for cell
+    psMetadata *header = psMetadataCopy(NULL, hdu->header); // Header to write
+    psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.X", PS_META_REPLACE,
+                     "Centre of covariance matrix in x", -xMinCovar);
+    psMetadataAddS32(header, PS_LIST_TAIL, "COVARIANCE.CENTRE.Y", PS_META_REPLACE,
+                     "Centre of covariance matrix in y", -yMinCovar);
+
+    // Write images
+    if (!psFitsWriteImageCube(fits, header, images, extname)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to write covariances from chip %s, cell %s to extension %s",
+                chipName, cellName, extname);
+        psFree(extname);
+        psFree(header);
+        psFree(images);
+        return 0;
+    }
+    psFree(extname);
+    psFree(header);
+    psFree(images);
+
+    return true;
+}
+
+
+bool pmChipWriteCovariance(psFits *fits, const pmChip *chip)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_PTR_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 (!pmCellWriteCovariance(fits, cell)) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+bool pmFPAWriteCovariance(psFits *fits, const pmFPA *fpa)
+{
+    PS_ASSERT_PTR_NON_NULL(fpa, false);
+    PS_ASSERT_PTR_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 (!pmChipWriteCovariance(fits, chip)) {
+            return false;
+        }
+    }
+
+    return true;
+}

trunk/psModules/src/camera/pmFPAWrite.h

@@ -4 +4 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-04 02:39:36 $
+ * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
+ *
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -24 +25 @@
                         psFits *fits,   ///<  FITS file to which to write
                         int z           ///<  Image plane to write
-                       );
+    );
 
 /// Write a cell to a FITS file
@@ -36 +37 @@
                  pmConfig *config,      ///<  Configuration
                  bool blank             ///<  Write a blank PHU?
-                );
+    );
 
 /// Write a chip to a FITS file
@@ -49 +50 @@
                  bool blank,            ///<  Write a blank PHU?
                  bool recurse           ///<  Recurse to lower levels?
-                );
+    );
 
 /// Write an FPA to a FITS file
@@ -62 +63 @@
                 bool blank,             ///<  Write a blank PHU?
                 bool recurse            ///<  Recurse to lower levels?
-               );
+    );
 
 /// Write a cell mask to a FITS file
@@ -74 +75 @@
                      pmConfig *config,  ///<  Configuration
                      bool blank         ///<  Write a blank PHU?
-                    );
+    );
 
 /// Write a chip mask to a FITS file
@@ -88 +89 @@
                      bool blank,        ///<  Write a blank PHU?
                      bool recurse       ///<  Recurse to lower levels?
-                    );
+    );
 
 /// Write an FPA mask to a FITS file
@@ -102 +103 @@
                     bool blank,         ///<  Write a blank PHU?
                     bool recurse        ///<  Recurse to lower levels?
-                   );
+    );
 
-/// Write a cell weight to a FITS file
+/// Write a cell variance to a FITS file
 ///
-/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
+/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
 /// header) is written only if specifically requested.  Writes the concepts to the various locations, and then
-/// the HDU weight to the FITS file.  This function should be called at the beginning of the output cell loop
-/// with blank=true in order to produce the correct file structure.
-bool pmCellWriteWeight(pmCell *cell,    ///<  Cell to write
-                       psFits *fits,    ///<  FITS file to which to write
-                       pmConfig *config, ///<  Configuration
-                       bool blank       ///<  Write a blank PHU?
-                      );
+/// the HDU variance to the FITS file.  This function should be called at the beginning of the output cell
+/// loop with blank=true in order to produce the correct file structure.
+bool pmCellWriteVariance(pmCell *cell,    ///<  Cell to write
+                         psFits *fits,    ///<  FITS file to which to write
+                         pmConfig *config, ///<  Configuration
+                         bool blank       ///<  Write a blank PHU?
+    );
 
-/// Write a chip weight to a FITS file
+/// Write a chip variance to a FITS file
 ///
-/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
+/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
 /// header) is written only if specifically requested.  Writes the concepts to the various locations, and then
-/// the HDU weight to the FITS file, optionally recursing to lower levels.  This function should be called at
-/// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct
+/// the HDU variance to the FITS file, optionally recursing to lower levels.  This function should be called
+/// at the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct
 /// file structure.
-bool pmChipWriteWeight(pmChip *chip,    ///<  Chip to write
-                       psFits *fits,    ///<  FITS file to which to write
-                       pmConfig *config, ///<  Configuration
-                       bool blank,      ///<  Write a blank PHU?
-                       bool recurse     ///<  Recurse to lower levels?
-                      );
+bool pmChipWriteVariance(pmChip *chip,    ///<  Chip to write
+                         psFits *fits,    ///<  FITS file to which to write
+                         pmConfig *config, ///<  Configuration
+                         bool blank,      ///<  Write a blank PHU?
+                         bool recurse     ///<  Recurse to lower levels?
+    );
 
-/// Write an FPA weight to a FITS file
+/// Write an FPA variance to a FITS file
 ///
-/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required.  A blank (i.e., image-less
+/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU variance pixels if required.  A blank (i.e., image-less
 /// header) is written only if specifically requested.  Writes the concepts to the various locations, and then
-/// the HDU weight to the FITS file, optionally recursing to lower levels.  This function should be called at
-/// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct
+/// the HDU variance to the FITS file, optionally recursing to lower levels.  This function should be called
+/// at the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct
 /// file structure.
-bool pmFPAWriteWeight(pmFPA *fpa,       ///<  FPA to write
-                      psFits *fits,     ///<  FITS file to which to write
-                      pmConfig *config, ///<  Configuration
-                      bool blank,       ///<  Write a blank PHU?
-                      bool recurse      ///<  Recurse to lower levels?
-                     );
+bool pmFPAWriteVariance(pmFPA *fpa,       ///<  FPA to write
+                        psFits *fits,     ///<  FITS file to which to write
+                        pmConfig *config, ///<  Configuration
+                        bool blank,       ///<  Write a blank PHU?
+                        bool recurse      ///<  Recurse to lower levels?
+    );
 
 
@@ -153 +154 @@
                      const pmCell *cell, ///< Cell containing FITS table in the analysis metadata
                      const char *name   ///< Name for the table data, and the extension name
-                    );
+    );
 
 int pmChipWriteTable(psFits *fits,      ///< FITS file to which to write
                      const pmChip *chip, ///< Chip containing cells with tables to write
                      const char *name   ///< Name for the table data, and the extension name
-                    );
+    );
 
 int pmFPAWriteTable(psFits *fits,       ///< FITS file to which to write
                     const pmFPA *fpa,   ///< FPA containing cells with tables to write
                     const char *name    ///< Name for the table data, and the extension name
-                   );
+    );
+
+/// Write covariance matrix to a FITS file
+///
+/// The covariance matrices for a cell are written to an independent extension, named after the chip and cell
+/// name.
+bool pmCellWriteCovariance(psFits *fits,///< FITS file to which to write
+                           const pmCell *cell ///< Cell for which to write covariance
+    );
+
+bool pmChipWriteCovariance(psFits *fits,///< FITS file to which to write
+                           const pmChip *chip ///< Chip for which to write covariance
+    );
+
+bool pmFPAWriteCovariance(psFits *fits,///< FITS file to which to write
+                          const pmFPA *fpa ///< FPA for which to write covariance
+    );
 
 // Update the header before writing to be consistent

trunk/psModules/src/camera/pmFPAfile.c

@@ -476 +476 @@
         return PM_FPA_FILE_MASK;
     }
-    if (!strcasecmp(type, "WEIGHT"))     {
-        return PM_FPA_FILE_WEIGHT;
+    if (!strcasecmp(type, "VARIANCE"))     {
+        return PM_FPA_FILE_VARIANCE;
     }
     if (!strcasecmp(type, "FRINGE")) {
@@ -530 +530 @@
       case PM_FPA_FILE_MASK:
         return ("MASK");
-      case PM_FPA_FILE_WEIGHT:
-        return ("WEIGHT");
+      case PM_FPA_FILE_VARIANCE:
+        return ("VARIANCE");
       case PM_FPA_FILE_FRINGE:
         return ("FRINGE");

trunk/psModules/src/camera/pmFPAfile.h

@@ -4 +4 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-04 02:39:36 $
+ * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
+ *
  * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
  */
@@ -37 +38 @@
     PM_FPA_FILE_IMAGE,
     PM_FPA_FILE_MASK,
-    PM_FPA_FILE_WEIGHT,
+    PM_FPA_FILE_VARIANCE,
     PM_FPA_FILE_FRINGE,
     PM_FPA_FILE_DARK,

trunk/psModules/src/camera/pmFPAfileFitsIO.c

@@ -148 +148 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_HEADER:
       case PM_FPA_FILE_FRINGE:
@@ -329 +329 @@
 }
 
-bool pmFPAviewReadFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
-{
-    PS_ASSERT_PTR_NON_NULL(view, false);
-    PS_ASSERT_PTR_NON_NULL(file, false);
-    return fpaViewReadFitsImage(view, file, config, pmFPAReadWeight, pmChipReadWeight, pmCellReadWeight);
+bool pmFPAviewReadFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+    return fpaViewReadFitsImage(view, file, config, pmFPAReadVariance, pmChipReadVariance,
+                                pmCellReadVariance);
 }
 
@@ -347 +348 @@
     PS_ASSERT_PTR_NON_NULL(view, false);
     PS_ASSERT_PTR_NON_NULL(file, false);
-    return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, pmCellReadHeaderSet);
+    return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet,
+                                pmCellReadHeaderSet);
 }
 
@@ -432 +434 @@
 }
 
-bool pmFPAviewWriteFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
-{
-    PS_ASSERT_PTR_NON_NULL(view, false);
-    PS_ASSERT_PTR_NON_NULL(file, false);
-    return fpaViewWriteFitsImage(view, file, config, pmFPAWriteWeight, pmChipWriteWeight, pmCellWriteWeight);
+bool pmFPAviewWriteFitsVariance(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+    return fpaViewWriteFitsImage(view, file, config, pmFPAWriteVariance, pmChipWriteVariance,
+                                 pmCellWriteVariance);
 }
 

trunk/psModules/src/camera/pmFPAfileFitsIO.h

@@ -5 +5 @@
  * @author PAP, IfA
  *
- * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-07-17 22:38:15 $
+ * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
  * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
  */
@@ -27 +27 @@
                             pmConfig *config
                           );
-/// Read a weight map into the current view
-bool pmFPAviewReadFitsWeight(const pmFPAview *view,  ///< View specifying level of interest
+/// Read a variance map into the current view
+bool pmFPAviewReadFitsVariance(const pmFPAview *view,  ///< View specifying level of interest
                              pmFPAfile *file, ///< FPA file into which to read
                             pmConfig *config
@@ -57 +57 @@
                            );
 
-/// Write the weight map for the specified view
-bool pmFPAviewWriteFitsWeight(const pmFPAview *view, ///< View specifying level of interest
+/// Write the variance map for the specified view
+bool pmFPAviewWriteFitsVariance(const pmFPAview *view, ///< View specifying level of interest
                               pmFPAfile *file, ///< FPA file to write
                               pmConfig *config ///< Configuration

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -173 +173 @@
         status = pmFPAviewReadFitsMask(view, file, config);
         break;
-      case PM_FPA_FILE_WEIGHT:
-        status = pmFPAviewReadFitsWeight(view, file, config);
+      case PM_FPA_FILE_VARIANCE:
+        status = pmFPAviewReadFitsVariance(view, file, config);
         break;
       case PM_FPA_FILE_HEADER:
@@ -262 +262 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_DARK: {
@@ -427 +427 @@
         status = pmFPAviewWriteFitsMask(view, file, config);
         break;
-      case PM_FPA_FILE_WEIGHT:
-        status = pmFPAviewWriteFitsWeight(view, file, config);
+      case PM_FPA_FILE_VARIANCE:
+        status = pmFPAviewWriteFitsVariance(view, file, config);
         break;
       case PM_FPA_FILE_HEADER:
@@ -524 +524 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_HEADER:
       case PM_FPA_FILE_FRINGE:
@@ -589 +589 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_HEADER:
       case PM_FPA_FILE_FRINGE:
@@ -730 +730 @@
     psString tmpName = pmConfigConvertFilename (file->filename, config, create, false);
     if (!tmpName) {
-        psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename);
-        return false;
+        psError(PS_ERR_IO, false, "failed to determine real name from template for %s\n", file->filename);
+        return false;
     }
     psFree (file->filename);
@@ -740 +740 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_HEADER:
       case PM_FPA_FILE_FRINGE:
@@ -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
           // XXX can I actually reach this with camera not specified??
-          psMetadata *camera = NULL;
-          psString formatName = NULL;
-          psString cameraName = NULL;
-          file->format = pmConfigCameraFormatFromHeader (&camera, &cameraName, &formatName, config, phu, true);
+          psMetadata *camera = NULL;
+          psString formatName = NULL;
+          psString cameraName = NULL;
+          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);
@@ -794 +785 @@
           psFree(phu);
 
-          pmFPA *newFPA = pmFPAConstruct (camera, formatName);
-          if (!newFPA) {
-              psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName);
-              psFree(camera);
-              psFree(formatName);
-              return NULL;
-          }
-          psFree(camera);
-          psFree(formatName);
-          psFree(cameraName);
-
-          // XXX this is really dangerous...
-          psFree (file->fpa);
-          file->fpa = newFPA;
+          pmFPA *newFPA = pmFPAConstruct (camera, formatName);
+          if (!newFPA) {
+              psError(PS_ERR_IO, false, "Failed to construct FPA from %s for %s", file->filename, formatName);
+              psFree(camera);
+              psFree(formatName);
+              return NULL;
+          }
+          psFree(camera);
+          psFree(formatName);
+          psFree(cameraName);
+
+          // XXX this is really dangerous...
+          psFree (file->fpa);
+          file->fpa = newFPA;
         }
 
@@ -938 +929 @@
       case PM_FPA_FILE_IMAGE:
       case PM_FPA_FILE_MASK:
-      case PM_FPA_FILE_WEIGHT:
+      case PM_FPA_FILE_VARIANCE:
       case PM_FPA_FILE_DARK:
       case PM_FPA_FILE_FRINGE:

trunk/psModules/src/camera/pmHDU.c

@@ -45 +45 @@
     psFree(hdu->header);
     psFree(hdu->images);
-    psFree(hdu->weights);
+    psFree(hdu->variances);
     psFree(hdu->masks);
 }
@@ -69 +69 @@
     hdu->header  = NULL;
     hdu->images  = NULL;
-    hdu->weights = NULL;
+    hdu->variances = NULL;
     hdu->masks   = NULL;
 
@@ -155 +155 @@
 }
 
-bool pmHDUReadWeight(pmHDU *hdu, psFits *fits)
-{
-    PS_ASSERT_PTR_NON_NULL(hdu, false);
-    PS_ASSERT_PTR_NON_NULL(fits, false);
-
-    return hduRead(hdu, &hdu->weights, fits);
+bool pmHDUReadVariance(pmHDU *hdu, psFits *fits)
+{
+    PS_ASSERT_PTR_NON_NULL(hdu, false);
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+
+    return hduRead(hdu, &hdu->variances, fits);
 }
 
@@ -235 +235 @@
 }
 
-bool pmHDUWriteWeight(pmHDU *hdu, psFits *fits, const pmConfig *config)
+bool pmHDUWriteVariance(pmHDU *hdu, psFits *fits, const pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(hdu, false);
@@ -241 +241 @@
 
     psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask
-    return hduWrite(hdu, hdu->weights, hdu->masks, maskVal, fits);
+    return hduWrite(hdu, hdu->variances, hdu->masks, maskVal, fits);
 }
 

trunk/psModules/src/camera/pmHDU.h

@@ -4 +4 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-04 02:39:36 $
+ * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-06 02:31:24 $
+ *
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */
@@ -18 +19 @@
 /// @{
 
+#define PM_HDU_COVARIANCE_KEYWORD "PS_COVAR" // FITS keyword to indicate presence of a covariance matrix
+
 /// An instance of the FITS Header Data Unit
 ///
-/// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weight data, but these are
+/// 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 *weights;                   ///< 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
+} pmHDU;
 
 
@@ -60 +61 @@
                   );
 
-/// Read the HDU header and weight map
+/// Read the HDU header and variance map
 ///
 /// Moves to the appropriate extension
-bool pmHDUReadWeight(pmHDU *hdu,        ///< HDU to read
-                     psFits *fits       ///< FITS file to read from
+bool pmHDUReadVariance(pmHDU *hdu,        ///< HDU to read
+                       psFits *fits       ///< FITS file to read from
     );
 
@@ -79 +80 @@
     );
 
-/// Write the HDU header and weight map
-bool pmHDUWriteWeight(pmHDU *hdu,       ///< HDU to write
-                      psFits *fits,     ///< FITS file to write to
-                      const pmConfig *config  ///< Configuration
+/// Write the HDU header and variance map
+bool pmHDUWriteVariance(pmHDU *hdu,       ///< HDU to write
+                        psFits *fits,     ///< FITS file to write to
+                        const pmConfig *config  ///< Configuration
     );
 

trunk/psModules/src/camera/pmHDUGenerate.c

@@ -274 +274 @@
         pmReadout *readout = cell->readouts->data[0]; // The first readout, as representative
         // The proper image, used to get the size
-        psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout->weight);
+        psImage *image = readout->image ? readout->image : (readout->mask ? readout->mask : readout->variance);
         if (!image) {
             continue;
@@ -283 +283 @@
                      image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows);
         }
-        if (readout->weight &&
-                (readout->weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) {
-            psWarning("Image and weight have different sizes (%dx%d vs %dx%d)!\n",
-                     image->numCols, image->numRows, readout->weight->numCols, readout->weight->numRows);
+        if (readout->variance &&
+                (readout->variance->numCols != image->numCols || readout->variance->numRows != image->numRows)) {
+            psWarning("Image and variance have different sizes (%dx%d vs %dx%d)!\n",
+                     image->numCols, image->numRows, readout->variance->numCols, readout->variance->numRows);
         }
         // New reference
@@ -355 +355 @@
     psElemType imageType = 0;           // Type of readout images
     psElemType maskType = 0;            // Type of readout masks
-    psElemType weightType = 0;          // Type of readout weights
+    psElemType varianceType = 0;        // Type of readout variances
     {
         psListIterator *iter = psListIteratorAlloc(cells, PS_LIST_HEAD, false); // Iterator for cells
@@ -381 +381 @@
                     maskType = checkTypes(maskType, readout->mask->type.type);
                 }
-                if (!hdu->weights && readout->weight) {
-                    weightType = checkTypes(weightType, readout->weight->type.type);
+                if (!hdu->variances && readout->variance) {
+                    varianceType = checkTypes(varianceType, readout->variance->type.type);
                 }
             }
@@ -388 +388 @@
         psFree(iter);
     }
-    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && weightType == 0)) {
+    if (numReadouts == 0 || (imageType == 0 && maskType == 0 && varianceType == 0)) {
         // Nothing from which to create an HDU
         psFree(cells);
@@ -421 +421 @@
         }
     }
-    if (weightType) {
-        hdu->weights = psArrayAlloc(numReadouts);
+    if (varianceType) {
+        hdu->variances = psArrayAlloc(numReadouts);
         for (int i = 0; i < numReadouts; i++) {
-            psImage *weight = psImageAlloc(xSize, ySize, weightType);
-            psImageInit(weight, 0.0);
-            hdu->weights->data[i] = weight;
+            psImage *variance = psImageAlloc(xSize, ySize, varianceType);
+            psImageInit(variance, 0.0);
+            hdu->variances->data[i] = variance;
         }
     }
@@ -448 +448 @@
 
             psArray *readouts = cell->readouts; // Array of readouts
+
             psArray *hduImages = hdu->images; // Array of images in the HDU
             psArray *hduMasks = hdu->masks; // Array of masks in the HDU
-            psArray *hduWeights = hdu->weights; // Array of weights in the HDU
+            psArray *hduVariances = hdu->variances; // Array of variances in the HDU
             for (int i = 0; i < readouts->n; i++) {
                 pmReadout *readout = readouts->data[i]; // The readout of interest
@@ -467 +468 @@
                     readout->mask = new;
                 }
-                if (readout->weight) {
-                    psImage *new = pasteImage(hduWeights->data[i], readout->weight, trimsec);
-                    psFree(readout->weight);
-                    readout->weight = new;
+                if (readout->variance) {
+                    psImage *new = pasteImage(hduVariances->data[i], readout->variance, trimsec);
+                    psFree(readout->variance);
+                    readout->variance = new;
                 }
 
@@ -504 +505 @@
 
 // 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 +513 @@
     }
     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 +582 @@
         return false;
     }
-    if (hdu->images && hdu->masks && hdu->weights) {
+    if (hdu->images && hdu->masks && hdu->variances) {
         // It's already here!
         return true;
@@ -631 +632 @@
         return generateForCells(chip);
     }
-    if (hdu->images && hdu->masks && hdu->weights) {
+    if (hdu->images && hdu->masks && hdu->variances) {
         // It's already here!
         return true;
@@ -680 +681 @@
         return generateForChips(fpa);
     }
-    if (hdu->images && hdu->masks && hdu->weights) {
+    if (hdu->images && hdu->masks && hdu->variances) {
         // It's already here!
         return true;

trunk/psModules/src/camera/pmHDUUtils.c

@@ -18 +18 @@
 
     for (int i = 0; i < fpa->chips->n; i++) {
-        pmChip *chip = fpa->chips->data[i];
-        if (!chip) continue;
-        if (chip->hdu) return chip->hdu;
-        if (!chip->cells) continue;
-        for (int j = 0; j < chip->cells->n; j++) {
-            pmCell *cell = chip->cells->data[j];
-            if (!cell) continue;
-            if (cell->hdu) return cell->hdu;
-        }
+        pmChip *chip = fpa->chips->data[i];
+        if (!chip) continue;
+        if (chip->hdu) return chip->hdu;
+        if (!chip->cells) continue;
+        for (int j = 0; j < chip->cells->n; j++) {
+            pmCell *cell = chip->cells->data[j];
+            if (!cell) continue;
+            if (cell->hdu) return cell->hdu;
+        }
     }
     return NULL;
@@ -158 +158 @@
 
     INDENT(fd, level + 1);
-    if (hdu->weights) {
-        fprintf(fd, "Weights:\n");
-        for (long i = 0; i < hdu->weights->n; i++) {
-            psImage *weight = hdu->weights->data[i]; // Weight image of interest
+    if (hdu->variances) {
+        fprintf(fd, "Variances:\n");
+        for (long i = 0; i < hdu->variances->n; i++) {
+            psImage *variance = hdu->variances->data[i]; // Variance image of interest
             INDENT(fd, level + 2);
-            fprintf(fd, "%ld: %dx%d\n", i, weight->numCols, weight->numRows);
+            fprintf(fd, "%ld: %dx%d\n", i, variance->numCols, variance->numRows);
         }
     } else {
-        fprintf(fd, "NO weights.\n");
+        fprintf(fd, "NO variances.\n");
     }
 

trunk/psModules/src/camera/pmReadoutStack.c

@@ -12 +12 @@
 // XXX should it be an error for the image to exist?
 psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image
-                                   const char *name, // Name of image in analysis metadata
-                                   int numCols, int numRows, // Expected size of image
-                                   psElemType type, // Expected type of image
-                                   double init // Initial value
+                                   const char *name, // Name of image in analysis metadata
+                                   int numCols, int numRows, // Expected size of image
+                                   psElemType type, // Expected type of image
+                                   double init // Initial value
     )
 {
@@ -24 +24 @@
 
     if (!psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image)) {
-        psAbort ("analysis image already exists");
+        psAbort ("analysis image already exists");
     }
     psImageInit(image, init);
@@ -35 +35 @@
 // XXX not sure why this should call psMemIncrRefCounter
 psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image
-                                   const char *name       // Name of image in analysis metadata
+                                   const char *name       // Name of image in analysis metadata
     )
 {
@@ -78 +78 @@
 
 // XXX for the moment, use col0, row0, numCols, numRows supplied from the outside
-bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool weight, psImageMaskType blank)
+bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool variance, psImageMaskType blank)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -92 +92 @@
 
     if (mask) {
-        // XXX is this an error?
+        // XXX is this an error?
         if (readout->mask) return false;
-        readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-        psImageInit(readout->mask, blank);
-    }
-
-    if (weight) {
-        // XXX is this an error?
-        if (readout->weight) return false;
-        readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        psImageInit(readout->weight, NAN);
+        readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+        psImageInit(readout->mask, blank);
+    }
+
+    if (variance) {
+        // XXX is this an error?
+        if (readout->variance) return false;
+        readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        psImageInit(readout->variance, NAN);
     }
 
@@ -162 +162 @@
 
 bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows,
-                         int numCols, int numRows, bool mask, bool weight,
+                         int numCols, int numRows, bool mask, bool variance,
                          psImageMaskType blank)
 {
@@ -203 +203 @@
     }
 
-    if (weight) {
-        if (!readout->weight) {
-            readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            psImageInit(readout->weight, NAN);
-        }
-        if (readout->weight->numCols < numCols || readout->weight->numRows < numRows) {
-            psImage *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            psImageInit(newWeight, NAN);
-            psImageOverlaySection(newWeight, readout->weight, readout->col0, readout->row0, "=");
-            psFree(readout->weight);
-            readout->weight = newWeight;
+    if (variance) {
+        if (!readout->variance) {
+            readout->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            psImageInit(readout->variance, NAN);
+        }
+        if (readout->variance->numCols < numCols || readout->variance->numRows < numRows) {
+            psImage *newVariance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            psImageInit(newVariance, NAN);
+            psImageOverlaySection(newVariance, readout->variance, readout->col0, readout->row0, "=");
+            psFree(readout->variance);
+            readout->variance = newVariance;
         }
     }