"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/imcompress.c" (20 May 2021, 380829 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "imcompress.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 3430_vs_3440.

    1 # include <stdio.h>
    2 # include <stdlib.h>
    3 # include <string.h>
    4 # include <math.h>
    5 # include <ctype.h>
    6 # include <time.h>
    7 # include "fitsio2.h"
    8 
    9 #define NULL_VALUE -2147483647 /* value used to represent undefined pixels */
   10 #define ZERO_VALUE -2147483646 /* value used to represent zero-valued pixels */
   11 
   12 /* nearest integer function */
   13 # define NINT(x)  ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5))
   14 
   15 /* special quantize level value indicates that floating point image pixels */
   16 /* should not be quantized and instead losslessly compressed (with GZIP) */
   17 #define NO_QUANTIZE 9999
   18 
   19 /* string array for storing the individual column compression stats */
   20 char results[999][30];
   21 
   22 float *fits_rand_value = 0;
   23 
   24 int imcomp_write_nocompress_tile(fitsfile *outfptr, long row, int datatype, 
   25     void *tiledata, long tilelen, int nullcheck, void *nullflagval, int *status);
   26 int imcomp_convert_tile_tshort(fitsfile *outfptr, void *tiledata, long tilelen,
   27     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   28     double zero, double actual_bzero, int *intlength, int *status);
   29 int imcomp_convert_tile_tushort(fitsfile *outfptr, void *tiledata, long tilelen,
   30     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   31     double zero, int *intlength, int *status);
   32 int imcomp_convert_tile_tint(fitsfile *outfptr, void *tiledata, long tilelen,
   33     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   34     double zero, int *intlength, int *status);
   35 int imcomp_convert_tile_tuint(fitsfile *outfptr, void *tiledata, long tilelen,
   36     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   37     double zero, int *intlength, int *status);
   38 int imcomp_convert_tile_tbyte(fitsfile *outfptr, void *tiledata, long tilelen,
   39     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   40     double zero, int *intlength, int *status);
   41 int imcomp_convert_tile_tsbyte(fitsfile *outfptr, void *tiledata, long tilelen,
   42     int nullcheck, void *nullflagval, int nullval, int zbitpix, double scale,
   43     double zero, int *intlength, int *status);
   44 int imcomp_convert_tile_tfloat(fitsfile *outfptr, long row, void *tiledata, long tilelen,
   45     long tilenx, long tileny, int nullcheck, void *nullflagval, int nullval, int zbitpix,
   46     double scale, double zero, int *intlength, int *flag, double *bscale, double *bzero,int *status);
   47 int imcomp_convert_tile_tdouble(fitsfile *outfptr, long row, void *tiledata, long tilelen,
   48     long tilenx, long tileny, int nullcheck, void *nullflagval, int nullval, int zbitpix, 
   49     double scale, double zero, int *intlength, int *flag, double *bscale, double *bzero, int *status);
   50 
   51 static int unquantize_i1r4(long row,
   52             unsigned char *input,         /* I - array of values to be converted     */
   53             long ntodo,           /* I - number of elements in the array     */
   54             double scale,         /* I - FITS TSCALn or BSCALE value         */
   55             double zero,          /* I - FITS TZEROn or BZERO  value         */
   56             int dither_method,    /* I - which subtractive dither method to use */
   57             int nullcheck,        /* I - null checking code; 0 = don't check */
   58                                   /*     1:set null pixels = nullval         */
   59                                   /*     2: if null pixel, set nullarray = 1 */
   60             unsigned char tnull,          /* I - value of FITS TNULLn keyword if any */
   61             float nullval,        /* I - set null pixels, if nullcheck = 1   */
   62             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
   63             int  *anynull,        /* O - set to 1 if any pixels are null     */
   64             float *output,        /* O - array of converted pixels           */
   65             int *status);          /* IO - error status                       */
   66 static int unquantize_i2r4(long row,
   67             short *input,         /* I - array of values to be converted     */
   68             long ntodo,           /* I - number of elements in the array     */
   69             double scale,         /* I - FITS TSCALn or BSCALE value         */
   70             double zero,          /* I - FITS TZEROn or BZERO  value         */
   71             int dither_method,    /* I - which subtractive dither method to use */
   72             int nullcheck,        /* I - null checking code; 0 = don't check */
   73                                   /*     1:set null pixels = nullval         */
   74                                   /*     2: if null pixel, set nullarray = 1 */
   75             short tnull,          /* I - value of FITS TNULLn keyword if any */
   76             float nullval,        /* I - set null pixels, if nullcheck = 1   */
   77             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
   78             int  *anynull,        /* O - set to 1 if any pixels are null     */
   79             float *output,        /* O - array of converted pixels           */
   80             int *status);          /* IO - error status                       */
   81 static int unquantize_i4r4(long row,
   82             INT32BIT *input,      /* I - array of values to be converted     */
   83             long ntodo,           /* I - number of elements in the array     */
   84             double scale,         /* I - FITS TSCALn or BSCALE value         */
   85             double zero,          /* I - FITS TZEROn or BZERO  value         */
   86             int dither_method,    /* I - which subtractive dither method to use */
   87             int nullcheck,        /* I - null checking code; 0 = don't check */
   88                                   /*     1:set null pixels = nullval         */
   89                                   /*     2: if null pixel, set nullarray = 1 */
   90             INT32BIT tnull,       /* I - value of FITS TNULLn keyword if any */
   91             float nullval,        /* I - set null pixels, if nullcheck = 1   */
   92             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
   93             int  *anynull,        /* O - set to 1 if any pixels are null     */
   94             float *output,        /* O - array of converted pixels           */
   95             int *status);          /* IO - error status                       */
   96 static int unquantize_i1r8(long row,
   97             unsigned char *input,         /* I - array of values to be converted     */
   98             long ntodo,           /* I - number of elements in the array     */
   99             double scale,         /* I - FITS TSCALn or BSCALE value         */
  100             double zero,          /* I - FITS TZEROn or BZERO  value         */
  101             int dither_method,    /* I - which subtractive dither method to use */
  102             int nullcheck,        /* I - null checking code; 0 = don't check */
  103                                   /*     1:set null pixels = nullval         */
  104                                   /*     2: if null pixel, set nullarray = 1 */
  105             unsigned char tnull,          /* I - value of FITS TNULLn keyword if any */
  106             double nullval,        /* I - set null pixels, if nullcheck = 1   */
  107             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
  108             int  *anynull,        /* O - set to 1 if any pixels are null     */
  109             double *output,        /* O - array of converted pixels           */
  110             int *status);          /* IO - error status                       */
  111 static int unquantize_i2r8(long row,
  112             short *input,         /* I - array of values to be converted     */
  113             long ntodo,           /* I - number of elements in the array     */
  114             double scale,         /* I - FITS TSCALn or BSCALE value         */
  115             double zero,          /* I - FITS TZEROn or BZERO  value         */
  116             int dither_method,    /* I - which subtractive dither method to use */
  117             int nullcheck,        /* I - null checking code; 0 = don't check */
  118                                   /*     1:set null pixels = nullval         */
  119                                   /*     2: if null pixel, set nullarray = 1 */
  120             short tnull,          /* I - value of FITS TNULLn keyword if any */
  121             double nullval,        /* I - set null pixels, if nullcheck = 1   */
  122             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
  123             int  *anynull,        /* O - set to 1 if any pixels are null     */
  124             double *output,        /* O - array of converted pixels           */
  125             int *status);          /* IO - error status                       */
  126 static int unquantize_i4r8(long row,
  127             INT32BIT *input,      /* I - array of values to be converted     */
  128             long ntodo,           /* I - number of elements in the array     */
  129             double scale,         /* I - FITS TSCALn or BSCALE value         */
  130             double zero,          /* I - FITS TZEROn or BZERO  value         */
  131             int dither_method,    /* I - which subtractive dither method to use */
  132             int nullcheck,        /* I - null checking code; 0 = don't check */
  133                                   /*     1:set null pixels = nullval         */
  134                                   /*     2: if null pixel, set nullarray = 1 */
  135             INT32BIT tnull,       /* I - value of FITS TNULLn keyword if any */
  136             double nullval,        /* I - set null pixels, if nullcheck = 1   */
  137             char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
  138             int  *anynull,        /* O - set to 1 if any pixels are null     */
  139             double *output,        /* O - array of converted pixels           */
  140             int *status);          /* IO - error status                       */
  141 static int imcomp_float2nan(float *indata, long tilelen, int *outdata,
  142     float nullflagval,  int *status);
  143 static int imcomp_double2nan(double *indata, long tilelen, LONGLONG *outdata,
  144     double nullflagval,  int *status);    
  145 static int fits_read_write_compressed_img(fitsfile *fptr,   /* I - FITS file pointer */
  146             int  datatype,  /* I - datatype of the array to be returned      */
  147             LONGLONG  *infpixel, /* I - 'bottom left corner' of the subsection    */
  148             LONGLONG  *inlpixel, /* I - 'top right corner' of the subsection      */
  149             long  *ininc,    /* I - increment to be applied in each dimension */
  150             int  nullcheck,  /* I - 0 for no null checking                   */
  151                               /*     1: set undefined pixels = nullval       */
  152             void *nullval,    /* I - value for undefined pixels              */
  153             int  *anynul,     /* O - set to 1 if any values are null; else 0 */
  154             fitsfile *outfptr,   /* I - FITS file pointer                    */
  155             int  *status);
  156 
  157 static int fits_shuffle_8bytes(char *heap, LONGLONG length, int *status);
  158 static int fits_shuffle_4bytes(char *heap, LONGLONG length, int *status);
  159 static int fits_shuffle_2bytes(char *heap, LONGLONG length, int *status);
  160 static int fits_unshuffle_8bytes(char *heap, LONGLONG length, int *status);
  161 static int fits_unshuffle_4bytes(char *heap, LONGLONG length, int *status);
  162 static int fits_unshuffle_2bytes(char *heap, LONGLONG length, int *status);
  163 
  164 static int fits_int_to_longlong_inplace(int *intarray, long length, int *status);
  165 static int fits_short_to_int_inplace(short *intarray, long length, int shift, int *status);
  166 static int fits_ushort_to_int_inplace(unsigned short *intarray, long length, int shift, int *status);
  167 static int fits_sbyte_to_int_inplace(signed char *intarray, long length, int *status);
  168 static int fits_ubyte_to_int_inplace(unsigned char *intarray, long length, int *status);
  169 
  170 static int fits_calc_tile_rows(long *tlpixel, long *tfpixel, int ndim, long *trowsize, long *ntrows, int *status); 
  171 
  172 /* only used for diagnoitic purposes */
  173 /* int fits_get_case(int *c1, int*c2, int*c3); */ 
  174 /*---------------------------------------------------------------------------*/
  175 int fits_init_randoms(void) {
  176 
  177 /* initialize an array of random numbers */
  178 
  179     int ii;
  180     double a = 16807.0;
  181     double m = 2147483647.0;
  182     double temp, seed;
  183 
  184     FFLOCK;
  185  
  186     if (fits_rand_value) {
  187        FFUNLOCK;
  188        return(0);  /* array is already initialized */
  189     }
  190 
  191     /* allocate array for the random number sequence */
  192     /* THIS MEMORY IS NEVER FREED */
  193     fits_rand_value = calloc(N_RANDOM, sizeof(float));
  194 
  195     if (!fits_rand_value) {
  196         FFUNLOCK;
  197     return(MEMORY_ALLOCATION);
  198     }
  199                
  200     /*  We need a portable algorithm that anyone can use to generate this
  201         exact same sequence of random number.  The C 'rand' function is not
  202     suitable because it is not available to Fortran or Java programmers.
  203     Instead, use a well known simple algorithm published here: 
  204     "Random number generators: good ones are hard to find", Communications of the ACM,
  205         Volume 31 ,  Issue 10  (October 1988) Pages: 1192 - 1201 
  206     */  
  207 
  208     /* initialize the random numbers */
  209     seed = 1;
  210     for (ii = 0; ii < N_RANDOM; ii++) {
  211         temp = a * seed;
  212     seed = temp -m * ((int) (temp / m) );
  213     fits_rand_value[ii] = (float) (seed / m);
  214     }
  215 
  216     FFUNLOCK;
  217 
  218     /* 
  219     IMPORTANT NOTE: the 10000th seed value must have the value 1043618065 if the 
  220        algorithm has been implemented correctly */
  221     
  222     if ( (int) seed != 1043618065) {
  223         ffpmsg("fits_init_randoms generated incorrect random number sequence");
  224     return(1);
  225     } else {
  226         return(0);
  227     }
  228 }
  229 /*--------------------------------------------------------------------------*/
  230 void bz_internal_error(int errcode)
  231 {
  232     /* external function declared by the bzip2 code in bzlib_private.h */
  233     ffpmsg("bzip2 returned an internal error");
  234     ffpmsg("This should never happen");
  235     return;
  236 }
  237 /*--------------------------------------------------------------------------*/
  238 int fits_set_compression_type(fitsfile *fptr,  /* I - FITS file pointer     */
  239        int ctype,    /* image compression type code;                        */
  240                      /* allowed values: RICE_1, GZIP_1, GZIP_2, PLIO_1,     */
  241                      /*  HCOMPRESS_1, BZIP2_1, and NOCOMPRESS               */
  242        int *status)  /* IO - error status                                   */
  243 {
  244 /*
  245    This routine specifies the image compression algorithm that should be
  246    used when writing a FITS image.  The image is divided into tiles, and
  247    each tile is compressed and stored in a row of at variable length binary
  248    table column.
  249 */
  250 
  251     if (ctype != RICE_1 && 
  252         ctype != GZIP_1 && 
  253         ctype != GZIP_2 && 
  254         ctype != PLIO_1 && 
  255         ctype != HCOMPRESS_1 && 
  256         ctype != BZIP2_1 && 
  257         ctype != NOCOMPRESS &&
  258     ctype != 0)
  259     {
  260     ffpmsg("unknown compression algorithm (fits_set_compression_type)");
  261     *status = DATA_COMPRESSION_ERR; 
  262     } else {
  263         (fptr->Fptr)->request_compress_type = ctype;
  264     }
  265     return(*status);
  266 }
  267 /*--------------------------------------------------------------------------*/
  268 int fits_set_tile_dim(fitsfile *fptr,  /* I - FITS file pointer             */
  269            int ndim,   /* number of dimensions in the compressed image      */
  270            long *dims, /* size of image compression tile in each dimension  */
  271                       /* default tile size = (NAXIS1, 1, 1, ...)            */
  272            int *status)         /* IO - error status                        */
  273 {
  274 /*
  275    This routine specifies the size (dimension) of the image
  276    compression  tiles that should be used when writing a FITS
  277    image.  The image is divided into tiles, and each tile is compressed
  278    and stored in a row of at variable length binary table column.
  279 */
  280     int ii;
  281 
  282     if (ndim < 0 || ndim > MAX_COMPRESS_DIM)
  283     {
  284         *status = BAD_DIMEN;
  285     ffpmsg("illegal number of tile dimensions (fits_set_tile_dim)");
  286         return(*status);
  287     }
  288 
  289     for (ii = 0; ii < ndim; ii++)
  290     {
  291         (fptr->Fptr)->request_tilesize[ii] = dims[ii];
  292     }
  293 
  294     return(*status);
  295 }
  296 /*--------------------------------------------------------------------------*/
  297 int fits_set_quantize_level(fitsfile *fptr,  /* I - FITS file pointer   */
  298            float qlevel,        /* floating point quantization level      */
  299            int *status)         /* IO - error status                */
  300 {
  301 /*
  302    This routine specifies the value of the quantization level, q,  that
  303    should be used when compressing floating point images.  The image is
  304    divided into tiles, and each tile is compressed and stored in a row
  305    of at variable length binary table column.
  306 */
  307     if (qlevel == 0.)
  308     {
  309         /* this means don't quantize the floating point values. Instead, */
  310     /* the floating point values will be losslessly compressed */
  311        (fptr->Fptr)->request_quantize_level = NO_QUANTIZE;
  312     } else {
  313 
  314         (fptr->Fptr)->request_quantize_level = qlevel;
  315     }
  316 
  317     return(*status);
  318 }
  319 /*--------------------------------------------------------------------------*/
  320 int fits_set_quantize_method(fitsfile *fptr,  /* I - FITS file pointer   */
  321            int method,          /* quantization method       */
  322            int *status)         /* IO - error status                */
  323 {
  324 /*
  325    This routine specifies what type of dithering (randomization) should
  326    be performed when quantizing floating point images to integer prior to
  327    compression.   A value of -1 means do no dithering.  A value of 0 means
  328    use the default SUBTRACTIVE_DITHER_1 (which is equivalent to dither = 1).
  329    A value of 2 means use SUBTRACTIVE_DITHER_2.
  330 */
  331 
  332     if (method < -1 || method > 2)
  333     {
  334     ffpmsg("illegal dithering value (fits_set_quantize_method)");
  335     *status = DATA_COMPRESSION_ERR; 
  336     } else {
  337        
  338         if (method == 0) method = 1;
  339         (fptr->Fptr)->request_quantize_method = method;
  340     }
  341 
  342     return(*status);
  343 }
  344 /*--------------------------------------------------------------------------*/
  345 int fits_set_quantize_dither(fitsfile *fptr,  /* I - FITS file pointer   */
  346            int dither,        /* dither type      */
  347            int *status)         /* IO - error status                */
  348 {
  349 /*
  350    the name of this routine has changed.  This is kept here only for backwards
  351    compatibility for any software that may be calling the old routine.
  352 */
  353 
  354     fits_set_quantize_method(fptr, dither, status);
  355     return(*status);
  356 }
  357 /*--------------------------------------------------------------------------*/
  358 int fits_set_dither_seed(fitsfile *fptr,  /* I - FITS file pointer   */
  359            int seed,        /* random dithering seed value (1 to 10000) */
  360            int *status)         /* IO - error status                */
  361 {
  362 /*
  363    This routine specifies the value of the offset that should be applied when
  364    calculating the random dithering when quantizing floating point iamges.
  365    A random offset should be applied to each image to avoid quantization 
  366    effects when taking the difference of 2 images, or co-adding a set of
  367    images.  Without this random offset, the corresponding pixel in every image
  368    will have exactly the same dithering.
  369    
  370    offset = 0 means use the default random dithering based on system time
  371    offset = negative means randomly chose dithering based on 1st tile checksum
  372    offset = [1 - 10000] means use that particular dithering pattern
  373 
  374 */
  375     /* if positive, ensure that the value is in the range 1 to 10000 */
  376     if (seed > 10000) {
  377     ffpmsg("illegal dithering seed value (fits_set_dither_seed)");
  378     *status = DATA_COMPRESSION_ERR;
  379     } else {
  380        (fptr->Fptr)->request_dither_seed = seed; 
  381     }
  382     
  383     return(*status);
  384 }
  385 /*--------------------------------------------------------------------------*/
  386 int fits_set_dither_offset(fitsfile *fptr,  /* I - FITS file pointer   */
  387            int offset,        /* random dithering offset value (1 to 10000) */
  388            int *status)         /* IO - error status                */
  389 {
  390 /*
  391     The name of this routine has changed.  This is kept just for
  392     backwards compatibility with any software that calls the old name
  393 */
  394 
  395     fits_set_dither_seed(fptr, offset, status);
  396     return(*status);
  397 }
  398 /*--------------------------------------------------------------------------*/
  399 int fits_set_noise_bits(fitsfile *fptr,  /* I - FITS file pointer   */
  400            int noisebits,       /* noise_bits parameter value       */
  401                                 /* (default = 4)                    */
  402            int *status)         /* IO - error status                */
  403 {
  404 /*
  405    ********************************************************************
  406    ********************************************************************
  407    THIS ROUTINE IS PROVIDED ONLY FOR BACKWARDS COMPATIBILITY;
  408    ALL NEW SOFTWARE SHOULD CALL fits_set_quantize_level INSTEAD
  409    ********************************************************************
  410    ********************************************************************
  411 
  412    This routine specifies the value of the noice_bits parameter that
  413    should be used when compressing floating point images.  The image is
  414    divided into tiles, and each tile is compressed and stored in a row
  415    of at variable length binary table column.
  416 
  417    Feb 2008:  the "noisebits" parameter has been replaced with the more
  418    general "quantize level" parameter.
  419 */
  420     float qlevel;
  421 
  422     if (noisebits < 1 || noisebits > 16)
  423     {
  424         *status = DATA_COMPRESSION_ERR;
  425     ffpmsg("illegal number of noise bits (fits_set_noise_bits)");
  426         return(*status);
  427     }
  428 
  429     qlevel = (float) pow (2., (double)noisebits);
  430     fits_set_quantize_level(fptr, qlevel, status);
  431     
  432     return(*status);
  433 }
  434 /*--------------------------------------------------------------------------*/
  435 int fits_set_hcomp_scale(fitsfile *fptr,  /* I - FITS file pointer   */
  436            float scale,       /* hcompress scale parameter value       */
  437                                 /* (default = 0.)                    */
  438            int *status)         /* IO - error status                */
  439 {
  440 /*
  441    This routine specifies the value of the hcompress scale parameter.
  442 */
  443     (fptr->Fptr)->request_hcomp_scale = scale;
  444     return(*status);
  445 }
  446 /*--------------------------------------------------------------------------*/
  447 int fits_set_hcomp_smooth(fitsfile *fptr,  /* I - FITS file pointer   */
  448            int smooth,       /* hcompress smooth parameter value       */
  449                                 /* if scale > 1 and smooth != 0, then */
  450                 /*  the image will be smoothed when it is */
  451                 /* decompressed to remove some of the */
  452                 /* 'blockiness' in the image produced */
  453                 /* by the lossy compression    */
  454            int *status)         /* IO - error status                */
  455 {
  456 /*
  457    This routine specifies the value of the hcompress scale parameter.
  458 */
  459 
  460     (fptr->Fptr)->request_hcomp_smooth = smooth;
  461     return(*status);
  462 }
  463 /*--------------------------------------------------------------------------*/
  464 int fits_set_lossy_int(fitsfile *fptr,  /* I - FITS file pointer   */
  465            int lossy_int,       /* I - True (!= 0) or False (0) */
  466            int *status)         /* IO - error status                */
  467 {
  468 /*
  469    This routine specifies whether images with integer pixel values should
  470    quantized and compressed the same way float images are compressed.
  471    The default is to not do this, and instead apply a lossless compression
  472    algorithm to integer images.
  473 */
  474 
  475     (fptr->Fptr)->request_lossy_int_compress = lossy_int;
  476     return(*status);
  477 }
  478 /*--------------------------------------------------------------------------*/
  479 int fits_set_huge_hdu(fitsfile *fptr,  /* I - FITS file pointer   */
  480            int huge,       /* I - True (!= 0) or False (0) */
  481            int *status)         /* IO - error status                */
  482 {
  483 /*
  484    This routine specifies whether the HDU that is being compressed is so large
  485    (i.e., > 4 GB) that the 'Q' type variable length array columns should be used
  486    rather than the normal 'P' type.  The allows the heap pointers to be stored
  487    as 64-bit quantities, rather than just 32-bits.
  488 */
  489 
  490     (fptr->Fptr)->request_huge_hdu = huge;
  491     return(*status);
  492 }
  493 /*--------------------------------------------------------------------------*/
  494 int fits_get_compression_type(fitsfile *fptr,  /* I - FITS file pointer     */
  495        int *ctype,   /* image compression type code;                        */
  496                      /* allowed values:                                     */
  497              /* RICE_1, GZIP_1, GZIP_2, PLIO_1, HCOMPRESS_1, BZIP2_1 */
  498        int *status)  /* IO - error status                                   */
  499 {
  500 /*
  501    This routine returns the image compression algorithm that should be
  502    used when writing a FITS image.  The image is divided into tiles, and
  503    each tile is compressed and stored in a row of at variable length binary
  504    table column.
  505 */
  506     *ctype = (fptr->Fptr)->request_compress_type;
  507 
  508     if (*ctype != RICE_1 && 
  509         *ctype != GZIP_1 && 
  510         *ctype != GZIP_2 && 
  511         *ctype != PLIO_1 && 
  512         *ctype != HCOMPRESS_1 && 
  513         *ctype != BZIP2_1 && 
  514         *ctype != NOCOMPRESS &&
  515     *ctype != 0   ) 
  516 
  517     {
  518     ffpmsg("unknown compression algorithm (fits_get_compression_type)");
  519     *status = DATA_COMPRESSION_ERR; 
  520     }
  521  
  522     return(*status);
  523 }
  524 /*--------------------------------------------------------------------------*/
  525 int fits_get_tile_dim(fitsfile *fptr,  /* I - FITS file pointer             */
  526            int ndim,   /* number of dimensions in the compressed image      */
  527            long *dims, /* size of image compression tile in each dimension  */
  528                        /* default tile size = (NAXIS1, 1, 1, ...)           */
  529            int *status)         /* IO - error status                        */
  530 {
  531 /*
  532    This routine returns the size (dimension) of the image
  533    compression  tiles that should be used when writing a FITS
  534    image.  The image is divided into tiles, and each tile is compressed
  535    and stored in a row of at variable length binary table column.
  536 */
  537     int ii;
  538 
  539     if (ndim < 0 || ndim > MAX_COMPRESS_DIM)
  540     {
  541         *status = BAD_DIMEN;
  542     ffpmsg("illegal number of tile dimensions (fits_get_tile_dim)");
  543         return(*status);
  544     }
  545 
  546     for (ii = 0; ii < ndim; ii++)
  547     {
  548         dims[ii] = (fptr->Fptr)->request_tilesize[ii];
  549     }
  550 
  551     return(*status);
  552 }
  553 /*--------------------------------------------------------------------------*/
  554 int fits_unset_compression_param(
  555       fitsfile *fptr,
  556       int *status) 
  557 {
  558     int ii;
  559 
  560     (fptr->Fptr)->compress_type = 0;
  561     (fptr->Fptr)->quantize_level = 0;
  562     (fptr->Fptr)->quantize_method = 0;
  563     (fptr->Fptr)->dither_seed = 0; 
  564     (fptr->Fptr)->hcomp_scale = 0;
  565 
  566     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
  567     {
  568         (fptr->Fptr)->tilesize[ii] = 0;
  569     }
  570 
  571     return(*status);
  572 }
  573 /*--------------------------------------------------------------------------*/
  574 int fits_unset_compression_request(
  575       fitsfile *fptr,
  576       int *status) 
  577 {
  578     int ii;
  579 
  580     (fptr->Fptr)->request_compress_type = 0;
  581     (fptr->Fptr)->request_quantize_level = 0;
  582     (fptr->Fptr)->request_quantize_method = 0;
  583     (fptr->Fptr)->request_dither_seed = 0; 
  584     (fptr->Fptr)->request_hcomp_scale = 0;
  585     (fptr->Fptr)->request_lossy_int_compress = 0;
  586     (fptr->Fptr)->request_huge_hdu = 0;
  587 
  588     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
  589     {
  590         (fptr->Fptr)->request_tilesize[ii] = 0;
  591     }
  592 
  593     return(*status);
  594 }
  595 /*--------------------------------------------------------------------------*/
  596 int fits_set_compression_pref(
  597       fitsfile *infptr,
  598       fitsfile *outfptr,
  599       int *status) 
  600 {
  601 /*
  602    Set the preference for various compression options, based
  603    on keywords in the input file that
  604    provide guidance about how the HDU should be compressed when written
  605    to the output file.
  606 */
  607 
  608     int ii, naxis, nkeys, comptype;
  609     int  ivalue;
  610     long tiledim[6]= {1,1,1,1,1,1};
  611     char card[FLEN_CARD], value[FLEN_VALUE];
  612     double  qvalue;
  613     float hscale;
  614     LONGLONG datastart, dataend; 
  615     if (*status > 0)
  616         return(*status);
  617 
  618     /* check the size of the HDU that is to be compressed */
  619     fits_get_hduaddrll(infptr, NULL, &datastart, &dataend, status);
  620     if ( (LONGLONG)(dataend - datastart) > UINT32_MAX) {
  621        /* use 64-bit '1Q' variable length columns instead of '1P' columns */
  622        /* for large files, in case the heap size becomes larger than 2**32 bytes*/
  623        fits_set_huge_hdu(outfptr, 1, status);
  624     }
  625 
  626     fits_get_hdrspace(infptr, &nkeys, NULL, status);
  627  
  628    /* look for a image compression directive keywords (begin with 'FZ') */
  629     for (ii = 2; ii <= nkeys; ii++) {
  630         
  631     fits_read_record(infptr, ii, card, status);
  632 
  633     if (!strncmp(card, "FZ", 2) ){
  634     
  635             /* get the keyword value string */
  636             fits_parse_value(card, value, NULL, status);
  637         
  638         if      (!strncmp(card+2, "ALGOR", 5) ) {
  639 
  640             /* set the desired compression algorithm */
  641                 /* allowed values: RICE_1, GZIP_1, GZIP_2, PLIO_1,     */
  642                 /*  HCOMPRESS_1, BZIP2_1, and NOCOMPRESS               */
  643 
  644                 if        (!fits_strncasecmp(value, "'RICE_1", 7) ) {
  645             comptype = RICE_1;
  646                 } else if (!fits_strncasecmp(value, "'GZIP_1", 7) ) {
  647             comptype = GZIP_1;
  648                 } else if (!fits_strncasecmp(value, "'GZIP_2", 7) ) {
  649             comptype = GZIP_2;
  650                 } else if (!fits_strncasecmp(value, "'PLIO_1", 7) ) {
  651             comptype = PLIO_1;
  652                 } else if (!fits_strncasecmp(value, "'HCOMPRESS_1", 12) ) {
  653             comptype = HCOMPRESS_1;
  654                 } else if (!fits_strncasecmp(value, "'NONE", 5) ) {
  655             comptype = NOCOMPRESS;
  656         } else {
  657             ffpmsg("Unknown FZALGOR keyword compression algorithm:");
  658             ffpmsg(value);
  659             return(*status = DATA_COMPRESSION_ERR);
  660         }  
  661 
  662             fits_set_compression_type (outfptr, comptype, status);
  663 
  664         } else if (!strncmp(card+2, "TILE  ", 6) ) {
  665 
  666                 if (!fits_strncasecmp(value, "'row", 4) ) {
  667                    tiledim[0] = -1;
  668         } else if (!fits_strncasecmp(value, "'whole", 6) ) {
  669                    tiledim[0] = -1;
  670                    tiledim[1] = -1;
  671                    tiledim[2] = -1;
  672                 } else {
  673            ffdtdm(infptr, value, 0,6, &naxis, tiledim, status);
  674                 }
  675 
  676             /* set the desired tile size */
  677         fits_set_tile_dim (outfptr, 6, tiledim, status);
  678 
  679         } else if (!strncmp(card+2, "QVALUE", 6) ) {
  680 
  681             /* set the desired Q quantization value */
  682         qvalue = atof(value);
  683         fits_set_quantize_level (outfptr, (float) qvalue, status);
  684 
  685         } else if (!strncmp(card+2, "QMETHD", 6) ) {
  686 
  687                     if (!fits_strncasecmp(value, "'no_dither", 10) ) {
  688                         ivalue = -1; /* just quantize, with no dithering */
  689             } else if (!fits_strncasecmp(value, "'subtractive_dither_1", 21) ) {
  690                         ivalue = SUBTRACTIVE_DITHER_1; /* use subtractive dithering */
  691             } else if (!fits_strncasecmp(value, "'subtractive_dither_2", 21) ) {
  692                         ivalue = SUBTRACTIVE_DITHER_2; /* dither, except preserve zero-valued pixels */
  693             } else {
  694                 ffpmsg("Unknown value for FZQUANT keyword: (set_compression_pref)");
  695             ffpmsg(value);
  696                         return(*status = DATA_COMPRESSION_ERR);
  697             }
  698 
  699             fits_set_quantize_method(outfptr, ivalue, status);
  700             
  701         } else if (!strncmp(card+2, "DTHRSD", 6) ) {
  702 
  703                 if (!fits_strncasecmp(value, "'checksum", 9) ) {
  704                     ivalue = -1; /* use checksum of first tile */
  705         } else if (!fits_strncasecmp(value, "'clock", 6) ) {
  706                     ivalue = 0; /* set dithering seed based on system clock */
  707         } else {  /* read integer value */
  708             if (*value == '\'')
  709                         ivalue = (int) atol(value+1); /* allow for leading quote character */
  710                     else 
  711                         ivalue = (int) atol(value); 
  712 
  713                     if (ivalue < 1 || ivalue > 10000) {
  714                 ffpmsg("Invalid value for FZDTHRSD keyword: (set_compression_pref)");
  715             ffpmsg(value);
  716                         return(*status = DATA_COMPRESSION_ERR);
  717                     }
  718         }
  719 
  720             /* set the desired dithering */
  721         fits_set_dither_seed(outfptr, ivalue, status);
  722 
  723         } else if (!strncmp(card+2, "I2F", 3) ) {
  724 
  725             /* set whether to convert integers to float then use lossy compression */
  726                 if (!fits_strcasecmp(value, "t") ) {
  727             fits_set_lossy_int (outfptr, 1, status);
  728         } else if (!fits_strcasecmp(value, "f") ) {
  729             fits_set_lossy_int (outfptr, 0, status);
  730         } else {
  731                 ffpmsg("Unknown value for FZI2F keyword: (set_compression_pref)");
  732             ffpmsg(value);
  733                         return(*status = DATA_COMPRESSION_ERR);
  734                 }
  735 
  736         } else if (!strncmp(card+2, "HSCALE ", 6) ) {
  737 
  738             /* set the desired Hcompress scale value */
  739         hscale = (float) atof(value);
  740         fits_set_hcomp_scale (outfptr, hscale, status);
  741             }
  742     }    
  743     }
  744     return(*status);
  745 }
  746 /*--------------------------------------------------------------------------*/
  747 int fits_get_noise_bits(fitsfile *fptr,  /* I - FITS file pointer   */
  748            int *noisebits,       /* noise_bits parameter value       */
  749                                 /* (default = 4)                    */
  750            int *status)         /* IO - error status                */
  751 {
  752 /*
  753    ********************************************************************
  754    ********************************************************************
  755    THIS ROUTINE IS PROVIDED ONLY FOR BACKWARDS COMPATIBILITY;
  756    ALL NEW SOFTWARE SHOULD CALL fits_set_quantize_level INSTEAD
  757    ********************************************************************
  758    ********************************************************************
  759 
  760 
  761    This routine returns the value of the noice_bits parameter that
  762    should be used when compressing floating point images.  The image is
  763    divided into tiles, and each tile is compressed and stored in a row
  764    of at variable length binary table column.
  765 
  766    Feb 2008: code changed to use the more general "quantize level" parameter
  767    rather than the "noise bits" parameter.  If quantize level is greater than
  768    zero, then the previous noisebits parameter is approximately given by
  769    
  770    noise bits = natural logarithm (quantize level) / natural log (2)
  771    
  772    This result is rounded to the nearest integer.
  773 */
  774     double qlevel;
  775 
  776     qlevel = (fptr->Fptr)->request_quantize_level;
  777 
  778     if (qlevel > 0. && qlevel < 65537. )
  779          *noisebits =  (int) ((log(qlevel) / log(2.0)) + 0.5);
  780     else 
  781         *noisebits = 0;
  782 
  783     return(*status);
  784 }
  785 /*--------------------------------------------------------------------------*/
  786 int fits_get_quantize_level(fitsfile *fptr,  /* I - FITS file pointer   */
  787            float *qlevel,       /* quantize level parameter value       */
  788            int *status)         /* IO - error status                */
  789 {
  790 /*
  791    This routine returns the value of the noice_bits parameter that
  792    should be used when compressing floating point images.  The image is
  793    divided into tiles, and each tile is compressed and stored in a row
  794    of at variable length binary table column.
  795 */
  796 
  797     if ((fptr->Fptr)->request_quantize_level == NO_QUANTIZE) {
  798       *qlevel = 0;
  799     } else {
  800       *qlevel = (fptr->Fptr)->request_quantize_level;
  801     }
  802 
  803     return(*status);
  804 }
  805 /*--------------------------------------------------------------------------*/
  806 int fits_get_dither_seed(fitsfile *fptr,  /* I - FITS file pointer   */
  807            int *offset,       /* dithering offset parameter value       */
  808            int *status)         /* IO - error status                */
  809 {
  810 /*
  811    This routine returns the value of the dithering offset parameter that
  812    is used when compressing floating point images.  The image is
  813    divided into tiles, and each tile is compressed and stored in a row
  814    of at variable length binary table column.
  815 */
  816 
  817     *offset = (fptr->Fptr)->request_dither_seed;
  818     return(*status);
  819 }/*--------------------------------------------------------------------------*/
  820 int fits_get_hcomp_scale(fitsfile *fptr,  /* I - FITS file pointer   */
  821            float *scale,          /* Hcompress scale parameter value       */
  822            int *status)         /* IO - error status                */
  823 
  824 {
  825 /*
  826    This routine returns the value of the noice_bits parameter that
  827    should be used when compressing floating point images.  The image is
  828    divided into tiles, and each tile is compressed and stored in a row
  829    of at variable length binary table column.
  830 */
  831 
  832     *scale = (fptr->Fptr)->request_hcomp_scale;
  833     return(*status);
  834 }
  835 /*--------------------------------------------------------------------------*/
  836 int fits_get_hcomp_smooth(fitsfile *fptr,  /* I - FITS file pointer   */
  837            int *smooth,          /* Hcompress smooth parameter value       */
  838            int *status)         /* IO - error status                */
  839 
  840 {
  841     *smooth = (fptr->Fptr)->request_hcomp_smooth;
  842     return(*status);
  843 }
  844 /*--------------------------------------------------------------------------*/
  845 int fits_img_compress(fitsfile *infptr, /* pointer to image to be compressed */
  846                  fitsfile *outfptr, /* empty HDU for output compressed image */
  847                  int *status)       /* IO - error status               */
  848 
  849 /*
  850    This routine initializes the output table, copies all the keywords,
  851    and  loops through the input image, compressing the data and
  852    writing the compressed tiles to the output table.
  853    
  854    This is a high level routine that is called by the fpack and funpack
  855    FITS compression utilities.
  856 */
  857 {
  858     int bitpix, naxis;
  859     long naxes[MAX_COMPRESS_DIM];
  860 /*    int c1, c2, c3; */
  861 
  862     if (*status > 0)
  863         return(*status);
  864 
  865 
  866     /* get datatype and size of input image */
  867     if (fits_get_img_param(infptr, MAX_COMPRESS_DIM, &bitpix, 
  868                        &naxis, naxes, status) > 0)
  869         return(*status);
  870 
  871     if (naxis < 1 || naxis > MAX_COMPRESS_DIM)
  872     {
  873         ffpmsg("Image cannot be compressed: NAXIS out of range");
  874         return(*status = BAD_NAXIS);
  875     }
  876 
  877     /* create a new empty HDU in the output file now, before setting the */
  878     /* compression preferences.  This HDU will become a binary table that */
  879     /* contains the compressed image.  If necessary, create a dummy primary */
  880     /* array, which much precede the binary table extension. */
  881     
  882     ffcrhd(outfptr, status);  /* this does nothing if the output file is empty */
  883 
  884     if ((outfptr->Fptr)->curhdu == 0)  /* have to create dummy primary array */
  885     {
  886        ffcrim(outfptr, 16, 0, NULL, status);
  887        ffcrhd(outfptr, status);
  888     } else {
  889         /* unset any compress parameter preferences that may have been
  890            set when closing the previous HDU in the output file */
  891         fits_unset_compression_param(outfptr, status);
  892     }
  893     
  894     /* set any compress parameter preferences as given in the input file */
  895     fits_set_compression_pref(infptr, outfptr, status);
  896 
  897     /* special case: the quantization level is not given by a keyword in  */
  898     /* the HDU header, so we have to explicitly copy the requested value */
  899     /* to the actual value */
  900 /* do this in imcomp_get_compressed_image_par, instead
  901     if ( (outfptr->Fptr)->request_quantize_level != 0.)
  902         (outfptr->Fptr)->quantize_level = (outfptr->Fptr)->request_quantize_level;
  903 */
  904     /* if requested, treat integer images same as a float image. */
  905     /* Then the pixels will be quantized (lossy algorithm) to achieve */
  906     /* higher amounts of compression than with lossless algorithms */
  907 
  908     if ( (outfptr->Fptr)->request_lossy_int_compress != 0  && bitpix > 0) 
  909     bitpix = FLOAT_IMG;  /* compress integer images as if float */
  910 
  911     /* initialize output table */
  912     if (imcomp_init_table(outfptr, bitpix, naxis, naxes, 0, status) > 0)
  913         return (*status);    
  914 
  915     /* Copy the image header keywords to the table header. */
  916     if (imcomp_copy_img2comp(infptr, outfptr, status) > 0)
  917         return (*status);
  918 
  919     /* turn off any intensity scaling (defined by BSCALE and BZERO */
  920     /* keywords) so that unscaled values will be read by CFITSIO */
  921     /* (except if quantizing an int image, same as a float image) */
  922     if ( (outfptr->Fptr)->request_lossy_int_compress == 0 && bitpix > 0) 
  923         ffpscl(infptr, 1.0, 0.0, status);
  924 
  925     /* force a rescan of the output file keywords, so that */
  926     /* the compression parameters will be copied to the internal */
  927     /* fitsfile structure used by CFITSIO */
  928     ffrdef(outfptr, status);
  929 
  930     /* turn off any intensity scaling (defined by BSCALE and BZERO */
  931     /* keywords) so that unscaled values will be written by CFITSIO */
  932     /* (except if quantizing an int image, same as a float image) */
  933     if ( (outfptr->Fptr)->request_lossy_int_compress == 0 && bitpix > 0) 
  934         ffpscl(outfptr, 1.0, 0.0, status);
  935 
  936     /* Read each image tile, compress, and write to a table row. */
  937     imcomp_compress_image (infptr, outfptr, status);
  938 
  939     /* force another rescan of the output file keywords, to */
  940     /* update PCOUNT and TFORMn = '1PB(iii)' keyword values. */
  941     ffrdef(outfptr, status);
  942 
  943     /* unset any previously set compress parameter preferences */
  944     fits_unset_compression_request(outfptr, status);
  945 
  946 /*
  947     fits_get_case(&c1, &c2, &c3);
  948     printf("c1, c2, c3 = %d, %d, %d\n", c1, c2, c3); 
  949 */
  950 
  951     return (*status);
  952 }
  953 /*--------------------------------------------------------------------------*/
  954 int imcomp_init_table(fitsfile *outfptr,
  955         int inbitpix,
  956         int naxis,
  957         long *naxes,
  958     int writebitpix,    /* write the ZBITPIX, ZNAXIS, and ZNAXES keyword? */
  959         int *status)
  960 /* 
  961   create a BINTABLE extension for the output compressed image.
  962 */
  963 {
  964     char keyname[FLEN_KEYWORD], zcmptype[12];
  965     int ii,  remain,  ndiv, addToDim, ncols, bitpix;
  966     long nrows;
  967     char *ttype[] = {"COMPRESSED_DATA", "ZSCALE", "ZZERO"};
  968     char *tform[3];
  969     char tf0[4], tf1[4], tf2[4];
  970     char *tunit[] = {"\0",            "\0",            "\0"  };
  971     char comm[FLEN_COMMENT];
  972     long actual_tilesize[MAX_COMPRESS_DIM]; /* Actual size to use for tiles */
  973     int is_primary=0; /* Is this attempting to write to the primary? */
  974     int nQualifyDims=0; /* For Hcompress, number of image dimensions with required pixels. */
  975     int noHigherDims=1; /* Set to true if all tile dims other than x are size 1. */
  976     int firstDim=-1, secondDim=-1; /* Indices of first and second tiles dimensions
  977                                 with width > 1 */
  978     
  979     if (*status > 0)
  980         return(*status);
  981 
  982     /* check for special case of losslessly compressing floating point */
  983     /* images.  Only compression algorithm that supports this is GZIP */
  984     if ( (inbitpix < 0) && ((outfptr->Fptr)->request_quantize_level == NO_QUANTIZE) ) {
  985        if (((outfptr->Fptr)->request_compress_type != GZIP_1) &&
  986            ((outfptr->Fptr)->request_compress_type != GZIP_2)) {
  987          ffpmsg("Lossless compression of floating point images must use GZIP (imcomp_init_table)");
  988          return(*status = DATA_COMPRESSION_ERR);
  989        }
  990     }
  991  
  992      /* set default compression parameter values, if undefined */
  993     
  994     if ( (outfptr->Fptr)->request_compress_type == 0) {
  995     /* use RICE_1 by default */
  996     (outfptr->Fptr)->request_compress_type = RICE_1;
  997     }
  998 
  999     if (inbitpix < 0 && (outfptr->Fptr)->request_quantize_level != NO_QUANTIZE) {  
 1000     /* set defaults for quantizing floating point images */
 1001     if ( (outfptr->Fptr)->request_quantize_method == 0) {
 1002           /* set default dithering method */
 1003               (outfptr->Fptr)->request_quantize_method = SUBTRACTIVE_DITHER_1;
 1004     }
 1005 
 1006     if ( (outfptr->Fptr)->request_quantize_level == 0) {
 1007         if ((outfptr->Fptr)->request_quantize_method == NO_DITHER) {
 1008             /* must use finer quantization if no dithering is done */
 1009             (outfptr->Fptr)->request_quantize_level = 16; 
 1010         } else {
 1011             (outfptr->Fptr)->request_quantize_level = 4; 
 1012         }
 1013         }
 1014     }
 1015 
 1016     /* special case: the quantization level is not given by a keyword in  */
 1017     /* the HDU header, so we have to explicitly copy the requested value */
 1018     /* to the actual value */
 1019 /* do this in imcomp_get_compressed_image_par, instead
 1020     if ( (outfptr->Fptr)->request_quantize_level != 0.)
 1021         (outfptr->Fptr)->quantize_level = (outfptr->Fptr)->request_quantize_level;
 1022 */
 1023     /* test for the 2 special cases that represent unsigned integers */
 1024     if (inbitpix == USHORT_IMG)
 1025         bitpix = SHORT_IMG;
 1026     else if (inbitpix == ULONG_IMG)
 1027         bitpix = LONG_IMG;
 1028     else if (inbitpix == SBYTE_IMG)
 1029         bitpix = BYTE_IMG;
 1030     else 
 1031         bitpix = inbitpix;
 1032 
 1033     /* reset default tile dimensions too if required */
 1034     memcpy(actual_tilesize, outfptr->Fptr->request_tilesize, MAX_COMPRESS_DIM * sizeof(long));
 1035 
 1036     if ((outfptr->Fptr)->request_compress_type == HCOMPRESS_1) {
 1037          
 1038          /* Tiles must ultimately have 2 (and only 2) dimensions, each with
 1039              at least 4 pixels. First catch the case where the image
 1040              itself won't allow this. */
 1041          if (naxis < 2 ) {
 1042             ffpmsg("Hcompress cannot be used with 1-dimensional images (imcomp_init_table)");
 1043             return(*status = DATA_COMPRESSION_ERR);
 1044      }
 1045          for (ii=0; ii<naxis; ii++)
 1046          {
 1047             if (naxes[ii] >= 4)
 1048                ++nQualifyDims;
 1049          }
 1050          if (nQualifyDims < 2)
 1051          {
 1052             ffpmsg("Hcompress minimum image dimension is 4 pixels (imcomp_init_table)");
 1053             return(*status = DATA_COMPRESSION_ERR);            
 1054          }
 1055 
 1056          /* Handle 2 special cases for backwards compatibility.
 1057             1) If both X and Y tile dims are set to full size, ignore
 1058                any other requested dimensions and just set their sizes to 1. 
 1059             2) If X is full size and all the rest are size 1, attempt to
 1060                find a reasonable size for Y. All other 1-D tile specifications
 1061                will be rejected. */
 1062          for (ii=1; ii<naxis; ++ii)
 1063             if (actual_tilesize[ii] != 0 && actual_tilesize[ii] != 1)
 1064             {
 1065                noHigherDims = 0;
 1066                break;
 1067             }
 1068 
 1069          if ((actual_tilesize[0] <= 0) &&
 1070              (actual_tilesize[1] == -1) ){
 1071          
 1072         /* compress the whole image as a single tile */
 1073              actual_tilesize[0] = naxes[0];
 1074              actual_tilesize[1] = naxes[1];
 1075 
 1076               for (ii = 2; ii < naxis; ii++) {
 1077                  /* set all higher tile dimensions = 1 */
 1078                      actual_tilesize[ii] = 1;
 1079           }
 1080 
 1081          } else if ((actual_tilesize[0] <= 0) && noHigherDims) {
 1082          
 1083              /*
 1084               The Hcompress algorithm is inherently 2D in nature, so the row by row
 1085           tiling that is used for other compression algorithms is not appropriate.
 1086           If the image has less than 30 rows, then the entire image will be compressed
 1087           as a single tile.  Otherwise the tiles will consist of 16 rows of the image. 
 1088           This keeps the tiles to a reasonable size, and it also includes enough rows
 1089           to allow good compression efficiency.  If the last tile of the image 
 1090           happens to contain less than 4 rows, then find another tile size with
 1091           between 14 and 30 rows (preferably even), so that the last tile has 
 1092           at least 4 rows
 1093          */ 
 1094           
 1095              /* 1st tile dimension is the row length of the image */
 1096              actual_tilesize[0] = naxes[0];
 1097 
 1098               if (naxes[1] <= 30) {  /* use whole image if it is small */
 1099                    actual_tilesize[1] = naxes[1];
 1100           } else {
 1101                 /* look for another good tile dimension */
 1102               if        (naxes[1] % 16 == 0 || naxes[1] % 16 > 3) {
 1103                       actual_tilesize[1] = 16;
 1104           } else if (naxes[1] % 24 == 0 || naxes[1] % 24 > 3) {
 1105                       actual_tilesize[1] = 24;
 1106           } else if (naxes[1] % 20 == 0 || naxes[1] % 20 > 3) {
 1107                       actual_tilesize[1] = 20;
 1108           } else if (naxes[1] % 30 == 0 || naxes[1] % 30 > 3) {
 1109                       actual_tilesize[1] = 30;
 1110           } else if (naxes[1] % 28 == 0 || naxes[1] % 28 > 3) {
 1111                       actual_tilesize[1] = 28;
 1112           } else if (naxes[1] % 26 == 0 || naxes[1] % 26 > 3) {
 1113                       actual_tilesize[1] = 26;
 1114           } else if (naxes[1] % 22 == 0 || naxes[1] % 22 > 3) {
 1115                       actual_tilesize[1] = 22;
 1116           } else if (naxes[1] % 18 == 0 || naxes[1] % 18 > 3) {
 1117                       actual_tilesize[1] = 18;
 1118           } else if (naxes[1] % 14 == 0 || naxes[1] % 14 > 3) {
 1119                       actual_tilesize[1] = 14;
 1120           } else  {
 1121                       actual_tilesize[1] = 17;
 1122           }
 1123           }
 1124         } else {
 1125            if (actual_tilesize[0] <= 0)
 1126               actual_tilesize[0] = naxes[0];
 1127            for (ii=1; ii<naxis; ++ii)
 1128            {
 1129               if (actual_tilesize[ii] < 0)
 1130                  actual_tilesize[ii] = naxes[ii];
 1131               else if (actual_tilesize[ii] == 0)
 1132                  actual_tilesize[ii] = 1;
 1133            }
 1134         }
 1135         
 1136         for (ii=0; ii<naxis; ++ii)
 1137         {
 1138            if (actual_tilesize[ii] > 1)
 1139            {
 1140               if (firstDim < 0)
 1141                  firstDim = ii;
 1142               else if (secondDim < 0)
 1143                  secondDim = ii;
 1144               else
 1145               {
 1146                  ffpmsg("Hcompress tiles can only have 2 dimensions (imcomp_init_table)");
 1147                  return(*status = DATA_COMPRESSION_ERR);
 1148               }
 1149            }
 1150         }
 1151         if (firstDim < 0 || secondDim < 0)
 1152         {
 1153             ffpmsg("Hcompress tiles must have 2 dimensions (imcomp_init_table)");
 1154             return(*status = DATA_COMPRESSION_ERR);
 1155         }
 1156         
 1157         if (actual_tilesize[firstDim] < 4 || actual_tilesize[secondDim] < 4)
 1158         {
 1159            ffpmsg("Hcompress minimum tile dimension is 4 pixels (imcomp_init_table)");
 1160            return (*status = DATA_COMPRESSION_ERR);
 1161         }
 1162     
 1163         /* check if requested tile size causes the last tile to to have less than 4 pixels */
 1164         remain = naxes[firstDim] % (actual_tilesize[firstDim]);  /* 1st dimension */
 1165         if (remain > 0 && remain < 4) {
 1166             ndiv = naxes[firstDim]/actual_tilesize[firstDim]; /* integer truncation is intentional */
 1167             addToDim = ceil((double)remain/ndiv);
 1168             (actual_tilesize[firstDim]) += addToDim; /* increase tile size */
 1169        
 1170             remain = naxes[firstDim] % (actual_tilesize[firstDim]);
 1171             if (remain > 0 && remain < 4) {
 1172                 ffpmsg("Last tile along 1st dimension has less than 4 pixels (imcomp_init_table)");
 1173                 return(*status = DATA_COMPRESSION_ERR); 
 1174             }        
 1175         }
 1176 
 1177         remain = naxes[secondDim] % (actual_tilesize[secondDim]);  /* 2nd dimension */
 1178         if (remain > 0 && remain < 4) {
 1179             ndiv = naxes[secondDim]/actual_tilesize[secondDim]; /* integer truncation is intentional */
 1180             addToDim = ceil((double)remain/ndiv);
 1181             (actual_tilesize[secondDim]) += addToDim; /* increase tile size */
 1182        
 1183             remain = naxes[secondDim] % (actual_tilesize[secondDim]);
 1184             if (remain > 0 && remain < 4) {
 1185                 ffpmsg("Last tile along 2nd dimension has less than 4 pixels (imcomp_init_table)");
 1186                 return(*status = DATA_COMPRESSION_ERR); 
 1187             }        
 1188         }
 1189 
 1190     } /* end, if HCOMPRESS_1 */
 1191     
 1192     for (ii = 0; ii < naxis; ii++) {
 1193     if (ii == 0) { /* first axis is different */
 1194         if (actual_tilesize[ii] <= 0) {
 1195                 actual_tilesize[ii] = naxes[ii]; 
 1196         }
 1197     } else {
 1198         if (actual_tilesize[ii] < 0) {
 1199                 actual_tilesize[ii] = naxes[ii];  /* negative value maean use whole length */
 1200         } else if (actual_tilesize[ii] == 0) {
 1201                 actual_tilesize[ii] = 1;  /* zero value means use default value = 1 */
 1202         }
 1203     }
 1204     }
 1205 
 1206     /* ---- set up array of TFORM strings -------------------------------*/
 1207     if ( (outfptr->Fptr)->request_huge_hdu != 0) {
 1208         strcpy(tf0, "1QB");
 1209     } else {
 1210         strcpy(tf0, "1PB");
 1211     }
 1212     strcpy(tf1, "1D");
 1213     strcpy(tf2, "1D");
 1214 
 1215     tform[0] = tf0;
 1216     tform[1] = tf1;
 1217     tform[2] = tf2;
 1218 
 1219     /* calculate number of rows in output table */
 1220     nrows = 1;
 1221     for (ii = 0; ii < naxis; ii++)
 1222     {
 1223         nrows = nrows * ((naxes[ii] - 1)/ (actual_tilesize[ii]) + 1);
 1224     }
 1225 
 1226     /* determine the default  number of columns in the output table */
 1227     if (bitpix < 0 && (outfptr->Fptr)->request_quantize_level != NO_QUANTIZE)  
 1228         ncols = 3;  /* quantized and scaled floating point image */
 1229     else
 1230         ncols = 1; /* default table has just one 'COMPRESSED_DATA' column */
 1231 
 1232     if ((outfptr->Fptr)->request_compress_type == RICE_1)
 1233     {
 1234         strcpy(zcmptype, "RICE_1");
 1235     }
 1236     else if ((outfptr->Fptr)->request_compress_type == GZIP_1)
 1237     {
 1238         strcpy(zcmptype, "GZIP_1");
 1239     }
 1240     else if ((outfptr->Fptr)->request_compress_type == GZIP_2)
 1241     {
 1242         strcpy(zcmptype, "GZIP_2");
 1243     }
 1244     else if ((outfptr->Fptr)->request_compress_type == BZIP2_1)
 1245     {
 1246         strcpy(zcmptype, "BZIP2_1");
 1247     }
 1248     else if ((outfptr->Fptr)->request_compress_type == PLIO_1)
 1249     {
 1250         strcpy(zcmptype, "PLIO_1");
 1251        /* the PLIO compression algorithm outputs short integers, not bytes */
 1252         if ( (outfptr->Fptr)->request_huge_hdu != 0) {
 1253             strcpy(tform[0], "1QI");
 1254         } else {
 1255             strcpy(tform[0], "1PI");
 1256         }
 1257     }
 1258     else if ((outfptr->Fptr)->request_compress_type == HCOMPRESS_1)
 1259     {
 1260         strcpy(zcmptype, "HCOMPRESS_1");
 1261     }
 1262     else if ((outfptr->Fptr)->request_compress_type == NOCOMPRESS)
 1263     {
 1264         strcpy(zcmptype, "NOCOMPRESS");
 1265     }    
 1266     else
 1267     {
 1268         ffpmsg("unknown compression type (imcomp_init_table)");
 1269         return(*status = DATA_COMPRESSION_ERR);
 1270     }
 1271 
 1272     /* If attempting to write compressed image to primary, the
 1273        call to ffcrtb will increment Fptr->curhdu to 1.  Therefore
 1274        we need to test now for setting is_primary */
 1275     is_primary = (outfptr->Fptr->curhdu == 0);
 1276     /* create the bintable extension to contain the compressed image */
 1277     ffcrtb(outfptr, BINARY_TBL, nrows, ncols, ttype, 
 1278                 tform, tunit, 0, status);
 1279 
 1280     /* Add standard header keywords. */
 1281     ffpkyl (outfptr, "ZIMAGE", 1, 
 1282            "extension contains compressed image", status);  
 1283 
 1284     if (writebitpix) {
 1285         /*  write the keywords defining the datatype and dimensions of */
 1286     /*  the uncompressed image.  If not, these keywords will be */
 1287         /*  copied later from the input uncompressed image  */
 1288     
 1289         if (is_primary)   
 1290             ffpkyl (outfptr, "ZSIMPLE", 1,
 1291             "file does conform to FITS standard", status);
 1292         ffpkyj (outfptr, "ZBITPIX", bitpix,
 1293             "data type of original image", status);
 1294         ffpkyj (outfptr, "ZNAXIS", naxis,
 1295             "dimension of original image", status);
 1296 
 1297         for (ii = 0;  ii < naxis;  ii++)
 1298         {
 1299             snprintf (keyname, FLEN_KEYWORD,"ZNAXIS%d", ii+1);
 1300             ffpkyj (outfptr, keyname, naxes[ii],
 1301             "length of original image axis", status);
 1302         }
 1303     }
 1304                       
 1305     for (ii = 0;  ii < naxis;  ii++)
 1306     {
 1307         snprintf (keyname, FLEN_KEYWORD,"ZTILE%d", ii+1);
 1308         ffpkyj (outfptr, keyname, actual_tilesize[ii],
 1309             "size of tiles to be compressed", status);
 1310     }
 1311 
 1312     if (bitpix < 0) {
 1313        
 1314     if ((outfptr->Fptr)->request_quantize_level == NO_QUANTIZE) {
 1315         ffpkys(outfptr, "ZQUANTIZ", "NONE", 
 1316           "Lossless compression without quantization", status);
 1317     } else {
 1318         
 1319         /* Unless dithering has been specifically turned off by setting */
 1320         /* request_quantize_method = -1, use dithering by default */
 1321         /* when quantizing floating point images. */
 1322     
 1323         if ( (outfptr->Fptr)->request_quantize_method == 0) 
 1324               (outfptr->Fptr)->request_quantize_method = SUBTRACTIVE_DITHER_1;
 1325        
 1326         if ((outfptr->Fptr)->request_quantize_method == SUBTRACTIVE_DITHER_1) {
 1327           ffpkys(outfptr, "ZQUANTIZ", "SUBTRACTIVE_DITHER_1", 
 1328             "Pixel Quantization Algorithm", status);
 1329 
 1330           /* also write the associated ZDITHER0 keyword with a default value */
 1331           /* which may get updated later. */
 1332               ffpky(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->request_dither_seed), 
 1333            "dithering offset when quantizing floats", status);
 1334  
 1335             } else if ((outfptr->Fptr)->request_quantize_method == SUBTRACTIVE_DITHER_2) {
 1336           ffpkys(outfptr, "ZQUANTIZ", "SUBTRACTIVE_DITHER_2", 
 1337             "Pixel Quantization Algorithm", status);
 1338 
 1339           /* also write the associated ZDITHER0 keyword with a default value */
 1340           /* which may get updated later. */
 1341               ffpky(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->request_dither_seed), 
 1342            "dithering offset when quantizing floats", status);
 1343 
 1344           if (!strcmp(zcmptype, "RICE_1"))  {
 1345             /* when using this new dithering method, change the compression type */
 1346         /* to an alias, so that old versions of funpack will not be able to */
 1347         /* created a corrupted uncompressed image. */
 1348         /* ******* can remove this cludge after about June 2015, after most old versions of fpack are gone */
 1349             strcpy(zcmptype, "RICE_ONE");
 1350           }
 1351 
 1352             } else if ((outfptr->Fptr)->request_quantize_method == NO_DITHER) {
 1353           ffpkys(outfptr, "ZQUANTIZ", "NO_DITHER", 
 1354             "No dithering during quantization", status);
 1355         }
 1356 
 1357     }
 1358     }
 1359 
 1360     ffpkys (outfptr, "ZCMPTYPE", zcmptype,
 1361               "compression algorithm", status);
 1362 
 1363     /* write any algorithm-specific keywords */
 1364     if ((outfptr->Fptr)->request_compress_type == RICE_1)
 1365     {
 1366         ffpkys (outfptr, "ZNAME1", "BLOCKSIZE",
 1367             "compression block size", status);
 1368 
 1369         /* for now at least, the block size is always 32 */
 1370         ffpkyj (outfptr, "ZVAL1", 32,
 1371             "pixels per block", status);
 1372 
 1373         ffpkys (outfptr, "ZNAME2", "BYTEPIX",
 1374             "bytes per pixel (1, 2, 4, or 8)", status);
 1375 
 1376         if (bitpix == BYTE_IMG)
 1377             ffpkyj (outfptr, "ZVAL2", 1,
 1378             "bytes per pixel (1, 2, 4, or 8)", status);
 1379         else if (bitpix == SHORT_IMG)
 1380             ffpkyj (outfptr, "ZVAL2", 2,
 1381             "bytes per pixel (1, 2, 4, or 8)", status);
 1382         else 
 1383             ffpkyj (outfptr, "ZVAL2", 4,
 1384             "bytes per pixel (1, 2, 4, or 8)", status);
 1385 
 1386     }
 1387     else if ((outfptr->Fptr)->request_compress_type == HCOMPRESS_1)
 1388     {
 1389         ffpkys (outfptr, "ZNAME1", "SCALE",
 1390             "HCOMPRESS scale factor", status);
 1391         ffpkye (outfptr, "ZVAL1", (outfptr->Fptr)->request_hcomp_scale,
 1392         7, "HCOMPRESS scale factor", status);
 1393 
 1394         ffpkys (outfptr, "ZNAME2", "SMOOTH",
 1395             "HCOMPRESS smooth option", status);
 1396         ffpkyj (outfptr, "ZVAL2", (long) (outfptr->Fptr)->request_hcomp_smooth,
 1397             "HCOMPRESS smooth option", status);
 1398     }
 1399 
 1400     /* Write the BSCALE and BZERO keywords, if an unsigned integer image */
 1401     if (inbitpix == USHORT_IMG)
 1402     {
 1403         strcpy(comm, "offset data range to that of unsigned short");
 1404         ffpkyg(outfptr, "BZERO", 32768., 0, comm, status);
 1405         strcpy(comm, "default scaling factor");
 1406         ffpkyg(outfptr, "BSCALE", 1.0, 0, comm, status);
 1407     }
 1408     else if (inbitpix == SBYTE_IMG)
 1409     {
 1410         strcpy(comm, "offset data range to that of signed byte");
 1411         ffpkyg(outfptr, "BZERO", -128., 0, comm, status);
 1412         strcpy(comm, "default scaling factor");
 1413         ffpkyg(outfptr, "BSCALE", 1.0, 0, comm, status);
 1414     }
 1415     else if (inbitpix == ULONG_IMG)
 1416     {
 1417         strcpy(comm, "offset data range to that of unsigned long");
 1418         ffpkyg(outfptr, "BZERO", 2147483648., 0, comm, status);
 1419         strcpy(comm, "default scaling factor");
 1420         ffpkyg(outfptr, "BSCALE", 1.0, 0, comm, status);
 1421     }
 1422 
 1423     return(*status);
 1424 }
 1425 /*--------------------------------------------------------------------------*/
 1426 int imcomp_calc_max_elem (int comptype, int nx, int zbitpix, int blocksize)
 1427 
 1428 /* This function returns the maximum number of bytes in a compressed
 1429    image line.
 1430 
 1431     nx = maximum number of pixels in a tile
 1432     blocksize is only relevant for RICE compression
 1433 */
 1434 {    
 1435     if (comptype == RICE_1)
 1436     {
 1437         if (zbitpix == 16)
 1438             return (sizeof(short) * nx + nx / blocksize + 2 + 4);
 1439     else
 1440             return (sizeof(float) * nx + nx / blocksize + 2 + 4);
 1441     }
 1442     else if ((comptype == GZIP_1) || (comptype == GZIP_2))
 1443     {
 1444         /* gzip usually compressed by at least a factor of 2 for I*4 images */
 1445         /* and somewhat less for I*2 images */
 1446         /* If this size turns out to be too small, then the gzip */
 1447         /* compression routine will allocate more space as required */
 1448         /* to be on the safe size, allocate buffer same size as input */
 1449     
 1450         if (zbitpix == 16)
 1451             return(nx * 2);
 1452     else if (zbitpix == 8)
 1453             return(nx);
 1454     else
 1455             return(nx * 4);
 1456     }
 1457     else if (comptype == BZIP2_1)
 1458     {
 1459         /* To guarantee that the compressed data will fit, allocate an output
 1460        buffer of size 1% larger than the uncompressed data, plus 600 bytes */
 1461 
 1462             return((int) (nx * 1.01 * zbitpix / 8. + 601.));
 1463     }
 1464      else if (comptype == HCOMPRESS_1)
 1465     {
 1466         /* Imperical evidence suggests in the worst case, 
 1467        the compressed stream could be up to 10% larger than the original
 1468        image.  Add 26 byte overhead, only significant for very small tiles
 1469        
 1470          Possible improvement: may need to allow a larger size for 32-bit images */
 1471 
 1472         if (zbitpix == 16 || zbitpix == 8)
 1473     
 1474             return( (int) (nx * 2.2 + 26));   /* will be compressing 16-bit int array */
 1475         else
 1476             return( (int) (nx * 4.4 + 26));   /* will be compressing 32-bit int array */
 1477     }
 1478     else
 1479         return(nx * sizeof(int));
 1480 }
 1481 /*--------------------------------------------------------------------------*/
 1482 int imcomp_compress_image (fitsfile *infptr, fitsfile *outfptr, int *status)
 1483 
 1484 /* This routine does the following:
 1485         - reads an image one tile at a time
 1486         - if it is a float or double image, then it tries to quantize the pixels
 1487           into scaled integers.
 1488         - it then compressess the integer pixels, or if the it was not
 1489       possible to quantize the floating point pixels, then it losslessly
 1490       compresses them with gzip
 1491     - writes the compressed byte stream to the output FITS file
 1492 */
 1493 {
 1494     double *tiledata;
 1495     int anynul, gotnulls = 0, datatype;
 1496     long ii, row;
 1497     int naxis;
 1498     double dummy = 0., dblnull = DOUBLENULLVALUE;
 1499     float fltnull = FLOATNULLVALUE;
 1500     long maxtilelen, tilelen, incre[] = {1, 1, 1, 1, 1, 1};
 1501     long naxes[MAX_COMPRESS_DIM], fpixel[MAX_COMPRESS_DIM];
 1502     long lpixel[MAX_COMPRESS_DIM], tile[MAX_COMPRESS_DIM];
 1503     long tilesize[MAX_COMPRESS_DIM];
 1504     long i0, i1, i2, i3, i4, i5, trowsize, ntrows;
 1505     char card[FLEN_CARD];
 1506 
 1507     if (*status > 0)
 1508         return(*status);
 1509 
 1510     maxtilelen = (outfptr->Fptr)->maxtilelen;
 1511 
 1512     /* 
 1513      Allocate buffer to hold 1 tile of data; size depends on which compression 
 1514      algorithm is used:
 1515 
 1516       Rice and GZIP will compress byte, short, or int arrays without conversion.
 1517       PLIO requires 4-byte int values, so byte and short arrays must be converted to int.
 1518       HCompress internally converts byte or short values to ints, and
 1519          converts int values to 8-byte longlong integers.
 1520     */
 1521     
 1522     if ((outfptr->Fptr)->zbitpix == FLOAT_IMG)
 1523     {
 1524         datatype = TFLOAT;
 1525 
 1526         if ( (outfptr->Fptr)->compress_type == HCOMPRESS_1) {
 1527         /* need twice as much scratch space (8 bytes per pixel) */
 1528             tiledata = (double*) malloc (maxtilelen * 2 *sizeof (float));   
 1529     } else {
 1530             tiledata = (double*) malloc (maxtilelen * sizeof (float));
 1531     }
 1532     }
 1533     else if ((outfptr->Fptr)->zbitpix == DOUBLE_IMG)
 1534     {
 1535         datatype = TDOUBLE;
 1536         tiledata = (double*) malloc (maxtilelen * sizeof (double));
 1537     }
 1538     else if ((outfptr->Fptr)->zbitpix == SHORT_IMG)
 1539     {
 1540         datatype = TSHORT;
 1541         if ( (outfptr->Fptr)->compress_type == RICE_1  ||
 1542          (outfptr->Fptr)->compress_type == GZIP_1  ||
 1543          (outfptr->Fptr)->compress_type == GZIP_2  ||
 1544          (outfptr->Fptr)->compress_type == BZIP2_1 ||
 1545              (outfptr->Fptr)->compress_type == NOCOMPRESS) {
 1546         /* only need  buffer of I*2 pixels for gzip, bzip2, and Rice */
 1547 
 1548             tiledata = (double*) malloc (maxtilelen * sizeof (short));  
 1549     } else {
 1550         /*  need  buffer of I*4 pixels for Hcompress and PLIO */
 1551             tiledata = (double*) malloc (maxtilelen * sizeof (int));
 1552         }
 1553     }
 1554     else if ((outfptr->Fptr)->zbitpix == BYTE_IMG)
 1555     {
 1556 
 1557         datatype = TBYTE;
 1558         if ( (outfptr->Fptr)->compress_type == RICE_1  ||
 1559          (outfptr->Fptr)->compress_type == BZIP2_1 ||
 1560          (outfptr->Fptr)->compress_type == GZIP_1  ||
 1561          (outfptr->Fptr)->compress_type == GZIP_2) {
 1562         /* only need  buffer of I*1 pixels for gzip, bzip2, and Rice */
 1563 
 1564             tiledata = (double*) malloc (maxtilelen);   
 1565     } else {
 1566         /*  need  buffer of I*4 pixels for Hcompress and PLIO */
 1567             tiledata = (double*) malloc (maxtilelen * sizeof (int));
 1568         }
 1569     }
 1570     else if ((outfptr->Fptr)->zbitpix == LONG_IMG)
 1571     {
 1572         datatype = TINT;
 1573         if ( (outfptr->Fptr)->compress_type == HCOMPRESS_1) {
 1574         /* need twice as much scratch space (8 bytes per pixel) */
 1575 
 1576             tiledata = (double*) malloc (maxtilelen * 2 * sizeof (int));    
 1577     } else {
 1578         /* only need  buffer of I*4 pixels for gzip, bzip2,  Rice, and PLIO */
 1579 
 1580             tiledata = (double*) malloc (maxtilelen * sizeof (int));
 1581         }
 1582     }
 1583     else
 1584     {
 1585     ffpmsg("Bad image datatype. (imcomp_compress_image)");
 1586     return (*status = MEMORY_ALLOCATION);
 1587     }
 1588     
 1589     if (tiledata == NULL)
 1590     {
 1591     ffpmsg("Out of memory. (imcomp_compress_image)");
 1592     return (*status = MEMORY_ALLOCATION);
 1593     }
 1594 
 1595     /*  calculate size of tile in each dimension */
 1596     naxis = (outfptr->Fptr)->zndim;
 1597     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 1598     {
 1599         if (ii < naxis)
 1600         {
 1601              naxes[ii] = (outfptr->Fptr)->znaxis[ii];
 1602              tilesize[ii] = (outfptr->Fptr)->tilesize[ii];
 1603         }
 1604         else
 1605         {
 1606             naxes[ii] = 1;
 1607             tilesize[ii] = 1;
 1608         }
 1609     }
 1610     row = 1;
 1611 
 1612     /* set up big loop over up to 6 dimensions */
 1613     for (i5 = 1; i5 <= naxes[5]; i5 += tilesize[5])
 1614     {
 1615      fpixel[5] = i5;
 1616      lpixel[5] = minvalue(i5 + tilesize[5] - 1, naxes[5]);
 1617      tile[5] = lpixel[5] - fpixel[5] + 1;
 1618      for (i4 = 1; i4 <= naxes[4]; i4 += tilesize[4])
 1619      {
 1620       fpixel[4] = i4;
 1621       lpixel[4] = minvalue(i4 + tilesize[4] - 1, naxes[4]);
 1622       tile[4] = lpixel[4] - fpixel[4] + 1;
 1623       for (i3 = 1; i3 <= naxes[3]; i3 += tilesize[3])
 1624       {
 1625        fpixel[3] = i3;
 1626        lpixel[3] = minvalue(i3 + tilesize[3] - 1, naxes[3]);
 1627        tile[3] = lpixel[3] - fpixel[3] + 1;
 1628        for (i2 = 1; i2 <= naxes[2]; i2 += tilesize[2])
 1629        {
 1630         fpixel[2] = i2;
 1631         lpixel[2] = minvalue(i2 + tilesize[2] - 1, naxes[2]);
 1632         tile[2] = lpixel[2] - fpixel[2] + 1;
 1633         for (i1 = 1; i1 <= naxes[1]; i1 += tilesize[1])
 1634         {
 1635          fpixel[1] = i1;
 1636          lpixel[1] = minvalue(i1 + tilesize[1] - 1, naxes[1]);
 1637          tile[1] = lpixel[1] - fpixel[1] + 1;
 1638          for (i0 = 1; i0 <= naxes[0]; i0 += tilesize[0])
 1639          {
 1640           fpixel[0] = i0;
 1641           lpixel[0] = minvalue(i0 + tilesize[0] - 1, naxes[0]);
 1642           tile[0] = lpixel[0] - fpixel[0] + 1;
 1643 
 1644           /* number of pixels in this tile */
 1645           tilelen = tile[0];
 1646           for (ii = 1; ii < naxis; ii++)
 1647           {
 1648              tilelen *= tile[ii];
 1649           }
 1650 
 1651           /* read next tile of data from image */
 1652       anynul = 0;
 1653           if (datatype == TFLOAT)
 1654           {
 1655               ffgsve(infptr, 1, naxis, naxes, fpixel, lpixel, incre, 
 1656                   FLOATNULLVALUE, (float *) tiledata,  &anynul, status);
 1657           }
 1658           else if (datatype == TDOUBLE)
 1659           {
 1660               ffgsvd(infptr, 1, naxis, naxes, fpixel, lpixel, incre, 
 1661                   DOUBLENULLVALUE, tiledata, &anynul, status);
 1662           }
 1663           else if (datatype == TINT)
 1664           {
 1665               ffgsvk(infptr, 1, naxis, naxes, fpixel, lpixel, incre, 
 1666                   0, (int *) tiledata,  &anynul, status);
 1667           }
 1668           else if (datatype == TSHORT)
 1669           {
 1670               ffgsvi(infptr, 1, naxis, naxes, fpixel, lpixel, incre, 
 1671                   0, (short *) tiledata,  &anynul, status);
 1672           }
 1673           else if (datatype == TBYTE)
 1674           {
 1675               ffgsvb(infptr, 1, naxis, naxes, fpixel, lpixel, incre, 
 1676                   0, (unsigned char *) tiledata,  &anynul, status);
 1677           }
 1678           else 
 1679           {
 1680               ffpmsg("Error bad datatype of image tile to compress");
 1681               free(tiledata);
 1682               return (*status);
 1683           }
 1684 
 1685           /* now compress the tile, and write to row of binary table */
 1686           /*   NOTE: we don't have to worry about the presence of null values in the
 1687            array if it is an integer array:  the null value is simply encoded
 1688            in the compressed array just like any other pixel value.  
 1689            
 1690            If it is a floating point array, then we need to check for null
 1691            only if the anynul parameter returned a true value when reading the tile
 1692       */
 1693           
 1694           /* Collapse sizes of higher dimension tiles into 2 dimensional
 1695              equivalents needed by the quantizing algorithms for
 1696              floating point types */
 1697           fits_calc_tile_rows(lpixel, fpixel, naxis, &trowsize,
 1698                             &ntrows, status);
 1699 
 1700           if (anynul && datatype == TFLOAT) {
 1701               imcomp_compress_tile(outfptr, row, datatype, tiledata, tilelen,
 1702                                trowsize, ntrows, 1, &fltnull, status);
 1703           } else if (anynul && datatype == TDOUBLE) {
 1704               imcomp_compress_tile(outfptr, row, datatype, tiledata, tilelen,
 1705                                trowsize, ntrows, 1, &dblnull, status);
 1706           } else {
 1707               imcomp_compress_tile(outfptr, row, datatype, tiledata, tilelen,
 1708                                trowsize, ntrows, 0, &dummy, status);
 1709           }
 1710 
 1711           /* set flag if we found any null values */
 1712           if (anynul)
 1713               gotnulls = 1;
 1714 
 1715           /* check for any error in the previous operations */
 1716           if (*status > 0)
 1717           {
 1718               ffpmsg("Error writing compressed image to table");
 1719               free(tiledata);
 1720               return (*status);
 1721           }
 1722 
 1723       row++;
 1724          }
 1725         }
 1726        }
 1727       }
 1728      }
 1729     }
 1730 
 1731     free (tiledata);  /* finished with this buffer */
 1732 
 1733     /* insert ZBLANK keyword if necessary; only for TFLOAT or TDOUBLE images */
 1734     if (gotnulls)
 1735     {
 1736           ffgcrd(outfptr, "ZCMPTYPE", card, status);
 1737           ffikyj(outfptr, "ZBLANK", COMPRESS_NULL_VALUE, 
 1738              "null value in the compressed integer array", status);
 1739     }
 1740 
 1741     return (*status);
 1742 }
 1743 /*--------------------------------------------------------------------------*/
 1744 int imcomp_compress_tile (fitsfile *outfptr,
 1745     long row,  /* tile number = row in the binary table that holds the compressed data */
 1746     int datatype, 
 1747     void *tiledata, 
 1748     long tilelen,
 1749     long tilenx,
 1750     long tileny,
 1751     int nullcheck,
 1752     void *nullflagval,
 1753     int *status)
 1754 
 1755 /*
 1756    This is the main compression routine.
 1757 
 1758    This routine does the following to the input tile of pixels:
 1759         - if it is a float or double image, then it quantizes the pixels
 1760         - compresses the integer pixel values
 1761         - writes the compressed byte stream to the FITS file.
 1762 
 1763    If the tile cannot be quantized than the raw float or double values
 1764    are losslessly compressed with gzip and then written to the output table.
 1765    
 1766    This input array may be modified by this routine.  If the array is of type TINT
 1767    or TFLOAT, and the compression type is HCOMPRESS, then it must have been 
 1768    allocated to be twice as large (8 bytes per pixel) to provide scratch space.
 1769 
 1770   Note that this routine does not fully support the implicit datatype conversion that
 1771   is supported when writing to normal FITS images.  The datatype of the input array
 1772   must have the same datatype (either signed or unsigned) as the output (compressed)
 1773   FITS image in some cases.
 1774 */
 1775 {
 1776     int *idata;     /* quantized integer data */
 1777     int cn_zblank, zbitpix, nullval;
 1778     int flag = 1;  /* true by default; only = 0 if float data couldn't be quantized */
 1779     int intlength;      /* size of integers to be compressed */
 1780     double scale, zero, actual_bzero;
 1781     long ii;
 1782     size_t clen;        /* size of cbuf */
 1783     short *cbuf;    /* compressed data */
 1784     int  nelem = 0;     /* number of bytes */
 1785     int tilecol;
 1786     size_t gzip_nelem = 0;
 1787     unsigned int bzlen;
 1788     int ihcompscale;
 1789     float hcompscale;
 1790     double noise2, noise3, noise5;
 1791     double bscale[1] = {1.}, bzero[1] = {0.};   /* scaling parameters */
 1792     long  hcomp_len;
 1793     LONGLONG *lldata;
 1794 
 1795     if (*status > 0)
 1796         return(*status);
 1797 
 1798     /* check for special case of losslessly compressing floating point */
 1799     /* images.  Only compression algorithm that supports this is GZIP */
 1800     if ( (outfptr->Fptr)->quantize_level == NO_QUANTIZE) {
 1801        if (((outfptr->Fptr)->compress_type != GZIP_1) &&
 1802            ((outfptr->Fptr)->compress_type != GZIP_2)) {
 1803            switch (datatype) {
 1804             case TFLOAT:
 1805             case TDOUBLE:
 1806             case TCOMPLEX:
 1807             case TDBLCOMPLEX:
 1808               ffpmsg("Lossless compression of floating point images must use GZIP (imcomp_compress_tile)");
 1809               return(*status = DATA_COMPRESSION_ERR);
 1810             default:
 1811               break;
 1812           }
 1813        }
 1814     }
 1815 
 1816     /* free the previously saved tile if the input tile is for the same row */
 1817     if ((outfptr->Fptr)->tilerow) {  /* has the tile cache been allocated? */
 1818 
 1819       /* calculate the column bin of the compressed tile */
 1820       tilecol = (row - 1) % ((long)(((outfptr->Fptr)->znaxis[0] - 1) / ((outfptr->Fptr)->tilesize[0])) + 1);
 1821       
 1822       if ((outfptr->Fptr)->tilerow[tilecol] == row) {
 1823         if (((outfptr->Fptr)->tiledata)[tilecol]) {
 1824             free(((outfptr->Fptr)->tiledata)[tilecol]);
 1825         }
 1826       
 1827         if (((outfptr->Fptr)->tilenullarray)[tilecol]) {
 1828             free(((outfptr->Fptr)->tilenullarray)[tilecol]);
 1829         }
 1830 
 1831         ((outfptr->Fptr)->tiledata)[tilecol] = 0;
 1832         ((outfptr->Fptr)->tilenullarray)[tilecol] = 0;
 1833         (outfptr->Fptr)->tilerow[tilecol] = 0;
 1834         (outfptr->Fptr)->tiledatasize[tilecol] = 0;
 1835         (outfptr->Fptr)->tiletype[tilecol] = 0;
 1836         (outfptr->Fptr)->tileanynull[tilecol] = 0;
 1837       }
 1838     }
 1839 
 1840     if ( (outfptr->Fptr)->compress_type == NOCOMPRESS) {
 1841          /* Special case when using NOCOMPRESS for diagnostic purposes in fpack */
 1842          if (imcomp_write_nocompress_tile(outfptr, row, datatype, tiledata, tilelen, 
 1843          nullcheck, nullflagval, status) > 0) {
 1844              return(*status);
 1845          }
 1846          return(*status);
 1847     }
 1848 
 1849     /* =========================================================================== */
 1850     /* initialize various parameters */
 1851     idata = (int *) tiledata;   /* may overwrite the input tiledata in place */
 1852 
 1853     /* zbitpix is the BITPIX keyword value in the uncompressed FITS image */
 1854     zbitpix = (outfptr->Fptr)->zbitpix;
 1855 
 1856     /* if the tile/image has an integer datatype, see if a null value has */
 1857     /* been defined (with the BLANK keyword in a normal FITS image).  */
 1858     /* If so, and if the input tile array also contains null pixels, */
 1859     /* (represented by pixels that have a value = nullflagval) then  */
 1860     /* any pixels whose value = nullflagval, must be set to the value = nullval */
 1861     /* before the pixel array is compressed.  These null pixel values must */
 1862     /* not be inverse scaled by the BSCALE/BZERO values, if present. */
 1863 
 1864     cn_zblank = (outfptr->Fptr)->cn_zblank;
 1865     nullval = (outfptr->Fptr)->zblank;
 1866 
 1867     if (zbitpix > 0 && cn_zblank != -1)  /* If the integer image has no defined null */
 1868         nullcheck = 0;    /* value, then don't bother checking input array for nulls. */
 1869 
 1870     /* if the BSCALE and BZERO keywords exist, then the input values must */
 1871     /* be inverse scaled by this factor, before the values are compressed. */
 1872     /* (The program may have turned off scaling, which over rides the keywords) */
 1873     
 1874     scale = (outfptr->Fptr)->cn_bscale;
 1875     zero  = (outfptr->Fptr)->cn_bzero;
 1876     actual_bzero = (outfptr->Fptr)->cn_actual_bzero;
 1877 
 1878     /* =========================================================================== */
 1879     /* prepare the tile of pixel values for compression */
 1880     if (datatype == TSHORT) {
 1881        imcomp_convert_tile_tshort(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1882            nullval, zbitpix, scale, zero, actual_bzero, &intlength, status);
 1883     } else if (datatype == TUSHORT) {
 1884        imcomp_convert_tile_tushort(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1885            nullval, zbitpix, scale, zero, &intlength, status);
 1886     } else if (datatype == TBYTE) {
 1887        imcomp_convert_tile_tbyte(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1888            nullval, zbitpix, scale, zero,  &intlength, status);
 1889     } else if (datatype == TSBYTE) {
 1890        imcomp_convert_tile_tsbyte(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1891            nullval, zbitpix, scale, zero,  &intlength, status);
 1892     } else if (datatype == TINT) {
 1893        imcomp_convert_tile_tint(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1894            nullval, zbitpix, scale, zero, &intlength, status);
 1895     } else if (datatype == TUINT) {
 1896        imcomp_convert_tile_tuint(outfptr, tiledata, tilelen, nullcheck, nullflagval,
 1897            nullval, zbitpix, scale, zero, &intlength, status);
 1898     } else if (datatype == TLONG && sizeof(long) == 8) {
 1899            ffpmsg("Integer*8 Long datatype is not supported when writing to compressed images");
 1900            return(*status = BAD_DATATYPE);
 1901     } else if (datatype == TULONG && sizeof(long) == 8) {
 1902            ffpmsg("Unsigned integer*8 datatype is not supported when writing to compressed images");
 1903            return(*status = BAD_DATATYPE);
 1904     } else if (datatype == TFLOAT) {
 1905         imcomp_convert_tile_tfloat(outfptr, row, tiledata, tilelen, tilenx, tileny, nullcheck,
 1906         nullflagval, nullval, zbitpix, scale, zero, &intlength, &flag, bscale, bzero, status);
 1907     } else if (datatype == TDOUBLE) {
 1908        imcomp_convert_tile_tdouble(outfptr, row, tiledata, tilelen, tilenx, tileny, nullcheck,
 1909        nullflagval, nullval, zbitpix, scale, zero, &intlength, &flag, bscale, bzero, status);
 1910     } else {
 1911           ffpmsg("unsupported image datatype (imcomp_compress_tile)");
 1912           return(*status = BAD_DATATYPE);
 1913     }
 1914 
 1915     if (*status > 0)
 1916       return(*status);      /* return if error occurs */
 1917 
 1918     /* =========================================================================== */
 1919     if (flag)   /* now compress the integer data array */
 1920     {
 1921         /* allocate buffer for the compressed tile bytes */
 1922         clen = (outfptr->Fptr)->maxelem;
 1923         cbuf = (short *) calloc (clen, sizeof (unsigned char));
 1924 
 1925         if (cbuf == NULL) {
 1926             ffpmsg("Memory allocation failure. (imcomp_compress_tile)");
 1927         return (*status = MEMORY_ALLOCATION);
 1928         }
 1929 
 1930         /* =========================================================================== */
 1931         if ( (outfptr->Fptr)->compress_type == RICE_1)
 1932         {
 1933             if (intlength == 2) {
 1934             nelem = fits_rcomp_short ((short *)idata, tilelen, (unsigned char *) cbuf,
 1935                        clen, (outfptr->Fptr)->rice_blocksize);
 1936             } else if (intlength == 1) {
 1937             nelem = fits_rcomp_byte ((signed char *)idata, tilelen, (unsigned char *) cbuf,
 1938                        clen, (outfptr->Fptr)->rice_blocksize);
 1939             } else {
 1940             nelem = fits_rcomp (idata, tilelen, (unsigned char *) cbuf,
 1941                        clen, (outfptr->Fptr)->rice_blocksize);
 1942             }
 1943 
 1944         if (nelem < 0)  /* data compression error condition */
 1945             {
 1946             free (cbuf);
 1947                 ffpmsg("error Rice compressing image tile (imcomp_compress_tile)");
 1948                 return (*status = DATA_COMPRESSION_ERR);
 1949             }
 1950 
 1951         /* Write the compressed byte stream. */
 1952             ffpclb(outfptr, (outfptr->Fptr)->cn_compressed, row, 1,
 1953                      nelem, (unsigned char *) cbuf, status);
 1954         }
 1955 
 1956         /* =========================================================================== */
 1957         else if ( (outfptr->Fptr)->compress_type == PLIO_1)
 1958         {
 1959               for (ii = 0; ii < tilelen; ii++)  {
 1960                 if (idata[ii] < 0 || idata[ii] > 16777215)
 1961                 {
 1962                    /* plio algorithn only supports positive 24 bit ints */
 1963                    ffpmsg("data out of range for PLIO compression (0 - 2**24)");
 1964                    return(*status = DATA_COMPRESSION_ERR);
 1965                 }
 1966               }
 1967 
 1968           nelem = pl_p2li (idata, 1, cbuf, tilelen);
 1969 
 1970           if (nelem < 0)  /* data compression error condition */
 1971               {
 1972             free (cbuf);
 1973                 ffpmsg("error PLIO compressing image tile (imcomp_compress_tile)");
 1974                 return (*status = DATA_COMPRESSION_ERR);
 1975               }
 1976 
 1977           /* Write the compressed byte stream. */
 1978               ffpcli(outfptr, (outfptr->Fptr)->cn_compressed, row, 1,
 1979                      nelem, cbuf, status);
 1980         }
 1981 
 1982         /* =========================================================================== */
 1983         else if ( ((outfptr->Fptr)->compress_type == GZIP_1) ||
 1984                   ((outfptr->Fptr)->compress_type == GZIP_2) )   {
 1985 
 1986         if ((outfptr->Fptr)->quantize_level == NO_QUANTIZE && datatype == TFLOAT) {
 1987           /* Special case of losslessly compressing floating point pixels with GZIP */
 1988           /* In this case we compress the input tile array directly */
 1989 
 1990 #if BYTESWAPPED
 1991                ffswap4((int*) tiledata, tilelen); 
 1992 #endif
 1993                if ( (outfptr->Fptr)->compress_type == GZIP_2 )
 1994             fits_shuffle_4bytes((char *) tiledata, tilelen, status);
 1995 
 1996                 compress2mem_from_mem((char *) tiledata, tilelen * sizeof(float),
 1997                     (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 1998 
 1999         } else if ((outfptr->Fptr)->quantize_level == NO_QUANTIZE && datatype == TDOUBLE) {
 2000           /* Special case of losslessly compressing double pixels with GZIP */
 2001           /* In this case we compress the input tile array directly */
 2002 
 2003 #if BYTESWAPPED
 2004                ffswap8((double *) tiledata, tilelen); 
 2005 #endif
 2006                if ( (outfptr->Fptr)->compress_type == GZIP_2 )
 2007             fits_shuffle_8bytes((char *) tiledata, tilelen, status);
 2008 
 2009                 compress2mem_from_mem((char *) tiledata, tilelen * sizeof(double),
 2010                     (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2011 
 2012         } else {
 2013 
 2014             /* compress the integer idata array */
 2015 
 2016 #if BYTESWAPPED
 2017            if (intlength == 2)
 2018                  ffswap2((short *) idata, tilelen); 
 2019            else if (intlength == 4)
 2020                  ffswap4(idata, tilelen); 
 2021 #endif
 2022 
 2023                if (intlength == 2) {
 2024 
 2025                   if ( (outfptr->Fptr)->compress_type == GZIP_2 )
 2026             fits_shuffle_2bytes((char *) tiledata, tilelen, status);
 2027 
 2028                   compress2mem_from_mem((char *) idata, tilelen * sizeof(short),
 2029                    (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2030 
 2031                } else if (intlength == 1) {
 2032 
 2033                   compress2mem_from_mem((char *) idata, tilelen * sizeof(unsigned char),
 2034                    (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2035 
 2036                } else {
 2037 
 2038                   if ( (outfptr->Fptr)->compress_type == GZIP_2 )
 2039             fits_shuffle_4bytes((char *) tiledata, tilelen, status);
 2040 
 2041                   compress2mem_from_mem((char *) idata, tilelen * sizeof(int),
 2042                    (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2043                }
 2044             }
 2045 
 2046         /* Write the compressed byte stream. */
 2047             ffpclb(outfptr, (outfptr->Fptr)->cn_compressed, row, 1,
 2048                      gzip_nelem, (unsigned char *) cbuf, status);
 2049 
 2050         /* =========================================================================== */
 2051         } else if ( (outfptr->Fptr)->compress_type == BZIP2_1) {
 2052 
 2053 #if BYTESWAPPED
 2054        if (intlength == 2)
 2055                ffswap2((short *) idata, tilelen); 
 2056        else if (intlength == 4)
 2057                ffswap4(idata, tilelen); 
 2058 #endif
 2059 
 2060            bzlen = (unsigned int) clen;
 2061        
 2062            /* call bzip2 with blocksize = 900K, verbosity = 0, and default workfactor */
 2063 
 2064 /*  bzip2 is not supported in the public release.  This is only for test purposes.
 2065            if (BZ2_bzBuffToBuffCompress( (char *) cbuf, &bzlen,
 2066              (char *) idata, (unsigned int) (tilelen * intlength), 9, 0, 0) ) 
 2067 */
 2068        {
 2069                    ffpmsg("bzip2 compression error");
 2070                    return(*status = DATA_COMPRESSION_ERR);
 2071            }
 2072 
 2073         /* Write the compressed byte stream. */
 2074             ffpclb(outfptr, (outfptr->Fptr)->cn_compressed, row, 1,
 2075                      bzlen, (unsigned char *) cbuf, status);
 2076 
 2077         /* =========================================================================== */
 2078         }  else if ( (outfptr->Fptr)->compress_type == HCOMPRESS_1)     {
 2079         /*
 2080           if hcompscale is positive, then we have to multiply
 2081           the value by the RMS background noise to get the 
 2082           absolute scale value.  If negative, then it gives the
 2083           absolute scale value directly.
 2084         */
 2085             hcompscale = (outfptr->Fptr)->hcomp_scale;
 2086 
 2087         if (hcompscale > 0.) {
 2088            fits_img_stats_int(idata, tilenx, tileny, nullcheck,
 2089                     nullval, 0,0,0,0,0,0,&noise2,&noise3,&noise5,status);
 2090 
 2091         /* use the minimum of the 3 noise estimates */
 2092         if (noise2 != 0. && noise2 < noise3) noise3 = noise2;
 2093         if (noise5 != 0. && noise5 < noise3) noise3 = noise5;
 2094         
 2095         hcompscale = (float) (hcompscale * noise3);
 2096 
 2097         } else if (hcompscale < 0.) {
 2098 
 2099         hcompscale = hcompscale * -1.0F;
 2100         }
 2101 
 2102         ihcompscale = (int) (hcompscale + 0.5);
 2103 
 2104             hcomp_len = clen;  /* allocated size of the buffer */
 2105         
 2106             if (zbitpix == BYTE_IMG || zbitpix == SHORT_IMG) {
 2107                 fits_hcompress(idata, tilenx, tileny, 
 2108           ihcompscale, (char *) cbuf, &hcomp_len, status);
 2109 
 2110             } else {
 2111                  /* have to convert idata to an I*8 array, in place */
 2112                  /* idata must have been allocated large enough to do this */
 2113 
 2114                 fits_int_to_longlong_inplace(idata, tilelen, status);
 2115                 lldata = (LONGLONG *) idata;        
 2116 
 2117                 fits_hcompress64(lldata, tilenx, tileny, 
 2118           ihcompscale, (char *) cbuf, &hcomp_len, status);
 2119             }
 2120 
 2121         /* Write the compressed byte stream. */
 2122             ffpclb(outfptr, (outfptr->Fptr)->cn_compressed, row, 1,
 2123                      hcomp_len, (unsigned char *) cbuf, status);
 2124         }
 2125 
 2126         /* =========================================================================== */
 2127         if ((outfptr->Fptr)->cn_zscale > 0)
 2128         {
 2129               /* write the linear scaling parameters for this tile */
 2130           ffpcld (outfptr, (outfptr->Fptr)->cn_zscale, row, 1, 1,
 2131                       bscale, status);
 2132           ffpcld (outfptr, (outfptr->Fptr)->cn_zzero,  row, 1, 1,
 2133                       bzero,  status);
 2134         }
 2135 
 2136         free(cbuf);  /* finished with this buffer */
 2137 
 2138     /* =========================================================================== */
 2139     } else {    /* if flag == 0., floating point data couldn't be quantized */
 2140 
 2141      /* losslessly compress the data with gzip. */
 2142 
 2143          /* if gzip2 compressed data column doesn't exist, create it */
 2144          if ((outfptr->Fptr)->cn_gzip_data < 1) {
 2145               if ( (outfptr->Fptr)->request_huge_hdu != 0) {
 2146                  fits_insert_col(outfptr, 999, "GZIP_COMPRESSED_DATA", "1QB", status);
 2147               } else {
 2148                  fits_insert_col(outfptr, 999, "GZIP_COMPRESSED_DATA", "1PB", status);
 2149               }
 2150 
 2151                  if (*status <= 0)  /* save the number of this column */
 2152                        ffgcno(outfptr, CASEINSEN, "GZIP_COMPRESSED_DATA",
 2153                                 &(outfptr->Fptr)->cn_gzip_data, status);
 2154          }
 2155 
 2156          if (datatype == TFLOAT)  {
 2157                /* allocate buffer for the compressed tile bytes */
 2158            /* make it 10% larger than the original uncompressed data */
 2159                clen = (size_t) (tilelen * sizeof(float) * 1.1);
 2160                cbuf = (short *) calloc (clen, sizeof (unsigned char));
 2161 
 2162                if (cbuf == NULL)
 2163                {
 2164                    ffpmsg("Memory allocation error. (imcomp_compress_tile)");
 2165                return (*status = MEMORY_ALLOCATION);
 2166                }
 2167 
 2168            /* convert null values to NaNs in place, if necessary */
 2169            if (nullcheck == 1) {
 2170                imcomp_float2nan((float *) tiledata, tilelen, (int *) tiledata,
 2171                    *(float *) (nullflagval), status);
 2172            }
 2173 
 2174 #if BYTESWAPPED
 2175                ffswap4((int*) tiledata, tilelen); 
 2176 #endif
 2177                compress2mem_from_mem((char *) tiledata, tilelen * sizeof(float),
 2178                     (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2179 
 2180          } else {  /* datatype == TDOUBLE */
 2181 
 2182                /* allocate buffer for the compressed tile bytes */
 2183            /* make it 10% larger than the original uncompressed data */
 2184                clen = (size_t) (tilelen * sizeof(double) * 1.1);
 2185                cbuf = (short *) calloc (clen, sizeof (unsigned char));
 2186 
 2187                if (cbuf == NULL)
 2188                {
 2189                    ffpmsg("Memory allocation error. (imcomp_compress_tile)");
 2190                return (*status = MEMORY_ALLOCATION);
 2191                }
 2192 
 2193            /* convert null values to NaNs in place, if necessary */
 2194            if (nullcheck == 1) {
 2195                imcomp_double2nan((double *) tiledata, tilelen, (LONGLONG *) tiledata,
 2196                    *(double *) (nullflagval), status);
 2197            }
 2198 
 2199 #if BYTESWAPPED
 2200                ffswap8((double*) tiledata, tilelen); 
 2201 #endif
 2202                compress2mem_from_mem((char *) tiledata, tilelen * sizeof(double),
 2203                     (char **) &cbuf,  &clen, realloc, &gzip_nelem, status);
 2204         }
 2205 
 2206     /* Write the compressed byte stream. */
 2207         ffpclb(outfptr, (outfptr->Fptr)->cn_gzip_data, row, 1,
 2208              gzip_nelem, (unsigned char *) cbuf, status);
 2209 
 2210         free(cbuf);  /* finished with this buffer */
 2211     }
 2212 
 2213     return(*status);
 2214 }
 2215 
 2216 /*--------------------------------------------------------------------------*/
 2217 int imcomp_write_nocompress_tile(fitsfile *outfptr,
 2218     long row,
 2219     int datatype, 
 2220     void *tiledata, 
 2221     long tilelen,
 2222     int nullcheck,
 2223     void *nullflagval,
 2224     int *status)
 2225 {
 2226     char coltype[4];
 2227 
 2228     /* Write the uncompressed image tile pixels to the tile-compressed image file. */
 2229     /* This is a special case when using NOCOMPRESS for diagnostic purposes in fpack. */ 
 2230     /* Currently, this only supports a limited number of data types and */
 2231     /* does not fully support null-valued pixels in the image. */
 2232 
 2233     if ((outfptr->Fptr)->cn_uncompressed < 1) {
 2234         /* uncompressed data column doesn't exist, so append new column to table */
 2235         if (datatype == TSHORT) {
 2236         strcpy(coltype, "1PI");
 2237     } else if (datatype == TINT) {
 2238         strcpy(coltype, "1PJ");
 2239     } else if (datatype == TFLOAT) {
 2240         strcpy(coltype, "1QE");
 2241         } else {
 2242         ffpmsg("NOCOMPRESSION option only supported for int*2, int*4, and float*4 images");
 2243             return(*status = DATA_COMPRESSION_ERR);
 2244         }
 2245 
 2246         fits_insert_col(outfptr, 999, "UNCOMPRESSED_DATA", coltype, status); /* create column */
 2247     }
 2248 
 2249     fits_get_colnum(outfptr, CASEINSEN, "UNCOMPRESSED_DATA",
 2250                     &(outfptr->Fptr)->cn_uncompressed, status);  /* save col. num. */
 2251     
 2252     fits_write_col(outfptr, datatype, (outfptr->Fptr)->cn_uncompressed, row, 1,
 2253                       tilelen, tiledata, status);  /* write the tile data */
 2254     return (*status);
 2255 }
 2256  /*--------------------------------------------------------------------------*/
 2257 int imcomp_convert_tile_tshort(
 2258     fitsfile *outfptr,
 2259     void *tiledata, 
 2260     long tilelen,
 2261     int nullcheck,
 2262     void *nullflagval,
 2263     int nullval,
 2264     int zbitpix,
 2265     double scale,
 2266     double zero,
 2267     double actual_bzero,
 2268     int *intlength,
 2269     int *status)
 2270 {
 2271     /*  Prepare the input tile array of pixels for compression. */
 2272     /*  Convert input integer*2 tile array in place to 4 or 8-byte ints for compression, */
 2273     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2274     /*  Note that the calling routine must have allocated the input array big enough */
 2275     /* to be able to do this.  */
 2276 
 2277     short *sbuff;
 2278     int flagval, *idata;
 2279     long ii;
 2280     
 2281        /* We only support writing this integer*2 tile data to a FITS image with 
 2282           BITPIX = 16 and with BZERO = 0 and BSCALE = 1.  */
 2283       
 2284        if (zbitpix != SHORT_IMG || scale != 1.0 || zero != 0.0) {
 2285            ffpmsg("Datatype conversion/scaling is not supported when writing to compressed images");
 2286            return(*status = DATA_COMPRESSION_ERR);
 2287        } 
 2288 
 2289        sbuff = (short *) tiledata;
 2290        idata = (int *) tiledata;
 2291        
 2292        if ( (outfptr->Fptr)->compress_type == RICE_1 || (outfptr->Fptr)->compress_type == GZIP_1
 2293          || (outfptr->Fptr)->compress_type == GZIP_2 || (outfptr->Fptr)->compress_type == BZIP2_1 ) 
 2294        {
 2295            /* don't have to convert to int if using gzip, bzip2 or Rice compression */
 2296            *intlength = 2;
 2297              
 2298            if (nullcheck == 1) {
 2299                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2300                flagval = *(short *) (nullflagval);
 2301                if (flagval != nullval) {
 2302                   for (ii = tilelen - 1; ii >= 0; ii--) {
 2303                 if (sbuff[ii] == (short) flagval)
 2304                sbuff[ii] = (short) nullval;
 2305                   }
 2306                }
 2307            }
 2308        } else if ((outfptr->Fptr)->compress_type == HCOMPRESS_1) {
 2309            /* have to convert to int if using HCOMPRESS */
 2310            *intlength = 4;
 2311 
 2312            if (nullcheck == 1) {
 2313                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2314                flagval = *(short *) (nullflagval);
 2315                for (ii = tilelen - 1; ii >= 0; ii--) {
 2316                 if (sbuff[ii] == (short) flagval)
 2317                idata[ii] = nullval;
 2318                     else 
 2319                        idata[ii] = (int) sbuff[ii];
 2320                }
 2321            } else {  /* just do the data type conversion to int */
 2322                  /* have to convert sbuff to an I*4 array, in place */
 2323                  /* sbuff must have been allocated large enough to do this */
 2324                  fits_short_to_int_inplace(sbuff, tilelen, 0, status);
 2325            }
 2326        } else {
 2327            /* have to convert to int if using PLIO */
 2328            *intlength = 4;
 2329            if (zero == 0. && actual_bzero == 32768.) {
 2330              /* Here we are compressing unsigned 16-bit integers that have */
 2331          /* been offset by -32768 using the standard FITS convention. */
 2332          /* Since PLIO cannot deal with negative values, we must apply */
 2333          /* the shift of 32786 to the values to make them all positive. */
 2334          /* The inverse negative shift will be applied in */
 2335          /* imcomp_decompress_tile when reading the compressed tile. */
 2336              if (nullcheck == 1) {
 2337                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2338                flagval = *(short *) (nullflagval);
 2339                for (ii = tilelen - 1; ii >= 0; ii--) {
 2340                 if (sbuff[ii] == (short) flagval)
 2341                idata[ii] = nullval;
 2342                     else
 2343                        idata[ii] = (int) sbuff[ii] + 32768;
 2344                }
 2345              } else {  
 2346                  /* have to convert sbuff to an I*4 array, in place */
 2347                  /* sbuff must have been allocated large enough to do this */
 2348                  fits_short_to_int_inplace(sbuff, tilelen, 32768, status);
 2349              }
 2350            } else {
 2351          /* This is not an unsigned 16-bit integer array, so process normally */
 2352              if (nullcheck == 1) {
 2353                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2354                flagval = *(short *) (nullflagval);
 2355                for (ii = tilelen - 1; ii >= 0; ii--) {
 2356                 if (sbuff[ii] == (short) flagval)
 2357                idata[ii] = nullval;
 2358                     else
 2359                        idata[ii] = (int) sbuff[ii];
 2360                }
 2361              } else {  /* just do the data type conversion to int */
 2362                  /* have to convert sbuff to an I*4 array, in place */
 2363                  /* sbuff must have been allocated large enough to do this */
 2364                  fits_short_to_int_inplace(sbuff, tilelen, 0, status);
 2365              }
 2366            }
 2367         }
 2368         return(*status);
 2369 }
 2370  /*--------------------------------------------------------------------------*/
 2371 int imcomp_convert_tile_tushort(
 2372     fitsfile *outfptr,
 2373     void *tiledata, 
 2374     long tilelen,
 2375     int nullcheck,
 2376     void *nullflagval,
 2377     int nullval,
 2378     int zbitpix,
 2379     double scale,
 2380     double zero,
 2381     int *intlength,
 2382     int *status)
 2383 {
 2384     /*  Prepare the input  tile array of pixels for compression. */
 2385     /*  Convert input unsigned integer*2 tile array in place to 4 or 8-byte ints for compression, */
 2386     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2387     /*  Note that the calling routine must have allocated the input array big enough */
 2388     /* to be able to do this.  */
 2389 
 2390     unsigned short *usbuff;
 2391     short *sbuff;
 2392     int flagval, *idata;
 2393     long ii;
 2394     
 2395        /* datatype of input array is unsigned short.  We only support writing this datatype
 2396           to a FITS image with BITPIX = 16 and with BZERO = 0 and BSCALE = 32768.  */
 2397 
 2398        if (zbitpix != SHORT_IMG || scale != 1.0 || zero != 32768.) {
 2399            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2400            return(*status = DATA_COMPRESSION_ERR);
 2401        } 
 2402 
 2403        usbuff = (unsigned short *) tiledata;
 2404        sbuff = (short *) tiledata;
 2405        idata = (int *) tiledata;
 2406 
 2407        if ((outfptr->Fptr)->compress_type == RICE_1 || (outfptr->Fptr)->compress_type == GZIP_1
 2408         || (outfptr->Fptr)->compress_type == GZIP_2 || (outfptr->Fptr)->compress_type == BZIP2_1) 
 2409        {
 2410            /* don't have to convert to int if using gzip, bzip2, or Rice compression */
 2411            *intlength = 2;
 2412 
 2413           /* offset the unsigned value by -32768 to a signed short value. */
 2414       /* It is more efficient to do this by just flipping the most significant of the 16 bits */
 2415 
 2416            if (nullcheck == 1) {
 2417                /* reset pixels equal to flagval to the FITS null value, prior to compression  */
 2418                flagval = *(unsigned short *) (nullflagval);
 2419                for (ii = tilelen - 1; ii >= 0; ii--) {
 2420                 if (usbuff[ii] == (unsigned short) flagval)
 2421                sbuff[ii] = (short) nullval;
 2422                     else
 2423                usbuff[ii] =  (usbuff[ii]) ^ 0x8000;
 2424                }
 2425            } else {
 2426                /* just offset the pixel values by 32768 (by flipping the MSB */
 2427                for (ii = tilelen - 1; ii >= 0; ii--)
 2428                usbuff[ii] =  (usbuff[ii]) ^ 0x8000;
 2429            }
 2430        } else {
 2431            /* have to convert to int if using HCOMPRESS or PLIO */
 2432            *intlength = 4;
 2433 
 2434            if (nullcheck == 1) {
 2435                /* offset the pixel values by 32768, and */
 2436                /* reset pixels equal to flagval to nullval */
 2437                flagval = *(unsigned short *) (nullflagval);
 2438                for (ii = tilelen - 1; ii >= 0; ii--) {
 2439                 if (usbuff[ii] == (unsigned short) flagval)
 2440                idata[ii] = nullval;
 2441                     else
 2442                idata[ii] = ((int) usbuff[ii]) - 32768;
 2443                }
 2444            } else {  /* just do the data type conversion to int */
 2445                /* for HCOMPRESS we need to simply subtract 32768 */
 2446                /* for PLIO, have to convert usbuff to an I*4 array, in place */
 2447                /* usbuff must have been allocated large enough to do this */
 2448 
 2449                if ((outfptr->Fptr)->compress_type == HCOMPRESS_1) {
 2450                     fits_ushort_to_int_inplace(usbuff, tilelen, -32768, status);
 2451                } else {
 2452                     fits_ushort_to_int_inplace(usbuff, tilelen, 0, status);
 2453                }
 2454            }
 2455         }
 2456 
 2457         return(*status);
 2458 }
 2459  /*--------------------------------------------------------------------------*/
 2460 int imcomp_convert_tile_tint(
 2461     fitsfile *outfptr,
 2462     void *tiledata, 
 2463     long tilelen,
 2464     int nullcheck,
 2465     void *nullflagval,
 2466     int nullval,
 2467     int zbitpix,
 2468     double scale,
 2469     double zero,
 2470     int *intlength,
 2471     int *status)
 2472 {
 2473     /*  Prepare the input tile array of pixels for compression. */
 2474     /*  Convert input integer tile array in place to 4 or 8-byte ints for compression, */
 2475     /*  If needed, do null value substitution. */
 2476    
 2477     int flagval, *idata;
 2478     long ii;
 2479     
 2480  
 2481         /* datatype of input array is int.  We only support writing this datatype
 2482            to a FITS image with BITPIX = 32 and with BZERO = 0 and BSCALE = 1.  */
 2483 
 2484        if (zbitpix != LONG_IMG || scale != 1.0 || zero != 0.) {
 2485            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2486            return(*status = DATA_COMPRESSION_ERR);
 2487        } 
 2488 
 2489        idata = (int *) tiledata;
 2490        *intlength = 4;
 2491 
 2492        if (nullcheck == 1) {
 2493                /* no datatype conversion is required for any of the compression algorithms,
 2494              except possibly for HCOMPRESS (to I*8), which is handled later.
 2495          Just reset pixels equal to flagval to the FITS null value */
 2496                flagval = *(int *) (nullflagval);
 2497                if (flagval != nullval) {
 2498                   for (ii = tilelen - 1; ii >= 0; ii--) {
 2499                 if (idata[ii] == flagval)
 2500                idata[ii] = nullval;
 2501                   }
 2502                }
 2503        }
 2504 
 2505        return(*status);
 2506 }
 2507  /*--------------------------------------------------------------------------*/
 2508 int imcomp_convert_tile_tuint(
 2509     fitsfile *outfptr,
 2510     void *tiledata, 
 2511     long tilelen,
 2512     int nullcheck,
 2513     void *nullflagval,
 2514     int nullval,
 2515     int zbitpix,
 2516     double scale,
 2517     double zero,
 2518     int *intlength,
 2519     int *status)
 2520 {
 2521     /*  Prepare the input tile array of pixels for compression. */
 2522     /*  Convert input unsigned integer tile array in place to 4 or 8-byte ints for compression, */
 2523     /*  If needed, do null value substitution. */
 2524 
 2525 
 2526     int *idata;
 2527     unsigned int *uintbuff, uintflagval;
 2528     long ii;
 2529  
 2530        /* datatype of input array is unsigned int.  We only support writing this datatype
 2531           to a FITS image with BITPIX = 32 and with BZERO = 0 and BSCALE = 2147483648.  */
 2532 
 2533        if (zbitpix != LONG_IMG || scale != 1.0 || zero != 2147483648.) {
 2534            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2535            return(*status = DATA_COMPRESSION_ERR);
 2536        } 
 2537 
 2538        *intlength = 4;
 2539        idata = (int *) tiledata;
 2540        uintbuff = (unsigned int *) tiledata;
 2541 
 2542        /* offset the unsigned value by -2147483648 to a signed int value. */
 2543        /* It is more efficient to do this by just flipping the most significant of the 32 bits */
 2544 
 2545        if (nullcheck == 1) {
 2546                /* reset pixels equal to flagval to nullval and */
 2547                /* offset the other pixel values (by flipping the MSB) */
 2548                uintflagval = *(unsigned int *) (nullflagval);
 2549                for (ii = tilelen - 1; ii >= 0; ii--) {
 2550                 if (uintbuff[ii] == uintflagval)
 2551                idata[ii] = nullval;
 2552                     else
 2553                uintbuff[ii] = (uintbuff[ii]) ^ 0x80000000;
 2554                }
 2555        } else {
 2556                /* just offset the pixel values (by flipping the MSB) */
 2557                for (ii = tilelen - 1; ii >= 0; ii--)
 2558                uintbuff[ii] = (uintbuff[ii]) ^ 0x80000000;
 2559        }
 2560 
 2561        return(*status);
 2562 }
 2563  /*--------------------------------------------------------------------------*/
 2564 int imcomp_convert_tile_tbyte(
 2565     fitsfile *outfptr,
 2566     void *tiledata, 
 2567     long tilelen,
 2568     int nullcheck,
 2569     void *nullflagval,
 2570     int nullval,
 2571     int zbitpix,
 2572     double scale,
 2573     double zero,
 2574     int *intlength,
 2575     int *status)
 2576 {
 2577     /*  Prepare the input tile array of pixels for compression. */
 2578     /*  Convert input unsigned integer*1 tile array in place to 4 or 8-byte ints for compression, */
 2579     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2580     /*  Note that the calling routine must have allocated the input array big enough */
 2581     /* to be able to do this.  */
 2582 
 2583     int flagval, *idata;
 2584     long ii;
 2585     unsigned char *usbbuff;
 2586         
 2587        /* datatype of input array is unsigned byte.  We only support writing this datatype
 2588           to a FITS image with BITPIX = 8 and with BZERO = 0 and BSCALE = 1.  */
 2589 
 2590        if (zbitpix != BYTE_IMG || scale != 1.0 || zero != 0.) {
 2591            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2592            return(*status = DATA_COMPRESSION_ERR);
 2593        } 
 2594 
 2595        idata = (int *) tiledata;
 2596        usbbuff = (unsigned char *) tiledata;
 2597 
 2598        if ( (outfptr->Fptr)->compress_type == RICE_1 || (outfptr->Fptr)->compress_type == GZIP_1
 2599          || (outfptr->Fptr)->compress_type == GZIP_2 || (outfptr->Fptr)->compress_type == BZIP2_1 ) 
 2600        {
 2601            /* don't have to convert to int if using gzip, bzip2, or Rice compression */
 2602            *intlength = 1;
 2603              
 2604            if (nullcheck == 1) {
 2605                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2606                flagval = *(unsigned char *) (nullflagval);
 2607                if (flagval != nullval) {
 2608                   for (ii = tilelen - 1; ii >= 0; ii--) {
 2609                 if (usbbuff[ii] == (unsigned char) flagval)
 2610                usbbuff[ii] = (unsigned char) nullval;
 2611                     }
 2612                }
 2613            }
 2614        } else {
 2615            /* have to convert to int if using HCOMPRESS or PLIO */
 2616            *intlength = 4;
 2617 
 2618            if (nullcheck == 1) {
 2619                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2620                flagval = *(unsigned char *) (nullflagval);
 2621                for (ii = tilelen - 1; ii >= 0; ii--) {
 2622                 if (usbbuff[ii] == (unsigned char) flagval)
 2623                idata[ii] = nullval;
 2624                     else
 2625                        idata[ii] = (int) usbbuff[ii];
 2626                }
 2627            } else {  /* just do the data type conversion to int */
 2628                  /* have to convert usbbuff to an I*4 array, in place */
 2629                  /* usbbuff must have been allocated large enough to do this */
 2630                  fits_ubyte_to_int_inplace(usbbuff, tilelen, status);
 2631            }
 2632        }
 2633 
 2634        return(*status);
 2635 }
 2636  /*--------------------------------------------------------------------------*/
 2637 int imcomp_convert_tile_tsbyte(
 2638     fitsfile *outfptr,
 2639     void *tiledata, 
 2640     long tilelen,
 2641     int nullcheck,
 2642     void *nullflagval,
 2643     int nullval,
 2644     int zbitpix,
 2645     double scale,
 2646     double zero,
 2647     int *intlength,
 2648     int *status)
 2649 {
 2650     /*  Prepare the input tile array of pixels for compression. */
 2651     /*  Convert input integer*1 tile array in place to 4 or 8-byte ints for compression, */
 2652     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2653     /*  Note that the calling routine must have allocated the input array big enough */
 2654     /* to be able to do this.  */
 2655 
 2656     int flagval, *idata;
 2657     long ii;
 2658     signed char *sbbuff;
 2659 
 2660        /* datatype of input array is signed byte.  We only support writing this datatype
 2661           to a FITS image with BITPIX = 8 and with BZERO = 0 and BSCALE = -128.  */
 2662 
 2663        if (zbitpix != BYTE_IMG|| scale != 1.0 || zero != -128.) {
 2664            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2665            return(*status = DATA_COMPRESSION_ERR);
 2666        }
 2667 
 2668        idata = (int *) tiledata;
 2669        sbbuff = (signed char *) tiledata;
 2670 
 2671        if ( (outfptr->Fptr)->compress_type == RICE_1 || (outfptr->Fptr)->compress_type == GZIP_1
 2672          || (outfptr->Fptr)->compress_type == GZIP_2 || (outfptr->Fptr)->compress_type == BZIP2_1 ) 
 2673        {
 2674            /* don't have to convert to int if using gzip, bzip2 or Rice compression */
 2675            *intlength = 1;
 2676              
 2677            if (nullcheck == 1) {
 2678                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2679                /* offset the other pixel values (by flipping the MSB) */
 2680 
 2681                flagval = *(signed char *) (nullflagval);
 2682                for (ii = tilelen - 1; ii >= 0; ii--) {
 2683                 if (sbbuff[ii] == (signed char) flagval)
 2684                sbbuff[ii] = (signed char) nullval;
 2685                     else
 2686                sbbuff[ii] = (sbbuff[ii]) ^ 0x80;               }
 2687            } else {  /* just offset the pixel values (by flipping the MSB) */
 2688                for (ii = tilelen - 1; ii >= 0; ii--) 
 2689                sbbuff[ii] = (sbbuff[ii]) ^ 0x80;
 2690            }
 2691 
 2692        } else {
 2693            /* have to convert to int if using HCOMPRESS or PLIO */
 2694            *intlength = 4;
 2695 
 2696            if (nullcheck == 1) {
 2697                /* reset pixels equal to flagval to the FITS null value, prior to compression */
 2698                flagval = *(signed char *) (nullflagval);
 2699                for (ii = tilelen - 1; ii >= 0; ii--) {
 2700                 if (sbbuff[ii] == (signed char) flagval)
 2701                idata[ii] = nullval;
 2702                     else
 2703                        idata[ii] = ((int) sbbuff[ii]) + 128;
 2704                }
 2705            } else {  /* just do the data type conversion to int */
 2706                  /* have to convert sbbuff to an I*4 array, in place */
 2707                  /* sbbuff must have been allocated large enough to do this */
 2708                  fits_sbyte_to_int_inplace(sbbuff, tilelen, status);
 2709            }
 2710        }
 2711  
 2712        return(*status);
 2713 }
 2714  /*--------------------------------------------------------------------------*/
 2715 int imcomp_convert_tile_tfloat(
 2716     fitsfile *outfptr,
 2717     long row,
 2718     void *tiledata, 
 2719     long tilelen,
 2720     long tilenx,
 2721     long tileny,
 2722     int nullcheck,
 2723     void *nullflagval,
 2724     int nullval,
 2725     int zbitpix,
 2726     double scale,
 2727     double zero,
 2728     int *intlength,
 2729     int *flag,
 2730     double *bscale,
 2731     double *bzero,
 2732     int *status)
 2733 {
 2734     /*  Prepare the input tile array of pixels for compression. */
 2735     /*  Convert input float tile array in place to 4 or 8-byte ints for compression, */
 2736     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2737     /*  Note that the calling routine must have allocated the input array big enough */
 2738     /* to be able to do this.  */
 2739 
 2740     int *idata;
 2741     long irow, ii;
 2742     float floatnull;
 2743     unsigned char *usbbuff;
 2744     unsigned long dithersum;
 2745     int iminval = 0, imaxval = 0;  /* min and max quantized integers */
 2746 
 2747         /* datatype of input array is double.  We only support writing this datatype
 2748            to a FITS image with BITPIX = -64 or -32, except we also support the special case where
 2749        BITPIX = 32 and BZERO = 0 and BSCALE = 1.  */
 2750 
 2751        if ((zbitpix != LONG_IMG && zbitpix != DOUBLE_IMG && zbitpix != FLOAT_IMG) || scale != 1.0 || zero != 0.) {
 2752            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2753            return(*status = DATA_COMPRESSION_ERR);
 2754        } 
 2755 
 2756            *intlength = 4;
 2757            idata = (int *) tiledata;
 2758 
 2759           /* if the tile-compressed table contains zscale and zzero columns */
 2760           /* then scale and quantize the input floating point data.    */
 2761 
 2762           if ((outfptr->Fptr)->cn_zscale > 0) {
 2763         /* quantize the float values into integers */
 2764 
 2765             if (nullcheck == 1)
 2766           floatnull = *(float *) (nullflagval);
 2767         else
 2768           floatnull = FLOATNULLVALUE;  /* NaNs are represented by this, by default */
 2769 
 2770             if ((outfptr->Fptr)->quantize_method == SUBTRACTIVE_DITHER_1  ||
 2771             (outfptr->Fptr)->quantize_method == SUBTRACTIVE_DITHER_2) {
 2772           
 2773               /* see if the dithering offset value needs to be initialized */                  
 2774               if ((outfptr->Fptr)->request_dither_seed == 0 && (outfptr->Fptr)->dither_seed == 0) {
 2775 
 2776              /* This means randomly choose the dithering offset based on the system time. */
 2777              /* The offset will have a value between 1 and 10000, inclusive. */
 2778              /* The time function returns an integer value that is incremented each second. */
 2779              /* The clock function returns the elapsed CPU time, in integer CLOCKS_PER_SEC units. */
 2780              /* The CPU time returned by clock is typically (on linux PC) only good to 0.01 sec */
 2781              /* Summing the 2 quantities may help avoid cases where 2 executions of the program */
 2782              /* (perhaps in a multithreaded environoment) end up with exactly the same dither seed */
 2783              /* value.  The sum is incremented by the current HDU number in the file to provide */
 2784              /* further randomization.  This randomization is desireable if multiple compressed */
 2785              /* images will be summed (or differenced). In such cases, the benefits of dithering */
 2786              /* may be lost if all the images use exactly the same sequence of random numbers when */
 2787              /* calculating the dithering offsets. */        
 2788              
 2789              (outfptr->Fptr)->dither_seed = 
 2790                (( (int)time(NULL) + ( (int) clock() / (int) (CLOCKS_PER_SEC / 100)) + (outfptr->Fptr)->curhdu) % 10000) + 1;
 2791              
 2792                      /* update the header keyword with this new value */
 2793              fits_update_key(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->dither_seed), 
 2794                             NULL, status);
 2795 
 2796               } else if ((outfptr->Fptr)->request_dither_seed < 0 && (outfptr->Fptr)->dither_seed < 0) {
 2797 
 2798              /* this means randomly choose the dithering offset based on some hash function */
 2799              /* of the first input tile of data to be quantized and compressed.  This ensures that */
 2800                      /* the same offset value is used for a given image every time it is compressed. */
 2801 
 2802              usbbuff = (unsigned char *) tiledata;
 2803              dithersum = 0;
 2804              for (ii = 0; ii < 4 * tilelen; ii++) {
 2805                  dithersum += usbbuff[ii];  /* doesn't matter if there is an integer overflow */
 2806                  }
 2807              (outfptr->Fptr)->dither_seed = ((int) (dithersum % 10000)) + 1;
 2808         
 2809                      /* update the header keyword with this new value */
 2810              fits_update_key(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->dither_seed), 
 2811                             NULL, status);
 2812           }
 2813 
 2814                   /* subtract 1 to convert from 1-based to 0-based element number */
 2815               irow = row + (outfptr->Fptr)->dither_seed - 1; /* dither the quantized values */
 2816 
 2817           } else if ((outfptr->Fptr)->quantize_method == -1) {
 2818               irow = 0;  /* do not dither the quantized values */
 2819               } else {
 2820                   ffpmsg("Unknown dithering method.");
 2821                   ffpmsg("May need to install a newer version of CFITSIO.");
 2822                   return(*status = DATA_COMPRESSION_ERR);
 2823               }
 2824 
 2825               *flag = fits_quantize_float (irow, (float *) tiledata, tilenx, tileny,
 2826                    nullcheck, floatnull, (outfptr->Fptr)->quantize_level, 
 2827            (outfptr->Fptr)->quantize_method, idata, bscale, bzero, &iminval, &imaxval);
 2828 
 2829               if (*flag > 1)
 2830            return(*status = *flag);
 2831           }
 2832           else if ((outfptr->Fptr)->quantize_level != NO_QUANTIZE)
 2833       {
 2834         /* if floating point pixels are not being losslessly compressed, then */
 2835         /* input float data is implicitly converted (truncated) to integers */
 2836             if ((scale != 1. || zero != 0.))  /* must scale the values */
 2837            imcomp_nullscalefloats((float *) tiledata, tilelen, idata, scale, zero,
 2838                nullcheck, *(float *) (nullflagval), nullval, status);
 2839              else
 2840            imcomp_nullfloats((float *) tiledata, tilelen, idata,
 2841                nullcheck, *(float *) (nullflagval), nullval,  status);
 2842           }
 2843           else if ((outfptr->Fptr)->quantize_level == NO_QUANTIZE)
 2844       {
 2845           /* just convert null values to NaNs in place, if necessary, then do lossless gzip compression */
 2846         if (nullcheck == 1) {
 2847                 imcomp_float2nan((float *) tiledata, tilelen, (int *) tiledata,
 2848                     *(float *) (nullflagval), status);
 2849         }
 2850           }
 2851 
 2852           return(*status);
 2853 }
 2854  /*--------------------------------------------------------------------------*/
 2855 int imcomp_convert_tile_tdouble(
 2856     fitsfile *outfptr,
 2857     long row,
 2858     void *tiledata, 
 2859     long tilelen,
 2860     long tilenx,
 2861     long tileny,
 2862     int nullcheck,
 2863     void *nullflagval,
 2864     int nullval,
 2865     int zbitpix,
 2866     double scale,
 2867     double zero,
 2868     int *intlength,
 2869     int *flag,
 2870     double *bscale,
 2871     double *bzero,
 2872     int *status)
 2873 {
 2874     /*  Prepare the input tile array of pixels for compression. */
 2875     /*  Convert input double tile array in place to 4-byte ints for compression, */
 2876     /*  If needed, convert 4 or 8-byte ints and do null value substitution. */
 2877     /*  Note that the calling routine must have allocated the input array big enough */
 2878     /* to be able to do this.  */
 2879 
 2880     int *idata;
 2881     long irow, ii;
 2882     double doublenull;
 2883     unsigned char *usbbuff;
 2884     unsigned long dithersum;
 2885     int iminval = 0, imaxval = 0;  /* min and max quantized integers */
 2886 
 2887         /* datatype of input array is double.  We only support writing this datatype
 2888            to a FITS image with BITPIX = -64 or -32, except we also support the special case where
 2889        BITPIX = 32 and BZERO = 0 and BSCALE = 1.  */
 2890 
 2891        if ((zbitpix != LONG_IMG && zbitpix != DOUBLE_IMG && zbitpix != FLOAT_IMG) || scale != 1.0 || zero != 0.) {
 2892            ffpmsg("Implicit datatype conversion is not supported when writing to compressed images");
 2893            return(*status = DATA_COMPRESSION_ERR);
 2894        } 
 2895 
 2896            *intlength = 4;
 2897            idata = (int *) tiledata;
 2898 
 2899           /* if the tile-compressed table contains zscale and zzero columns */
 2900           /* then scale and quantize the input floating point data.    */
 2901           /* Otherwise, just truncate the floats to integers.          */
 2902 
 2903           if ((outfptr->Fptr)->cn_zscale > 0)
 2904           {
 2905             if (nullcheck == 1)
 2906           doublenull = *(double *) (nullflagval);
 2907         else
 2908           doublenull = DOUBLENULLVALUE;
 2909 
 2910             /* quantize the double values into integers */
 2911               if ((outfptr->Fptr)->quantize_method == SUBTRACTIVE_DITHER_1 ||
 2912               (outfptr->Fptr)->quantize_method == SUBTRACTIVE_DITHER_2) {
 2913 
 2914               /* see if the dithering offset value needs to be initialized (see above) */                  
 2915               if ((outfptr->Fptr)->request_dither_seed == 0 && (outfptr->Fptr)->dither_seed == 0) {
 2916 
 2917              (outfptr->Fptr)->dither_seed = 
 2918                (( (int)time(NULL) + ( (int) clock() / (int) (CLOCKS_PER_SEC / 100)) + (outfptr->Fptr)->curhdu) % 10000) + 1;
 2919              
 2920                      /* update the header keyword with this new value */
 2921              fits_update_key(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->dither_seed), 
 2922                             NULL, status);
 2923 
 2924               } else if ((outfptr->Fptr)->request_dither_seed < 0 && (outfptr->Fptr)->dither_seed < 0) {
 2925 
 2926              usbbuff = (unsigned char *) tiledata;
 2927              dithersum = 0;
 2928              for (ii = 0; ii < 8 * tilelen; ii++) {
 2929                  dithersum += usbbuff[ii];
 2930                  }
 2931              (outfptr->Fptr)->dither_seed = ((int) (dithersum % 10000)) + 1;
 2932         
 2933                      /* update the header keyword with this new value */
 2934              fits_update_key(outfptr, TINT, "ZDITHER0", &((outfptr->Fptr)->dither_seed), 
 2935                             NULL, status);
 2936           }
 2937 
 2938               irow = row + (outfptr->Fptr)->dither_seed - 1; /* dither the quantized values */
 2939 
 2940           } else if ((outfptr->Fptr)->quantize_method == -1) {
 2941               irow = 0;  /* do not dither the quantized values */
 2942               } else {
 2943                   ffpmsg("Unknown subtractive dithering method.");
 2944                   ffpmsg("May need to install a newer version of CFITSIO.");
 2945                   return(*status = DATA_COMPRESSION_ERR);
 2946               }
 2947 
 2948             *flag = fits_quantize_double (irow, (double *) tiledata, tilenx, tileny,
 2949                nullcheck, doublenull, (outfptr->Fptr)->quantize_level, 
 2950            (outfptr->Fptr)->quantize_method, idata,
 2951                bscale, bzero, &iminval, &imaxval);
 2952 
 2953             if (*flag > 1)
 2954         return(*status = *flag);
 2955           }
 2956           else if ((outfptr->Fptr)->quantize_level != NO_QUANTIZE)
 2957       {
 2958         /* if floating point pixels are not being losslessly compressed, then */
 2959         /* input float data is implicitly converted (truncated) to integers */
 2960              if ((scale != 1. || zero != 0.))  /* must scale the values */
 2961            imcomp_nullscaledoubles((double *) tiledata, tilelen, idata, scale, zero,
 2962                nullcheck, *(double *) (nullflagval), nullval, status);
 2963              else
 2964            imcomp_nulldoubles((double *) tiledata, tilelen, idata,
 2965                nullcheck, *(double *) (nullflagval), nullval,  status);
 2966           }
 2967           else if ((outfptr->Fptr)->quantize_level == NO_QUANTIZE)
 2968       {
 2969           /* just convert null values to NaNs in place, if necessary, then do lossless gzip compression */
 2970         if (nullcheck == 1) {
 2971                 imcomp_double2nan((double *) tiledata, tilelen, (LONGLONG *) tiledata,
 2972                     *(double *) (nullflagval), status);
 2973         }
 2974           }
 2975  
 2976           return(*status);
 2977 }
 2978 /*---------------------------------------------------------------------------*/
 2979 int imcomp_nullscale(
 2980      int *idata, 
 2981      long tilelen,
 2982      int nullflagval,
 2983      int nullval,
 2984      double scale,
 2985      double zero,
 2986      int *status)
 2987 /*
 2988    do null value substitution AND scaling of the integer array.
 2989    If array value = nullflagval, then set the value to nullval.
 2990    Otherwise, inverse scale the integer value.
 2991 */
 2992 {
 2993     long ii;
 2994     double dvalue;
 2995     
 2996     for (ii=0; ii < tilelen; ii++)
 2997     {
 2998         if (idata[ii] == nullflagval)
 2999         idata[ii] = nullval;
 3000     else 
 3001     {
 3002             dvalue = (idata[ii] - zero) / scale;
 3003 
 3004             if (dvalue < DINT_MIN)
 3005             {
 3006                 *status = OVERFLOW_ERR;
 3007                 idata[ii] = INT32_MIN;
 3008             }
 3009             else if (dvalue > DINT_MAX)
 3010             {
 3011                 *status = OVERFLOW_ERR;
 3012                 idata[ii] = INT32_MAX;
 3013             }
 3014             else
 3015             {
 3016                 if (dvalue >= 0)
 3017                     idata[ii] = (int) (dvalue + .5);
 3018                 else
 3019                     idata[ii] = (int) (dvalue - .5);
 3020             }
 3021         }
 3022     }
 3023     return(*status);
 3024 }
 3025 /*---------------------------------------------------------------------------*/
 3026 int imcomp_nullvalues(
 3027      int *idata, 
 3028      long tilelen,
 3029      int nullflagval,
 3030      int nullval,
 3031      int *status)
 3032 /*
 3033    do null value substitution.
 3034    If array value = nullflagval, then set the value to nullval.
 3035 */
 3036 {
 3037     long ii;
 3038     
 3039     for (ii=0; ii < tilelen; ii++)
 3040     {
 3041         if (idata[ii] == nullflagval)
 3042         idata[ii] = nullval;
 3043     }
 3044     return(*status);
 3045 }
 3046 /*---------------------------------------------------------------------------*/
 3047 int imcomp_scalevalues(
 3048      int *idata, 
 3049      long tilelen,
 3050      double scale,
 3051      double zero,
 3052      int *status)
 3053 /*
 3054    do inverse scaling the integer values.
 3055 */
 3056 {
 3057     long ii;
 3058     double dvalue;
 3059     
 3060     for (ii=0; ii < tilelen; ii++)
 3061     {
 3062             dvalue = (idata[ii] - zero) / scale;
 3063 
 3064             if (dvalue < DINT_MIN)
 3065             {
 3066                 *status = OVERFLOW_ERR;
 3067                 idata[ii] = INT32_MIN;
 3068             }
 3069             else if (dvalue > DINT_MAX)
 3070             {
 3071                 *status = OVERFLOW_ERR;
 3072                 idata[ii] = INT32_MAX;
 3073             }
 3074             else
 3075             {
 3076                 if (dvalue >= 0)
 3077                     idata[ii] = (int) (dvalue + .5);
 3078                 else
 3079                     idata[ii] = (int) (dvalue - .5);
 3080             }
 3081     }
 3082     return(*status);
 3083 }
 3084 /*---------------------------------------------------------------------------*/
 3085 int imcomp_nullscalei2(
 3086      short *idata, 
 3087      long tilelen,
 3088      short nullflagval,
 3089      short nullval,
 3090      double scale,
 3091      double zero,
 3092      int *status)
 3093 /*
 3094    do null value substitution AND scaling of the integer array.
 3095    If array value = nullflagval, then set the value to nullval.
 3096    Otherwise, inverse scale the integer value.
 3097 */
 3098 {
 3099     long ii;
 3100     double dvalue;
 3101     
 3102     for (ii=0; ii < tilelen; ii++)
 3103     {
 3104         if (idata[ii] == nullflagval)
 3105         idata[ii] = nullval;
 3106     else 
 3107     {
 3108             dvalue = (idata[ii] - zero) / scale;
 3109 
 3110             if (dvalue < DSHRT_MIN)
 3111             {
 3112                 *status = OVERFLOW_ERR;
 3113                 idata[ii] = SHRT_MIN;
 3114             }
 3115             else if (dvalue > DSHRT_MAX)
 3116             {
 3117                 *status = OVERFLOW_ERR;
 3118                 idata[ii] = SHRT_MAX;
 3119             }
 3120             else
 3121             {
 3122                 if (dvalue >= 0)
 3123                     idata[ii] = (int) (dvalue + .5);
 3124                 else
 3125                     idata[ii] = (int) (dvalue - .5);
 3126             }
 3127         }
 3128     }
 3129     return(*status);
 3130 }
 3131 /*---------------------------------------------------------------------------*/
 3132 int imcomp_nullvaluesi2(
 3133      short *idata, 
 3134      long tilelen,
 3135      short nullflagval,
 3136      short nullval,
 3137      int *status)
 3138 /*
 3139    do null value substitution.
 3140    If array value = nullflagval, then set the value to nullval.
 3141 */
 3142 {
 3143     long ii;
 3144     
 3145     for (ii=0; ii < tilelen; ii++)
 3146     {
 3147         if (idata[ii] == nullflagval)
 3148         idata[ii] = nullval;
 3149     }
 3150     return(*status);
 3151 }
 3152 /*---------------------------------------------------------------------------*/
 3153 int imcomp_scalevaluesi2(
 3154      short *idata, 
 3155      long tilelen,
 3156      double scale,
 3157      double zero,
 3158      int *status)
 3159 /*
 3160    do inverse scaling the integer values.
 3161 */
 3162 {
 3163     long ii;
 3164     double dvalue;
 3165     
 3166     for (ii=0; ii < tilelen; ii++)
 3167     {
 3168             dvalue = (idata[ii] - zero) / scale;
 3169 
 3170             if (dvalue < DSHRT_MIN)
 3171             {
 3172                 *status = OVERFLOW_ERR;
 3173                 idata[ii] = SHRT_MIN;
 3174             }
 3175             else if (dvalue > DSHRT_MAX)
 3176             {
 3177                 *status = OVERFLOW_ERR;
 3178                 idata[ii] = SHRT_MAX;
 3179             }
 3180             else
 3181             {
 3182                 if (dvalue >= 0)
 3183                     idata[ii] = (int) (dvalue + .5);
 3184                 else
 3185                     idata[ii] = (int) (dvalue - .5);
 3186             }
 3187     }
 3188     return(*status);
 3189 }
 3190 /*---------------------------------------------------------------------------*/
 3191 int imcomp_nullfloats(
 3192      float *fdata,
 3193      long tilelen,
 3194      int *idata, 
 3195      int nullcheck,
 3196      float nullflagval,
 3197      int nullval,
 3198      int *status)
 3199 /*
 3200    do null value substitution  of the float array.
 3201    If array value = nullflagval, then set the output value to FLOATNULLVALUE.
 3202 */
 3203 {
 3204     long ii;
 3205     double dvalue;
 3206     
 3207     if (nullcheck == 1) /* must check for null values */
 3208     {
 3209       for (ii=0; ii < tilelen; ii++)
 3210       {
 3211         if (fdata[ii] == nullflagval)
 3212         idata[ii] = nullval;
 3213     else 
 3214     {
 3215             dvalue = fdata[ii];
 3216 
 3217             if (dvalue < DINT_MIN)
 3218             {
 3219                 *status = OVERFLOW_ERR;
 3220                 idata[ii] = INT32_MIN;
 3221             }
 3222             else if (dvalue > DINT_MAX)
 3223             {
 3224                 *status = OVERFLOW_ERR;
 3225                 idata[ii] = INT32_MAX;
 3226             }
 3227             else
 3228             {
 3229                 if (dvalue >= 0)
 3230                     idata[ii] = (int) (dvalue + .5);
 3231                 else
 3232                     idata[ii] = (int) (dvalue - .5);
 3233             }
 3234         }
 3235       }
 3236     }
 3237     else  /* don't have to worry about null values */
 3238     {
 3239       for (ii=0; ii < tilelen; ii++)
 3240       {
 3241             dvalue = fdata[ii];
 3242 
 3243             if (dvalue < DINT_MIN)
 3244             {
 3245                 *status = OVERFLOW_ERR;
 3246                 idata[ii] = INT32_MIN;
 3247             }
 3248             else if (dvalue > DINT_MAX)
 3249             {
 3250                 *status = OVERFLOW_ERR;
 3251                 idata[ii] = INT32_MAX;
 3252             }
 3253             else
 3254             {
 3255                 if (dvalue >= 0)
 3256                     idata[ii] = (int) (dvalue + .5);
 3257                 else
 3258                     idata[ii] = (int) (dvalue - .5);
 3259             }
 3260       }
 3261     }
 3262     return(*status);
 3263 }
 3264 /*---------------------------------------------------------------------------*/
 3265 int imcomp_nullscalefloats(
 3266      float *fdata,
 3267      long tilelen,
 3268      int *idata, 
 3269      double scale,
 3270      double zero,
 3271      int nullcheck,
 3272      float nullflagval,
 3273      int nullval,
 3274      int *status)
 3275 /*
 3276    do null value substitution  of the float array.
 3277    If array value = nullflagval, then set the output value to FLOATNULLVALUE.
 3278    Otherwise, inverse scale the integer value.
 3279 */
 3280 {
 3281     long ii;
 3282     double dvalue;
 3283     
 3284     if (nullcheck == 1) /* must check for null values */
 3285     {
 3286       for (ii=0; ii < tilelen; ii++)
 3287       {
 3288         if (fdata[ii] == nullflagval)
 3289         idata[ii] = nullval;
 3290     else 
 3291     {
 3292             dvalue = (fdata[ii] - zero) / scale;
 3293 
 3294             if (dvalue < DINT_MIN)
 3295             {
 3296                 *status = OVERFLOW_ERR;
 3297                 idata[ii] = INT32_MIN;
 3298             }
 3299             else if (dvalue > DINT_MAX)
 3300             {
 3301                 *status = OVERFLOW_ERR;
 3302                 idata[ii] = INT32_MAX;
 3303             }
 3304             else
 3305             {
 3306                 if (dvalue >= 0.)
 3307                     idata[ii] = (int) (dvalue + .5);
 3308                 else
 3309                     idata[ii] = (int) (dvalue - .5);
 3310             }
 3311         }
 3312       }
 3313     }
 3314     else  /* don't have to worry about null values */
 3315     {
 3316       for (ii=0; ii < tilelen; ii++)
 3317       {
 3318             dvalue = (fdata[ii] - zero) / scale;
 3319 
 3320             if (dvalue < DINT_MIN)
 3321             {
 3322                 *status = OVERFLOW_ERR;
 3323                 idata[ii] = INT32_MIN;
 3324             }
 3325             else if (dvalue > DINT_MAX)
 3326             {
 3327                 *status = OVERFLOW_ERR;
 3328                 idata[ii] = INT32_MAX;
 3329             }
 3330             else
 3331             {
 3332                 if (dvalue >= 0.)
 3333                     idata[ii] = (int) (dvalue + .5);
 3334                 else
 3335                     idata[ii] = (int) (dvalue - .5);
 3336             }
 3337       }
 3338     }
 3339     return(*status);
 3340 }
 3341 /*---------------------------------------------------------------------------*/
 3342 int imcomp_nulldoubles(
 3343      double *fdata,
 3344      long tilelen,
 3345      int *idata, 
 3346      int nullcheck,
 3347      double nullflagval,
 3348      int nullval,
 3349      int *status)
 3350 /*
 3351    do null value substitution  of the float array.
 3352    If array value = nullflagval, then set the output value to FLOATNULLVALUE.
 3353    Otherwise, inverse scale the integer value.
 3354 */
 3355 {
 3356     long ii;
 3357     double dvalue;
 3358     
 3359     if (nullcheck == 1) /* must check for null values */
 3360     {
 3361       for (ii=0; ii < tilelen; ii++)
 3362       {
 3363         if (fdata[ii] == nullflagval)
 3364         idata[ii] = nullval;
 3365     else 
 3366     {
 3367             dvalue = fdata[ii];
 3368 
 3369             if (dvalue < DINT_MIN)
 3370             {
 3371                 *status = OVERFLOW_ERR;
 3372                 idata[ii] = INT32_MIN;
 3373             }
 3374             else if (dvalue > DINT_MAX)
 3375             {
 3376                 *status = OVERFLOW_ERR;
 3377                 idata[ii] = INT32_MAX;
 3378             }
 3379             else
 3380             {
 3381                 if (dvalue >= 0.)
 3382                     idata[ii] = (int) (dvalue + .5);
 3383                 else
 3384                     idata[ii] = (int) (dvalue - .5);
 3385             }
 3386         }
 3387       }
 3388     }
 3389     else  /* don't have to worry about null values */
 3390     {
 3391       for (ii=0; ii < tilelen; ii++)
 3392       {
 3393             dvalue = fdata[ii];
 3394 
 3395             if (dvalue < DINT_MIN)
 3396             {
 3397                 *status = OVERFLOW_ERR;
 3398                 idata[ii] = INT32_MIN;
 3399             }
 3400             else if (dvalue > DINT_MAX)
 3401             {
 3402                 *status = OVERFLOW_ERR;
 3403                 idata[ii] = INT32_MAX;
 3404             }
 3405             else
 3406             {
 3407                 if (dvalue >= 0.)
 3408                     idata[ii] = (int) (dvalue + .5);
 3409                 else
 3410                     idata[ii] = (int) (dvalue - .5);
 3411             }
 3412       }
 3413     }
 3414     return(*status);
 3415 }
 3416 /*---------------------------------------------------------------------------*/
 3417 int imcomp_nullscaledoubles(
 3418      double *fdata,
 3419      long tilelen,
 3420      int *idata, 
 3421      double scale,
 3422      double zero,
 3423      int nullcheck,
 3424      double nullflagval,
 3425      int nullval,
 3426      int *status)
 3427 /*
 3428    do null value substitution  of the float array.
 3429    If array value = nullflagval, then set the output value to FLOATNULLVALUE.
 3430    Otherwise, inverse scale the integer value.
 3431 */
 3432 {
 3433     long ii;
 3434     double dvalue;
 3435     
 3436     if (nullcheck == 1) /* must check for null values */
 3437     {
 3438       for (ii=0; ii < tilelen; ii++)
 3439       {
 3440         if (fdata[ii] == nullflagval)
 3441         idata[ii] = nullval;
 3442     else 
 3443     {
 3444             dvalue = (fdata[ii] - zero) / scale;
 3445 
 3446             if (dvalue < DINT_MIN)
 3447             {
 3448                 *status = OVERFLOW_ERR;
 3449                 idata[ii] = INT32_MIN;
 3450             }
 3451             else if (dvalue > DINT_MAX)
 3452             {
 3453                 *status = OVERFLOW_ERR;
 3454                 idata[ii] = INT32_MAX;
 3455             }
 3456             else
 3457             {
 3458                 if (dvalue >= 0.)
 3459                     idata[ii] = (int) (dvalue + .5);
 3460                 else
 3461                     idata[ii] = (int) (dvalue - .5);
 3462             }
 3463         }
 3464       }
 3465     }
 3466     else  /* don't have to worry about null values */
 3467     {
 3468       for (ii=0; ii < tilelen; ii++)
 3469       {
 3470             dvalue = (fdata[ii] - zero) / scale;
 3471 
 3472             if (dvalue < DINT_MIN)
 3473             {
 3474                 *status = OVERFLOW_ERR;
 3475                 idata[ii] = INT32_MIN;
 3476             }
 3477             else if (dvalue > DINT_MAX)
 3478             {
 3479                 *status = OVERFLOW_ERR;
 3480                 idata[ii] = INT32_MAX;
 3481             }
 3482             else
 3483             {
 3484                 if (dvalue >= 0.)
 3485                     idata[ii] = (int) (dvalue + .5);
 3486                 else
 3487                     idata[ii] = (int) (dvalue - .5);
 3488             }
 3489       }
 3490     }
 3491     return(*status);
 3492 }
 3493 /*---------------------------------------------------------------------------*/
 3494 int fits_write_compressed_img(fitsfile *fptr,   /* I - FITS file pointer     */
 3495             int  datatype,   /* I - datatype of the array to be written      */
 3496             long  *infpixel, /* I - 'bottom left corner' of the subsection   */
 3497             long  *inlpixel, /* I - 'top right corner' of the subsection     */
 3498             int  nullcheck,  /* I - 0 for no null checking                   */
 3499                              /*     1: pixels that are = nullval will be     */
 3500                              /*     written with the FITS null pixel value   */
 3501                              /*     (floating point arrays only)             */
 3502             void *array,     /* I - array of values to be written            */
 3503             void *nullval,   /* I - undefined pixel value                    */
 3504             int  *status)    /* IO - error status                            */
 3505 /*
 3506    Write a section of a compressed image.
 3507 */
 3508 {
 3509     int  tiledim[MAX_COMPRESS_DIM];
 3510     long naxis[MAX_COMPRESS_DIM];
 3511     long tilesize[MAX_COMPRESS_DIM], thistilesize[MAX_COMPRESS_DIM];
 3512     long ftile[MAX_COMPRESS_DIM], ltile[MAX_COMPRESS_DIM];
 3513     long tfpixel[MAX_COMPRESS_DIM], tlpixel[MAX_COMPRESS_DIM];
 3514     long rowdim[MAX_COMPRESS_DIM], offset[MAX_COMPRESS_DIM],ntemp;
 3515     long fpixel[MAX_COMPRESS_DIM], lpixel[MAX_COMPRESS_DIM];
 3516     long i5, i4, i3, i2, i1, i0, irow, trowsize, ntrows;
 3517     int ii, ndim, pixlen, tilenul;
 3518     int  tstatus, buffpixsiz;
 3519     void *buffer;
 3520     char *bnullarray = 0, card[FLEN_CARD];
 3521 
 3522     if (*status > 0) 
 3523         return(*status);
 3524 
 3525     if (!fits_is_compressed_image(fptr, status) )
 3526     {
 3527         ffpmsg("CHDU is not a compressed image (fits_write_compressed_img)");
 3528         return(*status = DATA_COMPRESSION_ERR);
 3529     }
 3530 
 3531     /* reset position to the correct HDU if necessary */
 3532     if (fptr->HDUposition != (fptr->Fptr)->curhdu)
 3533         ffmahd(fptr, (fptr->HDUposition) + 1, NULL, status);
 3534 
 3535     /* rescan header if data structure is undefined */
 3536     else if ((fptr->Fptr)->datastart == DATA_UNDEFINED)
 3537         if ( ffrdef(fptr, status) > 0)               
 3538             return(*status);
 3539 
 3540 
 3541     /* ===================================================================== */
 3542 
 3543 
 3544     if (datatype == TSHORT || datatype == TUSHORT)
 3545     {
 3546        pixlen = sizeof(short);
 3547     }
 3548     else if (datatype == TINT || datatype == TUINT)
 3549     {
 3550        pixlen = sizeof(int);
 3551     }
 3552     else if (datatype == TBYTE || datatype == TSBYTE)
 3553     {
 3554        pixlen = 1;
 3555     }
 3556     else if (datatype == TLONG || datatype == TULONG)
 3557     {
 3558        pixlen = sizeof(long);
 3559     }
 3560     else if (datatype == TFLOAT)
 3561     {
 3562        pixlen = sizeof(float);
 3563     }
 3564     else if (datatype == TDOUBLE)
 3565     {
 3566        pixlen = sizeof(double);
 3567     }
 3568     else
 3569     {
 3570         ffpmsg("unsupported datatype for compressing image");
 3571         return(*status = BAD_DATATYPE);
 3572     }
 3573 
 3574     /* ===================================================================== */
 3575 
 3576     /* allocate scratch space for processing one tile of the image */
 3577     buffpixsiz = pixlen;  /* this is the minimum pixel size */
 3578     
 3579     if ( (fptr->Fptr)->compress_type == HCOMPRESS_1) { /* need 4 or 8 bytes per pixel */
 3580         if ((fptr->Fptr)->zbitpix == BYTE_IMG ||
 3581         (fptr->Fptr)->zbitpix == SHORT_IMG )
 3582                 buffpixsiz = maxvalue(buffpixsiz, 4);
 3583         else
 3584             buffpixsiz = 8;
 3585     }
 3586     else if ( (fptr->Fptr)->compress_type == PLIO_1) { /* need 4 bytes per pixel */
 3587                 buffpixsiz = maxvalue(buffpixsiz, 4);
 3588     }
 3589     else if ( (fptr->Fptr)->compress_type == RICE_1  ||
 3590               (fptr->Fptr)->compress_type == GZIP_1 ||
 3591               (fptr->Fptr)->compress_type == GZIP_2 ||
 3592               (fptr->Fptr)->compress_type == BZIP2_1) {  /* need 1, 2, or 4 bytes per pixel */
 3593         if ((fptr->Fptr)->zbitpix == BYTE_IMG)
 3594             buffpixsiz = maxvalue(buffpixsiz, 1);
 3595         else if ((fptr->Fptr)->zbitpix == SHORT_IMG)
 3596             buffpixsiz = maxvalue(buffpixsiz, 2);
 3597         else 
 3598             buffpixsiz = maxvalue(buffpixsiz, 4);
 3599     }
 3600     else
 3601     {
 3602         ffpmsg("unsupported image compression algorithm");
 3603         return(*status = BAD_DATATYPE);
 3604     }
 3605     
 3606     /* cast to double to force alignment on 8-byte addresses */
 3607     buffer = (double *) calloc ((fptr->Fptr)->maxtilelen, buffpixsiz);
 3608 
 3609     if (buffer == NULL)
 3610     {
 3611         ffpmsg("Out of memory (fits_write_compress_img)");
 3612         return (*status = MEMORY_ALLOCATION);
 3613     }
 3614 
 3615     /* ===================================================================== */
 3616 
 3617     /* initialize all the arrays */
 3618     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 3619     {
 3620         naxis[ii] = 1;
 3621         tiledim[ii] = 1;
 3622         tilesize[ii] = 1;
 3623         ftile[ii] = 1;
 3624         ltile[ii] = 1;
 3625         rowdim[ii] = 1;
 3626     }
 3627 
 3628     ndim = (fptr->Fptr)->zndim;
 3629     ntemp = 1;
 3630     for (ii = 0; ii < ndim; ii++)
 3631     {
 3632         fpixel[ii] = infpixel[ii];
 3633         lpixel[ii] = inlpixel[ii];
 3634 
 3635         /* calc number of tiles in each dimension, and tile containing */
 3636         /* the first and last pixel we want to read in each dimension  */
 3637         naxis[ii] = (fptr->Fptr)->znaxis[ii];
 3638         if (fpixel[ii] < 1)
 3639         {
 3640             free(buffer);
 3641             return(*status = BAD_PIX_NUM);
 3642         }
 3643 
 3644         tilesize[ii] = (fptr->Fptr)->tilesize[ii];
 3645         tiledim[ii] = (naxis[ii] - 1) / tilesize[ii] + 1;
 3646         ftile[ii]   = (fpixel[ii] - 1)   / tilesize[ii] + 1;
 3647         ltile[ii]   = minvalue((lpixel[ii] - 1) / tilesize[ii] + 1, 
 3648                                 tiledim[ii]);
 3649         rowdim[ii]  = ntemp;  /* total tiles in each dimension */
 3650         ntemp *= tiledim[ii];
 3651     }
 3652 
 3653     /* support up to 6 dimensions for now */
 3654     /* tfpixel and tlpixel are the first and last image pixels */
 3655     /* along each dimension of the compression tile */
 3656     for (i5 = ftile[5]; i5 <= ltile[5]; i5++)
 3657     {
 3658      tfpixel[5] = (i5 - 1) * tilesize[5] + 1;
 3659      tlpixel[5] = minvalue(tfpixel[5] + tilesize[5] - 1, 
 3660                             naxis[5]);
 3661      thistilesize[5] = tlpixel[5] - tfpixel[5] + 1;
 3662      offset[5] = (i5 - 1) * rowdim[5];
 3663      for (i4 = ftile[4]; i4 <= ltile[4]; i4++)
 3664      {
 3665       tfpixel[4] = (i4 - 1) * tilesize[4] + 1;
 3666       tlpixel[4] = minvalue(tfpixel[4] + tilesize[4] - 1, 
 3667                             naxis[4]);
 3668       thistilesize[4] = thistilesize[5] * (tlpixel[4] - tfpixel[4] + 1);
 3669       offset[4] = (i4 - 1) * rowdim[4] + offset[5];
 3670       for (i3 = ftile[3]; i3 <= ltile[3]; i3++)
 3671       {
 3672         tfpixel[3] = (i3 - 1) * tilesize[3] + 1;
 3673         tlpixel[3] = minvalue(tfpixel[3] + tilesize[3] - 1, 
 3674                               naxis[3]);
 3675         thistilesize[3] = thistilesize[4] * (tlpixel[3] - tfpixel[3] + 1);
 3676         offset[3] = (i3 - 1) * rowdim[3] + offset[4];
 3677         for (i2 = ftile[2]; i2 <= ltile[2]; i2++)
 3678         {
 3679           tfpixel[2] = (i2 - 1) * tilesize[2] + 1;
 3680           tlpixel[2] = minvalue(tfpixel[2] + tilesize[2] - 1, 
 3681                                 naxis[2]);
 3682           thistilesize[2] = thistilesize[3] * (tlpixel[2] - tfpixel[2] + 1);
 3683           offset[2] = (i2 - 1) * rowdim[2] + offset[3];
 3684           for (i1 = ftile[1]; i1 <= ltile[1]; i1++)
 3685           {
 3686             tfpixel[1] = (i1 - 1) * tilesize[1] + 1;
 3687             tlpixel[1] = minvalue(tfpixel[1] + tilesize[1] - 1, 
 3688                                   naxis[1]);
 3689             thistilesize[1] = thistilesize[2] * (tlpixel[1] - tfpixel[1] + 1);
 3690             offset[1] = (i1 - 1) * rowdim[1] + offset[2];
 3691             for (i0 = ftile[0]; i0 <= ltile[0]; i0++)
 3692             {
 3693               tfpixel[0] = (i0 - 1) * tilesize[0] + 1;
 3694               tlpixel[0] = minvalue(tfpixel[0] + tilesize[0] - 1, 
 3695                                     naxis[0]);
 3696               thistilesize[0] = thistilesize[1] * (tlpixel[0] - tfpixel[0] + 1);
 3697               /* calculate row of table containing this tile */
 3698               irow = i0 + offset[1];
 3699 
 3700               /* read and uncompress this row (tile) of the table */
 3701               /* also do type conversion and undefined pixel substitution */
 3702               /* at this point */
 3703               imcomp_decompress_tile(fptr, irow, thistilesize[0],
 3704                     datatype, nullcheck, nullval, buffer, bnullarray, &tilenul,
 3705                      status);
 3706 
 3707               if (*status == NO_COMPRESSED_TILE)
 3708               {
 3709                    /* tile doesn't exist, so initialize to zero */
 3710                    memset(buffer, 0, pixlen * thistilesize[0]);
 3711                    *status = 0;
 3712               }
 3713 
 3714               /* copy the intersecting pixels to this tile from the input */
 3715               imcomp_merge_overlap(buffer, pixlen, ndim, tfpixel, tlpixel, 
 3716                      bnullarray, array, fpixel, lpixel, nullcheck, status);
 3717                      
 3718              /* Collapse sizes of higher dimension tiles into 2 dimensional
 3719                 equivalents needed by the quantizing algorithms for
 3720                 floating point types */
 3721               fits_calc_tile_rows(tlpixel, tfpixel, ndim, &trowsize,
 3722                               &ntrows, status);
 3723 
 3724               /* compress the tile again, and write it back to the FITS file */
 3725               imcomp_compress_tile (fptr, irow, datatype, buffer, 
 3726                                     thistilesize[0],
 3727                     trowsize,
 3728                     ntrows,
 3729                     nullcheck, nullval, 
 3730                     status);
 3731             }
 3732           }
 3733         }
 3734       }
 3735      }
 3736     }
 3737     free(buffer);
 3738     
 3739 
 3740     if ((fptr->Fptr)->zbitpix < 0 && nullcheck != 0) { 
 3741 /*
 3742      This is a floating point FITS image with possible null values.
 3743      It is too messy to test if any null values are actually written, so 
 3744      just assume so.  We need to make sure that the
 3745      ZBLANK keyword is present in the compressed image header.  If it is not
 3746      there then we need to insert the keyword. 
 3747 */   
 3748         tstatus = 0;
 3749         ffgcrd(fptr, "ZBLANK", card, &tstatus);
 3750 
 3751     if (tstatus) {   /* have to insert the ZBLANK keyword */
 3752            ffgcrd(fptr, "ZCMPTYPE", card, status);
 3753            ffikyj(fptr, "ZBLANK", COMPRESS_NULL_VALUE, 
 3754                 "null value in the compressed integer array", status);
 3755     
 3756            /* set this value into the internal structure; it is used if */
 3757        /* the program reads back the values from the array */
 3758      
 3759           (fptr->Fptr)->zblank = COMPRESS_NULL_VALUE;
 3760           (fptr->Fptr)->cn_zblank = -1;  /* flag for a constant ZBLANK */
 3761         }  
 3762     }  
 3763     
 3764     return(*status);
 3765 }
 3766 /*--------------------------------------------------------------------------*/
 3767 int fits_write_compressed_pixels(fitsfile *fptr, /* I - FITS file pointer   */
 3768             int  datatype,  /* I - datatype of the array to be written      */
 3769             LONGLONG   fpixel,  /* I - 'first pixel to write          */
 3770             LONGLONG   npixel,  /* I - number of pixels to write      */
 3771             int  nullcheck,  /* I - 0 for no null checking                   */
 3772                              /*     1: pixels that are = nullval will be     */
 3773                              /*     written with the FITS null pixel value   */
 3774                              /*     (floating point arrays only)             */
 3775             void *array,      /* I - array of values to write                */
 3776             void *nullval,    /* I - value used to represent undefined pixels*/
 3777             int  *status)     /* IO - error status                           */
 3778 /*
 3779    Write a consecutive set of pixels to a compressed image.  This routine
 3780    interpretes the n-dimensional image as a long one-dimensional array. 
 3781    This is actually a rather inconvenient way to write compressed images in
 3782    general, and could be rather inefficient if the requested pixels to be
 3783    written are located in many different image compression tiles.    
 3784 
 3785    The general strategy used here is to write the requested pixels in blocks
 3786    that correspond to rectangular image sections.  
 3787 */
 3788 {
 3789     int naxis, ii, bytesperpixel;
 3790     long naxes[MAX_COMPRESS_DIM], nread;
 3791     LONGLONG tfirst, tlast, last0, last1, dimsize[MAX_COMPRESS_DIM];
 3792     long nplane, firstcoord[MAX_COMPRESS_DIM], lastcoord[MAX_COMPRESS_DIM];
 3793     char *arrayptr;
 3794 
 3795     if (*status > 0)
 3796         return(*status);
 3797 
 3798     arrayptr = (char *) array;
 3799 
 3800     /* get size of array pixels, in bytes */
 3801     bytesperpixel = ffpxsz(datatype);
 3802 
 3803     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 3804     {
 3805         naxes[ii] = 1;
 3806         firstcoord[ii] = 0;
 3807         lastcoord[ii] = 0;
 3808     }
 3809 
 3810     /*  determine the dimensions of the image to be written */
 3811     ffgidm(fptr, &naxis, status);
 3812     ffgisz(fptr, MAX_COMPRESS_DIM, naxes, status);
 3813 
 3814     /* calc the cumulative number of pixels in each successive dimension */
 3815     dimsize[0] = 1;
 3816     for (ii = 1; ii < MAX_COMPRESS_DIM; ii++)
 3817          dimsize[ii] = dimsize[ii - 1] * naxes[ii - 1];
 3818 
 3819     /*  determine the coordinate of the first and last pixel in the image */
 3820     /*  Use zero based indexes here */
 3821     tfirst = fpixel - 1;
 3822     tlast = tfirst + npixel - 1;
 3823     for (ii = naxis - 1; ii >= 0; ii--)
 3824     {
 3825         firstcoord[ii] = (long) (tfirst / dimsize[ii]);
 3826         lastcoord[ii]  = (long) (tlast / dimsize[ii]);
 3827         tfirst = tfirst - firstcoord[ii] * dimsize[ii];
 3828         tlast = tlast - lastcoord[ii] * dimsize[ii];
 3829     }
 3830 
 3831     /* to simplify things, treat 1-D, 2-D, and 3-D images as separate cases */
 3832 
 3833     if (naxis == 1)
 3834     {
 3835         /* Simple: just write the requested range of pixels */
 3836 
 3837         firstcoord[0] = firstcoord[0] + 1;
 3838         lastcoord[0] = lastcoord[0] + 1;
 3839         fits_write_compressed_img(fptr, datatype, firstcoord, lastcoord,
 3840             nullcheck, array, nullval, status);
 3841         return(*status);
 3842     }
 3843     else if (naxis == 2)
 3844     {
 3845         nplane = 0;  /* write 1st (and only) plane of the image */
 3846         fits_write_compressed_img_plane(fptr, datatype, bytesperpixel,
 3847           nplane, firstcoord, lastcoord, naxes, nullcheck,
 3848           array, nullval, &nread, status);
 3849     }
 3850     else if (naxis == 3)
 3851     {
 3852         /* test for special case: writing an integral number of planes */
 3853         if (firstcoord[0] == 0 && firstcoord[1] == 0 &&
 3854             lastcoord[0] == naxes[0] - 1 && lastcoord[1] == naxes[1] - 1)
 3855         {
 3856             for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 3857             {
 3858                 /* convert from zero base to 1 base */
 3859                 (firstcoord[ii])++;
 3860                 (lastcoord[ii])++;
 3861             }
 3862 
 3863             /* we can write the contiguous block of pixels in one go */
 3864             fits_write_compressed_img(fptr, datatype, firstcoord, lastcoord,
 3865                 nullcheck, array, nullval, status);
 3866             return(*status);
 3867         }
 3868 
 3869         /* save last coordinate in temporary variables */
 3870         last0 = lastcoord[0];
 3871         last1 = lastcoord[1];
 3872 
 3873         if (firstcoord[2] < lastcoord[2])
 3874         {
 3875             /* we will write up to the last pixel in all but the last plane */
 3876             lastcoord[0] = naxes[0] - 1;
 3877             lastcoord[1] = naxes[1] - 1;
 3878         }
 3879 
 3880         /* write one plane of the cube at a time, for simplicity */
 3881         for (nplane = firstcoord[2]; nplane <= lastcoord[2]; nplane++)
 3882         {
 3883             if (nplane == lastcoord[2])
 3884             {
 3885                 lastcoord[0] = (long) last0;
 3886                 lastcoord[1] = (long) last1;
 3887             }
 3888 
 3889             fits_write_compressed_img_plane(fptr, datatype, bytesperpixel,
 3890               nplane, firstcoord, lastcoord, naxes, nullcheck,
 3891               arrayptr, nullval, &nread, status);
 3892 
 3893             /* for all subsequent planes, we start with the first pixel */
 3894             firstcoord[0] = 0;
 3895             firstcoord[1] = 0;
 3896 
 3897             /* increment pointers to next elements to be written */
 3898             arrayptr = arrayptr + nread * bytesperpixel;
 3899         }
 3900     }
 3901     else
 3902     {
 3903         ffpmsg("only 1D, 2D, or 3D images are currently supported");
 3904         return(*status = DATA_COMPRESSION_ERR);
 3905     }
 3906 
 3907     return(*status);
 3908 }
 3909 /*--------------------------------------------------------------------------*/
 3910 int fits_write_compressed_img_plane(fitsfile *fptr, /* I - FITS file    */
 3911             int  datatype,  /* I - datatype of the array to be written    */
 3912             int  bytesperpixel, /* I - number of bytes per pixel in array */
 3913             long   nplane,  /* I - which plane of the cube to write      */
 3914             long *firstcoord, /* I coordinate of first pixel to write */
 3915             long *lastcoord,  /* I coordinate of last pixel to write */
 3916             long *naxes,     /* I size of each image dimension */
 3917             int  nullcheck,  /* I - 0 for no null checking                   */
 3918                              /*     1: pixels that are = nullval will be     */
 3919                              /*     written with the FITS null pixel value   */
 3920                              /*     (floating point arrays only)             */
 3921             void *array,      /* I - array of values that are written        */
 3922             void *nullval,    /* I - value for undefined pixels              */
 3923             long *nread,      /* O - total number of pixels written          */
 3924             int  *status)     /* IO - error status                           */
 3925 
 3926    /*
 3927            in general we have to write the first partial row of the image,
 3928            followed by the middle complete rows, followed by the last
 3929            partial row of the image.  If the first or last rows are complete,
 3930            then write them at the same time as all the middle rows.
 3931     */
 3932 {
 3933     /* bottom left coord. and top right coord. */
 3934     long blc[MAX_COMPRESS_DIM], trc[MAX_COMPRESS_DIM]; 
 3935     char *arrayptr;
 3936 
 3937     *nread = 0;
 3938 
 3939     arrayptr = (char *) array;
 3940 
 3941     blc[2] = nplane + 1;
 3942     trc[2] = nplane + 1;
 3943 
 3944     if (firstcoord[0] != 0)
 3945     { 
 3946             /* have to read a partial first row */
 3947             blc[0] = firstcoord[0] + 1;
 3948             blc[1] = firstcoord[1] + 1;
 3949             trc[1] = blc[1];  
 3950             if (lastcoord[1] == firstcoord[1])
 3951                trc[0] = lastcoord[0] + 1; /* 1st and last pixels in same row */
 3952             else
 3953                trc[0] = naxes[0];  /* read entire rest of the row */
 3954 
 3955             fits_write_compressed_img(fptr, datatype, blc, trc,
 3956                 nullcheck, arrayptr, nullval, status);
 3957 
 3958             *nread = *nread + trc[0] - blc[0] + 1;
 3959 
 3960             if (lastcoord[1] == firstcoord[1])
 3961             {
 3962                return(*status);  /* finished */
 3963             }
 3964 
 3965             /* set starting coord to beginning of next line */
 3966             firstcoord[0] = 0;
 3967             firstcoord[1] += 1;
 3968             arrayptr = arrayptr + (trc[0] - blc[0] + 1) * bytesperpixel;
 3969     }
 3970 
 3971     /* write contiguous complete rows of the image, if any */
 3972     blc[0] = 1;
 3973     blc[1] = firstcoord[1] + 1;
 3974     trc[0] = naxes[0];
 3975 
 3976     if (lastcoord[0] + 1 == naxes[0])
 3977     {
 3978             /* can write the last complete row, too */
 3979             trc[1] = lastcoord[1] + 1;
 3980     }
 3981     else
 3982     {
 3983             /* last row is incomplete; have to read it separately */
 3984             trc[1] = lastcoord[1];
 3985     }
 3986 
 3987     if (trc[1] >= blc[1])  /* must have at least one whole line to read */
 3988     {
 3989         fits_write_compressed_img(fptr, datatype, blc, trc,
 3990                 nullcheck, arrayptr, nullval, status);
 3991 
 3992         *nread = *nread + (trc[1] - blc[1] + 1) * naxes[0];
 3993 
 3994         if (lastcoord[1] + 1 == trc[1])
 3995                return(*status);  /* finished */
 3996 
 3997         /* increment pointers for the last partial row */
 3998         arrayptr = arrayptr + (trc[1] - blc[1] + 1) * naxes[0] * bytesperpixel;
 3999 
 4000      }
 4001 
 4002     if (trc[1] == lastcoord[1] + 1)
 4003         return(*status);           /* all done */
 4004 
 4005     /* set starting and ending coord to last line */
 4006 
 4007     trc[0] = lastcoord[0] + 1;
 4008     trc[1] = lastcoord[1] + 1;
 4009     blc[1] = trc[1];
 4010 
 4011     fits_write_compressed_img(fptr, datatype, blc, trc,
 4012                 nullcheck, arrayptr, nullval, status);
 4013 
 4014     *nread = *nread + trc[0] - blc[0] + 1;
 4015 
 4016     return(*status);
 4017 }
 4018 
 4019 /* ######################################################################## */
 4020 /* ###                 Image Decompression Routines                     ### */
 4021 /* ######################################################################## */
 4022 
 4023 /*--------------------------------------------------------------------------*/
 4024 int fits_img_decompress (fitsfile *infptr, /* image (bintable) to uncompress */
 4025               fitsfile *outfptr,   /* empty HDU for output uncompressed image */
 4026               int *status)         /* IO - error status               */
 4027 
 4028 /* 
 4029   This routine decompresses the whole image and writes it to the output file.
 4030 */
 4031 
 4032 {
 4033     int ii, datatype = 0;
 4034     int nullcheck, anynul;
 4035     LONGLONG fpixel[MAX_COMPRESS_DIM], lpixel[MAX_COMPRESS_DIM];
 4036     long inc[MAX_COMPRESS_DIM];
 4037     long imgsize;
 4038     float *nulladdr, fnulval;
 4039     double dnulval;
 4040 
 4041     if (fits_img_decompress_header(infptr, outfptr, status) > 0)
 4042     {
 4043         return (*status);
 4044     }
 4045 
 4046     /* force a rescan of the output header keywords, then reset the scaling */
 4047     /* in case the BSCALE and BZERO keywords are present, so that the       */
 4048     /* decompressed values won't be scaled when written to the output image */
 4049     ffrdef(outfptr, status);
 4050     ffpscl(outfptr, 1.0, 0.0, status);
 4051     ffpscl(infptr, 1.0, 0.0, status);
 4052 
 4053     /* initialize; no null checking is needed for integer images */
 4054     nullcheck = 0;
 4055     nulladdr =  &fnulval;
 4056 
 4057     /* determine datatype for image */
 4058     if ((infptr->Fptr)->zbitpix == BYTE_IMG)
 4059     {
 4060         datatype = TBYTE;
 4061     }
 4062     else if ((infptr->Fptr)->zbitpix == SHORT_IMG)
 4063     {
 4064         datatype = TSHORT;
 4065     }
 4066     else if ((infptr->Fptr)->zbitpix == LONG_IMG)
 4067     {
 4068         datatype = TINT;
 4069     }
 4070     else if ((infptr->Fptr)->zbitpix == FLOAT_IMG)
 4071     {
 4072         /* In the case of float images we must check for NaNs  */
 4073         nullcheck = 1;
 4074         fnulval = FLOATNULLVALUE;
 4075         nulladdr =  &fnulval;
 4076         datatype = TFLOAT;
 4077     }
 4078     else if ((infptr->Fptr)->zbitpix == DOUBLE_IMG)
 4079     {
 4080         /* In the case of double images we must check for NaNs  */
 4081         nullcheck = 1;
 4082         dnulval = DOUBLENULLVALUE;
 4083         nulladdr = (float *) &dnulval;
 4084         datatype = TDOUBLE;
 4085     }
 4086 
 4087     /* calculate size of the image (in pixels) */
 4088     imgsize = 1;
 4089     for (ii = 0; ii < (infptr->Fptr)->zndim; ii++)
 4090     {
 4091         imgsize *= (infptr->Fptr)->znaxis[ii];
 4092         fpixel[ii] = 1;              /* Set first and last pixel to */
 4093         lpixel[ii] = (infptr->Fptr)->znaxis[ii]; /* include the entire image. */
 4094         inc[ii] = 1;
 4095     }
 4096 
 4097     /* uncompress the input image and write to output image, one tile at a time */
 4098 
 4099     fits_read_write_compressed_img(infptr, datatype, fpixel, lpixel, inc,  
 4100             nullcheck, nulladdr, &anynul, outfptr, status);
 4101 
 4102     return (*status);
 4103 }
 4104 /*--------------------------------------------------------------------------*/
 4105 int fits_decompress_img (fitsfile *infptr, /* image (bintable) to uncompress */
 4106               fitsfile *outfptr,   /* empty HDU for output uncompressed image */
 4107               int *status)         /* IO - error status               */
 4108 
 4109 /* 
 4110   THIS IS AN OBSOLETE ROUTINE.  USE fits_img_decompress instead!!!
 4111   
 4112   This routine decompresses the whole image and writes it to the output file.
 4113 */
 4114 
 4115 {
 4116     double *data;
 4117     int ii, datatype = 0, byte_per_pix = 0;
 4118     int nullcheck, anynul;
 4119     LONGLONG fpixel[MAX_COMPRESS_DIM], lpixel[MAX_COMPRESS_DIM];
 4120     long inc[MAX_COMPRESS_DIM];
 4121     long imgsize, memsize;
 4122     float *nulladdr, fnulval;
 4123     double dnulval;
 4124 
 4125     if (*status > 0)
 4126         return(*status);
 4127 
 4128     if (!fits_is_compressed_image(infptr, status) )
 4129     {
 4130         ffpmsg("CHDU is not a compressed image (fits_decompress_img)");
 4131         return(*status = DATA_DECOMPRESSION_ERR);
 4132     }
 4133 
 4134     /* create an empty output image with the correct dimensions */
 4135     if (ffcrim(outfptr, (infptr->Fptr)->zbitpix, (infptr->Fptr)->zndim, 
 4136        (infptr->Fptr)->znaxis, status) > 0)
 4137     {
 4138         ffpmsg("error creating output decompressed image HDU");
 4139         return (*status);
 4140     }
 4141     /* Copy the table header to the image header. */
 4142     if (imcomp_copy_imheader(infptr, outfptr, status) > 0)
 4143     {
 4144         ffpmsg("error copying header of compressed image");
 4145         return (*status);
 4146     }
 4147 
 4148     /* force a rescan of the output header keywords, then reset the scaling */
 4149     /* in case the BSCALE and BZERO keywords are present, so that the       */
 4150     /* decompressed values won't be scaled when written to the output image */
 4151     ffrdef(outfptr, status);
 4152     ffpscl(outfptr, 1.0, 0.0, status);
 4153     ffpscl(infptr, 1.0, 0.0, status);
 4154 
 4155     /* initialize; no null checking is needed for integer images */
 4156     nullcheck = 0;
 4157     nulladdr =  &fnulval;
 4158 
 4159     /* determine datatype for image */
 4160     if ((infptr->Fptr)->zbitpix == BYTE_IMG)
 4161     {
 4162         datatype = TBYTE;
 4163         byte_per_pix = 1;
 4164     }
 4165     else if ((infptr->Fptr)->zbitpix == SHORT_IMG)
 4166     {
 4167         datatype = TSHORT;
 4168         byte_per_pix = sizeof(short);
 4169     }
 4170     else if ((infptr->Fptr)->zbitpix == LONG_IMG)
 4171     {
 4172         datatype = TINT;
 4173         byte_per_pix = sizeof(int);
 4174     }
 4175     else if ((infptr->Fptr)->zbitpix == FLOAT_IMG)
 4176     {
 4177         /* In the case of float images we must check for NaNs  */
 4178         nullcheck = 1;
 4179         fnulval = FLOATNULLVALUE;
 4180         nulladdr =  &fnulval;
 4181         datatype = TFLOAT;
 4182         byte_per_pix = sizeof(float);
 4183     }
 4184     else if ((infptr->Fptr)->zbitpix == DOUBLE_IMG)
 4185     {
 4186         /* In the case of double images we must check for NaNs  */
 4187         nullcheck = 1;
 4188         dnulval = DOUBLENULLVALUE;
 4189         nulladdr = (float *) &dnulval;
 4190         datatype = TDOUBLE;
 4191         byte_per_pix = sizeof(double);
 4192     }
 4193 
 4194     /* calculate size of the image (in pixels) */
 4195     imgsize = 1;
 4196     for (ii = 0; ii < (infptr->Fptr)->zndim; ii++)
 4197     {
 4198         imgsize *= (infptr->Fptr)->znaxis[ii];
 4199         fpixel[ii] = 1;              /* Set first and last pixel to */
 4200         lpixel[ii] = (infptr->Fptr)->znaxis[ii]; /* include the entire image. */
 4201         inc[ii] = 1;
 4202     }
 4203     /* Calc equivalent number of double pixels same size as whole the image. */
 4204     /* We use double datatype to force the memory to be aligned properly */
 4205     memsize = ((imgsize * byte_per_pix) - 1) / sizeof(double) + 1;
 4206 
 4207     /* allocate memory for the image */
 4208     data = (double*) calloc (memsize, sizeof(double));
 4209     if (!data)
 4210     { 
 4211         ffpmsg("Couldn't allocate memory for the uncompressed image");
 4212         return(*status = MEMORY_ALLOCATION);
 4213     }
 4214 
 4215     /* uncompress the entire image into memory */
 4216     /* This routine should be enhanced sometime to only need enough */
 4217     /* memory to uncompress one tile at a time.  */
 4218     fits_read_compressed_img(infptr, datatype, fpixel, lpixel, inc,  
 4219             nullcheck, nulladdr, data, NULL, &anynul, status);
 4220 
 4221     /* write the image to the output file */
 4222     if (anynul)
 4223         fits_write_imgnull(outfptr, datatype, 1, imgsize, data, nulladdr, 
 4224                           status);
 4225     else
 4226         fits_write_img(outfptr, datatype, 1, imgsize, data, status);
 4227 
 4228     free(data);
 4229     return (*status);
 4230 }
 4231 /*--------------------------------------------------------------------------*/
 4232 int fits_img_decompress_header(fitsfile *infptr, /* image (bintable) to uncompress */
 4233               fitsfile *outfptr,   /* empty HDU for output uncompressed image */
 4234               int *status)         /* IO - error status               */
 4235 
 4236 /* 
 4237   This routine reads the header of the input tile compressed image and 
 4238   converts it to that of a standard uncompress FITS image.
 4239 */
 4240 
 4241 {
 4242     int writeprime = 0;
 4243     int hdupos, inhdupos, numkeys;
 4244     int nullprime = 0, copyprime = 0, norec = 0, tstatus;
 4245     char card[FLEN_CARD];
 4246     int ii, naxis, bitpix;
 4247     long naxes[MAX_COMPRESS_DIM];
 4248 
 4249     if (*status > 0)
 4250         return(*status);
 4251     else if (*status == -1) {
 4252         *status = 0;
 4253     writeprime = 1;
 4254     }
 4255 
 4256     if (!fits_is_compressed_image(infptr, status) )
 4257     {
 4258         ffpmsg("CHDU is not a compressed image (fits_img_decompress)");
 4259         return(*status = DATA_DECOMPRESSION_ERR);
 4260     }
 4261 
 4262     /* get information about the state of the output file; does it already */
 4263     /* contain any keywords and HDUs?  */
 4264     fits_get_hdu_num(infptr, &inhdupos);  /* Get the current output HDU position */
 4265     fits_get_hdu_num(outfptr, &hdupos);  /* Get the current output HDU position */
 4266     fits_get_hdrspace(outfptr, &numkeys, 0, status);
 4267 
 4268     /* Was the input compressed HDU originally the primary array image? */
 4269     tstatus = 0;
 4270     if (!fits_read_card(infptr, "ZSIMPLE", card, &tstatus)) { 
 4271       /* yes, input HDU was a primary array (not an IMAGE extension) */
 4272       /* Now determine if we can uncompress it into the primary array of */
 4273       /* the output file.  This is only possible if the output file */
 4274       /* currently only contains a null primary array, with no addition */
 4275       /* header keywords and with no following extension in the FITS file. */
 4276       
 4277       if (hdupos == 1) {  /* are we positioned at the primary array? */
 4278             if (numkeys == 0) { /* primary HDU is completely empty */
 4279             nullprime = 1;
 4280             } else {
 4281                 fits_get_img_param(outfptr, MAX_COMPRESS_DIM, &bitpix, &naxis, naxes, status);
 4282     
 4283             if (naxis == 0) { /* is this a null image? */
 4284                    nullprime = 1;
 4285 
 4286            if (inhdupos == 2)  /* must be at the first extension */
 4287               copyprime = 1;
 4288         }
 4289            }
 4290       }
 4291     } 
 4292 
 4293     if (nullprime) {  
 4294        /* We will delete the existing keywords in the null primary array
 4295           and uncompress the input image into the primary array of the output.
 4296       Some of these keywords may be added back to the uncompressed image
 4297       header later.
 4298        */
 4299 
 4300        for (ii = numkeys; ii > 0; ii--)
 4301           fits_delete_record(outfptr, ii, status);
 4302 
 4303     } else  {
 4304 
 4305        /* if the ZTENSION keyword doesn't exist, then we have to 
 4306           write the required keywords manually */
 4307        tstatus = 0;
 4308        if (fits_read_card(infptr, "ZTENSION", card, &tstatus)) {
 4309 
 4310           /* create an empty output image with the correct dimensions */
 4311           if (ffcrim(outfptr, (infptr->Fptr)->zbitpix, (infptr->Fptr)->zndim, 
 4312              (infptr->Fptr)->znaxis, status) > 0)
 4313           {
 4314              ffpmsg("error creating output decompressed image HDU");
 4315              return (*status);
 4316           }
 4317 
 4318       norec = 1;  /* the required keywords have already been written */
 4319 
 4320        } else {  /* the input compressed image does have ZTENSION keyword */
 4321        
 4322           if (writeprime) {  /* convert the image extension to a primary array */
 4323           /* have to write the required keywords manually */
 4324 
 4325               /* create an empty output image with the correct dimensions */
 4326               if (ffcrim(outfptr, (infptr->Fptr)->zbitpix, (infptr->Fptr)->zndim, 
 4327                  (infptr->Fptr)->znaxis, status) > 0)
 4328               {
 4329                  ffpmsg("error creating output decompressed image HDU");
 4330                  return (*status);
 4331               }
 4332 
 4333           norec = 1;  /* the required keywords have already been written */
 4334 
 4335           } else {  /* write the input compressed image to an image extension */
 4336 
 4337               if (numkeys == 0) {  /* the output file is currently completely empty */
 4338       
 4339              /* In this case, the input is a compressed IMAGE extension. */
 4340              /* Since the uncompressed output file is currently completely empty, */
 4341              /* we need to write a null primary array before uncompressing the */
 4342                  /* image extension */
 4343          
 4344                  ffcrim(outfptr, 8, 0, naxes, status); /* naxes is not used */
 4345          
 4346              /* now create the empty extension to uncompress into */
 4347                  if (fits_create_hdu(outfptr, status) > 0)
 4348                  {
 4349                       ffpmsg("error creating output decompressed image HDU");
 4350                       return (*status);
 4351                  }
 4352       
 4353           } else {
 4354                   /* just create a new empty extension, then copy all the required */
 4355               /* keywords into it.  */
 4356                  fits_create_hdu(outfptr, status);
 4357           }
 4358            }
 4359        }
 4360 
 4361     }
 4362 
 4363     if (*status > 0)  {
 4364         ffpmsg("error creating output decompressed image HDU");
 4365         return (*status);
 4366     }
 4367 
 4368     /* Copy the table header to the image header. */
 4369 
 4370     if (imcomp_copy_comp2img(infptr, outfptr, norec, status) > 0)
 4371     {
 4372         ffpmsg("error copying header keywords from compressed image");
 4373     }
 4374 
 4375     if (copyprime) {  
 4376     /* append any unexpected keywords from the primary array.
 4377        This includes any keywords except SIMPLE, BITPIX, NAXIS,
 4378        EXTEND, COMMENT, HISTORY, CHECKSUM, and DATASUM.
 4379     */
 4380 
 4381         fits_movabs_hdu(infptr, 1, NULL, status);  /* move to primary array */
 4382     
 4383         /* do this so that any new keywords get written before any blank
 4384        keywords that may have been appended by imcomp_copy_comp2img  */
 4385         fits_set_hdustruc(outfptr, status);
 4386 
 4387         if (imcomp_copy_prime2img(infptr, outfptr, status) > 0)
 4388         {
 4389             ffpmsg("error copying primary keywords from compressed file");
 4390         }
 4391 
 4392         fits_movabs_hdu(infptr, 2, NULL, status); /* move back to where we were */
 4393     }
 4394 
 4395     return (*status);
 4396 }
 4397 /*---------------------------------------------------------------------------*/
 4398 int fits_read_compressed_img(fitsfile *fptr,   /* I - FITS file pointer      */
 4399             int  datatype,  /* I - datatype of the array to be returned      */
 4400             LONGLONG  *infpixel, /* I - 'bottom left corner' of the subsection    */
 4401             LONGLONG  *inlpixel, /* I - 'top right corner' of the subsection      */
 4402             long  *ininc,    /* I - increment to be applied in each dimension */
 4403             int  nullcheck,  /* I - 0 for no null checking                   */
 4404                               /*     1: set undefined pixels = nullval       */
 4405                               /*     2: set nullarray=1 for undefined pixels */
 4406             void *nullval,    /* I - value for undefined pixels              */
 4407             void *array,      /* O - array of values that are returned       */
 4408             char *nullarray,  /* O - array of flags = 1 if nullcheck = 2     */
 4409             int  *anynul,     /* O - set to 1 if any values are null; else 0 */
 4410             int  *status)     /* IO - error status                           */
 4411 /*
 4412    Read a section of a compressed image;  Note: lpixel may be larger than the 
 4413    size of the uncompressed image.  Only the pixels within the image will be
 4414    returned.
 4415 */
 4416 {
 4417     long naxis[MAX_COMPRESS_DIM], tiledim[MAX_COMPRESS_DIM];
 4418     long tilesize[MAX_COMPRESS_DIM], thistilesize[MAX_COMPRESS_DIM];
 4419     long ftile[MAX_COMPRESS_DIM], ltile[MAX_COMPRESS_DIM];
 4420     long tfpixel[MAX_COMPRESS_DIM], tlpixel[MAX_COMPRESS_DIM];
 4421     long rowdim[MAX_COMPRESS_DIM], offset[MAX_COMPRESS_DIM],ntemp;
 4422     long fpixel[MAX_COMPRESS_DIM], lpixel[MAX_COMPRESS_DIM];
 4423     long inc[MAX_COMPRESS_DIM];
 4424     long i5, i4, i3, i2, i1, i0, irow;
 4425     int ii, ndim, pixlen, tilenul=0;
 4426     void *buffer;
 4427     char *bnullarray = 0;
 4428     double testnullval = 0.;
 4429 
 4430     if (*status > 0) 
 4431         return(*status);
 4432 
 4433     if (!fits_is_compressed_image(fptr, status) )
 4434     {
 4435         ffpmsg("CHDU is not a compressed image (fits_read_compressed_img)");
 4436         return(*status = DATA_DECOMPRESSION_ERR);
 4437     }
 4438 
 4439     /* get temporary space for uncompressing one image tile */
 4440     if (datatype == TSHORT)
 4441     {
 4442        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (short)); 
 4443        pixlen = sizeof(short);
 4444        if (nullval)
 4445            testnullval = *(short *) nullval;
 4446     }
 4447     else if (datatype == TINT)
 4448     {
 4449        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (int));
 4450        pixlen = sizeof(int);
 4451        if (nullval)
 4452            testnullval = *(int *) nullval;
 4453     }
 4454     else if (datatype == TLONG)
 4455     {
 4456        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (long));
 4457        pixlen = sizeof(long);
 4458        if (nullval)
 4459            testnullval = *(long *) nullval;
 4460     }
 4461     else if (datatype == TFLOAT)
 4462     {
 4463        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (float));
 4464        pixlen = sizeof(float);
 4465        if (nullval)
 4466            testnullval = *(float *) nullval;
 4467     }
 4468     else if (datatype == TDOUBLE)
 4469     {
 4470        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (double));
 4471        pixlen = sizeof(double);
 4472        if (nullval)
 4473            testnullval = *(double *) nullval;
 4474     }
 4475     else if (datatype == TUSHORT)
 4476     {
 4477        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned short));
 4478        pixlen = sizeof(short);
 4479        if (nullval)
 4480            testnullval = *(unsigned short *) nullval;
 4481     }
 4482     else if (datatype == TUINT)
 4483     {
 4484        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned int));
 4485        pixlen = sizeof(int);
 4486        if (nullval)
 4487            testnullval = *(unsigned int *) nullval;
 4488     }
 4489     else if (datatype == TULONG)
 4490     {
 4491        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned long));
 4492        pixlen = sizeof(long);
 4493        if (nullval)
 4494            testnullval = *(unsigned long *) nullval;
 4495     }
 4496     else if (datatype == TBYTE || datatype == TSBYTE)
 4497     {
 4498        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (char));
 4499        pixlen = 1;
 4500        if (nullval)
 4501            testnullval = *(unsigned char *) nullval;
 4502     }
 4503     else
 4504     {
 4505         ffpmsg("unsupported datatype for uncompressing image");
 4506         return(*status = BAD_DATATYPE);
 4507     }
 4508 
 4509     /* If nullcheck ==1 and nullval == 0, then this means that the */
 4510     /* calling routine does not want to check for null pixels in the array */
 4511     if (nullcheck == 1 && testnullval == 0.)
 4512         nullcheck = 0;
 4513 
 4514     if (buffer == NULL)
 4515     {
 4516         ffpmsg("Out of memory (fits_read_compress_img)");
 4517         return (*status = MEMORY_ALLOCATION);
 4518     }
 4519     
 4520     /* allocate memory for a null flag array, if needed */
 4521     if (nullcheck == 2)
 4522     {
 4523         bnullarray = calloc ((fptr->Fptr)->maxtilelen, sizeof (char));
 4524 
 4525         if (bnullarray == NULL)
 4526         {
 4527         ffpmsg("Out of memory (fits_read_compress_img)");
 4528             free(buffer);
 4529         return (*status = MEMORY_ALLOCATION);
 4530         }
 4531     }
 4532 
 4533     /* initialize all the arrays */
 4534     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 4535     {
 4536         naxis[ii] = 1;
 4537         tiledim[ii] = 1;
 4538         tilesize[ii] = 1;
 4539         ftile[ii] = 1;
 4540         ltile[ii] = 1;
 4541         rowdim[ii] = 1;
 4542     }
 4543 
 4544     ndim = (fptr->Fptr)->zndim;
 4545     ntemp = 1;
 4546     for (ii = 0; ii < ndim; ii++)
 4547     {
 4548         /* support for mirror-reversed image sections */
 4549         if (infpixel[ii] <= inlpixel[ii])
 4550         {
 4551            fpixel[ii] = (long) infpixel[ii];
 4552            lpixel[ii] = (long) inlpixel[ii];
 4553            inc[ii]    = ininc[ii];
 4554         }
 4555         else
 4556         {
 4557            fpixel[ii] = (long) inlpixel[ii];
 4558            lpixel[ii] = (long) infpixel[ii];
 4559            inc[ii]    = -ininc[ii];
 4560         }
 4561 
 4562         /* calc number of tiles in each dimension, and tile containing */
 4563         /* the first and last pixel we want to read in each dimension  */
 4564         naxis[ii] = (fptr->Fptr)->znaxis[ii];
 4565         if (fpixel[ii] < 1)
 4566         {
 4567             if (nullcheck == 2)
 4568             {
 4569                 free(bnullarray);
 4570             }
 4571             free(buffer);
 4572             return(*status = BAD_PIX_NUM);
 4573         }
 4574 
 4575         tilesize[ii] = (fptr->Fptr)->tilesize[ii];
 4576         tiledim[ii] = (naxis[ii] - 1) / tilesize[ii] + 1;
 4577         ftile[ii]   = (fpixel[ii] - 1)   / tilesize[ii] + 1;
 4578         ltile[ii]   = minvalue((lpixel[ii] - 1) / tilesize[ii] + 1, 
 4579                                 tiledim[ii]);
 4580         rowdim[ii]  = ntemp;  /* total tiles in each dimension */
 4581         ntemp *= tiledim[ii];
 4582     }
 4583 
 4584     if (anynul)
 4585        *anynul = 0;  /* initialize */
 4586 
 4587     /* support up to 6 dimensions for now */
 4588     /* tfpixel and tlpixel are the first and last image pixels */
 4589     /* along each dimension of the compression tile */
 4590     for (i5 = ftile[5]; i5 <= ltile[5]; i5++)
 4591     {
 4592      tfpixel[5] = (i5 - 1) * tilesize[5] + 1;
 4593      tlpixel[5] = minvalue(tfpixel[5] + tilesize[5] - 1, 
 4594                             naxis[5]);
 4595      thistilesize[5] = tlpixel[5] - tfpixel[5] + 1;
 4596      offset[5] = (i5 - 1) * rowdim[5];
 4597      for (i4 = ftile[4]; i4 <= ltile[4]; i4++)
 4598      {
 4599       tfpixel[4] = (i4 - 1) * tilesize[4] + 1;
 4600       tlpixel[4] = minvalue(tfpixel[4] + tilesize[4] - 1, 
 4601                             naxis[4]);
 4602       thistilesize[4] = thistilesize[5] * (tlpixel[4] - tfpixel[4] + 1);
 4603       offset[4] = (i4 - 1) * rowdim[4] + offset[5];
 4604       for (i3 = ftile[3]; i3 <= ltile[3]; i3++)
 4605       {
 4606         tfpixel[3] = (i3 - 1) * tilesize[3] + 1;
 4607         tlpixel[3] = minvalue(tfpixel[3] + tilesize[3] - 1, 
 4608                               naxis[3]);
 4609         thistilesize[3] = thistilesize[4] * (tlpixel[3] - tfpixel[3] + 1);
 4610         offset[3] = (i3 - 1) * rowdim[3] + offset[4];
 4611         for (i2 = ftile[2]; i2 <= ltile[2]; i2++)
 4612         {
 4613           tfpixel[2] = (i2 - 1) * tilesize[2] + 1;
 4614           tlpixel[2] = minvalue(tfpixel[2] + tilesize[2] - 1, 
 4615                                 naxis[2]);
 4616           thistilesize[2] = thistilesize[3] * (tlpixel[2] - tfpixel[2] + 1);
 4617           offset[2] = (i2 - 1) * rowdim[2] + offset[3];
 4618           for (i1 = ftile[1]; i1 <= ltile[1]; i1++)
 4619           {
 4620             tfpixel[1] = (i1 - 1) * tilesize[1] + 1;
 4621             tlpixel[1] = minvalue(tfpixel[1] + tilesize[1] - 1, 
 4622                                   naxis[1]);
 4623             thistilesize[1] = thistilesize[2] * (tlpixel[1] - tfpixel[1] + 1);
 4624             offset[1] = (i1 - 1) * rowdim[1] + offset[2];
 4625             for (i0 = ftile[0]; i0 <= ltile[0]; i0++)
 4626             {
 4627              tfpixel[0] = (i0 - 1) * tilesize[0] + 1;
 4628              tlpixel[0] = minvalue(tfpixel[0] + tilesize[0] - 1, 
 4629                                     naxis[0]);
 4630               thistilesize[0] = thistilesize[1] * (tlpixel[0] - tfpixel[0] + 1);
 4631               /* calculate row of table containing this tile */
 4632               irow = i0 + offset[1];
 4633 
 4634 /*
 4635 printf("row %d, %d %d, %d %d, %d %d; %d\n",
 4636               irow, tfpixel[0],tlpixel[0],tfpixel[1],tlpixel[1],tfpixel[2],tlpixel[2],
 4637           thistilesize[0]);
 4638 */   
 4639               /* test if there are any intersecting pixels in this tile and the output image */
 4640               if (imcomp_test_overlap(ndim, tfpixel, tlpixel, 
 4641                       fpixel, lpixel, inc, status)) {
 4642                   /* read and uncompress this row (tile) of the table */
 4643                   /* also do type conversion and undefined pixel substitution */
 4644                   /* at this point */
 4645 
 4646                   imcomp_decompress_tile(fptr, irow, thistilesize[0],
 4647                     datatype, nullcheck, nullval, buffer, bnullarray, &tilenul,
 4648                      status);
 4649 
 4650                   if (tilenul && anynul)
 4651                       *anynul = 1;  /* there are null pixels */
 4652 /*
 4653 printf(" pixlen=%d, ndim=%d, %d %d %d, %d %d %d, %d %d %d\n",
 4654      pixlen, ndim, fpixel[0],lpixel[0],inc[0],fpixel[1],lpixel[1],inc[1],
 4655      fpixel[2],lpixel[2],inc[2]);
 4656 */
 4657                   /* copy the intersecting pixels from this tile to the output */
 4658                   imcomp_copy_overlap(buffer, pixlen, ndim, tfpixel, tlpixel, 
 4659                      bnullarray, array, fpixel, lpixel, inc, nullcheck, 
 4660                      nullarray, status);
 4661                }
 4662             }
 4663           }
 4664         }
 4665       }
 4666      }
 4667     }
 4668     if (nullcheck == 2)
 4669     {
 4670         free(bnullarray);
 4671     }
 4672     free(buffer);
 4673 
 4674     return(*status);
 4675 }
 4676 /*---------------------------------------------------------------------------*/
 4677 int fits_read_write_compressed_img(fitsfile *fptr,   /* I - FITS file pointer      */
 4678             int  datatype,  /* I - datatype of the array to be returned      */
 4679             LONGLONG  *infpixel, /* I - 'bottom left corner' of the subsection    */
 4680             LONGLONG  *inlpixel, /* I - 'top right corner' of the subsection      */
 4681             long  *ininc,    /* I - increment to be applied in each dimension */
 4682             int  nullcheck,  /* I - 0 for no null checking                   */
 4683                               /*     1: set undefined pixels = nullval       */
 4684             void *nullval,    /* I - value for undefined pixels              */
 4685             int  *anynul,     /* O - set to 1 if any values are null; else 0 */
 4686             fitsfile *outfptr,   /* I - FITS file pointer                    */
 4687             int  *status)     /* IO - error status                           */
 4688 /*
 4689    This is similar to fits_read_compressed_img, except that it writes
 4690    the pixels to the output image, on a tile by tile basis instead of returning
 4691    the array.
 4692 */
 4693 {
 4694     long naxis[MAX_COMPRESS_DIM], tiledim[MAX_COMPRESS_DIM];
 4695     long tilesize[MAX_COMPRESS_DIM], thistilesize[MAX_COMPRESS_DIM];
 4696     long ftile[MAX_COMPRESS_DIM], ltile[MAX_COMPRESS_DIM];
 4697     long tfpixel[MAX_COMPRESS_DIM], tlpixel[MAX_COMPRESS_DIM];
 4698     long rowdim[MAX_COMPRESS_DIM], offset[MAX_COMPRESS_DIM],ntemp;
 4699     long fpixel[MAX_COMPRESS_DIM], lpixel[MAX_COMPRESS_DIM];
 4700     long inc[MAX_COMPRESS_DIM];
 4701     long i5, i4, i3, i2, i1, i0, irow;
 4702     int ii, ndim, tilenul;
 4703     void *buffer;
 4704     char *bnullarray = 0, *cnull;
 4705     LONGLONG firstelem;
 4706 
 4707     if (*status > 0) 
 4708         return(*status);
 4709 
 4710     if (!fits_is_compressed_image(fptr, status) )
 4711     {
 4712         ffpmsg("CHDU is not a compressed image (fits_read_compressed_img)");
 4713         return(*status = DATA_DECOMPRESSION_ERR);
 4714     }
 4715 
 4716     cnull = (char *) nullval;  /* used to test if the nullval = 0 */
 4717     
 4718     /* get temporary space for uncompressing one image tile */
 4719     /* If nullval == 0, then this means that the */
 4720     /* calling routine does not want to check for null pixels in the array */
 4721     if (datatype == TSHORT)
 4722     {
 4723        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (short)); 
 4724        if (cnull) {
 4725          if (cnull[0] == 0 && cnull[1] == 0 ) {
 4726            nullcheck = 0;
 4727      }
 4728        }
 4729     }
 4730     else if (datatype == TINT)
 4731     {
 4732        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (int));
 4733        if (cnull) {
 4734          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0 ) {
 4735            nullcheck = 0;
 4736      }
 4737        }
 4738     }
 4739     else if (datatype == TLONG)
 4740     {
 4741        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (long));
 4742        if (cnull) {
 4743          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0 ) {
 4744            nullcheck = 0;
 4745      }
 4746        }
 4747     }
 4748     else if (datatype == TFLOAT)
 4749     {
 4750        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (float));
 4751        if (cnull) {
 4752          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0  ) {
 4753            nullcheck = 0;
 4754      }
 4755        }
 4756     }
 4757     else if (datatype == TDOUBLE)
 4758     {
 4759        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (double));
 4760        if (cnull) {
 4761          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0 &&
 4762          cnull[4] == 0 && cnull[5] == 0 && cnull[6] == 0 && cnull[7] == 0 ) {
 4763            nullcheck = 0;
 4764      }
 4765        }
 4766     }
 4767     else if (datatype == TUSHORT)
 4768     {
 4769        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned short));
 4770        if (cnull) {
 4771          if (cnull[0] == 0 && cnull[1] == 0 ){
 4772            nullcheck = 0;
 4773      }
 4774        }
 4775     }
 4776     else if (datatype == TUINT)
 4777     {
 4778        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned int));
 4779        if (cnull) {
 4780          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0 ){
 4781            nullcheck = 0;
 4782      }
 4783        }
 4784     }
 4785     else if (datatype == TULONG)
 4786     {
 4787        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (unsigned long));
 4788        if (cnull) {
 4789          if (cnull[0] == 0 && cnull[1] == 0 && cnull[2] == 0 && cnull[3] == 0 ){
 4790            nullcheck = 0;
 4791      }
 4792        }
 4793     }
 4794     else if (datatype == TBYTE || datatype == TSBYTE)
 4795     {
 4796        buffer =  malloc ((fptr->Fptr)->maxtilelen * sizeof (char));
 4797        if (cnull) {
 4798          if (cnull[0] == 0){
 4799            nullcheck = 0;
 4800      }
 4801        }
 4802     }
 4803     else
 4804     {
 4805         ffpmsg("unsupported datatype for uncompressing image");
 4806         return(*status = BAD_DATATYPE);
 4807     }
 4808 
 4809     if (buffer == NULL)
 4810     {
 4811         ffpmsg("Out of memory (fits_read_compress_img)");
 4812         return (*status = MEMORY_ALLOCATION);
 4813     }
 4814 
 4815     /* initialize all the arrays */
 4816     for (ii = 0; ii < MAX_COMPRESS_DIM; ii++)
 4817     {
 4818         naxis[ii] = 1;
 4819         tiledim[ii] = 1;
 4820         tilesize[ii] = 1;
 4821         ftile[ii] = 1;
 4822         ltile[ii] = 1;
 4823         rowdim[ii] = 1;
 4824     }
 4825 
 4826     ndim = (fptr->Fptr)->zndim;
 4827     ntemp = 1;
 4828     for (ii = 0; ii < ndim; ii++)
 4829     {
 4830         /* support for mirror-reversed image sections */
 4831         if (infpixel[ii] <= inlpixel[ii])
 4832         {
 4833            fpixel[ii] = (long) infpixel[ii];
 4834            lpixel[ii] = (long) inlpixel[ii];
 4835            inc[ii]    = ininc[ii];
 4836         }
 4837         else
 4838         {
 4839            fpixel[ii] = (long) inlpixel[ii];
 4840            lpixel[ii] = (long) infpixel[ii];
 4841            inc[ii]    = -ininc[ii];
 4842         }
 4843 
 4844         /* calc number of tiles in each dimension, and tile containing */
 4845         /* the first and last pixel we want to read in each dimension  */
 4846         naxis[ii] = (fptr->Fptr)->znaxis[ii];
 4847         if (fpixel[ii] < 1)
 4848         {
 4849             free(buffer);
 4850             return(*status = BAD_PIX_NUM);
 4851         }
 4852 
 4853         tilesize[ii] = (fptr->Fptr)->tilesize[ii];
 4854         tiledim[ii] = (naxis[ii] - 1) / tilesize[ii] + 1;
 4855         ftile[ii]   = (fpixel[ii] - 1)   / tilesize[ii] + 1;
 4856         ltile[ii]   = minvalue((lpixel[ii] - 1) / tilesize[ii] + 1, 
 4857                                 tiledim[ii]);
 4858         rowdim[ii]  = ntemp;  /* total tiles in each dimension */
 4859         ntemp *= tiledim[ii];
 4860     }
 4861 
 4862     if (anynul)
 4863        *anynul = 0;  /* initialize */
 4864 
 4865     firstelem = 1;
 4866 
 4867     /* support up to 6 dimensions for now */
 4868     /* tfpixel and tlpixel are the first and last image pixels */
 4869     /* along each dimension of the compression tile */
 4870     for (i5 = ftile[5]; i5 <= ltile[5]; i5++)
 4871     {
 4872      tfpixel[5] = (i5 - 1) * tilesize[5] + 1;
 4873      tlpixel[5] = minvalue(tfpixel[5] + tilesize[5] - 1, 
 4874                             naxis[5]);
 4875      thistilesize[5] = tlpixel[5] - tfpixel[5] + 1;
 4876      offset[5] = (i5 - 1) * rowdim[5];
 4877      for (i4 = ftile[4]; i4 <= ltile[4]; i4++)
 4878      {
 4879       tfpixel[4] = (i4 - 1) * tilesize[4] + 1;
 4880       tlpixel[4] = minvalue(tfpixel[4] + tilesize[4] - 1, 
 4881                             naxis[4]);
 4882       thistilesize[4] = thistilesize[5] * (tlpixel[4] - tfpixel[4] + 1);
 4883       offset[4] = (i4 - 1) * rowdim[4] + offset[5];
 4884       for (i3 = ftile[3]; i3 <= ltile[3]; i3++)
 4885       {
 4886         tfpixel[3] = (i3 - 1) * tilesize[3] + 1;
 4887         tlpixel[3] = minvalue(tfpixel[3] + tilesize[3] - 1, 
 4888                               naxis[3]);
 4889         thistilesize[3] = thistilesize[4] * (tlpixel[3] - tfpixel[3] + 1);
 4890         offset[3] = (i3 - 1) * rowdim[3] + offset[4];
 4891         for (i2 = ftile[2]; i2 <= ltile[2]; i2++)
 4892         {
 4893           tfpixel[2] = (i2 - 1) * tilesize[2] + 1;
 4894           tlpixel[2] = minvalue(tfpixel[2] + tilesize[2] - 1, 
 4895                                 naxis[2]);
 4896           thistilesize[2] = thistilesize[3] * (tlpixel[2] - tfpixel[2] + 1);
 4897           offset[2] = (i2 - 1) * rowdim[2] + offset[3];
 4898           for (i1 = ftile[1]; i1 <= ltile[1]; i1++)
 4899           {
 4900             tfpixel[1] = (i1 - 1) * tilesize[1] + 1;
 4901             tlpixel[1] = minvalue(tfpixel[1] + tilesize[1] - 1, 
 4902                                   naxis[1]);
 4903             thistilesize[1] = thistilesize[2] * (tlpixel[1] - tfpixel[1] + 1);
 4904             offset[1] = (i1 - 1) * rowdim[1] + offset[2];
 4905             for (i0 = ftile[0]; i0 <= ltile[0]; i0++)
 4906             {
 4907               tfpixel[0] = (i0 - 1) * tilesize[0] + 1;
 4908               tlpixel[0] = minvalue(tfpixel[0] + tilesize[0] - 1, 
 4909                                     naxis[0]);
 4910               thistilesize[0] = thistilesize[1] * (tlpixel[0] - tfpixel[0] + 1);
 4911               /* calculate row of table containing this tile */
 4912               irow = i0 + offset[1];
 4913  
 4914               /* read and uncompress this row (tile) of the table */
 4915               /* also do type conversion and undefined pixel substitution */
 4916               /* at this point */
 4917 
 4918               imcomp_decompress_tile(fptr, irow, thistilesize[0],
 4919                     datatype, nullcheck, nullval, buffer, bnullarray, &tilenul,
 4920                      status);
 4921 
 4922                /* write the image to the output file */
 4923 
 4924               if (tilenul && anynul) {     
 4925                    /* this assumes that the tiled pixels are in the same order
 4926               as in the uncompressed FITS image.  This is not necessarily
 4927               the case, but it almost alway is in practice.  
 4928               Note that null checking is not performed for integer images,
 4929               so this could only be a problem for tile compressed floating
 4930               point images that use an unconventional tiling pattern.
 4931            */
 4932                    fits_write_imgnull(outfptr, datatype, firstelem, thistilesize[0],
 4933               buffer, nullval, status);
 4934               } else {
 4935                   fits_write_subset(outfptr, datatype, tfpixel, tlpixel, 
 4936               buffer, status);
 4937               }
 4938 
 4939               firstelem += thistilesize[0];
 4940 
 4941             }
 4942           }
 4943         }
 4944       }
 4945      }
 4946     }
 4947 
 4948     free(buffer);
 4949 
 4950     return(*status);
 4951 }
 4952 /*--------------------------------------------------------------------------*/
 4953 int fits_read_compressed_pixels(fitsfile *fptr, /* I - FITS file pointer    */
 4954             int  datatype,  /* I - datatype of the array to be returned     */
 4955             LONGLONG   fpixel, /* I - 'first pixel to read          */
 4956             LONGLONG   npixel,  /* I - number of pixels to read      */
 4957             int  nullcheck,  /* I - 0 for no null checking                   */
 4958                               /*     1: set undefined pixels = nullval       */
 4959                               /*     2: set nullarray=1 for undefined pixels */
 4960             void *nullval,    /* I - value for undefined pixels              */
 4961             void *array,      /* O - array of values that are returned       */
 4962             char *nullarray,  /* O - array of flags = 1 if nullcheck = 2     */
 4963             int  *anynul,     /* O - set to 1 if any values are null; else 0 */
 4964             int  *status)     /* IO - error status                           */
 4965 /*
 4966    Read a consecutive set of pixels from a compressed image.  This routine
 4967    interpretes the n-dimensional image as a long one-dimensional array. 
 4968    This is actually a rather inconvenient way to read compressed images in
 4969    general, and could be rather inefficient if the requested pixels to be
 4970    read are located in many different image compression tiles.    
 4971 
 4972    The general strategy used here is to read the requested pixels in blocks
 4973    that correspond to rectangular image sections.  
 4974 */
 4975 {
 4976     int naxis, ii, bytesperpixel, planenul;
 4977     long naxes[MAX_COMPRESS_DIM], nread;
 4978     long nplane, inc[MAX_COMPRESS_DIM];
 4979     LONGLONG tfirst, tlast, last0, last1, dimsize[MAX_COMPRESS_DIM];
 4980     LONGLONG firstcoord[MAX_COMPRESS_DIM], lastcoord[MAX_COMPRESS_DIM];