Context Navigation

source: trunk/Mars/mcore/ofits.h@ 17224

Visit:

Last change on this file since 17224 was 17224, checked in by tbretz, 11 years ago
Added missing include; setup a reasonable scheme for calling AddColumn for potentally compressed data
File size: 31.3 KB

Line
1	#ifndef MARS_ofits
2	#define MARS_ofits
3
4	#include <string>
5	#include <string.h>
6	#include <algorithm>
7	#include <sstream>
8	#include <iostream>
9	#include <fstream>
10	#include <iomanip>
11	#include <vector>
12	#include <algorithm>
13	#include <stdexcept>
14
15	#ifdef __CINT__
16	#define off_t size_t
17	#endif
18
19	#ifndef __MARS__
20	#define gLog cerr
21	#define ___err___ ""
22	#define ___warn___ ""
23	#define ___all___ ""
24	#else
25	#include "MLog.h"
26	#include "MLogManip.h"
27	#define ___err___ err
28	#define ___warn___ warn
29	#define ___all___ all
30	#endif
31
32	#include "FITS.h"
33	#include "checksum.h"
34
35	#ifndef __MARS__
36	namespace std
37	{
38	#else
39	using namespace std;
40	#endif
41
42	// Sloppy: allow / <--- left
43	// allow all characters (see specs for what is possible)
44
45	// units: m kg s rad sr K A mol cd Hz J W V N Pa C Ohm S F Wb T Hlm lx
46
47	class ofits : public ofstream
48	{
49	public:
50	struct Key
51	{
52	string key;
53	bool delim;
54	string value;
55	string comment;
56	string fitsString;
57
58	off_t offset; // offset in file
59
60	bool changed; // For closing the file
61
62	Key(const string &k="") : key(k), delim(false), fitsString(""), offset(0), changed(true) { }
63
64	string Trim(const string &str)
65	{
66	// Trim Both leading and trailing spaces
67	const size_t first = str.find_first_not_of(' '); // Find the first character position after excluding leading blank spaces
68	const size_t last = str.find_last_not_of(' '); // Find the first character position from reverse af
69
70	// if all spaces or empty return an empty string
71	if (string::npos==first \|\| string::npos==last)
72	return string();
73
74	return str.substr(first, last-first+1);
75	}
76
77	bool FormatKey()
78	{
79	key = Trim(key);
80	if (key.empty())
81	{
82	#ifdef __EXCEPTIONS
83	throw runtime_error("Key name empty.");
84	#else
85	gLog << ___err___ << "ERROR - Key name empty." << endl;
86	return false;
87	#endif
88	}
89	if (key.size()>8)
90	{
91	ostringstream sout;
92	sout << "Key '" << key << "' exceeds 8 bytes.";
93	#ifdef __EXCEPTIONS
94	throw runtime_error(sout.str());
95	#else
96	gLog << ___err___ << "ERROR - " << sout.str() << endl;
97	return false;
98	#endif
99	}
100
101	//transform(key.begin(), key.end(), key.begin(), toupper);
102
103	for (string::const_iterator c=key.begin(); c<key.end(); c++)
104	if ((c<'A' \|\| c>'Z') && (c<'0' \|\| c>'9') && c!='-' && c!='_')
105	{
106	ostringstream sout;
107	sout << "Invalid character '" << *c << "' found in key '" << key << "'";
108	#ifdef __EXCEPTIONS
109	throw runtime_error(sout.str());
110	#else
111	gLog << ___err___ << "ERROR - " << sout.str() << endl;
112	return false;
113	#endif
114	}
115
116	return true;
117	}
118
119	bool FormatComment()
120	{
121	comment = Trim(comment);
122
123	for (string::const_iterator c=key.begin(); c<key.end(); c++)
124	if (c<32 \|\| c>126)
125	{
126	ostringstream sout;
127	sout << "Invalid character '" << c << "' [" << int(c) << "] found in comment '" << comment << "'";
128	#ifdef __EXCEPTIONS
129	throw runtime_error(sout.str());
130	#else
131	gLog << ___err___ << "ERROR - " << sout.str() << endl;
132	return false;
133	#endif
134	}
135
136	return true;
137	}
138
139	bool check(bool trim=false)
140	{
141	if (!FormatKey())
142	return false;
143
144	if (!FormatComment())
145	return false;
146
147	size_t sz = CalcSize();
148	if (sz<=80)
149	return true;
150
151	if (!trim)
152	{
153	ostringstream sout;
154	sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters.";
155	#ifdef __EXCEPTIONS
156	throw runtime_error(sout.str());
157	#else
158	gLog << ___err___ << "ERROR - " << sout.str() << endl;
159	#endif
160	return false;
161	}
162
163	//looks like something went wrong. Maybe entry is too long ?
164	//try to remove the comment
165	comment = "";
166
167	sz = CalcSize();
168	if (sz<=80)
169	{
170	#ifndef __EXCEPTIONS
171	ostringstream sout;
172	sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters... removed comment.";
173	gLog << ___warn___ << "WARNING - " << sout.str() << endl;
174	#endif
175	return true;
176	}
177
178	ostringstream sout;
179	sout << "Size " << sz << " of entry for key '" << key << "' exceeds 80 characters even without comment.";
180	#ifdef __EXCEPTIONS
181	throw runtime_error(sout.str());
182	#else
183	gLog << ___err___ << "ERROR - " << sout.str() << endl;
184	return false;
185	#endif
186	}
187
188	size_t CalcSize() const
189	{
190	if (!delim)
191	return 10+comment.size();
192
193	return 10 + (value.size()<20?20:value.size()) + 3 + comment.size();
194	}
195
196	string Compile()
197	{
198
199	if (fitsString != "")
200	return fitsString;
201
202	ostringstream sout;
203	sout << std::left << setw(8) << key;
204
205	if (!delim)
206	{
207	sout << " " << comment;
208	return sout.str();
209	}
210
211	sout << "= ";
212	sout << (!value.empty() && value[0]=='\''?std::left:std::right);
213	sout << setw(20) << value << std::left;
214
215	if (!comment.empty())
216	sout << " / " << comment;
217
218	return sout.str();
219	}
220
221	Checksum checksum;
222
223	void Out(ofstream &fout)
224	{
225	if (!changed)
226	return;
227
228	string str = Compile();
229	str.insert(str.end(), 80-str.size(), ' ');
230
231	if (offset==0)
232	offset = fout.tellp();
233
234	//cout << "Write[" << offset << "]: " << key << "/" << value << endl;
235
236	fout.seekp(offset);
237	fout << str;
238
239	checksum.reset();
240	checksum.add(str.c_str(), 80);
241
242	changed = false;
243	}
244	/*
245	void Out(ostream &out)
246	{
247	string str = Compile();
248
249	str.insert(str.end(), 80-str.size(), ' ');
250
251	out << str;
252	changed = false;
253	}*/
254	};
255
256	private:
257	vector<Key> fKeys;
258
259	vector<Key>::iterator findkey(const string &key)
260	{
261	for (auto it=fKeys.begin(); it!=fKeys.end(); it++)
262	if (key==it->key)
263	return it;
264
265	return fKeys.end();
266	}
267
268	bool Set(const string &key="", bool delim=false, const string &value="", const string &comment="")
269	{
270	// If no delimit add the row no matter if it alread exists
271	if (delim)
272	{
273	// if the row already exists: update it
274	auto it = findkey(key);
275	if (it!=fKeys.end())
276	{
277	it->value = value;
278	it->changed = true;
279	return true;
280	}
281	}
282
283	if (fTable.num_rows>0)
284	{
285	ostringstream sout;
286	sout << "No new header key can be defined, rows were already written to the file... ignoring new key '" << key << "'";
287	#ifdef __EXCEPTIONS
288	throw runtime_error(sout.str());
289	#else
290	gLog << ___err___ << "ERROR - " << sout.str() << endl;
291	return false;
292	#endif
293	}
294
295	Key entry;
296
297	entry.key = key;
298	entry.delim = delim;
299	entry.value = value;
300	entry.comment = comment;
301	entry.offset = 0;
302	entry.changed = true;
303
304	if (!entry.check(fCommentTrimming))
305	return false;
306
307	fKeys.push_back(entry);
308	return true;
309	}
310
311	protected:
312	struct Table
313	{
314	off_t offset;
315
316	size_t bytes_per_row;
317	size_t num_rows;
318	size_t num_cols;
319
320	struct Column
321	{
322	string name;
323	size_t offset;
324	size_t num;
325	size_t size;
326	char type;
327	};
328
329	vector<Column> cols;
330
331	Table() : offset(0), bytes_per_row(0), num_rows(0), num_cols(0)
332	{
333	}
334	};
335
336
337	Table fTable;
338
339	vector<char> fOutputBuffer;
340
341	vector<Table::Column>::iterator findcol(const string &name)
342	{
343	for (auto it=fTable.cols.begin(); it!=fTable.cols.end(); it++)
344	if (name==it->name)
345	return it;
346
347	return fTable.cols.end();
348	}
349
350	Checksum fDataSum;
351	Checksum fHeaderSum;
352
353	bool fCommentTrimming;
354	bool fManualExtName;
355
356	public:
357	ofits() : fCommentTrimming(false),
358	fManualExtName(false)
359	{
360	}
361	ofits(const char *fname) : ofstream(),
362	fCommentTrimming(false),
363	fManualExtName(false)
364	{
365	this->open(fname);
366	}
367	virtual ~ofits() { if (is_open()) close(); }
368
369	virtual void open(const char * filename, bool addEXTNAMEKey=true)
370	{
371	fDataSum = 0;
372	fHeaderSum = 0;
373
374	fTable = Table();
375	fKeys.clear();
376
377	SetStr("XTENSION", "BINTABLE", "binary table extension");
378	SetInt("BITPIX", 8, "8-bit bytes");
379	SetInt("NAXIS", 2, "2-dimensional binary table");
380	SetInt("NAXIS1", 0, "width of table in bytes");
381	SetInt("NAXIS2", 0, "number of rows in table");
382	SetInt("PCOUNT", 0, "size of special data area");
383	SetInt("GCOUNT", 1, "one data group (required keyword)");
384	SetInt("TFIELDS", 0, "number of fields in each row");
385	if (addEXTNAMEKey)
386	SetStr("EXTNAME", "", "name of extension table");
387	else
388	fManualExtName = true;
389	SetStr("CHECKSUM", "0000000000000000", "Checksum for the whole HDU");
390	SetStr("DATASUM", " 0", "Checksum for the data block");
391
392	ofstream::open(filename);
393	}
394
395	void AllowCommentsTrimming(bool allow)
396	{
397	fCommentTrimming = allow;
398	}
399	//Etienne: required to enable 1 to 1 reconstruction of files
400	bool SetKeyComment(const string& key, const string& comment)
401	{
402	auto it = findkey(key);
403	if (it==fKeys.end())
404	return false;
405	it->comment = comment;
406	it->changed = true;
407	return true;
408	}
409	bool SetKeyFromFitsString(const string& fitsString)
410	{
411	if (fTable.num_rows>0)
412	{
413	ostringstream sout;
414	sout << "No new header key can be defined, rows were already written to the file... ignoring new key '" << fitsString << "'";
415	#ifdef __EXCEPTIONS
416	throw runtime_error(sout.str());
417	#else
418	gLog << ___err___ << "ERROR - " << sout.str() << endl;
419	return false;
420	#endif
421	}
422
423	Key entry;
424	entry.fitsString = fitsString;
425	entry.changed = true;
426	fKeys.push_back(entry);
427	return true;
428	}
429	bool SetRaw(const string &key, const string &val, const string &comment)
430	{
431	return Set(key, true, val, comment);
432	}
433
434	bool SetBool(const string &key, bool b, const string &comment="")
435	{
436	return Set(key, true, b?"T":"F", comment);
437	}
438
439	bool AddEmpty(const string &key, const string &comment="")
440	{
441	return Set(key, true, "", comment);
442	}
443
444	bool SetStr(const string &key, string s, const string &comment="")
445	{
446	for (uint i=0; i<s.length(); i++)
447	if (s[i]=='\'')
448	s.insert(i++, "\'");
449
450	return Set(key, true, "'"+s+"'", comment);
451	}
452
453	bool SetInt(const string &key, int64_t i, const string &comment="")
454	{
455	ostringstream sout;
456	sout << i;
457
458	return Set(key, true, sout.str(), comment);
459	}
460
461	bool SetFloat(const string &key, double f, int p, const string &comment="")
462	{
463	ostringstream sout;
464
465	if (p<0)
466	sout << setprecision(-p) << fixed;
467	if (p>0)
468	sout << setprecision(p);
469	if (p==0)
470	sout << setprecision(f>1e-100 && f<1e100 ? 15 : 14);
471
472	sout << f;
473
474	string str = sout.str();
475
476	replace(str.begin(), str.end(), 'e', 'E');
477
478	if (str.find_first_of('E')==string::npos && str.find_first_of('.')==string::npos)
479	str += ".";
480
481	return Set(key, true, str, comment);
482	}
483
484	bool SetFloat(const string &key, double f, const string &comment="")
485	{
486	return SetFloat(key, f, 0, comment);
487	}
488
489	bool SetHex(const string &key, uint64_t i, const string &comment="")
490	{
491	ostringstream sout;
492	sout << std::hex << "0x" << i;
493	return SetStr(key, sout.str(), comment);
494	}
495
496	bool AddComment(const string &comment)
497	{
498	return Set("COMMENT", false, "", comment);
499	}
500
501	bool AddHistory(const string &comment)
502	{
503	return Set("HISTORY", false, "", comment);
504	}
505
506	void End()
507	{
508	Set("END");
509	while (fKeys.size()%36!=0)
510	fKeys.emplace_back();
511	}
512
513	string CommentFromType(char type)
514	{
515	string comment;
516
517	switch (type)
518	{
519	case 'L': comment = "1-byte BOOL]"; break;
520	case 'A': comment = "1-byte CHAR]"; break;
521	case 'B': comment = "1-byte BOOL]"; break;
522	case 'I': comment = "2-byte INT]"; break;
523	case 'J': comment = "4-byte INT]"; break;
524	case 'K': comment = "8-byte INT]"; break;
525	case 'E': comment = "4-byte FLOAT]"; break;
526	case 'D': comment = "8-byte FLOAT]"; break;
527	case 'Q': comment = "var. Length]"; break;
528	}
529
530	return comment;
531	}
532
533	uint32_t SizeFromType(char type)
534	{
535	size_t size = 0;
536
537	switch (type)
538	{
539	case 'L': size = 1; break;
540	case 'A': size = 1; break;
541	case 'B': size = 1; break;
542	case 'I': size = 2; break;
543	case 'J': size = 4; break;
544	case 'K': size = 8; break;
545	case 'E': size = 4; break;
546	case 'D': size = 8; break;
547	case 'Q': size = 16; break;
548	}
549
550	return size;
551	}
552	//ETIENNE to be able to restore the file 1 to 1, I must restore the header keys myself
553	virtual bool AddColumn(uint32_t cnt, char typechar, const string &name, const string &unit, const string &comment="", bool addHeaderKeys=true)
554	{
555	if (tellp()<0)
556	{
557	ostringstream sout;
558	sout << "File not open... ignoring column '" << name << "'";
559	#ifdef __EXCEPTIONS
560	throw runtime_error(sout.str());
561	#else
562	gLog << ___err___ << "ERROR - " << sout.str() << endl;
563	return false;
564	#endif
565	}
566
567	if (tellp()>0)
568	{
569	ostringstream sout;
570	sout << "Header already written, no new column can be defined... ignoring column '" << name << "'";
571	#ifdef __EXCEPTIONS
572	throw runtime_error(sout.str());
573	#else
574	gLog << ___err___ << "ERROR - " << sout.str() << endl;
575	return false;
576	#endif
577	}
578
579	if (findcol(name)!=fTable.cols.end())
580	{
581	ostringstream sout;
582	sout << "A column with the name '" << name << "' already exists.";
583	#ifdef __EXCEPTIONS
584	throw runtime_error(sout.str());
585	#else
586	gLog << ___err___ << "ERROR - " << sout.str() << endl;
587	return false;
588	#endif
589	}
590
591	typechar = toupper(typechar);
592
593	static const string allow("LABIJKEDQ");
594	if (std::find(allow.begin(), allow.end(), typechar)==allow.end())
595	{
596	ostringstream sout;
597	sout << "Column type '" << typechar << "' not supported.";
598	#ifdef __EXCEPTIONS
599	throw runtime_error(sout.str());
600	#else
601	gLog << ___err___ << "ERROR - " << sout.str() << endl;
602	return false;
603	#endif
604	}
605
606	ostringstream type;
607	type << cnt;
608	if (typechar=='Q')
609	type << "QB";
610	else
611	type << typechar;
612
613	fTable.num_cols++;
614
615	ostringstream typekey, formkey, unitkey, unitcom, typecom;
616	typekey << "TTYPE" << fTable.num_cols;
617	formkey << "TFORM" << fTable.num_cols;
618	unitkey << "TUNIT" << fTable.num_cols;
619	unitcom << "unit of " << name;
620
621	typecom << "format of " << name << " [";
622
623	typecom << CommentFromType(typechar);
624
625	if (addHeaderKeys)
626	{
627	SetStr(formkey.str(), type.str(), typecom.str());
628	SetStr(typekey.str(), name, comment);
629	if (!unit.empty())
630	SetStr(unitkey.str(), unit, unitcom.str());
631	}
632	size_t size = SizeFromType(typechar);
633
634	Table::Column col;
635
636	col.name = name;
637	col.type = typechar;
638	col.num = cnt;
639	col.size = size;
640	col.offset = fTable.bytes_per_row;
641
642	fTable.cols.push_back(col);
643
644	fTable.bytes_per_row += size*cnt;
645
646	// Align to four bytes
647	fOutputBuffer.resize(fTable.bytes_per_row + 4-fTable.bytes_per_row%4);
648
649	return true;
650	}
651
652	virtual bool AddColumn(const FITS::Compression&, uint32_t cnt, char typechar, const string& name, const string& unit, const string& comment="", bool addHeaderKeys=true)
653	{
654	return AddColumn(cnt, typechar, name, unit, comment, addHeaderKeys);
655	}
656
657	bool AddColumnShort(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
658	{ return AddColumn(cnt, 'I', name, unit, comment); }
659	bool AddColumnInt(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
660	{ return AddColumn(cnt, 'J', name, unit, comment); }
661	bool AddColumnLong(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
662	{ return AddColumn(cnt, 'K', name, unit, comment); }
663	bool AddColumnFloat(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
664	{ return AddColumn(cnt, 'E', name, unit, comment); }
665	bool AddColumnDouble(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
666	{ return AddColumn(cnt, 'D', name, unit, comment); }
667	bool AddColumnChar(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
668	{ return AddColumn(cnt, 'A', name, unit, comment); }
669	bool AddColumnByte(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
670	{ return AddColumn(cnt, 'B', name, unit, comment); }
671	bool AddColumnBool(uint32_t cnt, const string &name, const string &unit="", const string &comment="")
672	{ return AddColumn(cnt, 'L', name, unit, comment); }
673
674	bool AddColumnShort(const string &name, const string &unit="", const string &comment="")
675	{ return AddColumn(1, 'I', name, unit, comment); }
676	bool AddColumnInt(const string &name, const string &unit="", const string &comment="")
677	{ return AddColumn(1, 'J', name, unit, comment); }
678	bool AddColumnLong(const string &name, const string &unit="", const string &comment="")
679	{ return AddColumn(1, 'K', name, unit, comment); }
680	bool AddColumnFloat(const string &name, const string &unit="", const string &comment="")
681	{ return AddColumn(1, 'E', name, unit, comment); }
682	bool AddColumnDouble(const string &name, const string &unit="", const string &comment="")
683	{ return AddColumn(1, 'D', name, unit, comment); }
684	bool AddColumnChar(const string &name, const string &unit="", const string &comment="")
685	{ return AddColumn(1, 'A', name, unit, comment); }
686	bool AddColumnByte(const string &name, const string &unit="", const string &comment="")
687	{ return AddColumn(1, 'B', name, unit, comment); }
688	bool AddColumnBool(const string &name, const string &unit="", const string &comment="")
689	{ return AddColumn(1, 'L', name, unit, comment); }
690
691	bool AddColumnShort(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
692	{ return AddColumn(comp, cnt, 'I', name, unit, comment); }
693	bool AddColumnInt(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
694	{ return AddColumn(comp, cnt, 'J', name, unit, comment); }
695	bool AddColumnLong(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
696	{ return AddColumn(comp, cnt, 'K', name, unit, comment); }
697	bool AddColumnFloat(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
698	{ return AddColumn(comp, cnt, 'E', name, unit, comment); }
699	bool AddColumnDouble(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
700	{ return AddColumn(comp, cnt, 'D', name, unit, comment); }
701	bool AddColumnChar(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
702	{ return AddColumn(comp, cnt, 'A', name, unit, comment); }
703	bool AddColumnByte(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
704	{ return AddColumn(comp, cnt, 'B', name, unit, comment); }
705	bool AddColumnBool(const FITS::Compression &comp, uint32_t cnt, const string &name, const string &unit="", const string &comment="")
706	{ return AddColumn(comp, cnt, 'L', name, unit, comment); }
707
708	bool AddColumnShort(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
709	{ return AddColumn(comp, 1, 'I', name, unit, comment); }
710	bool AddColumnInt(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
711	{ return AddColumn(comp, 1, 'J', name, unit, comment); }
712	bool AddColumnLong(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
713	{ return AddColumn(comp, 1, 'K', name, unit, comment); }
714	bool AddColumnFloat(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
715	{ return AddColumn(comp, 1, 'E', name, unit, comment); }
716	bool AddColumnDouble(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
717	{ return AddColumn(comp, 1, 'D', name, unit, comment); }
718	bool AddColumnChar(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
719	{ return AddColumn(comp, 1, 'A', name, unit, comment); }
720	bool AddColumnByte(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
721	{ return AddColumn(comp, 1, 'B', name, unit, comment); }
722	bool AddColumnBool(const FITS::Compression &comp, const string &name, const string &unit="", const string &comment="")
723	{ return AddColumn(comp, 1, 'L', name, unit, comment); }
724
725	/*
726	bool AddKey(string key, double d, const string &comment)
727	{
728	ostringstream out;
729	out << d;
730
731	string s = out.str();
732
733	replace(s.begin(), s.end(), "e", "E");
734
735	return AddKey(key, s, comment);
736	}*/
737
738
739	Checksum WriteHeader(ofstream &fout)
740	{
741	Checksum sum;
742	uint32_t count=0;
743	for (auto it=fKeys.begin(); it!=fKeys.end(); it++)
744	{
745	it->Out(fout);
746	sum += it->checksum;
747	count++;
748	}
749	fout.flush();
750
751	return sum;
752	}
753
754	Checksum WriteHeader()
755	{
756	return WriteHeader(*this);
757	}
758
759	void FlushHeader()
760	{
761	const off_t pos = tellp();
762	WriteHeader();
763	seekp(pos);
764	}
765
766	Checksum WriteFitsHeader()
767	{
768	ofits h;
769
770	h.SetBool("SIMPLE", true, "file does conform to FITS standard");
771	h.SetInt ("BITPIX", 8, "number of bits per data pixel");
772	h.SetInt ("NAXIS", 0, "number of data axes");
773	h.SetBool("EXTEND", true, "FITS dataset may contain extensions");
774	h.SetStr ("CHECKSUM","0000000000000000", "Checksum for the whole HDU");
775	h.SetStr ("DATASUM", " 0", "Checksum for the data block");
776	h.AddComment("FITS (Flexible Image Transport System) format is defined in 'Astronomy");
777	h.AddComment("and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H");
778	h.End();
779
780	const Checksum sum = h.WriteHeader(*this);
781
782	h.SetStr("CHECKSUM", sum.str());
783
784	const size_t offset = tellp();
785	h.WriteHeader(*this);
786	seekp(offset);
787
788	return sum;
789	}
790
791	virtual bool WriteDrsOffsetsTable ()
792	{
793	return true;
794	}
795
796	virtual bool WriteCatalog()
797	{
798	return true;
799	}
800
801	virtual bool WriteTableHeader(const char *name="DATA")
802	{
803	if (tellp()>0)
804	{
805	#ifdef __EXCEPTIONS
806	throw runtime_error("File not empty anymore.");
807	#else
808	gLog << ___err___ << "ERROR - File not empty anymore." << endl;
809	return false;
810	#endif
811	}
812
813	fHeaderSum = WriteFitsHeader();
814
815	WriteDrsOffsetsTable();
816
817	if (!fManualExtName)
818	SetStr("EXTNAME", name);
819	SetInt("NAXIS1", fTable.bytes_per_row);
820	SetInt("TFIELDS", fTable.cols.size());
821
822	End();
823
824	WriteHeader();
825
826	WriteCatalog();
827
828	return good();
829	}
830
831	template<size_t N>
832	void revcpy(char dest, const char src, int num)
833	{
834	const char pend = src + numN;
835	for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
836	reverse_copy(ptr, ptr+N, dest);
837	}
838
839	virtual uint32_t GetBytesPerRow() const { return fTable.bytes_per_row; }
840
841	virtual bool WriteRow(const void *ptr, size_t cnt, bool byte_swap=true)
842	{
843	// FIXME: Make sure that header was already written
844	// or write header now!
845	if (cnt!=fTable.bytes_per_row)
846	{
847	ostringstream sout;
848	sout << "WriteRow - Size " << cnt << " does not match expected size " << fTable.bytes_per_row;
849	#ifdef __EXCEPTIONS
850	throw runtime_error(sout.str());
851	#else
852	gLog << ___err___ << "ERROR - " << sout.str() << endl;
853	return false;
854	#endif
855	}
856
857	// For the checksum we need everything to be correctly aligned
858	const uint8_t offset = fTable.offset%4;
859
860	char *buffer = fOutputBuffer.data() + offset;
861
862	auto ib = fOutputBuffer.begin();
863	auto ie = fOutputBuffer.rbegin();
864	*ib++ = 0;
865	*ib++ = 0;
866	*ib++ = 0;
867	*ib = 0;
868
869	*ie++ = 0;
870	*ie++ = 0;
871	*ie++ = 0;
872	*ie = 0;
873
874	if (!byte_swap)
875	memcpy(buffer, ptr, cnt);
876	else
877	{
878	for (auto it=fTable.cols.begin(); it!=fTable.cols.end(); it++)
879	{
880	const char src = reinterpret_cast<const char>(ptr) + it->offset;
881	char *dest = buffer + it->offset;
882
883	// Let the compiler do some optimization by
884	// knowing the we only have 1, 2, 4 and 8
885	switch (it->size)
886	{
887	case 1: memcpy (dest, src, it->num*it->size); break;
888	case 2: revcpy<2>(dest, src, it->num); break;
889	case 4: revcpy<4>(dest, src, it->num); break;
890	case 8: revcpy<8>(dest, src, it->num); break;
891	}
892	}
893	}
894
895	write(buffer, cnt);
896	fDataSum.add(fOutputBuffer);
897
898	fTable.num_rows++;
899	fTable.offset += cnt;
900	return good();
901	}
902
903	template<typename N>
904	bool WriteRow(const vector<N> &vec)
905	{
906	return WriteRow(vec.data(), vec.size()*sizeof(N));
907	}
908
909	// Flushes the number of rows to the header on disk
910	virtual void FlushNumRows()
911	{
912	SetInt("NAXIS2", fTable.num_rows);
913	FlushHeader();
914	}
915
916	size_t GetNumRows() const { return fTable.num_rows; }
917
918	void AlignTo2880Bytes()
919	{
920	if (tellp()%(80*36)>0)
921	{
922	vector<char> filler(8036-tellp()%(8036));
923	write(filler.data(), filler.size());
924	}
925	}
926
927	Checksum UpdateHeaderChecksum()
928	{
929
930	ostringstream dataSumStr;
931	dataSumStr << fDataSum.val();
932	SetStr("DATASUM", dataSumStr.str());
933
934	const Checksum sum = WriteHeader();
935
936	//sum += headersum;
937
938	SetStr("CHECKSUM", (sum+fDataSum).str());
939
940	return WriteHeader();
941	}
942	virtual bool close()
943	{
944	if (tellp()<0)
945	return false;
946
947	AlignTo2880Bytes();
948
949	// We don't have to jump back to the end of the file
950	SetInt("NAXIS2", fTable.num_rows);
951
952
953	const Checksum chk = UpdateHeaderChecksum();
954
955	ofstream::close();
956
957	if ((chk+fDataSum).valid())
958	return true;
959
960	ostringstream sout;
961	sout << "Checksum (" << std::hex << chk.val() << ") invalid.";
962	#ifdef __EXCEPTIONS
963	throw runtime_error(sout.str());
964	#else
965	gLog << ___err___ << "ERROR - " << sout.str() << endl;
966	return false;
967	#endif
968	}
969
970	pair<string, int> GetChecksumData()
971	{
972	string datasum;
973	string checksum;
974	//cannot directly use the Get methods, because they are only in fits.h
975	for (vector<Key>::const_iterator it=fKeys.begin(); it!= fKeys.end(); it++)
976	{
977	if (it->key == "CHECKSUM") checksum = it->value;
978	if (it->key == "DATASUM") datasum = it->value;
979	}
980	if (checksum[0] == '\'')
981	checksum = checksum.substr(1,checksum.size()-2);
982	if (datasum[0] == '\'')
983	datasum = datasum.substr(1, datasum.size()-2);
984	return make_pair(checksum, atoi(datasum.c_str()));
985	}
986	};
987
988	#ifndef __MARS__
989	};
990	#endif
991
992	#if 0
993	#include "fits.h"
994
995	int main()
996	{
997	using namespace std;
998
999	ofits h2("delme.fits");
1000
1001	h2.SetInt("KEY1", 1, "comment 1");
1002	h2.AddColumnInt(2, "testcol1", "counts", "My comment");
1003	h2.AddColumnInt("testcol2", "counts", "My comment");
1004	//h2.AddColumnInt("testcol2", "counts", "My comment");
1005	h2.SetInt("KEY2", 2, "comment 2");
1006
1007	/*
1008	AddFloat("X0", 0.000123456, "number of fields in each row");
1009	AddFloat("X1", 0, "number of fields in each row");
1010	AddFloat("X2", 12345, "number of fields in each row");
1011	AddFloat("X3", 123456.67890, "number of fields in each row");
1012	AddFloat("X4", 1234567890123456789.12345678901234567890, "number of fields in each row");
1013	AddFloat("X5", 1234567890.1234567890e20, "number of fields in each row");
1014	AddFloat("X6", 1234567890.1234567890e-20, "number of fields in each row");
1015	AddFloat("XB", 1234567890.1234567890e-111, "number of fields in each row");
1016	AddFloat("X7", 1e-5, "number of fields in each row");
1017	AddFloat("X8", 1e-6, "number of fields in each row");
1018	//AddStr("12345678", "123456789012345678", "12345678901234567890123456789012345678901234567");
1019	*/
1020	// -
1021
1022	h2.WriteTableHeader("TABLE_NAME");
1023
1024	for (int i=0; i<10; i++)
1025	{
1026	int j[3] = { i+10, i10, i100 };
1027	h2.WriteRow(j, 3*sizeof(i));
1028	}
1029
1030	//h2.AddColumnInt("testcol2", "counts", "My comment");
1031	//h2.SetInt("KEY3", 2, "comment 2");
1032	h2.SetInt("KEY2", 2, "comment 2xxx");
1033	h2.SetInt("KEY1", 11);
1034
1035	h2.close();
1036
1037	cout << "---" << endl;
1038
1039	fits f("delme.fits");
1040	if (!f)
1041	throw runtime_error("xxx");
1042
1043	cout << "Header is valid: " << f.IsHeaderOk() << endl;
1044
1045	while (f.GetNextRow());
1046
1047	cout << "File is valid: " << f.IsFileOk() << endl;
1048
1049	cout << "---" << endl;
1050
1051	return 0;
1052	}
1053	#endif
1054
1055	#endif

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

Download in other formats: