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

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