Index: trunk/psLib/src/fits/psFitsImage.c
===================================================================
--- trunk/psLib/src/fits/psFitsImage.c	(revision 16095)
+++ trunk/psLib/src/fits/psFitsImage.c	(revision 16185)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-01-16 20:10:35 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-01-23 03:08:03 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -14,5 +14,5 @@
 
 #ifdef HAVE_CONFIG_H
-# include "config.h"
+#include "config.h"
 #endif
 
@@ -36,4 +36,5 @@
 #include "psFitsFloatFile.h"
 #include "psFitsHeader.h"
+#include "psFitsScale.h"
 
 #include "psMemory.h"
@@ -241,264 +242,4 @@
 }
 
-# if (0)
-// XXX this needs to be optional (eg, invalid for a mask)
-
-// Apply the BSCALE and BZERO for an image with a "fuzz", so that we get the image as it should be written to
-// disk.
-// The idea is that the "fuzz" (adding a random number between 0 and 1) preserves the expectation value of
-// the image (e.g., a value of 0.1 will get translated to zero 90% of the time, and unity 10% of the time),
-// though at the cost of adding an additional variance of 1/12 (a standard deviation of ~0.29).
-static psImage *scaleImageForDisk(psImage *image, // Image to which to apply BSCALE and BZERO
-                                  int bitpix, // Output BITPIX
-                                  double bscale, // Scaling
-                                  double bzero, // Zero point
-                                  psRandom *rng // Random number generator (for the "fuzz"), or NULL
-                                  )
-{
-    assert(image);
-
-    if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) {
-        return psMemIncrRefCounter(image);
-    }
-
-    psElemType outType;                 // Type for output image
-    // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS
-    // file
-    switch (bitpix) {
-      case 8:
-        outType = PS_TYPE_S8;
-        break;
-      case 16:
-        outType = PS_TYPE_S16;
-        break;
-      case 32:
-        outType = PS_TYPE_S32;
-        break;
-      case 64:
-        outType = PS_TYPE_S64;
-        break;
-      default:
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
-        return NULL;
-    }
-
-    if (bscale == 1.0 && bzero == 0.0) {
-        return psImageCopy(NULL, image, outType);
-    }
-
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
-    psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
-
-    if (!psMemIncrRefCounter(rng)) {
-        // Don't blab about which seed we're going to get --- it's not necessary for this purpose
-        psU64 seed = p_psRandomGetSystemSeed(false);
-        rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
-    }
-
-
-#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
-    case PS_TYPE_##OUTTYPE: { \
-        ps##INTYPE scale = 1.0 / bscale; \
-        ps##INTYPE zero = bzero; \
-        for (int y = 0; y < numRows; y++) { \
-            for (int x = 0; x < numCols; x++) { \
-                /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \
-                ps##INTYPE random = psRandomUniform(rng); \
-                (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
-            } \
-        } \
-        break; \
-    }
-
-#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
-    case PS_TYPE_##INTYPE: { \
-        switch (outType) { \
-            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \
-            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \
-            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \
-            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \
-          default: \
-            psAbort("Should be unreachable."); \
-        } \
-        break; \
-    }
-
-    switch (image->type.type) {
-        SCALE_WRITE_IN_CASE(image, F32, out);
-        SCALE_WRITE_IN_CASE(image, F64, out);
-      default:
-        psAbort("Should be unreachable.");
-    }
-
-    psFree(rng);
-
-    return out;
-}
-# endif
-
-# if (0)
-// XXX supporting code needs to make this an optional operation
-// Determine BSCALE and BZERO for an image, and generate a new image with it applied
-// TRUE = BZERO + BSCALE * FITS
-static psImage *scaleImageDetermine(double *bscale, // Scaling, to return
-                                    double *bzero, // Zero point, to return
-                                    psImage *image, // Image to scale
-                                    int bitpix, // Desired bits per pixel
-                                    psRandom *rng // Random number generator for scaleImageForDisk
-                                    )
-{
-    PS_ASSERT_PTR_NON_NULL(bscale, NULL);
-    PS_ASSERT_PTR_NON_NULL(bzero, NULL);
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
-
-    *bscale = 0.0;
-    *bzero = 0.0;
-
-    switch (bitpix) {
-      case 0:
-        // No scaling applied
-        return psMemIncrRefCounter(image);
-      case 8:
-      case 16:
-      case 32:
-      case 64:
-        // Nothing to do; just allowing these values to pass through
-        break;
-      default:
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
-        return NULL;
-    }
-
-    int numCols = image->numCols, numRows = image->numRows;
-    double range = pow(2.0, bitpix); // Range of values for target BITPIX
-
-#define SCALE_DETERMINE_CASE(IN, INTYPE) \
-    case PS_TYPE_##INTYPE: { \
-        ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \
-        for (int y = 0; y < numRows; y++) { \
-            for (int x = 0; x < numCols; x++) { \
-                ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \
-                if (isfinite(value)) { \
-                    if (value < min) { \
-                        min = value; \
-                    } \
-                    if (value > max) { \
-                        max = value; \
-                    } \
-                } \
-            } \
-        } \
-        if (!isfinite(min) || !isfinite(max)) { \
-            psWarning("No valid values in image to derive BSCALE,BZERO --- using original image."); \
-            *bscale = 1.0; \
-            *bzero = 0.0; \
-            return psMemIncrRefCounter(image); \
-        } \
-        if (min == max) { \
-            *bscale = 1.0; \
-            *bzero = min; \
-        } else { \
-            *bscale = (max - min) / (range - 1.0); \
-            *bzero = min + 0.5 * range * (*bscale); \
-        } \
-        break; \
-    }
-
-    switch (image->type.type) {
-        SCALE_DETERMINE_CASE(image, F32);
-        SCALE_DETERMINE_CASE(image, F64);
-      default:
-        psAbort("Should be unreachable.");
-    }
-    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero);
-
-    return scaleImageForDisk(image, bitpix, *bscale, *bzero, rng);
-}
-# endif
-
-#if 0
-// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
-// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
-// later desire it.
-static psImage *scaleImageFromDisk(psFits *fits, psImage *image)
-{
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-
-    if (bscale == 0.0) {
-        // BSCALE = 0 means don't apply anything
-        return psMemIncrRefCounter(image);
-    }
-
-    psElemType inType = image->type.type; // Type for input image
-    psElemType outType;                 // Type for output image
-    switch (inType) {
-      case PS_TYPE_S8:
-      case PS_TYPE_S16:
-      case PS_TYPE_S32:
-      case PS_TYPE_U8:
-      case PS_TYPE_U16:
-        outType = PS_TYPE_F32;
-        break;
-      case PS_TYPE_S64:
-      case PS_TYPE_U32:
-      case PS_TYPE_U64:
-        outType = PS_TYPE_F64;
-        break;
-        // Including floating-point types just in case someone wants to apply a BSCALE and BZERO to them.
-      case PS_TYPE_F32:
-        outType = PS_TYPE_F32;
-        break;
-      case PS_TYPE_F64:
-        outType = PS_TYPE_F64;
-        break;
-      default:
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unsupported image type: %x", inType);
-        return NULL;
-    }
-
-    int numCols = image->numCols, numRows = image->numRows;
-    psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
-
-
-#define SCALE_READ_OUT_CASE(INTYPE, OUTTYPE) \
-  case PS_TYPE_##OUTTYPE: { \
-      for (int y = 0; y < numRows; y++) { \
-          for (int x = 0; x < numCols; x++) { \
-              out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
-          } \
-      } \
-      break; \
-  }
-
-#define SCALE_READ_IN_CASE(INTYPE) \
-  case PS_TYPE_##INTYPE: { \
-      switch (outType) { \
-          SCALE_READ_OUT_CASE(INTYPE, F32); \
-          SCALE_READ_OUT_CASE(INTYPE, F64); \
-        default: \
-          psAbort("Should never get here: type %x should be F32 or F64", outType); \
-      } \
-      break; \
-  }
-
-    switch (inType) {
-        SCALE_READ_IN_CASE(S8);
-        SCALE_READ_IN_CASE(S16);
-        SCALE_READ_IN_CASE(S32);
-        SCALE_READ_IN_CASE(S64);
-        SCALE_READ_IN_CASE(U8);
-        SCALE_READ_IN_CASE(U16);
-        SCALE_READ_IN_CASE(U32);
-        SCALE_READ_IN_CASE(U64);
-        SCALE_READ_IN_CASE(F32);
-        SCALE_READ_IN_CASE(F64);
-      default:
-          psAbort("Should never get here: type %x should be integer", inType);
-    }
-
-    return out;
-}
-#endif
 
 // Convert an image to the desired BITPIX, i.e., the desired disk representation
@@ -522,43 +263,51 @@
     *floatType = PS_FITS_FLOAT_NONE;
 
+    psFitsOptions *options = fits->options;
+    if (!options) {
+        return psMemIncrRefCounter((psImage*)image); // Casting away const
+    }
+
     // Custom floating-point
-    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix &&
-        fits->floatType != PS_FITS_FLOAT_NONE) {
-        *floatType = fits->floatType;
-        return psFitsFloatImageToDisk(image, fits->floatType);
+    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->conventions.psBitpix &&
+        options->floatType != PS_FITS_FLOAT_NONE) {
+        *floatType = options->floatType;
+        return psFitsFloatImageToDisk(image, options->floatType);
     }
 
     // Automatically select what we're given
-    if (fits->bitpix == 0) {
+    if (options->bitpix == 0) {
         return psMemIncrRefCounter((psImage*)image); // Casting away const
     }
 
     // Quantise floating-point images
-    // XXX this needs to be more controlled: certainly not valid for output masks!
-    # if (0)
-    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) {
+    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->bitpix > 0) {
         if (newScaleZero) {
-            return scaleImageDetermine(bscale, bzero, (psImage*)image, fits->bitpix, rng);
-        }
-        // Get the current BSCALE and BZERO
-        int status = 0;                 // Status of cfitsio
-        if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
-            psFitsError(status, true, "Unable to read header.");
-            return NULL;
-        }
-        status = 0;
-        if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
-            psFitsError(status, true, "Unable to read header.");
-            return NULL;
-        }
-        status = 0;
-        if (*bscale == 0.0) {
-            psError(PS_ERR_IO, true,
-                    "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
-            return NULL;
-        }
-        return scaleImageForDisk((psImage*)image, fits->bitpix, *bscale, *bzero, rng);
-    }
-    # endif
+            // Choose an appropriate BSCALE and BZERO
+            if (!psFitsScaleDetermine(bscale, bzero, image, fits)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine BSCALE and BZERO for image.");
+                return NULL;
+            }
+        } else {
+            // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
+            int status = 0;                 // Status of cfitsio
+            if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
+                psFitsError(status, true, "Unable to read header.");
+                return NULL;
+            }
+            status = 0;
+            if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
+                psFitsError(status, true, "Unable to read header.");
+                return NULL;
+            }
+            status = 0;
+            if (*bscale == 0.0) {
+                psError(PS_ERR_IO, true,
+                        "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
+                return NULL;
+            }
+        }
+
+        return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng);
+    }
 
     // Choose the appropriate output type, given the input type and desired bits per pixel
@@ -576,5 +325,5 @@
     psElemType inType = image->type.type; // Type for input image
     psElemType outType;                 // Type for output image
-    switch (fits->bitpix) {
+    switch (options->bitpix) {
         CONVERT_TYPE_INT_CASE(outType, inType, 8);
         CONVERT_TYPE_INT_CASE(outType, inType, 16);
@@ -585,5 +334,5 @@
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64",
-                fits->bitpix);
+                options->bitpix);
         return NULL;
     }
@@ -592,4 +341,9 @@
         return psMemIncrRefCounter((psImage*)image);
     }
+
+    if (PSELEMTYPE_SIZEOF(inType) > PSELEMTYPE_SIZEOF(outType)) {
+        psWarning("Truncating image pixels to write to disk.");
+    }
+
     return psImageCopy(NULL, image, outType);
 }
@@ -765,5 +519,6 @@
         bzero = cfitsioBzero;
     }
-    assert(bitPix == fits->bitpix || fits->bitpix == 0);
+    psFitsOptions *options = fits->options; // FITS I/O options
+    assert(!options || bitPix == options->bitpix || options->bitpix == 0);
 
     int naxis = 3;                      // Number of axes
@@ -807,5 +562,5 @@
     // an unsigned integer type).  In all other cases, we have already converted the image to use the
     // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
-    fits_set_bscale(fits->fd, 1.0, bzero, &status);
+    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
 
     if (bscale != 0.0) {
@@ -895,5 +650,6 @@
         bzero = cfitsioBzero;
     }
-    assert(bitPix == fits->bitpix || fits->bitpix == 0);
+    psFitsOptions *options = fits->options; // FITS I/O options
+    assert(!options || bitPix == options->bitpix || options->bitpix == 0);
 
     //check to see if the HDU has the same datatype