Index: trunk/psLib/src/fits/psFitsImage.c
===================================================================
--- trunk/psLib/src/fits/psFitsImage.c	(revision 19270)
+++ trunk/psLib/src/fits/psFitsImage.c	(revision 19383)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-08-29 01:01:39 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-09-05 07:51:09 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -28,6 +28,6 @@
 #include "psLogMsg.h"
 #include "psTrace.h"
-#include "psVector.h"
 #include "psRandom.h"
+#include "psImage.h"
 #include "psImageStructManip.h"
 
@@ -249,5 +249,7 @@
                                           psFitsFloat *floatType, // Type of custom floating-point
                                           psFits *fits, // FITS file pointer
-                                          const psImage *image, // Current type
+                                          const psImage *image, // Image to convert
+                                          const psImage *mask, // Mask image, or NULL
+                                          psMaskType maskVal, // Value to mask
                                           psRandom *rng, // Random number generator
                                           bool newScaleZero // Determine a new BSCALE and BZERO?
@@ -285,5 +287,5 @@
         if (newScaleZero) {
             // Choose an appropriate BSCALE and BZERO
-            if (!psFitsScaleDetermine(bscale, bzero, blank, image, fits)) {
+            if (!psFitsScaleDetermine(bscale, bzero, blank, image, mask, maskVal, fits)) {
                 // We can't have the write dying for this reason --- try to save it somehow!
                 psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
@@ -404,5 +406,5 @@
 }
 
-psImage* psFitsReadImage(const psFits *fits, // the psFits object
+psImage *psFitsReadImage(const psFits *fits, // the psFits object
                          psRegion region, // the region in the FITS image to read
                          int z          // the z-plane in the FITS image cube to read
@@ -445,5 +447,5 @@
 }
 
-psImage* psFitsReadImageBuffer(psImage *outImage, // Output image buffer
+psImage *psFitsReadImageBuffer(psImage *outImage, // Output image buffer
                                const psFits *fits, // the psFits object
                                psRegion region, // the region in the FITS image to read
@@ -497,9 +499,13 @@
 }
 
-bool psFitsWriteImage(psFits* fits,
-                      psMetadata* header,
-                      const psImage* input,
-                      int numZPlanes,
-                      const char* extname)
+bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input,
+                      int numZPlanes, const char *extname)
+{
+    return psFitsWriteImageWithMask(fits, header, input, NULL, 0, numZPlanes, extname);
+}
+
+bool psFitsWriteImageWithMask(psFits *fits, psMetadata *header, const psImage *input,
+                              const psImage *mask, psMaskType maskVal, int numZPlanes,
+                              const char *extname)
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -509,13 +515,24 @@
 
     psFitsMoveLast(fits);
-    return psFitsInsertImage(fits,header,input,numZPlanes,extname,true);
-}
-
-bool psFitsInsertImage(psFits* fits, psMetadata* header, const psImage* image, int numZPlanes,
-                       const char* extname, bool after)
+    return psFitsInsertImageWithMask(fits, header, input, mask, maskVal, numZPlanes, extname, true);
+}
+
+bool psFitsInsertImage(psFits *fits, psMetadata *header, const psImage *image, int numZPlanes,
+                       const char *extname, bool after)
+{
+    return psFitsInsertImageWithMask(fits, header, image, NULL, 0, numZPlanes, extname, after);
+}
+
+bool psFitsInsertImageWithMask(psFits *fits, psMetadata *header, const psImage *image,
+                               const psImage *mask, psMaskType maskVal, int numZPlanes,
+                               const char *extname, bool after)
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
     PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_IMAGE_NON_NULL(image, false);
+    if (mask) {
+        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false);
+    }
 
     int numCols = image->numCols;       // Number of columns for image
@@ -529,5 +546,5 @@
     psFitsFloat floatType;              // Custom floating-point convention type
     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, image,
-                                                   NULL, true); // Image to write out
+                                                   mask, maskVal, NULL, true); // Image to write out
     if (!diskImage) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
@@ -557,6 +574,6 @@
     }
     if (cfitsioBzero != 0.0) {
-        // p_psFitsTypeToCfitsio and imageToDiskRepresentation must not clash!
-        psAssert(bzero == 0.0 && bscale == 1.0, "impossible");
+        psAssert(bzero == 0.0 && bscale == 1.0,
+                 "p_psFitsTypeToCfitsio and imageToDiskRepresentation must not clash!");
         bscale = 1.0;
         bzero = cfitsioBzero;
@@ -564,5 +581,6 @@
 
     psFitsOptions *options = fits->options; // FITS I/O options
-    psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0, "impossible");
+    psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0,
+             "Something's not consistent");
 
     int naxis = 3;                      // Number of axes
@@ -630,4 +648,6 @@
         // own quantisation, the correct keyword is always "BLANK" instead of "ZBLANK".
         fits_write_key_lng(fits->fd, "BLANK", blank, "Value for undefined pixels", &status);
+        // But it seems that cfitsio uses ZBLANK
+        fits_write_key_lng(fits->fd, "ZBLANK", blank, "Value for undefined pixels", &status);
         fits_set_imgnull(fits->fd, blank, &status);
         if (psFitsError(status, true, "Could not write BLANK header to file.")) {
@@ -681,9 +701,19 @@
 }
 
-bool psFitsUpdateImage(psFits* fits, const psImage* input, int x0, int y0, int z)
+bool psFitsUpdateImage(psFits *fits, const psImage *input, int x0, int y0, int z)
+{
+    return psFitsUpdateImageWithMask(fits, input, NULL, 0, x0, y0, z);
+}
+
+bool psFitsUpdateImageWithMask(psFits *fits, const psImage *input, const psImage *mask, psMaskType maskVal,
+                               int x0, int y0, int z)
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
     PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_IMAGE_NON_NULL(input, false);
+    if (mask) {
+        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, input, false);
+    }
 
     int status = 0;
@@ -709,5 +739,5 @@
     psFitsFloat floatType;              // Custom floating-point convention type
     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, input,
-                                                   NULL, false); // Image to write out
+                                                   mask, maskVal, NULL, false); // Image to write out
     if (!diskImage) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
@@ -737,11 +767,12 @@
     }
     if (cfitsioBzero != 0.0) {
-        // p_psFitsTypeToCfitsio and imageToDiskRepresentation must not clash!
-        psAssert(bzero == 0.0 && bscale == 1.0, "impossible");
+        psAssert(bzero == 0.0 && bscale == 1.0,
+                 "p_psFitsTypeToCfitsio and imageToDiskRepresentation must not clash!");
         bscale = 1.0;
         bzero = cfitsioBzero;
     }
     psFitsOptions *options = fits->options; // FITS I/O options
-    psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0, "impossible");
+    psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0,
+             "Something's not consistent");
 
     // Check to see if the HDU has the same datatype