Index: /trunk/psLib/src/fits/Makefile.am =================================================================== --- /trunk/psLib/src/fits/Makefile.am (revision 16184) +++ /trunk/psLib/src/fits/Makefile.am (revision 16185) @@ -10,5 +10,6 @@ psFitsTable.c \ psFitsFloat.c \ - psFitsFloatFile.c + psFitsFloatFile.c \ + psFitsScale.c EXTRA_DIST = fits.i @@ -20,5 +21,6 @@ psFitsTable.h \ psFitsFloat.h \ - psFitsFloatFile.h + psFitsFloatFile.h \ + psFitsScale.h CLEANFILES = *~ *.bb *.bbg *.da Index: /trunk/psLib/src/fits/psFits.c =================================================================== --- /trunk/psLib/src/fits/psFits.c (revision 16184) +++ /trunk/psLib/src/fits/psFits.c (revision 16185) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.76 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-25 01:29:42 $ + * @version $Revision: 1.77 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-23 03:08:03 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -78,5 +78,5 @@ fitsClose(fits); } - psFree (fits->extword); + psFree(fits->options); } @@ -160,14 +160,45 @@ psFits* fits = psAlloc(sizeof(psFits)); + psMemSetDeallocator(fits, (psFreeFunc)fitsFree); + fits->fd = fptr; fits->writable = (iomode == READWRITE); - fits->extword = NULL; - fits->conventions.compression = true; - fits->conventions.psBitpix = true; - fits->bitpix = 0; - fits->floatType = PS_FITS_FLOAT_NONE; - psMemSetDeallocator(fits,(psFreeFunc)fitsFree); + + fits->options = NULL; return fits; +} + + +static void fitsOptionsFree(psFitsOptions *options) +{ + psFree(options->extword); +} + + +psFitsOptions *psFitsOptionsAlloc(void) +{ + psFitsOptions *options = psAlloc(sizeof(psFitsOptions)); // Options, to return + psMemSetDeallocator(options, (psFreeFunc)fitsOptionsFree); + + options->extword = NULL; + + options->conventions.compression = true; + options->conventions.psBitpix = true; + + options->floatType = PS_FITS_FLOAT_NONE; + + options->bitpix = 0; + + options->scaling = PS_FITS_SCALE_NONE; + options->fuzz = true; + options->bscale = 1.0; + options->bzero = 0.0; + options->mean = NAN; + options->stdev = NAN; + options->stdevBits = 4; + options->stdevNum = 5.0; + + return options; } @@ -227,11 +258,15 @@ } -bool psFitsSetExtnameWord (psFits *fits, const char *extword) +bool psFitsSetExtnameWord(psFits *fits, const char *extword) { PS_ASSERT_PTR_NON_NULL(fits, false); PS_ASSERT_PTR_NON_NULL(extword, false); - psFree (fits->extword); - fits->extword = psStringCopy (extword); + if (!fits->options) { + fits->options = psFitsOptionsAlloc(); + } + + psFree(fits->options->extword); + fits->options->extword = psStringCopy(extword); return true; } @@ -251,5 +286,6 @@ int status = 0; - if (!fits->conventions.compression && !fits->extword) { + psFitsOptions *options = fits->options; // FITS options + if (options && !options->conventions.compression && !options->extword) { // User wants to use cfitsio. Good luck to them! if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) { @@ -260,7 +296,7 @@ } - bool ignoreCI = (fits->conventions.compression && + bool ignoreCI = (options->conventions.compression && (strcmp(extname, "COMPRESSED_IMAGE") != 0)); // Ignore COMPRESSED_IMAGE extension name? - char *extword = (fits->extword ? fits->extword : "EXTNAME"); // Word to use as extension name + char *extword = (options->extword ? options->extword : "EXTNAME"); // Word to use as extension name #if 0 @@ -365,9 +401,9 @@ char name[MAX_STRING_LENGTH]; - char *extword = (fits->extword == NULL) ? defaultExtword : fits->extword; + psFitsOptions *options = fits->options; // FITS options + char *extword = (!options || !options->extword) ? defaultExtword : options->extword; if (fits_read_key_str(fits->fd, extword, name, NULL, &status) != 0) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Header keyword %s is not found"), extword); + psError(PS_ERR_BAD_PARAMETER_NULL, true, _("Header keyword %s is not found"), extword); return NULL; } @@ -382,12 +418,11 @@ int status = 0; - char *extword = (fits->extword == NULL) ? defaultExtword : fits->extword; + psFitsOptions *options = fits->options; // FITS options + char *extword = (!options || !options->extword) ? defaultExtword : options->extword; if (fits_update_key_str(fits->fd, extword, (char*)name, NULL, &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); + psError(PS_ERR_IO, true, _("Could not write data to file. CFITSIO Error: %s"), fitsErr); return false; } @@ -404,8 +439,6 @@ // move to the specified HDU - if (! psFitsMoveExtNum(fits,extnum,relative) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "Failed to delete HDU #%d", - extnum); + if (!psFitsMoveExtNum(fits, extnum, relative) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Failed to delete HDU #%d", extnum); return false; } @@ -417,7 +450,5 @@ 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); + psError(PS_ERR_IO, true, _("Could not write data to file. CFITSIO Error: %s"), fitsErr); return false; } @@ -598,12 +629,12 @@ // that is the restriction; data must be 32 or 64 bit for noise bits to be valid. if (type != PS_FITS_COMPRESS_PLIO) { - if (fits_set_noise_bits(fits->fd, noisebits, &status)) { - fits_set_compression_type(fits->fd, 0x0, &status); - char fitsErr[MAX_STRING_LENGTH]; - fits_get_errstatus(status, fitsErr); - psError(PS_ERR_BAD_FITS, true, - "Error while configuring compression. CFITSIO error: %s", fitsErr); - return false; - } + if (fits_set_noise_bits(fits->fd, noisebits, &status)) { + fits_set_compression_type(fits->fd, 0x0, &status); + char fitsErr[MAX_STRING_LENGTH]; + fits_get_errstatus(status, fitsErr); + psError(PS_ERR_BAD_FITS, true, + "Error while configuring compression. CFITSIO error: %s", fitsErr); + return false; + } } Index: /trunk/psLib/src/fits/psFits.h =================================================================== --- /trunk/psLib/src/fits/psFits.h (revision 16184) +++ /trunk/psLib/src/fits/psFits.h (revision 16185) @@ -4,6 +4,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-16 01:04:56 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-23 03:08:03 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -25,9 +25,5 @@ #include "psErrorCodes.h" -/** FITS HDU type. - * - * Enumeration for FITS HDU type. - * - */ +/// FITS HDU type. typedef enum { PS_FITS_TYPE_NONE = -1, ///< Unknown HDU type @@ -38,21 +34,25 @@ } psFitsType; +/// FITS compression type typedef enum { - PS_FITS_COMPRESS_NONE = 0, - PS_FITS_COMPRESS_GZIP, - PS_FITS_COMPRESS_RICE, - PS_FITS_COMPRESS_HCOMPRESS, - PS_FITS_COMPRESS_PLIO + PS_FITS_COMPRESS_NONE = 0, ///< No compression + PS_FITS_COMPRESS_GZIP, ///< GZIP compression (of the pixels only) + PS_FITS_COMPRESS_RICE, ///< RICE compression (of the pixels only) + PS_FITS_COMPRESS_HCOMPRESS, ///< HCOMPRESS compression (of the pixels only) + PS_FITS_COMPRESS_PLIO ///< PLIO compression (of the pixels only; appropriate for masks) } psFitsCompressionType; -/** FITS file object. - * - * This object should be considered opaque to the user; no item in this - * struct should be accessed directly. - * - */ +/// FITS scaling method: how to set BSCALE and BZERO +typedef enum { + PS_FITS_SCALE_NONE, ///< No auto-scaling to be applied (BSCALE = 1, BZERO = 0) + PS_FITS_SCALE_RANGE, ///< Auto-scale to preserve dynamic range + PS_FITS_SCALE_STDEV_POSITIVE, ///< Auto-scale to sample stdev, place mean at lower limit + PS_FITS_SCALE_STDEV_NEGATIVE, ///< Auto-scale to sample stdev, place mean at upper limit + PS_FITS_SCALE_STDEV_BOTH, ///< Auto-scale to sample stdev, place mean at middle + PS_FITS_SCALE_MANUAL ///< Manual scaling (use specified BSCALE and BZERO) +} psFitsScaling; + +/// Options for FITS I/O typedef struct { - fitsfile* fd; ///< the CFITSIO fits files handle. - bool writable; ///< Is the file writable? char *extword; ///< user-specified word to name extensions (NULL implies EXTNAME) struct { @@ -60,7 +60,23 @@ bool psBitpix; ///< Custom floating-point image } conventions; ///< Conventions to honour + // The following options are particular to writing images; they needn't be set for anything else. + psFitsFloat floatType; ///< Desired custom floating-point for output images int bitpix; ///< Desired BITPIX for output images; 0 to use as provided - psFitsFloat floatType; ///< Desired custom floating-point for output images + psFitsScaling scaling; ///< Scaling scheme to use when quantising floating-point values + bool fuzz; ///< Fuzz the values when quantising floating-point values? + double bscale, bzero; ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL) + double mean, stdev; ///< Mean and standard deviation of image + int stdevBits; ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*) + float stdevNum; ///< Number of standard deviations to pad off the edge +} psFitsOptions; + + +/// FITS file +typedef struct { + fitsfile* fd; ///< the CFITSIO fits files handle. + bool writable; ///< Is the file writable? + psFitsOptions *options; ///< Options for FITS I/O, or NULL } psFits; + /** FITS compression settings. */ @@ -140,4 +156,7 @@ ); +/// Allocator for options +psFitsOptions *psFitsOptionsAlloc(void); + /** Enables/configures FITS compression. * Index: /trunk/psLib/src/fits/psFitsFloatFile.c =================================================================== --- /trunk/psLib/src/fits/psFitsFloatFile.c (revision 16184) +++ /trunk/psLib/src/fits/psFitsFloatFile.c (revision 16185) @@ -44,5 +44,7 @@ PS_ASSERT_FITS_NON_NULL(fits, PS_FITS_FLOAT_NONE); - if (!fits->conventions.psBitpix) { + psFitsOptions *options = fits->options; // FITS I/O options + + if (!options || !options->conventions.psBitpix) { return PS_FITS_FLOAT_NONE; } Index: /trunk/psLib/src/fits/psFitsHeader.c =================================================================== --- /trunk/psLib/src/fits/psFitsHeader.c (revision 16184) +++ /trunk/psLib/src/fits/psFitsHeader.c (revision 16185) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-16 01:04:56 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-23 03:08:03 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -112,5 +112,5 @@ PS_ASSERT_FITS_NON_NULL(fits, false); - if (!fits->conventions.compression) { + if (fits->options && !fits->options->conventions.compression) { // User has turned off compression conventions; doesn't want any nasty surprises return false; @@ -220,5 +220,6 @@ bool compressed = false; // Is this a compressed image? - if (fits->conventions.compression && fits_is_compressed_image(fits->fd, &status)) { + if ((!fits->options || fits->options->conventions.compression) && + fits_is_compressed_image(fits->fd, &status)) { compressed = true; } Index: /trunk/psLib/src/fits/psFitsImage.c =================================================================== --- /trunk/psLib/src/fits/psFitsImage.c (revision 16184) +++ /trunk/psLib/src/fits/psFitsImage.c (revision 16185) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-16 20:10:35 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-23 03:08:03 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -14,5 +14,5 @@ #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif @@ -36,4 +36,5 @@ #include "psFitsFloatFile.h" #include "psFitsHeader.h" +#include "psFitsScale.h" #include "psMemory.h" @@ -241,264 +242,4 @@ } -# if (0) -// XXX this needs to be optional (eg, invalid for a mask) - -// Apply the BSCALE and BZERO for an image with a "fuzz", so that we get the image as it should be written to -// disk. -// The idea is that the "fuzz" (adding a random number between 0 and 1) preserves the expectation value of -// the image (e.g., a value of 0.1 will get translated to zero 90% of the time, and unity 10% of the time), -// though at the cost of adding an additional variance of 1/12 (a standard deviation of ~0.29). -static psImage *scaleImageForDisk(psImage *image, // Image to which to apply BSCALE and BZERO - int bitpix, // Output BITPIX - double bscale, // Scaling - double bzero, // Zero point - psRandom *rng // Random number generator (for the "fuzz"), or NULL - ) -{ - assert(image); - - if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) { - return psMemIncrRefCounter(image); - } - - psElemType outType; // Type for output image - // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS - // file - switch (bitpix) { - case 8: - outType = PS_TYPE_S8; - break; - case 16: - outType = PS_TYPE_S16; - break; - case 32: - outType = PS_TYPE_S32; - break; - case 64: - outType = PS_TYPE_S64; - break; - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix); - return NULL; - } - - if (bscale == 1.0 && bzero == 0.0) { - return psImageCopy(NULL, image, outType); - } - - int numCols = image->numCols, numRows = image->numRows; // Size of image - psImage *out = psImageAlloc(numCols, numRows, outType); // Output image - - if (!psMemIncrRefCounter(rng)) { - // Don't blab about which seed we're going to get --- it's not necessary for this purpose - psU64 seed = p_psRandomGetSystemSeed(false); - rng = psRandomAlloc(PS_RANDOM_TAUS, seed); - } - - -#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \ - case PS_TYPE_##OUTTYPE: { \ - ps##INTYPE scale = 1.0 / bscale; \ - ps##INTYPE zero = bzero; \ - for (int y = 0; y < numRows; y++) { \ - for (int x = 0; x < numCols; x++) { \ - /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \ - ps##INTYPE random = psRandomUniform(rng); \ - (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \ - } \ - } \ - break; \ - } - -#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \ - case PS_TYPE_##INTYPE: { \ - switch (outType) { \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \ - default: \ - psAbort("Should be unreachable."); \ - } \ - break; \ - } - - switch (image->type.type) { - SCALE_WRITE_IN_CASE(image, F32, out); - SCALE_WRITE_IN_CASE(image, F64, out); - default: - psAbort("Should be unreachable."); - } - - psFree(rng); - - return out; -} -# endif - -# if (0) -// XXX supporting code needs to make this an optional operation -// Determine BSCALE and BZERO for an image, and generate a new image with it applied -// TRUE = BZERO + BSCALE * FITS -static psImage *scaleImageDetermine(double *bscale, // Scaling, to return - double *bzero, // Zero point, to return - psImage *image, // Image to scale - int bitpix, // Desired bits per pixel - psRandom *rng // Random number generator for scaleImageForDisk - ) -{ - PS_ASSERT_PTR_NON_NULL(bscale, NULL); - PS_ASSERT_PTR_NON_NULL(bzero, NULL); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL); - - *bscale = 0.0; - *bzero = 0.0; - - switch (bitpix) { - case 0: - // No scaling applied - return psMemIncrRefCounter(image); - case 8: - case 16: - case 32: - case 64: - // Nothing to do; just allowing these values to pass through - break; - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix); - return NULL; - } - - int numCols = image->numCols, numRows = image->numRows; - double range = pow(2.0, bitpix); // Range of values for target BITPIX - -#define SCALE_DETERMINE_CASE(IN, INTYPE) \ - case PS_TYPE_##INTYPE: { \ - ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \ - for (int y = 0; y < numRows; y++) { \ - for (int x = 0; x < numCols; x++) { \ - ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \ - if (isfinite(value)) { \ - if (value < min) { \ - min = value; \ - } \ - if (value > max) { \ - max = value; \ - } \ - } \ - } \ - } \ - if (!isfinite(min) || !isfinite(max)) { \ - psWarning("No valid values in image to derive BSCALE,BZERO --- using original image."); \ - *bscale = 1.0; \ - *bzero = 0.0; \ - return psMemIncrRefCounter(image); \ - } \ - if (min == max) { \ - *bscale = 1.0; \ - *bzero = min; \ - } else { \ - *bscale = (max - min) / (range - 1.0); \ - *bzero = min + 0.5 * range * (*bscale); \ - } \ - break; \ - } - - switch (image->type.type) { - SCALE_DETERMINE_CASE(image, F32); - SCALE_DETERMINE_CASE(image, F64); - default: - psAbort("Should be unreachable."); - } - psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero); - - return scaleImageForDisk(image, bitpix, *bscale, *bzero, rng); -} -# endif - -#if 0 -// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at -// the present time, since cfitsio should apply the scaling itself in the process of reading. However, we may -// later desire it. -static psImage *scaleImageFromDisk(psFits *fits, psImage *image) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - - if (bscale == 0.0) { - // BSCALE = 0 means don't apply anything - return psMemIncrRefCounter(image); - } - - psElemType inType = image->type.type; // Type for input image - psElemType outType; // Type for output image - switch (inType) { - case PS_TYPE_S8: - case PS_TYPE_S16: - case PS_TYPE_S32: - case PS_TYPE_U8: - case PS_TYPE_U16: - outType = PS_TYPE_F32; - break; - case PS_TYPE_S64: - case PS_TYPE_U32: - case PS_TYPE_U64: - outType = PS_TYPE_F64; - break; - // Including floating-point types just in case someone wants to apply a BSCALE and BZERO to them. - case PS_TYPE_F32: - outType = PS_TYPE_F32; - break; - case PS_TYPE_F64: - outType = PS_TYPE_F64; - break; - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unsupported image type: %x", inType); - return NULL; - } - - int numCols = image->numCols, numRows = image->numRows; - psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image - - -#define SCALE_READ_OUT_CASE(INTYPE, OUTTYPE) \ - case PS_TYPE_##OUTTYPE: { \ - for (int y = 0; y < numRows; y++) { \ - for (int x = 0; x < numCols; x++) { \ - out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \ - } \ - } \ - break; \ - } - -#define SCALE_READ_IN_CASE(INTYPE) \ - case PS_TYPE_##INTYPE: { \ - switch (outType) { \ - SCALE_READ_OUT_CASE(INTYPE, F32); \ - SCALE_READ_OUT_CASE(INTYPE, F64); \ - default: \ - psAbort("Should never get here: type %x should be F32 or F64", outType); \ - } \ - break; \ - } - - switch (inType) { - SCALE_READ_IN_CASE(S8); - SCALE_READ_IN_CASE(S16); - SCALE_READ_IN_CASE(S32); - SCALE_READ_IN_CASE(S64); - SCALE_READ_IN_CASE(U8); - SCALE_READ_IN_CASE(U16); - SCALE_READ_IN_CASE(U32); - SCALE_READ_IN_CASE(U64); - SCALE_READ_IN_CASE(F32); - SCALE_READ_IN_CASE(F64); - default: - psAbort("Should never get here: type %x should be integer", inType); - } - - return out; -} -#endif // Convert an image to the desired BITPIX, i.e., the desired disk representation @@ -522,43 +263,51 @@ *floatType = PS_FITS_FLOAT_NONE; + psFitsOptions *options = fits->options; + if (!options) { + return psMemIncrRefCounter((psImage*)image); // Casting away const + } + // Custom floating-point - if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix && - fits->floatType != PS_FITS_FLOAT_NONE) { - *floatType = fits->floatType; - return psFitsFloatImageToDisk(image, fits->floatType); + if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->conventions.psBitpix && + options->floatType != PS_FITS_FLOAT_NONE) { + *floatType = options->floatType; + return psFitsFloatImageToDisk(image, options->floatType); } // Automatically select what we're given - if (fits->bitpix == 0) { + if (options->bitpix == 0) { return psMemIncrRefCounter((psImage*)image); // Casting away const } // Quantise floating-point images - // XXX this needs to be more controlled: certainly not valid for output masks! - # if (0) - if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) { + if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->bitpix > 0) { if (newScaleZero) { - return scaleImageDetermine(bscale, bzero, (psImage*)image, fits->bitpix, rng); - } - // Get the current BSCALE and BZERO - int status = 0; // Status of cfitsio - if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) { - psFitsError(status, true, "Unable to read header."); - return NULL; - } - status = 0; - if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) { - psFitsError(status, true, "Unable to read header."); - return NULL; - } - status = 0; - if (*bscale == 0.0) { - psError(PS_ERR_IO, true, - "Supposed to use old values of BSCALE and BZERO, but they don't exist."); - return NULL; - } - return scaleImageForDisk((psImage*)image, fits->bitpix, *bscale, *bzero, rng); - } - # endif + // Choose an appropriate BSCALE and BZERO + if (!psFitsScaleDetermine(bscale, bzero, image, fits)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine BSCALE and BZERO for image."); + return NULL; + } + } else { + // Don't want to muck with the current BSCALE and BZERO. Get the current values and use those. + int status = 0; // Status of cfitsio + if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + return NULL; + } + status = 0; + if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + return NULL; + } + status = 0; + if (*bscale == 0.0) { + psError(PS_ERR_IO, true, + "Supposed to use old values of BSCALE and BZERO, but they don't exist."); + return NULL; + } + } + + return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng); + } // Choose the appropriate output type, given the input type and desired bits per pixel @@ -576,5 +325,5 @@ psElemType inType = image->type.type; // Type for input image psElemType outType; // Type for output image - switch (fits->bitpix) { + switch (options->bitpix) { CONVERT_TYPE_INT_CASE(outType, inType, 8); CONVERT_TYPE_INT_CASE(outType, inType, 16); @@ -585,5 +334,5 @@ default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", - fits->bitpix); + options->bitpix); return NULL; } @@ -592,4 +341,9 @@ return psMemIncrRefCounter((psImage*)image); } + + if (PSELEMTYPE_SIZEOF(inType) > PSELEMTYPE_SIZEOF(outType)) { + psWarning("Truncating image pixels to write to disk."); + } + return psImageCopy(NULL, image, outType); } @@ -765,5 +519,6 @@ bzero = cfitsioBzero; } - assert(bitPix == fits->bitpix || fits->bitpix == 0); + psFitsOptions *options = fits->options; // FITS I/O options + assert(!options || bitPix == options->bitpix || options->bitpix == 0); int naxis = 3; // Number of axes @@ -807,5 +562,5 @@ // an unsigned integer type). In all other cases, we have already converted the image to use the // appropriate scale and zero (because we want to apply a randomiser to the quantisation). - fits_set_bscale(fits->fd, 1.0, bzero, &status); + fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status); if (bscale != 0.0) { @@ -895,5 +650,6 @@ bzero = cfitsioBzero; } - assert(bitPix == fits->bitpix || fits->bitpix == 0); + psFitsOptions *options = fits->options; // FITS I/O options + assert(!options || bitPix == options->bitpix || options->bitpix == 0); //check to see if the HDU has the same datatype Index: /trunk/psLib/src/fits/psFitsScale.c =================================================================== --- /trunk/psLib/src/fits/psFitsScale.c (revision 16185) +++ /trunk/psLib/src/fits/psFitsScale.c (revision 16185) @@ -0,0 +1,464 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "psAbort.h" +#include "psAssert.h" +#include "psError.h" +#include "psTrace.h" + +#include "psFits.h" +#include "psImageBackground.h" +#include "psStats.h" +#include "psImageStructManip.h" + +#include "psMemory.h" + +#include "psFitsScale.h" + +// Remember: +// TRUE(i.e., value in memory) = BZERO + BSCALE * FITS(i.e., value on disk) + + +#define MEAN_STAT PS_STAT_ROBUST_MEDIAN // Statistic to use for mean +#define STDEV_STAT PS_STAT_ROBUST_STDEV // Statistic to use for stdev + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Determine appropriate BSCALE and BZERO for an image, preserving dynamic range +static bool scaleRange(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + const psImage *image, // Image to scale + const psFitsOptions *options // FITS options + ) +{ + assert(bscale); + assert(bzero); + assert(image); + assert(options); + + double range = pow(2.0, options->bitpix); // Range of values for target BITPIX + double min = INFINITY, max = -INFINITY; // Minimum and maximum values + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Determine the minimum and maximum values +#define SCALE_DETERMINE_CASE(IN, INTYPE) \ + case PS_TYPE_##INTYPE: { \ + 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; \ + } \ + } \ + } \ + } \ + break; \ + } + + switch (image->type.type) { + SCALE_DETERMINE_CASE(image, F32); + SCALE_DETERMINE_CASE(image, F64); + default: + psAbort("Should be unreachable."); + } + + if (!isfinite(min) || !isfinite(max)) { + psWarning("No valid values in image to derive BSCALE,BZERO"); + *bscale = 1.0; + *bzero = 0.0; + return false; + } + if (min == max) { + *bscale = 1.0; + *bzero = min; + } else { + *bscale = (max - min) / (range - 1.0); + *bzero = min + 0.5 * range * (*bscale); + } + + return true; +} + +// Determine appropriate BSCALE and BZERO for an image, mapping the standard deviation to the nominated number +// of bits +static bool scaleStdev(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + const psImage *image, // Image to scale + const psFitsOptions *options // FITS options + ) +{ + assert(bscale); + assert(bzero); + assert(image); + assert(options); + + double mean = options->mean, stdev = options->stdev; // Mean and standard deviation + if (!isfinite(mean) || !isfinite(stdev)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mean (%f) or stdev (%f) is non-finite. Did you forget to call psFitsScaleMeasure()?", + mean, stdev); + return false; + } + + long range = 1 << options->stdevBits; // Range of values to carry standard deviation + *bscale = stdev / (double) range; + + double imageVal; // Value on image + long diskVal; // Corresponding value on disk + switch (options->scaling) { + case PS_FITS_SCALE_STDEV_POSITIVE: + // Put (mean - N sigma) at the lowest possible value: predominantly positive images + imageVal = mean - options->stdevNum * stdev; + diskVal = - (1 << (options->bitpix - 1)); + break; + case PS_FITS_SCALE_STDEV_NEGATIVE: + // Put (mean + N sigma) at the highest possible value: predominantly negative images + imageVal = mean + options->stdevNum * stdev; + diskVal = (1 << (options->bitpix - 1)) - 1; + break; + case PS_FITS_SCALE_STDEV_BOTH: + // Put mean right in the middle: images with an equal abundance of positive and negative values + imageVal = mean; + diskVal = 0; + break; + default: + psAbort("Should never get here."); + } + + *bzero = imageVal - *bscale * diskVal; + + return true; + + + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + +// Measure statistics on the image for the PS_FITS_SCALE_STDEV_* scaling options +bool psFitsScaleMeasure(psFits *fits, + const psImage *image, + const psImage *mask, + psMaskType maskVal, + psRandom *rng // Random number generator, or NULL + ) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + // Let psImageBackground worry about the mask + // If we don't have a RNG, we will allocate one + + psFitsOptions *options = fits->options; // FITS options + if (!options) { + // No scaling desired + return true; + } + + if (options->scaling != PS_FITS_SCALE_STDEV_POSITIVE && + options->scaling != PS_FITS_SCALE_STDEV_NEGATIVE && + options->scaling != PS_FITS_SCALE_STDEV_BOTH) { + // No need to do anything + return true; + } + + if (!psMemIncrRefCounter(rng)) { + psU64 seed = p_psRandomGetSystemSeed(false); + rng = psRandomAlloc(PS_RANDOM_TAUS, seed); + } + + psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on image"); + psFree(rng); + psFree(stats); + return false; + } + psFree(rng); + + options->mean = psStatsGetValue(stats, MEAN_STAT); // Mean + options->stdev = psStatsGetValue(stats, STDEV_STAT); // Standard deviation + psFree(stats); + + return true; +} + + +bool psFitsScaleDetermine(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + const psImage *image, // Image to scale + const psFits *fits // FITS options + ) +{ + PS_ASSERT_PTR_NON_NULL(bscale, false); + PS_ASSERT_PTR_NON_NULL(bzero, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_FITS_NON_NULL(fits, false); + + *bscale = 1.0; + *bzero = 0.0; + + psFitsOptions *options = fits->options; // FITS options + if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || !options) { + return true; + } + + switch (options->bitpix) { + case 0: + // No scaling applied + return true; + 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 0,8,16,32,64", + options->bitpix); + return false; + } + + switch (options->scaling) { + case PS_FITS_SCALE_NONE: + // No scaling applied + break; + case PS_FITS_SCALE_RANGE: + if (!scaleRange(bscale, bzero, image, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_STDEV_POSITIVE: + case PS_FITS_SCALE_STDEV_NEGATIVE: + case PS_FITS_SCALE_STDEV_BOTH: + if (!scaleStdev(bscale, bzero, image, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + case PS_FITS_SCALE_MANUAL: + *bscale = options->bscale; + *bzero = options->bzero; + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised FITS scaling method: %x", + options->scaling); + return false; + } + + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero); + return true; +} + + +psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, + psRandom *rng) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_FITS_NON_NULL(fits, NULL); + + psFitsOptions *options = fits->options; // FITS options + if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || !options || options->bitpix == 0) { + // No scaling desired + return psMemIncrRefCounter((psImage*)image); // Casting away "const" + } + + int bitpix = options->bitpix; // Bits per pixel + 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 && !options->fuzz) { + 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) && options->fuzz) { + // 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; \ + ps##INTYPE min = - (1 << (options->bitpix - 1)); \ + ps##INTYPE max = (1 << (options->bitpix - 1)) - 1; \ + for (int y = 0; y < numRows; y++) { \ + for (int x = 0; x < numCols; x++) { \ + ps##INTYPE value = ((IN)->data.INTYPE[y][x] - zero) * scale; \ + if (options->fuzz) { \ + /* Add random factor [0,1): adds a variance of 1/12, */ \ + /* but preserves the expectation value */ \ + value += psRandomUniform(rng); \ + } \ + /* Check for underflow and overflow */ \ + (OUT)->data.OUTTYPE[y][x] = (value < min ? min : (value > max ? max : value)); \ + } \ + } \ + 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; +} + + +#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 (e.g., if we ever make our own FITS implementation). +psImage *psFitsScaleFromDisk(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 + + + +psFitsScaling psFitsScalingFromString(const char *string) +{ + PS_ASSERT_STRING_NON_EMPTY(string, PS_FITS_SCALE_NONE); + + if (strcasecmp(string, "RANGE") == 0) return PS_FITS_SCALE_RANGE; + if (strcasecmp(string, "STDEV_POSITIVE") == 0) return PS_FITS_SCALE_STDEV_POSITIVE; + if (strcasecmp(string, "STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_STDEV_NEGATIVE; + if (strcasecmp(string, "STDEV_BOTH") == 0) return PS_FITS_SCALE_STDEV_BOTH; + if (strcasecmp(string, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL; + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to interpret FITS scaling: %s", string); + return PS_FITS_SCALE_NONE; +} + + Index: /trunk/psLib/src/fits/psFitsScale.h =================================================================== --- /trunk/psLib/src/fits/psFitsScale.h (revision 16185) +++ /trunk/psLib/src/fits/psFitsScale.h (revision 16185) @@ -0,0 +1,47 @@ +#ifndef PS_FITS_SCALE_H +#define PS_FITS_SCALE_H + +#include +#include +#include +#include + +/// Measure statistics on the image for the PS_FITS_SCALE_STDEV_ options +/// +/// NB: It is the user's responsibility to call this function if there is the possibility of scaling being +/// performed. The reason this is required is that this function takes a mask (required for measuring the +/// statistics properly), whereas we would like to avoid having to pass a mask to all the other FITS +/// functions. Therefore, since the other FITS functions don't have the mask, while the user may, it is the +/// user's responsibility to call this function. +bool psFitsScaleMeasure(psFits *fits, ///< FITS file + const psImage *image, ///< Image to measure + const psImage *mask, ///< Mask image, or NULL + psMaskType maskVal, ///< Value to mask + psRandom *rng ///< Random number generator, or NULL + ); + +/// Determine BSCALE and BZERO for an image +bool psFitsScaleDetermine(double *bscale, ///< Scaling, to return + double *bzero, ///< Zero point, to return + const psImage *image, ///< Image to scale + const psFits *fits ///< FITS options + ); + +/// Apply the BSCALE and BZERO for an image, so that we get the image as it should be written to disk. +/// +/// "Fuzz" may be optionally added. 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). +psImage *psFitsScaleForDisk(const psImage *image, ///< Image to which to apply BSCALE and BZERO + const psFits *fits, ///< FITS file + double bscale, ///< Scaling + double bzero, ///< Zero point + psRandom *rng ///< Random number generator (for the "fuzz"), or NULL + ); + +/// Interpret a string as a scaling method +psFitsScaling psFitsScalingFromString(const char *string ///< String to interpret + ); + +#endif Index: /trunk/psLib/src/pslib_strict.h =================================================================== --- /trunk/psLib/src/pslib_strict.h (revision 16184) +++ /trunk/psLib/src/pslib_strict.h (revision 16185) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-12-15 01:20:03 $ +* @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ +* @date $Date: 2008-01-23 03:08:03 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -42,4 +42,5 @@ #include "psFitsFloat.h" #include "psFitsFloatFile.h" +#include "psFitsScale.h" //#include "psXML.h"