Context Navigation

source: trunk/Mars/mcore/fits.h@ 12551

Visit:

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: