Index: /branches/pap_branch_080117/psLib/src/fits/Makefile.am =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/Makefile.am (revision 16121) +++ /branches/pap_branch_080117/psLib/src/fits/Makefile.am (revision 16122) @@ -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: /branches/pap_branch_080117/psLib/src/fits/psFits.c =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/psFits.c (revision 16121) +++ /branches/pap_branch_080117/psLib/src/fits/psFits.c (revision 16122) @@ -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.76.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-18 01:04:21 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -165,6 +165,16 @@ fits->conventions.compression = true; fits->conventions.psBitpix = true; + + fits->floatType = PS_FITS_FLOAT_NONE; fits->bitpix = 0; - fits->floatType = PS_FITS_FLOAT_NONE; + fits->scaling = PS_FITS_SCALE_NONE; + fits->fuzz = true; + fits->bscale = 1.0; + fits->bzero = 0.0; + fits->mean = NAN; + fits->stdev = NAN; + fits->stdevBits = 8; + fits->stdevNum = 5.0; + psMemSetDeallocator(fits,(psFreeFunc)fitsFree); @@ -598,12 +608,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: /branches/pap_branch_080117/psLib/src/fits/psFits.h =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/psFits.h (revision 16121) +++ /branches/pap_branch_080117/psLib/src/fits/psFits.h (revision 16122) @@ -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.35.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-18 01:04:21 $ * 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,11 +34,22 @@ } 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 scaling: 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_LOWER, ///< Auto-scale to sample stdev, place mean at lower limit + PS_FITS_SCALE_STDEV_UPPER, ///< Auto-scale to sample stdev, place mean at upper limit + PS_FITS_SCALE_STDEV_MIDDLE, ///< Auto-scale to sample stdev, place mean at middle + PS_FITS_SCALE_MANUAL ///< Manual scaling (use specified BSCALE and BZERO) +} psFitsScaling; /** FITS file object. @@ -60,6 +67,13 @@ 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 (SCALE_MANUAL) + double mean, stdev; ///< Mean and standard deviation of image + int stdevBits; ///< Number of bits to sample a standard deviation (SCALE_STDEV) + float stdevNum; ///< Number of standard deviations to pad off the edge } psFits; Index: /branches/pap_branch_080117/psLib/src/fits/psFitsImage.c =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/psFitsImage.c (revision 16121) +++ /branches/pap_branch_080117/psLib/src/fits/psFitsImage.c (revision 16122) @@ -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.23.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-01-18 01:04:21 $ * * 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 @@ -535,30 +276,33 @@ // 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 (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->bitpix, *bscale, *bzero, fits->fuzz, rng); + } // Choose the appropriate output type, given the input type and desired bits per pixel @@ -592,4 +336,5 @@ return psMemIncrRefCounter((psImage*)image); } + return psImageCopy(NULL, image, outType); } Index: /branches/pap_branch_080117/psLib/src/fits/psFitsScale.c =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/psFitsScale.c (revision 16122) +++ /branches/pap_branch_080117/psLib/src/fits/psFitsScale.c (revision 16122) @@ -0,0 +1,436 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#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(MEMORY) = BZERO + BSCALE * FITS(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 psFits *fits // FITS options + ) +{ + assert(bscale); + assert(bzero); + assert(image); + assert(fits); + + double range = pow(2.0, fits->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 psFits *fits // FITS options + ) +{ + assert(bscale); + assert(bzero); + assert(image); + assert(fits); + + double mean = fits->mean, stdev = fits->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 << fits->stdevBits; // Range of values to carry standard deviation + *bscale = fits->stdev / (double) range; + + int bitpix = fits->bitpix; // Bits per pixel + double imageVal; // Value on image + long diskVal; // Corresponding value on disk + switch (fits->scaling) { + case PS_FITS_SCALE_STDEV_LOWER: + // Put (mean - N sigma) at the lowest possible value: predominantly positive images + imageVal = mean - fits->stdevNum * stdev; + diskVal = - (1 << (bitpix - 1)); + break; + case PS_FITS_SCALE_STDEV_UPPER: + // Put (mean + N sigma) at the highest possible value: predominantly negative images + imageVal = mean + fits->stdevNum * stdev; + diskVal = (1 << (bitpix - 1)) - 1; + break; + case PS_FITS_SCALE_STDEV_MIDDLE: + // 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 + + if (fits->scaling != PS_FITS_SCALE_STDEV_LOWER && + fits->scaling != PS_FITS_SCALE_STDEV_UPPER && + fits->scaling != PS_FITS_SCALE_STDEV_MIDDLE) { + // 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); + + fits->mean = psStatsGetValue(stats, MEAN_STAT); // Mean + fits->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; + + if (!PS_IS_PSELEMTYPE_REAL(image->type.type)) { + return true; + } + + switch (fits->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", + fits->bitpix); + return false; + } + + switch (fits->scaling) { + case PS_FITS_SCALE_NONE: + // No scaling applied + break; + case PS_FITS_SCALE_RANGE: + if (!scaleRange(bscale, bzero, image, fits)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_STDEV_LOWER: + case PS_FITS_SCALE_STDEV_UPPER: + case PS_FITS_SCALE_STDEV_MIDDLE: + if (!scaleStdev(bscale, bzero, image, fits)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + case PS_FITS_SCALE_MANUAL: + *bscale = fits->bscale; + *bzero = fits->bzero; + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised FITS scaling option: %x", fits->scaling); + return false; + } + + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero); + return true; +} + + +psImage *psFitsScaleForDisk(const psImage *image, int bitpix, double bscale, double bzero, bool fuzz, + psRandom *rng) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + + if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) { + return psMemIncrRefCounter((psImage*)image); // Casting away "const" + } + + 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) && 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; \ + if (fuzz) { \ + /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \ + for (int y = 0; y < numRows; y++) { \ + for (int x = 0; x < numCols; x++) { \ + ps##INTYPE random = psRandomUniform(rng); \ + (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \ + } \ + } \ + } else { \ + /* No random factor added: data preserved, but expectation value is not */ \ + for (int y = 0; y < numRows; y++) { \ + for (int x = 0; x < numCols; x++) { \ + (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale; \ + } \ + } \ + } \ + 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 Index: /branches/pap_branch_080117/psLib/src/fits/psFitsScale.h =================================================================== --- /branches/pap_branch_080117/psLib/src/fits/psFitsScale.h (revision 16122) +++ /branches/pap_branch_080117/psLib/src/fits/psFitsScale.h (revision 16122) @@ -0,0 +1,45 @@ +#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 options + 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 + int bitpix, // Output BITPIX + double bscale, // Scaling + double bzero, // Zero point + bool fuzz, // Apply fuzz? + psRandom *rng // Random number generator (for the "fuzz"), or NULL + ); + + +#endif