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

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