IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38355


Ignore:
Timestamp:
Jun 3, 2015, 4:10:26 AM (11 years ago)
Author:
eugene
Message:

clarify image vs compression buffer vs compressed pixel sizes to allow hcompress and rice

Location:
branches/eam_branches/ohana.20150429/src/libfits
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ohana.20150429/src/libfits/include/gfitsio.h

    r38350 r38355  
    192192int     gfits_uncompress_image         PROTO((Header *header, Matrix *matrix, FTable *ftable));
    193193int     gfits_uncompress_data          PROTO((char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int out_pixsize));
    194 int     gfits_distribute_data          PROTO((Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero));
    195 
    196 int     gfits_compress_data            PROTO((char *zdata, int *Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *rawdata, int Nraw, int raw_pixsize));
     194int     gfits_distribute_data          PROTO((Matrix *matrix, char *data, int Ndata, int bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero));
     195
     196int     gfits_compress_data            PROTO((char *zdata, int *Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *rawdata, int Nraw, int raw_pixsize, int Nx, int Ny));
    197197int     gfits_compress_image           PROTO((Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype));
    198198int     gfits_collect_data             PROTO((Matrix *matrix, char *raw, int Nraw, int raw_pixsize, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero));
    … …  
    205205int     gfits_imtile_next              PROTO((Matrix *matrix, int *ztile, int *ntile, int *otile));
    206206
    207 int     gfits_byteswap_zdata           PROTO((char *zdata, int Nzdata, int bitpix));
    208207int     gfits_extension_is_compressed_image  PROTO((Header *header));
    209208int     gfits_extension_is_compressed_table  PROTO((Header *header));
     209int     gfits_byteswap_zdata            PROTO((char *zdata, int Nzdata, int pixsize));
     210int     gfits_compressed_data_pixsize   PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions));
    210211int     gfits_uncompressed_data_pixsize PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions));
    211 int     gfits_vartable_heap_pixsize    PROTO((char format));
     212int     gfits_uncompressed_data_bitpix  PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions));
     213int     gfits_vartable_heap_pixsize     PROTO((char format));
    212214
    213215int     gfits_varlength_column_add_data PROTO((FTable *table, char *data, off_t Ndata, int row, VarLengthColumn *column));

  • branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c

    r38340 r38355  
    1818// we do not try to support compression of an image with an associated table (not valid)
    1919
    20 # define ESCAPE(A) { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; }
     20# define ESCAPE { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; }
    2121
    2222int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype) {
    … …  
    4646
    4747  // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_IMAGE
    48   if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE(A);
     48  if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE;
    4949
    5050  // by using "P" format, we are limited to 32bit pointers (2GB heap)
    5151  // XXX how is the "B" format used?
    52   if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE (A);
    53   if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE(A);
     52  if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE;
     53  if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE;
    5454
    5555  // allocates the default data array (ftable->buffer)
    56   if (!gfits_create_table (theader, ftable)) ESCAPE (A);
     56  if (!gfits_create_table (theader, ftable)) ESCAPE;
    5757
    5858  // copy original header to output header (XXX this needs to be finished)
    59   if (!gfits_copy_keywords_compress (header, theader)) ESCAPE (A);
     59  if (!gfits_copy_keywords_compress (header, theader)) ESCAPE;
    6060 
    6161  // this allocates the pointer data array for the full set of tiles (NOT the heap)
    6262  char *tmpbuffer = NULL;
    6363  ALLOCATE_ZERO (tmpbuffer, char, Ntile);
    64   if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE (A);
     64  if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE;
    6565  free (tmpbuffer);
    6666
    6767  // add ZIMAGE from output header
    68   if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE (A);
     68  if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE;
    6969
    7070  // define compression-specific keywords, update header as needed.
    71   if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE(A);
    72 
    73   if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE (A);
     71  if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE;
     72
     73  if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE;
    7474  // gfits_modify (header, "ZTENSION", "%s", 1, exttype);
    7575
    7676  // supply the Z* header values from the source image header
    77   if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE (A);
    78 
    79   if (!gfits_modify (theader, "ZNAXIS",  "%d", 1, header[0].Naxes)) ESCAPE (A);
     77  if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE;
     78
     79  if (!gfits_modify (theader, "ZNAXIS",  "%d", 1, header[0].Naxes)) ESCAPE;
    8080
    8181  char keyword[11];
    8282  for (i = 0; i < header[0].Naxes; i++) {
    8383    snprintf (keyword, 10, "ZNAXIS%d", i + 1);
    84     if (!gfits_modify (theader, keyword, OFF_T_FMT, 1,  header[0].Naxis[i])) ESCAPE (A);
     84    if (!gfits_modify (theader, keyword, OFF_T_FMT, 1,  header[0].Naxis[i])) ESCAPE;
    8585  }
    8686
    8787  for (i = 0; i < header->Naxes; i++) {
    8888    snprintf (keyword, 10, "ZTILE%d", i + 1);
    89     if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE(A);
    90   }
    91 
    92   if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE (A);
    93   if (!gfits_modify (theader, "BZERO",  "%lf", 1, header[0].bzero)) ESCAPE (A);
     89    if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE;
     90  }
     91
     92  if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE;
     93  if (!gfits_modify (theader, "BZERO",  "%lf", 1, header[0].bzero)) ESCAPE;
    9494
    9595  // other keywords to define:
    … …  
    128128
    129129  VarLengthColumn zdef;
    130   if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE(A);
    131 
    132   // allocate the working buffers to the max tile size
     130  if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE;
     131
     132  // pixel sizes and tile sizes:
     133  // * We have Ntile tiles, each of Nx * Ny * Nz ... tile pixels
     134  // -> A pixel in the tile (in the image) has size tile_pixsize = f(BITPIX)
     135  // * We copy the pixels in the tile to a continuous buffer to be compressed.
     136  //   During this copy, we can scale the data (eg, floats -> ints), in which
     137  //   case the raw buffer may have a different pixel size (and type) than the
     138  //   image tile
     139  // -> A pixel in the raw buffer is f(zcmptype,BITPIX):
     140  // -- GZIP_1 : raw_pixsize = tile_pixsize (no scaling is performed)
     141  // -- GZIP_2 : raw_pixsize = tile_pixsize (no scaling is performed)
     142  // -- RICE_1 : raw pixels must be integers. raw_pixsize is user specified? (ZNAMEn,ZVALn : BYTEPIX = 1,2,4,8)
     143  // -- PLIO_1 : raw_pixsize = 2 bytes
     144  // -- HCOMPRESS_1 : raw pixels must be integers. 
     145
     146  // size (in pixels) of the largest tile
    133147  int max_tile_size = gfits_imtile_maxsize (matrix, ztile);
    134148
    135   // size of a pixel in the raw image tile:
    136   int raw_pixsize = abs(header[0].bitpix) / 8;
    137 
    138   // size of a pixel in the raw image tile:
    139   int tile_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
    140 
     149  // size of a pixel in the image tile -- this is probably not needed
     150  int tile_pixsize = abs(header[0].bitpix) / 8;
     151
     152  // size of a pixel in the raw pixel buffer (before compression)
     153  int raw_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
     154  int raw_bitpix  = gfits_uncompressed_data_bitpix (zcmptype, header[0].bitpix, NULL, NULL, 0);
     155  int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
     156
     157  fprintf (stderr, "raw_pixsize: %d, cmp_pixsize: %d, tile_pixsize: %d, raw_bitpix: %d\n",
     158           raw_pixsize, cmp_pixsize, tile_pixsize, raw_bitpix);
     159
     160  // allocate the buffer for compression work
    141161  int Nzdata_alloc = raw_pixsize*max_tile_size + 100;
    142162  ALLOCATE (raw,   char, raw_pixsize*max_tile_size);
    … …  
    144164
    145165  // init the otile[] counters
    146   if (!gfits_imtile_start (matrix, otile)) ESCAPE (A);
     166  if (!gfits_imtile_start (matrix, otile)) ESCAPE;
    147167
    148168  // compress the data : copy into a tile, compress the tile, then add to the output table
    … …  
    150170  for (i = 0; i < Ntile; i++) {
    151171
    152     // size of the current tile
     172    // size of the current tile in pixels
    153173    int Nraw = gfits_tile_size (matrix, otile, ztile);
    154174
    155175    // copy the raw pixels from their native matrix locations to the temporary output buffer
    156     if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE(A);
    157 
     176    // for float -> int scaling by zscale, zzero may be applied
     177    if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     178
     179    if (0 && (i == 0)) {
     180      int k;
     181      fprintf (stderr, "cmp mat: ");
     182      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", matrix->buffer[k]); }
     183      fprintf (stderr, "\n");
     184      fprintf (stderr, "cmp raw: ");
     185      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); }
     186      fprintf (stderr, "\n");
     187    }
     188
     189    // gzip compresses bytes which are in BIG-ENDIAN order
    158190    if (!strcasecmp(zcmptype, "GZIP_1")) {
    159       if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE(A);
    160     }
    161 
    162     // XXX the tile must not be > 2GB
    163    
    164     // optname, optvalue = NULL, Noptions = 0
     191      if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE;
     192    }
     193    if (0 && (i == 0)) {
     194      int k;
     195      fprintf (stderr, "cmp swp: ");
     196      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); }
     197      fprintf (stderr, "\n");
     198    }
     199
     200    // optname, optvalue = NULL, Noptions = 0 : defaults for RICE_1, HCOMPRESS_1
    165201
    166202    int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression
    167     if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, tile_pixsize)) ESCAPE(A);
    168 
    169     if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE(A);
    170 
    171     // XXX need to activate
    172     // if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) ESCAPE(A);
     203    if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, raw_pixsize, ztile[1], ztile[0])) ESCAPE;
     204    if (0 && (i == 0)) {
     205      int k;
     206      fprintf (stderr, "cmp dat: ");
     207      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); }
     208      fprintf (stderr, "\n");
     209    }
     210
     211    if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE;
     212
     213    // all other compression modes require swapping after compression
     214    // (compresssion & decompression operate on native ENDIAN)
     215    if (strcasecmp(zcmptype, "GZIP_1")) {
     216      if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
     217    }
    173218
    174219    // update the otile[] counters, carrying to the next dimension if needed
    … …  
    176221  }
    177222
    178   if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE(A);
     223  if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE;
     224  if (0) {
     225    int k;
     226    fprintf (stderr, "cmp tbl: ");
     227    for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", ftable->buffer[k]); }
     228    fprintf (stderr, "\n");
     229  }
    179230
    180231  free (raw);
    … …  
    209260// raw_pixsize is the bytes / pixel for the input matrix
    210261// place the raw image bytes for the current tile into the tile buffer
    211 int gfits_collect_data (Matrix *matrix, char *raw, int Nraw, int raw_pixsize, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
     262// * otile is the counter for the current output tile
     263// * ztile gives the size of the i-th dimension of the current tile
     264// * raw_bitpix defines the size and type of the raw buffer
     265int gfits_collect_data (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
    212266
    213267  int i, j, start, offset, coord, Nline;
    … …  
    246300  // we need to set up switches for all of the possible combinations:
    247301  // this macro is used at the inner switch to run the actual loop
    248   // note we copy Optr to Iptr
    249 # define SCALE_AND_DIST(TYPE, SIZE) {                                   \
    250     TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];       \
    251     for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) {                    \
    252       if (*Optr == oblank) {                                            \
    253         *Iptr = zblank;                                                 \
     302# define SCALE_AND_DIST_INT(TYPE, SIZE) {                               \
     303    TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];    \
     304    for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) {               \
     305      if (*TILEptr == oblank) {                                         \
     306        *RAWptr = zblank;                                               \
    254307      } else {                                                          \
    255         *Iptr = (*Optr - zzero) / zscale;                               \
     308        *RAWptr = (*TILEptr - zzero) / zscale;                          \
    256309      } } }
    257310
    258311  // we need to set up switches for all of the possible combinations:
    259312  // this macro is used at the inner switch to run the actual loop
    260   // note we copy Optr to Iptr
    261 # define SCALE_AND_DIST_FLOAT(TYPE, SIZE) {                                     \
    262     TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];       \
    263     for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) {                    \
    264       if (!isfinite(*Optr)) {                                           \
    265         *Iptr = zblank;                                                 \
     313# define SCALE_AND_DIST_FLOAT(TYPE, SIZE) {                             \
     314    TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];    \
     315    for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) {               \
     316      if (!isfinite(*TILEptr)) {                                        \
     317        *RAWptr = zblank;                                               \
    266318      } else {                                                          \
    267         *Iptr = (*Optr - zzero) / zscale;                               \
     319        *RAWptr = (*TILEptr - zzero) / zscale;                          \
    268320      } } }
    269321
    270   // this macro sets up the outer switch and calls above macro with all output bitpix options
    271 # define SETUP_INSIZE(TYPE, SIZE) {                     \
    272     TYPE *Iptr = (TYPE *) &raw[i*SIZE*Ztile[0]];        \
     322  // this macro sets up the outer switch and calls above the macro with all RAW buffer bitpix options
     323  //
     324# define SETUP_RAWSIZE(TYPE, SIZE) {                    \
     325    TYPE *RAWptr = (TYPE *) &raw[i*SIZE*Ztile[0]];      \
    273326    switch (matrix->bitpix) {                           \
    274       case 8:                                           \
    275         SCALE_AND_DIST (char, 1);                       \
    276         break;                                          \
    277       case 16:                                          \
    278         SCALE_AND_DIST (short, 2);                      \
    279         break;                                          \
    280       case 32:                                          \
    281         SCALE_AND_DIST (int, 4);                        \
    282         break;                                          \
    283       case -32:                                         \
    284         SCALE_AND_DIST_FLOAT (float, 4);                        \
    285         break;                                          \
    286       case -64:                                         \
    287         SCALE_AND_DIST_FLOAT (double, 8);                       \
    288         break;                                          \
    289       default:                                          \
    290         abort();                                        \
     327      case   8: SCALE_AND_DIST_INT   (char,   1); break;        \
     328      case  16: SCALE_AND_DIST_INT   (short,  2); break;        \
     329      case  32: SCALE_AND_DIST_INT   (int,    4); break;        \
     330      case -32: SCALE_AND_DIST_FLOAT (float,  4); break;        \
     331      case -64: SCALE_AND_DIST_FLOAT (double, 8); break;                                                \
     332      default: abort();                                 \
    291333    } }
    292334
    293335  // loop over lines in the tile
    294336  for (i = 0; i < Nline; i++) {
    295     switch (raw_pixsize) {
    296       case 1:
    297         SETUP_INSIZE (char, 1);
    298         break;
    299       case 2:
    300         SETUP_INSIZE (short, 2);
    301         break;
    302       case 4:
    303         SETUP_INSIZE (int, 4);
    304         break;
    305       case 8:
    306         SETUP_INSIZE (double, 8);
    307         break;
    308       default:
    309         abort();
     337    switch (raw_bitpix) {
     338      case   8: SETUP_RAWSIZE (char,   1); break;
     339      case  16: SETUP_RAWSIZE (short,  2); break;
     340      case  32: SETUP_RAWSIZE (int,    4); break;
     341      case -32: SETUP_RAWSIZE (float,  4); break;
     342      case -64: SETUP_RAWSIZE (double, 8); break;
     343      default: abort();
    310344    }
    311345

  • branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_data.c

    r38340 r38355  
    55/* functions defined in ricecomp.c */
    66int fits_rcomp(int a[], int nx, unsigned char *c, int clen, int nblock);
     7int fits_rcomp_short(short a[], int nx, unsigned char *c, int clen, int nblock);
     8int fits_rcomp_byte(char a[], int nx, unsigned char *c, int clen, int nblock);
    79int fits_rdecomp (unsigned char *c, int clen, unsigned int array[], int nx, int nblock);
    810int fits_rdecomp_short (unsigned char *c, int clen, unsigned short array[], int nx, int nblock);
    … …  
    2931# define ESCAPE(RET) { fprintf (stderr, "gzip error in %s @ %s:%d\n", __func__, __FILE__, __LINE__); return (RET); }
    3032
     33// XXX I'm putting Nx and Ny on the argument list -- only used by hcompress
     34// with a little work, these could go on the option list
    3135int gfits_compress_data (char *zdata, int *Nzdata, char *cmptype,
    3236                         char **optname, char **optvalue, int Nopt,
    33                          char *rawdata, int Nrawpix, int rawpix_size) {
     37                         char *rawdata, int Nrawpix, int rawpix_size,
     38                         int Nx, int Ny) {
    3439
     40  int Nout;
    3541  int status;
    3642
    … …  
    5965  }
    6066
    61 # if (0)
    6267  if (!strcasecmp(cmptype, "RICE_1")) {
    6368    int i, blocksize;
    … …  
    7479    }
    7580
    76     int Npix;
    77     // Npix = *Nout * (out_pixsize / 4.0);
    78     Npix = *Nout;
    79 
    80     Ninsum += Nzdata;
    81     Noutsum += *Nout;
     81    // Ninsum += Nzdata;
     82    // Noutsum += *Nout;
    8283    // fprintf (stderr, "%d comp bytes; %d uncomp 'pixels', totals: %d %d\n", Nzdata, Npix, Ninsum, Noutsum);
    8384
    84     switch (out_pixsize) {
     85    switch (rawpix_size) {
    8586      case 4:
    86         // rice decompression from the CFITSIO source tree : we need to tell it the expected number of pixels
    87         // is also REQUIRES 4byte output, which is fairly stupid.
    88         if (fits_rdecomp ((unsigned char *) zdata, Nzdata, (unsigned int *) outdata, Npix, blocksize)) {
    89           fprintf (stderr, "error in rice decompression\n");
    90           return (FALSE);
    91         }
    92         return (TRUE);
     87        // rice compression from the CFITSIO source tree
     88        Nout = fits_rcomp ((int *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize);
     89        break;
    9390
    9491      case 2:
    95         if (fits_rdecomp_short ((unsigned char *) zdata, Nzdata, (unsigned short *) outdata, Npix, blocksize)) {
    96           fprintf (stderr, "error in rice decompression\n");
    97           return (FALSE);
    98         }
    99         return (TRUE);
     92        Nout = fits_rcomp_short ((short *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize);
     93        break;
    10094
    10195      case 1:
    102         if (fits_rdecomp_byte ((unsigned char *) zdata, Nzdata, (unsigned char *) outdata, Npix, blocksize)) {
    103           fprintf (stderr, "error in rice decompression\n");
    104           return (FALSE);
    105         }
    106         return (TRUE);
     96        Nout = fits_rcomp_byte ((char *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize);
     97        break;
    10798       
    10899      default:
    109         fprintf (stderr, "invalid output pixel size %d\n", out_pixsize);
     100        fprintf (stderr, "invalid output pixel size %d\n", rawpix_size);
     101        *Nzdata = 0;
    110102        return (FALSE);
    111103    }
    112    
     104    if (Nout < 0) {
     105      fprintf (stderr, "error in rice decompression\n");
     106      *Nzdata = 0;
     107      return (FALSE);
     108    }
     109    *Nzdata = Nout;
     110    return TRUE;
    113111  }
    114112 
    115113  if (!strcasecmp(cmptype, "PLIO_1")) {
    116     int Npix;
    117     Npix = pl_l2pi ((short *) zdata, 1, (int *) outdata, *Nout);
    118     if (Npix != *Nout) {
    119       fprintf (stderr, "error in plio decompression\n");
    120       return (FALSE);
    121     }
     114    // note the fortan starting point: zdata is decremented at start
     115    Nout = pl_p2li ((int *) rawdata, 1, (short *) zdata, Nrawpix);
     116    *Nzdata = Nout;
     117    if (!Nout) return (FALSE);
    122118    return (TRUE);
    123119  }
    124120
    125121  if (!strcasecmp(cmptype, "HCOMPRESS_1")) {
    126     int Nx, Ny, scale;
     122    long Nbytes = Nrawpix * rawpix_size;
     123    int scale = 0;
    127124    status = 0;
    128125    // call hdecompress without smoothing
    129     fits_hdecompress ((unsigned char *) zdata, FALSE, (int *) outdata, &Nx, &Ny, &scale, &status);
     126
     127    if (rawpix_size == 4) {
     128      fits_hcompress ((int *) rawdata, Nx, Ny, scale, (char *) zdata, &Nbytes, &status);
     129    }
     130    if (rawpix_size == 8) {
     131      fits_hcompress64 ((long long *) rawdata, Nx, Ny, scale, (char *) zdata, &Nbytes, &status);
     132    }
     133
    130134    if (status) {
    131135      fprintf (stderr, "error in hdecompress\n");
    … …  
    134138    // fprintf (stderr, "decompression yields image %d x %d (scale: %d)\n", Nx, Ny, scale);
    135139   
    136     if (Nx * Ny != *Nout) {
    137       fprintf (stderr, "error in hdecompress: mismatched output size\n");
    138       return (FALSE);
    139     }
     140    *Nzdata = Nbytes / rawpix_size;
    140141    return (TRUE);
    141142  }
    142 # endif
    143143
    144144  fprintf (stderr, "unknown compression %s\n", cmptype);

  • branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_utils.c

    r38340 r38355  
    174174  int i, Nbyte;
    175175
    176   if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "NONE")) {
    177     if (out_bitpix == 8)  return (1);
    178     if (out_bitpix == 16) return (2);
    179     if (out_bitpix == 32) return (4);
     176  // GZIP_1, GZIP_2 use matched input & output pixels
     177  if (!strcasecmp(cmptype, "GZIP_1") ||
     178      !strcasecmp(cmptype, "GZIP_2") ||
     179      !strcasecmp(cmptype, "NONE")) {
     180    if (out_bitpix ==   8) return (1);
     181    if (out_bitpix ==  16) return (2);
     182    if (out_bitpix ==  32) return (4);
    180183    if (out_bitpix == -32) return (4);
    181184    if (out_bitpix == -64) return (8);
    182185    return (1);
    183186  }
    184   if (!strcasecmp(cmptype, "RICE_1")) {
    185 
     187
     188  // RICE_1 must operate on integer pixels
     189  if (!strcasecmp(cmptype, "RICE_1") ||
     190      !strcasecmp(cmptype, "RICE_ONE")) {
    186191    // if BYTEPIX option is specified, use that for Nbyte
    187192    for (i = 0; i < Noptions; i++) {
    … …  
    191196      }
    192197    }
     198    if (out_bitpix ==   8) return (1);
     199    if (out_bitpix ==  16) return (2);
     200    if (out_bitpix ==  32) return (4);
    193201    return (4);
    194202  }
     203
     204  // PLIO_1 always results in 4-byte ints (not 2 byte ints?)
    195205  if (!strcasecmp(cmptype, "PLIO_1")) {
    196206    return (4);
    197207  }
     208
    198209  if (!strcasecmp(cmptype, "HCOMPRESS_1")) {
    199210    if (out_bitpix == 8)  return (4);
    … …  
    204215}
    205216
     217int gfits_compressed_data_pixsize (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) {
     218
     219  int i, Nbyte;
     220
     221  // GZIP_1, GZIP_2 should not be swapped after compression
     222  if (!strcasecmp(cmptype, "GZIP_1") ||
     223      !strcasecmp(cmptype, "GZIP_2") ||
     224      !strcasecmp(cmptype, "NONE")) {
     225    if (out_bitpix ==   8) return (1);
     226    if (out_bitpix ==  16) return (1);
     227    if (out_bitpix ==  32) return (1);
     228    if (out_bitpix == -32) return (1);
     229    if (out_bitpix == -64) return (1);
     230    return (1);
     231  }
     232
     233  // RICE_1 must operate on integer pixels
     234  if (!strcasecmp(cmptype, "RICE_1") ||
     235      !strcasecmp(cmptype, "RICE_ONE")) {
     236    // if BYTEPIX option is specified, use that for Nbyte
     237    for (i = 0; i < Noptions; i++) {
     238      if (!strcmp(optname[i], "BYTEPIX")) {
     239        Nbyte = atoi (optvalue[i]);
     240        return (Nbyte);
     241      }
     242    }
     243    if (out_bitpix ==   8) return (1);
     244    if (out_bitpix ==  16) return (2);
     245    if (out_bitpix ==  32) return (4);
     246    return (4);
     247  }
     248
     249  // PLIO_1 always results in 4-byte ints (not 2 byte ints?)
     250  if (!strcasecmp(cmptype, "PLIO_1")) {
     251    return (4);
     252  }
     253
     254  if (!strcasecmp(cmptype, "HCOMPRESS_1")) {
     255    if (out_bitpix == 8)  return (4);
     256    if (out_bitpix == 16) return (4);
     257    return (8);
     258  }
     259  return (0);
     260}
     261
     262int gfits_uncompressed_data_bitpix (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) {
     263
     264  int i, Nbyte;
     265
     266  // GZIP_1, GZIP_2 use matched input & output pixels
     267  if (!strcasecmp(cmptype, "GZIP_1") ||
     268      !strcasecmp(cmptype, "GZIP_2") ||
     269      !strcasecmp(cmptype, "NONE")) {
     270    return (out_bitpix);
     271  }
     272
     273  // RICE_1 must operate on integer pixels
     274  if (!strcasecmp(cmptype, "RICE_1") ||
     275      !strcasecmp(cmptype, "RICE_ONE")) {
     276    // if BYTEPIX option is specified, use that for Nbyte
     277    for (i = 0; i < Noptions; i++) {
     278      if (!strcmp(optname[i], "BYTEPIX")) {
     279        Nbyte = atoi (optvalue[i]);
     280        return (Nbyte*8);
     281      }
     282    }
     283    if (out_bitpix ==   8) return (out_bitpix);
     284    if (out_bitpix ==  16) return (out_bitpix);
     285    if (out_bitpix ==  32) return (out_bitpix);
     286    return (32);
     287  }
     288
     289  // PLIO_1 always results in 4-byte ints (not 2 byte ints?)
     290  if (!strcasecmp(cmptype, "PLIO_1")) {
     291    return (32);
     292  }
     293
     294  if (!strcasecmp(cmptype, "HCOMPRESS_1")) {
     295    if (out_bitpix == 8)  return (32);
     296    if (out_bitpix == 16) return (32);
     297    return (64);
     298  }
     299  return (0);
     300}
     301
    206302int gfits_vartable_heap_pixsize (char format) {
    207303

  • branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_M.c

    r38340 r38355  
    4646  float zscale, zzero;
    4747
    48   int zdata_pixsize, odata_pixsize;
    49 
    5048  int *ztile = NULL;
    5149  int *otile = NULL;
    … …  
    242240
    243241  // size of an element in the vartable heap section
    244   zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format);
     242  // zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); XXX
    245243
    246244  // size of a pixel in the final image (not needed)
    … …  
    248246
    249247  // size of a pixel in the output from the decompression routine
    250   odata_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions);
    251 
    252   int Nout_alloc = odata_pixsize*max_tile_size;
     248  int cmp_pixsize = gfits_compressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions);
     249  int raw_bitpix  = gfits_uncompressed_data_bitpix (cmptype, header[0].bitpix, optname, optvalue, Noptions);
     250  int raw_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions);
     251
     252  int Nout_alloc = raw_pixsize*max_tile_size;
    253253  ALLOCATE (out, char, Nout_alloc);
    254254
    255255  off_t row;
     256
     257  if (0) {
     258    int k;
     259    fprintf (stderr, "unc tbl: ");
     260    for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", ftable->buffer[k]); }
     261    fprintf (stderr, "\n");
     262  }
    256263
    257264  // uncompress the data
    … …  
    266273    if (!zdata) ESCAPE;
    267274
    268     // XXX not certain this is the correct place to do the swap (or if zdata_pixsize is
    269     // the correct size to guide the swap)
    270 
    271275    if (strcasecmp(cmptype, "GZIP_1")) {
    272       if (!gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize)) ESCAPE;
     276      if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
     277    }
     278
     279    if (0 && (row == 0)) {
     280      int k;
     281      fprintf (stderr, "unc dat: ");
     282      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); }
     283      fprintf (stderr, "\n");
    273284    }
    274285
    275286    // gfits_uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated
    276287    // XXX the tile must not be > 2GB
    277     if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, odata_pixsize)) ESCAPE;
    278 
     288    if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, raw_pixsize)) ESCAPE;
     289
     290    if (0 && (row == 0)) {
     291      int k;
     292      fprintf (stderr, "unc swp: ");
     293      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", out[k]); }
     294      fprintf (stderr, "\n");
     295    }
     296   
    279297    if (!strcasecmp(cmptype, "GZIP_1")) {
    280       if (!gfits_byteswap_zdata (out, Nout, odata_pixsize)) ESCAPE;
     298      if (!gfits_byteswap_zdata (out, Nout, raw_pixsize)) ESCAPE;
     299    }
     300   
     301    if (0 && (row == 0)) {
     302      int k;
     303      fprintf (stderr, "unc raw: ");
     304      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", out[k]); }
     305      fprintf (stderr, "\n");
    281306    }
    282307
    283308    // copy the uncompressed pixels into their correct locations           
    284     if (!gfits_distribute_data (matrix, odata_pixsize, out, Nout, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     309    if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     310
     311    if (0 && (row == 0)) {
     312      int k;
     313      fprintf (stderr, "unc mat: ");
     314      for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", matrix->buffer[k]); }
     315      fprintf (stderr, "\n");
     316    }
    285317
    286318    // update the tile counters, carrying to the next dimension if needed
    … …  
    311343
    312344// bitpix is the input data size/type
    313 int gfits_distribute_data (Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
     345int gfits_distribute_data (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
    314346
    315347  int i, j, start, offset, coord, Nline;
    … …  
    334366
    335367  // double check reported size (needs the Isize)
    336   assert (Ndata == Ztile[0]*Nline);
     368  assert (Nraw == Ztile[0]*Nline);
    337369
    338370  // set the starting point of the tile:
    … …  
    350382  // we need to set up switches for all of the possible combinations:
    351383  // this macro is used at the inner switch to run the actual loop
    352 # define SCALE_AND_DIST(TYPE, SIZE) { \
    353     TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \
    354     for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \
    355       if (*Iptr == zblank) { \
    356           *Optr = oblank;    \
    357       } else { \
    358           *Optr = *Iptr * zscale + zzero;       \
     384# define SCALE_AND_DIST_INT(TYPE, SIZE) {                               \
     385    TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];    \
     386    for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) {               \
     387      if (*TILEptr == zblank) {                                         \
     388        *TILEptr = oblank;                                              \
     389      } else {                                                          \
     390        *TILEptr = *RAWptr * zscale + zzero;                            \
    359391      } } }
    360392
     393  // we need to set up switches for all of the possible combinations:
     394  // this macro is used at the inner switch to run the actual loop
     395# define SCALE_AND_DIST_FLOAT(TYPE, SIZE) {                             \
     396    TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];    \
     397    for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) {               \
     398      if (!isfinite(*TILEptr)) {                                        \
     399        *TILEptr = oblank;                                              \
     400      } else {                                                          \
     401        *TILEptr = *RAWptr * zscale + zzero;                            \
     402      } } }
     403 
    361404  // this macro sets up the outer switch and calls above macro with all output bitpix options
    362 # define SETUP_INSIZE(TYPE, SIZE) { \
    363     TYPE *Iptr = (TYPE *) &data[i*SIZE*Ztile[0]]; \
    364     switch (matrix->bitpix) { \
    365       case 8: \
    366         SCALE_AND_DIST (char, 1); \
    367         break; \
    368       case 16: \
    369         SCALE_AND_DIST (short, 2); \
    370         break; \
    371       case 32: \
    372         SCALE_AND_DIST (int, 4); \
    373         break; \
    374       case -32: \
    375         SCALE_AND_DIST (float, 4); \
    376         break; \
    377       case -64: \
    378         SCALE_AND_DIST (double, 8); \
    379         break; \
    380       default: \
    381         abort(); \
     405# define SETUP_RAWSIZE(TYPE, SIZE) {               \
     406    TYPE *RAWptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \
     407    switch (matrix->bitpix) {                      \
     408      case   8: SCALE_AND_DIST_INT   (char,   1); break; \
     409      case  16: SCALE_AND_DIST_INT   (short,  2); break; \
     410      case  32: SCALE_AND_DIST_INT   (int,    4); break; \
     411      case -32: SCALE_AND_DIST_FLOAT (float,  4); break; \
     412      case -64: SCALE_AND_DIST_FLOAT (double, 8); break; \
     413      default:  abort(); \
    382414    } }
    383415
    384416  // loop over lines in the tile
    385417  for (i = 0; i < Nline; i++) {
    386     switch (bitpix) {
    387       case 1:
    388         SETUP_INSIZE (char, 1);
    389         break;
    390       case 2:
    391         SETUP_INSIZE (short, 2);
    392         break;
    393       case 4:
    394         SETUP_INSIZE (int, 4);
    395         break;
    396       case 8:
    397         SETUP_INSIZE (double, 8);
    398         break;
    399       default:
    400         abort();
     418    switch (raw_bitpix) {
     419      case   8: SETUP_RAWSIZE (char,   1); break;
     420      case  16: SETUP_RAWSIZE (short,  2); break;
     421      case  32: SETUP_RAWSIZE (int,    4); break;
     422      case -32: SETUP_RAWSIZE (float,  4); break;
     423      case -64: SETUP_RAWSIZE (double, 8); break;
     424      default:  abort();
    401425    }
    402426

  • branches/eam_branches/ohana.20150429/src/libfits/table/F_compress_T.c

    r38348 r38355  
    161161
    162162      int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression
    163       if (!gfits_compress_data (zdata, &Nzdata, fields[j].zctype, NULL, NULL, 0, raw, Nraw, fields[j].pixsize)) ESCAPE(A);
     163      if (!gfits_compress_data (zdata, &Nzdata, fields[j].zctype, NULL, NULL, 0, raw, Nraw, fields[j].pixsize, 0, 0)) ESCAPE(A);
    164164     
    165165      if (!gfits_varlength_column_add_data (tgttable, zdata, Nzdata, i, &fields[j].zdef)) ESCAPE(A);
Note: See TracChangeset for help on using the changeset viewer.