#ifndef MARS_fits #define MARS_fits #include #include #include #include #include #include #include #ifdef __CINT__ #define off_t size_t #endif #if !defined(__MARS__) && !defined(__CINT__) #include #endif #ifndef __MARS__ #include #include #include #define gLog cerr #define ___err___ "" #define ___warn___ "" #define ___all___ "" #else #include "MLog.h" #include "MLogManip.h" #define ___err___ err #define ___warn___ warn #define ___all___ all #endif #if defined(HAVE_ZLIB) || defined(__CINT__) #include "izstream.h" #else #include #define izstream ifstream #warning Support for zipped FITS files disabled. #endif #include "checksum.h" #ifndef __MARS__ namespace std { #else using namespace std; #endif 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; string fitsString; template T Get() const { T t; istringstream str(value); str >> t; return t; } }; struct Table { off_t offset; bool is_compressed; string name; size_t bytes_per_row; size_t num_rows; size_t num_cols; size_t total_bytes; // NAXIS1*NAXIS2 struct Column { size_t offset; size_t num; size_t size; size_t bytes; // num*size char type; string unit; Compression_t comp; }; typedef map Keys; typedef map Columns; typedef vector SortedColumns; Columns cols; SortedColumns sorted_cols; Keys keys; int64_t datasum; string Trim(const string &str, char c=' ') const { // Trim Both leading and trailing spaces const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces const size_t pend = 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==pstart || string::npos==pend) return string(); return str.substr(pstart, pend-pstart+1); } bool Check(const string &key, char type, const string &value="") const { const Keys::const_iterator it = keys.find(key); if (it==keys.end()) { ostringstream str; str << "Key '" << key << "' not found."; #ifdef __EXCEPTIONS throw runtime_error(str.str()); #else gLog << ___err___ << "ERROR - " << str.str() << endl; return false; #endif } if (it->second.type!=type) { 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; return false; #endif } if (!value.empty() && it->second.value!=value) { 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; return false; #endif } return true; } Keys ParseBlock(const vector &vec) const { map rc; for (unsigned int i=0; i &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("GCOUNT", 'I', "1") || !Check("EXTNAME", 'T') || !Check("NAXIS1", 'I') || !Check("NAXIS2", 'I') || !Check("TFIELDS", 'I')) return; 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("NAXIS1")*Get("NAXIS2"); bytes_per_row = is_compressed ? Get("ZNAXIS1") : Get("NAXIS1"); num_rows = is_compressed ? Get("ZNAXIS2") : Get("NAXIS2"); num_cols = Get("TFIELDS"); datasum = is_compressed ? Get("ZDATASUM", -1) : Get("DATASUM", -1); size_t bytes = 0; const string tFormName = is_compressed ? "ZFORM" : "TFORM"; for (long long unsigned int i=1; i<=num_cols; i++) { const string num(to_string(i)); if (!Check("TTYPE"+num, 'T') || !Check(tFormName+num, 'T')) return; const string id = Get("TTYPE"+num); const string fmt = Get(tFormName+num); const string unit = Get("TUNIT"+num, ""); const string comp = Get("ZCTYP"+num, ""); Compression_t compress = kCompUnknown; if (comp == "FACT") compress = kCompFACT; istringstream sin(fmt); int n = 0; sin >> n; if (!sin) 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 'L': // logical case 'A': // char 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 // 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 // case 'P': size = 8; break; // array descriptor (32bit) // case 'Q': size = 16; break; // array descriptor (64bit) default: { ostringstream str; str << "FITS format TFORM='" << fmt << "' not yet supported."; #ifdef __EXCEPTIONS throw runtime_error(str.str()); #else gLog << ___err___ << "ERROR - " << str.str() << endl; return; #endif } } const Table::Column col = { bytes, n, size, n*size, type, unit, compress}; cols[id] = col; sorted_cols.push_back(col); bytes += n*size; } if (bytes!=bytes_per_row) { #ifdef __EXCEPTIONS throw runtime_error("Column size mismatch"); #else gLog << ___err___ << "ERROR - Column size mismatch" << endl; return; #endif } name = Get("EXTNAME"); } void PrintKeys(bool display_all=false) const { for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++) { if (!display_all && (it->first.substr(0, 6)=="TTYPE" || it->first.substr(0, 6)=="TFORM" || it->first.substr(0, 6)=="TUNIT" || it->first=="TFIELDS" || it->first=="XTENSION" || it->first=="NAXIS" || it->first=="BITPIX" || it->first=="PCOUNT" || it->first=="GCOUNT") ) continue; gLog << ___all___ << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl; }} void PrintColumns() const { typedef map, 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 << "| "; 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; case 'I': gLog << "short(16)"; break; case 'J': gLog << "int(32)"; break; case 'K': gLog << "int(64)"; break; case 'E': gLog << "float(32)"; break; case 'D': gLog << "double(64)"; break; } gLog << ": " << it->first.second << " [" << it->second.unit << "]" << endl; } } operator bool() const { return !name.empty(); } bool HasKey(const string &key) const { return keys.find(key)!=keys.end(); } bool HasColumn(const 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 T Get(const string &key) const { const map::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; return T(); #endif } return it->second.Get(); } // Values of keys are always signed template T Get(const string &key, const T &deflt) const { const map::const_iterator it = keys.find(key); return it==keys.end() ? deflt :it->second.Get(); } size_t GetN(const 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("THEAP")) return 0; const size_t shift = Get("THEAP"); 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("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)); } }; protected: ofstream fCopy; Table fTable; typedef pair Address; typedef vector
Addresses; //map fAddresses; Addresses fAddresses; #if defined(__MARS__) || defined(__CINT__) typedef map Pointers; #else typedef unordered_map Pointers; #endif Pointers fPointers; vector> fGarbage; vector fBufferRow; vector fBufferDat; size_t fRow; Checksum fChkHeader; Checksum fChkData; bool ReadBlock(vector &vec) { int endtag = 0; for (int i=0; i<36; i++) { char c[81]; c[80] = 0; read(c, 80); if (!good()) break; fChkHeader.add(c, 80); // if (c[0]==0) // return vector(); string str(c); // if (!str.empty()) // cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl; 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); } // Make sure that no empty vector is returned if (endtag && vec.size()%36==0) vec.emplace_back("END = '' / "); return endtag==2; } string Compile(const string &key, int16_t i=-1) const { if (i<0) return key; ostringstream str; str << key << i; return str.str(); } void Constructor(const string &fname, string fout, const string& tableName, bool force) { char simple[10]; read(simple, 10); if (!good()) return; 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; return; #endif } seekg(0); while (good()) { vector 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 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; return; #endif } if (block.size()%36) { if (!rc && !force) { clear(rdstate()|ios::badbit); #ifdef __EXCEPTIONS throw runtime_error("END keyword missing in FITS header."); #else gLog << ___err___ << "ERROR - END keyword missing in FITS file... file might be corrupted." << endl; return; #endif } break; } } if (block.empty()) break; if (block[0].substr(0, 9)=="SIMPLE =") { fChkHeader.reset(); continue; } if (block[0].substr(0, 9)=="XTENSION=") { fTable = Table(block, tellg()); fRow = (size_t)-1; if (!fTable) { clear(rdstate()|ios::badbit); return; } // Check for table name. Skip until eof or requested table are found. // skip the current table? if ((!tableName.empty() && tableName!=fTable.Get("EXTNAME")) || ( tableName.empty() && "ZDrsCellOffsets"==fTable.Get("EXTNAME"))) { const streamoff skip = fTable.GetTotalBytes(); seekg(skip, ios_base::cur); fChkHeader.reset(); continue; } fBufferRow.resize(fTable.bytes_per_row + 8-fTable.bytes_per_row%4); fBufferDat.resize(fTable.bytes_per_row); break; } } 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()|ios::badbit); #ifdef __EXCEPTIONS throw runtime_error("Could not open output file."); #else gLog << ___err___ << "ERROR - Failed to open output file." << endl; #endif } const streampos p = tellg(); seekg(0); vector buf(p); read(buf.data(), p); fCopy.write(buf.data(), p); if (!fCopy) clear(rdstate()|ios::badbit); } public: fits(const string &fname, const string& tableName="", bool force=false) : izstream(fname.c_str()) { Constructor(fname, "", tableName, force); if ((fTable.is_compressed && !force) || (fTable.name == "ZDrsCellOffsets" && !force)) { #ifdef __EXCEPTIONS throw 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." << endl; #endif clear(rdstate()|ios::badbit); } } fits(const string &fname, const string &fout, const string& tableName, bool force=false) : izstream(fname.c_str()) { Constructor(fname, fout, tableName, force); if ((fTable.is_compressed && !force) || (fTable.name == "ZDrsCellOffsets" && !force)) { #ifdef __EXCEPTIONS throw 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." << endl; #endif clear(rdstate()|ios::badbit); } } fits() : izstream() { } ~fits() { copy(istreambuf_iterator(*this), istreambuf_iterator(), ostreambuf_iterator(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()|ios::badbit); } else if (fCopy.is_open()) clear(rdstate()|ios::badbit); } void ZeroBufferForChecksum(vector& 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 void revcpy(char *dest, const char *src, const int &num) { const char *pend = src + num*N; for (const char *ptr = src; ptr(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__) virtual bool GetRow(size_t row, bool check=true) #else virtual bool GetRowNum(size_t row, bool check=true) #endif { if (check && row>=fTable.num_rows) return false; const uint8_t offset = ReadRow(row); if (!good()) return good(); const char *ptr = fBufferRow.data() + offset; for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++) { const Table::Column &c = it->second; const char *src = ptr + c.offset; char *dest = reinterpret_cast(it->first); MoveColumnDataToUserSpace(dest, src, c); } return good(); } bool GetNextRow(bool check=true) { #if !defined(__MARS__) && !defined(__CINT__) return GetRow(fRow+1, check); #else return GetRowNum(fRow+1, check); #endif } virtual bool SkipNextRow() { seekg(fTable.offset+(++fRow)*fTable.bytes_per_row); return good(); } static bool Compare(const Address &p1, const Address &p2) { return p1.first>p2.first; } template const T &GetAs(const string &name) { return *reinterpret_cast(fPointers[name]); } void *SetPtrAddress(const string &name) { 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; 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 bool SetPtrAddress(const 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; return false; #endif } if (sizeof(T)!=fTable.cols[name].size) { 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; return false; #endif } if (cnt!=fTable.cols[name].num) { 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; return false; #endif } // if (fAddresses.count(ptr)>0) // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; //fAddresses[ptr] = fTable.cols[name]; fPointers[name] = ptr; fAddresses.emplace_back(ptr, fTable.cols[name]); sort(fAddresses.begin(), fAddresses.end(), Compare); return true; } template bool SetRefAddress(const string &name, T &ptr) { return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T)); } template bool SetVecAddress(const string &name, vector &vec) { return SetPtrAddress(name, vec.data(), vec.size()); } template T Get(const string &key) const { return fTable.Get(key); } template T Get(const string &key, const string &deflt) const { return fTable.Get(key, deflt); } bool SetPtrAddress(const 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; return false; #endif } if (cnt && cnt!=fTable.cols[name].num) { 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; return false; #endif } // if (fAddresses.count(ptr)>0) // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; //fAddresses[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); } 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.sorted_cols;} const Table::Keys &GetKeys() const { return fTable.GetKeys();} int64_t GetInt(const string &key) const { return fTable.Get(key); } uint64_t GetUInt(const string &key) const { return fTable.Get(key); } double GetFloat(const string &key) const { return fTable.Get(key); } string GetStr(const string &key) const { return fTable.Get(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==(size_t)-1 ? 0 : fRow; } operator bool() const { return fTable && fTable.offset!=0; } void PrintKeys() const { fTable.PrintKeys(); } void PrintColumns() const { fTable.PrintColumns(); } 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("NAXIS2"); } virtual size_t GetBytesPerRow() const { return fTable.Get("NAXIS1"); } }; #ifndef __MARS__ }; #endif #endif