Index: trunk/Mars/mbase/MFits.h
===================================================================
--- trunk/Mars/mbase/MFits.h	(revision 11436)
+++ trunk/Mars/mbase/MFits.h	(revision 11436)
@@ -0,0 +1,532 @@
+#ifndef MARS_MFits
+#define MARS_MFits
+
+#include <map>
+#include <string>
+#include <sstream>
+#include <algorithm>
+
+#include <stdint.h>
+
+#ifdef __CINT__
+#define off_t size_t
+#endif
+
+#ifndef __MARS__
+#include <vector>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#define MZlib ifstream
+#define gLog cerr
+#define err ""
+#define all ""
+#else
+#include "MZlib.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+#endif
+
+using namespace std;
+
+class MFits : public MZlib
+{
+    struct Entry
+    {
+        char   type;
+        string value;
+        string comment;
+
+        template<typename T>
+            T Get() const
+        {
+            T t;
+
+            istringstream str(value);
+            str >> t;
+
+            return t;
+        }
+    };
+
+    vector<string> ReadBlock(vector<string> &vec)
+    {
+        bool end = false;
+        for (int i=0; i<36; i++)
+        {
+            char c[81];
+            c[80] = 0;
+            read(c, 80);
+            if (!good())
+                break;
+
+            if (c[0]==0)
+                return vector<string>();
+
+            string str(c);
+
+            //if (!str.empty())
+            //    cout << setw(2) << i << "|" << str << "|" << endl;
+
+            if (str=="END                                                                             ")
+                end = true;
+
+            if (end)
+                continue;
+
+            vec.push_back(str);
+        }
+
+        return vec;
+    }
+
+public:
+    struct Table
+    {
+        off_t offset;
+
+        string name;
+        size_t bytes_per_row;
+        size_t num_rows;
+        size_t num_cols;
+
+        struct Column
+        {
+            size_t offset;
+            size_t num;
+            size_t size;
+            char   type;
+        };
+
+        typedef map<string, Entry>  Keys;
+        typedef map<string, Column> Columns;
+
+        Columns cols;
+        Keys    keys;
+
+        string Trim(const string &str, char c=' ') const
+        {
+            // Trim Both leading and trailing spaces
+            const size_t start = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces
+            const size_t end   = str.find_last_not_of(c);  // Find the first character position from reverse af
+
+            // if all spaces or empty return an empty string
+            if (string::npos==start || string::npos==end)
+                return string();
+
+            return str.substr(start, end-start+1);
+        }
+
+        bool Check(const string &key, char type, const string &value="") const
+        {
+            const Keys::const_iterator it = keys.find(key);
+            if (it==keys.end())
+            {
+                gLog << err << "ERROR - Key '" << key << "' not found." << endl;
+                return false;
+            }
+
+            if (it->second.type!=type)
+            {
+                gLog << err << "ERROR - Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << "." << endl;
+                return false;
+            }
+
+            if (!value.empty() && it->second.value!=value)
+            {
+                gLog << err << "ERROR - Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << "." << endl;
+                return false;
+            }
+
+            return true;
+        }
+
+        Keys ParseBlock(const vector<string> &vec) const
+        {
+            map<string,Entry> rc;
+
+            for (unsigned int i=0; i<vec.size(); i++)
+            {
+                const string key = Trim(vec[i].substr(0,8));
+                // Keywords without a value, like COMMENT / HISTORY
+                if (vec[i].substr(8,2)!="= ")
+                    continue;
+
+                char type = 0;
+
+                string com;
+                string val = Trim(vec[i].substr(10));
+                if (val[0]=='\'')
+                {
+                    // First skip all '' in the string
+                    size_t p = 1;
+                    while (1)
+                    {
+                        const size_t pp = val.find_first_of('\'', p);
+                        if (pp==string::npos)
+                            break;
+
+                        p = val[pp+1]=='\'' ? pp+2 : pp+1;
+                    }
+
+                    // Now find the comment
+                    const size_t ppp = val.find_first_of('/', p);
+
+                    // Set value, comment and type
+                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
+                    val  = Trim(val.substr(1, p-2));
+                    type = 'T';
+                }
+                else
+                {
+                    const size_t p = val.find_first_of('/');
+
+                    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)
+                        type = 'B';
+                    else
+                        type = val.find_last_of('.')==string::npos ? 'I' : 'F';
+                }
+
+                const Entry e = { type, val, com };
+                rc[key] = e;
+            }
+
+            return rc;
+        }
+
+        Table() : offset(0) { }
+        Table(const vector<string> &vec, off_t off) :
+            offset(off), keys(ParseBlock(vec))
+        {
+            if (!Check("XTENSION", 'T', "BINTABLE") ||
+                !Check("NAXIS",    'I', "2")        ||
+                !Check("BITPIX",   'I', "8")        ||
+                !Check("PCOUNT",   'I', "0")        ||
+                !Check("GCOUNT",   'I', "1")        ||
+                !Check("EXTNAME",  'T')             ||
+                !Check("NAXIS1",   'I')             ||
+                !Check("NAXIS2",   'I')             ||
+                !Check("TFIELDS",  'I'))
+                return;
+
+            bytes_per_row = Get<size_t>("NAXIS1");
+            num_rows      = Get<size_t>("NAXIS2");
+            num_cols      = Get<size_t>("TFIELDS");
+
+            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'))
+                    return;
+
+                const string name = Get<string>("TTYPE"+num.str());
+                const string fmt  = Get<string>("TFORM"+num.str());
+
+                istringstream in(fmt);
+                int n = 0;
+                in >> n;
+                if (!in)
+                    n = 1;
+
+                const char type = fmt[fmt.length()-1];
+
+                size_t size = 0;
+                switch (type)
+                {
+                    // We could use negative values to mark floats
+                    // otheriwse we could just cast them to int64_t?
+                case 'B': size = 1; break; // byte
+                case 'I': size = 2; break; // short
+                case 'J': size = 4; break; // int
+                case 'K': size = 8; break; // long long
+                case 'E': size = 4; break; // float
+                case 'D': size = 8; break; // double
+                default:
+                    gLog << err << "FITS format '" << fmt << "' not yet supported." << endl;
+                    continue;
+                }
+
+                const Table::Column col = { bytes, n, size, type };
+
+                cols[name] = col;
+                bytes  += n*size;
+            }
+
+            if (bytes!=bytes_per_row)
+            {
+                gLog << err << "Column size mismatch" << endl;
+                return;
+            }
+
+            name = Get<string>("EXTNAME");
+        }
+
+        void PrintKeys() const
+        {
+            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+                gLog << all << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl;
+        }
+
+        void PrintColumns() const
+        {
+            typedef map<pair<size_t, 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 << "| " << it->first.second << " [" << it->second.type << '=' << it->second.num << 'x' << it->second.size << "byte]" << endl;
+        }
+
+        operator bool() const { return !name.empty(); }
+
+        bool HasKey(const string &key) const
+        {
+            return keys.find(key)!=keys.end();
+        }
+
+            // 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);
+            if (it==keys.end())
+            {
+                gLog << err << "Key '" << key << "' not found." << endl;
+                return 0;
+            }
+            return it->second.Get<T>();
+        }
+
+        size_t GetN(const string &key) const
+        {
+            const Columns::const_iterator it = cols.find(key);
+            return it==cols.end() ? 0 : it->second.num;
+        }
+    };
+
+private:
+    Table fTable;
+
+    map<void*, Table::Column> fAddresses;
+
+    vector<char> fBufferRow;
+    vector<char> fBufferDat;
+
+    size_t fRow;
+
+    string Compile(const string &key, int16_t i=-1)
+    {
+        if (i<0)
+            return key;
+
+        ostringstream str;
+        str << key << i;
+        return str.str();
+    }
+
+public:
+    MFits(const string &fname) : MZlib(fname.c_str())
+    {
+        while (good())
+        {
+            vector<string> block;
+            while (1)
+            {
+                // FIXME: Set limit on memory consumption
+                ReadBlock(block);
+                if (!good())
+                {
+                    gLog << err << "FITS file corrupted." << endl;
+                    clear(rdstate()|ios::badbit);
+                    return;
+                }
+
+                if (block.size()%36)
+                    break;
+            }
+
+            if (block.size()==0)
+                break;
+
+            if (block[0].substr(0, 9)=="SIMPLE  =")
+                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)
+                    return;
+
+                if (!fTable.Check("TELESCOP", 'T', "FACT"))
+                {
+                    gLog << err << "Not a valid FACT fits file." << endl;
+                    return;
+                }
+
+                fBufferRow.resize(fTable.bytes_per_row);
+                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 (row!=fRow+1) // Fast seeking is ensured by MZlib
+        seekg(fTable.offset+row*fTable.bytes_per_row);
+
+        fRow = row;
+
+        read(fBufferRow.data(), fBufferRow.size());
+        //fin.clear(fin.rdstate()&~ios::eofbit);
+    }
+
+    bool GetRow(size_t row)
+    {
+        ReadRow(row);
+        if (!good())
+            return good();
+
+        for (map<void*,Table::Column>::const_iterator it=fAddresses.begin();
+             it!=fAddresses.end(); it++)
+        {
+            const Table::Column &c = it->second;
+
+            const char *src = fBufferRow.data() + c.offset;
+            const char *end = src + c.num*c.size;
+
+            char *dest = reinterpret_cast<char*>(it->first);
+            for (const char *ptr = src; ptr<end; ptr+=c.size, dest+=c.size)
+                reverse_copy(ptr, ptr+c.size, dest);
+        }
+
+        return good();
+    }
+
+    bool GetNextRow()
+    {
+        return GetRow(fRow+1);
+    }
+
+    bool SkipNextRow()
+    {
+        seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
+        return good();
+    }
+
+    template<typename T>
+    bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+            return false;
+        }
+
+        if (sizeof(T)!=fTable.cols[name].size)
+        {
+            gLog << err << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
+                << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
+
+            return false;
+        }
+
+        if (cnt!=fTable.cols[name].num)
+        {
+            gLog << err << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
+                << fTable.cols[name].num << " from header, got " << cnt << endl;
+
+            return false;
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        fAddresses[ptr] = fTable.cols[name];
+
+        return true;
+    }
+
+    template<class T>
+    bool SetRefAddress(const string &name, T &ptr)
+    {
+        return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
+    }
+
+    template<typename T>
+    bool SetVecAddress(const string &name, vector<T> &vec)
+    {
+        return SetPtrAddress(name, vec.data(), vec.size());
+    }
+
+    template<typename T>
+        T Get(const string &key) const
+    {
+        return fTable.Get<T>(key);
+    }
+
+    bool SetPtrAddress(const string &name, void *ptr)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+            return false;
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        fAddresses[ptr] = fTable.cols[name];
+
+        return true;
+    }
+
+    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); }
+
+    size_t GetN(const string &key) const
+    {
+        return fTable.GetN(key);
+    }
+
+    size_t GetNumRows() const { return fTable.num_rows; }
+    size_t GetRow() const { return fRow; }
+
+    operator bool() const { return fTable && fTable.offset!=0; }
+
+    void PrintKeys() const { fTable.PrintKeys(); }
+    void PrintColumns() const { fTable.PrintColumns(); }
+};
+#endif