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

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