Changeset 17660 for trunk/psLib

Timestamp:

May 13, 2008, 3:27:25 PM (18 years ago)

Author:

Paul Price

Message:

Some fixes to image I/O: (1) when BITPIX=8, the FITS standard says unsigned ints (cf two's complement for BITPIX=16,32,64). Added some special cases for quantising to BITPIX=8 to accomodate this. (2) Adding support for marking bad pixels when doing compression. Out-of-range and non-finite pixels are written with value 2^{(BITPIX-1)-1, which is recovered (to NAN) when reading (through use of (Z)BLANK keyword). Added test for this. It doesn't pass yet --- the highest pixel in the image is being interpreted as NAN when it shouldn't be. Solution probably involves tweaking scaleRange() in psFitsScale.c, but not worrying about that now.}

Location:

trunk/psLib

Files:

• src/fits/psFitsImage.c (modified) (7 diffs)
• src/fits/psFitsScale.c (modified) (12 diffs)
• src/fits/psFitsScale.h (modified) (1 diff)
• test/fits (modified) (1 prop)
• test/fits/.cvsignore (modified) (1 diff)
• test/fits/Makefile.am (modified) (1 diff)
• test/fits/tap_psFitsImage.c (added)

trunk/psLib/src/fits/psFitsImage.c

@@ -7 +7 @@
  *  @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 +246 @@
 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 +285 @@
         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 +366 @@
     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 +523 @@
 
     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 +602 @@
     }
 
+    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 +668 @@
 
     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) {

trunk/psLib/src/fits/psFitsScale.c

@@ -46 +46 @@
     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 +86 @@
         *bzero = min;
     } else {
-        *bscale = (max - min) / (range - 1.0);
+        *bscale = (max - min) / range ;
         *bzero = min + 0.5 * range * (*bscale);
     }
@@ -157 +157 @@
 
     return true;
-
-
-    return true;
 }
 
@@ -168 +165 @@
 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 +172 @@
     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 +178 @@
     *bscale = 1.0;
     *bzero = 0.0;
+    *blank = 0;
 
     psFitsOptions *options = fits->options; // FITS options
@@ -201 +201 @@
     }
 
+    *blank = (1 << (options->bitpix - 1)) - 1;
+
     switch (options->scaling) {
       case PS_FITS_SCALE_NONE:
@@ -229 +231 @@
     }
 
-    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 +261 @@
     switch (bitpix) {
       case 8:
-        outType = PS_TYPE_S8;
+        // Note: Use unsigned integer for BITPIX=8
+        outType = PS_TYPE_U8;
         break;
       case 16:
@@ -283 +293 @@
 #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 +312 @@
                         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 +323 @@
     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); \

trunk/psLib/src/fits/psFitsScale.h

@@ -10 +10 @@
 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

trunk/psLib/test/fits

@@ -13 +13 @@
 gmon.out
 tap_psFits
+tap_psFitsImage

@@ -13 +13 @@
 gmon.out
 tap_psFits
+tap_psFitsImage

trunk/psLib/test/fits/.cvsignore

@@ -13 +13 @@
 gmon.out
 tap_psFits
+tap_psFitsImage

trunk/psLib/test/fits/Makefile.am

@@ -11 +11 @@
 
 TEST_PROGS = \
-        tap_psFits
-        tap_psFitsBlank_00
+        tap_psFits \
+        tap_psFitsBlank_00 \
+        tap_psFitsImage
 
 if BUILD_TESTS

@@ -13 +13 @@
 gmon.out
 tap_psFits
+tap_psFitsImage