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

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