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

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