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 Index: trunk/psLib/test/fits/.cvsignore =================================================================== --- trunk/psLib/test/fits/.cvsignore (revision 17656) +++ trunk/psLib/test/fits/.cvsignore (revision 17660) @@ -13,2 +13,3 @@ gmon.out tap_psFits +tap_psFitsImage Index: trunk/psLib/test/fits/Makefile.am =================================================================== --- trunk/psLib/test/fits/Makefile.am (revision 17656) +++ trunk/psLib/test/fits/Makefile.am (revision 17660) @@ -11,6 +11,7 @@ TEST_PROGS = \ - tap_psFits - tap_psFitsBlank_00 + tap_psFits \ + tap_psFitsBlank_00 \ + tap_psFitsImage if BUILD_TESTS Index: trunk/psLib/test/fits/tap_psFitsImage.c =================================================================== --- trunk/psLib/test/fits/tap_psFitsImage.c (revision 17660) +++ trunk/psLib/test/fits/tap_psFitsImage.c (revision 17660) @@ -0,0 +1,95 @@ +#include +#include "tap.h" +#include "pstap.h" +#include + +#define NUMCOLS 12 // Number of columns for image +#define NUMROWS 34 // Number of rows for image +#define BITPIX 32 // Bits per pixel +#define OUTTYPE PS_TYPE_F64 // Expected output type + +// Generate an image +static psImage *generateImage(void) +{ + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 12345); // Random number generator + psImage *image = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); // Generated image + for (int y = 0; y < NUMROWS; y++) { + for (int x = 0; x < NUMCOLS; x++) { + image->data.F32[y][x] = ((x + y) % 2) ? NAN : psRandomUniform(rng); + } + } + psFree(rng); + return image; +} + + +int main(int argc, const char *argv[]) +{ + plan_tests(7); + + { + psMemId id = psMemGetId(); + + psString name = NULL; // Name of FITS image + psStringAppend(&name, "%s.fits", argv[0]); + + psImage *original = generateImage(); // Image to test + psMetadata *header = psMetadataAlloc(); // FITS header + + psFits *fits = psFitsOpen(name, "w"); + + // Set compression + psVector *tiles = psVectorAlloc(2, PS_TYPE_S32); + tiles->data.S32[0] = 0; + tiles->data.S32[1] = 1; + ok(psFitsSetCompression(fits, PS_FITS_COMPRESS_RICE, tiles, 4, 0, 0), "Set compression"); + psFree(tiles); + + // Set quantisation + fits->options = psFitsOptionsAlloc(); + fits->options->bitpix = BITPIX; + fits->options->scaling = PS_FITS_SCALE_RANGE; + + bool write = psFitsWriteImage(fits, header, original, 0, NULL); + psFitsClose(fits); + + ok(write, "Write image"); + skip_start(!fits || !write, 4, "Unable to write image"); + { + psFits *fits = psFitsOpen(name, "r"); + psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); + psFitsClose(fits); + + ok(image, "Read image"); + skip_start(!image, 3, "Unable to read image"); + { + ok(image->type.type == OUTTYPE, "Image type"); + ok(image->numCols == NUMCOLS && image->numRows == NUMROWS, "Image size"); + bool consistent = true; // Images are consistent? + double tol = 1.0 / (double)(1 << (BITPIX - 1)) + FLT_EPSILON; // Expected tolerance + for (int y = 0; y < NUMROWS; y++) { + for (int x = 0; x < NUMCOLS; x++) { + double imageVal = psImageGet(image, x, y); + if ((!isfinite(original->data.F32[y][x]) && isfinite(imageVal)) || + (isfinite(original->data.F32[y][x]) && !isfinite(imageVal) && + original->data.F32[y][x] < 1.0 - 0.5 * tol) || + (fabs(original->data.F32[y][x] - imageVal) > tol)) { + consistent = false; + diag("Inconsistent at %d,%d: %lf vs %f --> %le", x, y, + imageVal, original->data.F32[y][x], imageVal - original->data.F32[y][x]); + } + } + } + ok(consistent, "Images consistent."); + psFree(image); + } + skip_end(); + } + skip_end(); + psFree(original); + psFree(header); + psFree(name); + + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } +}