IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16185


Ignore:
Timestamp:
Jan 22, 2008, 5:08:03 PM (18 years ago)
Author:
Paul Price
Message:

Merging pap_branch_080117 onto mainline.

Location:
trunk/psLib/src
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fits/Makefile.am

    r15630 r16185  
    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

  • trunk/psLib/src/fits/psFits.c

    r15919 r16185  
    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.77 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-01-23 03:08:03 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    7878        fitsClose(fits);
    7979    }
    80     psFree (fits->extword);
     80    psFree(fits->options);
    8181}
    8282
    … …  
    160160
    161161    psFits* fits = psAlloc(sizeof(psFits));
     162    psMemSetDeallocator(fits, (psFreeFunc)fitsFree);
     163
    162164    fits->fd = fptr;
    163165    fits->writable = (iomode == READWRITE);
    164     fits->extword = NULL;
    165     fits->conventions.compression = true;
    166     fits->conventions.psBitpix = true;
    167     fits->bitpix = 0;
    168     fits->floatType = PS_FITS_FLOAT_NONE;
    169     psMemSetDeallocator(fits,(psFreeFunc)fitsFree);
     166
     167    fits->options = NULL;
    170168
    171169    return fits;
     170}
     171
     172
     173static void fitsOptionsFree(psFitsOptions *options)
     174{
     175    psFree(options->extword);
     176}
     177
     178
     179psFitsOptions *psFitsOptionsAlloc(void)
     180{
     181    psFitsOptions *options = psAlloc(sizeof(psFitsOptions)); // Options, to return
     182    psMemSetDeallocator(options, (psFreeFunc)fitsOptionsFree);
     183
     184    options->extword = NULL;
     185
     186    options->conventions.compression = true;
     187    options->conventions.psBitpix = true;
     188
     189    options->floatType = PS_FITS_FLOAT_NONE;
     190
     191    options->bitpix = 0;
     192
     193    options->scaling = PS_FITS_SCALE_NONE;
     194    options->fuzz = true;
     195    options->bscale = 1.0;
     196    options->bzero = 0.0;
     197    options->mean = NAN;
     198    options->stdev = NAN;
     199    options->stdevBits = 4;
     200    options->stdevNum = 5.0;
     201
     202    return options;
    172203}
    173204
    … …  
    227258}
    228259
    229 bool psFitsSetExtnameWord (psFits *fits, const char *extword)
     260bool psFitsSetExtnameWord(psFits *fits, const char *extword)
    230261{
    231262    PS_ASSERT_PTR_NON_NULL(fits,    false);
    232263    PS_ASSERT_PTR_NON_NULL(extword, false);
    233264
    234     psFree (fits->extword);
    235     fits->extword = psStringCopy (extword);
     265    if (!fits->options) {
     266        fits->options = psFitsOptionsAlloc();
     267    }
     268
     269    psFree(fits->options->extword);
     270    fits->options->extword = psStringCopy(extword);
    236271    return true;
    237272}
    … …  
    251286    int status = 0;
    252287
    253     if (!fits->conventions.compression && !fits->extword) {
     288    psFitsOptions *options = fits->options; // FITS options
     289    if (options && !options->conventions.compression && !options->extword) {
    254290        // User wants to use cfitsio.  Good luck to them!
    255291        if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
    … …  
    260296    }
    261297
    262     bool ignoreCI = (fits->conventions.compression &&
     298    bool ignoreCI = (options->conventions.compression &&
    263299                     (strcmp(extname, "COMPRESSED_IMAGE") != 0)); // Ignore COMPRESSED_IMAGE extension name?
    264     char *extword = (fits->extword ? fits->extword : "EXTNAME"); // Word to use as extension name
     300    char *extword = (options->extword ? options->extword : "EXTNAME"); // Word to use as extension name
    265301
    266302#if 0
    … …  
    365401    char name[MAX_STRING_LENGTH];
    366402
    367     char *extword = (fits->extword == NULL) ? defaultExtword : fits->extword;
     403    psFitsOptions *options = fits->options; // FITS options
     404    char *extword = (!options || !options->extword) ? defaultExtword : options->extword;
    368405
    369406    if (fits_read_key_str(fits->fd, extword, name, NULL, &status) != 0) {
    370         psError(PS_ERR_BAD_PARAMETER_NULL, true,
    371                 _("Header keyword %s is not found"), extword);
     407        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("Header keyword %s is not found"), extword);
    372408        return NULL;
    373409    }
    … …  
    382418    int status = 0;
    383419
    384     char *extword = (fits->extword == NULL) ? defaultExtword : fits->extword;
     420    psFitsOptions *options = fits->options; // FITS options
     421    char *extword = (!options || !options->extword) ? defaultExtword : options->extword;
    385422
    386423    if (fits_update_key_str(fits->fd, extword, (char*)name, NULL, &status) != 0) {
    387424        char fitsErr[MAX_STRING_LENGTH];
    388425        (void)fits_get_errstatus(status, fitsErr);
    389         psError(PS_ERR_IO, true,
    390                 _("Could not write data to file. CFITSIO Error: %s"),
    391                 fitsErr);
     426        psError(PS_ERR_IO, true, _("Could not write data to file. CFITSIO Error: %s"), fitsErr);
    392427        return false;
    393428    }
    … …  
    404439
    405440    // move to the specified HDU
    406     if (! psFitsMoveExtNum(fits,extnum,relative) ) {
    407         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    408                 "Failed to delete HDU #%d",
    409                 extnum);
     441    if (!psFitsMoveExtNum(fits, extnum, relative) ) {
     442        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Failed to delete HDU #%d", extnum);
    410443        return false;
    411444    }
    … …  
    417450        char fitsErr[MAX_STRING_LENGTH];
    418451        (void)fits_get_errstatus(status, fitsErr);
    419         psError(PS_ERR_IO, true,
    420                 _("Could not write data to file. CFITSIO Error: %s"),
    421                 fitsErr);
     452        psError(PS_ERR_IO, true, _("Could not write data to file. CFITSIO Error: %s"), fitsErr);
    422453        return false;
    423454    }
    … …  
    598629    // that is the restriction; data must be 32 or 64 bit for noise bits to be valid.
    599630    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         }
     631        if (fits_set_noise_bits(fits->fd, noisebits, &status)) {
     632            fits_set_compression_type(fits->fd, 0x0, &status);
     633            char fitsErr[MAX_STRING_LENGTH];
     634            fits_get_errstatus(status, fitsErr);
     635            psError(PS_ERR_BAD_FITS, true,
     636                    "Error while configuring compression. CFITSIO error: %s", fitsErr);
     637            return false;
     638        }
    608639    }
    609640

  • trunk/psLib/src/fits/psFits.h

    r15630 r16185  
    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.36 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2008-01-23 03:08:03 $
    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;
    4744
    48 /** FITS file object.
    49  *
    50  *  This object should be considered opaque to the user; no item in this
    51  *  struct should be accessed directly.
    52  *
    53  */
     45/// FITS scaling method: 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_POSITIVE,       ///< Auto-scale to sample stdev, place mean at lower limit
     50    PS_FITS_SCALE_STDEV_NEGATIVE,       ///< Auto-scale to sample stdev, place mean at upper limit
     51    PS_FITS_SCALE_STDEV_BOTH,           ///< Auto-scale to sample stdev, place mean at middle
     52    PS_FITS_SCALE_MANUAL                ///< Manual scaling (use specified BSCALE and BZERO)
     53} psFitsScaling;
     54
     55/// Options for FITS I/O
    5456typedef struct {
    55     fitsfile* fd;                       ///< the CFITSIO fits files handle.
    56     bool writable;                      ///< Is the file writable?
    5757    char *extword;                      ///< user-specified word to name extensions (NULL implies EXTNAME)
    5858    struct {
    … …  
    6060        bool psBitpix;                  ///< Custom floating-point image
    6161    } conventions;                      ///< Conventions to honour
     62    // The following options are particular to writing images; they needn't be set for anything else.
     63    psFitsFloat floatType;              ///< Desired custom floating-point for output images
    6264    int bitpix;                         ///< Desired BITPIX for output images; 0 to use as provided
    63     psFitsFloat floatType;              ///< Desired custom floating-point for output images
     65    psFitsScaling scaling;              ///< Scaling scheme to use when quantising floating-point values
     66    bool fuzz;                          ///< Fuzz the values when quantising floating-point values?
     67    double bscale, bzero;               ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL)
     68    double mean, stdev;                 ///< Mean and standard deviation of image
     69    int stdevBits;                      ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*)
     70    float stdevNum;                     ///< Number of standard deviations to pad off the edge
     71} psFitsOptions;
     72
     73
     74/// FITS file
     75typedef struct {
     76    fitsfile* fd;                       ///< the CFITSIO fits files handle.
     77    bool writable;                      ///< Is the file writable?
     78    psFitsOptions *options;             ///< Options for FITS I/O, or NULL
    6479} psFits;
     80
    6581
    6682/** FITS compression settings. */
    … …  
    140156);
    141157
     158/// Allocator for options
     159psFitsOptions *psFitsOptionsAlloc(void);
     160
    142161/** Enables/configures FITS compression.
    143162 *

  • trunk/psLib/src/fits/psFitsFloatFile.c

    r15630 r16185  
    4444    PS_ASSERT_FITS_NON_NULL(fits, PS_FITS_FLOAT_NONE);
    4545
    46     if (!fits->conventions.psBitpix) {
     46    psFitsOptions *options = fits->options; // FITS I/O options
     47
     48    if (!options || !options->conventions.psBitpix) {
    4749        return PS_FITS_FLOAT_NONE;
    4850    }

  • trunk/psLib/src/fits/psFitsHeader.c

    r15630 r16185  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-11-16 01:04:56 $
     9 *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-01-23 03:08:03 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    112112    PS_ASSERT_FITS_NON_NULL(fits, false);
    113113
    114     if (!fits->conventions.compression) {
     114    if (fits->options && !fits->options->conventions.compression) {
    115115        // User has turned off compression conventions; doesn't want any nasty surprises
    116116        return false;
    … …  
    220220
    221221    bool compressed = false;            // Is this a compressed image?
    222     if (fits->conventions.compression && fits_is_compressed_image(fits->fd, &status)) {
     222    if ((!fits->options || fits->options->conventions.compression) &&
     223        fits_is_compressed_image(fits->fd, &status)) {
    223224        compressed = true;
    224225    }

  • trunk/psLib/src/fits/psFitsImage.c

    r16095 r16185  
    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.24 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-01-23 03:08:03 $
    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
    … …  
    522263    *floatType = PS_FITS_FLOAT_NONE;
    523264
     265    psFitsOptions *options = fits->options;
     266    if (!options) {
     267        return psMemIncrRefCounter((psImage*)image); // Casting away const
     268    }
     269
    524270    // Custom floating-point
    525     if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix &&
    526         fits->floatType != PS_FITS_FLOAT_NONE) {
    527         *floatType = fits->floatType;
    528         return psFitsFloatImageToDisk(image, fits->floatType);
     271    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->conventions.psBitpix &&
     272        options->floatType != PS_FITS_FLOAT_NONE) {
     273        *floatType = options->floatType;
     274        return psFitsFloatImageToDisk(image, options->floatType);
    529275    }
    530276
    531277    // Automatically select what we're given
    532     if (fits->bitpix == 0) {
     278    if (options->bitpix == 0) {
    533279        return psMemIncrRefCounter((psImage*)image); // Casting away const
    534280    }
    535281
    536282    // Quantise floating-point images
    537     // XXX this needs to be more controlled: certainly not valid for output masks!
    538     # if (0)
    539     if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) {
     283    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && options->bitpix > 0) {
    540284        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
     285            // Choose an appropriate BSCALE and BZERO
     286            if (!psFitsScaleDetermine(bscale, bzero, image, fits)) {
     287                psError(PS_ERR_UNKNOWN, false, "Unable to determine BSCALE and BZERO for image.");
     288                return NULL;
     289            }
     290        } else {
     291            // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
     292            int status = 0;                 // Status of cfitsio
     293            if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
     294                psFitsError(status, true, "Unable to read header.");
     295                return NULL;
     296            }
     297            status = 0;
     298            if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
     299                psFitsError(status, true, "Unable to read header.");
     300                return NULL;
     301            }
     302            status = 0;
     303            if (*bscale == 0.0) {
     304                psError(PS_ERR_IO, true,
     305                        "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
     306                return NULL;
     307            }
     308        }
     309
     310        return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng);
     311    }
    563312
    564313    // Choose the appropriate output type, given the input type and desired bits per pixel
    … …  
    576325    psElemType inType = image->type.type; // Type for input image
    577326    psElemType outType;                 // Type for output image
    578     switch (fits->bitpix) {
     327    switch (options->bitpix) {
    579328        CONVERT_TYPE_INT_CASE(outType, inType, 8);
    580329        CONVERT_TYPE_INT_CASE(outType, inType, 16);
    … …  
    585334      default:
    586335        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64",
    587                 fits->bitpix);
     336                options->bitpix);
    588337        return NULL;
    589338    }
    … …  
    592341        return psMemIncrRefCounter((psImage*)image);
    593342    }
     343
     344    if (PSELEMTYPE_SIZEOF(inType) > PSELEMTYPE_SIZEOF(outType)) {
     345        psWarning("Truncating image pixels to write to disk.");
     346    }
     347
    594348    return psImageCopy(NULL, image, outType);
    595349}
    … …  
    765519        bzero = cfitsioBzero;
    766520    }
    767     assert(bitPix == fits->bitpix || fits->bitpix == 0);
     521    psFitsOptions *options = fits->options; // FITS I/O options
     522    assert(!options || bitPix == options->bitpix || options->bitpix == 0);
    768523
    769524    int naxis = 3;                      // Number of axes
    … …  
    807562    // an unsigned integer type).  In all other cases, we have already converted the image to use the
    808563    // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
    809     fits_set_bscale(fits->fd, 1.0, bzero, &status);
     564    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
    810565
    811566    if (bscale != 0.0) {
    … …  
    895650        bzero = cfitsioBzero;
    896651    }
    897     assert(bitPix == fits->bitpix || fits->bitpix == 0);
     652    psFitsOptions *options = fits->options; // FITS I/O options
     653    assert(!options || bitPix == options->bitpix || options->bitpix == 0);
    898654
    899655    //check to see if the HDU has the same datatype

  • trunk/psLib/src/pslib_strict.h

    r15841 r16185  
    99*  @author Eric Van Alst, MHPCC
    1010*
    11 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2007-12-15 01:20:03 $
     11*  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2008-01-23 03:08:03 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4242#include "psFitsFloat.h"
    4343#include "psFitsFloatFile.h"
     44#include "psFitsScale.h"
    4445
    4546//#include "psXML.h"
Note: See TracChangeset for help on using the changeset viewer.