IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38362


Ignore:
Timestamp:
Jun 3, 2015, 8:20:19 PM (11 years ago)
Author:
eugene
Message:

adding GZIP_2 and getting RICE_1 to work (mostly)

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

Legend:

Unmodified
Added
Removed

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

    r38356 r38362  
    33# include <zlib.h>
    44# define VERBOSE_DUMP 0
     5
     6int gfits_collect_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero);
    57
    68// the user needs to specify the various compression options:
    … …  
    7375  // define compression-specific keywords, update header as needed.
    7476  if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE;
     77
     78  // define compression-specific keywords, update header as needed.
     79  if (!strcasecmp(zcmptype, "GZIP_1") || !strcasecmp(zcmptype, "GZIP_2")) {
     80    if ((header[0].bitpix == -32) || (header[0].bitpix == -64)) {
     81      if (!gfits_modify (theader, "ZQUANTIZ", "%s", 1, "NONE")) ESCAPE;
     82    }
     83  }
    7584
    7685  if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE;
    … …  
    155164  int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
    156165
     166  // define compression-specific keywords, update header as needed.
     167  if (!strcasecmp(zcmptype, "RICE_1")) {
     168    if (!gfits_modify (theader, "ZNAME1", "%s", 1, "BLOCKSIZE")) ESCAPE;
     169    if (!gfits_modify (theader, "ZVAL1", "%d", 1, 32)) ESCAPE;
     170    if (!gfits_modify (theader, "ZNAME2", "%s", 1, "BYTEPIX")) ESCAPE;
     171    if (!gfits_modify (theader, "ZVAL2", "%d", 1, raw_pixsize)) ESCAPE;
     172  }
     173
    157174  // size of a pixel in the image tile -- this is probably not needed
    158175  // int tile_pixsize = abs(header[0].bitpix) / 8;
    … …  
    176193    // copy the raw pixels from their native matrix locations to the temporary output buffer
    177194    // for float -> int scaling by zscale, zzero may be applied
    178     if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     195    if (!strcasecmp (zcmptype, "GZIP_2")) {
     196      if (!gfits_collect_gzp2 (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     197    } else {
     198      if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     199    }
    179200
    180201    if (VERBOSE_DUMP && (i == 0)) {
    … …  
    214235    }
    215236
    216     // all other compression modes require swapping after compression
    217     // (compresssion & decompression operate on native ENDIAN)
    218     if (strcasecmp(zcmptype, "GZIP_1")) {
     237    // GZIP_1 is swapped before compression; GZIP_2 implies swapping.  All other
     238    // compression modes require swapping after compression (compresssion & decompression
     239    // operate on native ENDIAN)
     240    if (strcasecmp(zcmptype, "GZIP_1") && strcasecmp(zcmptype, "GZIP_2")) {
    219241      // Nzdata is number of bytes
    220242      if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
    … …  
    383405  return (TRUE);
    384406}
     407
     408// GZIP_2 organizes data in bytes by significance: (A1 A2 A3 A4), (B1 B2 B3 B4) -> (A1 B1), (A2, B2), (A3, B3), (A4, B4)
     409// NOTE that this organization depends on the machine endian state.
     410// raw_pixsize is the bytes / pixel for the input matrix
     411// place the raw image bytes for the current tile into the tile buffer
     412// * otile is the counter for the current output tile
     413// * ztile gives the size of the i-th dimension of the current tile
     414// * raw_bitpix defines the size and type of the raw buffer
     415int gfits_collect_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
     416
     417  int i, j, k;
     418  int *counter = NULL;
     419  int *Ztile = NULL;
     420
     421  ALLOCATE (counter, int, matrix->Naxes); // counter for current row in tile to copy
     422  ALLOCATE (Ztile, int, matrix->Naxes); // true sizes of this tile (in pixels)
     423 
     424  for (i = 0; i < matrix->Naxes; i++) {
     425    counter[i] = 0;
     426    Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]);
     427  }
     428
     429  // number of lines in the tile (in pixels)
     430  int Nline = 1;
     431  int Npix = matrix->Naxis[0];
     432  for (i = 1; i < matrix->Naxes; i++) {
     433    Nline *= Ztile[i];
     434    Npix  *= matrix->Naxis[i];
     435  }
     436
     437  // confirm the tile size:
     438  assert (Nraw == Ztile[0]*Nline);
     439
     440  // set the starting point of the tile in the matrix buffer
     441  // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...;
     442  // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...));
     443  int start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
     444  for (i = matrix->Naxes - 2; i >= 0; i--) {
     445    int coord = otile[i]*ztile[i];
     446    start = start*matrix->Naxis[i] + coord;
     447  }
     448 
     449  int size = 0;
     450  switch (raw_bitpix) {
     451    case   8: size = 1; break;
     452    case  16: size = 2; break;
     453    case  32: size = 4; break;
     454    case -32: size = 4; break;
     455    case -64: size = 8; break;
     456    default: abort();
     457  }
     458
     459  // pixel offset in output array relative to tile start
     460  int offset = 0;
     461 
     462  static int pass = 0;
     463  for (k = 0; k < size; k++) {
     464    for (i = 0; i < Nline; i++) {
     465# ifdef BYTE_SWAP     
     466      char *srcptr = &matrix->buffer[size*(offset + start) + (size - k - 1)];
     467# else
     468      char *srcptr = &matrix->buffer[size*(offset + start) + k];
     469# endif
     470      char *rawptr = &raw[i*size*Ztile[0] + k*Nraw];   
     471      for (j = 0; j < Ztile[0]; j++, srcptr += size, rawptr ++) {               
     472        if (FALSE && (i == 0) && (j < 4) && (pass == 0)) {
     473          fprintf (stderr, "0x%02hhx -> 0x%02hhx\n", (int)(srcptr - matrix->buffer), (int)(rawptr - raw));
     474        }
     475        *rawptr = *srcptr;
     476        myAssert (srcptr - matrix->buffer < Npix*size, "memory overrun");
     477        myAssert (rawptr - raw < Nraw*size, "memory overrun");
     478      }
     479
     480      // update the counters, carrying to the next dimension if needed
     481      for (j = 1; j < matrix->Naxes; j++) {
     482        counter[j] ++;
     483        if (counter[j] == Ztile[j]) {
     484          counter[j] = 0;
     485        } else {
     486          break;
     487        }
     488      }
     489      if (j == matrix->Naxes) assert (i == Nline - 1); // we should be done here...
     490
     491      // Naxes = 3
     492      // offset = counter[1]*matrix->Naxis[0] + counter[2]*matrix->Naxis[0]*matrix->Naxis[1] +
     493      // offset = matrix->Naxis[0]*(counter[1] + matrix->Naxis[1]*(counter[2] + matrix->Naxis[2]*...))
     494
     495      // determine the offset of the next line relative to the start position
     496      offset = counter[matrix->Naxes - 1];
     497      for (j = matrix->Naxes - 2; j >= 0; j--) {
     498        offset = offset*matrix->Naxis[j] + counter[j];
     499      }     
     500    }
     501  }
     502  pass ++;
     503
     504  free (counter);
     505  free (Ztile);
     506
     507  return (TRUE);
     508}

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

    r38356 r38362  
    5050  }
    5151
    52   if (!strcasecmp(cmptype, "GZIP_1")) {
     52  // GZIP_1 uses inflate / deflate, GZIP_2 does as well, but operates on a buffer with
     53  // a different byte organization (gfits_distribute_gzp2)
     54  if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "GZIP_2")) {
    5355    unsigned long Nbytes = Nrawpix * rawpix_size;
    5456

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

    r38356 r38362  
    33# include <zlib.h>
    44# define VERBOSE_DUMP 0
     5
     6int gfits_distribute_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero);
    57
    68# define ESCAPE { \
    … …  
    321323
    322324    // copy the uncompressed pixels into their correct locations           
    323     if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     325    if (!strcasecmp(cmptype, "GZIP_2")) {
     326      if (!gfits_distribute_gzp2 (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     327    } else {
     328      if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
     329    }
    324330
    325331    if (VERBOSE_DUMP && (row == 0)) {
    … …  
    359365int gfits_distribute_data (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
    360366
    361   int i, j, start, offset, coord, Nline;
     367  int i, j;
    362368  int *counter = NULL;
    363369  int *Ztile = NULL;
    … …  
    374380
    375381  // number of lines in the tile (in pixels)
    376   Nline = 1;
     382  int Nline = 1;
    377383  for (i = 1; i < matrix->Naxes; i++) {
    378384    Nline *= Ztile[i];
    … …  
    385391  // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...;
    386392  // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...));
    387   start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
     393  int start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
    388394  for (i = matrix->Naxes - 2; i >= 0; i--) {
    389     coord = otile[i]*ztile[i];
     395    int coord = otile[i]*ztile[i];
    390396    start = start*matrix->Naxis[i] + coord;
    391397  }
    392398 
    393399  // pixel offset in output array relative to tile start
    394   offset = 0;
     400  int offset = 0;
    395401
    396402  // we need to set up switches for all of the possible combinations:
    … …  
    466472}
    467473
     474// bitpix is the input data size/type
     475int gfits_distribute_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
     476
     477  int i, j, k;
     478  int *counter = NULL;
     479  int *Ztile = NULL;
     480
     481  ALLOCATE (counter, int, matrix->Naxes);
     482  ALLOCATE (Ztile, int, matrix->Naxes);
     483
     484  // counter for current row in tile to copy
     485  // true sizes of this tile (in pixels)
     486  for (i = 0; i < matrix->Naxes; i++) {
     487    counter[i] = 0;
     488    Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]);
     489  }
     490
     491  // number of lines in the tile (in pixels)
     492  int Nline = 1;
     493  int Npix = matrix->Naxis[0];
     494  for (i = 1; i < matrix->Naxes; i++) {
     495    Nline *= Ztile[i];
     496    Npix  *= matrix->Naxis[i];
     497  }
     498
     499  // double check reported size (needs the Isize)
     500  assert (Nraw == Ztile[0]*Nline);
     501
     502  // set the starting point of the tile:
     503  // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...;
     504  // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...));
     505  int start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
     506  for (i = matrix->Naxes - 2; i >= 0; i--) {
     507    int coord = otile[i]*ztile[i];
     508    start = start*matrix->Naxis[i] + coord;
     509  }
     510 
     511  int size = 0;
     512  switch (raw_bitpix) {
     513    case   8: size = 1; break;
     514    case  16: size = 2; break;
     515    case  32: size = 4; break;
     516    case -32: size = 4; break;
     517    case -64: size = 8; break;
     518    default: abort();
     519  }
     520
     521  // pixel offset in output array relative to tile start
     522  int offset = 0;
     523
     524  static int pass = 0;
     525  for (k = 0; k < size; k++) {
     526    for (i = 0; i < Nline; i++) {
     527# ifdef BYTE_SWAP     
     528      char *srcptr = &matrix->buffer[size*(offset + start) + (size - k - 1)];
     529# else
     530      char *srcptr = &matrix->buffer[size*(offset + start) + k];
     531# endif
     532      char *rawptr = &raw[i*size*Ztile[0] + k*Nraw];   
     533      for (j = 0; j < Ztile[0]; j++, srcptr += size, rawptr ++) {               
     534        if (FALSE && (i == 0) && (j < 4) && (pass == 0)) {
     535          fprintf (stderr, "0x%02hhx -> 0x%02hhx\n", (int)(srcptr - matrix->buffer), (int)(rawptr - raw));
     536        }
     537        *srcptr = *rawptr;
     538        myAssert (srcptr - matrix->buffer < Npix*size, "memory overrun");
     539        myAssert (rawptr - raw < Nraw*size, "memory overrun");
     540      }
     541
     542      // update the counters, carrying to the next dimension if needed
     543      for (j = 1; j < matrix->Naxes; j++) {
     544        counter[j] ++;
     545        if (counter[j] == Ztile[j]) {
     546          counter[j] = 0;
     547        } else {
     548          break;
     549        }
     550      }
     551      if (j == matrix->Naxes) assert (i == Nline - 1); // we should be done here...
     552
     553      // Naxes = 3
     554      // offset = counter[1]*matrix->Naxis[0] + counter[2]*matrix->Naxis[0]*matrix->Naxis[1] +
     555      // offset = matrix->Naxis[0]*(counter[1] + matrix->Naxis[1]*(counter[2] + matrix->Naxis[2]*...))
     556
     557      // determine the offset of the next line relative to the start position
     558      offset = counter[matrix->Naxes - 1];
     559      for (j = matrix->Naxes - 2; j >= 0; j--) {
     560        offset = offset*matrix->Naxis[j] + counter[j];
     561      }     
     562    }
     563  }
     564  pass ++;
     565
     566  free (counter);
     567  free (Ztile);
     568  return (TRUE);
     569}
     570

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

    r38356 r38362  
    3939  }
    4040
    41   if (!strcasecmp(cmptype, "GZIP_1")) {
     41  if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "GZIP_2")) {
    4242    // unsigned long tNout = *Nout * out_pixsize;
    4343    // input value of *Nout is number of BYTES

  • branches/eam_branches/ohana.20150429/src/libfits/test/compress.sh

    r38356 r38362  
    22macro test_image_all
    33
    4   foreach mode NONE GZIP_1 RICE_1
     4  # foreach mode NONE GZIP_1 RICE_1
     5  # foreach mode GZIP_1
     6  foreach mode GZIP_2
     7  # foreach mode RICE_1
    58    foreach bitpix 8 16 32 -32 -64
     9    # foreach bitpix 32
    610      echo ===== bitpix = $bitpix =====
    711      test_image_single $bitpix $mode
    … …  
    3034 if ($cmpmode == NONE) return
    3135
    32  if ($bitpix == -64)
    33    echo "*** bitpix $bitpix fails fpack / funpack ***"
    34  else
    35    exec fpack -g -S test.raw.fits > test.fpk.fits
    36    echo "load fpack version"
    37    rd c test.fpk.fits -x 0
    38    set d = x - c
    39    stat d
    40  end
    41 
    42  if ($cmpmode == RICE_1)
     36 if ($cmpmode == GZIP_1)
     37   exec fpack -q 0 -g -S test.raw.fits > test.fpk.fits
     38 end
     39 if ($cmpmode == GZIP_2)
     40   exec fpack -q 0 -g2 -S test.raw.fits > test.fpk.fits
     41 end
     42 if ($cmpmode == RICE_1)
     43   exec fpack -r -S test.raw.fits > test.fpk.fits
     44 end
     45 echo "load fpack version"
     46 rd c test.fpk.fits -x 0
     47 set d = x - c
     48 stat d
     49
     50 if ($cmpmode == RICE_1_0)
    4351   echo "*** funpack fails on RICE_1 ****"
    4452   return
    4553 end
    4654
    47  if ($bitpix == -64)
    48    echo "*** bitpix $bitpix fails fpack / funpack ***"
    49  else
    50    echo "run funpack on our version"
    51    exec funpack -S test.cmp.fits > test.fun.fits
    52    rd e test.fun.fits
    53    set d = x - e
    54    stat d
    55  end
     55 echo "run funpack on our version"
     56 exec funpack -S test.cmp.fits > test.fun.fits
     57 rd e test.fun.fits
     58 set d = x - e
     59 stat d
    5660end
    5761

  • branches/eam_branches/ohana.20150429/src/libfits/test/imagecomp.c

    r38356 r38362  
    66
    77// char *cmptype[] = {"NONE", "GZIP_1", "GZIP_2", "PLIO_1", "RICE_1", "RICE_ONE", "HCOMPRESS_1", NULL};
    8 char *cmptype[] = {"NONE", "GZIP_1", "RICE_1", "RICE_ONE", NULL};
    9 int bitpix[] = {8, 16, 32, -32, -64, 0};
     8// char *cmptype[] = {"NONE", "GZIP_1", "RICE_1", "RICE_ONE", NULL};
     9// int bitpix[] = {8, 16, 32, -32, -64, 0};
    1010
    1111// char *cmptype[] = {"NONE", "GZIP_1", NULL};
    … …  
    1313// int bitpix[] = {8, 16, 32, 0};
    1414
    15 // char *cmptype[] = {"RICE_1", NULL};
    16 // int bitpix[] = {-32, -64, 0};
     15char *cmptype[] = {"GZIP_2", NULL};
     16int bitpix[] = {-32, 0};
    1717
    1818int main (int argc, char **argv) {
    … …  
    4343}
    4444
    45 # define NX 100
     45# define NX 10
    4646# define NY 10
    4747
Note: See TracChangeset for help on using the changeset viewer.