Changeset 12723

Timestamp:

12/15/11 10:39:26 (13 years ago)

Author:

lyard

Message:

replaced cfitsio by fits.h

File:

• trunk/FACT++/src/fitsdump.cc (modified) (35 diffs)

trunk/FACT++/src/fitsdump.cc

@@ -11 +11 @@
 #include <fstream>
 
-#include <CCfits/CCfits>
-
-//#define PLOTTING_PLEASE
+#include "externals/fits.h"
+
+#define PLOTTING_PLEASE
 
 #ifdef PLOTTING_PLEASE
@@ -32 +32 @@
 using namespace std;
 
-class MyColumn : public CCfits::Column
-{
-public:
-    const string &comment() const { return CCfits::Column::comment(); }
-    long width() const
-    {//the width() method seems to be buggy (or empty ?) as it always return 1... redo it ourselves.
-        string inter = format();
-        inter = inter.substr(0, inter.size()-1);
-        return atoi(inter.c_str());
-    }
-};
 
 #ifdef PLOTTING_PLEASE
@@ -68 +57 @@
 
 private:
-    
-    CCfits::FITS  *fFile;                 /// FITS pointer
-    CCfits::Table *fTable;                /// Table pointer
-    map<string, MyColumn*> fColMap; /// map between the column names and their CCfits objects
+    bool fDotsPlease;
+    fits* fFile;
+    string fFilename;
+
+    fits::Table::Columns fColMap;
+    fits::Table::Keys fKeyMap;
 
     // Convert CCfits::ValueType into a human readable string
-    string ValueTypeToStr(CCfits::ValueType type) const;
+    string ValueTypeToStr(char type) const;
 
     // Convert CCfits::ValueType into a number of associated bytes
-    int    ValueTypeToSize(CCfits::ValueType type) const;
+    int    ValueTypeToSize(char type) const;
 
     /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
-    vector<int> CalculateOffsets() const;
+//    vector<int> CalculateOffsets() const;
 
     template<class T>
@@ -89 +80 @@
 
     /// Write a single row of the selected data
-    int  WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
+ //   int  WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
 
     bool OpenFile(const string &);        /// Open a file
@@ -108 +99 @@
     int doCurvesDisplay( const vector<string> &list, const string& tableName);
 #endif
+    int doStatsPlease(const string &filename, const vector<string>& list, int precision);
     //    bool Plot(const vector<string> &list);
 
@@ -121 +113 @@
 //!        the ostream where to redirect the outputs
 //
-FitsDumper::FitsDumper() : fFile(0), fTable(0)
+FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false)
 {
 }
@@ -136 +128 @@
 
 
-string FitsDumper::ValueTypeToStr(CCfits::ValueType type) const
+string FitsDumper::ValueTypeToStr(char type) const
 {
     switch (type)
     {
-    case CCfits::Tbyte:     return "uint8_t";
-    case CCfits::Tushort:   return "uint16_t";
-    case CCfits::Tshort:    return "int16_t";
-    case CCfits::Tuint:     return "uint32_t";
-    case CCfits::Tint:      return "int32_t";
-    case CCfits::Tulong:    return "uint32_t";
-    case CCfits::Tlong:     return "int32_t";
-    case CCfits::Tlonglong: return "int64_t";
-    case CCfits::Tfloat:    return "float";
-    case CCfits::Tdouble:   return "double";
-
+        case 'L': return "bool(8)";
+        case 'B': return "byte(8)";
+        case 'I': return "short(16)";
+        case 'J': return "int(32)";
+        case 'K': return "int(64)";
+        case 'E': return "float(32)";
+        case 'D': return "double(64)";
     default:
         return "unknwown";
@@ -156 +144 @@
 }
 
-int FitsDumper::ValueTypeToSize(CCfits::ValueType type) const
+int FitsDumper::ValueTypeToSize(char type) const
 {
     switch (type)
     {
-    case CCfits::Tbyte:     return sizeof(uint8_t);
-    case CCfits::Tushort:
-    case CCfits::Tshort:    return sizeof(uint16_t);
-    case CCfits::Tuint:
-    case CCfits::Tint:
-    case CCfits::Tulong:
-    case CCfits::Tlong:     return sizeof(uint32_t);
-    case CCfits::Tlonglong: return sizeof(uint64_t);
-    case CCfits::Tfloat:    return sizeof(float);
-    case CCfits::Tdouble:   return sizeof(double);
-
+        case 'L': return sizeof(uint8_t);
+        case 'B': return sizeof(int8_t);
+        case 'I': return sizeof(int16_t);
+        case 'J': return sizeof(int32_t);
+        case 'K': return sizeof(int64_t);
+        case 'E': return sizeof(float);
+        case 'D': return sizeof(double);
     default:
         return 0;
@@ -185 +169 @@
     return t;
 }
-/*
-template<class T>
-double FitsDumper::PtrToDouble(const unsigned char *ptr) const
-{
-    T t;
-    reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
-    return t;
-}
-
-double FitsDumper::PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const
-{
-    switch (type)
-    {
-    case CCfits::Tbyte:     return PtrToDouble<uint8_t> (ptr);
-    case CCfits::Tushort:   return PtrToDouble<uint16_t>(ptr);
-    case CCfits::Tuint:
-    case CCfits::Tulong:    return PtrToDouble<uint32_t>(ptr);
-    case CCfits::Tshort:    return PtrToDouble<int16_t> (ptr);
-    case CCfits::Tint:
-    case CCfits::Tlong:     return PtrToDouble<int32_t> (ptr);
-    case CCfits::Tlonglong: return PtrToDouble<int64_t> (ptr);
-    case CCfits::Tfloat:    return PtrToDouble<float>   (ptr);
-    case CCfits::Tdouble:   return PtrToDouble<double>  (ptr);
-
-    default:
-        throw runtime_error("Data type not implemented yet.");
-    }
-}
-*/
-
-// --------------------------------------------------------------------------
-//
-//! Writes a single row of the selected FITS data to the output file.
-//!
-//! Data type \b not yet implemented:
-//!   CCfits::Tnull, CCfits::Tbit, CCfits::Tlogical, CCfits::Tstring,
-//!   CCfits::Tcomplex, CCfits::Tdblcomplex, CCfits::VTbit,
-//!   CCfits::VTbyte, CCfits::VTlogical, CCfits::VTushort,
-//!   CCfits::VTshort, CCfits::VTuint, CCfits::VTint, CCfits::VTulong,
-//!   CCfits::VTlong, CCfits::VTlonglong, CCfits::VTfloat,
-//!   CCfits::VTdouble, CCfits::VTcomplex, CCfits::VTdblcomplex
-//!
-//! @param offsets
-//!         a vector containing the offsets to the columns (in bytes)
-//! @param targetFile
-//!         the ofstream where to write to
-//! @param fitsBuffer
-//!         the memory were the row has been loaded by cfitsio
-//
-int FitsDumper::WriteRow(ostream &out, const vector<MyColumn*> &list, const vector<int> &offsets, unsigned char* fitsBuffer, const vector<pair<int, int> >& ranges) const
-{
-    int cnt = 0;
-    vector<pair<int, int> >::const_iterator jt = ranges.begin();
-    for (vector<MyColumn*>::const_iterator it=list.begin(); it!=list.end(); it++, jt++)
-    {
-        if (jt == ranges.end())
-        {
-            cout << "ERROR: END OF RANGE POINTER REACHED" << endl;
-            return false;
-        }
-        const MyColumn *col = *it;
-
-        // CCfits starts counting at 1 not 0
-        const int offset = offsets[col->index()-1];
-
-        // Get the pointer to the array which we are supposed to print
-        const unsigned char *ptr = fitsBuffer + offset;
-
-        // Loop over all array entries
-        int sizeToSkip = 0;
-        switch (col->type())
-        {
-        case CCfits::Tbyte:     sizeToSkip = sizeof(uint8_t); break;
-        case CCfits::Tushort:   sizeToSkip = sizeof(uint16_t); break;
-        case CCfits::Tuint:
-        case CCfits::Tulong:    sizeToSkip = sizeof(uint32_t); break;
-        case CCfits::Tshort:    sizeToSkip = sizeof(int16_t); break;
-        case CCfits::Tint:
-        case CCfits::Tlong:     sizeToSkip = sizeof(int32_t); break;
-        case CCfits::Tlonglong: sizeToSkip = sizeof(int64_t); break;
-        case CCfits::Tfloat:    sizeToSkip = sizeof(float); break;
-        case CCfits::Tdouble:   sizeToSkip = sizeof(double); break;
-        default:
-            cerr << "Data type not implemented yet." << endl;
-            return 0;
-        }
-        ptr += sizeToSkip*jt->first;
-        for (int width=jt->first; width<jt->second; width++)
-        {
-            switch (col->type())
-            {
-            case CCfits::Tbyte:     out << PtrToValue<uint8_t> (ptr); break;
-            case CCfits::Tushort:   out << PtrToValue<uint16_t>(ptr); break;
-            case CCfits::Tuint:
-            case CCfits::Tulong:    out << PtrToValue<uint32_t>(ptr); break;
-            case CCfits::Tshort:    out << PtrToValue<int16_t> (ptr); break;
-            case CCfits::Tint:
-            case CCfits::Tlong:     out << PtrToValue<int32_t> (ptr); break;
-            case CCfits::Tlonglong: out << PtrToValue<int64_t> (ptr); break;
-            case CCfits::Tfloat:    out << PtrToValue<float>   (ptr); break;
-            case CCfits::Tdouble:   out << PtrToValue<double>  (ptr); break;
-
-            default:
-                cerr << "Data type not implemented yet." << endl;
-                return 0;
-            }
-
-            out << " ";
-            cnt++;
-        }
-    }
-
-    if (cnt>0)
-        out << endl;
-
-    if (out.fail())
-    {
-        cerr << "ERROR - writing output: " << strerror(errno) << endl;
-        return -1;
-    }
-
-    return cnt;
-}
-
-// --------------------------------------------------------------------------
-//
-//! Calculates the required buffer size for reading one row of the current
-//! table. Also calculates the offsets to all the columns
-//
-vector<int> FitsDumper::CalculateOffsets() const
-{
-    map<int,int> sizes;
-
-    for (map<string, MyColumn*>::const_iterator it=fColMap.begin();
-         it!=fColMap.end(); it++)
-    {
-        const int &width = it->second->width();
-        const int &idx   = it->second->index();
-
-        const int size = ValueTypeToSize(it->second->type());
-        if (size==0)
-        {
-            cerr << "Data type " << (int)it->second->type() << " not implemented yet." << endl;
-            return vector<int>();
-        }
-//cout << "data size: " << size << " width: " << width << " index: " << idx << endl;
-        int realwidth = (width == 0)? 1 : width;
-        sizes[idx] = size*realwidth;
-    }
-
-    //calculate the offsets in the vector.
-    vector<int> result(1, 0);
-
-    int size = 0;
-    int idx  = 0;
-
-    for (map<int,int>::const_iterator it=sizes.begin(); it!=sizes.end(); it++)
-    {
-        size += it->second;
-        result.push_back(size);
-
-        if (it->first == ++idx)
-            continue;
-
-        cerr << "Expected index " << idx << ", but found " << it->first << endl;
-        return vector<int>();
-    }
-
-    return result;
-}
-
 // --------------------------------------------------------------------------
 //
@@ -363 +176 @@
 {
     if (fFile)
+    {
+        fFile->close();
         delete fFile;
-
-    ostringstream str;
-    try
-    {
-        fFile = new CCfits::FITS(filename);
-    }
-    catch (CCfits::FitsException e)
-    {
-        cerr << "Could not open FITS file " << filename << " reason: " << e.message() << endl;
+    }
+
+    try {
+        fFile = new fits(filename);
+    }
+    catch (std::runtime_error e)
+    {
+        cout << "Something went wrong while trying to open " << filename;
+        cout << ": " << e.what() << " Aborting dump." << endl;
         return false;
     }
+    fFilename = filename;
+
+    const fits::Table::Columns& tCols = fFile->getColumns();
+
+    for (auto it=tCols.begin(); it != tCols.end(); it++)
+        fColMap.insert(*it);
+
+    const fits::Table::Keys& tkeys = fFile->getKeys();
+
+    for (auto it=tkeys.begin(); it != tkeys.end(); it++)
+        fKeyMap.insert(*it);
 
     return true;
@@ -387 +213 @@
     }
 
-    const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension();
-    if (extMap.find(tablename) == extMap.end())
-    {
-        cerr << "Table '" << tablename << "' not found." << endl;
-        return false;
-    }
-
-    fTable  = dynamic_cast<CCfits::Table*>(extMap.find(tablename)->second);
-    if (!fTable)
-    {
-        cerr << "Object '" << tablename << "' returned not a CCfits::Table." << endl;
-        return false;
-    }
-
-    map<string, CCfits::Column*>& tCols = fTable->column();
-    for (map<string, CCfits::Column*>::const_iterator it = tCols.begin(); it != tCols.end(); it++)
-    {
-        fColMap.insert(make_pair(it->first, static_cast<MyColumn*>(it->second)));
-    }
-//    fColMap = fTable->column();
-
-    fTable->makeThisCurrent();
-
-    fTable->getComments();
-    fTable->getHistory();
-    fTable->readAllKeys();
 
     return true;
@@ -426 +226 @@
     }
 
-    cout << "\nFile: " << fFile->name() << "\n";
-
-    const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension();
-    for (std::multimap<string, CCfits::ExtHDU*>::const_iterator it=extMap.begin(); it != extMap.end(); it++)
-    {
-
-        CCfits::Table *table = dynamic_cast<CCfits::Table*>(extMap.find(it->first)->second);
-
-        table->makeThisCurrent();
-        table->readData();
-
-        cout << " " << it->first << " [" << table->rows() << "]\n";
-
-        const map<string, CCfits::Column*> &cols = table->column();
-
-//        for (map<string, CCfits::Column*>::const_iterator id = cols.begin(); ib != cols.end(); ib++)
-//        {
-//            TFORM
-//        }
-
-        for (map<string, CCfits::Column*>::const_iterator ic=cols.begin();
-             ic != cols.end(); ic++)
-        {
-            const MyColumn *col = static_cast<MyColumn*>(ic->second);
-
-            cout << "   " << col->name() << "[" ;
-            if (col->width() == 0)
-                cout << 1;
-            else
-                cout << col->width();
-            cout  << "] * " << col->scale() << " (" << col->unit() << ":" <<  ValueTypeToStr(col->type())<<  ") " << col->comment() << "\n";
-            /*
-             inline size_t Column::repeat () const
-             inline bool   Column::varLength () const
-             inline double Column::zero () const
-             inline const  String& Column::display () const
-             inline const  String& Column::dimen () const
-             inline const  String& Column::TBCOL ()
-             inline const  String& Column::TTYPE ()
-             inline const  String& Column::TFORM ()
-             inline const  String& Column::TDISP ()
-             inline const  String& Column::TZERO ()
-             inline const  String& Column::TDIM  ()
-             inline const  String& Column::TNULL ()
-             inline const  String& Column::TLMIN ()
-             inline const  String& Column::TLMAX ()
-             inline const  String& Column::TDMAX ()
-             inline const  String& Column::TDMIN ()
-            */
-        }
-        cout << '\n';
-    }
+    cout << "\nFile: " << fFilename << "\n";
+
+    cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
+    cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
+
+    for (auto it = fColMap.begin(); it != fColMap.end(); it++)
+    {
+        cout << "   " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
+        for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
+            if (jt->second.value == it->first)
+                cout << jt->second.comment << endl;
+    }
+
+    cout << endl;
     cout << flush;
 }
 
-class MyKeyword : public CCfits::Keyword
-{
-public:
-    CCfits::ValueType keytype() const { return CCfits::Keyword::keytype(); }
-};
-
 void FitsDumper::ListKeywords(ostream &out)
 {
-    map<string,CCfits::Keyword*> keys = fTable->keyWord();
-    for (map<string,CCfits::Keyword*>::const_iterator it=keys.begin();
-         it!=keys.end(); it++)
-    {
-        string str;
-        double d;
-        int l;
-
-        const MyKeyword *kw = static_cast<MyKeyword*>(it->second);
-        kw->keytype();
-        out << "## " << setw(8) << kw->name() << " = " << setw(10);
-        if (kw->keytype()==16)
-            out << ("'"+kw->value(str)+"'");
-        if (kw->keytype()==31)
-            out << kw->value(l);
-        if (kw->keytype()==82)
-            out << kw->value(d);
-        out << " / " << kw->comment() << endl;
+    for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) {
+        out << "## " << setw(8) << it->first << " = " << setw(10);
+        out << "'" << it->second.value << "'" << " / " << it->second.comment << endl;
     }
 }
@@ -512 +253 @@
 void FitsDumper::ListHeader()
 {
-    if (!fTable)
+    if (!fFile)
     {
         cerr << "No table open." << endl;
@@ -518 +259 @@
     }
 
-    cout << "\nTable: " << fTable->name() << " (rows=" << fTable->rows() << ")\n";
-    if (!fTable->comment().empty())
-        cout << "Comment: \t" << fTable->comment() << '\n';
-    if (!fTable->history().empty())
-        cout << "History: \t" << fTable->history() << '\n';
+    cout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
+    if (fKeyMap.find("COMMENT") != fKeyMap.end())
+        cout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
 
     ListKeywords(cout);
     cout << endl;
+
 }
 
@@ -538 +278 @@
         unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
         unsigned long colonIndex = columnNameOnly.find_first_of(':');
-//        cout << bracketIndex0 << " " << bracketIndex1 << " " << colonIndex << endl;
         int columnStart = -1;
         int columnEnd = -1;
@@ -553 +292 @@
                 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
                 columnEnd = columnStart+1;
-//                cout << "Cstart " << columnStart  << " end: " << columnEnd << endl;
             }
             columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
@@ -563 +301 @@
             return false;
         }
-//        cout << "The column name is: " << columnNameOnly << endl;
-        MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(columnNameOnly)->second);
+        fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second;
         if (bracketIndex0 == string::npos)
         {//no range given: use the full range
-            ranges.push_back(make_pair(0, cCol->width()));
+            ranges.push_back(make_pair(0, cCol.num));
             columnStart = 0;
             columnEnd = 1;
@@ -578 +315 @@
                 return false;
             }
-            if (columnEnd>1 && columnEnd > cCol->width())
+            if (columnEnd>1 && columnEnd > (int)(cCol.num))
             {
-                cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol->width() << " vs " << columnEnd << "). Aborting" << endl;
+                cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl;
                 return false;
             }
@@ -607 +344 @@
     }
 
-    // FIXME: Maybe do this when opening a table?
-    const vector<int> offsets = CalculateOffsets();
- //   for (int i=0;i<offsets.size(); i++)
- //       cout << offsets[i] << " ";
-//    cout << endl;
-    if (offsets.size()==0)
-        return false;
-
     ofstream out(filename=="-"?"/dev/stdout":filename);
     if (!out)
@@ -627 +356 @@
     if (filename!="-")
         out << "## File:    \t" << filename << '\n';
-    out << "## Table:   \t" << fTable->name() << '\n';
-    if (!fTable->comment().empty())
-        out << "## Comment: \t" << fTable->comment() << '\n';
-    if (!fTable->history().empty())
-        out << "## History: \t" << fTable->history() << '\n';
-    out << "## NumRows: \t" << fTable->rows() << '\n';
+    out << "## Table:   \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
+        out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
+
+    out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n';
     out << "## --------------------------------------------------------------------------\n";
     ListKeywords(out);
     out << "## --------------------------------------------------------------------------\n";
     out << "#\n";
-    //vector<pair<int, int> >::const_iterator rangesIt = ranges.begin();
-    //the above is soooo yesterday ;) let's try the new auto feature of c++0x
+
     auto rangesIt = ranges.begin();
     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
     {
-        const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
-        if (rangesIt->first != 0 || rangesIt->second != col->width())
+        const fits::Table::Column& col = fColMap[*it];
+//        const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
+        if (rangesIt->first != 0 || rangesIt->second != col.num)
         {
             out << "#";
             for (int i=rangesIt->first; i<rangesIt->second; i++)
-                out << " " << col->name() << "[" << i << "]";
-            out << ": " << col->unit();
+                out << " " << *it << "[" << i << "]";
+            out << ": " << col.unit;
         }
         else
-            out << "# " << col->name() << "[" << col->width() << "]: " << col->unit();
-
-        if (!col->comment().empty())
-            out << " (" <<col->comment() << ")";
+            out << "# " << *it << "[" << col.num << "]: " << col.unit;
+
+//        FIXME: retrive the column comment
+//        if (!col->comment().empty())
+//            out << " (" <<col->comment() << ")";
         out << '\n';
     }
@@ -661 +389 @@
 
     // Loop over all columns in our list of requested columns
-    vector<MyColumn*> columns;
+    vector<pair<char, char*> > columnsData;
+
     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
     {
-        MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(*it)->second);
-        columns.push_back(cCol);
-    }
-    const int size = offsets[offsets.size()-1];
-    unsigned char* fitsBuffer = new unsigned char[size];
-
-    int status = 0;
-    int endIndex = (fColMap.find("Time") == fColMap.end()) ? fTable->rows()-1 : fTable->rows();
-    for (int i=1; i<=endIndex; i++)
-    {
-        fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
-        if (status)
+        fits::Table::Column& cCol = fColMap.find(*it)->second;
+        columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
+        fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
+    }
+    int numRows = fFile->GetInt("NAXIS2");
+    int row = 0;
+    while (fFile->GetNextRow() && row < numRows)
+    {
+        row++;
+        rangesIt = ranges.begin();
+        for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
         {
-            char text[30];//max length of cfitsio error strings (from doc)
-            fits_get_errstatus(status, text);
-
-            cerr << "Reading row " << i << " failed: " << text << " (fits_insert_rows,rc=" << status << ")" << endl;
-            break;
-        }
-        if (WriteRow(out, columns, offsets, fitsBuffer, ranges)<0)
-        {
-            status=1;
-            break;
-        }
-    }
-    delete[] fitsBuffer;
-
-    return status==0;
-}
-
-// --------------------------------------------------------------------------
-//
-//! Perform the actual dump, based on the current parameters
-//
-/*
-bool FitsDumper::Plot(const vector<string> &list)
-{
-   for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
-        if (fColMap.find(*it) == fColMap.end())
-        {
-            cerr << "WARNING - Column '" << *it << "' not found in table." << endl;
-            return false;
-        }
-
-    // FIXME: Maybe do this when opening a table?
-    const vector<int> offsets = CalculateOffsets();
-    if (offsets.size()==0)
-        return false;
-
-    // Loop over all columns in our list of requested columns
-    const CCfits::Column *col[3] =
-    {
-        list.size()>0 ? fColMap.find(list[0])->second : 0,
-        list.size()>1 ? fColMap.find(list[1])->second : 0,
-        list.size()>2 ? fColMap.find(list[2])->second : 0
-    };
-
-    const int size = offsets[offsets.size()-1];
-    unsigned char* fitsBuffer = new unsigned char[size];
-
-    const int idx = 0;
-
-    // CCfits starts counting at 1 not 0
-    const size_t pos[3] =
-    {
-        col[0] ? offsets[col[0]->index()-1] + idx*col[0]->width()*ValueTypeToSize(col[0]->type()) : 0,
-        col[1] ? offsets[col[1]->index()-1] + idx*col[1]->width()*ValueTypeToSize(col[1]->type()) : 0,
-        col[2] ? offsets[col[2]->index()-1] + idx*col[2]->width()*ValueTypeToSize(col[2]->type()) : 0
-    };
-
-    int status = 0;
-    for (int i=1; i<=fTable->rows(); i++)
-    {
-        fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
-        if (status)
-        {
-            cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl;
-            break;
-        }
-
-        try
-        {
-            const double x[3] =
+            for (int i=rangesIt->first; i<rangesIt->second; i++)
             {
-                col[0] ? PtrToDouble(fitsBuffer+pos[0], col[0]->type()) : 0,
-                col[1] ? PtrToDouble(fitsBuffer+pos[1], col[1]->type()) : 0,
-                col[2] ? PtrToDouble(fitsBuffer+pos[2], col[2]->type()) : 0
-            };
-        }
-        catch (const runtime_error &e)
-        {
-            cerr << e.what() << endl;
-            status=1;
-            break;
-        }
-    }
-    delete[] fitsBuffer;
-
-    return status==0;
-}
-*/
+            switch (it->first) {
+                case 'L':
+                        out << reinterpret_cast<bool*>(it->second)[i] << " ";
+                        break;
+                case 'B':
+                        out << reinterpret_cast<char*>(it->second)[i] << " ";
+                        break;
+                case 'I':
+                        out << reinterpret_cast<int16_t*>(it->second)[i] << " ";
+                        break;
+                case 'J':
+                        out << reinterpret_cast<int32_t*>(it->second)[i] << " ";
+                        break;
+                case 'K':
+                        out << reinterpret_cast<int64_t*>(it->second)[i] << " ";
+                        break;
+                case 'E':
+                        out << reinterpret_cast<float*>(it->second)[i] << " ";
+                        break;
+                case 'D':
+                        out << reinterpret_cast<double*>(it->second)[i] << " ";
+                        break;
+                default:
+                    ;
+            }
+            }
+        }
+    }
+    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
+        delete[] it->second;
+    return true;
+}
 
 // --------------------------------------------------------------------------
@@ -780 +454 @@
     }
 
+    if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
+    {
+        cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl;
+        return -1;
+    }
+
     if (conf.Get<bool>("list"))
         List();
@@ -789 +469 @@
     }
 
-#ifdef PLOTTING_PLEASE
-    if (conf.Get<bool>("graph"))
+    if (conf.Get<bool>("stat"))
     {
         if (!conf.Has("col"))
@@ -797 +476 @@
             return 0;
         }
+        doStatsPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
+        return 0;
+    }
+
+#ifdef PLOTTING_PLEASE
+    if (conf.Get<bool>("graph"))
+    {
+        if (!conf.Has("col"))
+        {
+            cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
+            return 0;
+        }
+        if (conf.Get<bool>("dots"))
+            fDotsPlease = true;
         doCurvesDisplay(conf.Get<vector<string>>("col"),
                         conf.Get<string>("tablename"));
@@ -839 +532 @@
 {
     // 
+}
+
+template<typename T>
+void displayStats(char* array, int numElems)
+{
+    sort(reinterpret_cast<T*>(array), &reinterpret_cast<T*>(array)[numElems]);
+    cout << "Min: " << reinterpret_cast<T*>(array)[0] << endl;
+    cout << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
+    if (numElems%2 == 0)
+        cout << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl;
+    else
+        cout << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;
+
+    double average = 0;
+    for (int i=0;i<numElems;i++)
+        average += reinterpret_cast<T*>(array)[i];
+    cout << "Mea: " << average/(double)numElems << endl;
+
+}
+int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
+{
+
+    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
+    vector<pair<int, int> > ranges;
+    vector<string> listNamesOnly;
+
+    if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
+    {
+        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
+        return false;
+    }
+
+    ofstream out(filename=="-"?"/dev/stdout":filename);
+    if (!out)
+    {
+        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
+        return false;
+    }
+
+    out.precision(precision);
+
+    // Loop over all columns in our list of requested columns
+    vector<pair<char, char*> > columnsData;
+    vector<char*> statData;
+    int numRows = fFile->GetInt("NAXIS2");
+    auto rangesIt = ranges.begin();
+    for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
+    {
+        fits::Table::Column& cCol = fColMap.find(*it)->second;
+        columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
+        statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]);
+        fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
+    }
+
+    int row = 0;
+    while (fFile->GetNextRow() && row < numRows)
+    {
+        rangesIt = ranges.begin();
+        auto statsIt = statData.begin();
+        for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
+        {
+            int span = rangesIt->second - rangesIt->first;
+            for (int i=rangesIt->first; i<rangesIt->second; i++)
+            {
+            switch (it->first) {
+                case 'L':
+                        reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<bool*>(it->second)[i];
+                        break;
+                case 'B':
+                        reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<char*>(it->second)[i];
+                        break;
+                case 'I':
+                        reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int16_t*>(it->second)[i];
+                        break;
+                case 'J':
+                        reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int32_t*>(it->second)[i];
+                        break;
+                case 'K':
+                        reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int64_t*>(it->second)[i];
+                        break;
+                case 'E':
+                        reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<float*>(it->second)[i];
+                        break;
+                case 'D':
+                        reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<double*>(it->second)[i];
+                        break;
+                default:
+                    ;
+            }
+            }
+        }
+        row++;
+    }
+    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
+        delete[] it->second;
+
+    //okay. So now I've got ALL the data, loaded.
+    //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
+    rangesIt = ranges.begin();
+    auto statsIt = statData.begin();
+    int round;
+    auto nameIt = listNamesOnly.begin();
+    for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
+    {
+        int span = rangesIt->second - rangesIt->first;
+        cout << *nameIt << ": " << endl;
+        switch (it->first) {
+            case 'L':
+                    displayStats<bool>(*statsIt, numRows*span);
+                    break;
+            case 'B':
+                    displayStats<char>(*statsIt, numRows*span);
+                    break;
+            case 'I':
+                    displayStats<int16_t>(*statsIt, numRows*span);
+                    break;
+            case 'J':
+                    displayStats<int32_t>(*statsIt, numRows*span);
+                    break;
+            case 'K':
+                    displayStats<int64_t>(*statsIt, numRows*span);
+                    break;
+            case 'E':
+                    displayStats<float>(*statsIt, numRows*span);
+                    break;
+            case 'D':
+                    displayStats<double>(*statsIt, numRows*span);
+                    break;
+            default:
+                ;
+        }
+    }
+    return true;
 }
 #ifdef PLOTTING_PLEASE
@@ -894 +720 @@
     for (unsigned int i=0;i<list.size();i++)
         totalSize += ranges[i].second - ranges[i].first;
- //   cout << "Total size: " << totalSize << endl;
+
     vector<QwtPlotCurve*> curves(totalSize);
     int ii=0;
@@ -924 +750 @@
         };
         ii++;
-        (*it)->setStyle(QwtPlotCurve::Lines);
+        if (fDotsPlease)
+            (*it)->setStyle(QwtPlotCurve::Dots);
+        else
+            (*it)->setStyle(QwtPlotCurve::Lines);
         (*it)->attach(plot);
     }
     plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
 
-     const vector<int> offsets = CalculateOffsets();
-     if (offsets.size()==0)
-         return false;
-
-
-     // Loop over all columns in our list of requested columns
-      vector<MyColumn*> columns;
+
+      vector<pair<char, char*> > columnsData;
+
       for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
       {
-          MyColumn *cCol = static_cast<MyColumn*>(fColMap.find(*it)->second);
-          columns.push_back(cCol);
+          fits::Table::Column& cCol = fColMap.find(*it)->second;
+          columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
+          fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
       }
+
       //add the time column to the given columns
       if (fColMap.find("Time") == fColMap.end() && fColMap.find("UnixTimeUTC") == fColMap.end())
@@ -947 +774 @@
           return false;
       }
-      MyColumn* timeCol;
-      bool unixTime = false;
-      int endIndex = 0;
-      if (fColMap.find("Time") != fColMap.end())
-      {
-          timeCol = static_cast<MyColumn*>(fColMap.find("Time")->second);
+      const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
+      bool unixTime = (fColMap.find("Time") == fColMap.end());
+      if (unixTime)
+          ranges.push_back(make_pair(0,2));
+      else
           ranges.push_back(make_pair(0,1));
-          endIndex = fTable->rows();
-      }
-      if (fColMap.find("UnixTimeUTC") != fColMap.end())
-      {
-          timeCol = static_cast<MyColumn*>(fColMap.find("UnixTimeUTC")->second);
-          ranges.push_back(make_pair(0,2));
-          endIndex = fTable->rows()-1;
-          unixTime = true;
-      }
-      columns.push_back(timeCol);
-      /////
-      const int size = offsets[offsets.size()-1];
-      unsigned char* fitsBuffer = new unsigned char[size];
+      columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
+      fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);
 
 //      stringstream str;
@@ -977 +792 @@
       str.precision(20);
       for (auto it=xValues.begin(); it!=xValues.end(); it++)
-          *it = new double[fTable->rows()];
-
-      yValues = new double[fTable->rows()];
+          *it = new double[fFile->GetInt("NAXIS2")];
+
+      yValues = new double[fFile->GetInt("NAXIS2")];
 
       cout.precision(3);
-      for (int i=1; i<=endIndex; i++)
+      int endIndex = 0;
+      int numRows = fFile->GetInt("NAXIS2");
+      for (int i=1;i<numRows;i++)
       {
-          cout << "/r" << "Constructing graph " << ((float)(i)/(float)(endIndex))*100.0 << "%";
-          fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status);
-          if (status)
+          fFile->GetNextRow();
+          cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
+          endIndex++;
+          auto rangesIt = ranges.begin();
+          for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
           {
-              cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl;
-              break;
+              for (int i=rangesIt->first; i<rangesIt->second; i++)
+              {
+              switch (it->first) {
+                  case 'L':
+                          str << reinterpret_cast<bool*>(it->second)[i] << " ";
+                          break;
+                  case 'B':
+                          str << reinterpret_cast<char*>(it->second)[i] << " ";
+                          break;
+                  case 'I':
+                          str << reinterpret_cast<int16_t*>(it->second)[i] << " ";
+                          break;
+                  case 'J':
+                          str << reinterpret_cast<int32_t*>(it->second)[i] << " ";
+                          break;
+                  case 'K':
+                          str << reinterpret_cast<int64_t*>(it->second)[i] << " ";
+                          break;
+                  case 'E':
+                          str << reinterpret_cast<float*>(it->second)[i] << " ";
+                          break;
+                  case 'D':
+                          str << reinterpret_cast<double*>(it->second)[i] << " ";
+                          break;
+                  default:
+                      ;
+              }
+              }
           }
-          if (WriteRow(str, columns, offsets, fitsBuffer, ranges)<0)
-          {
-              status=1;
-              cerr << "An Error occured while reading the fits row " << i << endl;
-              return -1;
-          }
-//          yValues[i-1] = i;
+
           for (auto it=xValues.begin(); it!= xValues.end(); it++)
           {
               str >> (*it)[i-1];
-//              cout << (*it)[i-1] << " ";
           }
           if (unixTime)
@@ -1007 +845 @@
               long u1, u2;
               str >> u1 >> u2;
- //             cout << u1 << " " << u2;
+
               boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
                       boost::posix_time::seconds(u1) +  boost::posix_time::microsec(u2));
@@ -1013 +851 @@
               Time mjdTime(unixTimeT);
               yValues[i-1] = mjdTime.Mjd();
- //             cout << " " << mjdTime.Mjd() << endl;
           }
           else
@@ -1020 +857 @@
               if (yValues[i-1] < 40587)
                   yValues[i-1] += 40587;
+
+              Time t(yValues[i-1]);
+              string time = t.GetAsStr("%H:%M:%S%F");
+              while (time[time.size()-1] == '0' && time.size() > 2)
+              {
+                  time = time.substr(0, time.size()-1);
+              }
           }
           if (i==1)
@@ -1027 +871 @@
               plot->setTitle(title.c_str());
           }
-//          cout << yValues[i-1] << " ";
-//          cout << endl;
       }
       //set the actual data.
       auto jt = xValues.begin();
       for (auto it=curves.begin(); it != curves.end(); it++, jt++)
-          (*it)->setRawData(yValues, *jt, endIndex);
+          (*it)->setRawData(yValues, *jt, endIndex-1);
 
       QStack<QRectF> stack;
@@ -1061 +903 @@
       zoom.setZoomStack(stack);
 
-      delete[] fitsBuffer;
+//      delete[] fitsBuffer;
+      for (auto it = columnsData.begin(); it != columnsData.end(); it++)
+          delete[] it->second;
     window.resize(600, 400);
     window.show();
@@ -1101 +945 @@
         ("list,l",      po_switch(),                "List all tables and columns in file")
         ("header,h",    po_switch(),                "Dump header of given table")
+        ("stat,s",      po_switch(),                "Perform statistics instead of dump")
 #ifdef PLOTTING_PLEASE
         ("graph,g",     po_switch(),                "Plot the columns instead of dumping them")
+        ("dots,d",      po_switch(),                "Plot using dots instead of lines")
 #endif
         ;