source: trunk/Mars/mcore/fits.h@ 15215

Last change on this file since 15215 was 15085, checked in by tbretz, 13 years ago
To 'solve' the problem with the unordered_map MARS will now use a std::map. Fixed a problem which caused the ReadRow to write below the buffer - weird that this never resulted in any problem.
File size: 24.4 KB
Line 
1#ifndef MARS_fits
2#define MARS_fits
3
4#ifdef __CINT__
5#define int8_t Char_t
6#define int16_t Short_t
7#define int32_t Int_t
8#define int64_t Long64_t
9#define uint8_t UChar_t
10#define uint16_t UShort_t
11#define uint32_t UInt_t
12#define uint64_t ULong64_t
13#else
14#include <stdint.h>
15#endif
16
17#include <map>
18#include <string>
19#include <sstream>
20#include <algorithm>
21
22#ifdef __EXCEPTIONS
23#include <stdexcept>
24#endif
25
26#ifdef __CINT__
27#define off_t size_t
28#else
29#include <unordered_map>
30#endif
31
32#ifndef __MARS__
33#include <vector>
34#include <iomanip>
35#include <iostream>
36#define gLog cerr
37#define ___err___ ""
38#define ___warn___ ""
39#define ___all___ ""
40#else
41#include "MLog.h"
42#include "MLogManip.h"
43#define ___err___ err
44#define ___warn___ warn
45#define ___all___ all
46#endif
47
48#ifdef HAVE_ZLIB
49#include "izstream.h"
50#else
51#include <fstream>
52#define izstream ifstream
53#warning Support for zipped FITS files disabled.
54#endif
55
56#include "checksum.h"
57
58#ifndef __MARS__
59namespace std
60{
61#else
62using namespace std;
63#endif
64
65class fits : public izstream
66{
67public:
68 struct Entry
69 {
70 char type;
71 string value;
72 string comment;
73
74 template<typename T>
75 T Get() const
76 {
77 T t;
78
79 istringstream str(value);
80 str >> t;
81
82 return t;
83 }
84 };
85
86 struct Table
87 {
88 off_t offset;
89
90 string name;
91 size_t bytes_per_row;
92 size_t num_rows;
93 size_t num_cols;
94
95 struct Column
96 {
97 size_t offset;
98 size_t num;
99 size_t size;
100 char type;
101 string unit;
102 };
103
104 typedef map<string, Entry> Keys;
105 typedef map<string, Column> Columns;
106
107 Columns cols;
108 Keys keys;
109
110 int64_t datasum;
111
112 string Trim(const string &str, char c=' ') const
113 {
114 // Trim Both leading and trailing spaces
115 const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces
116 const size_t pend = str.find_last_not_of(c); // Find the first character position from reverse af
117
118 // if all spaces or empty return an empty string
119 if (string::npos==pstart || string::npos==pend)
120 return string();
121
122 return str.substr(pstart, pend-pstart+1);
123 }
124
125 bool Check(const string &key, char type, const string &value="") const
126 {
127 const Keys::const_iterator it = keys.find(key);
128 if (it==keys.end())
129 {
130 ostringstream str;
131 str << "Key '" << key << "' not found.";
132#ifdef __EXCEPTIONS
133 throw runtime_error(str.str());
134#else
135 gLog << ___err___ << "ERROR - " << str.str() << endl;
136 return false;
137#endif
138 }
139
140 if (it->second.type!=type)
141 {
142 ostringstream str;
143 str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << ".";
144#ifdef __EXCEPTIONS
145 throw runtime_error(str.str());
146#else
147 gLog << ___err___ << "ERROR - " << str.str() << endl;
148 return false;
149#endif
150 }
151
152 if (!value.empty() && it->second.value!=value)
153 {
154 ostringstream str;
155 str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << ".";
156#ifdef __EXCEPTIONS
157 throw runtime_error(str.str());
158#else
159 gLog << ___err___ << "ERROR - " << str.str() << endl;
160 return false;
161#endif
162 }
163
164 return true;
165 }
166
167 Keys ParseBlock(const vector<string> &vec) const
168 {
169 map<string,Entry> rc;
170
171 for (unsigned int i=0; i<vec.size(); i++)
172 {
173 const string key = Trim(vec[i].substr(0,8));
174 // Keywords without a value, like COMMENT / HISTORY
175 if (vec[i].substr(8,2)!="= ")
176 continue;
177
178 char type = 0;
179
180 string com;
181 string val = Trim(vec[i].substr(10));
182 if (val[0]=='\'')
183 {
184 // First skip all '' in the string
185 size_t p = 1;
186 while (1)
187 {
188 const size_t pp = val.find_first_of('\'', p);
189 if (pp==string::npos)
190 break;
191
192 p = val[pp+1]=='\'' ? pp+2 : pp+1;
193 }
194
195 // Now find the comment
196 const size_t ppp = val.find_first_of('/', p);
197
198 // Set value, comment and type
199 com = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
200 val = Trim(val.substr(1, p-2));
201 type = 'T';
202 }
203 else
204 {
205 const size_t p = val.find_first_of('/');
206
207 com = Trim(val.substr(p+2));
208 val = Trim(val.substr(0, p));
209
210 if (val.empty() || val.find_first_of('T')!=string::npos || val.find_first_of('F')!=string::npos)
211 type = 'B';
212 else
213 type = val.find_last_of('.')==string::npos ? 'I' : 'F';
214 }
215
216 const Entry e = { type, val, com };
217 rc[key] = e;
218 }
219
220 return rc;
221 }
222
223 Table() : offset(0) { }
224 Table(const vector<string> &vec, off_t off) :
225 offset(off), keys(ParseBlock(vec))
226 {
227 if (!Check("XTENSION", 'T', "BINTABLE") ||
228 !Check("NAXIS", 'I', "2") ||
229 !Check("BITPIX", 'I', "8") ||
230 !Check("PCOUNT", 'I', "0") ||
231 !Check("GCOUNT", 'I', "1") ||
232 !Check("EXTNAME", 'T') ||
233 !Check("NAXIS1", 'I') ||
234 !Check("NAXIS2", 'I') ||
235 !Check("TFIELDS", 'I'))
236 return;
237
238 bytes_per_row = Get<size_t>("NAXIS1");
239 num_rows = Get<size_t>("NAXIS2");
240 num_cols = Get<size_t>("TFIELDS");
241 datasum = Get<int64_t>("DATASUM", -1);
242
243 size_t bytes = 0;
244 for (size_t i=1; i<=num_cols; i++)
245 {
246 ostringstream num;
247 num << i;
248
249 if (!Check("TTYPE"+num.str(), 'T') ||
250 !Check("TFORM"+num.str(), 'T'))
251 return;
252
253 const string id = Get<string>("TTYPE"+num.str());
254 const string fmt = Get<string>("TFORM"+num.str());
255 const string unit = Get<string>("TUNIT"+num.str(), "");
256
257 istringstream sin(fmt);
258 int n = 0;
259 sin >> n;
260 if (!sin)
261 n = 1;
262
263 const char type = fmt[fmt.length()-1];
264
265 size_t size = 0;
266 switch (type)
267 {
268 // We could use negative values to mark floats
269 // otheriwse we could just cast them to int64_t?
270 case 'L': // logical
271 case 'A': // char
272 case 'B': size = 1; break; // byte
273 case 'I': size = 2; break; // short
274 case 'J': size = 4; break; // int
275 case 'K': size = 8; break; // long long
276 case 'E': size = 4; break; // float
277 case 'D': size = 8; break; // double
278 // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits)
279 // case 'C': size = 8; break; // complex float
280 // case 'M': size = 16; break; // complex double
281 // case 'P': size = 8; break; // array descriptor (32bit)
282 // case 'Q': size = 16; break; // array descriptor (64bit)
283 default:
284 {
285 ostringstream str;
286 str << "FITS format TFORM='" << fmt << "' not yet supported.";
287#ifdef __EXCEPTIONS
288 throw runtime_error(str.str());
289#else
290 gLog << ___err___ << "ERROR - " << str.str() << endl;
291 return;
292#endif
293 }
294 }
295
296 const Table::Column col = { bytes, n, size, type, unit };
297
298 cols[id] = col;
299 bytes += n*size;
300 }
301
302 if (bytes!=bytes_per_row)
303 {
304#ifdef __EXCEPTIONS
305 throw runtime_error("Column size mismatch");
306#else
307 gLog << ___err___ << "ERROR - Column size mismatch" << endl;
308 return;
309#endif
310 }
311
312 name = Get<string>("EXTNAME");
313 }
314
315 void PrintKeys(bool display_all=false) const
316 {
317 for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
318 {
319 if (!display_all &&
320 (it->first.substr(0, 6)=="TTYPE" ||
321 it->first.substr(0, 6)=="TFORM" ||
322 it->first.substr(0, 6)=="TUNIT" ||
323 it->first=="TFIELDS" ||
324 it->first=="XTENSION" ||
325 it->first=="NAXIS" ||
326 it->first=="BITPIX" ||
327 it->first=="PCOUNT" ||
328 it->first=="GCOUNT")
329 )
330 continue;
331
332 gLog << ___all___ << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl;
333 }}
334
335 void PrintColumns() const
336 {
337 typedef map<pair<size_t, string>, Column> Sorted;
338
339 Sorted sorted;
340
341 for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++)
342 sorted[make_pair(it->second.offset, it->first)] = it->second;
343
344 for (Sorted::const_iterator it=sorted.begin(); it!=sorted.end(); it++)
345 {
346 gLog << ___all___ << setw(6) << it->second.offset << "| ";
347 gLog << it->second.num << 'x';
348 switch (it->second.type)
349 {
350 case 'A': gLog << "char(8)"; break;
351 case 'L': gLog << "bool(8)"; break;
352 case 'B': gLog << "byte(8)"; break;
353 case 'I': gLog << "short(16)"; break;
354 case 'J': gLog << "int(32)"; break;
355 case 'K': gLog << "int(64)"; break;
356 case 'E': gLog << "float(32)"; break;
357 case 'D': gLog << "double(64)"; break;
358 }
359 gLog << ": " << it->first.second << " [" << it->second.unit << "]" << endl;
360 }
361 }
362
363 operator bool() const { return !name.empty(); }
364
365 bool HasKey(const string &key) const
366 {
367 return keys.find(key)!=keys.end();
368 }
369
370 bool HasColumn(const string& col) const
371 {
372 return cols.find(col)!=cols.end();
373 }
374
375 const Columns &GetColumns() const
376 {
377 return cols;
378 }
379
380 const Keys &GetKeys() const
381 {
382 return keys;
383 }
384
385 // Values of keys are always signed
386 template<typename T>
387 T Get(const string &key) const
388 {
389 const map<string,Entry>::const_iterator it = keys.find(key);
390 if (it==keys.end())
391 {
392 ostringstream str;
393 str << "Key '" << key << "' not found." << endl;
394#ifdef __EXCEPTIONS
395 throw runtime_error(str.str());
396#else
397 gLog << ___err___ << "ERROR - " << str.str() << endl;
398 return T();
399#endif
400 }
401 return it->second.Get<T>();
402 }
403
404 // Values of keys are always signed
405 template<typename T>
406 T Get(const string &key, const T &deflt) const
407 {
408 const map<string,Entry>::const_iterator it = keys.find(key);
409 return it==keys.end() ? deflt :it->second.Get<T>();
410 }
411
412 size_t GetN(const string &key) const
413 {
414 const Columns::const_iterator it = cols.find(key);
415 return it==cols.end() ? 0 : it->second.num;
416 }
417 };
418
419private:
420 Table fTable;
421
422 typedef pair<void*, Table::Column> Address;
423 typedef vector<Address> Addresses;
424 //map<void*, Table::Column> fAddresses;
425 Addresses fAddresses;
426
427#ifdef __MARS__
428 typedef map<string, void*> Pointers;
429#else
430 typedef unordered_map<string, void*> Pointers;
431#endif
432 Pointers fPointers;
433
434 vector<vector<char>> fGarbage;
435
436 vector<char> fBufferRow;
437 vector<char> fBufferDat;
438
439 size_t fRow;
440
441 Checksum fChkHeader;
442 Checksum fChkData;
443
444 bool ReadBlock(vector<string> &vec)
445 {
446 int endtag = 0;
447 for (int i=0; i<36; i++)
448 {
449 char c[81];
450 c[80] = 0;
451 read(c, 80);
452 if (!good())
453 break;
454
455 fChkHeader.add(c, 80);
456
457// if (c[0]==0)
458// return vector<string>();
459
460 string str(c);
461
462// if (!str.empty())
463// cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl;
464
465 if (endtag==2 || str=="END ")
466 {
467 endtag = 2; // valid END tag found
468 continue;
469 }
470
471 if (endtag==1 || str==" ")
472 {
473 endtag = 1; // end tag not found, but expected to be there
474 continue;
475 }
476
477 vec.push_back(str);
478 }
479
480 // Make sure that no empty vector is returned
481 if (endtag && vec.size()%36==0)
482 vec.push_back(string("END = '' / "));
483
484 return endtag==2;
485 }
486
487 string Compile(const string &key, int16_t i=-1)
488 {
489 if (i<0)
490 return key;
491
492 ostringstream str;
493 str << key << i;
494 return str.str();
495 }
496
497public:
498 fits(const string &fname, bool force=false) : izstream(fname.c_str())
499 {
500 char simple[10];
501 read(simple, 10);
502 if (!good())
503 return;
504
505 if (memcmp(simple, "SIMPLE = ", 10))
506 {
507 clear(rdstate()|ios::badbit);
508#ifdef __EXCEPTIONS
509 throw runtime_error("File is not a FITS file.");
510#else
511 gLog << ___err___ << "ERROR - File is not a FITS file." << endl;
512 return;
513#endif
514 }
515
516 seekg(0);
517
518 while (good())
519 {
520 vector<string> block;
521 while (1)
522 {
523 // FIXME: Set limit on memory consumption
524 const int rc = ReadBlock(block);
525 if (!good())
526 {
527 clear(rdstate()|ios::badbit);
528#ifdef __EXCEPTIONS
529 throw runtime_error("FITS file corrupted.");
530#else
531 gLog << ___err___ << "ERROR - FITS file corrupted." << endl;
532 return;
533#endif
534 }
535
536 if (block.size()%36)
537 {
538 if (!rc && !force)
539 {
540 clear(rdstate()|ios::badbit);
541#ifdef __EXCEPTIONS
542 throw runtime_error("END keyword missing in FITS header.");
543#else
544 gLog << ___err___ << "ERROR - END keyword missing in FITS file... file might be corrupted." << endl;
545 return;
546#endif
547 }
548 break;
549 }
550 }
551
552 if (block.size()==0)
553 break;
554
555 if (block[0].substr(0, 9)=="SIMPLE =")
556 {
557 fChkHeader.reset();
558 continue;
559 }
560
561 if (block[0].substr(0, 9)=="XTENSION=")
562 {
563 // FIXME: Check for table name
564
565 fTable = Table(block, tellg());
566 fRow = (size_t)-1;
567
568 //fTable.PrintKeys();
569
570 if (!fTable)
571 {
572 clear(rdstate()|ios::badbit);
573 return;
574 }
575
576 fBufferRow.resize(fTable.bytes_per_row + 8-fTable.bytes_per_row%4);
577 fBufferDat.resize(fTable.bytes_per_row);
578
579 /*
580 // Next table should start at:
581 const size_t size = fTable.bytes_per_row*fTable.num_rows;
582 const size_t blks = size/(36*80);
583 const size_t rest = size%(36*80);
584
585 seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
586 if (!good())
587 gLog << ___err___ << "File seems to be incomplete (less data than expected from header)." << endl;
588
589 fRow = fTable.num_rows;
590 */
591
592 break;
593 }
594 }
595 }
596
597 uint8_t ReadRow(size_t row)
598 {
599 // if (row!=fRow+1) // Fast seeking is ensured by izstream
600 seekg(fTable.offset+row*fTable.bytes_per_row);
601
602 // For the checksum we need everything to be correctly aligned
603 const uint8_t offset = (row*fTable.bytes_per_row)%4;
604
605 auto ib = fBufferRow.begin();
606 auto ie = fBufferRow.end();
607 *ib++ = 0;
608 *ib++ = 0;
609 *ib++ = 0;
610 *ib = 0;
611
612 *--ie = 0;
613 *--ie = 0;
614 *--ie = 0;
615 *--ie = 0;
616 *--ie = 0;
617 *--ie = 0;
618 *--ie = 0;
619 *--ie = 0;
620
621 read(fBufferRow.data()+offset, fTable.bytes_per_row);
622 //fin.clear(fin.rdstate()&~ios::eofbit);
623
624 if (row==fRow+1)
625 fChkData.add(fBufferRow);
626
627 fRow = row;
628
629 return offset;
630 }
631
632 template<size_t N>
633 void revcpy(char *dest, const char *src, int num)
634 {
635 const char *pend = src + num*N;
636 for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
637 reverse_copy(ptr, ptr+N, dest);
638 }
639
640 bool GetRow(size_t row, bool check=true)
641 {
642 if (check && row>=fTable.num_rows)
643 return false;
644
645 const uint8_t offset = ReadRow(row);
646 if (!good())
647 return good();
648
649 const char *ptr = fBufferRow.data() + offset;
650
651 for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++)
652 {
653 const Table::Column &c = it->second;
654
655 const char *src = ptr + c.offset;
656 char *dest = reinterpret_cast<char*>(it->first);
657
658 // Let the compiler do some optimization by
659 // knowing the we only have 1, 2, 4 and 8
660 switch (c.size)
661 {
662 case 1: memcpy (dest, src, c.num*c.size); break;
663 case 2: revcpy<2>(dest, src, c.num); break;
664 case 4: revcpy<4>(dest, src, c.num); break;
665 case 8: revcpy<8>(dest, src, c.num); break;
666 }
667 }
668
669 return good();
670 }
671
672 bool GetNextRow(bool check=true)
673 {
674 return GetRow(fRow+1, check);
675 }
676
677 bool SkipNextRow()
678 {
679 seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
680 return good();
681 }
682
683 static bool Compare(const Address &p1, const Address &p2)
684 {
685 return p1.first>p2.first;
686 }
687
688 template<class T, class S>
689 const T &GetAs(const string &name)
690 {
691 return *reinterpret_cast<S*>(fPointers[name]);
692 }
693
694 void *SetPtrAddress(const string &name)
695 {
696 if (fTable.cols.count(name)==0)
697 {
698 ostringstream str;
699 str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
700#ifdef __EXCEPTIONS
701 throw runtime_error(str.str());
702#else
703 gLog << ___err___ << "ERROR - " << str.str() << endl;
704 return NULL;
705#endif
706 }
707
708 Pointers::const_iterator it = fPointers.find(name);
709 if (it!=fPointers.end())
710 return it->second;
711
712 fGarbage.push_back(vector<char>(fTable.cols[name].size*fTable.cols[name].num));
713
714 void *ptr = fGarbage.back().data();
715
716 fPointers[name] = ptr;
717 fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
718 sort(fAddresses.begin(), fAddresses.end(), Compare);
719 return ptr;
720 }
721
722 template<typename T>
723 bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
724 {
725 if (fTable.cols.count(name)==0)
726 {
727 ostringstream str;
728 str << "SetPtrAddress('" << name << "') - Column not found." << endl;
729#ifdef __EXCEPTIONS
730 throw runtime_error(str.str());
731#else
732 gLog << ___err___ << "ERROR - " << str.str() << endl;
733 return false;
734#endif
735 }
736
737 if (sizeof(T)!=fTable.cols[name].size)
738 {
739 ostringstream str;
740 str << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
741 << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
742#ifdef __EXCEPTIONS
743 throw runtime_error(str.str());
744#else
745 gLog << ___err___ << "ERROR - " << str.str() << endl;
746 return false;
747#endif
748 }
749
750 if (cnt!=fTable.cols[name].num)
751 {
752 ostringstream str;
753 str << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
754 << fTable.cols[name].num << " from header, got " << cnt << endl;
755#ifdef __EXCEPTIONS
756 throw runtime_error(str.str());
757#else
758 gLog << ___err___ << "ERROR - " << str.str() << endl;
759 return false;
760#endif
761 }
762
763 // if (fAddresses.count(ptr)>0)
764 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
765
766 //fAddresses[ptr] = fTable.cols[name];
767 fPointers[name] = ptr;
768 fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
769 sort(fAddresses.begin(), fAddresses.end(), Compare);
770 return true;
771 }
772
773 template<class T>
774 bool SetRefAddress(const string &name, T &ptr)
775 {
776 return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
777 }
778
779 template<typename T>
780 bool SetVecAddress(const string &name, vector<T> &vec)
781 {
782 return SetPtrAddress(name, vec.data(), vec.size());
783 }
784
785 template<typename T>
786 T Get(const string &key) const
787 {
788 return fTable.Get<T>(key);
789 }
790
791 template<typename T>
792 T Get(const string &key, const string &deflt) const
793 {
794 return fTable.Get<T>(key, deflt);
795 }
796
797 bool SetPtrAddress(const string &name, void *ptr)
798 {
799 if (fTable.cols.count(name)==0)
800 {
801 ostringstream str;
802 str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
803#ifdef __EXCEPTIONS
804 throw runtime_error(str.str());
805#else
806 gLog << ___err___ << "ERROR - " << str.str() << endl;
807 return false;
808#endif
809 }
810
811 // if (fAddresses.count(ptr)>0)
812 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
813
814 //fAddresses[ptr] = fTable.cols[name];
815 fPointers[name] = ptr;
816 fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
817 sort(fAddresses.begin(), fAddresses.end(), Compare);
818 return true;
819 }
820
821 bool HasKey(const string &key) const { return fTable.HasKey(key); }
822 bool HasColumn(const string& col) const { return fTable.HasColumn(col);}
823 const Table::Columns &GetColumns() const { return fTable.GetColumns();}
824 const Table::Keys &GetKeys() const { return fTable.GetKeys();}
825
826 int64_t GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
827 uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
828 double GetFloat(const string &key) const { return fTable.Get<double>(key); }
829 string GetStr(const string &key) const { return fTable.Get<string>(key); }
830
831 size_t GetN(const string &key) const
832 {
833 return fTable.GetN(key);
834 }
835
836 size_t GetNumRows() const { return fTable.num_rows; }
837 size_t GetRow() const { return fRow==(size_t)-1 ? 0 : fRow; }
838
839 operator bool() const { return fTable && fTable.offset!=0; }
840
841 void PrintKeys() const { fTable.PrintKeys(); }
842 void PrintColumns() const { fTable.PrintColumns(); }
843
844 bool IsHeaderOk() const { return fTable.datasum<0?false:(fChkHeader+Checksum(fTable.datasum)).valid(); }
845 bool IsFileOk() const { return (fChkHeader+fChkData).valid(); }
846
847};
848
849#ifndef __MARS__
850};
851#endif
852#endif
Note: See TracBrowser for help on using the repository browser.