source: trunk/Mars/mbase/MFits.h@ 11486

Last change on this file since 11486 was 11461, checked in by tbretz, 14 years ago
Made also Entry public -- for newer root versions and moved ReadBlock down.
File size: 15.2 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 map<void*, Table::Column> fAddresses;
297
298 vector<char> fBufferRow;
299 vector<char> fBufferDat;
300
301 size_t fRow;
302
303 vector<string> ReadBlock(vector<string> &vec)
304 {
305 bool end = false;
306 for (int i=0; i<36; i++)
307 {
308 char c[81];
309 c[80] = 0;
310 read(c, 80);
311 if (!good())
312 break;
313
314 if (c[0]==0)
315 return vector<string>();
316
317 string str(c);
318
319 //if (!str.empty())
320 // cout << setw(2) << i << "|" << str << "|" << endl;
321
322 if (str=="END ")
323 end = true;
324
325 if (end)
326 continue;
327
328 vec.push_back(str);
329 }
330
331 return vec;
332 }
333
334 string Compile(const string &key, int16_t i=-1)
335 {
336 if (i<0)
337 return key;
338
339 ostringstream str;
340 str << key << i;
341 return str.str();
342 }
343
344public:
345 MFits(const string &fname) : MZlib(fname.c_str())
346 {
347 while (good())
348 {
349 vector<string> block;
350 while (1)
351 {
352 // FIXME: Set limit on memory consumption
353 ReadBlock(block);
354 if (!good())
355 {
356 gLog << err << "FITS file corrupted." << endl;
357 clear(rdstate()|ios::badbit);
358 return;
359 }
360
361 if (block.size()%36)
362 break;
363 }
364
365 if (block.size()==0)
366 break;
367
368 if (block[0].substr(0, 9)=="SIMPLE =")
369 continue;
370
371 if (block[0].substr(0, 9)=="XTENSION=")
372 {
373 // FIXME: Check for table name
374
375 fTable = Table(block, tellg());
376 fRow = (size_t)-1;
377
378 fTable.PrintKeys();
379
380 if (!fTable)
381 return;
382
383 if (!fTable.Check("TELESCOP", 'T', "FACT"))
384 {
385 gLog << err << "Not a valid FACT fits file." << endl;
386 return;
387 }
388
389 fBufferRow.resize(fTable.bytes_per_row);
390 fBufferDat.resize(fTable.bytes_per_row);
391
392 /*
393 // Next table should start at:
394 const size_t size = fTable.bytes_per_row*fTable.num_rows;
395 const size_t blks = size/(36*80);
396 const size_t rest = size%(36*80);
397
398 seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
399 if (!good())
400 gLog << err << "File seems to be incomplete (less data than expected from header)." << endl;
401
402 fRow = fTable.num_rows;
403 */
404
405 break;
406 }
407 }
408 }
409
410 void ReadRow(size_t row)
411 {
412 // if (row!=fRow+1) // Fast seeking is ensured by MZlib
413 seekg(fTable.offset+row*fTable.bytes_per_row);
414
415 fRow = row;
416
417 read(fBufferRow.data(), fBufferRow.size());
418 //fin.clear(fin.rdstate()&~ios::eofbit);
419 }
420
421 bool GetRow(size_t row)
422 {
423 ReadRow(row);
424 if (!good())
425 return good();
426
427 for (map<void*,Table::Column>::const_iterator it=fAddresses.begin();
428 it!=fAddresses.end(); it++)
429 {
430 const Table::Column &c = it->second;
431
432 const char *src = fBufferRow.data() + c.offset;
433 const char *end = src + c.num*c.size;
434
435 char *dest = reinterpret_cast<char*>(it->first);
436 for (const char *ptr = src; ptr<end; ptr+=c.size, dest+=c.size)
437 reverse_copy(ptr, ptr+c.size, dest);
438 }
439
440 return good();
441 }
442
443 bool GetNextRow()
444 {
445 return GetRow(fRow+1);
446 }
447
448 bool SkipNextRow()
449 {
450 seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
451 return good();
452 }
453
454 template<typename T>
455 bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
456 {
457 if (fTable.cols.count(name)==0)
458 {
459 gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
460 return false;
461 }
462
463 if (sizeof(T)!=fTable.cols[name].size)
464 {
465 gLog << err << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
466 << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
467
468 return false;
469 }
470
471 if (cnt!=fTable.cols[name].num)
472 {
473 gLog << err << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
474 << fTable.cols[name].num << " from header, got " << cnt << endl;
475
476 return false;
477 }
478
479 // if (fAddresses.count(ptr)>0)
480 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
481
482 fAddresses[ptr] = fTable.cols[name];
483
484 return true;
485 }
486
487 template<class T>
488 bool SetRefAddress(const string &name, T &ptr)
489 {
490 return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
491 }
492
493 template<typename T>
494 bool SetVecAddress(const string &name, vector<T> &vec)
495 {
496 return SetPtrAddress(name, vec.data(), vec.size());
497 }
498
499 template<typename T>
500 T Get(const string &key) const
501 {
502 return fTable.Get<T>(key);
503 }
504
505 bool SetPtrAddress(const string &name, void *ptr)
506 {
507 if (fTable.cols.count(name)==0)
508 {
509 gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
510 return false;
511 }
512
513 // if (fAddresses.count(ptr)>0)
514 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
515
516 fAddresses[ptr] = fTable.cols[name];
517
518 return true;
519 }
520
521 int64_t GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
522 uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
523 double GetFloat(const string &key) const { return fTable.Get<double>(key); }
524
525 size_t GetN(const string &key) const
526 {
527 return fTable.GetN(key);
528 }
529
530 size_t GetNumRows() const { return fTable.num_rows; }
531 size_t GetRow() const { return fRow; }
532
533 operator bool() const { return fTable && fTable.offset!=0; }
534
535 void PrintKeys() const { fTable.PrintKeys(); }
536 void PrintColumns() const { fTable.PrintColumns(); }
537};
538#endif
Note: See TracBrowser for help on using the repository browser.