IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15630


Ignore:
Timestamp:
Nov 15, 2007, 3:04:56 PM (18 years ago)
Author:
Paul Price
Message:

Integrating the code to choose the output BITPIX (and quantising floating-point images if necessary) and custom floating-point schemes with the functions to read and write images. More detailed description, from my e-mail: psLib has always read any of the (legal) FITS types, and written a FITS image with the appropriate BITPIX for the image->type.type. Now the output BITPIX can be forced by setting psFits.bitpix to the desired BITPIX. Floating point images being quantised to integer types have BSCALE and BZERO set in such a fashion as to maintain the dynamic range of the image; quantisation after applying the scale and zero includes adding a random number [0,1) to preserve the expectation of the floating-point value. On reading, the BSCALE and BZERO are applied by cfitsio without any assistance by psLib. I've also added Gene's custom 16-bit floats (1S 5E 10M). In anticipation of further flavours down the road and to maintain compatibility, I've labelled this as FLOAT_16_0. An F32 image may be written out with this scheme by setting psFits.floatType = PS_FITS_FLOAT_16_0. The image is converted and written out with BITPIX=16, and the header keyword PSBITPIX added with the custom floating-point scheme name. Such images are automatically converted to F32 on reading (if psFits.conventions.psBitpix is true, which it is by default). One thing to keep in mind is that I'm not sure how NANs and INFs are handled in these processes. These additions are made on the top of the compression feature --- quantised or custom floating-point images may also be compressed (though the custom floating-point images don't gain much from it).

Location:
trunk/psLib/src
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fits/Makefile.am

    r15335 r15630  
    88        psFitsHeader.c \
    99        psFitsImage.c \
    10         psFitsTable.c
    11 #       psFitsFloat.c
     10        psFitsTable.c \
     11        psFitsFloat.c \
     12        psFitsFloatFile.c
    1213
    1314EXTRA_DIST = fits.i
    … …  
    1718        psFitsHeader.h \
    1819        psFitsImage.h \
    19         psFitsTable.h
    20 #       psFitsFloat.h
     20        psFitsTable.h \
     21        psFitsFloat.h \
     22        psFitsFloatFile.h
    2123
    2224CLEANFILES = *~ *.bb *.bbg *.da

  • trunk/psLib/src/fits/psFits.h

    r15335 r15630  
    44 * @author Robert DeSonia, MHPCC
    55 *
    6  * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-10-19 23:52:39 $
     6 * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2007-11-16 01:04:56 $
    88 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    99 */
    … …  
    2121#include "psMetadata.h"
    2222#include "psImage.h"
     23#include "psFitsFloat.h"
    2324
    2425#include "psErrorCodes.h"
    … …  
    5758    struct {
    5859        bool compression;               ///< Compression convention: handling of compressed images
    59         bool psBitpix;                  ///< Pan-STARRS BITPIX for floating-point conversion
     60        bool psBitpix;                  ///< Custom floating-point image
    6061    } conventions;                      ///< Conventions to honour
    61     double bscale;                      ///< Scaling to apply when reading FITS data; 0 for auto
    62     double bzero;                       ///< Zero point to apply when reading FITS data; auto if bscale = 0
     62    int bitpix;                         ///< Desired BITPIX for output images; 0 to use as provided
     63    psFitsFloat floatType;              ///< Desired custom floating-point for output images
    6364} psFits;
    6465

  • trunk/psLib/src/fits/psFitsFloat.c

    r15335 r15630  
    44
    55#include <ieee754.h>
     6#include <string.h>
     7#include <assert.h>
    68
    79#include "psAbort.h"
    … …  
    911#include "psError.h"
    1012#include "psImage.h"
    11 #include "psMetadata.h"
    12 #include "psRandom.h"
    1313#include "psFits.h"
     14#include "psTrace.h"
     15#include "psMemory.h"
    1416
    1517#include "psFitsFloat.h"
    1618
    17 #define BIAS_FLOAT_16_0 10              // Exponent bias
     19#define BIAS_FLOAT_16_0              10 // Exponent bias for FLOAT_16_0
     20
     21//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     22// Private functions
     23//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1824
    1925union float_16_0 {
    … …  
    2733    } f16;
    2834};
    29 
    3035
    3136
    … …  
    5762}
    5863
     64//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     65// Public functions
     66//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    5967
    60 psImage *psFitsScaleGenerate(psFits *fits, const psImage *image, int bitpix, psRandom *rng)
    61 {
    62     PS_ASSERT_FITS_NON_NULL(fits, NULL);
    63     PS_ASSERT_FITS_WRITABLE(fits, NULL);
    64     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    65     PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
    66     PS_ASSERT_RANDOM_NON_NULL(rng, NULL);
    6768
    68     psElemType outType;                 // Type for output image
    69     // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS
    70     // file.
    71     switch (bitpix) {
    72       case 8:
    73         outType = PS_TYPE_S8;
    74         break;
    75       case 16:
    76         outType = PS_TYPE_S16;
    77         break;
    78       case 32:
    79         outType = PS_TYPE_S32;
    80         break;
    81       case 64:
    82         outType = PS_TYPE_S64;
    83         break;
    84       default:
    85         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
    86         return NULL;
    87     }
    88 
    89     int numCols = image->numCols, numRows = image->numRows;
    90     double range = pow(2.0, bitpix);    // Range of values for target BITPIX
    91     psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
    92 
    93 #define SCALE_GENERATE_CASE(IN, INTYPE, OUT, OUTTYPE) \
    94     case PS_TYPE_##OUTTYPE: { \
    95         for (int y = 0; y < numRows; y++) { \
    96             for (int x = 0; x < numCols; x++) { \
    97                 ps##INTYPE random = psRandomUniform(rng) - 0.5; /* Randomise the quantisation */ \
    98                 (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
    99             } \
    100         } \
    101         break; \
    102     }
    103 
    104 #define SCALE_DETERMINE_CASE(IN, INTYPE, OUT) \
    105     case PS_TYPE_##INTYPE: { \
    106         ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \
    107         for (int y = 0; y < numRows; y++) { \
    108             for (int x = 0; x < numCols; x++) { \
    109                 ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \
    110                 if (isfinite(value)) { \
    111                     if (value < min) { \
    112                         min = value; \
    113                     } \
    114                     if (value > max) { \
    115                         max = value; \
    116                     } \
    117                 } \
    118             } \
    119         } \
    120         ps##INTYPE scale = range / (max - min); \
    121         ps##INTYPE zero = min - 0.5 * range; /* Minimum value for twos-complement int is - 2^(bitpix-1) */ \
    122         fits->bscale = 1.0 / scale; \
    123         fits->bzero = zero; \
    124         switch (outType) { \
    125             SCALE_GENERATE_CASE(IN, INTYPE, OUT, S8); \
    126             SCALE_GENERATE_CASE(IN, INTYPE, OUT, S16); \
    127             SCALE_GENERATE_CASE(IN, INTYPE, OUT, S32); \
    128             SCALE_GENERATE_CASE(IN, INTYPE, OUT, S64); \
    129           default: \
    130             psAbort("Should be unreachable."); \
    131         } \
    132         break; \
    133     }
    134 
    135     switch (image->type.type) {
    136         SCALE_DETERMINE_CASE(image, F32, out);
    137         SCALE_DETERMINE_CASE(image, F64, out);
    138       default:
    139         psAbort("Should be unreachable.");
    140     }
    141 
    142     return out;
    143 }
    144 
    145 psImage *psFitsScaleApply(const psImage *image, double bscale, double bzero, psElemType type)
     69psImage *psFitsFloatImageWrite(const psImage *image, psFitsFloat type)
    14670{
    14771    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    148 
    149     if (bscale == 0.0) {
    150         return psMemIncrRefCounter(image);
    151     }
    152 
    153     switch (
    154 
    155 
    156 psImage *psFitsFloatImageWrite(const psImage *image,
    157                                psMetadata *header,
    158                                psFitsFloat type
    159     )
    160 {
    161     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    162     PS_ASSERT_METADATA_NON_NULL(header, NULL); // Need the header, so we can mark it as compressed
     72    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
    16373
    16474    psImage *output = NULL;             // Output image, to return
    165     psString convName = NULL;           // Convention name
    16675
    16776    switch (type) {
     77      case PS_FITS_FLOAT_NONE:
     78        // No conversion to be performed
     79        return psMemIncrRefCounter((psImage*)image); // Casting away "const"
    16880      case PS_FITS_FLOAT_16_0:
    16981        output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S16); // Output image
    … …  
    17385            }
    17486        }
    175         convName = psStringCopy("FLOAT_16_0");
    17687        break;
    17788      default:
    … …  
    18091    }
    18192
    182     psMetadataAddStr(header, PS_LIST_TAIL, "PSBITPIX", PS_META_REPLACE, "Pan-STARRS bit encoding", convName);
    183 
    18493    return output;
    18594}
    18695
    187 psImage *psFitsFloatImageRead(const psFits *fits,
    188                               const psImage *image
    189     )
     96
     97psImage *psFitsFloatImageRead(psImage *out, const psImage *in, psFitsFloat type)
    19098{
    191     PS_ASSERT_FITS_NON_NULL(fits, NULL);
    192     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     99    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
    193100
    194     bool mdok;                          // Status of MD lookup
    195     const char *convName = psMetadataLookupStr(&mdok, header, "PSBITPIX"); // Convention used
    196     if (!mdok || !convName) {
    197         // It's not a custom floating-point image
    198         return psMemIncrRefCounter(image);
     101    psElemType elem = psFitsFloatImageType(type); // Type for elements
     102    if (elem == 0) {
     103        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise convention: %x", type);
     104        return NULL;
    199105    }
    200106
    201     psImage *output = NULL;             // Output image, to return
     107    int numCols = in->numCols, numRows = in->numRows; // Size of image
    202108
    203     if (strcmp(convName, "FLOAT_16_0") == 0 && image->type.type == PS_TYPE_S16) {
    204         output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); // Output image
    205         for (int y = 0; y < image->numRows; y++) {
    206             for (int x = 0; x < image->numCols; x++) {
    207                 output->data.F32[y][x] = convertF32fromFloat16_0(image->data.S16[y][x]);
     109    out = psImageRecycle(out, numCols, numRows, elem);
     110
     111    switch (type) {
     112      case PS_FITS_FLOAT_16_0:
     113        if (in->type.type != PS_TYPE_S16) {
     114            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     115                    "Convention claims to be FLOAT_16_0, but image type is not S16.");
     116            return NULL;
     117        }
     118        assert(out->type.type == PS_TYPE_F32);
     119        for (int y = 0; y < numRows; y++) {
     120            for (int x = 0; x < numCols; x++) {
     121                out->data.F32[y][x] = convertF32fromFloat16_0(in->data.S16[y][x]);
    208122            }
    209123        }
    210         return output;
     124        return out;
     125      case PS_FITS_FLOAT_NONE:
     126      default:
     127        psAbort("Should be unreachable");
     128    }
     129    return NULL;
     130}
     131
     132psElemType psFitsFloatImageType(psFitsFloat type)
     133{
     134    switch (type) {
     135      case PS_FITS_FLOAT_16_0:
     136        return PS_TYPE_F32;
     137      case PS_FITS_FLOAT_NONE:
     138      default:
     139        return 0;                       // Doesn't correspond to ANY type --- should flag a real error
     140    }
     141}
     142
     143
     144psFitsFloat psFitsFloatTypeFromString(const char *string)
     145{
     146    PS_ASSERT_STRING_NON_EMPTY(string, PS_FITS_FLOAT_NONE);
     147
     148    if (strcmp(string, "FLOAT_16_0") == 0) {
     149        return PS_FITS_FLOAT_16_0;
    211150    }
    212151
    213     psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise PSBITPIX convention name: %s", convName);
    214     return NULL;
    215 
    216     return output;
     152    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     153            "Unable to recognise custom floating-point convention name: %s", string);
     154    return PS_FITS_FLOAT_NONE;
    217155}

  • trunk/psLib/src/fits/psFitsFloat.h

    r15335 r15630  
    22#define PS_FITS_FLOAT_H
    33
    4 
     4#include <psFits.h>
     5#include <psType.h>
    56#include <psImage.h>
    6 #include <psMetadata.h>
    77
    88/// Type of custom floating point
    99typedef enum {
     10    PS_FITS_FLOAT_NONE,                 ///< No conversion to be performed
    1011    PS_FITS_FLOAT_16_0,                 ///< Original F16 proposal: 1S 5E 10M
    1112} psFitsFloat;
    1213
    1314/// Convert an image to custom floating-point in preparation for writing as FITS
    14 psImage *psFitsConvertFloatImage(const psImage *image, ///< Image to convert
    15                                  psMetadata *header, ///< Header to update
    16                                  psFitsFloat type ///< Custom floating point type
     15psImage *psFitsFloatImageWrite(const psImage *image, ///< Image to convert
     16                               psFitsFloat type ///< Custom floating point type
    1717    );
    1818
     19/// Convert the custom floating-point image to a floating-point image
     20psImage *psFitsFloatImageRead(psImage *out, ///< Output image, or NULL
     21                              const psImage *in, ///< Image to convert
     22                              psFitsFloat type ///< Custom floating point type
     23    );
     24
     25/// Return the appropriate element type for a custom floating-point
     26psElemType psFitsFloatImageType(psFitsFloat type ///< Custom floating-point type
     27    );
     28
     29/// Return the custom floating-point type from a string description
     30psFitsFloat psFitsFloatTypeFromString(const char *string ///< String with name of type
     31    );
    1932
    2033#endif

  • trunk/psLib/src/fits/psFitsHeader.c

    r15335 r15630  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-10-19 23:52:39 $
     9 *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-11-16 01:04:56 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4444static const char *noWriteFitsKeys[] = { "SIMPLE", "XTENSION", "BITPIX", "NAXIS", "EXTNAME", "BSCALE",
    4545                                         "BZERO", "TFIELDS", "PCOUNT", "GCOUNT", "ZIMAGE", "ZBITPIX",
    46                                          "ZCMPTYPE", NULL};
     46                                         "ZCMPTYPE", "PSBITPIX", NULL};
    4747
    4848// List of the start of FITS header keys not to write (handled by cfitsio); NULL-terminated

  • trunk/psLib/src/fits/psFitsImage.c

    r15335 r15630  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-10-19 23:52:39 $
     9 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-11-16 01:04:56 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2121#include <string.h>
    2222
    23 #include "psFits.h"
     23#include "psAbort.h"
     24#include "psType.h"
    2425#include "psAssert.h"
    25 #include "psFitsImage.h"
    2626#include "psError.h"
    27 
    28 #include "psImageStructManip.h"
    29 #include "psMemory.h"
    3027#include "psString.h"
    3128#include "psLogMsg.h"
    3229#include "psTrace.h"
    3330#include "psVector.h"
    34 
     31#include "psRandom.h"
     32#include "psImageStructManip.h"
     33
     34#include "psFits.h"
     35#include "psFitsFloat.h"
     36#include "psFitsFloatFile.h"
    3537#include "psFitsHeader.h"
    3638
    37 #define MAX_STRING_LENGTH 256  // maximum length string for FITS routines
     39#include "psMemory.h"
     40
     41#include "psFitsImage.h"
     42
     43#define MAX_STRING_LENGTH 256           // maximum length string for FITS routines
    3844
    3945// Information required to read a FITS file
    … …  
    4955} p_psFitsReadInfo;
    5056
    51 static p_psFitsReadInfo *p_psFitsReadInfoAlloc(const psFits *fits, // The FITS file handle
    52         psRegion region, // Region to read
    53         int z // z-plane to read in cube
    54                                               )
     57// Read the vital statistics of this FITS image, in preparation for reading the image
     58static p_psFitsReadInfo *p_psFitsReadInfoAlloc(
     59    const psFits *fits, // The FITS file handle
     60    psRegion region, // Region to read
     61    int z // z-plane to read in cube
     62    )
    5563{
    5664    PS_ASSERT_FITS_NON_NULL(fits, NULL);
    … …  
    6169
    6270    int status = 0;                     // CFITSIO status
    63     char fitsErr[80];                   // CFITSIO error message string
    6471
    6572    // check to see if we even are positioned on an image HDU
    6673    int hduType;                        // Type of HDU
    6774    if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
    68         (void)fits_get_errstatus(status, fitsErr);
    69         psError(PS_ERR_IO, true,
    70                 _("Could not determine the HDU type. CFITSIO Error: %s"),
    71                 fitsErr);
     75        psFitsError(status, true, "Could not determine the HDU type.");
    7276        goto bad;
    7377    }
    7478    if (hduType != IMAGE_HDU) {
    75         psError(PS_ERR_IO, true,
    76                 _("Current FITS HDU type must be an image."));
     79        psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image."));
    7780        goto bad;
    7881    }
    … …  
    8083    // Get the data type 'bitPix' from the FITS image
    8184    if (fits_get_img_equivtype(fits->fd, &info->bitPix, &status) != 0) {
    82         fits_get_errstatus(status, fitsErr);
    83         psError(PS_ERR_IO, true,
    84                 _("Could not determine image data type. CFITSIO Error: %s"),
    85                 fitsErr);
     85        psFitsError(status, true, "Could not determine image data type.");
    8686        goto bad;
    8787    }
    … …  
    8989    /* Get the dimensions 'nAxis' from the FITS image */
    9090    if (fits_get_img_dim(fits->fd, &info->nAxis, &status) != 0) {
    91         (void)fits_get_errstatus(status, fitsErr);
    92         psError(PS_ERR_IO, true,
    93                 _("Could not determine image dimensions. CFITSIO Error: %s"),
    94                 fitsErr);
     91        psFitsError(status, true, "Could not determine image dimensions.");
    9592        goto bad;
    9693    }
    … …  
    9996    if ((info->nAxis < 2) || (info->nAxis > 3)) {
    10097        psError(PS_ERR_IO, true,
    101                 _("Image number of dimensions, %d, is not valid.  Only two or three dimensions supported for FITS I/O."),
    102                 info->nAxis);
     98                _("Image number of dimensions, %d, is not supported."), info->nAxis);
    10399        goto bad;
    104100    }
    … …  
    106102    /* Get the Image size from the FITS file */
    107103    if (fits_get_img_size(fits->fd, info->nAxis, info->nAxes, &status) != 0) {
    108         (void)fits_get_errstatus(status, fitsErr);
    109         psError(PS_ERR_IO, true,
    110                 _("Could not determine image size. CFITSIO Error: %s"),
    111                 fitsErr);
     104        psFitsError(status, true, "Could not determine image size.");
    112105        goto bad;
    113106    }
    … …  
    133126    info->increment[2] = 1;
    134127
    135     switch (info->bitPix) {
    136     case BYTE_IMG:
    137         info->psDatatype = PS_TYPE_U8;
    138         info->fitsDatatype = TBYTE;
    139         break;
    140     case SBYTE_IMG:
    141         info->psDatatype = PS_TYPE_S8;
    142         info->fitsDatatype = TSBYTE;
    143         break;
    144     case USHORT_IMG:
    145         info->psDatatype = PS_TYPE_U16;
    146         info->fitsDatatype = TUSHORT;
    147         break;
    148     case SHORT_IMG:
    149         info->psDatatype = PS_TYPE_S16;
    150         info->fitsDatatype = TSHORT;
    151         break;
    152     case ULONG_IMG:
    153         info->psDatatype = PS_TYPE_U32;
    154         info->fitsDatatype = TUINT;
    155         break;
    156     case LONG_IMG:
    157         info->psDatatype = PS_TYPE_S32;
    158         info->fitsDatatype = TINT;
    159         break;
    160     case LONGLONG_IMG:
    161         info->psDatatype = PS_TYPE_S64;
    162         info->fitsDatatype = TLONGLONG;
    163         break;
    164     case FLOAT_IMG:
    165         info->psDatatype = PS_TYPE_F32;
    166         info->fitsDatatype = TFLOAT;
    167         break;
    168     case DOUBLE_IMG:
    169         info->psDatatype = PS_TYPE_F64;
    170         info->fitsDatatype = TDOUBLE;
    171         break;
    172     default:
    173         psError(PS_ERR_IO, true,
    174                 _("FITS image type, BITPIX=%d, is not supported."),
    175                 info->bitPix);
     128    // Check scale and zero
     129    double bscale = 0.0, bzero = 0.0;    // Scale and zero point
     130    if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) {
     131        psFitsError(status, true, "Unable to read header.");
    176132        goto bad;
    177133    }
    178 
     134    status = 0;
     135    if (fits_read_key_dbl(fits->fd, "BZERO", &bzero, NULL, &status) && status != KEY_NO_EXIST) {
     136        psFitsError(status, true, "Unable to read header.");
     137        goto bad;
     138    }
     139    status = 0;
     140
     141    if ((bscale != 0.0 && bscale != 1.0) || bzero != (int)bzero) {
     142        // It's a floating-point image that's been quantised
     143        // cfitsio will apply the scale and zero point for us if we choose the correct data type
     144        switch (info->bitPix) {
     145          case BYTE_IMG:
     146          case SBYTE_IMG:
     147          case USHORT_IMG:
     148          case SHORT_IMG:
     149          case ULONG_IMG:
     150          case LONG_IMG:
     151          case FLOAT_IMG:
     152            info->psDatatype = PS_TYPE_F32;
     153            info->fitsDatatype = TFLOAT;
     154            break;
     155          case LONGLONG_IMG:
     156          case DOUBLE_IMG:
     157            info->psDatatype = PS_TYPE_F64;
     158            info->fitsDatatype = TDOUBLE;
     159            break;
     160          default:
     161            psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix);
     162            goto bad;
     163        }
     164    } else {
     165        switch (info->bitPix) {
     166          case BYTE_IMG:
     167            info->psDatatype = PS_TYPE_U8;
     168            info->fitsDatatype = TBYTE;
     169            break;
     170          case SBYTE_IMG:
     171            info->psDatatype = PS_TYPE_S8;
     172            info->fitsDatatype = TSBYTE;
     173            break;
     174          case USHORT_IMG:
     175            info->psDatatype = PS_TYPE_U16;
     176            info->fitsDatatype = TUSHORT;
     177            break;
     178          case SHORT_IMG:
     179            info->psDatatype = PS_TYPE_S16;
     180            info->fitsDatatype = TSHORT;
     181            break;
     182          case ULONG_IMG:
     183            info->psDatatype = PS_TYPE_U32;
     184            info->fitsDatatype = TUINT;
     185            break;
     186          case LONG_IMG:
     187            info->psDatatype = PS_TYPE_S32;
     188            info->fitsDatatype = TINT;
     189            break;
     190          case LONGLONG_IMG:
     191            info->psDatatype = PS_TYPE_S64;
     192            info->fitsDatatype = TLONGLONG;
     193            break;
     194          case FLOAT_IMG:
     195            info->psDatatype = PS_TYPE_F32;
     196            info->fitsDatatype = TFLOAT;
     197            break;
     198          case DOUBLE_IMG:
     199            info->psDatatype = PS_TYPE_F64;
     200            info->fitsDatatype = TDOUBLE;
     201            break;
     202          default:
     203            psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix);
     204            goto bad;
     205        }
     206    }
    179207    return info;
    180208
    … …  
    211239
    212240    return true;
     241}
     242
     243// Apply the BSCALE and BZERO for an image with a "fuzz"
     244// The idea is that the "fuzz" (adding a random number between 0 and 1) preserves the expectation value of
     245// the image (e.g., a value of 0.1 will get translated to zero 90% of the time, and unity 10% of the time),
     246// though at the cost of adding an additional variance of 1/12 (a standard deviation of ~0.29).
     247static psImage *scaleImageWrite(psImage *image, // Image to which to apply BSCALE and BZERO
     248                                int bitpix, // Output BITPIX
     249                                double bscale, // Scaling
     250                                double bzero, // Zero point
     251                                psRandom *rng // Random number generator (for the "fuzz"), or NULL
     252    )
     253{
     254    assert(image);
     255
     256    if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) {
     257        return psMemIncrRefCounter(image);
     258    }
     259
     260    psElemType outType;                 // Type for output image
     261    // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS
     262    // file
     263    switch (bitpix) {
     264      case 8:
     265        outType = PS_TYPE_S8;
     266        break;
     267      case 16:
     268        outType = PS_TYPE_S16;
     269        break;
     270      case 32:
     271        outType = PS_TYPE_S32;
     272        break;
     273      case 64:
     274        outType = PS_TYPE_S64;
     275        break;
     276      default:
     277        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
     278        return NULL;
     279    }
     280
     281    if (bscale == 1.0 && bzero == 0.0) {
     282        return psImageCopy(NULL, image, outType);
     283    }
     284
     285    int numCols = image->numCols, numRows = image->numRows; // Size of image
     286    psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
     287
     288    if (!psMemIncrRefCounter(rng)) {
     289        // Don't blab about which seed we're going to get --- it's not necessary for this purpose
     290        psU64 seed = p_psRandomGetSystemSeed(false);
     291        rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
     292    }
     293
     294
     295#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
     296    case PS_TYPE_##OUTTYPE: { \
     297        ps##INTYPE scale = 1.0 / bscale; \
     298        ps##INTYPE zero = bzero; \
     299        for (int y = 0; y < numRows; y++) { \
     300            for (int x = 0; x < numCols; x++) { \
     301                /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \
     302                ps##INTYPE random = psRandomUniform(rng); \
     303                (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
     304            } \
     305        } \
     306        break; \
     307    }
     308
     309#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
     310    case PS_TYPE_##INTYPE: { \
     311        switch (outType) { \
     312            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \
     313            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \
     314            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \
     315            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \
     316          default: \
     317            psAbort("Should be unreachable."); \
     318        } \
     319        break; \
     320    }
     321
     322    switch (image->type.type) {
     323        SCALE_WRITE_IN_CASE(image, F32, out);
     324        SCALE_WRITE_IN_CASE(image, F64, out);
     325      default:
     326        psAbort("Should be unreachable.");
     327    }
     328
     329    psFree(rng);
     330
     331    return out;
     332}
     333
     334
     335
     336// Determine BSCALE and BZERO for an image, and generate a new image with it applied
     337// TRUE = BZERO + BSCALE * FITS
     338static psImage *scaleImageDetermineWrite(double *bscale, // Scaling, to return
     339                                         double *bzero, // Zero point, to return
     340                                         psImage *image, // Image to scale
     341                                         int bitpix, // Desired bits per pixel
     342                                         psRandom *rng // Random number generator for scaleImageWrite
     343    )
     344{
     345    PS_ASSERT_PTR_NON_NULL(bscale, NULL);
     346    PS_ASSERT_PTR_NON_NULL(bzero, NULL);
     347    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     348    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
     349
     350    *bscale = 0.0;
     351    *bzero = 0.0;
     352
     353    switch (bitpix) {
     354      case 0:
     355        // No scaling applied
     356        return psMemIncrRefCounter(image);
     357      case 8:
     358      case 16:
     359      case 32:
     360      case 64:
     361        // Nothing to do; just allowing these values to pass through
     362        break;
     363      default:
     364        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
     365        return NULL;
     366    }
     367
     368    int numCols = image->numCols, numRows = image->numRows;
     369    double range = pow(2.0, bitpix); // Range of values for target BITPIX
     370
     371#define SCALE_DETERMINE_CASE(IN, INTYPE) \
     372    case PS_TYPE_##INTYPE: { \
     373        ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \
     374        for (int y = 0; y < numRows; y++) { \
     375            for (int x = 0; x < numCols; x++) { \
     376                ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \
     377                if (isfinite(value)) { \
     378                    if (value < min) { \
     379                        min = value; \
     380                    } \
     381                    if (value > max) { \
     382                        max = value; \
     383                    } \
     384                } \
     385            } \
     386        } \
     387        if (!isfinite(min) || !isfinite(max)) { \
     388            psWarning("No valid values in image to derive BSCALE,BZERO --- using original image."); \
     389            *bscale = 1.0; \
     390            *bzero = 0.0; \
     391            return psMemIncrRefCounter(image); \
     392        } \
     393        if (min == max) { \
     394            *bscale = 1.0; \
     395            *bzero = min; \
     396        } else { \
     397            *bscale = (max - min) / (range - 1.0); \
     398            *bzero = min + 0.5 * range * (*bscale); \
     399        } \
     400        break; \
     401    }
     402
     403    switch (image->type.type) {
     404        SCALE_DETERMINE_CASE(image, F32);
     405        SCALE_DETERMINE_CASE(image, F64);
     406      default:
     407        psAbort("Should be unreachable.");
     408    }
     409    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero);
     410
     411    return scaleImageWrite(image, bitpix, *bscale, *bzero, rng);
     412}
     413
     414
     415#if 0
     416// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
     417// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
     418// later desire it.
     419static psImage *scaleImageRead(psFits *fits, psImage *image)
     420{
     421    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     422
     423    if (bscale == 0.0) {
     424        // BSCALE = 0 means don't apply anything
     425        return psMemIncrRefCounter(image);
     426    }
     427
     428    psElemType inType = image->type.type; // Type for input image
     429    psElemType outType;                 // Type for output image
     430    switch (inType) {
     431      case PS_TYPE_S8:
     432      case PS_TYPE_S16:
     433      case PS_TYPE_S32:
     434      case PS_TYPE_U8:
     435      case PS_TYPE_U16:
     436        outType = PS_TYPE_F32;
     437        break;
     438      case PS_TYPE_S64:
     439      case PS_TYPE_U32:
     440      case PS_TYPE_U64:
     441        outType = PS_TYPE_F64;
     442        break;
     443        // Including floating-point types just in case someone wants to apply a BSCALE and BZERO to them.
     444      case PS_TYPE_F32:
     445        outType = PS_TYPE_F32;
     446        break;
     447      case PS_TYPE_F64:
     448        outType = PS_TYPE_F64;
     449        break;
     450      default:
     451        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unsupported image type: %x", inType);
     452        return NULL;
     453    }
     454
     455    int numCols = image->numCols, numRows = image->numRows;
     456    psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
     457
     458
     459#define SCALE_READ_OUT_CASE(INTYPE, OUTTYPE) \
     460  case PS_TYPE_##OUTTYPE: { \
     461      for (int y = 0; y < numRows; y++) { \
     462          for (int x = 0; x < numCols; x++) { \
     463              out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
     464          } \
     465      } \
     466      break; \
     467  }
     468
     469#define SCALE_READ_IN_CASE(INTYPE) \
     470  case PS_TYPE_##INTYPE: { \
     471      switch (outType) { \
     472          SCALE_READ_OUT_CASE(INTYPE, F32); \
     473          SCALE_READ_OUT_CASE(INTYPE, F64); \
     474        default: \
     475          psAbort("Should never get here: type %x should be F32 or F64", outType); \
     476      } \
     477      break; \
     478  }
     479
     480    switch (inType) {
     481        SCALE_READ_IN_CASE(S8);
     482        SCALE_READ_IN_CASE(S16);
     483        SCALE_READ_IN_CASE(S32);
     484        SCALE_READ_IN_CASE(S64);
     485        SCALE_READ_IN_CASE(U8);
     486        SCALE_READ_IN_CASE(U16);
     487        SCALE_READ_IN_CASE(U32);
     488        SCALE_READ_IN_CASE(U64);
     489        SCALE_READ_IN_CASE(F32);
     490        SCALE_READ_IN_CASE(F64);
     491      default:
     492          psAbort("Should never get here: type %x should be integer", inType);
     493    }
     494
     495    return out;
     496}
     497#endif
     498
     499// Convert an image to the desired BITPIX
     500static psImage *convertImageWrite(double *bscale, // Scaling applied
     501                                  double *bzero, // Zero point applied
     502                                  psFitsFloat *floatType, // Type of custom floating-point
     503                                  psFits *fits, // FITS file pointer
     504                                  const psImage *image, // Current type
     505                                  psRandom *rng, // Random number generator
     506                                  bool newScaleZero // Determine a new BSCALE and BZERO?
     507    )
     508{
     509    assert(bscale);
     510    assert(bzero);
     511    assert(floatType);
     512    assert(fits);
     513    assert(image);
     514
     515    *bscale = 0.0;
     516    *bzero = 0.0;
     517    *floatType = PS_FITS_FLOAT_NONE;
     518
     519    // Custom floating-point
     520    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix &&
     521        fits->floatType != PS_FITS_FLOAT_NONE) {
     522        *floatType = fits->floatType;
     523        return psFitsFloatImageWrite(image, fits->floatType);
     524    }
     525
     526    // Automatically select what we're given
     527    if (fits->bitpix == 0) {
     528        return psMemIncrRefCounter((psImage*)image); // Casting away const
     529    }
     530
     531    // Quantise floating-point images
     532    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) {
     533        if (newScaleZero) {
     534            return scaleImageDetermineWrite(bscale, bzero, (psImage*)image, fits->bitpix, rng);
     535        }
     536        // Get the current BSCALE and BZERO
     537        int status = 0;                 // Status of cfitsio
     538        if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
     539            psFitsError(status, true, "Unable to read header.");
     540            return NULL;
     541        }
     542        status = 0;
     543        if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
     544            psFitsError(status, true, "Unable to read header.");
     545            return NULL;
     546        }
     547        status = 0;
     548        if (*bscale == 0.0) {
     549            psError(PS_ERR_IO, true,
     550                    "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
     551            return NULL;
     552        }
     553        return scaleImageWrite((psImage*)image, fits->bitpix, *bscale, *bzero, rng);
     554    }
     555
     556    // Choose the appropriate output type, given the input type and desired bits per pixel
     557#define CONVERT_TYPE_INT_CASE(OUTTYPE, INTYPE, BITPIX) \
     558  case BITPIX: \
     559    OUTTYPE = PS_IS_PSELEMTYPE_UNSIGNED(INTYPE) ? PS_TYPE_U##BITPIX : PS_TYPE_S##BITPIX; \
     560    break;
     561#define CONVERT_TYPE_FLOAT_CASE(OUTTYPE, BITPIX) \
     562  case -BITPIX: /* Note the use of the negative sign */ \
     563    OUTTYPE = PS_TYPE_F##BITPIX; \
     564    break;
     565
     566    *bscale = 1.0;
     567    *bzero = 0.0;
     568    psElemType inType = image->type.type; // Type for input image
     569    psElemType outType;                 // Type for output image
     570    switch (fits->bitpix) {
     571        CONVERT_TYPE_INT_CASE(outType, inType, 8);
     572        CONVERT_TYPE_INT_CASE(outType, inType, 16);
     573        CONVERT_TYPE_INT_CASE(outType, inType, 32);
     574        CONVERT_TYPE_INT_CASE(outType, inType, 64);
     575        CONVERT_TYPE_FLOAT_CASE(outType, 32);
     576        CONVERT_TYPE_FLOAT_CASE(outType, 64);
     577      default:
     578        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64",
     579                fits->bitpix);
     580        return NULL;
     581    }
     582
     583    if (outType == inType) {
     584        return psMemIncrRefCounter((psImage*)image);
     585    }
     586    return psImageCopy(NULL, image, outType);
    213587}
    214588
    … …  
    232606    if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel,
    233607                         info->increment, NULL, output->data.V[0], &anynull, &status) != 0) {
    234         char fitsErr[80];               // CFITSIO error message string
    235         (void)fits_get_errstatus(status, fitsErr);
    236         psError(PS_ERR_IO, true,
    237                 _("Reading FITS file failed. CFITSIO Error: %s"),
    238                 fitsErr);
    239         return false;
    240     }
     608        psFitsError(status, true, "Reading FITS file failed.");
     609        return false;
     610    }
     611
     612    // No need to apply the BSCALE, BZERO because cfitsio does this for us
    241613
    242614    return true;
    … …  
    259631    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z);
    260632
    261     psImage *output = psImageAlloc(info->lastPixel[0] - info->firstPixel[0] + 1,
    262                                    info->lastPixel[1] - info->firstPixel[1] + 1,
    263                                    info->psDatatype);
    264 
    265     if (!fitsReadImage(output, fits, info)) {
     633    // Size of image
     634    int numCols = info->lastPixel[0] - info->firstPixel[0] + 1;
     635    int numRows = info->lastPixel[1] - info->firstPixel[1] + 1;
     636
     637    psImage *inImage = psImageAlloc(numCols, numRows, info->psDatatype); // Image to read in
     638
     639    psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point
     640    psImage *outImage = (floatType == PS_FITS_FLOAT_NONE ? psMemIncrRefCounter(inImage) :
     641                         psImageAlloc(numCols, numRows, psFitsFloatImageType(floatType))); // Output image
     642
     643    if (!fitsReadImage(inImage, fits, info)) {
    266644        psFree(info);
    267         psFree(output);
     645        psFree(inImage);
    268646        return NULL;
    269647    }
    270 
    271648    psFree(info);
    272     return output;
    273 }
    274 
    275 psImage* psFitsReadImageBuffer(psImage *output, // Output image buffer
    276                                const psFits *fits,    // the psFits object
     649
     650    if (floatType != PS_FITS_FLOAT_NONE) {
     651        outImage = psFitsFloatImageRead(outImage, inImage, floatType);
     652    }
     653    psFree(inImage);
     654
     655    return outImage;
     656}
     657
     658psImage* psFitsReadImageBuffer(psImage *outImage, // Output image buffer
     659                               const psFits *fits, // the psFits object
    277660                               psRegion region, // the region in the FITS image to read
    278661                               int z           // the z-plane in the FITS image cube to read
    … …  
    282665    PS_ASSERT_INT_NONNEGATIVE(z, NULL);
    283666
    284     if (output && output->parent) {
     667    if (outImage && outImage->parent) {
    285668        psError(PS_ERR_IO, true, "Unable to read into a buffer for a child image.\n");
    286669        return NULL;
    … …  
    295678    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z);
    296679
    297     output = psImageRecycle(output, info->lastPixel[0] - info->firstPixel[0] + 1,
    298                             info->lastPixel[1] - info->firstPixel[1] + 1,
    299                             info->psDatatype);
    300 
    301     if (!fitsReadImage(output, fits, info)) {
     680    // Size of image
     681    int numCols = info->lastPixel[0] - info->firstPixel[0] + 1;
     682    int numRows = info->lastPixel[1] - info->firstPixel[1] + 1;
     683
     684    psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point
     685    psImage *inImage;                   // Image to read in
     686    if (floatType == PS_FITS_FLOAT_NONE) {
     687        inImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype);
     688        outImage = psMemIncrRefCounter(inImage);
     689    } else {
     690        inImage = psImageAlloc(numCols, numRows, info->psDatatype);
     691        outImage = psImageRecycle(outImage, numCols, numRows, psFitsFloatImageType(floatType));
     692    }
     693
     694    if (!fitsReadImage(inImage, fits, info)) {
    302695        psFree(info);
    303         psFree(output);
     696        psFree(inImage);
     697        psFree(outImage);
    304698        return NULL;
    305699    }
    306 
    307700    psFree(info);
    308     return output;
     701
     702    if (floatType != PS_FITS_FLOAT_NONE) {
     703        outImage = psFitsFloatImageRead(outImage, inImage, floatType);
     704    }
     705    psFree(inImage);
     706
     707    return outImage;
    309708}
    310709
    311710bool psFitsWriteImage(psFits* fits,
    312                       const psMetadata* header,
     711                      psMetadata* header,
    313712                      const psImage* input,
    314713                      int numZPlanes,
    … …  
    324723}
    325724
    326 bool psFitsInsertImage(psFits* fits,
    327                        const psMetadata* header,
    328                        const psImage* input,
    329                        int numZPlanes,
    330                        const char* extname,
    331                        bool after)
     725bool psFitsInsertImage(psFits* fits, psMetadata* header, const psImage* image, int numZPlanes,
     726                       const char* extname, bool after)
    332727{
    333728    PS_ASSERT_FITS_NON_NULL(fits, false);
    334729    PS_ASSERT_FITS_WRITABLE(fits, false);
    335     PS_ASSERT_IMAGE_NON_NULL(input, false);
    336 
    337     int numCols = input->numCols;       // Number of columns for image
    338     int numRows = input->numRows;       // Number of rows for image
     730    PS_ASSERT_IMAGE_NON_NULL(image, false);
     731
     732    int numCols = image->numCols;       // Number of columns for image
     733    int numRows = image->numRows;       // Number of rows for image
    339734    int status = 0;                     // Status from cfitsio
     735
     736    double bscale = 0.0, bzero = 0.0;   // Scale and zero point to put in header (*already* applied to data)
     737    psFitsFloat floatType;              // Custom floating-point convention type
     738    psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, image,
     739                                           NULL, true); // Image to write out
     740    if (!diskImage) {
     741        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
     742        return false;
     743    }
    340744
    341745    // determine the FITS-equivalent parameters
    342746    int bitPix;                         // Bits per pixel
    343     double bZero;                       // Zero offset
     747    double cfitsioBzero = 0.0;          // Zero point for cfitsio to apply
    344748    int dataType;                       // cfitsio data type
    345     if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) {
    346         return false;
    347     }
     749    if (!p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) {
     750        psFree(diskImage);
     751        return false;
     752    }
     753    if (cfitsioBzero != 0.0) {
     754        assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash
     755        bscale = 1.0;
     756        bzero = cfitsioBzero;
     757    }
     758    assert(bitPix == fits->bitpix || fits->bitpix == 0);
    348759
    349760    int naxis = 3;                      // Number of axes
    … …  
    377788    }
    378789
     790    // write the header, if any.
    379791    if (header && !psFitsWriteHeader(fits, header)) {
    380792        psError(PS_ERR_IO, false, "Unable to write FITS header.\n");
    381         return false;
    382     }
    383 
    384     if (bZero != 0) {        // set the bscale/bzero
    385         fits_write_key_dbl(fits->fd, "BZERO", bZero, 12, "Pixel Value Offset", &status);
    386         fits_write_key_dbl(fits->fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status);
    387         fits_set_bscale(fits->fd, 1.0, bZero, &status);
    388     }
    389 
    390     // write the header, if any.
     793        psFree(diskImage);
     794        return false;
     795    }
     796
     797    // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
     798    // an unsigned integer type).  In all other cases, we have already converted the image to use the
     799    // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
     800    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
     801
     802    if (bscale != 0.0) {
     803        fits_write_key_dbl(fits->fd, "BZERO", bzero, 12,
     804                           "Scaling: TRUE = BZERO + BSCALE * DISK", &status);
     805        fits_write_key_dbl(fits->fd, "BSCALE", bscale, 12,
     806                           "Scaling: TRUE = BZERO + BSCALE * DISK", &status);
     807        if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) {
     808            psFree(diskImage);
     809            return false;
     810        }
     811    }
     812
     813    if (floatType != PS_FITS_FLOAT_NONE) {
     814        psFitsFloatImageSet(fits, floatType);
     815    }
     816
    391817    if (extname && strlen(extname) > 0) {
    392         psFitsSetExtName(fits,extname);
    393     }
    394 
    395     if (input->parent == NULL) { // if no parent, assume that the image data is contiguous
    396         fits_write_img(fits->fd,
    397                        dataType,              // datatype
    398                        1,                     // writing to the first z-plane
    399                        numCols*numRows,       // number of elements to write, i.e., the whole image
    400                        input->data.V[0],      // the data
    401                        &status);
    402     } else { // image data may not be contiguous; write one row at a time
     818        psFitsSetExtName(fits, extname);
     819    }
     820
     821    if (image->parent == NULL) {
     822        // if no parent, assume that the image data is contiguous
     823        fits_write_img(fits->fd, dataType, 1, numCols*numRows, diskImage->data.V[0], &status);
     824    } else {
     825        // image data may not be contiguous; write one row at a time
    403826        int firstPixel = 1;
    404827        for (int row = 0; row < numRows; row++) {
    405             fits_write_img(fits->fd,
    406                            dataType,          // datatype
    407                            firstPixel,
    408                            numCols,           // number of elements to write, i.e., one row's worth
    409                            input->data.V[row],// the raw row data
    410                            &status);
    411             firstPixel += numCols;  // move to next row
    412         }
    413     }
    414 
    415     if (status != 0) {
    416         char fitsErr[MAX_STRING_LENGTH];
    417         (void)fits_get_errstatus(status, fitsErr);
    418         psError(PS_ERR_IO, true,
    419                 _("Could not write data to file. CFITSIO Error: %s"),
    420                 fitsErr);
     828            fits_write_img(fits->fd, dataType, firstPixel, numCols, diskImage->data.V[row], &status);
     829            firstPixel += numCols;
     830        }
     831    }
     832
     833    psFree(diskImage);
     834
     835    if (psFitsError(status, true, "Could not write image to file.")) {
    421836        return false;
    422837    }
    … …  
    426841}
    427842
    428 bool psFitsUpdateImage(psFits* fits,
    429                        const psImage* input,
    430                        int x0,
    431                        int y0,
    432                        int z)
     843bool psFitsUpdateImage(psFits* fits, const psImage* input, int x0, int y0, int z)
    433844{
    434845    PS_ASSERT_FITS_NON_NULL(fits, false);
    … …  
    440851    // check to see if we are positioned on an image HDU
    441852    int hduType;
    442     if ( fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
    443         char fitsErr[MAX_STRING_LENGTH];
    444         (void)fits_get_errstatus(status, fitsErr);
    445         psError(PS_ERR_IO, true,
    446                 _("Could not determine the HDU type. CFITSIO Error: %s"),
    447                 fitsErr);
     853    if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
     854        psFitsError(status, true, "Could not determine the HDU type.");
    448855        return NULL;
    449856    }
    450857    if (hduType != IMAGE_HDU) {
    451         psError(PS_ERR_IO, true,
    452                 _("Current FITS HDU type must be an image."));
     858        psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image."));
    453859        return NULL;
    454860    }
    … …  
    457863    int numRows = input->numRows;
    458864
     865    double bscale = 0.0, bzero = 0.0;   // Scale and zero point to put in header (*already* applied to data)
     866    psFitsFloat floatType;              // Custom floating-point convention type
     867    psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, input,
     868                                           NULL, false); // Image to write out
     869    if (!diskImage) {
     870        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
     871        return false;
     872    }
     873
    459874    // determine the FITS-equivalent parameters
    460     int bitPix;
    461     double bZero;
    462     int dataType;
    463     if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) {
    464         return false;
    465     }
     875    int bitPix;                         // Bits per pixel
     876    double cfitsioBzero = 0.0;          // Zero point for cfitsio to apply
     877    int dataType;                       // cfitsio data type
     878    if (! p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) {
     879        psFree(diskImage);
     880        return false;
     881    }
     882    if (cfitsioBzero != 0.0) {
     883        assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash
     884        bscale = 1.0;
     885        bzero = cfitsioBzero;
     886    }
     887    assert(bitPix == fits->bitpix || fits->bitpix == 0);
    466888
    467889    //check to see if the HDU has the same datatype
    … …  
    476898        char* fitsTypeStr;
    477899        char* imageTypeStr;
    478         PS_TYPE_NAME(fitsTypeStr,fileBitpix);
    479         PS_TYPE_NAME(imageTypeStr,input->type.type);
    480         psError(PS_ERR_IO, true,
    481                 _("Can not update a %s image given a %s image."),
    482                 fitsTypeStr, imageTypeStr);
     900        PS_TYPE_NAME(fitsTypeStr, fileBitpix);
     901        PS_TYPE_NAME(imageTypeStr, input->type.type);
     902        psError(PS_ERR_IO, true, _("Can not update a %s image given a %s image."), fitsTypeStr, imageTypeStr);
     903        psFree(diskImage);
    483904        return false;
    484905    }
    … …  
    487908    if (z >= nAxes[2]) {
    488909        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    489                 _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."),
    490                 nAxes[2], z);
     910                _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."), nAxes[2], z);
     911        psFree(diskImage);
    491912        return false;
    492913    }
    … …  
    505926
    506927    if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] ||
    507             firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
    508             lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
    509             lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
     928        firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
     929        lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
     930        lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
    510931        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    511932                "Input image [size of %ix%i] at position (%i,%i) does not all lay in the %lix%li FITS image.",
    512                 numCols, numRows,
    513                 x0, y0,
    514                 nAxes[0], nAxes[1]);
    515         return false;
    516     }
    517 
    518     fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, input->data.V[0], &status);
    519 
    520     if ( status != 0) {
    521         char fitsErr[MAX_STRING_LENGTH];
    522         (void)fits_get_errstatus(status, fitsErr);
    523         psError(PS_ERR_IO, true,
    524                 _("Could not write data to file. CFITSIO Error: %s"),
    525                 fitsErr);
     933                numCols, numRows, x0, y0, nAxes[0], nAxes[1]);
     934        psFree(diskImage);
     935        return false;
     936    }
     937
     938    // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
     939    // an unsigned integer type).  In all other cases, we have already converted the image to use the
     940    // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
     941    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
     942
     943    fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, diskImage->data.V[0], &status);
     944
     945    psFree(diskImage);
     946
     947    if (psFitsError(status, true, "Could not write data to file.")) {
    526948        return false;
    527949    }
    … …  
    537959    long nAxes[3];                      // Number of pixels on each axis
    538960    int status = 0;                     // cfitsio status value
    539     char fitsErr[80] = "";              // CFITSIO error message string
    540961
    541962    // Some of this replicates what is in psFitsReadImage, so it's a little inefficient.  But it saves
    … …  
    549970
    550971    if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) {
    551         (void)fits_get_errstatus(status, fitsErr);
    552         psError(PS_ERR_IO, true,
    553                 _("Could not determine image dimensions. CFITSIO Error: %s"),
    554                 fitsErr);
     972        psFitsError(status, true, "Could not determine image dimensions.");
    555973        return NULL;
    556974    }
    … …  
    563981    if (nAxis == 3) {
    564982        if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) {
    565             (void)fits_get_errstatus(status, fitsErr);
    566             psError(PS_ERR_IO, true,
    567                     _("Could not determine image size. CFITSIO Error: %s"),
    568                     fitsErr);
     983            psFitsError(status, true, "Could not determine image size.");
    569984            return NULL;
    570985        }
    … …  
    579994
    580995    // Bad dimensionality
    581     psError(PS_ERR_IO, true,
    582             _("Image number of dimensions, %d, is not valid."
    583               " Only two or three dimensions supported for FITS I/O."), nAxis);
     996    psError(PS_ERR_IO, true, _("Image number of dimensions, %d, is not supported."), nAxis);
    584997    return NULL;
    585998}

  • trunk/psLib/src/fits/psFitsImage.h

    r15335 r15630  
    44 * @author Robert DeSonia, MHPCC
    55 *
    6  * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-10-19 23:52:39 $
     6 * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2007-11-16 01:04:56 $
    88 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    99 */
    … …  
    5252bool psFitsWriteImage(
    5353    psFits* fits,                      ///< the psFits object
    54     const psMetadata* header,           ///< header items for the new HDU.  Can be NULL.
     54    psMetadata* header,                 ///< header items for the new HDU.  Can be NULL.
    5555    const psImage* input,              ///< the image to output
    5656    int depth,                         ///< the number of z-planes of the FITS image data cube
    … …  
    6565bool psFitsInsertImage(
    6666    psFits* fits,                      ///< the psFits object
    67     const psMetadata* header,           ///< header items for the new HDU.  Can be NULL.
     67    psMetadata* header,                 ///< header items for the new HDU.  Can be NULL.
    6868    const psImage* input,              ///< the image to output
    6969    int depth,                         ///< the number of z-planes of the FITS image data cube

  • trunk/psLib/src/pslib_strict.h

    r14929 r15630  
    99*  @author Eric Van Alst, MHPCC
    1010*
    11 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2007-09-20 23:57:31 $
     11*  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2007-11-16 01:04:56 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4040#include "psFitsImage.h"
    4141#include "psFitsTable.h"
     42#include "psFitsFloat.h"
     43#include "psFitsFloatFile.h"
    4244
    4345//#include "psXML.h"
Note: See TracChangeset for help on using the changeset viewer.