Index: trunk/psLib/src/fits/Makefile.am =================================================================== --- trunk/psLib/src/fits/Makefile.am (revision 15627) +++ trunk/psLib/src/fits/Makefile.am (revision 15630) @@ -8,6 +8,7 @@ psFitsHeader.c \ psFitsImage.c \ - psFitsTable.c -# psFitsFloat.c + psFitsTable.c \ + psFitsFloat.c \ + psFitsFloatFile.c EXTRA_DIST = fits.i @@ -17,6 +18,7 @@ psFitsHeader.h \ psFitsImage.h \ - psFitsTable.h -# psFitsFloat.h + psFitsTable.h \ + psFitsFloat.h \ + psFitsFloatFile.h CLEANFILES = *~ *.bb *.bbg *.da Index: trunk/psLib/src/fits/psFits.h =================================================================== --- trunk/psLib/src/fits/psFits.h (revision 15627) +++ trunk/psLib/src/fits/psFits.h (revision 15630) @@ -4,6 +4,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-19 23:52:39 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-16 01:04:56 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -21,4 +21,5 @@ #include "psMetadata.h" #include "psImage.h" +#include "psFitsFloat.h" #include "psErrorCodes.h" @@ -57,8 +58,8 @@ struct { bool compression; ///< Compression convention: handling of compressed images - bool psBitpix; ///< Pan-STARRS BITPIX for floating-point conversion + bool psBitpix; ///< Custom floating-point image } 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 + int bitpix; ///< Desired BITPIX for output images; 0 to use as provided + psFitsFloat floatType; ///< Desired custom floating-point for output images } psFits; Index: trunk/psLib/src/fits/psFitsFloat.c =================================================================== --- trunk/psLib/src/fits/psFitsFloat.c (revision 15627) +++ trunk/psLib/src/fits/psFitsFloat.c (revision 15630) @@ -4,4 +4,6 @@ #include +#include +#include #include "psAbort.h" @@ -9,11 +11,15 @@ #include "psError.h" #include "psImage.h" -#include "psMetadata.h" -#include "psRandom.h" #include "psFits.h" +#include "psTrace.h" +#include "psMemory.h" #include "psFitsFloat.h" -#define BIAS_FLOAT_16_0 10 // Exponent bias +#define BIAS_FLOAT_16_0 10 // Exponent bias for FLOAT_16_0 + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// union float_16_0 { @@ -27,5 +33,4 @@ } f16; }; - @@ -57,113 +62,20 @@ } +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// -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) +psImage *psFitsFloatImageWrite(const psImage *image, psFitsFloat 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 + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); psImage *output = NULL; // Output image, to return - psString convName = NULL; // Convention name switch (type) { + case PS_FITS_FLOAT_NONE: + // No conversion to be performed + return psMemIncrRefCounter((psImage*)image); // Casting away "const" case PS_FITS_FLOAT_16_0: output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S16); // Output image @@ -173,5 +85,4 @@ } } - convName = psStringCopy("FLOAT_16_0"); break; default: @@ -180,38 +91,65 @@ } - psMetadataAddStr(header, PS_LIST_TAIL, "PSBITPIX", PS_META_REPLACE, "Pan-STARRS bit encoding", convName); - return output; } -psImage *psFitsFloatImageRead(const psFits *fits, - const psImage *image - ) + +psImage *psFitsFloatImageRead(psImage *out, const psImage *in, psFitsFloat type) { - PS_ASSERT_FITS_NON_NULL(fits, NULL); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_NON_NULL(in, 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); + psElemType elem = psFitsFloatImageType(type); // Type for elements + if (elem == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise convention: %x", type); + return NULL; } - psImage *output = NULL; // Output image, to return + int numCols = in->numCols, numRows = in->numRows; // Size of image - 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]); + out = psImageRecycle(out, numCols, numRows, elem); + + switch (type) { + case PS_FITS_FLOAT_16_0: + if (in->type.type != PS_TYPE_S16) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Convention claims to be FLOAT_16_0, but image type is not S16."); + return NULL; + } + assert(out->type.type == PS_TYPE_F32); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + out->data.F32[y][x] = convertF32fromFloat16_0(in->data.S16[y][x]); } } - return output; + return out; + case PS_FITS_FLOAT_NONE: + default: + psAbort("Should be unreachable"); + } + return NULL; +} + +psElemType psFitsFloatImageType(psFitsFloat type) +{ + switch (type) { + case PS_FITS_FLOAT_16_0: + return PS_TYPE_F32; + case PS_FITS_FLOAT_NONE: + default: + return 0; // Doesn't correspond to ANY type --- should flag a real error + } +} + + +psFitsFloat psFitsFloatTypeFromString(const char *string) +{ + PS_ASSERT_STRING_NON_EMPTY(string, PS_FITS_FLOAT_NONE); + + if (strcmp(string, "FLOAT_16_0") == 0) { + return PS_FITS_FLOAT_16_0; } - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise PSBITPIX convention name: %s", convName); - return NULL; - - return output; + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Unable to recognise custom floating-point convention name: %s", string); + return PS_FITS_FLOAT_NONE; } Index: trunk/psLib/src/fits/psFitsFloat.h =================================================================== --- trunk/psLib/src/fits/psFitsFloat.h (revision 15627) +++ trunk/psLib/src/fits/psFitsFloat.h (revision 15630) @@ -2,19 +2,32 @@ #define PS_FITS_FLOAT_H - +#include +#include #include -#include /// Type of custom floating point typedef enum { + PS_FITS_FLOAT_NONE, ///< No conversion to be performed 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 +psImage *psFitsFloatImageWrite(const psImage *image, ///< Image to convert + psFitsFloat type ///< Custom floating point type ); +/// Convert the custom floating-point image to a floating-point image +psImage *psFitsFloatImageRead(psImage *out, ///< Output image, or NULL + const psImage *in, ///< Image to convert + psFitsFloat type ///< Custom floating point type + ); + +/// Return the appropriate element type for a custom floating-point +psElemType psFitsFloatImageType(psFitsFloat type ///< Custom floating-point type + ); + +/// Return the custom floating-point type from a string description +psFitsFloat psFitsFloatTypeFromString(const char *string ///< String with name of type + ); #endif Index: trunk/psLib/src/fits/psFitsFloatFile.c =================================================================== --- trunk/psLib/src/fits/psFitsFloatFile.c (revision 15630) +++ trunk/psLib/src/fits/psFitsFloatFile.c (revision 15630) @@ -0,0 +1,62 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "psError.h" +#include "psFits.h" +#include "psFitsFloat.h" +#include "psMemory.h" + +#include "psFitsFloatFile.h" + +bool psFitsFloatImageSet(const psFits *fits, psFitsFloat type) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_FITS_WRITABLE(fits, false); + + char *convName; // Convention name + + switch (type) { + case PS_FITS_FLOAT_NONE: + return true; + case PS_FITS_FLOAT_16_0: + convName = "FLOAT_16_0"; + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise convention: %x", type); + return false; + } + + int status = 0; // Status from cfitsio + fits_write_key_str(fits->fd, "PSBITPIX", convName, "Custom floating point convention name", &status); + if (psFitsError(status, true, "Could not write PSBITPIX header to file.")) { + return false; + } + return true; +} + + +psFitsFloat psFitsFloatImageCheck(const psFits *fits) +{ + PS_ASSERT_FITS_NON_NULL(fits, PS_FITS_FLOAT_NONE); + + if (!fits->conventions.psBitpix) { + return PS_FITS_FLOAT_NONE; + } + + int status = 0; // Status of CFITSIO calls + char convName[FLEN_CARD]; // Convention name for custom floating-point + if (fits_read_key_str(fits->fd, "PSBITPIX", convName, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + return PS_FITS_FLOAT_NONE; + } + + if (!convName || status == KEY_NO_EXIST) { + return PS_FITS_FLOAT_NONE; + } + + return psFitsFloatTypeFromString(convName); +} Index: trunk/psLib/src/fits/psFitsFloatFile.h =================================================================== --- trunk/psLib/src/fits/psFitsFloatFile.h (revision 15630) +++ trunk/psLib/src/fits/psFitsFloatFile.h (revision 15630) @@ -0,0 +1,20 @@ +// This file exists to resolve the co-dependency problem of placing these function definitions in +// psFitsFloat.h --- because they refer to psFits, and psFits refers to psFitsFloat, we would end up with a +// co-dependency problem if these function declarations were to be placed in psFitsFloat.h. +#ifndef PS_FITS_FLOAT_FILE_H +#define PS_FITS_FLOAT_FILE_H + +#include +#include + +/// Set a flag in the FITS header of the current extension that the image is a custom floating-point +bool psFitsFloatImageSet(const psFits *fits, ///< FITS file + psFitsFloat type ///< Custom floating-point type + ); + +/// Check if the current extension contains a custom floating-point, returning the appropriate type +psFitsFloat psFitsFloatImageCheck(const psFits *fits ///< FITS file + ); + + +#endif Index: trunk/psLib/src/fits/psFitsHeader.c =================================================================== --- trunk/psLib/src/fits/psFitsHeader.c (revision 15627) +++ trunk/psLib/src/fits/psFitsHeader.c (revision 15630) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-19 23:52:39 $ + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-16 01:04:56 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -44,5 +44,5 @@ static const char *noWriteFitsKeys[] = { "SIMPLE", "XTENSION", "BITPIX", "NAXIS", "EXTNAME", "BSCALE", "BZERO", "TFIELDS", "PCOUNT", "GCOUNT", "ZIMAGE", "ZBITPIX", - "ZCMPTYPE", NULL}; + "ZCMPTYPE", "PSBITPIX", NULL}; // List of the start of FITS header keys not to write (handled by cfitsio); NULL-terminated Index: trunk/psLib/src/fits/psFitsImage.c =================================================================== --- trunk/psLib/src/fits/psFitsImage.c (revision 15627) +++ trunk/psLib/src/fits/psFitsImage.c (revision 15630) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-19 23:52:39 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-16 01:04:56 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -21,19 +21,25 @@ #include -#include "psFits.h" +#include "psAbort.h" +#include "psType.h" #include "psAssert.h" -#include "psFitsImage.h" #include "psError.h" - -#include "psImageStructManip.h" -#include "psMemory.h" #include "psString.h" #include "psLogMsg.h" #include "psTrace.h" #include "psVector.h" - +#include "psRandom.h" +#include "psImageStructManip.h" + +#include "psFits.h" +#include "psFitsFloat.h" +#include "psFitsFloatFile.h" #include "psFitsHeader.h" -#define MAX_STRING_LENGTH 256 // maximum length string for FITS routines +#include "psMemory.h" + +#include "psFitsImage.h" + +#define MAX_STRING_LENGTH 256 // maximum length string for FITS routines // Information required to read a FITS file @@ -49,8 +55,10 @@ } p_psFitsReadInfo; -static p_psFitsReadInfo *p_psFitsReadInfoAlloc(const psFits *fits, // The FITS file handle - psRegion region, // Region to read - int z // z-plane to read in cube - ) +// Read the vital statistics of this FITS image, in preparation for reading the image +static p_psFitsReadInfo *p_psFitsReadInfoAlloc( + const psFits *fits, // The FITS file handle + psRegion region, // Region to read + int z // z-plane to read in cube + ) { PS_ASSERT_FITS_NON_NULL(fits, NULL); @@ -61,18 +69,13 @@ int status = 0; // CFITSIO status - char fitsErr[80]; // CFITSIO error message string // check to see if we even are positioned on an image HDU int hduType; // Type of HDU if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) { - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine the HDU type. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine the HDU type."); goto bad; } if (hduType != IMAGE_HDU) { - psError(PS_ERR_IO, true, - _("Current FITS HDU type must be an image.")); + psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image.")); goto bad; } @@ -80,8 +83,5 @@ // Get the data type 'bitPix' from the FITS image if (fits_get_img_equivtype(fits->fd, &info->bitPix, &status) != 0) { - fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine image data type. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine image data type."); goto bad; } @@ -89,8 +89,5 @@ /* Get the dimensions 'nAxis' from the FITS image */ if (fits_get_img_dim(fits->fd, &info->nAxis, &status) != 0) { - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine image dimensions. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine image dimensions."); goto bad; } @@ -99,6 +96,5 @@ if ((info->nAxis < 2) || (info->nAxis > 3)) { psError(PS_ERR_IO, true, - _("Image number of dimensions, %d, is not valid. Only two or three dimensions supported for FITS I/O."), - info->nAxis); + _("Image number of dimensions, %d, is not supported."), info->nAxis); goto bad; } @@ -106,8 +102,5 @@ /* Get the Image size from the FITS file */ if (fits_get_img_size(fits->fd, info->nAxis, info->nAxes, &status) != 0) { - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine image size. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine image size."); goto bad; } @@ -133,48 +126,83 @@ info->increment[2] = 1; - switch (info->bitPix) { - case BYTE_IMG: - info->psDatatype = PS_TYPE_U8; - info->fitsDatatype = TBYTE; - break; - case SBYTE_IMG: - info->psDatatype = PS_TYPE_S8; - info->fitsDatatype = TSBYTE; - break; - case USHORT_IMG: - info->psDatatype = PS_TYPE_U16; - info->fitsDatatype = TUSHORT; - break; - case SHORT_IMG: - info->psDatatype = PS_TYPE_S16; - info->fitsDatatype = TSHORT; - break; - case ULONG_IMG: - info->psDatatype = PS_TYPE_U32; - info->fitsDatatype = TUINT; - break; - case LONG_IMG: - info->psDatatype = PS_TYPE_S32; - info->fitsDatatype = TINT; - break; - case LONGLONG_IMG: - info->psDatatype = PS_TYPE_S64; - info->fitsDatatype = TLONGLONG; - break; - case FLOAT_IMG: - info->psDatatype = PS_TYPE_F32; - info->fitsDatatype = TFLOAT; - break; - case DOUBLE_IMG: - info->psDatatype = PS_TYPE_F64; - info->fitsDatatype = TDOUBLE; - break; - default: - psError(PS_ERR_IO, true, - _("FITS image type, BITPIX=%d, is not supported."), - info->bitPix); + // Check scale and zero + double bscale = 0.0, bzero = 0.0; // Scale and zero point + if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); goto bad; } - + status = 0; + if (fits_read_key_dbl(fits->fd, "BZERO", &bzero, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + goto bad; + } + status = 0; + + if ((bscale != 0.0 && bscale != 1.0) || bzero != (int)bzero) { + // It's a floating-point image that's been quantised + // cfitsio will apply the scale and zero point for us if we choose the correct data type + switch (info->bitPix) { + case BYTE_IMG: + case SBYTE_IMG: + case USHORT_IMG: + case SHORT_IMG: + case ULONG_IMG: + case LONG_IMG: + case FLOAT_IMG: + info->psDatatype = PS_TYPE_F32; + info->fitsDatatype = TFLOAT; + break; + case LONGLONG_IMG: + case DOUBLE_IMG: + info->psDatatype = PS_TYPE_F64; + info->fitsDatatype = TDOUBLE; + break; + default: + psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix); + goto bad; + } + } else { + switch (info->bitPix) { + case BYTE_IMG: + info->psDatatype = PS_TYPE_U8; + info->fitsDatatype = TBYTE; + break; + case SBYTE_IMG: + info->psDatatype = PS_TYPE_S8; + info->fitsDatatype = TSBYTE; + break; + case USHORT_IMG: + info->psDatatype = PS_TYPE_U16; + info->fitsDatatype = TUSHORT; + break; + case SHORT_IMG: + info->psDatatype = PS_TYPE_S16; + info->fitsDatatype = TSHORT; + break; + case ULONG_IMG: + info->psDatatype = PS_TYPE_U32; + info->fitsDatatype = TUINT; + break; + case LONG_IMG: + info->psDatatype = PS_TYPE_S32; + info->fitsDatatype = TINT; + break; + case LONGLONG_IMG: + info->psDatatype = PS_TYPE_S64; + info->fitsDatatype = TLONGLONG; + break; + case FLOAT_IMG: + info->psDatatype = PS_TYPE_F32; + info->fitsDatatype = TFLOAT; + break; + case DOUBLE_IMG: + info->psDatatype = PS_TYPE_F64; + info->fitsDatatype = TDOUBLE; + break; + default: + psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix); + goto bad; + } + } return info; @@ -211,4 +239,350 @@ return true; +} + +// Apply the BSCALE and BZERO for an image with a "fuzz" +// 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 *scaleImageWrite(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; +} + + + +// Determine BSCALE and BZERO for an image, and generate a new image with it applied +// TRUE = BZERO + BSCALE * FITS +static psImage *scaleImageDetermineWrite(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 scaleImageWrite + ) +{ + 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 scaleImageWrite(image, bitpix, *bscale, *bzero, rng); +} + + +#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 *scaleImageRead(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 +static psImage *convertImageWrite(double *bscale, // Scaling applied + double *bzero, // Zero point applied + psFitsFloat *floatType, // Type of custom floating-point + psFits *fits, // FITS file pointer + const psImage *image, // Current type + psRandom *rng, // Random number generator + bool newScaleZero // Determine a new BSCALE and BZERO? + ) +{ + assert(bscale); + assert(bzero); + assert(floatType); + assert(fits); + assert(image); + + *bscale = 0.0; + *bzero = 0.0; + *floatType = PS_FITS_FLOAT_NONE; + + // Custom floating-point + if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix && + fits->floatType != PS_FITS_FLOAT_NONE) { + *floatType = fits->floatType; + return psFitsFloatImageWrite(image, fits->floatType); + } + + // Automatically select what we're given + if (fits->bitpix == 0) { + return psMemIncrRefCounter((psImage*)image); // Casting away const + } + + // Quantise floating-point images + if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) { + if (newScaleZero) { + return scaleImageDetermineWrite(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 scaleImageWrite((psImage*)image, fits->bitpix, *bscale, *bzero, rng); + } + + // Choose the appropriate output type, given the input type and desired bits per pixel +#define CONVERT_TYPE_INT_CASE(OUTTYPE, INTYPE, BITPIX) \ + case BITPIX: \ + OUTTYPE = PS_IS_PSELEMTYPE_UNSIGNED(INTYPE) ? PS_TYPE_U##BITPIX : PS_TYPE_S##BITPIX; \ + break; +#define CONVERT_TYPE_FLOAT_CASE(OUTTYPE, BITPIX) \ + case -BITPIX: /* Note the use of the negative sign */ \ + OUTTYPE = PS_TYPE_F##BITPIX; \ + break; + + *bscale = 1.0; + *bzero = 0.0; + psElemType inType = image->type.type; // Type for input image + psElemType outType; // Type for output image + switch (fits->bitpix) { + CONVERT_TYPE_INT_CASE(outType, inType, 8); + CONVERT_TYPE_INT_CASE(outType, inType, 16); + CONVERT_TYPE_INT_CASE(outType, inType, 32); + CONVERT_TYPE_INT_CASE(outType, inType, 64); + CONVERT_TYPE_FLOAT_CASE(outType, 32); + CONVERT_TYPE_FLOAT_CASE(outType, 64); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", + fits->bitpix); + return NULL; + } + + if (outType == inType) { + return psMemIncrRefCounter((psImage*)image); + } + return psImageCopy(NULL, image, outType); } @@ -232,11 +606,9 @@ if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel, info->increment, NULL, output->data.V[0], &anynull, &status) != 0) { - char fitsErr[80]; // CFITSIO error message string - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Reading FITS file failed. CFITSIO Error: %s"), - fitsErr); - return false; - } + psFitsError(status, true, "Reading FITS file failed."); + return false; + } + + // No need to apply the BSCALE, BZERO because cfitsio does this for us return true; @@ -259,20 +631,31 @@ p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z); - psImage *output = psImageAlloc(info->lastPixel[0] - info->firstPixel[0] + 1, - info->lastPixel[1] - info->firstPixel[1] + 1, - info->psDatatype); - - if (!fitsReadImage(output, fits, info)) { + // Size of image + int numCols = info->lastPixel[0] - info->firstPixel[0] + 1; + int numRows = info->lastPixel[1] - info->firstPixel[1] + 1; + + psImage *inImage = psImageAlloc(numCols, numRows, info->psDatatype); // Image to read in + + psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point + psImage *outImage = (floatType == PS_FITS_FLOAT_NONE ? psMemIncrRefCounter(inImage) : + psImageAlloc(numCols, numRows, psFitsFloatImageType(floatType))); // Output image + + if (!fitsReadImage(inImage, fits, info)) { psFree(info); - psFree(output); + psFree(inImage); return NULL; } - psFree(info); - return output; -} - -psImage* psFitsReadImageBuffer(psImage *output, // Output image buffer - const psFits *fits, // the psFits object + + if (floatType != PS_FITS_FLOAT_NONE) { + outImage = psFitsFloatImageRead(outImage, inImage, floatType); + } + psFree(inImage); + + return outImage; +} + +psImage* psFitsReadImageBuffer(psImage *outImage, // Output image buffer + const psFits *fits, // the psFits object psRegion region, // the region in the FITS image to read int z // the z-plane in the FITS image cube to read @@ -282,5 +665,5 @@ PS_ASSERT_INT_NONNEGATIVE(z, NULL); - if (output && output->parent) { + if (outImage && outImage->parent) { psError(PS_ERR_IO, true, "Unable to read into a buffer for a child image.\n"); return NULL; @@ -295,20 +678,36 @@ p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z); - output = psImageRecycle(output, info->lastPixel[0] - info->firstPixel[0] + 1, - info->lastPixel[1] - info->firstPixel[1] + 1, - info->psDatatype); - - if (!fitsReadImage(output, fits, info)) { + // Size of image + int numCols = info->lastPixel[0] - info->firstPixel[0] + 1; + int numRows = info->lastPixel[1] - info->firstPixel[1] + 1; + + psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point + psImage *inImage; // Image to read in + if (floatType == PS_FITS_FLOAT_NONE) { + inImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype); + outImage = psMemIncrRefCounter(inImage); + } else { + inImage = psImageAlloc(numCols, numRows, info->psDatatype); + outImage = psImageRecycle(outImage, numCols, numRows, psFitsFloatImageType(floatType)); + } + + if (!fitsReadImage(inImage, fits, info)) { psFree(info); - psFree(output); + psFree(inImage); + psFree(outImage); return NULL; } - psFree(info); - return output; + + if (floatType != PS_FITS_FLOAT_NONE) { + outImage = psFitsFloatImageRead(outImage, inImage, floatType); + } + psFree(inImage); + + return outImage; } bool psFitsWriteImage(psFits* fits, - const psMetadata* header, + psMetadata* header, const psImage* input, int numZPlanes, @@ -324,26 +723,38 @@ } -bool psFitsInsertImage(psFits* fits, - const psMetadata* header, - const psImage* input, - int numZPlanes, - const char* extname, - bool after) +bool psFitsInsertImage(psFits* fits, psMetadata* header, const psImage* image, int numZPlanes, + const char* extname, bool after) { PS_ASSERT_FITS_NON_NULL(fits, false); PS_ASSERT_FITS_WRITABLE(fits, false); - PS_ASSERT_IMAGE_NON_NULL(input, false); - - int numCols = input->numCols; // Number of columns for image - int numRows = input->numRows; // Number of rows for image + PS_ASSERT_IMAGE_NON_NULL(image, false); + + int numCols = image->numCols; // Number of columns for image + int numRows = image->numRows; // Number of rows for image int status = 0; // Status from cfitsio + + double bscale = 0.0, bzero = 0.0; // Scale and zero point to put in header (*already* applied to data) + psFitsFloat floatType; // Custom floating-point convention type + psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, image, + NULL, true); // Image to write out + if (!diskImage) { + psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format."); + return false; + } // determine the FITS-equivalent parameters int bitPix; // Bits per pixel - double bZero; // Zero offset + double cfitsioBzero = 0.0; // Zero point for cfitsio to apply int dataType; // cfitsio data type - if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) { - return false; - } + if (!p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) { + psFree(diskImage); + return false; + } + if (cfitsioBzero != 0.0) { + assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash + bscale = 1.0; + bzero = cfitsioBzero; + } + assert(bitPix == fits->bitpix || fits->bitpix == 0); int naxis = 3; // Number of axes @@ -377,46 +788,50 @@ } + // write the header, if any. if (header && !psFitsWriteHeader(fits, header)) { psError(PS_ERR_IO, false, "Unable to write FITS header.\n"); - return false; - } - - if (bZero != 0) { // set the bscale/bzero - fits_write_key_dbl(fits->fd, "BZERO", bZero, 12, "Pixel Value Offset", &status); - fits_write_key_dbl(fits->fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status); - fits_set_bscale(fits->fd, 1.0, bZero, &status); - } - - // write the header, if any. + psFree(diskImage); + return false; + } + + // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is + // 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, cfitsioBzero, &status); + + if (bscale != 0.0) { + fits_write_key_dbl(fits->fd, "BZERO", bzero, 12, + "Scaling: TRUE = BZERO + BSCALE * DISK", &status); + fits_write_key_dbl(fits->fd, "BSCALE", bscale, 12, + "Scaling: TRUE = BZERO + BSCALE * DISK", &status); + if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) { + psFree(diskImage); + return false; + } + } + + if (floatType != PS_FITS_FLOAT_NONE) { + psFitsFloatImageSet(fits, floatType); + } + if (extname && strlen(extname) > 0) { - psFitsSetExtName(fits,extname); - } - - if (input->parent == NULL) { // if no parent, assume that the image data is contiguous - fits_write_img(fits->fd, - dataType, // datatype - 1, // writing to the first z-plane - numCols*numRows, // number of elements to write, i.e., the whole image - input->data.V[0], // the data - &status); - } else { // image data may not be contiguous; write one row at a time + psFitsSetExtName(fits, extname); + } + + if (image->parent == NULL) { + // if no parent, assume that the image data is contiguous + fits_write_img(fits->fd, dataType, 1, numCols*numRows, diskImage->data.V[0], &status); + } else { + // image data may not be contiguous; write one row at a time int firstPixel = 1; for (int row = 0; row < numRows; row++) { - fits_write_img(fits->fd, - dataType, // datatype - firstPixel, - numCols, // number of elements to write, i.e., one row's worth - input->data.V[row],// the raw row data - &status); - firstPixel += numCols; // move to next row - } - } - - if (status != 0) { - char fitsErr[MAX_STRING_LENGTH]; - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not write data to file. CFITSIO Error: %s"), - fitsErr); + fits_write_img(fits->fd, dataType, firstPixel, numCols, diskImage->data.V[row], &status); + firstPixel += numCols; + } + } + + psFree(diskImage); + + if (psFitsError(status, true, "Could not write image to file.")) { return false; } @@ -426,9 +841,5 @@ } -bool psFitsUpdateImage(psFits* fits, - const psImage* input, - int x0, - int y0, - int z) +bool psFitsUpdateImage(psFits* fits, const psImage* input, int x0, int y0, int z) { PS_ASSERT_FITS_NON_NULL(fits, false); @@ -440,15 +851,10 @@ // check to see if we are positioned on an image HDU int hduType; - if ( fits_get_hdu_type(fits->fd, &hduType, &status) != 0) { - char fitsErr[MAX_STRING_LENGTH]; - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine the HDU type. CFITSIO Error: %s"), - fitsErr); + if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) { + psFitsError(status, true, "Could not determine the HDU type."); return NULL; } if (hduType != IMAGE_HDU) { - psError(PS_ERR_IO, true, - _("Current FITS HDU type must be an image.")); + psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image.")); return NULL; } @@ -457,11 +863,27 @@ int numRows = input->numRows; + double bscale = 0.0, bzero = 0.0; // Scale and zero point to put in header (*already* applied to data) + psFitsFloat floatType; // Custom floating-point convention type + psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, input, + NULL, false); // Image to write out + if (!diskImage) { + psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format."); + return false; + } + // determine the FITS-equivalent parameters - int bitPix; - double bZero; - int dataType; - if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) { - return false; - } + int bitPix; // Bits per pixel + double cfitsioBzero = 0.0; // Zero point for cfitsio to apply + int dataType; // cfitsio data type + if (! p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) { + psFree(diskImage); + return false; + } + if (cfitsioBzero != 0.0) { + assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash + bscale = 1.0; + bzero = cfitsioBzero; + } + assert(bitPix == fits->bitpix || fits->bitpix == 0); //check to see if the HDU has the same datatype @@ -476,9 +898,8 @@ char* fitsTypeStr; char* imageTypeStr; - PS_TYPE_NAME(fitsTypeStr,fileBitpix); - PS_TYPE_NAME(imageTypeStr,input->type.type); - psError(PS_ERR_IO, true, - _("Can not update a %s image given a %s image."), - fitsTypeStr, imageTypeStr); + PS_TYPE_NAME(fitsTypeStr, fileBitpix); + PS_TYPE_NAME(imageTypeStr, input->type.type); + psError(PS_ERR_IO, true, _("Can not update a %s image given a %s image."), fitsTypeStr, imageTypeStr); + psFree(diskImage); return false; } @@ -487,6 +908,6 @@ if (z >= nAxes[2]) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, - _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."), - nAxes[2], z); + _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."), nAxes[2], z); + psFree(diskImage); return false; } @@ -505,23 +926,24 @@ if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] || - firstPixel[1] < 1 || firstPixel[1] > nAxes[1] || - lastPixel[0] < 1 || lastPixel[0] > nAxes[0] || - lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) { + firstPixel[1] < 1 || firstPixel[1] > nAxes[1] || + lastPixel[0] < 1 || lastPixel[0] > nAxes[0] || + lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input image [size of %ix%i] at position (%i,%i) does not all lay in the %lix%li FITS image.", - numCols, numRows, - x0, y0, - nAxes[0], nAxes[1]); - return false; - } - - fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, input->data.V[0], &status); - - if ( status != 0) { - char fitsErr[MAX_STRING_LENGTH]; - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not write data to file. CFITSIO Error: %s"), - fitsErr); + numCols, numRows, x0, y0, nAxes[0], nAxes[1]); + psFree(diskImage); + return false; + } + + // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is + // 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, cfitsioBzero, &status); + + fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, diskImage->data.V[0], &status); + + psFree(diskImage); + + if (psFitsError(status, true, "Could not write data to file.")) { return false; } @@ -537,5 +959,4 @@ long nAxes[3]; // Number of pixels on each axis int status = 0; // cfitsio status value - char fitsErr[80] = ""; // CFITSIO error message string // Some of this replicates what is in psFitsReadImage, so it's a little inefficient. But it saves @@ -549,8 +970,5 @@ if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) { - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine image dimensions. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine image dimensions."); return NULL; } @@ -563,8 +981,5 @@ if (nAxis == 3) { if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) { - (void)fits_get_errstatus(status, fitsErr); - psError(PS_ERR_IO, true, - _("Could not determine image size. CFITSIO Error: %s"), - fitsErr); + psFitsError(status, true, "Could not determine image size."); return NULL; } @@ -579,7 +994,5 @@ // Bad dimensionality - psError(PS_ERR_IO, true, - _("Image number of dimensions, %d, is not valid." - " Only two or three dimensions supported for FITS I/O."), nAxis); + psError(PS_ERR_IO, true, _("Image number of dimensions, %d, is not supported."), nAxis); return NULL; } Index: trunk/psLib/src/fits/psFitsImage.h =================================================================== --- trunk/psLib/src/fits/psFitsImage.h (revision 15627) +++ trunk/psLib/src/fits/psFitsImage.h (revision 15630) @@ -4,6 +4,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-19 23:52:39 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-16 01:04:56 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -52,5 +52,5 @@ bool psFitsWriteImage( psFits* fits, ///< the psFits object - const psMetadata* header, ///< header items for the new HDU. Can be NULL. + psMetadata* header, ///< header items for the new HDU. Can be NULL. const psImage* input, ///< the image to output int depth, ///< the number of z-planes of the FITS image data cube @@ -65,5 +65,5 @@ bool psFitsInsertImage( psFits* fits, ///< the psFits object - const psMetadata* header, ///< header items for the new HDU. Can be NULL. + psMetadata* header, ///< header items for the new HDU. Can be NULL. const psImage* input, ///< the image to output int depth, ///< the number of z-planes of the FITS image data cube Index: trunk/psLib/src/pslib_strict.h =================================================================== --- trunk/psLib/src/pslib_strict.h (revision 15627) +++ trunk/psLib/src/pslib_strict.h (revision 15630) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-09-20 23:57:31 $ +* @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-11-16 01:04:56 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -40,4 +40,6 @@ #include "psFitsImage.h" #include "psFitsTable.h" +#include "psFitsFloat.h" +#include "psFitsFloatFile.h" //#include "psXML.h"