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

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