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

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