IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16122


Ignore:
Timestamp:
Jan 17, 2008, 3:04:21 PM (18 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.

Location:
branches/pap_branch_080117/psLib/src/fits
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_080117/psLib/src/fits/Makefile.am

    r15630 r16122  
    1010        psFitsTable.c \
    1111        psFitsFloat.c \
    12         psFitsFloatFile.c
     12        psFitsFloatFile.c \
     13        psFitsScale.c
    1314
    1415EXTRA_DIST = fits.i
    … …  
    2021        psFitsTable.h \
    2122        psFitsFloat.h \
    22         psFitsFloatFile.h
     23        psFitsFloatFile.h \
     24        psFitsScale.h
    2325
    2426CLEANFILES = *~ *.bb *.bbg *.da

  • branches/pap_branch_080117/psLib/src/fits/psFits.c

    r15919 r16122  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.76 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-12-25 01:29:42 $
     9 *  @version $Revision: 1.76.2.1 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-01-18 01:04:21 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    165165    fits->conventions.compression = true;
    166166    fits->conventions.psBitpix = true;
     167
     168    fits->floatType = PS_FITS_FLOAT_NONE;
    167169    fits->bitpix = 0;
    168     fits->floatType = PS_FITS_FLOAT_NONE;
     170    fits->scaling = PS_FITS_SCALE_NONE;
     171    fits->fuzz = true;
     172    fits->bscale = 1.0;
     173    fits->bzero = 0.0;
     174    fits->mean = NAN;
     175    fits->stdev = NAN;
     176    fits->stdevBits = 8;
     177    fits->stdevNum = 5.0;
     178
    169179    psMemSetDeallocator(fits,(psFreeFunc)fitsFree);
    170180
    … …  
    598608    // that is the restriction; data must be 32 or 64 bit for noise bits to be valid.
    599609    if (type != PS_FITS_COMPRESS_PLIO) {
    600         if (fits_set_noise_bits(fits->fd, noisebits, &status)) {
    601             fits_set_compression_type(fits->fd, 0x0, &status);
    602             char fitsErr[MAX_STRING_LENGTH];
    603             fits_get_errstatus(status, fitsErr);
    604             psError(PS_ERR_BAD_FITS, true,
    605                     "Error while configuring compression. CFITSIO error: %s", fitsErr);
    606             return false;
    607         }
     610        if (fits_set_noise_bits(fits->fd, noisebits, &status)) {
     611            fits_set_compression_type(fits->fd, 0x0, &status);
     612            char fitsErr[MAX_STRING_LENGTH];
     613            fits_get_errstatus(status, fitsErr);
     614            psError(PS_ERR_BAD_FITS, true,
     615                    "Error while configuring compression. CFITSIO error: %s", fitsErr);
     616            return false;
     617        }
    608618    }
    609619

  • branches/pap_branch_080117/psLib/src/fits/psFits.h

    r15630 r16122  
    44 * @author Robert DeSonia, MHPCC
    55 *
    6  * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-11-16 01:04:56 $
     6 * @version $Revision: 1.35.2.1 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2008-01-18 01:04:21 $
    88 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    99 */
    … …  
    2525#include "psErrorCodes.h"
    2626
    27 /** FITS HDU type.
    28  *
    29  *  Enumeration for FITS HDU type.
    30  *
    31  */
     27/// FITS HDU type.
    3228typedef enum {
    3329    PS_FITS_TYPE_NONE = -1,            ///< Unknown HDU type
    … …  
    3834} psFitsType;
    3935
     36/// FITS compression type
    4037typedef enum {
    41     PS_FITS_COMPRESS_NONE = 0,
    42     PS_FITS_COMPRESS_GZIP,
    43     PS_FITS_COMPRESS_RICE,
    44     PS_FITS_COMPRESS_HCOMPRESS,
    45     PS_FITS_COMPRESS_PLIO
     38    PS_FITS_COMPRESS_NONE = 0,          ///< No compression
     39    PS_FITS_COMPRESS_GZIP,              ///< GZIP compression (of the pixels only)
     40    PS_FITS_COMPRESS_RICE,              ///< RICE compression (of the pixels only)
     41    PS_FITS_COMPRESS_HCOMPRESS,         ///< HCOMPRESS compression (of the pixels only)
     42    PS_FITS_COMPRESS_PLIO               ///< PLIO compression (of the pixels only; appropriate for masks)
    4643} psFitsCompressionType;
     44
     45/// FITS scaling: how to set BSCALE and BZERO
     46typedef enum {
     47    PS_FITS_SCALE_NONE,                 ///< No auto-scaling to be applied (BSCALE = 1, BZERO = 0)
     48    PS_FITS_SCALE_RANGE,                ///< Auto-scale to preserve dynamic range
     49    PS_FITS_SCALE_STDEV_LOWER,          ///< Auto-scale to sample stdev, place mean at lower limit
     50    PS_FITS_SCALE_STDEV_UPPER,          ///< Auto-scale to sample stdev, place mean at upper limit
     51    PS_FITS_SCALE_STDEV_MIDDLE,         ///< Auto-scale to sample stdev, place mean at middle
     52    PS_FITS_SCALE_MANUAL                ///< Manual scaling (use specified BSCALE and BZERO)
     53} psFitsScaling;
    4754
    4855/** FITS file object.
    … …  
    6067        bool psBitpix;                  ///< Custom floating-point image
    6168    } conventions;                      ///< Conventions to honour
     69    // The following options are particular to writing images; they needn't be set for anything else.
     70    psFitsFloat floatType;              ///< Desired custom floating-point for output images
    6271    int bitpix;                         ///< Desired BITPIX for output images; 0 to use as provided
    63     psFitsFloat floatType;              ///< Desired custom floating-point for output images
     72    psFitsScaling scaling;              ///< Scaling scheme to use when quantising floating-point values
     73    bool fuzz;                          ///< Fuzz the values when quantising floating-point values?
     74    double bscale, bzero;               ///< Manually specified BSCALE and BZERO (SCALE_MANUAL)
     75    double mean, stdev;                 ///< Mean and standard deviation of image
     76    int stdevBits;                      ///< Number of bits to sample a standard deviation (SCALE_STDEV)
     77    float stdevNum;                     ///< Number of standard deviations to pad off the edge
    6478} psFits;
    6579

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

    r16095 r16122  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2008-01-16 20:10:35 $
     9 *  @version $Revision: 1.23.2.1 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-01-18 01:04:21 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1414
    1515#ifdef HAVE_CONFIG_H
    16 # include "config.h"
     16#include "config.h"
    1717#endif
    1818
    … …  
    3636#include "psFitsFloatFile.h"
    3737#include "psFitsHeader.h"
     38#include "psFitsScale.h"
    3839
    3940#include "psMemory.h"
    … …  
    241242}
    242243
    243 # if (0)
    244 // XXX this needs to be optional (eg, invalid for a mask)
    245 
    246 // Apply the BSCALE and BZERO for an image with a "fuzz", so that we get the image as it should be written to
    247 // disk.
    248 // The idea is that the "fuzz" (adding a random number between 0 and 1) preserves the expectation value of
    249 // the image (e.g., a value of 0.1 will get translated to zero 90% of the time, and unity 10% of the time),
    250 // though at the cost of adding an additional variance of 1/12 (a standard deviation of ~0.29).
    251 static psImage *scaleImageForDisk(psImage *image, // Image to which to apply BSCALE and BZERO
    252                                   int bitpix, // Output BITPIX
    253                                   double bscale, // Scaling
    254                                   double bzero, // Zero point
    255                                   psRandom *rng // Random number generator (for the "fuzz"), or NULL
    256                                   )
    257 {
    258     assert(image);
    259 
    260     if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) {
    261         return psMemIncrRefCounter(image);
    262     }
    263 
    264     psElemType outType;                 // Type for output image
    265     // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS
    266     // file
    267     switch (bitpix) {
    268       case 8:
    269         outType = PS_TYPE_S8;
    270         break;
    271       case 16:
    272         outType = PS_TYPE_S16;
    273         break;
    274       case 32:
    275         outType = PS_TYPE_S32;
    276         break;
    277       case 64:
    278         outType = PS_TYPE_S64;
    279         break;
    280       default:
    281         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
    282         return NULL;
    283     }
    284 
    285     if (bscale == 1.0 && bzero == 0.0) {
    286         return psImageCopy(NULL, image, outType);
    287     }
    288 
    289     int numCols = image->numCols, numRows = image->numRows; // Size of image
    290     psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
    291 
    292     if (!psMemIncrRefCounter(rng)) {
    293         // Don't blab about which seed we're going to get --- it's not necessary for this purpose
    294         psU64 seed = p_psRandomGetSystemSeed(false);
    295         rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
    296     }
    297 
    298 
    299 #define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
    300     case PS_TYPE_##OUTTYPE: { \
    301         ps##INTYPE scale = 1.0 / bscale; \
    302         ps##INTYPE zero = bzero; \
    303         for (int y = 0; y < numRows; y++) { \
    304             for (int x = 0; x < numCols; x++) { \
    305                 /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \
    306                 ps##INTYPE random = psRandomUniform(rng); \
    307                 (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
    308             } \
    309         } \
    310         break; \
    311     }
    312 
    313 #define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
    314     case PS_TYPE_##INTYPE: { \
    315         switch (outType) { \
    316             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \
    317             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \
    318             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \
    319             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \
    320           default: \
    321             psAbort("Should be unreachable."); \
    322         } \
    323         break; \
    324     }
    325 
    326     switch (image->type.type) {
    327         SCALE_WRITE_IN_CASE(image, F32, out);
    328         SCALE_WRITE_IN_CASE(image, F64, out);
    329       default:
    330         psAbort("Should be unreachable.");
    331     }
    332 
    333     psFree(rng);
    334 
    335     return out;
    336 }
    337 # endif
    338 
    339 # if (0)
    340 // XXX supporting code needs to make this an optional operation
    341 // Determine BSCALE and BZERO for an image, and generate a new image with it applied
    342 // TRUE = BZERO + BSCALE * FITS
    343 static psImage *scaleImageDetermine(double *bscale, // Scaling, to return
    344                                     double *bzero, // Zero point, to return
    345                                     psImage *image, // Image to scale
    346                                     int bitpix, // Desired bits per pixel
    347                                     psRandom *rng // Random number generator for scaleImageForDisk
    348                                     )
    349 {
    350     PS_ASSERT_PTR_NON_NULL(bscale, NULL);
    351     PS_ASSERT_PTR_NON_NULL(bzero, NULL);
    352     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    353     PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
    354 
    355     *bscale = 0.0;
    356     *bzero = 0.0;
    357 
    358     switch (bitpix) {
    359       case 0:
    360         // No scaling applied
    361         return psMemIncrRefCounter(image);
    362       case 8:
    363       case 16:
    364       case 32:
    365       case 64:
    366         // Nothing to do; just allowing these values to pass through
    367         break;
    368       default:
    369         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
    370         return NULL;
    371     }
    372 
    373     int numCols = image->numCols, numRows = image->numRows;
    374     double range = pow(2.0, bitpix); // Range of values for target BITPIX
    375 
    376 #define SCALE_DETERMINE_CASE(IN, INTYPE) \
    377     case PS_TYPE_##INTYPE: { \
    378         ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \
    379         for (int y = 0; y < numRows; y++) { \
    380             for (int x = 0; x < numCols; x++) { \
    381                 ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \
    382                 if (isfinite(value)) { \
    383                     if (value < min) { \
    384                         min = value; \
    385                     } \
    386                     if (value > max) { \
    387                         max = value; \
    388                     } \
    389                 } \
    390             } \
    391         } \
    392         if (!isfinite(min) || !isfinite(max)) { \
    393             psWarning("No valid values in image to derive BSCALE,BZERO --- using original image."); \
    394             *bscale = 1.0; \
    395             *bzero = 0.0; \
    396             return psMemIncrRefCounter(image); \
    397         } \
    398         if (min == max) { \
    399             *bscale = 1.0; \
    400             *bzero = min; \
    401         } else { \
    402             *bscale = (max - min) / (range - 1.0); \
    403             *bzero = min + 0.5 * range * (*bscale); \
    404         } \
    405         break; \
    406     }
    407 
    408     switch (image->type.type) {
    409         SCALE_DETERMINE_CASE(image, F32);
    410         SCALE_DETERMINE_CASE(image, F64);
    411       default:
    412         psAbort("Should be unreachable.");
    413     }
    414     psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero);
    415 
    416     return scaleImageForDisk(image, bitpix, *bscale, *bzero, rng);
    417 }
    418 # endif
    419 
    420 #if 0
    421 // This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
    422 // the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
    423 // later desire it.
    424 static psImage *scaleImageFromDisk(psFits *fits, psImage *image)
    425 {
    426     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    427 
    428     if (bscale == 0.0) {
    429         // BSCALE = 0 means don't apply anything
    430         return psMemIncrRefCounter(image);
    431     }
    432 
    433     psElemType inType = image->type.type; // Type for input image
    434     psElemType outType;                 // Type for output image
    435     switch (inType) {
    436       case PS_TYPE_S8:
    437       case PS_TYPE_S16:
    438       case PS_TYPE_S32:
    439       case PS_TYPE_U8:
    440       case PS_TYPE_U16:
    441         outType = PS_TYPE_F32;
    442         break;
    443       case PS_TYPE_S64:
    444       case PS_TYPE_U32:
    445       case PS_TYPE_U64:
    446         outType = PS_TYPE_F64;
    447         break;
    448         // Including floating-point types just in case someone wants to apply a BSCALE and BZERO to them.
    449       case PS_TYPE_F32:
    450         outType = PS_TYPE_F32;
    451         break;
    452       case PS_TYPE_F64:
    453         outType = PS_TYPE_F64;
    454         break;
    455       default:
    456         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unsupported image type: %x", inType);
    457         return NULL;
    458     }
    459 
    460     int numCols = image->numCols, numRows = image->numRows;
    461     psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
    462 
    463 
    464 #define SCALE_READ_OUT_CASE(INTYPE, OUTTYPE) \
    465   case PS_TYPE_##OUTTYPE: { \
    466       for (int y = 0; y < numRows; y++) { \
    467           for (int x = 0; x < numCols; x++) { \
    468               out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
    469           } \
    470       } \
    471       break; \
    472   }
    473 
    474 #define SCALE_READ_IN_CASE(INTYPE) \
    475   case PS_TYPE_##INTYPE: { \
    476       switch (outType) { \
    477           SCALE_READ_OUT_CASE(INTYPE, F32); \
    478           SCALE_READ_OUT_CASE(INTYPE, F64); \
    479         default: \
    480           psAbort("Should never get here: type %x should be F32 or F64", outType); \
    481       } \
    482       break; \
    483   }
    484 
    485     switch (inType) {
    486         SCALE_READ_IN_CASE(S8);
    487         SCALE_READ_IN_CASE(S16);
    488         SCALE_READ_IN_CASE(S32);
    489         SCALE_READ_IN_CASE(S64);
    490         SCALE_READ_IN_CASE(U8);
    491         SCALE_READ_IN_CASE(U16);
    492         SCALE_READ_IN_CASE(U32);
    493         SCALE_READ_IN_CASE(U64);
    494         SCALE_READ_IN_CASE(F32);
    495         SCALE_READ_IN_CASE(F64);
    496       default:
    497           psAbort("Should never get here: type %x should be integer", inType);
    498     }
    499 
    500     return out;
    501 }
    502 #endif
    503244
    504245// Convert an image to the desired BITPIX, i.e., the desired disk representation
    … …  
    535276
    536277    // Quantise floating-point images
    537     // XXX this needs to be more controlled: certainly not valid for output masks!
    538     # if (0)
    539278    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) {
    540279        if (newScaleZero) {
    541             return scaleImageDetermine(bscale, bzero, (psImage*)image, fits->bitpix, rng);
    542         }
    543         // Get the current BSCALE and BZERO
    544         int status = 0;                 // Status of cfitsio
    545         if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
    546             psFitsError(status, true, "Unable to read header.");
    547             return NULL;
    548         }
    549         status = 0;
    550         if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
    551             psFitsError(status, true, "Unable to read header.");
    552             return NULL;
    553         }
    554         status = 0;
    555         if (*bscale == 0.0) {
    556             psError(PS_ERR_IO, true,
    557                     "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
    558             return NULL;
    559         }
    560         return scaleImageForDisk((psImage*)image, fits->bitpix, *bscale, *bzero, rng);
    561     }
    562     # endif
     280            // Choose an appropriate BSCALE and BZERO
     281            if (!psFitsScaleDetermine(bscale, bzero, image, fits)) {
     282                psError(PS_ERR_UNKNOWN, false, "Unable to determine BSCALE and BZERO for image.");
     283                return NULL;
     284            }
     285        } else {
     286            // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
     287            int status = 0;                 // Status of cfitsio
     288            if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
     289                psFitsError(status, true, "Unable to read header.");
     290                return NULL;
     291            }
     292            status = 0;
     293            if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
     294                psFitsError(status, true, "Unable to read header.");
     295                return NULL;
     296            }
     297            status = 0;
     298            if (*bscale == 0.0) {
     299                psError(PS_ERR_IO, true,
     300                        "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
     301                return NULL;
     302            }
     303        }
     304
     305        return psFitsScaleForDisk(image, fits->bitpix, *bscale, *bzero, fits->fuzz, rng);
     306    }
    563307
    564308    // Choose the appropriate output type, given the input type and desired bits per pixel
    … …  
    592336        return psMemIncrRefCounter((psImage*)image);
    593337    }
     338
    594339    return psImageCopy(NULL, image, outType);
    595340}
Note: See TracChangeset for help on using the changeset viewer.