Index: trunk/psLib/src/fits/psFitsImage.c =================================================================== --- trunk/psLib/src/fits/psFitsImage.c (revision 15335) +++ 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; }