Changeset 16122 for branches/pap_branch_080117/psLib/src/fits/psFitsImage.c

Timestamp:

Jan 17, 2008, 3:04:21 PM (19 years ago)

Author:

Paul Price

Message:

Updating FITS scaling capability (quantising floating-point images). See bug 1032. Moving scaling functions into their own file, and providing more flexibility.

File:

• branches/pap_branch_080117/psLib/src/fits/psFitsImage.c (modified) (6 diffs)

branches/pap_branch_080117/psLib/src/fits/psFitsImage.c

@@ -7 +7 @@
  *  @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 +14 @@
 
 #ifdef HAVE_CONFIG_H
-# include "config.h"
+#include "config.h"
 #endif
 
@@ -36 +36 @@
 #include "psFitsFloatFile.h"
 #include "psFitsHeader.h"
+#include "psFitsScale.h"
 
 #include "psMemory.h"
@@ -241 +242 @@
 }
 
-# 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 +276 @@
 
     // 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 +336 @@
         return psMemIncrRefCounter((psImage*)image);
     }
+
     return psImageCopy(NULL, image, outType);
 }