summaryrefslogtreecommitdiffstats
path: root/tksao/fitsy++/gzip.C
diff options
context:
space:
mode:
Diffstat (limited to 'tksao/fitsy++/gzip.C')
-rw-r--r--tksao/fitsy++/gzip.C339
1 files changed, 339 insertions, 0 deletions
diff --git a/tksao/fitsy++/gzip.C b/tksao/fitsy++/gzip.C
new file mode 100644
index 0000000..6f83af6
--- /dev/null
+++ b/tksao/fitsy++/gzip.C
@@ -0,0 +1,339 @@
+// Copyright (C) 1999-2016
+// 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(int);
+ 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::NODITHER:
+ val = FitsCompressm<T>::getValue((float*)obuf+ll,zs,zz,blank);
+ break;
+ 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::NODITHER:
+ val = FitsCompressm<T>::getValue((double*)obuf+ll,zs,zz,blank);
+ break;
+ 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>;
+