Context Navigation

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

Visit:

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

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

Download in other formats: