Index: /branches/rel10_ifa/psLib/src/fits/psFitsImage.c
===================================================================
--- /branches/rel10_ifa/psLib/src/fits/psFitsImage.c	(revision 6511)
+++ /branches/rel10_ifa/psLib/src/fits/psFitsImage.c	(revision 6512)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-14 22:18:30 $
+ *  @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-04 01:23:17 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -119,8 +119,8 @@
 }
 
-psImage* psFitsReadImage(psImage* output, // a psImage to recycle.
-                         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
+psImage* psFitsReadImage(//psImage* output, // a psImage to recycle.
+    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
 {
     psS32 status = 0;           /* CFITSIO file vars */
@@ -139,5 +139,4 @@
         psError(PS_ERR_BAD_PARAMETER_NULL, true,
                 PS_ERRORTEXT_psFits_NULL);
-        psFree(output);
         return NULL;
     }
@@ -165,5 +164,4 @@
                 PS_ERRORTEXT_psFits_DATATYPE_UNKNOWN,
                 fitsErr);
-        psFree(output);
         return NULL;
     }
@@ -175,5 +173,4 @@
                 PS_ERRORTEXT_psFits_IMAGE_DIM_UNKNOWN,
                 fitsErr);
-        psFree(output);
         return NULL;
     }
@@ -184,5 +181,4 @@
                 PS_ERRORTEXT_psFits_IMAGE_DIMENSION_UNSUPPORTED,
                 nAxis);
-        psFree(output);
         return NULL;
     }
@@ -194,5 +190,4 @@
                 PS_ERRORTEXT_psFits_IMAGE_SIZE_UNKNOWN,
                 fitsErr);
-        psFree(output);
         return NULL;
     }
@@ -259,18 +254,8 @@
                 PS_ERRORTEXT_psFits_FITS_TYPE_UNSUPPORTED,
                 bitPix);
-        psFree(output);
-        return NULL;
-    }
-
-    output = psImageRecycle(output,
-                            lastPixel[0]-firstPixel[0]+1,
-                            lastPixel[1]-firstPixel[1]+1,
-                            datatype);
-
-    if (output == NULL) {
-        psError(PS_ERR_UNKNOWN, false,
-                "Failed to allocate a properly sized image.");
-        return false;
-    }
+        return NULL;
+    }
+
+    output = psImageAlloc(lastPixel[0]-firstPixel[0]+1, lastPixel[1]-firstPixel[1]+1, datatype);
 
     // n.b., this assumes contiguous image buffer
@@ -288,4 +273,54 @@
 
 }
+
+psArray *psFitsReadImageCube(const psFits *fits, psRegion region)
+{
+    int nAxis = 0;                      // Number of axes
+    long nAxes[3];                      // Number of pixels on each axis
+
+    if (fits == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                PS_ERRORTEXT_psFits_NULL);
+        return NULL;
+    }
+
+    // Some of this replicates what is in psFitsReadImage, so it's a little inefficient.  But it saves
+    // code replication, and should be sufficient for our needs.
+
+    if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) {
+        (void)fits_get_errstatus(status, fitsErr);
+        psError(PS_ERR_IO, true,
+                PS_ERRORTEXT_psFits_IMAGE_DIM_UNKNOWN,
+                fitsErr);
+        return NULL;
+    }
+
+    if (nAxis == 2) {
+        psArray *images = psArrayAlloc(1); // Single image plane
+        images->data[0] = psFitsReadImage(fits, region, 0);
+        return images;
+    }
+    if (nAxis == 3) {
+        if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) {
+            (void)fits_get_errstatus(status, fitsErr);
+            psError(PS_ERR_IO, true,
+                    PS_ERRORTEXT_psFits_IMAGE_SIZE_UNKNOWN,
+                    fitsErr);
+            return NULL;
+        }
+
+        psArray *images = psArrayAlloc(nAxes[2]); // Array of image planes
+        for (int i = 0; i < nAxes[2]; i++) {
+            images->data[i] = psFitsReadImage(fits, region, i);
+        }
+
+        return images;
+    }
+
+    // Bad dimensionality
+    psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_IMAGE_DIMENSION_UNSUPPORTED, nAxis);
+    return NULL;
+}
+
 
 bool psFitsWriteImage(psFits* fits,
@@ -340,4 +375,9 @@
     // write the header, if any.
     if (header != NULL) {
+        // Update with BITPIX, NAXIS[12]
+        psMetadataAddS32(header, PS_LIST_HEAD, "BITPIX", PS_META_REPLACE, "Bits per pixel", bitPix);
+        psMetadataAddS32(header, PS_LIST_HEAD, "NAXIS1", PS_META_REPLACE, "Bits per pixel", numCols);
+        psMetadataAddS32(header, PS_LIST_HEAD, "NAXIS2", PS_META_REPLACE, "Bits per pixel", numRows);
+
         psFitsWriteHeader(header, (psPtr)fits);
     }
@@ -376,4 +416,77 @@
 }
 
+bool psFitsWriteImageCube(psFits *fits, psMetadata *header, const psArray *input, const char *extname)
+{
+    if (fits == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                PS_ERRORTEXT_psFits_NULL);
+        return false;
+    }
+
+    if (input == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                PS_ERRORTEXT_psFits_IMAGE_NULL);
+        return false;
+    }
+
+    if (input->n == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psFits_ARRAY_EMPTY);
+        return falase;
+    }
+
+    if (input->n == 1) {
+        // The problem reduces to one already solved
+        return psFitsWriteImage(fits, header, input->data[0], extname);
+    }
+
+    // Check that all images are of the same size
+    {
+        psImage *image = input->data[0];// First image off the array
+        int numCols = image->numCols;   // Number of columns
+        int numRows = image->numRows;   // Number of rows
+        for (int i = 1; i < input->n; i++)
+        {
+            image = input->data[i];
+            if (image->numCols != numCols || image->numRows != numRows) {
+                psError(PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psFits_ARRAY_SIZE_DIFFER);
+                return false;
+            }
+        }
+    }
+
+    // Need to check the header to make sure NAXIS and NAXIS[1-3] are correct
+    psMetadata *headerCopy = NULL;      // Copy of header
+    if (header) {
+        headerCopy = psMemIncrRefCounter(header);
+    } else {
+        headerCopy = psMetadataAlloc();
+    }
+    bool update = psMetadataAddS32(headerCopy, PS_LIST_HEAD, "NAXIS", PS_META_REPLACE, "Dimensionality", 3) &&
+                  psMetadataAddS32(headerCopy, PS_LIST_HEAD, "NAXIS1", PS_META_REPLACE, "Number of columns", numCols) &&
+                  psMetadataAddS32(headerCopy, PS_LIST_HEAD, "NAXIS2", PS_META_REPLACE, "Number of rows", numRows) &&
+                  psMetadataAddS32(headerCopy, PS_LIST_HEAD, "NAXIS3", PS_META_REPLACE, "Number of image planes",
+                                   images->n);
+    if (! update) {
+        psError(PS_ERR_UNKNOWN, false, psFits_METADATA_ADD_FAILED, "NAXIS, NAXIS1, NAXIS2, NAXIS3");
+        psFree(headerCopy);
+        return false;
+    }
+
+    // Now we can safely write the images out
+    for (int i = 0; i < images->n; i++) {
+        if (! psFitsWriteImage(fits, headerCopy, images->data[i], i, extname)) {
+            psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psFits_WRITE_PLANE_FAILED, i);
+            psFree(headerCopy);
+            return false;
+        }
+    }
+
+    psFree(headerCopy);                 // Free, or drop reference
+
+    return true;
+}
+
+
+
 bool psFitsUpdateImage(psFits* fits,
                        const psImage* input,