Context Navigation

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

Visit:

Last change on this file since 14625 was 14484, checked in by tbretz, 13 years ago
Make cint happy for newer root versions.
File size: 24.1 KB

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

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

Download in other formats: