source: trunk/Mars/mbase/fits.h@ 11559

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