Index: trunk/psLib/src/fits/psFitsImage.c =================================================================== --- trunk/psLib/src/fits/psFitsImage.c (revision 17656) +++ trunk/psLib/src/fits/psFitsImage.c (revision 17660) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-04-17 23:43:02 $ + * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-05-14 01:27:25 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -246,4 +246,5 @@ static psImage *imageToDiskRepresentation(double *bscale, // Scaling applied double *bzero, // Zero point applied + long *blank, // Blank value (integer data) psFitsFloat *floatType, // Type of custom floating-point psFits *fits, // FITS file pointer @@ -284,5 +285,5 @@ if (newScaleZero) { // Choose an appropriate BSCALE and BZERO - if (!psFitsScaleDetermine(bscale, bzero, image, fits)) { + if (!psFitsScaleDetermine(bscale, bzero, blank, image, fits)) { // We can't have the write dying for this reason --- try to save it somehow! psWarning("Unable to determine BSCALE and BZERO for image --- setting to 1.0, 0.0"); @@ -365,8 +366,32 @@ psAssert(!output->parent, "impossible"); // No parents means the buffer is contiguous + void *nullValue = NULL; // Null value for data + float nullFloat = NAN; // Null value for floating point + double nullDouble = NAN; // Null value for double + switch (info->psDatatype) { + case PS_TYPE_F32: + nullValue = &nullFloat; + break; + case PS_TYPE_F64: + nullValue = &nullDouble; + break; + case PS_TYPE_U8: + case PS_TYPE_U16: + case PS_TYPE_U32: + case PS_TYPE_U64: + case PS_TYPE_S8: + case PS_TYPE_S16: + case PS_TYPE_S32: + case PS_TYPE_S64: + // Can't mark bad pixels any further than what is in the FITS image + break; + default: + psAbort("Unknown type: %x", info->psDatatype); + } + int anynull = 0; // Are there any NULLs in the data? int status = 0; // cfitsio status if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel, - info->increment, NULL, output->data.V[0], &anynull, &status) != 0) { + info->increment, nullValue, output->data.V[0], &anynull, &status) != 0) { psFitsError(status, true, "Reading FITS file failed."); return false; @@ -498,6 +523,7 @@ double bscale = 0.0, bzero = 0.0; // Scale and zero point to put in header (*already* applied to data) + long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &floatType, fits, image, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, image, NULL, true); // Image to write out if (!diskImage) { @@ -576,4 +602,17 @@ } + if (blank != 0) { + // Some quantisation has taken place --- record the blank ("magic") pixel value + char *keyword = (psFitsCompressionGetType(fits) != PS_FITS_COMPRESS_NONE && + (!fits->options || fits->options->conventions.compression)) ? + "ZBLANK" : "BLANK"; // Keyword for recording blank pixel value + fits_write_key_lng(fits->fd, keyword, blank, "Value for undefined pixels", &status); + fits_set_imgnull(fits->fd, blank, &status); + if (psFitsError(status, true, "Could not write BLANK header to file.")) { + psFree(diskImage); + return false; + } + } + if (floatType != PS_FITS_FLOAT_NONE) { psFitsFloatImageSet(fits, floatType); @@ -629,6 +668,7 @@ double bscale = 0.0, bzero = 0.0; // Scale and zero point to put in header (*already* applied to data) + long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &floatType, fits, input, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, input, NULL, false); // Image to write out if (!diskImage) { Index: trunk/psLib/src/fits/psFitsScale.c =================================================================== --- trunk/psLib/src/fits/psFitsScale.c (revision 17656) +++ trunk/psLib/src/fits/psFitsScale.c (revision 17660) @@ -46,5 +46,5 @@ psAssert(options, "impossible"); - double range = pow(2.0, options->bitpix); // Range of values for target BITPIX + double range = pow(2.0, options->bitpix) - 1.0; // Range of values for target BITPIX, reduced by the BLANK double min = INFINITY, max = -INFINITY; // Minimum and maximum values int numCols = image->numCols, numRows = image->numRows; // Size of image @@ -86,5 +86,5 @@ *bzero = min; } else { - *bscale = (max - min) / (range - 1.0); + *bscale = (max - min) / range ; *bzero = min + 0.5 * range * (*bscale); } @@ -157,7 +157,4 @@ return true; - - - return true; } @@ -168,4 +165,5 @@ bool psFitsScaleDetermine(double *bscale, // Scaling, to return double *bzero, // Zero point, to return + long *blank, // Blank value, to return const psImage *image, // Image to scale const psFits *fits // FITS options @@ -174,4 +172,5 @@ PS_ASSERT_PTR_NON_NULL(bscale, false); PS_ASSERT_PTR_NON_NULL(bzero, false); + PS_ASSERT_PTR_NON_NULL(blank, false); PS_ASSERT_IMAGE_NON_NULL(image, false); PS_ASSERT_FITS_NON_NULL(fits, false); @@ -179,4 +178,5 @@ *bscale = 1.0; *bzero = 0.0; + *blank = 0; psFitsOptions *options = fits->options; // FITS options @@ -201,4 +201,6 @@ } + *blank = (1 << (options->bitpix - 1)) - 1; + switch (options->scaling) { case PS_FITS_SCALE_NONE: @@ -229,5 +231,12 @@ } - psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero); + if (options->bitpix == 8) { + // FITS standard wants unsigned for BITPIX=8, two's-complement for BITPIX=16,32,64. + *bzero -= *bscale * 128; + *blank = 255; + } + + + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BLANK = %ld\n", *bscale, *bzero, *blank); return true; } @@ -252,5 +261,6 @@ switch (bitpix) { case 8: - outType = PS_TYPE_S8; + // Note: Use unsigned integer for BITPIX=8 + outType = PS_TYPE_U8; break; case 16: @@ -283,8 +293,9 @@ #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; \ + double scale = 1.0 / bscale; \ + double zero = bzero; \ + /* Note: BITPIX=8 treated differently, since it uses unsigned values; the rest use signed */ \ + double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1); \ + double max = bitpix == 8 ? 255 : (pow(2.0, options->bitpix - 1) - 1.0); \ for (int y = 0; y < numRows; y++) { \ for (int x = 0; x < numCols; x++) { \ @@ -301,6 +312,6 @@ value += psRandomUniform(rng); \ } \ - /* Check for underflow and overflow */ \ - (OUT)->data.OUTTYPE[y][x] = (value < min ? min : (value > max ? max : value)); \ + /* Check for underflow and overflow; set either to max */ \ + (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : value); \ } \ } \ @@ -312,5 +323,5 @@ case PS_TYPE_##INTYPE: { \ switch (outType) { \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8); \ SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \ SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \ Index: trunk/psLib/src/fits/psFitsScale.h =================================================================== --- trunk/psLib/src/fits/psFitsScale.h (revision 17656) +++ trunk/psLib/src/fits/psFitsScale.h (revision 17660) @@ -10,4 +10,5 @@ bool psFitsScaleDetermine(double *bscale, ///< Scaling, to return double *bzero, ///< Zero point, to return + long *blank, ///< Blank value, to return const psImage *image, ///< Image to scale const psFits *fits ///< FITS options