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

Last change on this file since 11491 was 11487, checked in by tbretz, 15 years ago
Optimized reverse copy procedure; added HasKey
File size: 16.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 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 while (good())
351 {
352 vector<string> block;
353 while (1)
354 {
355 // FIXME: Set limit on memory consumption
356 ReadBlock(block);
357 if (!good())
358 {
359 gLog << err << "FITS file corrupted." << endl;
360 clear(rdstate()|ios::badbit);
361 return;
362 }
363
364 if (block.size()%36)
365 break;
366 }
367
368 if (block.size()==0)
369 break;
370
371 if (block[0].substr(0, 9)=="SIMPLE =")
372 continue;
373
374 if (block[0].substr(0, 9)=="XTENSION=")
375 {
376 // FIXME: Check for table name
377
378 fTable = Table(block, tellg());
379 fRow = (size_t)-1;
380
381 fTable.PrintKeys();
382
383 if (!fTable)
384 return;
385
386 if (!fTable.Check("TELESCOP", 'T', "FACT"))
387 {
388 gLog << err << "Not a valid FACT fits file." << endl;
389 return;
390 }
391
392 fBufferRow.resize(fTable.bytes_per_row);
393 fBufferDat.resize(fTable.bytes_per_row);
394
395 /*
396 // Next table should start at:
397 const size_t size = fTable.bytes_per_row*fTable.num_rows;
398 const size_t blks = size/(36*80);
399 const size_t rest = size%(36*80);
400
401 seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
402 if (!good())
403 gLog << err << "File seems to be incomplete (less data than expected from header)." << endl;
404
405 fRow = fTable.num_rows;
406 */
407
408 break;
409 }
410 }
411 }
412
413 void ReadRow(size_t row)
414 {
415 // if (row!=fRow+1) // Fast seeking is ensured by MZlib
416 seekg(fTable.offset+row*fTable.bytes_per_row);
417
418 fRow = row;
419
420 read(fBufferRow.data(), fBufferRow.size());
421 //fin.clear(fin.rdstate()&~ios::eofbit);
422 }
423
424 template<size_t N>
425 void revcpy(char *dest, const char *src, int num)
426 {
427 const char *end = src + num*N;
428 for (const char *ptr = src; ptr<end; ptr+=N, dest+=N)
429 reverse_copy(ptr, ptr+N, dest);
430 }
431
432 bool GetRow(size_t row)
433 {
434 ReadRow(row);
435 if (!good())
436 return good();
437
438 for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++)
439 {
440 const Table::Column &c = it->second;
441
442 const char *src = fBufferRow.data() + c.offset;
443 char *dest = reinterpret_cast<char*>(it->first);
444
445 // Let the compiler do some optimization by
446 // knowing the we only have 1, 2, 4 and 8
447 switch (c.size)
448 {
449 case 1: memcpy (dest, src, c.num*c.size); break;
450 case 2: revcpy<2>(dest, src, c.num); break;
451 case 4: revcpy<4>(dest, src, c.num); break;
452 case 8: revcpy<8>(dest, src, c.num); break;
453 }
454 }
455
456 return good();
457 }
458
459 bool GetNextRow()
460 {
461 return GetRow(fRow+1);
462 }
463
464 bool SkipNextRow()
465 {
466 seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
467 return good();
468 }
469
470 static bool Compare(const Address &p1, const Address &p2)
471 {
472 return p1.first>p2.first;
473 }
474
475 template<typename T>
476 bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
477 {
478 if (fTable.cols.count(name)==0)
479 {
480 gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
481 return false;
482 }
483
484 if (sizeof(T)!=fTable.cols[name].size)
485 {
486 gLog << err << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
487 << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
488
489 return false;
490 }
491
492 if (cnt!=fTable.cols[name].num)
493 {
494 gLog << err << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
495 << fTable.cols[name].num << " from header, got " << cnt << endl;
496
497 return false;
498 }
499
500 // if (fAddresses.count(ptr)>0)
501 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
502
503 //fAddresses[ptr] = fTable.cols[name];
504 fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
505 sort(fAddresses.begin(), fAddresses.end(), Compare);
506 return true;
507 }
508
509 template<class T>
510 bool SetRefAddress(const string &name, T &ptr)
511 {
512 return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
513 }
514
515 template<typename T>
516 bool SetVecAddress(const string &name, vector<T> &vec)
517 {
518 return SetPtrAddress(name, vec.data(), vec.size());
519 }
520
521 template<typename T>
522 T Get(const string &key) const
523 {
524 return fTable.Get<T>(key);
525 }
526
527 bool SetPtrAddress(const string &name, void *ptr)
528 {
529 if (fTable.cols.count(name)==0)
530 {
531 gLog << err <<"SetPtrAddress('" << name << "') - Column not found." << endl;
532 return false;
533 }
534
535 // if (fAddresses.count(ptr)>0)
536 // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
537
538 //fAddresses[ptr] = fTable.cols[name];
539 fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
540 sort(fAddresses.begin(), fAddresses.end(), Compare);
541 return true;
542 }
543
544 bool HasKey(const string &key) const { return fTable.HasKey(key); }
545 int64_t GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
546 uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
547 double GetFloat(const string &key) const { return fTable.Get<double>(key); }
548
549 size_t GetN(const string &key) const
550 {
551 return fTable.GetN(key);
552 }
553
554 size_t GetNumRows() const { return fTable.num_rows; }
555 size_t GetRow() const { return fRow; }
556
557 operator bool() const { return fTable && fTable.offset!=0; }
558
559 void PrintKeys() const { fTable.PrintKeys(); }
560 void PrintColumns() const { fTable.PrintColumns(); }
561};
562#endif
Note: See TracBrowser for help on using the repository browser.