IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15335 r15630  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-10-19 23:52:39 $
     9 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-11-16 01:04:56 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2121#include <string.h>
    2222
    23 #include "psFits.h"
     23#include "psAbort.h"
     24#include "psType.h"
    2425#include "psAssert.h"
    25 #include "psFitsImage.h"
    2626#include "psError.h"
    27 
    28 #include "psImageStructManip.h"
    29 #include "psMemory.h"
    3027#include "psString.h"
    3128#include "psLogMsg.h"
    3229#include "psTrace.h"
    3330#include "psVector.h"
    34 
     31#include "psRandom.h"
     32#include "psImageStructManip.h"
     33
     34#include "psFits.h"
     35#include "psFitsFloat.h"
     36#include "psFitsFloatFile.h"
    3537#include "psFitsHeader.h"
    3638
    37 #define MAX_STRING_LENGTH 256  // maximum length string for FITS routines
     39#include "psMemory.h"
     40
     41#include "psFitsImage.h"
     42
     43#define MAX_STRING_LENGTH 256           // maximum length string for FITS routines
    3844
    3945// Information required to read a FITS file
    … …  
    4955} p_psFitsReadInfo;
    5056
    51 static p_psFitsReadInfo *p_psFitsReadInfoAlloc(const psFits *fits, // The FITS file handle
    52         psRegion region, // Region to read
    53         int z // z-plane to read in cube
    54                                               )
     57// Read the vital statistics of this FITS image, in preparation for reading the image
     58static p_psFitsReadInfo *p_psFitsReadInfoAlloc(
     59    const psFits *fits, // The FITS file handle
     60    psRegion region, // Region to read
     61    int z // z-plane to read in cube
     62    )
    5563{
    5664    PS_ASSERT_FITS_NON_NULL(fits, NULL);
    … …  
    6169
    6270    int status = 0;                     // CFITSIO status
    63     char fitsErr[80];                   // CFITSIO error message string
    6471
    6572    // check to see if we even are positioned on an image HDU
    6673    int hduType;                        // Type of HDU
    6774    if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
    68         (void)fits_get_errstatus(status, fitsErr);
    69         psError(PS_ERR_IO, true,
    70                 _("Could not determine the HDU type. CFITSIO Error: %s"),
    71                 fitsErr);
     75        psFitsError(status, true, "Could not determine the HDU type.");
    7276        goto bad;
    7377    }
    7478    if (hduType != IMAGE_HDU) {
    75         psError(PS_ERR_IO, true,
    76                 _("Current FITS HDU type must be an image."));
     79        psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image."));
    7780        goto bad;
    7881    }
    … …  
    8083    // Get the data type 'bitPix' from the FITS image
    8184    if (fits_get_img_equivtype(fits->fd, &info->bitPix, &status) != 0) {
    82         fits_get_errstatus(status, fitsErr);
    83         psError(PS_ERR_IO, true,
    84                 _("Could not determine image data type. CFITSIO Error: %s"),
    85                 fitsErr);
     85        psFitsError(status, true, "Could not determine image data type.");
    8686        goto bad;
    8787    }
    … …  
    8989    /* Get the dimensions 'nAxis' from the FITS image */
    9090    if (fits_get_img_dim(fits->fd, &info->nAxis, &status) != 0) {
    91         (void)fits_get_errstatus(status, fitsErr);
    92         psError(PS_ERR_IO, true,
    93                 _("Could not determine image dimensions. CFITSIO Error: %s"),
    94                 fitsErr);
     91        psFitsError(status, true, "Could not determine image dimensions.");
    9592        goto bad;
    9693    }
    … …  
    9996    if ((info->nAxis < 2) || (info->nAxis > 3)) {
    10097        psError(PS_ERR_IO, true,
    101                 _("Image number of dimensions, %d, is not valid.  Only two or three dimensions supported for FITS I/O."),
    102                 info->nAxis);
     98                _("Image number of dimensions, %d, is not supported."), info->nAxis);
    10399        goto bad;
    104100    }
    … …  
    106102    /* Get the Image size from the FITS file */
    107103    if (fits_get_img_size(fits->fd, info->nAxis, info->nAxes, &status) != 0) {
    108         (void)fits_get_errstatus(status, fitsErr);
    109         psError(PS_ERR_IO, true,
    110                 _("Could not determine image size. CFITSIO Error: %s"),
    111                 fitsErr);
     104        psFitsError(status, true, "Could not determine image size.");
    112105        goto bad;
    113106    }
    … …  
    133126    info->increment[2] = 1;
    134127
    135     switch (info->bitPix) {
    136     case BYTE_IMG:
    137         info->psDatatype = PS_TYPE_U8;
    138         info->fitsDatatype = TBYTE;
    139         break;
    140     case SBYTE_IMG:
    141         info->psDatatype = PS_TYPE_S8;
    142         info->fitsDatatype = TSBYTE;
    143         break;
    144     case USHORT_IMG:
    145         info->psDatatype = PS_TYPE_U16;
    146         info->fitsDatatype = TUSHORT;
    147         break;
    148     case SHORT_IMG:
    149         info->psDatatype = PS_TYPE_S16;
    150         info->fitsDatatype = TSHORT;
    151         break;
    152     case ULONG_IMG:
    153         info->psDatatype = PS_TYPE_U32;
    154         info->fitsDatatype = TUINT;
    155         break;
    156     case LONG_IMG:
    157         info->psDatatype = PS_TYPE_S32;
    158         info->fitsDatatype = TINT;
    159         break;
    160     case LONGLONG_IMG:
    161         info->psDatatype = PS_TYPE_S64;
    162         info->fitsDatatype = TLONGLONG;
    163         break;
    164     case FLOAT_IMG:
    165         info->psDatatype = PS_TYPE_F32;
    166         info->fitsDatatype = TFLOAT;
    167         break;
    168     case DOUBLE_IMG:
    169         info->psDatatype = PS_TYPE_F64;
    170         info->fitsDatatype = TDOUBLE;
    171         break;
    172     default:
    173         psError(PS_ERR_IO, true,
    174                 _("FITS image type, BITPIX=%d, is not supported."),
    175                 info->bitPix);
     128    // Check scale and zero
     129    double bscale = 0.0, bzero = 0.0;    // Scale and zero point
     130    if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) {
     131        psFitsError(status, true, "Unable to read header.");
    176132        goto bad;
    177133    }
    178 
     134    status = 0;
     135    if (fits_read_key_dbl(fits->fd, "BZERO", &bzero, NULL, &status) && status != KEY_NO_EXIST) {
     136        psFitsError(status, true, "Unable to read header.");
     137        goto bad;
     138    }
     139    status = 0;
     140
     141    if ((bscale != 0.0 && bscale != 1.0) || bzero != (int)bzero) {
     142        // It's a floating-point image that's been quantised
     143        // cfitsio will apply the scale and zero point for us if we choose the correct data type
     144        switch (info->bitPix) {
     145          case BYTE_IMG:
     146          case SBYTE_IMG:
     147          case USHORT_IMG:
     148          case SHORT_IMG:
     149          case ULONG_IMG:
     150          case LONG_IMG:
     151          case FLOAT_IMG:
     152            info->psDatatype = PS_TYPE_F32;
     153            info->fitsDatatype = TFLOAT;
     154            break;
     155          case LONGLONG_IMG:
     156          case DOUBLE_IMG:
     157            info->psDatatype = PS_TYPE_F64;
     158            info->fitsDatatype = TDOUBLE;
     159            break;
     160          default:
     161            psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix);
     162            goto bad;
     163        }
     164    } else {
     165        switch (info->bitPix) {
     166          case BYTE_IMG:
     167            info->psDatatype = PS_TYPE_U8;
     168            info->fitsDatatype = TBYTE;
     169            break;
     170          case SBYTE_IMG:
     171            info->psDatatype = PS_TYPE_S8;
     172            info->fitsDatatype = TSBYTE;
     173            break;
     174          case USHORT_IMG:
     175            info->psDatatype = PS_TYPE_U16;
     176            info->fitsDatatype = TUSHORT;
     177            break;
     178          case SHORT_IMG:
     179            info->psDatatype = PS_TYPE_S16;
     180            info->fitsDatatype = TSHORT;
     181            break;
     182          case ULONG_IMG:
     183            info->psDatatype = PS_TYPE_U32;
     184            info->fitsDatatype = TUINT;
     185            break;
     186          case LONG_IMG:
     187            info->psDatatype = PS_TYPE_S32;
     188            info->fitsDatatype = TINT;
     189            break;
     190          case LONGLONG_IMG:
     191            info->psDatatype = PS_TYPE_S64;
     192            info->fitsDatatype = TLONGLONG;
     193            break;
     194          case FLOAT_IMG:
     195            info->psDatatype = PS_TYPE_F32;
     196            info->fitsDatatype = TFLOAT;
     197            break;
     198          case DOUBLE_IMG:
     199            info->psDatatype = PS_TYPE_F64;
     200            info->fitsDatatype = TDOUBLE;
     201            break;
     202          default:
     203            psError(PS_ERR_IO, true, _("FITS image type, BITPIX=%d, is not supported."), info->bitPix);
     204            goto bad;
     205        }
     206    }
    179207    return info;
    180208
    … …  
    211239
    212240    return true;
     241}
     242
     243// Apply the BSCALE and BZERO for an image with a "fuzz"
     244// The idea is that the "fuzz" (adding a random number between 0 and 1) preserves the expectation value of
     245// the image (e.g., a value of 0.1 will get translated to zero 90% of the time, and unity 10% of the time),
     246// though at the cost of adding an additional variance of 1/12 (a standard deviation of ~0.29).
     247static psImage *scaleImageWrite(psImage *image, // Image to which to apply BSCALE and BZERO
     248                                int bitpix, // Output BITPIX
     249                                double bscale, // Scaling
     250                                double bzero, // Zero point
     251                                psRandom *rng // Random number generator (for the "fuzz"), or NULL
     252    )
     253{
     254    assert(image);
     255
     256    if (!PS_IS_PSELEMTYPE_REAL(image->type.type) || bitpix == 0) {
     257        return psMemIncrRefCounter(image);
     258    }
     259
     260    psElemType outType;                 // Type for output image
     261    // Choosing to use signed types because those don't require BSCALE,BZERO to represent them in the FITS
     262    // file
     263    switch (bitpix) {
     264      case 8:
     265        outType = PS_TYPE_S8;
     266        break;
     267      case 16:
     268        outType = PS_TYPE_S16;
     269        break;
     270      case 32:
     271        outType = PS_TYPE_S32;
     272        break;
     273      case 64:
     274        outType = PS_TYPE_S64;
     275        break;
     276      default:
     277        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
     278        return NULL;
     279    }
     280
     281    if (bscale == 1.0 && bzero == 0.0) {
     282        return psImageCopy(NULL, image, outType);
     283    }
     284
     285    int numCols = image->numCols, numRows = image->numRows; // Size of image
     286    psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
     287
     288    if (!psMemIncrRefCounter(rng)) {
     289        // Don't blab about which seed we're going to get --- it's not necessary for this purpose
     290        psU64 seed = p_psRandomGetSystemSeed(false);
     291        rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
     292    }
     293
     294
     295#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
     296    case PS_TYPE_##OUTTYPE: { \
     297        ps##INTYPE scale = 1.0 / bscale; \
     298        ps##INTYPE zero = bzero; \
     299        for (int y = 0; y < numRows; y++) { \
     300            for (int x = 0; x < numCols; x++) { \
     301                /* Add random factor [0,1): adds a variance of 1/12, but preserves the expectation value */ \
     302                ps##INTYPE random = psRandomUniform(rng); \
     303                (OUT)->data.OUTTYPE[y][x] = ((IN)->data.INTYPE[y][x] - zero) * scale + random; \
     304            } \
     305        } \
     306        break; \
     307    }
     308
     309#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
     310    case PS_TYPE_##INTYPE: { \
     311        switch (outType) { \
     312            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S8); \
     313            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \
     314            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \
     315            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \
     316          default: \
     317            psAbort("Should be unreachable."); \
     318        } \
     319        break; \
     320    }
     321
     322    switch (image->type.type) {
     323        SCALE_WRITE_IN_CASE(image, F32, out);
     324        SCALE_WRITE_IN_CASE(image, F64, out);
     325      default:
     326        psAbort("Should be unreachable.");
     327    }
     328
     329    psFree(rng);
     330
     331    return out;
     332}
     333
     334
     335
     336// Determine BSCALE and BZERO for an image, and generate a new image with it applied
     337// TRUE = BZERO + BSCALE * FITS
     338static psImage *scaleImageDetermineWrite(double *bscale, // Scaling, to return
     339                                         double *bzero, // Zero point, to return
     340                                         psImage *image, // Image to scale
     341                                         int bitpix, // Desired bits per pixel
     342                                         psRandom *rng // Random number generator for scaleImageWrite
     343    )
     344{
     345    PS_ASSERT_PTR_NON_NULL(bscale, NULL);
     346    PS_ASSERT_PTR_NON_NULL(bzero, NULL);
     347    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     348    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(image, NULL);
     349
     350    *bscale = 0.0;
     351    *bzero = 0.0;
     352
     353    switch (bitpix) {
     354      case 0:
     355        // No scaling applied
     356        return psMemIncrRefCounter(image);
     357      case 8:
     358      case 16:
     359      case 32:
     360      case 64:
     361        // Nothing to do; just allowing these values to pass through
     362        break;
     363      default:
     364        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64", bitpix);
     365        return NULL;
     366    }
     367
     368    int numCols = image->numCols, numRows = image->numRows;
     369    double range = pow(2.0, bitpix); // Range of values for target BITPIX
     370
     371#define SCALE_DETERMINE_CASE(IN, INTYPE) \
     372    case PS_TYPE_##INTYPE: { \
     373        ps##INTYPE min = INFINITY, max = -INFINITY; /* Minimum and maximum values */ \
     374        for (int y = 0; y < numRows; y++) { \
     375            for (int x = 0; x < numCols; x++) { \
     376                ps##INTYPE value = (IN)->data.INTYPE[y][x]; /* Value of interest */ \
     377                if (isfinite(value)) { \
     378                    if (value < min) { \
     379                        min = value; \
     380                    } \
     381                    if (value > max) { \
     382                        max = value; \
     383                    } \
     384                } \
     385            } \
     386        } \
     387        if (!isfinite(min) || !isfinite(max)) { \
     388            psWarning("No valid values in image to derive BSCALE,BZERO --- using original image."); \
     389            *bscale = 1.0; \
     390            *bzero = 0.0; \
     391            return psMemIncrRefCounter(image); \
     392        } \
     393        if (min == max) { \
     394            *bscale = 1.0; \
     395            *bzero = min; \
     396        } else { \
     397            *bscale = (max - min) / (range - 1.0); \
     398            *bzero = min + 0.5 * range * (*bscale); \
     399        } \
     400        break; \
     401    }
     402
     403    switch (image->type.type) {
     404        SCALE_DETERMINE_CASE(image, F32);
     405        SCALE_DETERMINE_CASE(image, F64);
     406      default:
     407        psAbort("Should be unreachable.");
     408    }
     409    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf\n", *bscale, *bzero);
     410
     411    return scaleImageWrite(image, bitpix, *bscale, *bzero, rng);
     412}
     413
     414
     415#if 0
     416// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
     417// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
     418// later desire it.
     419static psImage *scaleImageRead(psFits *fits, psImage *image)
     420{
     421    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     422
     423    if (bscale == 0.0) {
     424        // BSCALE = 0 means don't apply anything
     425        return psMemIncrRefCounter(image);
     426    }
     427
     428    psElemType inType = image->type.type; // Type for input image
     429    psElemType outType;                 // Type for output image
     430    switch (inType) {
     431      case PS_TYPE_S8:
     432      case PS_TYPE_S16:
     433      case PS_TYPE_S32:
     434      case PS_TYPE_U8:
     435      case PS_TYPE_U16:
     436        outType = PS_TYPE_F32;
     437        break;
     438      case PS_TYPE_S64:
     439      case PS_TYPE_U32:
     440      case PS_TYPE_U64:
     441        outType = PS_TYPE_F64;
     442        break;
     443        // Including floating-point types just in case someone wants to apply a BSCALE and BZERO to them.
     444      case PS_TYPE_F32:
     445        outType = PS_TYPE_F32;
     446        break;
     447      case PS_TYPE_F64:
     448        outType = PS_TYPE_F64;
     449        break;
     450      default:
     451        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unsupported image type: %x", inType);
     452        return NULL;
     453    }
     454
     455    int numCols = image->numCols, numRows = image->numRows;
     456    psImage *out = psImageAlloc(numCols, numRows, outType); // Output scaled image
     457
     458
     459#define SCALE_READ_OUT_CASE(INTYPE, OUTTYPE) \
     460  case PS_TYPE_##OUTTYPE: { \
     461      for (int y = 0; y < numRows; y++) { \
     462          for (int x = 0; x < numCols; x++) { \
     463              out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
     464          } \
     465      } \
     466      break; \
     467  }
     468
     469#define SCALE_READ_IN_CASE(INTYPE) \
     470  case PS_TYPE_##INTYPE: { \
     471      switch (outType) { \
     472          SCALE_READ_OUT_CASE(INTYPE, F32); \
     473          SCALE_READ_OUT_CASE(INTYPE, F64); \
     474        default: \
     475          psAbort("Should never get here: type %x should be F32 or F64", outType); \
     476      } \
     477      break; \
     478  }
     479
     480    switch (inType) {
     481        SCALE_READ_IN_CASE(S8);
     482        SCALE_READ_IN_CASE(S16);
     483        SCALE_READ_IN_CASE(S32);
     484        SCALE_READ_IN_CASE(S64);
     485        SCALE_READ_IN_CASE(U8);
     486        SCALE_READ_IN_CASE(U16);
     487        SCALE_READ_IN_CASE(U32);
     488        SCALE_READ_IN_CASE(U64);
     489        SCALE_READ_IN_CASE(F32);
     490        SCALE_READ_IN_CASE(F64);
     491      default:
     492          psAbort("Should never get here: type %x should be integer", inType);
     493    }
     494
     495    return out;
     496}
     497#endif
     498
     499// Convert an image to the desired BITPIX
     500static psImage *convertImageWrite(double *bscale, // Scaling applied
     501                                  double *bzero, // Zero point applied
     502                                  psFitsFloat *floatType, // Type of custom floating-point
     503                                  psFits *fits, // FITS file pointer
     504                                  const psImage *image, // Current type
     505                                  psRandom *rng, // Random number generator
     506                                  bool newScaleZero // Determine a new BSCALE and BZERO?
     507    )
     508{
     509    assert(bscale);
     510    assert(bzero);
     511    assert(floatType);
     512    assert(fits);
     513    assert(image);
     514
     515    *bscale = 0.0;
     516    *bzero = 0.0;
     517    *floatType = PS_FITS_FLOAT_NONE;
     518
     519    // Custom floating-point
     520    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->conventions.psBitpix &&
     521        fits->floatType != PS_FITS_FLOAT_NONE) {
     522        *floatType = fits->floatType;
     523        return psFitsFloatImageWrite(image, fits->floatType);
     524    }
     525
     526    // Automatically select what we're given
     527    if (fits->bitpix == 0) {
     528        return psMemIncrRefCounter((psImage*)image); // Casting away const
     529    }
     530
     531    // Quantise floating-point images
     532    if (PS_IS_PSELEMTYPE_REAL(image->type.type) && fits->bitpix > 0) {
     533        if (newScaleZero) {
     534            return scaleImageDetermineWrite(bscale, bzero, (psImage*)image, fits->bitpix, rng);
     535        }
     536        // Get the current BSCALE and BZERO
     537        int status = 0;                 // Status of cfitsio
     538        if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
     539            psFitsError(status, true, "Unable to read header.");
     540            return NULL;
     541        }
     542        status = 0;
     543        if (fits_read_key_dbl(fits->fd, "BZERO", bzero, NULL, &status) && status != KEY_NO_EXIST) {
     544            psFitsError(status, true, "Unable to read header.");
     545            return NULL;
     546        }
     547        status = 0;
     548        if (*bscale == 0.0) {
     549            psError(PS_ERR_IO, true,
     550                    "Supposed to use old values of BSCALE and BZERO, but they don't exist.");
     551            return NULL;
     552        }
     553        return scaleImageWrite((psImage*)image, fits->bitpix, *bscale, *bzero, rng);
     554    }
     555
     556    // Choose the appropriate output type, given the input type and desired bits per pixel
     557#define CONVERT_TYPE_INT_CASE(OUTTYPE, INTYPE, BITPIX) \
     558  case BITPIX: \
     559    OUTTYPE = PS_IS_PSELEMTYPE_UNSIGNED(INTYPE) ? PS_TYPE_U##BITPIX : PS_TYPE_S##BITPIX; \
     560    break;
     561#define CONVERT_TYPE_FLOAT_CASE(OUTTYPE, BITPIX) \
     562  case -BITPIX: /* Note the use of the negative sign */ \
     563    OUTTYPE = PS_TYPE_F##BITPIX; \
     564    break;
     565
     566    *bscale = 1.0;
     567    *bzero = 0.0;
     568    psElemType inType = image->type.type; // Type for input image
     569    psElemType outType;                 // Type for output image
     570    switch (fits->bitpix) {
     571        CONVERT_TYPE_INT_CASE(outType, inType, 8);
     572        CONVERT_TYPE_INT_CASE(outType, inType, 16);
     573        CONVERT_TYPE_INT_CASE(outType, inType, 32);
     574        CONVERT_TYPE_INT_CASE(outType, inType, 64);
     575        CONVERT_TYPE_FLOAT_CASE(outType, 32);
     576        CONVERT_TYPE_FLOAT_CASE(outType, 64);
     577      default:
     578        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target bitpix (%d) is not one of 8,16,32,64",
     579                fits->bitpix);
     580        return NULL;
     581    }
     582
     583    if (outType == inType) {
     584        return psMemIncrRefCounter((psImage*)image);
     585    }
     586    return psImageCopy(NULL, image, outType);
    213587}
    214588
    … …  
    232606    if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel,
    233607                         info->increment, NULL, output->data.V[0], &anynull, &status) != 0) {
    234         char fitsErr[80];               // CFITSIO error message string
    235         (void)fits_get_errstatus(status, fitsErr);
    236         psError(PS_ERR_IO, true,
    237                 _("Reading FITS file failed. CFITSIO Error: %s"),
    238                 fitsErr);
    239         return false;
    240     }
     608        psFitsError(status, true, "Reading FITS file failed.");
     609        return false;
     610    }
     611
     612    // No need to apply the BSCALE, BZERO because cfitsio does this for us
    241613
    242614    return true;
    … …  
    259631    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z);
    260632
    261     psImage *output = psImageAlloc(info->lastPixel[0] - info->firstPixel[0] + 1,
    262                                    info->lastPixel[1] - info->firstPixel[1] + 1,
    263                                    info->psDatatype);
    264 
    265     if (!fitsReadImage(output, fits, info)) {
     633    // Size of image
     634    int numCols = info->lastPixel[0] - info->firstPixel[0] + 1;
     635    int numRows = info->lastPixel[1] - info->firstPixel[1] + 1;
     636
     637    psImage *inImage = psImageAlloc(numCols, numRows, info->psDatatype); // Image to read in
     638
     639    psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point
     640    psImage *outImage = (floatType == PS_FITS_FLOAT_NONE ? psMemIncrRefCounter(inImage) :
     641                         psImageAlloc(numCols, numRows, psFitsFloatImageType(floatType))); // Output image
     642
     643    if (!fitsReadImage(inImage, fits, info)) {
    266644        psFree(info);
    267         psFree(output);
     645        psFree(inImage);
    268646        return NULL;
    269647    }
    270 
    271648    psFree(info);
    272     return output;
    273 }
    274 
    275 psImage* psFitsReadImageBuffer(psImage *output, // Output image buffer
    276                                const psFits *fits,    // the psFits object
     649
     650    if (floatType != PS_FITS_FLOAT_NONE) {
     651        outImage = psFitsFloatImageRead(outImage, inImage, floatType);
     652    }
     653    psFree(inImage);
     654
     655    return outImage;
     656}
     657
     658psImage* psFitsReadImageBuffer(psImage *outImage, // Output image buffer
     659                               const psFits *fits, // the psFits object
    277660                               psRegion region, // the region in the FITS image to read
    278661                               int z           // the z-plane in the FITS image cube to read
    … …  
    282665    PS_ASSERT_INT_NONNEGATIVE(z, NULL);
    283666
    284     if (output && output->parent) {
     667    if (outImage && outImage->parent) {
    285668        psError(PS_ERR_IO, true, "Unable to read into a buffer for a child image.\n");
    286669        return NULL;
    … …  
    295678    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z);
    296679
    297     output = psImageRecycle(output, info->lastPixel[0] - info->firstPixel[0] + 1,
    298                             info->lastPixel[1] - info->firstPixel[1] + 1,
    299                             info->psDatatype);
    300 
    301     if (!fitsReadImage(output, fits, info)) {
     680    // Size of image
     681    int numCols = info->lastPixel[0] - info->firstPixel[0] + 1;
     682    int numRows = info->lastPixel[1] - info->firstPixel[1] + 1;
     683
     684    psFitsFloat floatType = psFitsFloatImageCheck(fits); // Type of custom floating-point
     685    psImage *inImage;                   // Image to read in
     686    if (floatType == PS_FITS_FLOAT_NONE) {
     687        inImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype);
     688        outImage = psMemIncrRefCounter(inImage);
     689    } else {
     690        inImage = psImageAlloc(numCols, numRows, info->psDatatype);
     691        outImage = psImageRecycle(outImage, numCols, numRows, psFitsFloatImageType(floatType));
     692    }
     693
     694    if (!fitsReadImage(inImage, fits, info)) {
    302695        psFree(info);
    303         psFree(output);
     696        psFree(inImage);
     697        psFree(outImage);
    304698        return NULL;
    305699    }
    306 
    307700    psFree(info);
    308     return output;
     701
     702    if (floatType != PS_FITS_FLOAT_NONE) {
     703        outImage = psFitsFloatImageRead(outImage, inImage, floatType);
     704    }
     705    psFree(inImage);
     706
     707    return outImage;
    309708}
    310709
    311710bool psFitsWriteImage(psFits* fits,
    312                       const psMetadata* header,
     711                      psMetadata* header,
    313712                      const psImage* input,
    314713                      int numZPlanes,
    … …  
    324723}
    325724
    326 bool psFitsInsertImage(psFits* fits,
    327                        const psMetadata* header,
    328                        const psImage* input,
    329                        int numZPlanes,
    330                        const char* extname,
    331                        bool after)
     725bool psFitsInsertImage(psFits* fits, psMetadata* header, const psImage* image, int numZPlanes,
     726                       const char* extname, bool after)
    332727{
    333728    PS_ASSERT_FITS_NON_NULL(fits, false);
    334729    PS_ASSERT_FITS_WRITABLE(fits, false);
    335     PS_ASSERT_IMAGE_NON_NULL(input, false);
    336 
    337     int numCols = input->numCols;       // Number of columns for image
    338     int numRows = input->numRows;       // Number of rows for image
     730    PS_ASSERT_IMAGE_NON_NULL(image, false);
     731
     732    int numCols = image->numCols;       // Number of columns for image
     733    int numRows = image->numRows;       // Number of rows for image
    339734    int status = 0;                     // Status from cfitsio
     735
     736    double bscale = 0.0, bzero = 0.0;   // Scale and zero point to put in header (*already* applied to data)
     737    psFitsFloat floatType;              // Custom floating-point convention type
     738    psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, image,
     739                                           NULL, true); // Image to write out
     740    if (!diskImage) {
     741        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
     742        return false;
     743    }
    340744
    341745    // determine the FITS-equivalent parameters
    342746    int bitPix;                         // Bits per pixel
    343     double bZero;                       // Zero offset
     747    double cfitsioBzero = 0.0;          // Zero point for cfitsio to apply
    344748    int dataType;                       // cfitsio data type
    345     if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) {
    346         return false;
    347     }
     749    if (!p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) {
     750        psFree(diskImage);
     751        return false;
     752    }
     753    if (cfitsioBzero != 0.0) {
     754        assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash
     755        bscale = 1.0;
     756        bzero = cfitsioBzero;
     757    }
     758    assert(bitPix == fits->bitpix || fits->bitpix == 0);
    348759
    349760    int naxis = 3;                      // Number of axes
    … …  
    377788    }
    378789
     790    // write the header, if any.
    379791    if (header && !psFitsWriteHeader(fits, header)) {
    380792        psError(PS_ERR_IO, false, "Unable to write FITS header.\n");
    381         return false;
    382     }
    383 
    384     if (bZero != 0) {        // set the bscale/bzero
    385         fits_write_key_dbl(fits->fd, "BZERO", bZero, 12, "Pixel Value Offset", &status);
    386         fits_write_key_dbl(fits->fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status);
    387         fits_set_bscale(fits->fd, 1.0, bZero, &status);
    388     }
    389 
    390     // write the header, if any.
     793        psFree(diskImage);
     794        return false;
     795    }
     796
     797    // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
     798    // an unsigned integer type).  In all other cases, we have already converted the image to use the
     799    // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
     800    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
     801
     802    if (bscale != 0.0) {
     803        fits_write_key_dbl(fits->fd, "BZERO", bzero, 12,
     804                           "Scaling: TRUE = BZERO + BSCALE * DISK", &status);
     805        fits_write_key_dbl(fits->fd, "BSCALE", bscale, 12,
     806                           "Scaling: TRUE = BZERO + BSCALE * DISK", &status);
     807        if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) {
     808            psFree(diskImage);
     809            return false;
     810        }
     811    }
     812
     813    if (floatType != PS_FITS_FLOAT_NONE) {
     814        psFitsFloatImageSet(fits, floatType);
     815    }
     816
    391817    if (extname && strlen(extname) > 0) {
    392         psFitsSetExtName(fits,extname);
    393     }
    394 
    395     if (input->parent == NULL) { // if no parent, assume that the image data is contiguous
    396         fits_write_img(fits->fd,
    397                        dataType,              // datatype
    398                        1,                     // writing to the first z-plane
    399                        numCols*numRows,       // number of elements to write, i.e., the whole image
    400                        input->data.V[0],      // the data
    401                        &status);
    402     } else { // image data may not be contiguous; write one row at a time
     818        psFitsSetExtName(fits, extname);
     819    }
     820
     821    if (image->parent == NULL) {
     822        // if no parent, assume that the image data is contiguous
     823        fits_write_img(fits->fd, dataType, 1, numCols*numRows, diskImage->data.V[0], &status);
     824    } else {
     825        // image data may not be contiguous; write one row at a time
    403826        int firstPixel = 1;
    404827        for (int row = 0; row < numRows; row++) {
    405             fits_write_img(fits->fd,
    406                            dataType,          // datatype
    407                            firstPixel,
    408                            numCols,           // number of elements to write, i.e., one row's worth
    409                            input->data.V[row],// the raw row data
    410                            &status);
    411             firstPixel += numCols;  // move to next row
    412         }
    413     }
    414 
    415     if (status != 0) {
    416         char fitsErr[MAX_STRING_LENGTH];
    417         (void)fits_get_errstatus(status, fitsErr);
    418         psError(PS_ERR_IO, true,
    419                 _("Could not write data to file. CFITSIO Error: %s"),
    420                 fitsErr);
     828            fits_write_img(fits->fd, dataType, firstPixel, numCols, diskImage->data.V[row], &status);
     829            firstPixel += numCols;
     830        }
     831    }
     832
     833    psFree(diskImage);
     834
     835    if (psFitsError(status, true, "Could not write image to file.")) {
    421836        return false;
    422837    }
    … …  
    426841}
    427842
    428 bool psFitsUpdateImage(psFits* fits,
    429                        const psImage* input,
    430                        int x0,
    431                        int y0,
    432                        int z)
     843bool psFitsUpdateImage(psFits* fits, const psImage* input, int x0, int y0, int z)
    433844{
    434845    PS_ASSERT_FITS_NON_NULL(fits, false);
    … …  
    440851    // check to see if we are positioned on an image HDU
    441852    int hduType;
    442     if ( fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
    443         char fitsErr[MAX_STRING_LENGTH];
    444         (void)fits_get_errstatus(status, fitsErr);
    445         psError(PS_ERR_IO, true,
    446                 _("Could not determine the HDU type. CFITSIO Error: %s"),
    447                 fitsErr);
     853    if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) {
     854        psFitsError(status, true, "Could not determine the HDU type.");
    448855        return NULL;
    449856    }
    450857    if (hduType != IMAGE_HDU) {
    451         psError(PS_ERR_IO, true,
    452                 _("Current FITS HDU type must be an image."));
     858        psError(PS_ERR_IO, true, _("Current FITS HDU type must be an image."));
    453859        return NULL;
    454860    }
    … …  
    457863    int numRows = input->numRows;
    458864
     865    double bscale = 0.0, bzero = 0.0;   // Scale and zero point to put in header (*already* applied to data)
     866    psFitsFloat floatType;              // Custom floating-point convention type
     867    psImage *diskImage = convertImageWrite(&bscale, &bzero, &floatType, fits, input,
     868                                           NULL, false); // Image to write out
     869    if (!diskImage) {
     870        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
     871        return false;
     872    }
     873
    459874    // determine the FITS-equivalent parameters
    460     int bitPix;
    461     double bZero;
    462     int dataType;
    463     if (! p_psFitsTypeToCfitsio(input->type.type, &bitPix, &bZero, &dataType) ) {
    464         return false;
    465     }
     875    int bitPix;                         // Bits per pixel
     876    double cfitsioBzero = 0.0;          // Zero point for cfitsio to apply
     877    int dataType;                       // cfitsio data type
     878    if (! p_psFitsTypeToCfitsio(diskImage->type.type, &bitPix, &cfitsioBzero, &dataType)) {
     879        psFree(diskImage);
     880        return false;
     881    }
     882    if (cfitsioBzero != 0.0) {
     883        assert(bzero == 0.0 && bscale == 1.0); // p_psFitsTypeToCfitsio and convertImageWrite must not clash
     884        bscale = 1.0;
     885        bzero = cfitsioBzero;
     886    }
     887    assert(bitPix == fits->bitpix || fits->bitpix == 0);
    466888
    467889    //check to see if the HDU has the same datatype
    … …  
    476898        char* fitsTypeStr;
    477899        char* imageTypeStr;
    478         PS_TYPE_NAME(fitsTypeStr,fileBitpix);
    479         PS_TYPE_NAME(imageTypeStr,input->type.type);
    480         psError(PS_ERR_IO, true,
    481                 _("Can not update a %s image given a %s image."),
    482                 fitsTypeStr, imageTypeStr);
     900        PS_TYPE_NAME(fitsTypeStr, fileBitpix);
     901        PS_TYPE_NAME(imageTypeStr, input->type.type);
     902        psError(PS_ERR_IO, true, _("Can not update a %s image given a %s image."), fitsTypeStr, imageTypeStr);
     903        psFree(diskImage);
    483904        return false;
    484905    }
    … …  
    487908    if (z >= nAxes[2]) {
    488909        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    489                 _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."),
    490                 nAxes[2], z);
     910                _("Current FITS HDU has %ld z-planes, but z-plane %d was specified."), nAxes[2], z);
     911        psFree(diskImage);
    491912        return false;
    492913    }
    … …  
    505926
    506927    if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] ||
    507             firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
    508             lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
    509             lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
     928        firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
     929        lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
     930        lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
    510931        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    511932                "Input image [size of %ix%i] at position (%i,%i) does not all lay in the %lix%li FITS image.",
    512                 numCols, numRows,
    513                 x0, y0,
    514                 nAxes[0], nAxes[1]);
    515         return false;
    516     }
    517 
    518     fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, input->data.V[0], &status);
    519 
    520     if ( status != 0) {
    521         char fitsErr[MAX_STRING_LENGTH];
    522         (void)fits_get_errstatus(status, fitsErr);
    523         psError(PS_ERR_IO, true,
    524                 _("Could not write data to file. CFITSIO Error: %s"),
    525                 fitsErr);
     933                numCols, numRows, x0, y0, nAxes[0], nAxes[1]);
     934        psFree(diskImage);
     935        return false;
     936    }
     937
     938    // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
     939    // an unsigned integer type).  In all other cases, we have already converted the image to use the
     940    // appropriate scale and zero (because we want to apply a randomiser to the quantisation).
     941    fits_set_bscale(fits->fd, 1.0, cfitsioBzero, &status);
     942
     943    fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, diskImage->data.V[0], &status);
     944
     945    psFree(diskImage);
     946
     947    if (psFitsError(status, true, "Could not write data to file.")) {
    526948        return false;
    527949    }
    … …  
    537959    long nAxes[3];                      // Number of pixels on each axis
    538960    int status = 0;                     // cfitsio status value
    539     char fitsErr[80] = "";              // CFITSIO error message string
    540961
    541962    // Some of this replicates what is in psFitsReadImage, so it's a little inefficient.  But it saves
    … …  
    549970
    550971    if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) {
    551         (void)fits_get_errstatus(status, fitsErr);
    552         psError(PS_ERR_IO, true,
    553                 _("Could not determine image dimensions. CFITSIO Error: %s"),
    554                 fitsErr);
     972        psFitsError(status, true, "Could not determine image dimensions.");
    555973        return NULL;
    556974    }
    … …  
    563981    if (nAxis == 3) {
    564982        if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) {
    565             (void)fits_get_errstatus(status, fitsErr);
    566             psError(PS_ERR_IO, true,
    567                     _("Could not determine image size. CFITSIO Error: %s"),
    568                     fitsErr);
     983            psFitsError(status, true, "Could not determine image size.");
    569984            return NULL;
    570985        }
    … …  
    579994
    580995    // Bad dimensionality
    581     psError(PS_ERR_IO, true,
    582             _("Image number of dimensions, %d, is not valid."
    583               " Only two or three dimensions supported for FITS I/O."), nAxis);
     996    psError(PS_ERR_IO, true, _("Image number of dimensions, %d, is not supported."), nAxis);
    584997    return NULL;
    585998}
Note: See TracChangeset for help on using the changeset viewer.