#ifndef MARS_ofits #define MARS_ofits #include #include #include #include #include #include #include #include #include #include #ifdef __EXCEPTIONS #include #endif #include "checksum.h" namespace std { // Sloppy: allow / <--- left // allow all characters (see specs for what is possible) // units: m kg s rad sr K A mol cd Hz J W V N Pa C Ohm S F Wb T Hlm lx class ofits : public ofstream { struct Key { string key; bool delim; string value; string comment; off_t offset; // offset in file bool changed; // For closing the file Key(const string &k="") : key(k), delim(false), offset(0), changed(true) { } string Trim(const string &str) { // Trim Both leading and trailing spaces const size_t first = str.find_first_not_of(' '); // Find the first character position after excluding leading blank spaces const size_t last = str.find_last_not_of(' '); // Find the first character position from reverse af // if all spaces or empty return an empty string if (string::npos==first || string::npos==last) return string(); return str.substr(first, last-first+1); } bool FormatKey() { key = Trim(key); if (key.size()==0) { #ifdef __EXCEPTIONS throw runtime_error("Key name empty."); #else gLog << ___err___ << "ERROR - Key name empty." << endl; return false; #endif } if (key.size()>8) { ostringstream sout; sout << "Key '" << key << "' exceeds 8 bytes."; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } //transform(key.begin(), key.end(), key.begin(), toupper); for (string::const_iterator c=key.begin(); c'Z') && (*c<'0' || *c>'9') && *c!='-' && *c!='_') { ostringstream sout; sout << "Invalid character '" << *c << "' found in key '" << key << "'"; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } return true; } bool FormatComment() { comment = Trim(comment); for (string::const_iterator c=key.begin(); c126) { ostringstream sout; sout << "Invalid character '" << *c << "' [" << int(*c) << "] found in comment '" << comment << "'"; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } return true; } bool check() { if (!FormatKey()) return false; if (!FormatComment()) return false; const size_t sz = CalcSize(); if (sz>80) { ostringstream sout; sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters."; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } return true; } size_t CalcSize() const { if (!delim) return 10+comment.size(); return 10 + (value.size()<20?20:value.size()) + 3 + comment.size(); } string Compile() { ostringstream sout; sout << std::left << setw(8) << key; if (!delim) { sout << " " << comment; return sout.str(); } sout << "= "; sout << (value[0]=='\''?std::left:std::right); sout << setw(20) << value << std::left; if (comment.size()>0) sout << " / " << comment; return sout.str(); } Checksum checksum; void Out(ofstream &fout) { if (!changed) return; string str = Compile(); str.insert(str.end(), 80-str.size(), ' '); if (offset==0) offset = fout.tellp(); //cout << "Write[" << offset << "]: " << key << "/" << value << endl; fout.seekp(offset); fout << str; checksum.reset(); checksum.add(str.c_str(), 80); changed = false; } /* void Out(ostream &out) { string str = Compile(); str.insert(str.end(), 80-str.size(), ' '); out << str; changed = false; }*/ }; vector fKeys; vector::iterator findkey(const string &key) { for (auto it=fKeys.begin(); it!=fKeys.end(); it++) if (key==it->key) return it; return fKeys.end(); } bool Set(const string &key="", bool delim=false, const string &value="", const string &comment="") { // If no delimit add the row no matter if it alread exists if (delim) { // if the row already exists: update it auto it = findkey(key); if (it!=fKeys.end()) { it->value = value; it->changed = true; return true; } } if (fTable.num_rows>0) { ostringstream sout; sout << "No new header key can be defined, rows were already written to the file... ignoring new key '" << key << "'"; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } Key entry; entry.key = key; entry.delim = delim; entry.value = value; entry.comment = comment; entry.offset = 0; entry.changed = true; if (!entry.check()) return false; fKeys.push_back(entry); return true; } struct Table { off_t offset; size_t bytes_per_row; size_t num_rows; size_t num_cols; struct Column { string name; size_t offset; size_t num; size_t size; char type; }; vector cols; Table() : offset(0), bytes_per_row(0), num_rows(0), num_cols(0) { } }; Table fTable; vector fOutputBuffer; vector::iterator findcol(const string &name) { for (auto it=fTable.cols.begin(); it!=fTable.cols.end(); it++) if (name==it->name) return it; return fTable.cols.end(); } Checksum fDataSum; Checksum fHeaderSum; public: ofits() { } ofits(const char *fname) : ofstream() { this->open(fname); } ~ofits() { close(); } void open(const char * filename) { fDataSum = 0; fHeaderSum = 0; fTable = Table(); fKeys.clear(); SetStr("XTENSION", "BINTABLE", "binary table extension"); SetInt("BITPIX", 8, "8-bit bytes"); SetInt("NAXIS", 2, "2-dimensional binary table"); SetInt("NAXIS1", 0, "width of table in bytes"); SetInt("NAXIS2", 0, "number of rows in table"); SetInt("PCOUNT", 0, "size of special data area"); SetInt("GCOUNT", 1, "one data group (required keyword)"); SetInt("TFIELDS", 0, "number of fields in each row"); SetStr("EXTNAME", "", "name of extension table"); SetStr("CHECKSUM", "0000000000000000", "Checksum for the whole HDU"); SetStr("DATASUM", " 0", "Checksum for the data block"); ofstream::open(filename); } bool SetBool(const string &key, bool b, const string &comment="") { return Set(key, true, b?"T":"F", comment); } bool AddEmpty(const string &key, const string &comment="") { return Set(key, true, "", comment); } bool SetStr(const string &key, string s, const string &comment="") { for (string::iterator c=s.begin(); c0) sout << setprecision(p); if (p==0) sout << setprecision(f>1e-100 && f<1e100 ? 15 : 14); sout << f; string str = sout.str(); replace(str.begin(), str.end(), 'e', 'E'); if (str.find_first_of('E')==string::npos && str.find_first_of('.')==string::npos) str += "."; return Set(key, true, str, comment); } bool SetFloat(const string &key, double f, const string &comment="") { return SetFloat(key, f, 0, comment); } bool SetHex(const string &key, uint64_t i, const string &comment="") { ostringstream sout; sout << std::hex << "0x" << i; return SetStr(key, sout.str(), comment); } bool AddComment(const string &comment) { return Set("COMMENT", false, "", comment); } bool AddHistory(const string &comment) { return Set("HISTORY", false, "", comment); } void End() { Set("END"); while (fKeys.size()%36!=0) fKeys.push_back(Key()); } bool AddColumn(uint32_t cnt, char typechar, const string &name, const string &unit, const string &comment="") { if (tellp()<0) { ostringstream sout; sout << "File not open... ignoring column '" << name << "'"; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } if (tellp()>0) { ostringstream sout; sout << "Header already writtenm, no new column can be defined... ignoring column '" << name << "'"; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } if (findcol(name)!=fTable.cols.end()) { ostringstream sout; sout << "A column with the name '" << name << "' already exists."; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } typechar = toupper(typechar); static const string allow("LABIJKED"); if (std::find(allow.begin(), allow.end(), typechar)==allow.end()) { ostringstream sout; sout << "Column type '" << typechar << "' not supported."; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } ostringstream type; type << cnt << typechar; fTable.num_cols++; ostringstream typekey, formkey, unitkey, unitcom, typecom; typekey << "TTYPE" << fTable.num_cols; formkey << "TFORM" << fTable.num_cols; unitkey << "TUNIT" << fTable.num_cols; unitcom << "unit of " << name; typecom << "format of " << name << " ["; switch (typechar) { case 'L': typecom << "1-byte BOOL]"; break; case 'A': typecom << "1-byte CHAR]"; break; case 'B': typecom << "1-byte BOOL]"; break; case 'I': typecom << "2-byte INT]"; break; case 'J': typecom << "4-byte INT]"; break; case 'K': typecom << "8-byte INT]"; break; case 'E': typecom << "4-byte FLOAT]"; break; case 'D': typecom << "8-byte FLOAT]"; break; } SetStr(formkey.str(), type.str(), typecom.str()); SetStr(typekey.str(), name, comment); if (!unit.empty()) SetStr(unitkey.str(), unit, unitcom.str()); size_t size = 0; switch (typechar) { case 'L': size = 1; break; case 'A': size = 1; break; case 'B': size = 1; break; case 'I': size = 2; break; case 'J': size = 4; break; case 'K': size = 8; break; case 'E': size = 4; break; case 'D': size = 8; break; } Table::Column col; col.name = name; col.type = typechar; col.num = cnt; col.size = size; col.offset = fTable.bytes_per_row; fTable.cols.push_back(col); fTable.bytes_per_row += size*cnt; // Align to four bytes fOutputBuffer.resize(fTable.bytes_per_row + 4-fTable.bytes_per_row%4); return true; } bool AddColumnShort(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'I', name, unit, comment); } bool AddColumnInt(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'J', name, unit, comment); } bool AddColumnLong(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'K', name, unit, comment); } bool AddColumnFloat(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'E', name, unit, comment); } bool AddColumnDouble(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'D', name, unit, comment); } bool AddColumnChar(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'A', name, unit, comment); } bool AddColumnByte(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'B', name, unit, comment); } bool AddColumnBool(uint32_t cnt, const string &name, const string &unit="", const string &comment="") { return AddColumn(cnt, 'L', name, unit, comment); } bool AddColumnShort(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'I', name, unit, comment); } bool AddColumnInt(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'J', name, unit, comment); } bool AddColumnLong(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'K', name, unit, comment); } bool AddColumnFloat(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'E', name, unit, comment); } bool AddColumnDouble(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'D', name, unit, comment); } bool AddColumnChar(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'A', name, unit, comment); } bool AddColumnByte(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'B', name, unit, comment); } bool AddColumnBool(const string &name, const string &unit="", const string &comment="") { return AddColumn(1, 'L', name, unit, comment); } /* bool AddKey(string key, double d, const string &comment) { ostringstream out; out << d; string s = out.str(); replace(s.begin(), s.end(), "e", "E"); return AddKey(key, s, comment); }*/ Checksum WriteHeader(ofstream &fout) { Checksum sum; for (auto it=fKeys.begin(); it!=fKeys.end(); it++) { it->Out(fout); sum += it->checksum; } fout.flush(); return sum; } Checksum WriteHeader() { return WriteHeader(*this); } void FlushHeader() { const off_t pos = tellp(); WriteHeader(); seekp(pos); } Checksum WriteFitsHeader() { ofits h; h.SetBool("SIMPLE", true, "file does conform to FITS standard"); h.SetInt ("BITPIX", 8, "number of bits per data pixel"); h.SetInt ("NAXIS", 0, "number of data axes"); h.SetBool("EXTEND", true, "FITS dataset may contain extensions"); h.SetStr ("CHECKSUM","0000000000000000", "Checksum for the whole HDU"); h.SetStr ("DATASUM", " 0", "Checksum for the data block"); h.AddComment("FITS (Flexible Image Transport System) format is defined in 'Astronomy"); h.AddComment("and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H"); h.End(); const Checksum sum = h.WriteHeader(*this); h.SetStr("CHECKSUM", sum.str()); const size_t offset = tellp(); h.WriteHeader(*this); seekp(offset); return sum; } bool WriteTableHeader(const char *name="DATA") { if (tellp()>0) { #ifdef __EXCEPTIONS throw runtime_error("File not empty anymore."); #else gLog << ___err___ << "ERROR - File not empty anymore." << endl; return false; #endif } fHeaderSum = WriteFitsHeader(); SetStr("EXTNAME", name); SetInt("NAXIS1", fTable.bytes_per_row); SetInt("TFIELDS", fTable.cols.size()); End(); WriteHeader(); return good(); } template void revcpy(char *dest, const char *src, int num) { const char *pend = src + num*N; for (const char *ptr = src; ptr(ptr) + it->offset; char *dest = buffer + it->offset; // Let the compiler do some optimization by // knowing the we only have 1, 2, 4 and 8 switch (it->size) { case 1: memcpy (dest, src, it->num*it->size); break; case 2: revcpy<2>(dest, src, it->num); break; case 4: revcpy<4>(dest, src, it->num); break; case 8: revcpy<8>(dest, src, it->num); break; } } } write(buffer, cnt); fDataSum.add(fOutputBuffer); fTable.num_rows++; fTable.offset += cnt; return good(); } template bool WriteRow(const vector &vec) { return WriteRow(vec.data(), vec.size()*sizeof(N)); } // Flushes the number of rows to the header on disk void FlushNumRows() { SetInt("NAXIS2", fTable.num_rows); FlushHeader(); } bool close() { if (tellp()<0) return false; if (tellp()%(80*36)>0) { const vector filler(80*36-tellp()%(80*36)); write(filler.data(), filler.size()); } // We don't have to jump back to the end of the file SetInt("NAXIS2", fTable.num_rows); ostringstream dataSumStr; dataSumStr << fDataSum.val(); SetStr("DATASUM", dataSumStr.str()); const Checksum sum = WriteHeader(); //sum += headersum; SetStr("CHECKSUM", (sum+fDataSum).str()); const Checksum chk = WriteHeader(); ofstream::close(); if ((chk+fDataSum).valid()) return true; ostringstream sout; sout << "Checksum (" << std::hex << chk.val() << ") invalid."; #ifdef __EXCEPTIONS throw runtime_error(sout.str()); #else gLog << ___err___ << "ERROR - " << sout.str() << endl; return false; #endif } }; } #if 0 #include "fits.h" int main() { using namespace std; ofits h2("delme.fits"); h2.SetInt("KEY1", 1, "comment 1"); h2.AddColumnInt(2, "testcol1", "counts", "My comment"); h2.AddColumnInt("testcol2", "counts", "My comment"); //h2.AddColumnInt("testcol2", "counts", "My comment"); h2.SetInt("KEY2", 2, "comment 2"); /* AddFloat("X0", 0.000123456, "number of fields in each row"); AddFloat("X1", 0, "number of fields in each row"); AddFloat("X2", 12345, "number of fields in each row"); AddFloat("X3", 123456.67890, "number of fields in each row"); AddFloat("X4", 1234567890123456789.12345678901234567890, "number of fields in each row"); AddFloat("X5", 1234567890.1234567890e20, "number of fields in each row"); AddFloat("X6", 1234567890.1234567890e-20, "number of fields in each row"); AddFloat("XB", 1234567890.1234567890e-111, "number of fields in each row"); AddFloat("X7", 1e-5, "number of fields in each row"); AddFloat("X8", 1e-6, "number of fields in each row"); //AddStr("12345678", "123456789012345678", "12345678901234567890123456789012345678901234567"); */ // - h2.WriteTableHeader("TABLE_NAME"); for (int i=0; i<10; i++) { int j[3] = { i+10, i*10, i*100 }; h2.WriteRow(j, 3*sizeof(i)); } //h2.AddColumnInt("testcol2", "counts", "My comment"); //h2.SetInt("KEY3", 2, "comment 2"); h2.SetInt("KEY2", 2, "comment 2xxx"); h2.SetInt("KEY1", 11); h2.close(); cout << "---" << endl; fits f("delme.fits"); if (!f) throw runtime_error("xxx"); cout << "Header is valid: " << f.IsHeaderOk() << endl; while (f.GetNextRow()); cout << "File is valid: " << f.IsFileOk() << endl; cout << "---" << endl; return 0; } #endif #endif