Changeset 16284 for trunk

Timestamp:

05/27/13 01:26:35 (12 years ago)

Author:

lyard

Message:

added fitsCompressor and compressed files read headers

File:

• trunk/Mars/mcore/fits.h (modified) (23 diffs)

trunk/Mars/mcore/fits.h

@@ -69 +69 @@
 {
 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.
+    typedef enum
+    {
+        UNCOMPRESSED,
+        SMOOTHMAN
+    } FitsCompression;
+
     struct Entry
     {
@@ -74 +82 @@
         string value;
         string comment;
+        string fitsString;
 
         template<typename T>
@@ -90 +99 @@
     {
         off_t offset;
+
+        bool isCompressed;
 
         string name;
@@ -103 +114 @@
             char   type;
             string unit;
+            FitsCompression comp;
         };
 
         typedef map<string, Entry>  Keys;
         typedef map<string, Column> Columns;
+        typedef vector<Column> SortedColumns;
 
         Columns cols;
+        SortedColumns sortedCols;
         Keys    keys;
 
@@ -200 +214 @@
 
                     // Set value, comment and type
-                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
+                    if (ppp==string::npos)
+                        com = "";
+                    else
+                    {// comments could be just spaces. take care of this.
+                        if (val.size() == ppp+1)
+                            com = "";
+                        else
+                            com = Trim(val.substr(ppp+1));
+                    }
+//                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
                     val  = Trim(val.substr(1, p-2));
                     type = 'T';
@@ -208 +231 @@
                     const size_t p = val.find_first_of('/');
 
-                    com = Trim(val.substr(p+2));
+                    if (val.size() == p+1)
+                        com = "";
+                    else
+                        com = Trim(val.substr(p+2));
                     val = Trim(val.substr(0, p));
 
@@ -217 +243 @@
                 }
 
-                const Entry e = { type, val, com };
+                const Entry e = { type, val, com, vec[i] };
                 rc[key] = e;
             }
@@ -226 +252 @@
         Table() : offset(0) { }
         Table(const vector<string> &vec, off_t off) :
-            offset(off), keys(ParseBlock(vec))
-        {
+            offset(off), isCompressed(false), keys(ParseBlock(vec))
+        {
+            if (HasKey("ZTABLE"))
+            if (Check("ZTABLE", 'B', "T"))
+            {
+                isCompressed = true;
+            }
+
             if (!Check("XTENSION", 'T', "BINTABLE") ||
                 !Check("NAXIS",    'I', "2")        ||
                 !Check("BITPIX",   'I', "8")        ||
-                !Check("PCOUNT",   'I', "0")        ||
                 !Check("GCOUNT",   'I', "1")        ||
                 !Check("EXTNAME",  'T')             ||
@@ -239 +270 @@
                 return;
 
-            bytes_per_row = Get<size_t>("NAXIS1");
-            num_rows      = Get<size_t>("NAXIS2");
+            if (isCompressed)
+            {
+                if (!Check("ZPCOUNT", 'I', "0"))
+                    return;
+            }
+            else
+            {
+                if (!Check("PCOUNT", 'I', "0"))
+                    return;
+            }
+
+            bytes_per_row = isCompressed ? Get<size_t>("ZNAXIS1") : Get<size_t>("NAXIS1");
+            num_rows      = isCompressed ? Get<size_t>("ZNAXIS2") : Get<size_t>("NAXIS2");
             num_cols      = Get<size_t>("TFIELDS");
-            datasum       = Get<int64_t>("DATASUM", -1);
+            datasum       = isCompressed ? Get<int64_t>("ZDATASUM", -1) : Get<int64_t>("DATASUM", -1);
 
             size_t bytes = 0;
+            string tFormName = isCompressed ? "ZFORM" : "TFORM";
             for (size_t i=1; i<=num_cols; i++)
             {
@@ -251 +294 @@
 
                 if (!Check("TTYPE"+num.str(), 'T') ||
-                    !Check("TFORM"+num.str(), 'T'))
+                    !Check(tFormName+num.str(), 'T'))
                     return;
 
                 const string id   = Get<string>("TTYPE"+num.str());
-                const string fmt  = Get<string>("TFORM"+num.str());
+                const string fmt  = Get<string>(tFormName+num.str());
                 const string unit = Get<string>("TUNIT"+num.str(), "");
+                const string comp = Get<string>("ZCTYP"+num.str(), "");
+
+                FitsCompression compress = UNCOMPRESSED;
+                if (comp == "SMOO")
+                    compress = SMOOTHMAN;
 
                 istringstream sin(fmt);
@@ -297 +345 @@
                 }
 
-                const Table::Column col = { bytes, n, size, type, unit };
+                const Table::Column col = { bytes, n, size, type, unit, compress};
 
                 cols[id] = col;
+                sortedCols.push_back(col);
                 bytes  += n*size;
             }
@@ -420 +469 @@
     };
 
-private:
+protected:
     ofstream fCopy;
 
@@ -500 +549 @@
     }
 
-    void Constructor(const string &fname, string fout, bool force)
+    //There may be a gap between the main table and the start of the heap: this computes the offset
+    size_t ComputeGapFromDataToHeap(const Table& table)
+    {
+        //get row width
+        const size_t rowWidth = table.Get<size_t>("NAXIS1");
+        const size_t numRows  = table.Get<size_t>("NAXIS2");
+
+        //get offset of special data area from start of main table
+        size_t heapShift = HasKey("THEAP") ? table.Get<size_t>("THEAP") : 0;
+        if (heapShift > rowWidth*numRows)
+            heapShift -= rowWidth*numRows;
+        else
+            heapShift = 0;
+
+        return heapShift;
+    }
+
+    //skip the current table
+    bool SkipCurrentTable(const vector<string>& vector)
+    {
+        Table currentTable(vector, tellg());
+        //get row width
+        const size_t rowWidth = currentTable.Get<size_t>("NAXIS1");
+        const size_t numRows  = currentTable.Get<size_t>("NAXIS2");
+
+        //get offset of special data area from start of main table
+        const off_t heapShift = ComputeGapFromDataToHeap(currentTable);
+
+        //and special data area size
+        const off_t heapSize  = HasKey("PCOUNT") ? currentTable.Get<size_t>("PCOUNT") : 0;
+
+        //skip the row of the table, padded to 2880 bytes
+        off_t paddingSize     = (rowWidth*numRows + heapSize + heapShift)%2880;
+        paddingSize += (paddingSize==0) ? 0 : 2880 - paddingSize;
+
+        const off_t whereDoIGo = tellg() + (off_t)(rowWidth*numRows) + heapSize + heapShift + paddingSize;
+        seekg(whereDoIGo);
+
+        return good();
+    }
+
+    void Constructor(const string &fname, string fout, const string& tableName, bool force, bool mightFail)
     {
         char simple[10];
@@ -529 +619 @@
                 if (!good())
                 {
+                    //we allow that the table is not found (to be checked later on with bool casting)
+                    if (mightFail) { return;}
+
                     clear(rdstate()|ios::badbit);
 #ifdef __EXCEPTIONS
@@ -565 +658 @@
             if (block[0].substr(0, 9)=="XTENSION=")
             {
-                // FIXME: Check for table name
-
+                //Check for table name. Skip until eof or requested table are found.
+                if (tableName != "")
+                {
+                    Table currentTable(block, tellg());
+                    string currentName = currentTable.Get<string>("EXTNAME");
+                    if (currentName != tableName)
+                    {
+                        SkipCurrentTable(block);
+                        fChkHeader.reset();
+                        continue;
+                    }
+                }
                 fTable = Table(block, tellg());
                 fRow   = (size_t)-1;
@@ -630 +733 @@
 
 public:
-    fits(const string &fname, bool force=false) : izstream(fname.c_str())
-    {
-        Constructor(fname, "", force);
-    }
-
-    fits(const string &fname, const string &fout, bool force=false) : izstream(fname.c_str())
-    {
-        Constructor(fname, fout, force);
+    fits(const string &fname, const string& tableName="", bool force=false, bool mightFail=false) : izstream(fname.c_str())
+    {
+        Constructor(fname, "", tableName, force, mightFail);
+    }
+
+    fits(const string &fname, const string &fout, const string& tableName="", bool force=false, bool mightFail=false) : izstream(fname.c_str())
+    {
+        Constructor(fname, fout, tableName, force, mightFail);
     }
 
@@ -647 +750 @@
     }
 
-    uint8_t ReadRow(size_t row)
+    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);
+    }
+
+    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;
@@ -671 +779 @@
         *--ie = 0;
 
-        read(fBufferRow.data()+offset, fTable.bytes_per_row);
-        //fin.clear(fin.rdstate()&~ios::eofbit);
+        StageRow(row, fBufferRow.data()+offset);
 
         if (row==fRow+1)
@@ -699 +806 @@
     }
 
+    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.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;
+        }
+    }
+
 #if !defined(__MARS__) && !defined(__CINT__)
-    bool GetRow(size_t row, bool check=true)
-#else
-    bool GetRowNum(size_t row, bool check=true)
+    virtual bool GetRow(size_t row, bool check=true)
+#else
+    virtual bool GetRowNum(size_t row, bool check=true)
 #endif
     {
@@ -721 +841 @@
             char *dest = reinterpret_cast<char*>(it->first);
 
-            // 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.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);
         }
 
@@ -744 +856 @@
     }
 
-    bool SkipNextRow()
+    virtual bool SkipNextRow()
     {
         seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
@@ -891 +1003 @@
     bool     HasColumn(const string& col) const { return fTable.HasColumn(col);}
     const Table::Columns &GetColumns() const { return fTable.GetColumns();}
+    const Table::SortedColumns& GetSortedColumns() const { return fTable.sortedCols;}
     const Table::Keys &GetKeys() const { return fTable.GetKeys();}
 
@@ -914 +1027 @@
     bool IsFileOk() const { return (fChkHeader+fChkData).valid(); }
 
+    bool IsCompressedFITS() const { return fTable.isCompressed;}
 };