#ifndef MARS_fits #define MARS_fits #ifdef __CINT__ #define int8_t Char_t #define int16_t Short_t #define int32_t Int_t #define int64_t Long64_t #define uint8_t UChar_t #define uint16_t UShort_t #define uint32_t UInt_t #define uint64_t ULong64_t #else #include #endif #include #include #include #include #include #ifdef __EXCEPTIONS #include #endif #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. typedef enum { UNKNOWN, FACT } FitsCompression; 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 isCompressed; 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; FitsCompression comp; }; typedef map Keys; typedef map Columns; typedef vector SortedColumns; Columns cols; SortedColumns sortedCols; 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), isCompressed(false), keys(ParseBlock(vec)) { if (HasKey("ZTABLE") && Check("ZTABLE", 'B', "T")) isCompressed = true; 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 (isCompressed) { 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 = isCompressed ? Get("ZNAXIS1") : Get("NAXIS1"); num_rows = isCompressed ? Get("ZNAXIS2") : Get("NAXIS2"); num_cols = Get("TFIELDS"); datasum = isCompressed ? Get("ZDATASUM", -1) : Get("DATASUM", -1); size_t bytes = 0; string tFormName = isCompressed ? "ZFORM" : "TFORM"; for (size_t i=1; i<=num_cols; i++) { ostringstream num; num << i; if (!Check("TTYPE"+num.str(), 'T') || !Check(tFormName+num.str(), 'T')) return; const string id = Get("TTYPE"+num.str()); const string fmt = Get(tFormName+num.str()); const string unit = Get("TUNIT"+num.str(), ""); const string comp = Get("ZCTYP"+num.str(), ""); FitsCompression compress = UNKNOWN; if (comp == "FACT") compress = FACT; 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; sortedCols.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")) { const streamoff skip = fTable.GetTotalBytes(); seekg(skip, ios_base::cur); fChkHeader.reset(); continue; } //fTable.PrintKeys(); 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); } fits(const string &fname, const string &fout, const string& tableName, bool force=false) : izstream(fname.c_str()) { Constructor(fname, fout, tableName, force); } ~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); } 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; auto ib = fBufferRow.begin(); auto ie = fBufferRow.end(); *ib++ = 0; *ib++ = 0; *ib++ = 0; *ib = 0; *--ie = 0; *--ie = 0; *--ie = 0; *--ie = 0; *--ie = 0; *--ie = 0; *--ie = 0; *--ie = 0; StageRow(row, fBufferRow.data()+offset); WriteRowToCopyFile(row); fRow = row; return offset; } template 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) { 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 (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.sortedCols;} 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(); } bool IsFileOk() const { return (fChkHeader+fChkData).valid(); } bool IsCompressedFITS() const { return fTable.isCompressed;} }; #ifndef __MARS__ }; #endif #endif