source: trunk/Mars/mcore/ofits.h@ 17242

Last change on this file since 17242 was 17240, checked in by tbretz, 12 years ago
Added SetDefaultKeys
File size: 32.3 KB
Line 
1#ifndef MARS_ofits
2#define MARS_ofits
3
4#include <string>
5#include <string.h>
6#include <algorithm>
7#include <sstream>
8#include <iostream>
9#include <fstream>
10#include <iomanip>
11#include <vector>
12#include <algorithm>
13#include <stdexcept>
14
15#ifdef __CINT__
16#define off_t size_t
17#endif
18
19#ifndef __MARS__
20#define gLog cerr
21#define ___err___ ""
22#define ___warn___ ""
23#define ___all___ ""
24#else
25#include "MLog.h"
26#include "MLogManip.h"
27#define ___err___ err
28#define ___warn___ warn
29#define ___all___ all
30#endif
31
32#include "FITS.h"
33#include "checksum.h"
34
35#ifndef __MARS__
36namespace std
37{
38#else
39using namespace std;
40#endif
41
42// Sloppy: allow / <--- left
43// allow all characters (see specs for what is possible)
44
45// units: m kg s rad sr K A mol cd Hz J W V N Pa C Ohm S F Wb T Hlm lx
46
47class ofits : public ofstream
48{
49public:
50 struct Key
51 {
52 string key;
53 bool delim;
54 string value;
55 string comment;
56 string fitsString;
57
58 off_t offset; // offset in file
59
60 bool changed; // For closing the file
61
62 Key(const string &k="") : key(k), delim(false), fitsString(""), offset(0), changed(true) { }
63
64 string Trim(const string &str)
65 {
66 // Trim Both leading and trailing spaces
67 const size_t first = str.find_first_not_of(' '); // Find the first character position after excluding leading blank spaces
68 const size_t last = str.find_last_not_of(' '); // Find the first character position from reverse af
69
70 // if all spaces or empty return an empty string
71 if (string::npos==first || string::npos==last)
72 return string();
73
74 return str.substr(first, last-first+1);
75 }
76
77 bool FormatKey()
78 {
79 key = Trim(key);
80 if (key.empty())
81 {
82#ifdef __EXCEPTIONS
83 throw runtime_error("Key name empty.");
84#else
85 gLog << ___err___ << "ERROR - Key name empty." << endl;
86 return false;
87#endif
88 }
89 if (key.size()>8)
90 {
91 ostringstream sout;
92 sout << "Key '" << key << "' exceeds 8 bytes.";
93#ifdef __EXCEPTIONS
94 throw runtime_error(sout.str());
95#else
96 gLog << ___err___ << "ERROR - " << sout.str() << endl;
97 return false;
98#endif
99 }
100
101 //transform(key.begin(), key.end(), key.begin(), toupper);
102
103 for (string::const_iterator c=key.begin(); c<key.end(); c++)
104 if ((*c<'A' || *c>'Z') && (*c<'0' || *c>'9') && *c!='-' && *c!='_')
105 {
106 ostringstream sout;
107 sout << "Invalid character '" << *c << "' found in key '" << key << "'";
108#ifdef __EXCEPTIONS
109 throw runtime_error(sout.str());
110#else
111 gLog << ___err___ << "ERROR - " << sout.str() << endl;
112 return false;
113#endif
114 }
115
116 return true;
117 }
118
119 bool FormatComment()
120 {
121 comment = Trim(comment);
122
123 for (string::const_iterator c=key.begin(); c<key.end(); c++)
124 if (*c<32 || *c>126)
125 {
126 ostringstream sout;
127 sout << "Invalid character '" << *c << "' [" << int(*c) << "] found in comment '" << comment << "'";
128#ifdef __EXCEPTIONS
129 throw runtime_error(sout.str());
130#else
131 gLog << ___err___ << "ERROR - " << sout.str() << endl;
132 return false;
133#endif
134 }
135
136 return true;
137 }
138
139 bool check(bool trim=false)
140 {
141 if (!FormatKey())
142 return false;
143
144 if (!FormatComment())
145 return false;
146
147 size_t sz = CalcSize();
148 if (sz<=80)
149 return true;
150
151 if (!trim)
152 {
153 ostringstream sout;
154 sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters.";
155#ifdef __EXCEPTIONS
156 throw runtime_error(sout.str());
157#else
158 gLog << ___err___ << "ERROR - " << sout.str() << endl;
159#endif
160 return false;
161 }
162
163 //looks like something went wrong. Maybe entry is too long ?
164 //try to remove the comment
165 comment = "";
166
167 sz = CalcSize();
168 if (sz<=80)
169 {
170#ifndef __EXCEPTIONS
171 ostringstream sout;
172 sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters... removed comment.";
173 gLog << ___warn___ << "WARNING - " << sout.str() << endl;
174#endif
175 return true;
176 }
177
178 ostringstream sout;
179 sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters even without comment.";
180#ifdef __EXCEPTIONS
181 throw runtime_error(sout.str());
182#else
183 gLog << ___err___ << "ERROR - " << sout.str() << endl;
184 return false;
185#endif
186 }
187
188 size_t CalcSize() const
189 {
190 if (!delim)
191 return 10+comment.size();
192
193 return 10 + (value.size()<20?20:value.size()) + 3 + comment.size();
194 }
195
196 string Compile()
197 {
198
199 if (fitsString != "")
200 return fitsString;
201
202 ostringstream sout;
203 sout << std::left << setw(8) << key;
204
205 if (!delim)
206 {
207 sout << " " << comment;
208 return sout.str();
209 }
210
211 sout << "= ";
212 sout << (!value.empty() && value[0]=='\''?std::left:std::right);
213 sout << setw(20) << value << std::left;
214
215 if (!comment.empty())
216 sout << " / " << comment;
217
218 return sout.str();
219 }
220
221 Checksum checksum;
222
223 void Out(ofstream &fout)
224 {
225 if (!changed)
226 return;
227
228 string str = Compile();
229 str.insert(str.end(), 80-str.size(), ' ');
230
231 if (offset==0)
232 offset = fout.tellp();
233
234 //cout << "Write[" << offset << "]: " << key << "/" << value << endl;
235
236 fout.seekp(offset);
237 fout << str;
238
239 checksum.reset();
240 checksum.add(str.c_str(), 80);
241
242 changed = false;
243 }
244 /*
245 void Out(ostream &out)
246 {
247 string str = Compile();
248
249 str.insert(str.end(), 80-str.size(), ' ');
250
251 out << str;
252 changed = false;
253 }*/
254 };
255
256private:
257 vector<Key> fKeys;
258
259 vector<Key>::iterator findkey(const string &key)
260 {
261 for (auto it=fKeys.begin(); it!=fKeys.end(); it++)
262 if (key==it->key)
263 return it;
264
265 return fKeys.end();
266 }
267
268 bool Set(const string &key="", bool delim=false, const string &value="", const string &comment="")
269 {
270 // If no delimit add the row no matter if it alread exists
271 if (delim)
272 {
273 // if the row already exists: update it
274 auto it = findkey(key);
275 if (it!=fKeys.end())
276 {
277 it->value = value;
278 it->changed = true;
279 return true;
280 }
281 }
282
283 if (fTable.num_rows>0)
284 {
285 ostringstream sout;
286 sout << "No new header key can be defined, rows were already written to the file... ignoring new key '" << key << "'";
287#ifdef __EXCEPTIONS
288 throw runtime_error(sout.str());
289#else
290 gLog << ___err___ << "ERROR - " << sout.str() << endl;
291 return false;
292#endif
293 }
294
295 Key entry;
296
297 entry.key = key;
298 entry.delim = delim;
299 entry.value = value;
300 entry.comment = comment;
301 entry.offset = 0;
302 entry.changed = true;
303
304 if (!entry.check(fCommentTrimming))
305 return false;
306
307 fKeys.push_back(entry);
308 return true;
309 }
310
311protected:
312 struct Table
313 {
314 off_t offset;
315
316 size_t bytes_per_row;
317 size_t num_rows;
318 size_t num_cols;
319
320 struct Column
321 {
322 string name;
323 size_t offset;
324 size_t num;
325 size_t size;
326 char type;
327 };
328
329 vector<Column> cols;
330
331 Table() : offset(0), bytes_per_row(0), num_rows(0), num_cols(0)
332 {
333 }
334 };
335
336
337 Table fTable;
338
339 vector<char> fOutputBuffer;
340
341 vector<Table::Column>::iterator findcol(const string &name)
342 {
343 for (auto it=fTable.cols.begin(); it!=fTable.cols.end(); it++)
344 if (name==it->name)
345 return it;
346
347 return fTable.cols.end();
348 }
349
350 Checksum fDataSum;
351 Checksum fHeaderSum;
352
353 bool fCommentTrimming;
354 bool fManualExtName;
355
356public:
357 ofits() : fCommentTrimming(false),
358 fManualExtName(false)
359 {
360 }
361 ofits(const char *fname) : ofstream(),
362 fCommentTrimming(false),
363 fManualExtName(false)
364 {
365 this->open(fname);
366 }
367 virtual ~ofits() { if (is_open()) close(); }
368
369 virtual void open(const char * filename, bool addEXTNAMEKey=true)
370 {
371 fDataSum = 0;
372 fHeaderSum = 0;
373
374 fTable = Table();
375 fKeys.clear();
376
377 SetStr("XTENSION", "BINTABLE", "binary table extension");
378 SetInt("BITPIX", 8, "8-bit bytes");
379 SetInt("NAXIS", 2, "2-dimensional binary table");
380 SetInt("NAXIS1", 0, "width of table in bytes");
381 SetInt("NAXIS2", 0, "number of rows in table");
382 SetInt("PCOUNT", 0, "size of special data area");
383 SetInt("GCOUNT", 1, "one data group (required keyword)");
384 SetInt("TFIELDS", 0, "number of fields in each row");
385 if (addEXTNAMEKey)
386 SetStr("EXTNAME", "", "name of extension table");
387 else
388 fManualExtName = true;
389 SetStr("CHECKSUM", "0000000000000000", "Checksum for the whole HDU");
390 SetStr("DATASUM", " 0", "Checksum for the data block");
391
392 ofstream::open(filename);
393 }
394
395 void AllowCommentsTrimming(bool allow)
396 {
397 fCommentTrimming = allow;
398 }
399 //Etienne: required to enable 1 to 1 reconstruction of files
400 bool SetKeyComment(const string& key, const string& comment)
401 {
402 auto it = findkey(key);
403 if (it==fKeys.end())
404 return false;
405 it->comment = comment;
406 it->changed = true;
407 return true;
408 }
409 bool SetKeyFromFitsString(const string& fitsString)
410 {
411 if (fTable.num_rows>0)
412 {
413 ostringstream sout;
414 sout << "No new header key can be defined, rows were already written to the file... ignoring new key '" << fitsString << "'";
415#ifdef __EXCEPTIONS
416 throw runtime_error(sout.str());
417#else
418 gLog << ___err___ << "ERROR - " << sout.str() << endl;
419 return false;
420#endif
421 }
422
423 Key entry;
424 entry.fitsString = fitsString;
425 entry.changed = true;
426 fKeys.push_back(entry);
427 return true;
428 }
429 bool SetRaw(const string &key, const string &val, const string &comment)
430 {
431 return Set(key, true, val, comment);
432 }
433
434 bool SetBool(const string &key, bool b, const string &comment="")
435 {
436 return Set(key, true, b?"T":"F", comment);
437 }
438
439 bool AddEmpty(const string &key, const string &comment="")
440 {
441 return Set(key, true, "", comment);
442 }
443
444 bool SetStr(const string &key, string s, const string &comment="")
445 {
446 for (uint i=0; i<s.length(); i++)
447 if (s[i]=='\'')
448 s.insert(i++, "\'");
449
450 return Set(key, true, "'"+s+"'", comment);
451 }
452
453 bool SetInt(const string &key, int64_t i, const string &comment="")
454 {
455 ostringstream sout;
456 sout << i;
457
458 return Set(key, true, sout.str(), comment);
459 }
460
461 bool SetFloat(const string &key, double f, int p, const string &comment="")
462 {
463 ostringstream sout;
464
465 if (p<0)
466 sout << setprecision(-p) << fixed;
467 if (p>0)
468 sout << setprecision(p);
469 if (p==0)
470 sout << setprecision(f>1e-100 && f<1e100 ? 15 : 14);
471
472 sout << f;
473
474 string str = sout.str();
475
476 replace(str.begin(), str.end(), 'e', 'E');
477
478 if (str.find_first_of('E')==string::npos && str.find_first_of('.')==string::npos)
479 str += ".";
480
481 return Set(key, true, str, comment);
482 }
483
484 bool SetFloat(const string &key, double f, const string &comment="")
485 {
486 return SetFloat(key, f, 0, comment);
487 }
488
489 bool SetHex(const string &key, uint64_t i, const string &comment="")
490 {
491 ostringstream sout;
492 sout << std::hex << "0x" << i;
493 return SetStr(key, sout.str(), comment);
494 }
495
496 bool AddComment(const string &comment)
497 {
498 return Set("COMMENT", false, "", comment);
499 }
500
501 bool AddHistory(const string &comment)
502 {
503 return Set("HISTORY", false, "", comment);
504 }
505
506 void End()
507 {
508 Set("END");
509 while (fKeys.size()%36!=0)
510 fKeys.emplace_back();
511 }
512
513 string CommentFromType(char type)
514 {
515 string comment;
516
517 switch (type)
518 {
519 case 'L': comment = "1-byte BOOL]"; break;
520 case 'A': comment = "1-byte CHAR]"; break;
521 case 'B': comment = "1-byte BOOL]"; break;
522 case 'I': comment = "2-byte INT]"; break;
523 case 'J': comment = "4-byte INT]"; break;
524 case 'K': comment = "8-byte INT]"; break;
525 case 'E': comment = "4-byte FLOAT]"; break;
526 case 'D': comment = "8-byte FLOAT]"; break;
527 case 'Q': comment = "var. Length]"; break;
528 }
529
530 return comment;
531 }
532
533 uint32_t SizeFromType(char type)
534 {
535 size_t size = 0;
536
537 switch (type)
538 {
539 case 'L': size = 1; break;
540 case 'A': size = 1; break;
541 case 'B': size = 1; break;
542 case 'I': size = 2; break;
543 case 'J': size = 4; break;
544 case 'K': size = 8; break;
545 case 'E': size = 4; break;
546 case 'D': size = 8; break;
547 case 'Q': size = 16; break;
548 }
549
550 return size;
551 }
552 //ETIENNE to be able to restore the file 1 to 1, I must restore the header keys myself
553 virtual bool AddColumn(uint32_t cnt, char typechar, const string &name, const string &unit, const string &comment="", bool addHeaderKeys=true)
554 {
555 if (tellp()<0)
556 {
557 ostringstream sout;
558 sout << "File not open... ignoring column '" << name << "'";
559#ifdef __EXCEPTIONS
560 throw runtime_error(sout.str());
561#else
562 gLog << ___err___ << "ERROR - " << sout.str() << endl;
563 return false;
564#endif
565 }
566
567 if (tellp()>0)
568 {
569 ostringstream sout;
570 sout << "Header already written, no new column can be defined... ignoring column '" << name << "'";
571#ifdef __EXCEPTIONS
572 throw runtime_error(sout.str());
573#else
574 gLog << ___err___ << "ERROR - " << sout.str() << endl;
575 return false;
576#endif
577 }
578
579 if (findcol(name)!=fTable.cols.end())
580 {
581 ostringstream sout;
582 sout << "A column with the name '" << name << "' already exists.";
583#ifdef __EXCEPTIONS
584 throw runtime_error(sout.str());
585#else
586 gLog << ___err___ << "ERROR - " << sout.str() << endl;
587 return false;
588#endif
589 }
590
591 typechar = toupper(typechar);
592
593 static const string allow("LABIJKEDQ");
594 if (std::find(allow.begin(), allow.end(), typechar)==allow.end())
595 {
596 ostringstream sout;
597 sout << "Column type '" << typechar << "' not supported.";
598#ifdef __EXCEPTIONS
599 throw runtime_error(sout.str());
600#else
601 gLog << ___err___ << "ERROR - " << sout.str() << endl;
602 return false;
603#endif
604 }
605
606 ostringstream type;
607 type << cnt;
608 if (typechar=='Q')
609 type << "QB";
610 else
611 type << typechar;
612
613 fTable.num_cols++;
614
615 ostringstream typekey, formkey, unitkey, unitcom, typecom;
616 typekey << "TTYPE" << fTable.num_cols;
617 formkey << "TFORM" << fTable.num_cols;
618 unitkey << "TUNIT" << fTable.num_cols;
619 unitcom << "unit of " << name;
620
621 typecom << "format of " << name << " [";
622
623 typecom << CommentFromType(typechar);
624
625 if (addHeaderKeys)
626 {
627 SetStr(formkey.str(), type.str(), typecom.str());
628 SetStr(typekey.str(), name, comment);
629 if (!unit.empty())
630 SetStr(unitkey.str(), unit, unitcom.str());
631 }
632 size_t size = SizeFromType(typechar);
633
634 Table::Column col;
635
636 col.name = name;
637 col.type = typechar;
638 col.num = cnt;
639 col.size = size;
640 col.offset = fTable.bytes_per_row;
641
642 fTable.cols.push_back(col);
643
644 fTable.bytes_per_row += size*cnt;
645
646 // Align to four bytes
647 fOutputBuffer.resize(fTable.bytes_per_row + 4-fTable.bytes_per_row%4);
648
649 return true;
650 }
651
652 virtual bool AddColumn(const FITS::Compression&, uint32_t cnt, char typechar, const string& name, const string& unit, const string& comment="", bool addHeaderKeys=true)
653 {
654 return AddColumn(cnt, typechar, name, unit, comment, addHeaderKeys);
655 }
656
657 bool AddColumnShort(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
658 { return AddColumn(cnt, 'I', name, unit, comment); }
659 bool AddColumnInt(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
660 { return AddColumn(cnt, 'J', name, unit, comment); }
661 bool AddColumnLong(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
662 { return AddColumn(cnt, 'K', name, unit, comment); }
663 bool AddColumnFloat(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
664 { return AddColumn(cnt, 'E', name, unit, comment); }
665 bool AddColumnDouble(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
666 { return AddColumn(cnt, 'D', name, unit, comment); }
667 bool AddColumnChar(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
668 { return AddColumn(cnt, 'A', name, unit, comment); }
669 bool AddColumnByte(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
670 { return AddColumn(cnt, 'B', name, unit, comment); }
671 bool AddColumnBool(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
672 { return AddColumn(cnt, 'L', name, unit, comment); }
673
674 bool AddColumnShort(const string &name, const string &unit="", const string &comment="")
675 { return AddColumn(1, 'I', name, unit, comment); }
676 bool AddColumnInt(const string &name, const string &unit="", const string &comment="")
677 { return AddColumn(1, 'J', name, unit, comment); }
678 bool AddColumnLong(const string &name, const string &unit="", const string &comment="")
679 { return AddColumn(1, 'K', name, unit, comment); }
680 bool AddColumnFloat(const string &name, const string &unit="", const string &comment="")
681 { return AddColumn(1, 'E', name, unit, comment); }
682 bool AddColumnDouble(const string &name, const string &unit="", const string &comment="")
683 { return AddColumn(1, 'D', name, unit, comment); }
684 bool AddColumnChar(const string &name, const string &unit="", const string &comment="")
685 { return AddColumn(1, 'A', name, unit, comment); }
686 bool AddColumnByte(const string &name, const string &unit="", const string &comment="")
687 { return AddColumn(1, 'B', name, unit, comment); }
688 bool AddColumnBool(const string &name, const string &unit="", const string &comment="")
689 { return AddColumn(1, 'L', name, unit, comment); }
690
691 bool AddColumnShort(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
692 { return AddColumn(comp, cnt, 'I', name, unit, comment); }
693 bool AddColumnInt(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
694 { return AddColumn(comp, cnt, 'J', name, unit, comment); }
695 bool AddColumnLong(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
696 { return AddColumn(comp, cnt, 'K', name, unit, comment); }
697 bool AddColumnFloat(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
698 { return AddColumn(comp, cnt, 'E', name, unit, comment); }
699 bool AddColumnDouble(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
700 { return AddColumn(comp, cnt, 'D', name, unit, comment); }
701 bool AddColumnChar(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
702 { return AddColumn(comp, cnt, 'A', name, unit, comment); }
703 bool AddColumnByte(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
704 { return AddColumn(comp, cnt, 'B', name, unit, comment); }
705 bool AddColumnBool(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
706 { return AddColumn(comp, cnt, 'L', name, unit, comment); }
707
708 bool AddColumnShort(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
709 { return AddColumn(comp, 1, 'I', name, unit, comment); }
710 bool AddColumnInt(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
711 { return AddColumn(comp, 1, 'J', name, unit, comment); }
712 bool AddColumnLong(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
713 { return AddColumn(comp, 1, 'K', name, unit, comment); }
714 bool AddColumnFloat(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
715 { return AddColumn(comp, 1, 'E', name, unit, comment); }
716 bool AddColumnDouble(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
717 { return AddColumn(comp, 1, 'D', name, unit, comment); }
718 bool AddColumnChar(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
719 { return AddColumn(comp, 1, 'A', name, unit, comment); }
720 bool AddColumnByte(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
721 { return AddColumn(comp, 1, 'B', name, unit, comment); }
722 bool AddColumnBool(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
723 { return AddColumn(comp, 1, 'L', name, unit, comment); }
724
725 /*
726 bool AddKey(string key, double d, const string &comment)
727 {
728 ostringstream out;
729 out << d;
730
731 string s = out.str();
732
733 replace(s.begin(), s.end(), "e", "E");
734
735 return AddKey(key, s, comment);
736 }*/
737
738
739 Checksum WriteHeader(ofstream &fout)
740 {
741 Checksum sum;
742 uint32_t count=0;
743 for (auto it=fKeys.begin(); it!=fKeys.end(); it++)
744 {
745 it->Out(fout);
746 sum += it->checksum;
747 count++;
748 }
749 fout.flush();
750
751 return sum;
752 }
753
754 Checksum WriteHeader()
755 {
756 return WriteHeader(*this);
757 }
758
759 void FlushHeader()
760 {
761 const off_t pos = tellp();
762 WriteHeader();
763 seekp(pos);
764 }
765
766 Checksum WriteFitsHeader()
767 {
768 ofits h;
769
770 h.SetBool("SIMPLE", true, "file does conform to FITS standard");
771 h.SetInt ("BITPIX", 8, "number of bits per data pixel");
772 h.SetInt ("NAXIS", 0, "number of data axes");
773 h.SetBool("EXTEND", true, "FITS dataset may contain extensions");
774 h.SetStr ("CHECKSUM","0000000000000000", "Checksum for the whole HDU");
775 h.SetStr ("DATASUM", " 0", "Checksum for the data block");
776 h.AddComment("FITS (Flexible Image Transport System) format is defined in 'Astronomy");
777 h.AddComment("and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H");
778 h.End();
779
780 const Checksum sum = h.WriteHeader(*this);
781
782 h.SetStr("CHECKSUM", sum.str());
783
784 const size_t offset = tellp();
785 h.WriteHeader(*this);
786 seekp(offset);
787
788 return sum;
789 }
790
791 virtual bool WriteDrsOffsetsTable ()
792 {
793 return true;
794 }
795
796 virtual bool WriteCatalog()
797 {
798 return true;
799 }
800
801 virtual bool WriteTableHeader(const char *name="DATA")
802 {
803 if (tellp()>0)
804 {
805#ifdef __EXCEPTIONS
806 throw runtime_error("File not empty anymore.");
807#else
808 gLog << ___err___ << "ERROR - File not empty anymore." << endl;
809 return false;
810#endif
811 }
812
813 fHeaderSum = WriteFitsHeader();
814
815 WriteDrsOffsetsTable();
816
817 if (!fManualExtName)
818 SetStr("EXTNAME", name);
819 SetInt("NAXIS1", fTable.bytes_per_row);
820 SetInt("TFIELDS", fTable.cols.size());
821
822 End();
823
824 WriteHeader();
825
826 WriteCatalog();
827
828 return good();
829 }
830
831 template<size_t N>
832 void revcpy(char *dest, const char *src, int num)
833 {
834 const char *pend = src + num*N;
835 for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
836 reverse_copy(ptr, ptr+N, dest);
837 }
838
839 virtual uint32_t GetBytesPerRow() const { return fTable.bytes_per_row; }
840
841 virtual bool WriteRow(const void *ptr, size_t cnt, bool byte_swap=true)
842 {
843 // FIXME: Make sure that header was already written
844 // or write header now!
845 if (cnt!=fTable.bytes_per_row)
846 {
847 ostringstream sout;
848 sout << "WriteRow - Size " << cnt << " does not match expected size " << fTable.bytes_per_row;
849#ifdef __EXCEPTIONS
850 throw runtime_error(sout.str());
851#else
852 gLog << ___err___ << "ERROR - " << sout.str() << endl;
853 return false;
854#endif
855 }
856
857 // For the checksum we need everything to be correctly aligned
858 const uint8_t offset = fTable.offset%4;
859
860 char *buffer = fOutputBuffer.data() + offset;
861
862 auto ib = fOutputBuffer.begin();
863 auto ie = fOutputBuffer.rbegin();
864 *ib++ = 0;
865 *ib++ = 0;
866 *ib++ = 0;
867 *ib = 0;
868
869 *ie++ = 0;
870 *ie++ = 0;
871 *ie++ = 0;
872 *ie = 0;
873
874 if (!byte_swap)
875 memcpy(buffer, ptr, cnt);
876 else
877 {
878 for (auto it=fTable.cols.begin(); it!=fTable.cols.end(); it++)
879 {
880 const char *src = reinterpret_cast<const char*>(ptr) + it->offset;
881 char *dest = buffer + it->offset;
882
883 // Let the compiler do some optimization by
884 // knowing the we only have 1, 2, 4 and 8
885 switch (it->size)
886 {
887 case 1: memcpy (dest, src, it->num*it->size); break;
888 case 2: revcpy<2>(dest, src, it->num); break;
889 case 4: revcpy<4>(dest, src, it->num); break;
890 case 8: revcpy<8>(dest, src, it->num); break;
891 }
892 }
893 }
894
895 write(buffer, cnt);
896 fDataSum.add(fOutputBuffer);
897
898 fTable.num_rows++;
899 fTable.offset += cnt;
900 return good();
901 }
902
903 template<typename N>
904 bool WriteRow(const vector<N> &vec)
905 {
906 return WriteRow(vec.data(), vec.size()*sizeof(N));
907 }
908
909 // Flushes the number of rows to the header on disk
910 virtual void FlushNumRows()
911 {
912 SetInt("NAXIS2", fTable.num_rows);
913 FlushHeader();
914 }
915
916 size_t GetNumRows() const { return fTable.num_rows; }
917
918 void AlignTo2880Bytes()
919 {
920 if (tellp()%(80*36)>0)
921 {
922 vector<char> filler(80*36-tellp()%(80*36));
923 write(filler.data(), filler.size());
924 }
925 }
926
927 Checksum UpdateHeaderChecksum()
928 {
929
930 ostringstream dataSumStr;
931 dataSumStr << fDataSum.val();
932 SetStr("DATASUM", dataSumStr.str());
933
934 const Checksum sum = WriteHeader();
935
936 //sum += headersum;
937
938 SetStr("CHECKSUM", (sum+fDataSum).str());
939
940 return WriteHeader();
941 }
942 virtual bool close()
943 {
944 if (tellp()<0)
945 return false;
946
947 AlignTo2880Bytes();
948
949 // We don't have to jump back to the end of the file
950 SetInt("NAXIS2", fTable.num_rows);
951
952
953 const Checksum chk = UpdateHeaderChecksum();
954
955 ofstream::close();
956
957 if ((chk+fDataSum).valid())
958 return true;
959
960 ostringstream sout;
961 sout << "Checksum (" << std::hex << chk.val() << ") invalid.";
962#ifdef __EXCEPTIONS
963 throw runtime_error(sout.str());
964#else
965 gLog << ___err___ << "ERROR - " << sout.str() << endl;
966 return false;
967#endif
968 }
969
970 pair<string, int> GetChecksumData()
971 {
972 string datasum;
973 string checksum;
974 //cannot directly use the Get methods, because they are only in fits.h
975 for (vector<Key>::const_iterator it=fKeys.begin(); it!= fKeys.end(); it++)
976 {
977 if (it->key == "CHECKSUM") checksum = it->value;
978 if (it->key == "DATASUM") datasum = it->value;
979 }
980 if (checksum[0] == '\'')
981 checksum = checksum.substr(1,checksum.size()-2);
982 if (datasum[0] == '\'')
983 datasum = datasum.substr(1, datasum.size()-2);
984 return make_pair(checksum, atoi(datasum.c_str()));
985 }
986
987 void SetDefaultKeys()
988 {
989 SetStr("TELESCOP", "FACT", "Telescope that acquired this data");
990 SetStr("CREATOR", typeid(*this).name(), "Class that wrote this file");
991 SetFloat("EXTREL", 1.0, "Release Number");
992 SetStr("COMPILED", __DATE__" "__TIME__, "Compile time");
993 SetStr("ORIGIN", "FACT", "Institution that wrote the file");
994 SetStr("TIMESYS", "UTC", "Time system");
995 SetStr("TIMEUNIT", "d", "Time given in days w.r.t. to MJDREF");
996 SetInt("MJDREF", 40587, "MJD to UNIX time (seconds since 1970/1/1)");
997 SetStr("PACKAGE", PACKAGE_NAME, "Package name");
998 SetStr("VERSION", PACKAGE_VERSION, "Package description");
999 SetStr("REVISION", REVISION, "SVN revision");
1000
1001 const time_t t0 = time(NULL);
1002 const struct tm *tmp1 = gmtime(&t0);
1003
1004 string str(21, '\0');
1005 if (tmp1 && strftime(const_cast<char*>(str.data()), 20, "%Y-%m-%dT%H:%M:%S", tmp1))
1006 SetStr("DATE", str, "File creation date");
1007 }
1008};
1009
1010#ifndef __MARS__
1011};
1012#endif
1013
1014#if 0
1015#include "fits.h"
1016
1017int main()
1018{
1019 using namespace std;
1020
1021 ofits h2("delme.fits");
1022
1023 h2.SetInt("KEY1", 1, "comment 1");
1024 h2.AddColumnInt(2, "testcol1", "counts", "My comment");
1025 h2.AddColumnInt("testcol2", "counts", "My comment");
1026 //h2.AddColumnInt("testcol2", "counts", "My comment");
1027 h2.SetInt("KEY2", 2, "comment 2");
1028
1029 /*
1030 AddFloat("X0", 0.000123456, "number of fields in each row");
1031 AddFloat("X1", 0, "number of fields in each row");
1032 AddFloat("X2", 12345, "number of fields in each row");
1033 AddFloat("X3", 123456.67890, "number of fields in each row");
1034 AddFloat("X4", 1234567890123456789.12345678901234567890, "number of fields in each row");
1035 AddFloat("X5", 1234567890.1234567890e20, "number of fields in each row");
1036 AddFloat("X6", 1234567890.1234567890e-20, "number of fields in each row");
1037 AddFloat("XB", 1234567890.1234567890e-111, "number of fields in each row");
1038 AddFloat("X7", 1e-5, "number of fields in each row");
1039 AddFloat("X8", 1e-6, "number of fields in each row");
1040 //AddStr("12345678", "123456789012345678", "12345678901234567890123456789012345678901234567");
1041 */
1042 // -
1043
1044 h2.WriteTableHeader("TABLE_NAME");
1045
1046 for (int i=0; i<10; i++)
1047 {
1048 int j[3] = { i+10, i*10, i*100 };
1049 h2.WriteRow(j, 3*sizeof(i));
1050 }
1051
1052 //h2.AddColumnInt("testcol2", "counts", "My comment");
1053 //h2.SetInt("KEY3", 2, "comment 2");
1054 h2.SetInt("KEY2", 2, "comment 2xxx");
1055 h2.SetInt("KEY1", 11);
1056
1057 h2.close();
1058
1059 cout << "---" << endl;
1060
1061 fits f("delme.fits");
1062 if (!f)
1063 throw runtime_error("xxx");
1064
1065 cout << "Header is valid: " << f.IsHeaderOk() << endl;
1066
1067 while (f.GetNextRow());
1068
1069 cout << "File is valid: " << f.IsFileOk() << endl;
1070
1071 cout << "---" << endl;
1072
1073 return 0;
1074}
1075#endif
1076
1077#endif
Note: See TracBrowser for help on using the repository browser.