Changeset 17687 for fact/tools/pyscripts/pyfact/fits.h

Timestamp:

04/24/14 15:58:56 (10 years ago)

Author:

dneise

Message:

differs from Mars/mcore only in special python accessors like GetPy_KeyKeys and similar.

File:

• fact/tools/pyscripts/pyfact/fits.h (modified) (38 diffs)

fact/tools/pyscripts/pyfact/fits.h

@@ -6 +6 @@
 #include <map>
 #include <string>
+#include <fstream>
 #include <sstream>
 #include <algorithm>
-
-#ifdef __EXCEPTIONS
 #include <stdexcept>
-#endif
-
-#ifdef __CINT__
+
+#define GCC_VERSION (__GNUC__ * 10000  + __GNUC_MINOR__ * 100  + __GNUC_PATCHLEVEL__)
+
+#ifndef __CINT__
+#include <unordered_map>
+#else
 #define off_t size_t
+namespace std
+{
+    template<class T, class S> class unordered_map<T, S>;
+}
 #endif
 
@@ -21 +27 @@
 #include <iomanip>
 #include <iostream>
-#define gLog cerr
-#define ___err___ ""
-#define ___all___ ""
+#ifndef gLog
+#define gLog std::cerr
+#define ___err___   ""
+#define ___warn___  ""
+#define ___all___   ""
+#endif
 #else
 #include "MLog.h"
 #include "MLogManip.h"
-#define ___err___ err
-#define ___all___ all
+#define ___err___   err
+#define ___warn___  warn
+#define ___all___   all
 #endif
 
 #define HAVE_ZLIB
-#ifdef HAVE_ZLIB
-#include "izstream.h"
+#if defined(HAVE_ZLIB) || defined(__CINT__)
+#include "extern_Mars_mcore/izstream.h"
 #else
 #include <fstream>
 #define izstream ifstream
-//#warning Support for zipped FITS files disabled.
-#endif
-
-#warning This file is considered DEPRICATED. Please use pyfits.h instead
-
-namespace std
-{
+#warning Support for zipped FITS files disabled.
+#endif
+
+#include "extern_Mars_mcore/checksum.h"
 
 class fits : public izstream
 {
 public:
+    //I know I know, you're going to yiell that this does not belong here.
+    //It will belong in the global scope eventually, and it makes the coding of zfits much simpler this way.
+    enum Compression_t
+    {
+        kCompUnknown,
+        kCompFACT
+    };
+
     struct Entry
     {
-        char   type;
-        string value;
-        string comment;
+        char        type;
+        std::string value;
+        std::string comment;
+        std::string fitsString;
 
         template<typename T>
@@ -59 +75 @@
             T t;
 
-            istringstream str(value);
+            std::istringstream str(value);
             str >> t;
 
@@ -70 +86 @@
         off_t offset;
 
-        string name;
+        bool is_compressed;
+
+        std::string name;
         size_t bytes_per_row;
         size_t num_rows;
         size_t num_cols;
+        size_t total_bytes; // NAXIS1*NAXIS2
 
         struct Column
@@ -80 +99 @@
             size_t num;
             size_t size;
+            size_t bytes;  // num*size
             char   type;
-            string unit;
+            std::string unit;
+            Compression_t comp;
         };
 
-        typedef map<string, Entry>  Keys;
-        typedef map<string, Column> Columns;
+        typedef std::map<std::string, Entry>  Keys;
+        typedef std::map<std::string, Column> Columns;
+        typedef std::vector<Column> SortedColumns;
 
         Columns cols;
+        SortedColumns sorted_cols;
         Keys    keys;
-
-        string Trim(const string &str, char c=' ') const
+        
+        //------------ vectors containing the same info as cols and keys------
+        // They are filled at the end of the contstuctor
+        // They are used py python, for retrieving the contents of the maps
+        // I did not find out how to return a map of <string, struct> via ROOT to python
+        // So I use this ugly hack. I guess its no problem, its just ugly
+        //
+        // btw: The way these variables are named, is imho as ugly as possible,
+        // but I will stick to it.
+        // usually a std::map contains key-value pairs
+        // that means the e.g. map<string, Entry> contains keys of type string
+        // and values of type Entry.
+        // but now the map itself was called 'keys', which means one should say:
+        // the keys of 'keys' are of type string
+        // and
+        // the value of 'keys' are of type Entry. 
+        // not to forget that:
+        // one field of a value of 'keys' is named 'value'
+        
+        // vectors for Keys keys;
+        vector<string> Py_KeyKeys;
+        vector<string> Py_KeyValues;
+        vector<string> Py_KeyComments;
+        vector<char>   Py_KeyTypes;
+        
+        // vectors for Columns cols;
+        vector<string> Py_ColumnKeys;
+        vector<size_t> Py_ColumnOffsets;
+        vector<size_t> Py_ColumnNums;
+        vector<size_t> Py_ColumnSizes;
+        vector<char>   Py_ColumnTypes;
+        vector<string> Py_ColumnUnits;
+        //-end of----- vectors containing the same info as cols and keys------
+        
+
+        int64_t datasum;
+
+        std::string Trim(const std::string &str, char c=' ') const
         {
             // Trim Both leading and trailing spaces
@@ -97 +156 @@
 
             // if all spaces or empty return an empty string
-            if (string::npos==pstart || string::npos==pend)
-                return string();
+            if (std::string::npos==pstart || std::string::npos==pend)
+                return std::string();
 
             return str.substr(pstart, pend-pstart+1);
         }
 
-        bool Check(const string &key, char type, const string &value="") const
+        bool Check(const std::string &key, char type, const std::string &value="") const
         {
             const Keys::const_iterator it = keys.find(key);
             if (it==keys.end())
             {
-                ostringstream str;
+                std::ostringstream str;
                 str << "Key '" << key << "' not found.";
 #ifdef __EXCEPTIONS
-                throw runtime_error(str.str());
-#else
-                gLog << ___err___ << "ERROR - " << str.str() << endl;
+                throw std::runtime_error(str.str());
+#else
+                gLog << ___err___ << "ERROR - " << str.str() << std::endl;
                 return false;
 #endif
@@ -120 +179 @@
             if (it->second.type!=type)
             {
-                ostringstream str;
+                std::ostringstream str;
                 str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << ".";
 #ifdef __EXCEPTIONS
-                throw runtime_error(str.str());
-#else
-                gLog << ___err___ << "ERROR - " << str.str() << endl;
+                throw std::runtime_error(str.str());
+#else
+                gLog << ___err___ << "ERROR - " << str.str() << std::endl;
                 return false;
 #endif
@@ -132 +191 @@
             if (!value.empty() && it->second.value!=value)
             {
-                ostringstream str;
+                std::ostringstream str;
                 str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << ".";
 #ifdef __EXCEPTIONS
-                throw runtime_error(str.str());
-#else
-                gLog << ___err___ << "ERROR - " << str.str() << endl;
+                throw std::runtime_error(str.str());
+#else
+                gLog << ___err___ << "ERROR - " << str.str() << std::endl;
                 return false;
 #endif
@@ -145 +204 @@
         }
 
-        Keys ParseBlock(const vector<string> &vec) const
-        {
-            map<string,Entry> rc;
+        Keys ParseBlock(const std::vector<std::string> &vec) const
+        {
+            std::map<std::string,Entry> rc;
 
             for (unsigned int i=0; i<vec.size(); i++)
             {
-                const string key = Trim(vec[i].substr(0,8));
+                const std::string key = Trim(vec[i].substr(0,8));
                 // Keywords without a value, like COMMENT / HISTORY
                 if (vec[i].substr(8,2)!="= ")
@@ -158 +217 @@
                 char type = 0;
 
-                string com;
-                string val = Trim(vec[i].substr(10));
+                std::string com;
+                std::string val = Trim(vec[i].substr(10));
+
                 if (val[0]=='\'')
                 {
@@ -167 +227 @@
                     {
                         const size_t pp = val.find_first_of('\'', p);
-                        if (pp==string::npos)
+                        if (pp==std::string::npos)
                             break;
 
@@ -177 +237 @@
 
                     // Set value, comment and type
-                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
+                    // comments could be just spaces. take care of this.
+                    if (ppp!=std::string::npos && val.size() != ppp+1)
+                        com = Trim(val.substr(ppp+1));
+
                     val  = Trim(val.substr(1, p-2));
                     type = 'T';
@@ -185 +248 @@
                     const size_t p = val.find_first_of('/');
 
-                    com = Trim(val.substr(p+2));
+                    if (val.size() != p+1)
+                        com = Trim(val.substr(p+2));
+
                     val = Trim(val.substr(0, p));
 
-                    if (val.empty() || val.find_first_of('T')!=string::npos || val.find_first_of('F')!=string::npos)
+                    if (val.empty() || val.find_first_of('T')!=std::string::npos || val.find_first_of('F')!=std::string::npos)
                         type = 'B';
                     else
-                        type = val.find_last_of('.')==string::npos ? 'I' : 'F';
+                        type = val.find_last_of('.')==std::string::npos ? 'I' : 'F';
                 }
 
-                const Entry e = { type, val, com };
+                const Entry e = { type, val, com, vec[i] };
+
                 rc[key] = e;
             }
@@ -201 +267 @@
         }
 
-        Table() : offset(0) { }
-        Table(const vector<string> &vec, off_t off) :
-            offset(off), keys(ParseBlock(vec))
-        {
+        Table() : offset(0), is_compressed(false) { }
+        Table(const std::vector<std::string> &vec, off_t off) : offset(off),
+            keys(ParseBlock(vec))
+        {
+            is_compressed = HasKey("ZTABLE") && Check("ZTABLE", 'B', "T");
+
             if (!Check("XTENSION", 'T', "BINTABLE") ||
                 !Check("NAXIS",    'I', "2")        ||
                 !Check("BITPIX",   'I', "8")        ||
-                !Check("PCOUNT",   'I', "0")        ||
                 !Check("GCOUNT",   'I', "1")        ||
                 !Check("EXTNAME",  'T')             ||
@@ -216 +283 @@
                 return;
 
-            bytes_per_row = Get<size_t>("NAXIS1");
-            num_rows      = Get<size_t>("NAXIS2");
+            if (is_compressed)
+            {
+                if (!Check("ZNAXIS1", 'I') ||
+                    !Check("ZNAXIS2", 'I') ||
+                    !Check("ZPCOUNT", 'I', "0"))
+                    return;
+            }
+            else
+            {
+                if (!Check("PCOUNT", 'I', "0"))
+                    return;
+            }
+
+            total_bytes   = Get<size_t>("NAXIS1")*Get<size_t>("NAXIS2");
+            bytes_per_row = is_compressed ? Get<size_t>("ZNAXIS1") : Get<size_t>("NAXIS1");
+            num_rows      = is_compressed ? Get<size_t>("ZNAXIS2") : Get<size_t>("NAXIS2");
             num_cols      = Get<size_t>("TFIELDS");
-
+            datasum       = is_compressed ? Get<int64_t>("DATASUM", -1) : Get<int64_t>("DATASUM", -1);
+//cout << "IS COMPRESSED =-========= " << is_compressed << " " << Get<size_t>("NAXIS1") << endl;
             size_t bytes = 0;
-            for (size_t i=1; i<=num_cols; i++)
-            {
-                ostringstream num;
-                num << i;
-
-                if (!Check("TTYPE"+num.str(), 'T') ||
-                    !Check("TFORM"+num.str(), 'T'))
+
+            const std::string tFormName = is_compressed ? "ZFORM" : "TFORM";
+            for (long long unsigned int i=1; i<=num_cols; i++)
+            {
+                const std::string num(std::to_string(i));
+
+                if (!Check("TTYPE"+num, 'T') ||
+                    !Check(tFormName+num, 'T'))
                     return;
 
-                const string id   = Get<string>("TTYPE"+num.str());
-                const string fmt  = Get<string>("TFORM"+num.str());
-                const string unit = Get<string>("TUNIT"+num.str(), "");
-
-                istringstream sin(fmt);
+                const std::string id   = Get<std::string>("TTYPE"+num);
+                const std::string fmt  = Get<std::string>(tFormName+num);
+                const std::string unit = Get<std::string>("TUNIT"+num, "");
+                const std::string comp = Get<std::string>("ZCTYP"+num, "");
+
+                Compression_t compress = kCompUnknown;
+                if (comp == "FACT")
+                    compress = kCompFACT;
+
+                std::istringstream sin(fmt);
                 int n = 0;
                 sin >> n;
@@ -247 +335 @@
                     // We could use negative values to mark floats
                     // otheriwse we could just cast them to int64_t?
-                case 'L': size = 1; break; // logical
-                // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits)
+                case 'L':                  // logical
+                case 'A':                  // char
                 case 'B': size = 1; break; // byte
                 case 'I': size = 2; break; // short
@@ -255 +343 @@
                 case 'E': size = 4; break; // float
                 case 'D': size = 8; break; // double
+                // case 'X': size =  n; break; // bits (n=number of bytes needed to contain all bits)
                 // case 'C': size =  8; break; // complex float
                 // case 'M': size = 16; break; // complex double
@@ -261 +350 @@
                 default:
                     {
-                        ostringstream str;
+                        std::ostringstream str;
                         str << "FITS format TFORM='" << fmt << "' not yet supported.";
 #ifdef __EXCEPTIONS
-                        throw runtime_error(str.str());
-#else
-                        gLog << ___err___ << "ERROR - " << str.str() << endl;
+                        throw std::runtime_error(str.str());
+#else
+                        gLog << ___err___ << "ERROR - " << str.str() << std::endl;
                         return;
 #endif
@@ -272 +361 @@
                 }
 
-                const Table::Column col = { bytes, n, size, type, unit };
+                const Table::Column col = { bytes, n, size, n*size, type, unit, compress};
 
                 cols[id] = col;
+                sorted_cols.push_back(col);
                 bytes  += n*size;
             }
@@ -281 +371 @@
             {
 #ifdef __EXCEPTIONS
-                throw runtime_error("Column size mismatch");
-#else
-                gLog << ___err___ << "ERROR - Column size mismatch" << endl;
+                throw std::runtime_error("Column size mismatch");
+#else
+                gLog << ___err___ << "ERROR - Column size mismatch" << std::endl;
                 return;
 #endif
             }
 
-            name = Get<string>("EXTNAME");
+            name = Get<std::string>("EXTNAME");
+            
+            Fill_Py_Vectors();
+        }
+        
+        void Fill_Py_Vectors(void)
+        {
+            // Fill the Py_ vectors see line: 90ff
+            // vectors for Keys keys;
+            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+            {
+                Py_KeyKeys.push_back(it->first);
+                Py_KeyValues.push_back(it->second.value);
+                Py_KeyComments.push_back(it->second.comment);
+                Py_KeyTypes.push_back(it->second.type);
+            }
+
+            // vectors for Columns cols;
+            for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++)
+            {
+                Py_ColumnKeys.push_back(it->first);
+                Py_ColumnTypes.push_back(it->second.type);
+                Py_ColumnOffsets.push_back(it->second.offset);
+                Py_ColumnNums.push_back(it->second.num);
+                Py_ColumnSizes.push_back(it->second.size);
+                Py_ColumnUnits.push_back(it->second.unit);
+            }
         }
 
         void PrintKeys(bool display_all=false) const
         {
-            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+            for (Keys::const_iterator it=keys.cbegin(); it!=keys.cend(); it++)
             {
                 if (!display_all &&
@@ -308 +424 @@
                     continue;
 
-                gLog << ___all___ << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl;
+                gLog << ___all___ << std::setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << std::endl;
             }}
 
         void PrintColumns() const
         {
-            typedef map<pair<size_t, string>, Column> Sorted;
+            typedef std::map<std::pair<size_t, std::string>, Column> Sorted;
 
             Sorted sorted;
 
-            for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++)
-                sorted[make_pair(it->second.offset, it->first)] = it->second;
-
-            for (Sorted::const_iterator it=sorted.begin(); it!=sorted.end(); it++)
-            {
-                gLog << ___all___ << setw(6) << it->second.offset << "| ";
+            for (Columns::const_iterator it=cols.cbegin(); it!=cols.cend(); it++)
+                sorted[std::make_pair(it->second.offset, it->first)] = it->second;
+
+            for (Sorted::const_iterator it=sorted.cbegin(); it!=sorted.cend(); it++)
+            {
+                gLog << ___all___ << std::setw(6) << it->second.offset << "| ";
                 gLog << it->second.num << 'x';
                 switch (it->second.type)
                 {
+                case 'A': gLog << "char(8)";    break;
                 case 'L': gLog << "bool(8)";    break;
                 case 'B': gLog << "byte(8)";    break;
@@ -334 +451 @@
                 case 'D': gLog << "double(64)"; break;
                 }
-                gLog << ": " << it->first.second << " [" << it->second.unit << "]" << endl;
+                gLog << ": " << it->first.second << " [" << it->second.unit << "]" << std::endl;
             }
         }
@@ -340 +457 @@
         operator bool() const { return !name.empty(); }
 
-        bool HasKey(const string &key) const
+        bool HasKey(const std::string &key) const
         {
             return keys.find(key)!=keys.end();
         }
 
-            // Values of keys are always signed
+        bool HasColumn(const std::string& col) const
+        {
+            return cols.find(col)!=cols.end();
+        }
+
+        const Columns &GetColumns() const
+        {
+            return cols;
+        }
+
+        const Keys &GetKeys() const
+        {
+            return keys;
+        }
+
+        // Values of keys are always signed
         template<typename T>
-            T Get(const string &key) const
-        {
-            const map<string,Entry>::const_iterator it = keys.find(key);
+            T Get(const std::string &key) const
+        {
+            const std::map<std::string,Entry>::const_iterator it = keys.find(key);
             if (it==keys.end())
             {
-                ostringstream str;
-                str << "Key '" << key << "' not found." << endl;
-#ifdef __EXCEPTIONS
-                throw runtime_error(str.str());
-#else
-                gLog << ___err___ << "ERROR - " << str.str() << endl;
+                std::ostringstream str;
+                str << "Key '" << key << "' not found.";
+#ifdef __EXCEPTIONS
+                throw std::runtime_error(str.str());
+#else
+                gLog << ___err___ << "ERROR - " << str.str() << std::endl;
+                return T();
+#endif
+            }
+            return it->second.Get<T>();
+        }
+
+        // Values of keys are always signed
+        template<typename T>
+            T Get(const std::string &key, const T &deflt) const
+        {
+            const std::map<std::string,Entry>::const_iterator it = keys.find(key);
+            return it==keys.end() ? deflt :it->second.Get<T>();
+        }
+
+        size_t GetN(const std::string &key) const
+        {
+            const Columns::const_iterator it = cols.find(key);
+            return it==cols.end() ? 0 : it->second.num;
+        }
+        
+
+
+        // There may be a gap between the main table and the start of the heap:
+        // this computes the offset
+        streamoff GetHeapShift() const
+        {
+            if (!HasKey("ZHEAPPTR"))
                 return 0;
-#endif
-            }
-            return it->second.Get<T>();
-        }
-
-            // Values of keys are always signed
-        template<typename T>
-            T Get(const string &key, const string &deflt) const
-        {
-            const map<string,Entry>::const_iterator it = keys.find(key);
-            return it==keys.end() ? deflt :it->second.Get<T>();
-        }
-
-        size_t GetN(const string &key) const
-        {
-            const Columns::const_iterator it = cols.find(key);
-            if (it==cols.end()){
-                ostringstream str;
-                str << "Key '" << key << "' not found." << endl;
-                str << "Possible keys are:" << endl;
-                for ( Columns::const_iterator it=cols.begin() ; it != cols.end(); ++it){
-                    str << it->first << endl;
-                }
-#ifdef __EXCEPTIONS
-                throw runtime_error(str.str());
-#else
-                gLog << ___err___ << "ERROR - " << str.str() << endl;
-                return 0;
-#endif
-            }
-            return it==cols.end() ? 0 : it->second.num;
+
+            const size_t shift = Get<size_t>("ZHEAPPTR");
+            return shift <= total_bytes ? 0 : shift - total_bytes;
+        }
+
+        // return total number of bytes 'all inclusive'
+        streamoff GetTotalBytes() const
+        {
+            //get offset of special data area from start of main table
+            const streamoff shift = GetHeapShift();
+
+            //and special data area size
+            const streamoff size  = HasKey("PCOUNT") ? Get<streamoff>("PCOUNT") : 0;
+
+            // Get the total size
+            const streamoff total = total_bytes + size + shift;
+
+            // check for padding
+            if (total%2880==0)
+                return total;
+
+            // padding necessary
+            return total + (2880 - (total%2880));
         }
     };
 
-private:
+protected:
+    std::ofstream fCopy;
+    std::vector<std::string> fListOfTables; // List of skipped tables. Last table is open table
+
     Table fTable;
 
-    typedef pair<void*, Table::Column> Address;
-    typedef vector<Address> Addresses;
+    typedef std::pair<void*, Table::Column> Address;
+    typedef std::vector<Address> Addresses;
     //map<void*, Table::Column> fAddresses;
     Addresses fAddresses;
 
-    vector<char> fBufferRow;
-    vector<char> fBufferDat;
+    typedef std::unordered_map<std::string, void*> Pointers;
+    Pointers fPointers;
+
+    std::vector<std::vector<char>> fGarbage;
+
+    std::vector<char> fBufferRow;
+    std::vector<char> fBufferDat;
 
     size_t fRow;
 
-    vector<string> ReadBlock(vector<string> &vec)
-    {
-        bool endtag = false;
+    Checksum fChkHeader;
+    Checksum fChkData;
+
+    bool ReadBlock(std::vector<std::string> &vec)
+    {
+        int endtag = 0;
         for (int i=0; i<36; i++)
         {
@@ -417 +579 @@
                 break;
 
-            if (c[0]==0)
-                return vector<string>();
-
-            string str(c);
+            fChkHeader.add(c, 80);
+
+//            if (c[0]==0)
+//                return vector<string>();
+
+            std::string str(c);
 
 //            if (!str.empty())
 //                cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl;
 
-            if (str=="END                                                                             ")
-            {
-                endtag = true;
-
-                // Make sure that no empty vector is returned
-                if (vec.size()%36==0)
-                    vec.push_back(string("END     = '' / "));
-            }
-
-            if (endtag)
+            if (endtag==2 || str=="END                                                                             ")
+            {
+                endtag = 2; // valid END tag found
                 continue;
+            }
+
+            if (endtag==1 || str=="                                                                                ")
+            {
+                endtag = 1; // end tag not found, but expected to be there
+                continue;
+            }
 
             vec.push_back(str);
         }
 
-        return vec;
-    }
-
-    string Compile(const string &key, int16_t i=-1)
-    {
-        if (i<0)
-            return key;
-
-        ostringstream str;
-        str << key << i;
-        return str.str();
-    }
-
-public:
-    fits(const string &fname) : izstream(fname.c_str())
+        // Make sure that no empty vector is returned
+        if (endtag && vec.size()%36==0)
+            vec.emplace_back("END     = '' / ");
+
+        return endtag==2;
+    }
+
+    std::string Compile(const std::string &key, int16_t i=-1) const
+    {
+#if GCC_VERSION < 40603
+        return i<0 ? key : key+std::to_string((long long int)(i));
+#else
+        return i<0 ? key : key+std::to_string(i);
+#endif
+    }
+
+    void Constructor(const std::string &fname, std::string fout="", const std::string& tableName="", bool force=false)
     {
         char simple[10];
@@ -463 +629 @@
         if (memcmp(simple, "SIMPLE  = ", 10))
         {
-            clear(rdstate()|ios::badbit);
-#ifdef __EXCEPTIONS
-            throw runtime_error("File is not a FITS file.");
-#else
-            gLog << ___err___ << "ERROR - File is not a FITS file." << endl;
+            clear(rdstate()|std::ios::badbit);
+#ifdef __EXCEPTIONS
+            throw std::runtime_error("File is not a FITS file.");
+#else
+            gLog << ___err___ << "ERROR - File is not a FITS file." << std::endl;
             return;
 #endif
@@ -476 +642 @@
         while (good())
         {
-            vector<string> block;
+            std::vector<std::string> block;
             while (1)
             {
+                // If we search for a table, we implicitly assume that
+                // not finding the table is not an error. The user
+                // can easily check that by eof() && !bad()
+                peek();
+                if (eof() && !bad() && !tableName.empty())
+                {
+                    break;
+                }
                 // FIXME: Set limit on memory consumption
-                ReadBlock(block);
+                const int rc = ReadBlock(block);
                 if (!good())
                 {
-                    clear(rdstate()|ios::badbit);
-#ifdef __EXCEPTIONS
-                    throw runtime_error("FITS file corrupted.");
-#else
-                    gLog << ___err___ << "ERROR - FITS file corrupted." << endl;
+                    clear(rdstate()|std::ios::badbit);
+#ifdef __EXCEPTIONS
+                    throw std::runtime_error("FITS file corrupted.");
+#else
+                    gLog << ___err___ << "ERROR - FITS file corrupted." << std::endl;
                     return;
 #endif
@@ -493 +667 @@
 
                 if (block.size()%36)
+                {
+                    if (!rc && !force)
+                    {
+                        clear(rdstate()|std::ios::badbit);
+#ifdef __EXCEPTIONS
+                        throw std::runtime_error("END keyword missing in FITS header.");
+#else
+                        gLog << ___err___ << "ERROR - END keyword missing in FITS file... file might be corrupted." << std::endl;
+                        return;
+#endif
+                    }
                     break;
-            }
-
-            if (block.size()==0)
+                }
+            }
+
+            if (block.empty())
                 break;
 
             if (block[0].substr(0, 9)=="SIMPLE  =")
+            {
+                fChkHeader.reset();
                 continue;
+            }
 
             if (block[0].substr(0, 9)=="XTENSION=")
             {
-                // FIXME: Check for table name
-
                 fTable = Table(block, tellg());
                 fRow   = (size_t)-1;
 
-                //fTable.PrintKeys();
-
                 if (!fTable)
                 {
-                    clear(rdstate()|ios::badbit);
+                    clear(rdstate()|std::ios::badbit);
                     return;
                 }
 
-                fBufferRow.resize(fTable.bytes_per_row);
+                const std::string &tname = fTable.Get<std::string>("EXTNAME");
+
+                fListOfTables.emplace_back(tname);
+
+                // Check for table name. Skip until eof or requested table are found.
+                // skip the current table?
+                if ((!tableName.empty() && tableName!=tname) || (tableName.empty() && "ZDrsCellOffsets"==tname))
+                {
+                    const streamoff skip = fTable.GetTotalBytes();
+                    seekg(skip, std::ios_base::cur);
+
+                    fChkHeader.reset();
+
+                    continue;
+                }
+
+                fBufferRow.resize(fTable.bytes_per_row + 8-fTable.bytes_per_row%4);
                 fBufferDat.resize(fTable.bytes_per_row);
 
-                /*
-                 // Next table should start at:
-                 const size_t size = fTable.bytes_per_row*fTable.num_rows;
-                 const size_t blks = size/(36*80);
-                 const size_t rest = size%(36*80);
-
-                seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
-                if (!good())
-                     gLog << ___err___ << "File seems to be incomplete (less data than expected from header)." << endl;
-
-                fRow   = fTable.num_rows;
-                 */
-
                 break;
             }
         }
-    }
-
-    void ReadRow(size_t row)
+
+        if (fout.empty())
+            return;
+
+        if (*fout.rbegin()=='/')
+        {
+            const size_t p = fname.find_last_of('/');
+            fout.append(fname.substr(p+1));
+        }
+
+        fCopy.open(fout);
+        if (!fCopy)
+        {
+            clear(rdstate()|std::ios::badbit);
+#ifdef __EXCEPTIONS
+            throw std::runtime_error("Could not open output file.");
+#else
+            gLog << ___err___ << "ERROR - Failed to open output file." << std::endl;
+            return;
+#endif
+        }
+
+        const streampos p = tellg();
+        seekg(0);
+
+        std::vector<char> buf(p);
+        read(buf.data(), p);
+
+        fCopy.write(buf.data(), p);
+        if (!fCopy)
+            clear(rdstate()|std::ios::badbit);
+    }
+
+public:
+    fits(const std::string &fname, const std::string& tableName="", bool force=false) : izstream(fname.c_str())
+    {
+        Constructor(fname, "", tableName, force);
+        if ((fTable.is_compressed ||fTable.name=="ZDrsCellOffsets") && !force)
+        {
+#ifdef __EXCEPTIONS
+            throw std::runtime_error("You are trying to read a compressed fits with the base fits class. Please use factfits instead.");
+#else
+            gLog << ___err___ << "ERROR - You are trying to read a compressed fits with the base fits class. Please use factfits instead." << std::endl;
+#endif
+            clear(rdstate()|std::ios::badbit);
+        }
+    }
+
+    fits(const std::string &fname, const std::string &fout, const std::string& tableName, bool force=false) : izstream(fname.c_str())
+    {
+        Constructor(fname, fout, tableName, force);
+        if ((fTable.is_compressed || fTable.name=="ZDrsCellOffsets") && !force)
+        {
+#ifdef __EXCEPTIONS
+            throw std::runtime_error("You are trying to read a compressed fits with the base fits class. Please use factfits instead.");
+#else
+            gLog << ___err___ << "ERROR - You are trying to read a compressed fits with the base fits class. Please use factfits instead." << std::endl;
+#endif
+            clear(rdstate()|std::ios::badbit);
+        }
+    }
+
+    fits() : izstream()
+    {
+
+    }
+
+    ~fits()
+    {
+        std::copy(std::istreambuf_iterator<char>(*this),
+                  std::istreambuf_iterator<char>(),
+                  std::ostreambuf_iterator<char>(fCopy));
+    }
+
+    virtual void StageRow(size_t row, char* dest)
     {
         // if (row!=fRow+1) // Fast seeking is ensured by izstream
         seekg(fTable.offset+row*fTable.bytes_per_row);
+        read(dest, fTable.bytes_per_row);
+        //fin.clear(fin.rdstate()&~ios::eofbit);
+    }
+
+    virtual void WriteRowToCopyFile(size_t row)
+    {
+        if (row==fRow+1)
+        {
+            const uint8_t offset = (row*fTable.bytes_per_row)%4;
+
+            fChkData.add(fBufferRow);
+            if (fCopy.is_open() && fCopy.good())
+                fCopy.write(fBufferRow.data()+offset, fTable.bytes_per_row);
+            if (!fCopy)
+                clear(rdstate()|std::ios::badbit);
+        }
+        else
+            if (fCopy.is_open())
+                clear(rdstate()|std::ios::badbit);
+    }
+
+    void ZeroBufferForChecksum(std::vector<char>& vec, const uint64_t extraZeros=0)
+    {
+        auto ib = vec.begin();
+        auto ie = vec.end();
+
+        *ib++ = 0;
+        *ib++ = 0;
+        *ib++ = 0;
+        *ib   = 0;
+
+        for (uint64_t i=0;i<extraZeros+8;i++)
+            *--ie = 0;
+    }
+
+    uint8_t ReadRow(size_t row)
+    {
+        // For the checksum we need everything to be correctly aligned
+        const uint8_t offset = (row*fTable.bytes_per_row)%4;
+
+        ZeroBufferForChecksum(fBufferRow);
+
+        StageRow(row, fBufferRow.data()+offset);
+
+        WriteRowToCopyFile(row);
 
         fRow = row;
 
-        read(fBufferRow.data(), fBufferRow.size());
-        //fin.clear(fin.rdstate()&~ios::eofbit);
+        return offset;
     }
 
     template<size_t N>
-        void revcpy(char *dest, const char *src, int num)
+        void revcpy(char *dest, const char *src, const int &num)
     {
         const char *pend = src + num*N;
         for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
-            reverse_copy(ptr, ptr+N, dest);
-    }
-
-    bool GetRow(size_t row)
-    {
-        ReadRow(row);
+            std::reverse_copy(ptr, ptr+N, dest);
+    }
+
+    virtual void MoveColumnDataToUserSpace(char *dest, const char *src, const Table::Column& c)
+    {
+        // Let the compiler do some optimization by
+        // knowing that we only have 1, 2, 4 and 8
+        switch (c.size)
+        {
+        case 1: memcpy   (dest, src, c.bytes); break;
+        case 2: revcpy<2>(dest, src, c.num);   break;
+        case 4: revcpy<4>(dest, src, c.num);   break;
+        case 8: revcpy<8>(dest, src, c.num);   break;
+        }
+    }
+
+    virtual bool GetRow(size_t row, bool check=true)
+    {
+        if (check && row>=fTable.num_rows)
+            return false;
+
+        const uint8_t offset = ReadRow(row);
         if (!good())
             return good();
 
-        for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++)
+        const char *ptr = fBufferRow.data() + offset;
+
+        for (Addresses::const_iterator it=fAddresses.cbegin(); it!=fAddresses.cend(); it++)
         {
             const Table::Column &c = it->second;
 
-            const char *src = fBufferRow.data() + c.offset;
+            const char *src = ptr + c.offset;
             char *dest = reinterpret_cast<char*>(it->first);
 
-            // Let the compiler do some optimization by
-            // knowing the we only have 1, 2, 4 and 8
-            switch (c.size)
-            {
-            case 1: memcpy   (dest, src, c.num*c.size); break;
-            case 2: revcpy<2>(dest, src, c.num);        break;
-            case 4: revcpy<4>(dest, src, c.num);        break;
-            case 8: revcpy<8>(dest, src, c.num);        break;
-            }
+            MoveColumnDataToUserSpace(dest, src, c);
         }
 
@@ -584 +898 @@
     }
 
-    bool GetNextRow()
-    {
-        return GetRow(fRow+1);
-    }
-
-    bool SkipNextRow()
+    bool GetNextRow(bool check=true)
+    {
+        return GetRow(fRow+1, check);
+    }
+
+    virtual bool SkipNextRow()
     {
         seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
@@ -600 +914 @@
     }
 
+    template<class T, class S>
+    const T &GetAs(const std::string &name)
+    {
+        return *reinterpret_cast<S*>(fPointers[name]);
+    }
+
+    void *SetPtrAddress(const std::string &name)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            std::ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Column not found.";
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
+            return NULL;
+#endif
+        }
+
+        Pointers::const_iterator it = fPointers.find(name);
+        if (it!=fPointers.end())
+            return it->second;
+
+        fGarbage.emplace_back(fTable.cols[name].bytes);
+
+        void *ptr = fGarbage.back().data();
+
+        fPointers[name] = ptr;
+        fAddresses.emplace_back(ptr, fTable.cols[name]);
+        sort(fAddresses.begin(), fAddresses.end(), Compare);
+        return ptr;
+    }
+
     template<typename T>
-    bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
+    bool SetPtrAddress(const std::string &name, T *ptr, size_t cnt)
     {
         if (fTable.cols.count(name)==0)
         {
-            ostringstream str;
-            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
-#ifdef __EXCEPTIONS
-            throw runtime_error(str.str());
-#else
-            gLog << ___err___ << "ERROR - " << str.str() << endl;
+            std::ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Column not found.";
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
             return false;
 #endif
@@ -617 +965 @@
         if (sizeof(T)!=fTable.cols[name].size)
         {
-            ostringstream str;
+            std::ostringstream str;
             str << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
-                << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
-#ifdef __EXCEPTIONS
-            throw runtime_error(str.str());
-#else
-            gLog << ___err___ << "ERROR - " << str.str() << endl;
+                << fTable.cols[name].size << " from header, got " << sizeof(T);
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
             return false;
 #endif
@@ -630 +978 @@
         if (cnt!=fTable.cols[name].num)
         {
-            ostringstream str;
+            std::ostringstream str;
             str << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
-                << fTable.cols[name].num << " from header, got " << cnt << endl;
-#ifdef __EXCEPTIONS
-            throw runtime_error(str.str());
-#else
-            gLog << ___err___ << "ERROR - " << str.str() << endl;
+                << fTable.cols[name].num << " from header, got " << cnt;
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
             return false;
 #endif
@@ -645 +993 @@
 
         //fAddresses[ptr] = fTable.cols[name];
-        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        fPointers[name] = ptr;
+        fAddresses.emplace_back(ptr, fTable.cols[name]);
         sort(fAddresses.begin(), fAddresses.end(), Compare);
         return true;
@@ -651 +1000 @@
 
     template<class T>
-    bool SetRefAddress(const string &name, T &ptr)
+    bool SetRefAddress(const std::string &name, T &ptr)
     {
         return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
@@ -657 +1006 @@
 
     template<typename T>
-    bool SetVecAddress(const string &name, vector<T> &vec)
+    bool SetVecAddress(const std::string &name, std::vector<T> &vec)
     {
         return SetPtrAddress(name, vec.data(), vec.size());
@@ -663 +1012 @@
 
     template<typename T>
-        T Get(const string &key) const
+        T Get(const std::string &key) const
     {
         return fTable.Get<T>(key);
@@ -669 +1018 @@
 
     template<typename T>
-        T Get(const string &key, const string &deflt) const
+        T Get(const std::string &key, const std::string &deflt) const
     {
         return fTable.Get<T>(key, deflt);
     }
 
-    bool SetPtrAddress(const string &name, void *ptr)
+    bool SetPtrAddress(const std::string &name, void *ptr, size_t cnt=0)
     {
         if (fTable.cols.count(name)==0)
         {
-            ostringstream str;
-            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
-#ifdef __EXCEPTIONS
-            throw runtime_error(str.str());
-#else
-            gLog << ___err___ << "ERROR - " << str.str() << endl;
+            std::ostringstream str;
+            str <<"SetPtrAddress('" << name << "') - Column not found.";
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
+            return false;
+#endif
+        }
+
+        if (cnt && cnt!=fTable.cols[name].num)
+        {
+            std::ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
+                << fTable.cols[name].num << " from header, got " << cnt;
+#ifdef __EXCEPTIONS
+            throw std::runtime_error(str.str());
+#else
+            gLog << ___err___ << "ERROR - " << str.str() << std::endl;
             return false;
 #endif
@@ -692 +1054 @@
 
         //fAddresses[ptr] = fTable.cols[name];
-        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        fPointers[name] = ptr;
+        fAddresses.emplace_back(ptr, fTable.cols[name]);
         sort(fAddresses.begin(), fAddresses.end(), Compare);
         return true;
     }
 
-    bool     HasKey(const string &key) const { return fTable.HasKey(key); }
-    int64_t  GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
-    uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
-    double   GetFloat(const string &key) const { return fTable.Get<double>(key); }
-    string   GetStr(const string &key) const { return fTable.Get<string>(key); }
-
-    size_t GetN(const string &key) const
+    bool     HasKey(const std::string &key) const { return fTable.HasKey(key); }
+    bool     HasColumn(const std::string& col) const { return fTable.HasColumn(col);}
+    const Table::Columns &GetColumns() const { return fTable.GetColumns();}
+    const Table::SortedColumns& GetSortedColumns() const { return fTable.sorted_cols;}
+    const Table::Keys &GetKeys() const { return fTable.GetKeys();}
+
+    int64_t     GetInt(const std::string &key) const { return fTable.Get<int64_t>(key); }
+    uint64_t    GetUInt(const std::string &key) const { return fTable.Get<uint64_t>(key); }
+    double      GetFloat(const std::string &key) const { return fTable.Get<double>(key); }
+    std::string GetStr(const std::string &key) const { return fTable.Get<std::string>(key); }
+
+    size_t GetN(const std::string &key) const
     {
         return fTable.GetN(key);
     }
 
-    size_t GetNumRows() const { return fTable.num_rows; }
-    size_t GetRow() const { return fRow; }
+//    size_t GetNumRows() const { return fTable.num_rows; }
+    size_t GetRow() const { return fRow==(size_t)-1 ? 0 : fRow; }
 
     operator bool() const { return fTable && fTable.offset!=0; }
@@ -715 +1083 @@
     void PrintKeys() const { fTable.PrintKeys(); }
     void PrintColumns() const { fTable.PrintColumns(); }
+    
+    // 'Wrappers' for the Table's Getters of the Py Vectors
+    vector<string> GetPy_KeyKeys()       {return fTable.Py_KeyKeys; }
+    vector<string> GetPy_KeyValues()     {return fTable.Py_KeyValues; }
+    vector<string> GetPy_KeyComments()   {return fTable.Py_KeyComments; }
+    vector<char>   GetPy_KeyTypes()      {return fTable.Py_KeyTypes; }
+    
+    vector<string> GetPy_ColumnKeys()    {return fTable.Py_ColumnKeys; }
+    vector<size_t> GetPy_ColumnOffsets() {return fTable.Py_ColumnOffsets; }
+    vector<size_t> GetPy_ColumnNums()    {return fTable.Py_ColumnNums; }
+    vector<size_t> GetPy_ColumnSizes()   {return fTable.Py_ColumnSizes;}
+    vector<char>   GetPy_ColumnTypes()   {return fTable.Py_ColumnTypes;}
+    vector<string> GetPy_ColumnUnits()   {return fTable.Py_ColumnUnits;}
+
+    bool IsHeaderOk() const { return fTable.datasum<0?false:(fChkHeader+Checksum(fTable.datasum)).valid(); }
+    virtual bool IsFileOk() const { return (fChkHeader+fChkData).valid(); }
+
+    bool IsCompressedFITS() const { return fTable.is_compressed;}
+
+    virtual size_t GetNumRows() const
+    {
+        return fTable.Get<size_t>("NAXIS2");
+    }
+
+    virtual size_t GetBytesPerRow() const
+    {
+        return fTable.Get<size_t>("NAXIS1");
+    }
+
+    const std::vector<std::string> &GetTables() const
+    {
+        return fListOfTables;
+    }
 };
 
-};
-#endif
+#endif