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

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