Changeset 15630 for trunk/psLib/src/fits/psFitsFloat.c

Timestamp:

Nov 15, 2007, 3:04:56 PM (18 years ago)

Author:

Paul Price

Message:

Integrating the code to choose the output BITPIX (and quantising floating-point images if necessary) and custom floating-point schemes with the functions to read and write images. More detailed description, from my e-mail: psLib has always read any of the (legal) FITS types, and written a FITS image with the appropriate BITPIX for the image->type.type. Now the output BITPIX can be forced by setting psFits.bitpix to the desired BITPIX. Floating point images being quantised to integer types have BSCALE and BZERO set in such a fashion as to maintain the dynamic range of the image; quantisation after applying the scale and zero includes adding a random number [0,1) to preserve the expectation of the floating-point value. On reading, the BSCALE and BZERO are applied by cfitsio without any assistance by psLib. I've also added Gene's custom 16-bit floats (1S 5E 10M). In anticipation of further flavours down the road and to maintain compatibility, I've labelled this as FLOAT_16_0. An F32 image may be written out with this scheme by setting psFits.floatType = PS_FITS_FLOAT_16_0. The image is converted and written out with BITPIX=16, and the header keyword PSBITPIX added with the custom floating-point scheme name. Such images are automatically converted to F32 on reading (if psFits.conventions.psBitpix is true, which it is by default). One thing to keep in mind is that I'm not sure how NANs and INFs are handled in these processes. These additions are made on the top of the compression feature --- quantised or custom floating-point images may also be compressed (though the custom floating-point images don't gain much from it).

File:

• trunk/psLib/src/fits/psFitsFloat.c (modified) (6 diffs)

trunk/psLib/src/fits/psFitsFloat.c

@@ -4 +4 @@
 
 #include <ieee754.h>
+#include <string.h>
+#include <assert.h>
 
 #include "psAbort.h"
@@ -9 +11 @@
 #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 +33 @@
     } f16;
 };
-
 
 
@@ -57 +62 @@
 }
 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// 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 +85 @@
             }
         }
-        convName = psStringCopy("FLOAT_16_0");
         break;
       default:
@@ -180 +91 @@
     }
 
-    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;
 }