diff options
Diffstat (limited to 'fitsy++/gzip.C')
-rw-r--r-- | fitsy++/gzip.C | 341 |
1 files changed, 341 insertions, 0 deletions
diff --git a/fitsy++/gzip.C b/fitsy++/gzip.C new file mode 100644 index 0000000..47ccb54 --- /dev/null +++ b/fitsy++/gzip.C @@ -0,0 +1,341 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#include <iostream> +#include <sstream> +#include <iomanip> +using namespace std; + +#include "gzip.h" +#include "zlib.h" +#include "util.h" + +template<class T> FitsGzipm<T>::FitsGzipm(FitsFile* fits) + : FitsCompressm<T>(fits) +{ + FitsCompressm<T>::uncompress(fits); +} + +template <class T> int FitsGzipm<T>::compressed(T* dest, char* sptr, + char* heap, + int kkstart, int kkstop, + int jjstart, int jjstop, + int iistart, int iistop) +{ + double zs = FitsCompressm<T>::bscale_; + if (FitsCompressm<T>::zscale_) + zs = FitsCompressm<T>::zscale_->value(sptr,0); + + double zz = FitsCompressm<T>::bzero_; + if (FitsCompressm<T>::zzero_) + zz = FitsCompressm<T>::zzero_->value(sptr,0); + + int blank = FitsCompressm<T>::blank_; + if (FitsCompressm<T>::zblank_) + blank = (int)FitsCompressm<T>::zblank_->value(sptr,0); + + int icnt=0; + unsigned char* ibuf = (unsigned char*)((FitsBinColumnArray*)FitsCompressm<T>::compress_)->get(heap, sptr, &icnt); + + // ibuf can be NULL + if (!ibuf || !icnt) + return 0; + + int ocnt = FitsCompressm<T>::tilesize_; + char* obuf = new char[ocnt*sizeof(long long)]; + if (!obuf) { + internalError("Fitsy++ gzip unable to alloc."); + return 0; + } + + z_stream zstrm; + zstrm.next_in = NULL; + zstrm.avail_in = 0; + zstrm.zalloc = NULL; + zstrm.zfree = NULL; + zstrm.opaque = NULL; + + // look for both zlib and gzip headers + if (inflateInit2(&zstrm, MAX_WBITS+32) != Z_OK) { + internalError("Fitsy++ gzip inflateInit error"); + if (obuf) + delete [] obuf; + return 0; + } + + zstrm.avail_in = icnt; + zstrm.next_in = ibuf; + zstrm.avail_out = ocnt*sizeof(long long); + zstrm.next_out = (Bytef*)obuf; + + if (DebugCompress) + cerr << " inflate START: avail_in " << zstrm.avail_in + << " avail_out " << zstrm.avail_out + << " total_in " << zstrm.total_in + << " total_out " << zstrm.total_out << endl; + + int result = ::inflate(&zstrm, Z_FINISH); + + switch (result) { + case Z_OK: + if (DebugCompress) + cerr << " inflate OK: avail_in " << zstrm.avail_in + << " avail_out " << zstrm.avail_out + << " total_in " << zstrm.total_in + << " total_out " << zstrm.total_out << endl; + break; + case Z_STREAM_END: + if (DebugCompress) + cerr << " inflate STREAM_END: avail_in " << zstrm.avail_in + << " avail_out " << zstrm.avail_out + << " total_in " << zstrm.total_in + << " total_out " << zstrm.total_out << endl; + break; + case Z_BUF_ERROR: + if (DebugCompress) + cerr << " inflate BUF_ERROR: avail_in " << zstrm.avail_in + << " avail_out " << zstrm.avail_out << endl; + if (obuf) + delete [] obuf; + return 0; + default: + internalError("Fitsy++ gzip inflate error"); + if (obuf) + delete [] obuf; + return 0; + } + + int bytepix = zstrm.total_out/FitsCompressm<T>::tilesize_; + + inflateEnd(&zstrm); + + // GZIP_2- unshuffle if needed + if (!strncmp(FitsCompressm<T>::type_,"GZIP_2",6)) { + switch (bytepix) { + case 1: + break; + case 2: + { + int ll = ocnt*sizeof(short); + char* nbuf = new char[ll]; + if (!nbuf) { + internalError("Fitsy++ gzip unable to alloc."); + if (obuf) + delete [] obuf; + return 0; + } + char* optr = obuf+ll-1; + char* nptr = nbuf+ll-1; + + for (int ii=0; ii<ocnt; ii++) { + *nptr = *optr; + nptr--; + *nptr = *(optr-ocnt); + nptr--; + + optr--; + } + memcpy(obuf,nbuf,ll); + delete [] nbuf; + } + break; + case 4: + { + int ll = ocnt*sizeof(int); + char* nbuf = new char[ll]; + if (!nbuf) { + internalError("Fitsy++ gzip unable to alloc."); + if (obuf) + delete [] obuf; + return 0; + } + char* optr = obuf+ll-1; + char* nptr = nbuf+ll-1; + + for (int ii=0; ii<ocnt; ii++) { + *nptr = *optr; + nptr--; + *nptr = *(optr-ocnt); + nptr--; + *nptr = *(optr-(2*ocnt)); + nptr--; + *nptr = *(optr-(3*ocnt)); + nptr--; + + optr--; + } + memcpy(obuf,nbuf,ll); + delete [] nbuf; + } + break; + case 8: + { + int ll = ocnt*sizeof(long long); + char* nbuf = new char[ll]; + if (!nbuf) { + internalError("Fitsy++ gzip unable to alloc."); + if (obuf) + delete [] obuf; + return 0; + } + char* optr = obuf+ll-1; + char* nptr = nbuf+ll-1; + + for (int ii=0; ii<ocnt; ii++) { + *nptr = *optr; + nptr--; + *nptr = *(optr-ocnt); + nptr--; + *nptr = *(optr-(2*ocnt)); + nptr--; + *nptr = *(optr-(3*ocnt)); + nptr--; + *nptr = *(optr-(4*ocnt)); + nptr--; + *nptr = *(optr-(5*ocnt)); + nptr--; + *nptr = *(optr-(6*ocnt)); + nptr--; + *nptr = *(optr-(7*ocnt)); + nptr--; + + optr--; + } + memcpy(obuf,nbuf,ll); + delete [] nbuf; + } + break; + } + } + + int ll=0; + switch (bytepix) { + case 1: + for (int kk=kkstart; kk<kkstop; kk++) + for (int jj=jjstart; jj<jjstop; jj++) + for (int ii=iistart; ii<iistop; ii++,ll++) { + size_t id = kk*FitsCompressm<T>::width_*FitsCompressm<T>::height_ + jj*FitsCompressm<T>::width_ + ii; + // very carefull about type conversions + T val = FitsCompressm<T>::getValue(obuf+ll,zs,zz,blank); + dest[id] = val; + } + break; + case 2: + for (int kk=kkstart; kk<kkstop; kk++) + for (int jj=jjstart; jj<jjstop; jj++) + for (int ii=iistart; ii<iistop; ii++,ll++) { + // swap if needed + if (FitsCompressm<T>::byteswap_) { + const char* p = (const char*)((short*)obuf+ll); + union { + char c[2]; + short s; + } u; + + u.c[1] = *p++; + u.c[0] = *p; + + *((short*)obuf+ll) = u.s; + } + // very carefull about type conversions + size_t id = kk*FitsCompressm<T>::width_*FitsCompressm<T>::height_ + jj*FitsCompressm<T>::width_ + ii; + T val = FitsCompressm<T>::getValue((short*)obuf+ll,zs,zz,blank); + dest[id] = val; + } + break; + case 4: + for (int kk=kkstart; kk<kkstop; kk++) + for (int jj=jjstart; jj<jjstop; jj++) + for (int ii=iistart; ii<iistop; ii++,ll++) { + // swap if needed + if (FitsCompressm<T>::byteswap_) { + const char* p = (const char*)((int*)obuf+ll); + union { + char c[4]; + int i; + } u; + + u.c[3] = *p++; + u.c[2] = *p++; + u.c[1] = *p++; + u.c[0] = *p; + + *((int*)obuf+ll) = u.i; + } + // very carefull about type conversions + size_t id = kk*FitsCompressm<T>::width_*FitsCompressm<T>::height_ + jj*FitsCompressm<T>::width_ + ii; + T val =0; + switch (FitsCompressm<T>::quantize_) { + case FitsCompress::NONE: + val = FitsCompressm<T>::getValue((float*)obuf+ll,zs,zz,blank); + break; + case FitsCompress::NODITHER: + case FitsCompress::SUBDITHER1: + case FitsCompress::SUBDITHER2: + val = FitsCompressm<T>::getValue((int*)obuf+ll,zs,zz,blank); + break; + } + dest[id] = val; + } + break; + case 8: + for (int kk=kkstart; kk<kkstop; kk++) + for (int jj=jjstart; jj<jjstop; jj++) + for (int ii=iistart; ii<iistop; ii++,ll++) { + // swap if needed + if (FitsCompressm<T>::byteswap_) { + const char* p = (const char*)((long long*)obuf+ll); + union { + char c[8]; + long long i; + } u; + + u.c[7] = *p++; + u.c[6] = *p++; + u.c[5] = *p++; + u.c[4] = *p++; + u.c[3] = *p++; + u.c[2] = *p++; + u.c[1] = *p++; + u.c[0] = *p; + + *((long long*)obuf+ll) = u.i; + } + // very carefull about type conversions + size_t id = kk*FitsCompressm<T>::width_*FitsCompressm<T>::height_ + jj*FitsCompressm<T>::width_ + ii; + T val =0; + switch (FitsCompressm<T>::quantize_) { + case FitsCompress::NONE: + val = FitsCompressm<T>::getValue((double*)obuf+ll,zs,zz,blank); + break; + case FitsCompress::NODITHER: + case FitsCompress::SUBDITHER1: + case FitsCompress::SUBDITHER2: + val = FitsCompressm<T>::getValue((long long*)obuf+ll,zs,zz,blank); + break; + } + dest[id] = val; + } + break; + + default: + internalError("Fitsy++ gzip illegal bytepix"); + if (obuf) + delete [] obuf; + return 0; + } + + if (obuf) + delete [] obuf; + return 1; +} + +template class FitsGzipm<unsigned char>; +template class FitsGzipm<short>; +template class FitsGzipm<unsigned short>; +template class FitsGzipm<int>; +template class FitsGzipm<long long>; +template class FitsGzipm<float>; +template class FitsGzipm<double>; + |