Index: /trunk/psLib/src/fits/Makefile.am
===================================================================
--- /trunk/psLib/src/fits/Makefile.am	(revision 15334)
+++ /trunk/psLib/src/fits/Makefile.am	(revision 15335)
@@ -9,4 +9,5 @@
 	psFitsImage.c \
 	psFitsTable.c
+#	psFitsFloat.c
 
 EXTRA_DIST = fits.i
@@ -17,4 +18,5 @@
 	psFitsImage.h \
 	psFitsTable.h
+#	psFitsFloat.h
 
 CLEANFILES = *~ *.bb *.bbg *.da
Index: /trunk/psLib/src/fits/psFits.c
===================================================================
--- /trunk/psLib/src/fits/psFits.c	(revision 15334)
+++ /trunk/psLib/src/fits/psFits.c	(revision 15335)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-03 21:27:21 $
+ *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-19 23:52:39 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -163,4 +163,7 @@
     fits->extword = NULL;
     fits->conventions.compression = true;
+    fits->conventions.psBitpix = true;
+    fits->bscale = 0.0;
+    fits->bzero = 0.0;
     psMemSetDeallocator(fits,(psFreeFunc)fitsFree);
 
@@ -397,10 +400,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
-
-    if (! fits->writable) {
-        psError(PS_ERR_IO, true,
-                _("The psFits object is not writable."));
-        return false;
-    }
+    PS_ASSERT_FITS_WRITABLE(fits, false);
 
     // move to the specified HDU
@@ -431,11 +429,6 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_STRING_NON_EMPTY(extname, false);
-
-    if (! fits->writable) {
-        psError(PS_ERR_IO, true,
-                _("The psFits object is not writable."));
-        return false;
-    }
 
     // move to the specified HDU
@@ -509,11 +502,6 @@
 bool psFitsTruncate(psFits* fits)
 {
-    PS_ASSERT_FITS_NON_NULL(fits, NULL);
-
-    if (! fits->writable) {
-        psError(PS_ERR_IO, true,
-                _("The psFits object is not writable."));
-        return PS_FITS_TYPE_NONE;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
 
     int newEnd = psFitsGetExtNum(fits);
@@ -547,4 +535,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
 
     // convert psFitsCompressionType to cfitsio compression types
Index: /trunk/psLib/src/fits/psFits.h
===================================================================
--- /trunk/psLib/src/fits/psFits.h	(revision 15334)
+++ /trunk/psLib/src/fits/psFits.h	(revision 15335)
@@ -4,6 +4,6 @@
  * @author Robert DeSonia, MHPCC
  *
- * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-10-09 03:00:05 $
+ * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-10-19 23:52:39 $
  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
@@ -57,5 +57,8 @@
     struct {
         bool compression;               ///< Compression convention: handling of compressed images
+        bool psBitpix;                  ///< Pan-STARRS BITPIX for floating-point conversion
     } conventions;                      ///< Conventions to honour
+    double bscale;                      ///< Scaling to apply when reading FITS data; 0 for auto
+    double bzero;                       ///< Zero point to apply when reading FITS data; auto if bscale = 0
 } psFits;
 
Index: /trunk/psLib/src/fits/psFitsFloat.c
===================================================================
--- /trunk/psLib/src/fits/psFitsFloat.c	(revision 15335)
+++ /trunk/psLib/src/fits/psFitsFloat.c	(revision 15335)
@@ -0,0 +1,217 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <ieee754.h>
+
+#include "psAbort.h"
+#include "psType.h"
+#include "psError.h"
+#include "psImage.h"
+#include "psMetadata.h"
+#include "psRandom.h"
+#include "psFits.h"
+
+#include "psFitsFloat.h"
+
+#define BIAS_FLOAT_16_0 10              // Exponent bias
+
+union float_16_0 {
+    psS16 s16;                      // 16-bit integer version
+
+    // Floating-point version
+    struct {
+        unsigned int negative:1;    // Sign bit
+        unsigned int exponent:5;    // Exponent bits
+        unsigned int mantissa:10;   // Mantissa bits
+    } f16;
+};
+
+
+
+static inline psS16 convertF32toFloat16_0(psF32 value)
+{
+    union ieee754_float in;             // Input value
+    union float_16_0 out;               // Output value
+
+    // XXX What happens to NAN and INF?
+    in.f = value;
+    out.f16.negative = in.ieee.negative;
+    out.f16.exponent = in.ieee.exponent - IEEE754_FLOAT_BIAS + BIAS_FLOAT_16_0;
+    out.f16.mantissa = in.ieee.mantissa >> 13;
+
+    return out.s16;
+}
+
+static inline psF32 convertF32fromFloat16_0(psS16 value)
+{
+    union float_16_0 in;                // Input value
+    union ieee754_float out;            // Output value
+
+    in.s16 = value;
+    out.ieee.negative = in.f16.negative;
+    out.ieee.exponent = in.f16.exponent + IEEE754_FLOAT_BIAS - BIAS_FLOAT_16_0;
+    out.ieee.mantissa = in.f16.mantissa << 13;
+
+    return out.f;
+}
+
+
+psImage *psFitsScaleGenerate(psFits *fits, const psImage *image, int bitpix, psRandom *rng)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+    PS_ASSERT_FITS_WRITABLE(fits, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
+    PS_ASSERT_RANDOM_NON_NULL(rng, NULL);
+
+    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;
+    }
+
+    int numCols = image->numCols, numRows = image->numRows;
+    double range = pow(2.0, bitpix);    // Range of values for target BITPIX
+    psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
+
+#define SCALE_GENERATE_CASE(IN, INTYPE, OUT, OUTTYPE) \
+    case PS_TYPE_##OUTTYPE: { \
+        for (int y = 0; y < numRows; y++) { \
+            for (int x = 0; x < numCols; x++) { \
+                ps##INTYPE random = psRandomUniform(rng) - 0.5; /* Randomise the quantisation */ \
+                (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
+            } \
+        } \
+        break; \
+    }
+
+#define SCALE_DETERMINE_CASE(IN, INTYPE, OUT) \
+    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; \
+                    } \
+                } \
+            } \
+        } \
+        ps##INTYPE scale = range / (max - min); \
+        ps##INTYPE zero = min - 0.5 * range; /* Minimum value for twos-complement int is - 2^(bitpix-1) */ \
+        fits->bscale = 1.0 / scale; \
+        fits->bzero = zero; \
+        switch (outType) { \
+            SCALE_GENERATE_CASE(IN, INTYPE, OUT, S8); \
+            SCALE_GENERATE_CASE(IN, INTYPE, OUT, S16); \
+            SCALE_GENERATE_CASE(IN, INTYPE, OUT, S32); \
+            SCALE_GENERATE_CASE(IN, INTYPE, OUT, S64); \
+          default: \
+            psAbort("Should be unreachable."); \
+        } \
+        break; \
+    }
+
+    switch (image->type.type) {
+        SCALE_DETERMINE_CASE(image, F32, out);
+        SCALE_DETERMINE_CASE(image, F64, out);
+      default:
+        psAbort("Should be unreachable.");
+    }
+
+    return out;
+}
+
+psImage *psFitsScaleApply(const psImage *image, double bscale, double bzero, psElemType type)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+
+    if (bscale == 0.0) {
+        return psMemIncrRefCounter(image);
+    }
+
+    switch (
+
+
+psImage *psFitsFloatImageWrite(const psImage *image,
+                               psMetadata *header,
+                               psFitsFloat type
+    )
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_METADATA_NON_NULL(header, NULL); // Need the header, so we can mark it as compressed
+
+    psImage *output = NULL;             // Output image, to return
+    psString convName = NULL;           // Convention name
+
+    switch (type) {
+      case PS_FITS_FLOAT_16_0:
+        output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S16); // Output image
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                output->data.S16[y][x] = convertF32toFloat16_0(image->data.F32[y][x]);
+            }
+        }
+        convName = psStringCopy("FLOAT_16_0");
+        break;
+      default:
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised custom floating-point type: %x", type);
+        return NULL;
+    }
+
+    psMetadataAddStr(header, PS_LIST_TAIL, "PSBITPIX", PS_META_REPLACE, "Pan-STARRS bit encoding", convName);
+
+    return output;
+}
+
+psImage *psFitsFloatImageRead(const psFits *fits,
+                              const psImage *image
+    )
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+
+    bool mdok;                          // Status of MD lookup
+    const char *convName = psMetadataLookupStr(&mdok, header, "PSBITPIX"); // Convention used
+    if (!mdok || !convName) {
+        // It's not a custom floating-point image
+        return psMemIncrRefCounter(image);
+    }
+
+    psImage *output = NULL;             // Output image, to return
+
+    if (strcmp(convName, "FLOAT_16_0") == 0 && image->type.type == PS_TYPE_S16) {
+        output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); // Output image
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                output->data.F32[y][x] = convertF32fromFloat16_0(image->data.S16[y][x]);
+            }
+        }
+        return output;
+    }
+
+    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise PSBITPIX convention name: %s", convName);
+    return NULL;
+
+    return output;
+}
Index: /trunk/psLib/src/fits/psFitsFloat.h
===================================================================
--- /trunk/psLib/src/fits/psFitsFloat.h	(revision 15335)
+++ /trunk/psLib/src/fits/psFitsFloat.h	(revision 15335)
@@ -0,0 +1,20 @@
+#ifndef PS_FITS_FLOAT_H
+#define PS_FITS_FLOAT_H
+
+
+#include <psImage.h>
+#include <psMetadata.h>
+
+/// Type of custom floating point
+typedef enum {
+    PS_FITS_FLOAT_16_0,                 ///< Original F16 proposal: 1S 5E 10M
+} psFitsFloat;
+
+/// Convert an image to custom floating-point in preparation for writing as FITS
+psImage *psFitsConvertFloatImage(const psImage *image, ///< Image to convert
+                                 psMetadata *header, ///< Header to update
+                                 psFitsFloat type ///< Custom floating point type
+    );
+
+
+#endif
Index: /trunk/psLib/src/fits/psFitsHeader.c
===================================================================
--- /trunk/psLib/src/fits/psFitsHeader.c	(revision 15334)
+++ /trunk/psLib/src/fits/psFitsHeader.c	(revision 15335)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-04 01:16:33 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-19 23:52:39 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -593,4 +593,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_METADATA_NON_NULL(output, false);
 
@@ -604,4 +605,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
 
     // We allow output == NULL in order to write a minimal header.
Index: /trunk/psLib/src/fits/psFitsImage.c
===================================================================
--- /trunk/psLib/src/fits/psFitsImage.c	(revision 15334)
+++ /trunk/psLib/src/fits/psFitsImage.c	(revision 15335)
@@ -7,6 +7,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-03 21:27:21 $
+ *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-19 23:52:39 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -185,4 +185,33 @@
 }
 
+
+bool psFitsImageSize(int *numCols, int *numRows, psElemType *type, const psFits *fits, psRegion region)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+
+    if (psFitsCheckSingleCompressedImagePHU(fits, NULL)) {
+        // This is really what we want, not the empty PHU
+        psTrace("psLib.fits", 1,
+                "This PHU should really be a compressed image --- reading that image instead.");
+    }
+
+    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, 0); // How big the region to read is
+
+    if (numCols) {
+        *numCols = info->lastPixel[0] - info->firstPixel[0] + 1;
+    }
+    if (numRows) {
+        *numRows = info->lastPixel[1] - info->firstPixel[1] + 1;
+    }
+    if (type) {
+        *type = info->psDatatype;
+    }
+
+    psFree(info);
+
+    return true;
+}
+
+
 // Read into an extant image of just the right size
 static bool fitsReadImage(psImage *output,   // Output image
@@ -287,4 +316,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_IMAGE_NON_NULL(input, false);
     // this is equivalent to insert after the last HDU
@@ -302,4 +332,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_IMAGE_NON_NULL(input, false);
 
@@ -402,4 +433,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_IMAGE_NON_NULL(input, false);
 
@@ -556,4 +588,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_ARRAY_NON_NULL(input, false);
 
@@ -622,4 +655,5 @@
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
     PS_ASSERT_ARRAY_NON_NULL(input, false);
 
Index: /trunk/psLib/src/fits/psFitsImage.h
===================================================================
--- /trunk/psLib/src/fits/psFitsImage.h	(revision 15334)
+++ /trunk/psLib/src/fits/psFitsImage.h	(revision 15335)
@@ -4,6 +4,6 @@
  * @author Robert DeSonia, MHPCC
  *
- * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-23 22:47:23 $
+ * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-10-19 23:52:39 $
  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
@@ -21,4 +21,11 @@
 #include "psMetadata.h"
 #include "psImage.h"
+
+/// Return the dimensions and type of the FITS image
+bool psFitsImageSize(int *numCols, int *numRows, ///< Size of image
+                     psElemType *type,  ///< Type of image
+                     const psFits *fits, ///< FITS file pointer
+                     psRegion region    ///< Region in the FITS image to read
+    );
 
 /** Reads an image, given the desired region and z-plane.