Index: trunk/psLib/src/fits/psFitsFloat.c =================================================================== --- trunk/psLib/src/fits/psFitsFloat.c (revision 15335) +++ trunk/psLib/src/fits/psFitsFloat.c (revision 15630) @@ -4,4 +4,6 @@ #include +#include +#include #include "psAbort.h" @@ -9,11 +11,15 @@ #include "psError.h" #include "psImage.h" -#include "psMetadata.h" -#include "psRandom.h" #include "psFits.h" +#include "psTrace.h" +#include "psMemory.h" #include "psFitsFloat.h" -#define BIAS_FLOAT_16_0 10 // Exponent bias +#define BIAS_FLOAT_16_0 10 // Exponent bias for FLOAT_16_0 + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// union float_16_0 { @@ -27,5 +33,4 @@ } f16; }; - @@ -57,113 +62,20 @@ } +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// -psImage *psFitsScaleGenerate(psFits *fits, const psImage *image, int bitpix, psRandom *rng) -{ - PS_ASSERT_FITS_NON_NULL(fits, NULL); - PS_ASSERT_FITS_WRITABLE(fits, NULL); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL); - PS_ASSERT_RANDOM_NON_NULL(rng, NULL); - 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; - } - - int numCols = image->numCols, numRows = image->numRows; - double range = pow(2.0, bitpix); // Range of values for target BITPIX - psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image - -#define SCALE_GENERATE_CASE(IN, INTYPE, OUT, OUTTYPE) \ - case PS_TYPE_##OUTTYPE: { \ - for (int y = 0; y < numRows; y++) { \ - for (int x = 0; x < numCols; x++) { \ - ps##INTYPE random = psRandomUniform(rng) - 0.5; /* Randomise the quantisation */ \ - (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \ - } \ - } \ - break; \ - } - -#define SCALE_DETERMINE_CASE(IN, INTYPE, OUT) \ - 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; \ - } \ - } \ - } \ - } \ - ps##INTYPE scale = range / (max - min); \ - ps##INTYPE zero = min - 0.5 * range; /* Minimum value for twos-complement int is - 2^(bitpix-1) */ \ - fits->bscale = 1.0 / scale; \ - fits->bzero = zero; \ - switch (outType) { \ - SCALE_GENERATE_CASE(IN, INTYPE, OUT, S8); \ - SCALE_GENERATE_CASE(IN, INTYPE, OUT, S16); \ - SCALE_GENERATE_CASE(IN, INTYPE, OUT, S32); \ - SCALE_GENERATE_CASE(IN, INTYPE, OUT, S64); \ - default: \ - psAbort("Should be unreachable."); \ - } \ - break; \ - } - - switch (image->type.type) { - SCALE_DETERMINE_CASE(image, F32, out); - SCALE_DETERMINE_CASE(image, F64, out); - default: - psAbort("Should be unreachable."); - } - - return out; -} - -psImage *psFitsScaleApply(const psImage *image, double bscale, double bzero, psElemType type) +psImage *psFitsFloatImageWrite(const psImage *image, psFitsFloat type) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); - - if (bscale == 0.0) { - return psMemIncrRefCounter(image); - } - - switch ( - - -psImage *psFitsFloatImageWrite(const psImage *image, - psMetadata *header, - psFitsFloat type - ) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_METADATA_NON_NULL(header, NULL); // Need the header, so we can mark it as compressed + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); psImage *output = NULL; // Output image, to return - psString convName = NULL; // Convention name switch (type) { + case PS_FITS_FLOAT_NONE: + // No conversion to be performed + return psMemIncrRefCounter((psImage*)image); // Casting away "const" case PS_FITS_FLOAT_16_0: output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S16); // Output image @@ -173,5 +85,4 @@ } } - convName = psStringCopy("FLOAT_16_0"); break; default: @@ -180,38 +91,65 @@ } - psMetadataAddStr(header, PS_LIST_TAIL, "PSBITPIX", PS_META_REPLACE, "Pan-STARRS bit encoding", convName); - return output; } -psImage *psFitsFloatImageRead(const psFits *fits, - const psImage *image - ) + +psImage *psFitsFloatImageRead(psImage *out, const psImage *in, psFitsFloat type) { - PS_ASSERT_FITS_NON_NULL(fits, NULL); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_NON_NULL(in, NULL); - bool mdok; // Status of MD lookup - const char *convName = psMetadataLookupStr(&mdok, header, "PSBITPIX"); // Convention used - if (!mdok || !convName) { - // It's not a custom floating-point image - return psMemIncrRefCounter(image); + psElemType elem = psFitsFloatImageType(type); // Type for elements + if (elem == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise convention: %x", type); + return NULL; } - psImage *output = NULL; // Output image, to return + int numCols = in->numCols, numRows = in->numRows; // Size of image - if (strcmp(convName, "FLOAT_16_0") == 0 && image->type.type == PS_TYPE_S16) { - output = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); // Output image - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - output->data.F32[y][x] = convertF32fromFloat16_0(image->data.S16[y][x]); + out = psImageRecycle(out, numCols, numRows, elem); + + switch (type) { + case PS_FITS_FLOAT_16_0: + if (in->type.type != PS_TYPE_S16) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Convention claims to be FLOAT_16_0, but image type is not S16."); + return NULL; + } + assert(out->type.type == PS_TYPE_F32); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + out->data.F32[y][x] = convertF32fromFloat16_0(in->data.S16[y][x]); } } - return output; + return out; + case PS_FITS_FLOAT_NONE: + default: + psAbort("Should be unreachable"); + } + return NULL; +} + +psElemType psFitsFloatImageType(psFitsFloat type) +{ + switch (type) { + case PS_FITS_FLOAT_16_0: + return PS_TYPE_F32; + case PS_FITS_FLOAT_NONE: + default: + return 0; // Doesn't correspond to ANY type --- should flag a real error + } +} + + +psFitsFloat psFitsFloatTypeFromString(const char *string) +{ + PS_ASSERT_STRING_NON_EMPTY(string, PS_FITS_FLOAT_NONE); + + if (strcmp(string, "FLOAT_16_0") == 0) { + return PS_FITS_FLOAT_16_0; } - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to recognise PSBITPIX convention name: %s", convName); - return NULL; - - return output; + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Unable to recognise custom floating-point convention name: %s", string); + return PS_FITS_FLOAT_NONE; }