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

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