File: fitswriter.h

package info (click to toggle)
wsclean 2.8-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 2,196 kB
  • sloc: cpp: 34,504; ansic: 234; python: 174; makefile: 10
file content (278 lines) | stat: -rw-r--r-- 7,982 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
#ifndef FITSWRITER_H
#define FITSWRITER_H

#include <fitsio.h>

#include <string>
#include <vector>
#include <map>
#include <cmath>

#include "polarization.h"
#include "fitsiochecker.h"

class FitsWriter : public FitsIOChecker
{
public:
	enum DimensionType { 
		FrequencyDimension, 
		PolarizationDimension, 
		AntennaDimension, 
		TimeDimension,
		MatrixDimension
	};
	
	FitsWriter() :
		_width(0), _height(0),
		_phaseCentreRA(0.0), _phaseCentreDec(0.0), _pixelSizeX(0.0), _pixelSizeY(0.0),
		_phaseCentreDL(0.0), _phaseCentreDM(0.0),
		_frequency(0.0), _bandwidth(0.0),
		_dateObs(0.0),
		_hasBeam(false),
		_beamMajorAxisRad(0.0), _beamMinorAxisRad(0.0), _beamPositionAngle(0.0),
		_polarization(Polarization::StokesI),
		_unit(JanskyPerBeam),
		_isUV(false),
		_telescopeName(), _observer(), _objectName(),
		_origin("AO/WSImager"), _originComment("Imager written by Andre Offringa"),
		_multiFPtr(nullptr)
	{
	}
	
	explicit FitsWriter(const class FitsReader& reader) :
		_width(0), _height(0),
		_phaseCentreRA(0.0), _phaseCentreDec(0.0), _pixelSizeX(0.0), _pixelSizeY(0.0),
		_phaseCentreDL(0.0), _phaseCentreDM(0.0),
		_frequency(0.0), _bandwidth(0.0),
		_dateObs(0.0),
		_hasBeam(false),
		_beamMajorAxisRad(0.0), _beamMinorAxisRad(0.0), _beamPositionAngle(0.0),
		_polarization(Polarization::StokesI),
		_unit(JanskyPerBeam),
		_isUV(false),
		_telescopeName(), _observer(), _objectName(),
		_origin("AO/WSImager"), _originComment("Imager written by Andre Offringa"),
		_multiFPtr(nullptr)
	{
		SetMetadata(reader);
	}
	
	~FitsWriter()
	{
		if(_multiFPtr != nullptr)
			FinishMulti();
	}
	
	template<typename NumType> void Write(const std::string& filename, const NumType* image) const;
	
	void WriteMask(const std::string& filename, const bool* mask) const;
	
	void StartMulti(const std::string& filename);
	
	template<typename NumType>
	void AddToMulti(const NumType* image)
	{
		if(_multiFPtr == 0)
			throw std::runtime_error("AddToMulti() called before StartMulti()");
		writeImage(_multiFPtr, _multiFilename, image, _currentPixel.data());
		size_t index = 2;
		_currentPixel[index]++;
		while(index < _currentPixel.size()-1 && _currentPixel[index] > long(_extraDimensions[index-2].size))
		{
			_currentPixel[index] = 1;
			++index;
			_currentPixel[index]++;
		}
	}
	
	void FinishMulti();
	
	void SetBeamInfo(double widthRad)
	{
		SetBeamInfo(widthRad, widthRad, 0.0);
	}
	void SetBeamInfo(double majorAxisRad, double minorAxisRad, double positionAngleRad)
	{
		_hasBeam = true;
		_beamMajorAxisRad = majorAxisRad;
		_beamMinorAxisRad = minorAxisRad;
		_beamPositionAngle = positionAngleRad;
	}
	void SetNoBeamInfo()
	{
		_hasBeam = false;
		_beamMajorAxisRad = 0.0;
		_beamMinorAxisRad = 0.0;
		_beamPositionAngle = 0.0;
	}
	void SetImageDimensions(size_t width, size_t height)
	{
		_width = width;
		_height = height;
	}
	void SetImageDimensions(size_t width, size_t height, double pixelSizeX, double pixelSizeY)
	{
		_width = width;
		_height = height;
		_pixelSizeX = pixelSizeX;
		_pixelSizeY = pixelSizeY;
	}
	void SetImageDimensions(size_t width, size_t height, double phaseCentreRA, double phaseCentreDec, double pixelSizeX, double pixelSizeY)
	{
		_width = width;
		_height = height;
		_phaseCentreRA = phaseCentreRA;
		_phaseCentreDec = phaseCentreDec;
		_pixelSizeX = pixelSizeX;
		_pixelSizeY = pixelSizeY;
	}
	void SetFrequency(double frequency, double bandwidth)
	{
		_frequency = frequency;
		_bandwidth = bandwidth;
	}
	void SetDate(double dateObs)
	{
		_dateObs = dateObs;
	}
	void SetPolarization(PolarizationEnum polarization)
	{
		_polarization = polarization;
	}
	Unit GetUnit() const { return _unit; }
	void SetUnit(Unit unit)
	{
		_unit = unit;
	}
	void SetIsUV(bool isUV)
	{
		_isUV = isUV;
	}
	void SetTelescopeName(const std::string& telescopeName)
	{
		_telescopeName = telescopeName;
	}
	void SetObserver(const std::string& observer)
	{
		_observer = observer;
	}
	void SetObjectName(const std::string& objectName)
	{
		_objectName = objectName;
	}
	void SetOrigin(const std::string& origin, const std::string& comment)
	{
		_origin = origin;
		_originComment = comment;
	}
	void SetHistory(const std::vector<std::string>& history)
	{
		_history = history;
	}
	void AddHistory(const std::string& historyLine)
	{
		_history.push_back(historyLine);
	}

	void SetMetadata(const class FitsReader& reader);
	
	double RA() const { return _phaseCentreRA; }
	double Dec() const { return _phaseCentreDec; }
	double Frequency() const { return _frequency; }
	double Bandwidth() const { return _bandwidth; }
	double BeamSizeMajorAxis() const { return _beamMajorAxisRad; }
	double BeamSizeMinorAxis() const { return _beamMinorAxisRad; }
	double BeamPositionAngle() const { return _beamPositionAngle; }
	
	void SetExtraKeyword(const std::string& name, const std::string& value)
	{
		if(_extraStringKeywords.count(name) != 0)
			_extraStringKeywords.erase(name);
		_extraStringKeywords.insert(std::make_pair(name, value));
	}
	void SetExtraKeyword(const std::string& name, double value)
	{
		if(_extraNumKeywords.count(name) != 0)
			_extraNumKeywords.erase(name);
		_extraNumKeywords.insert(std::make_pair(name, value));
	}
	void RemoveExtraKeyword(const std::string& name)
	{
		if(_extraNumKeywords.count(name) != 0)
			_extraNumKeywords.erase(name);
		if(_extraStringKeywords.count(name) != 0)
			_extraStringKeywords.erase(name);
	}
	void SetExtraStringKeywords(const std::map<std::string, std::string>& keywords)
	{
		_extraStringKeywords = keywords;
	}
	void SetExtraNumKeywords(const std::map<std::string, double>& keywords)
	{
		_extraNumKeywords = keywords;
	}
	void SetPhaseCentreShift(double dl, double dm)
	{
		_phaseCentreDL = dl;
		_phaseCentreDM = dm;
	}
	size_t Width() const { return _width; }
	size_t Height() const { return _height; }
	double PhaseCentreDL() const { return _phaseCentreDL; }
	double PhaseCentreDM() const { return _phaseCentreDM; }
	
	void CopyDoubleKeywordIfExists(class FitsReader& reader, const char* keywordName);
	void CopyStringKeywordIfExists(class FitsReader& reader, const char* keywordName);
	
	static void MJDToHMS(double mjd, int& hour, int& minutes, int& seconds, int& deciSec);
	
	void AddExtraDimension(enum DimensionType type, size_t size)
	{
		_extraDimensions.emplace_back(Dimension{type, size});
	}
	void SetTimeDirectionStart(double time) { _timeDirectionStart = time; }
	void SetTimeDirectionInc(double dTime) { _timeDirectionInc = dTime; }
private:
	struct Dimension
	{
		DimensionType type;
		size_t size;
	};
	
	template<typename T>
	static T setNotFiniteToZero(T num)
	{
		return std::isfinite(num) ? num : 0.0;
	}
	std::size_t _width, _height;
	double _phaseCentreRA, _phaseCentreDec, _pixelSizeX, _pixelSizeY;
	double _phaseCentreDL, _phaseCentreDM;
	double _frequency, _bandwidth;
	double _dateObs;
	bool _hasBeam;
	double _beamMajorAxisRad, _beamMinorAxisRad, _beamPositionAngle;
	PolarizationEnum _polarization;
	Unit _unit;
	bool _isUV;
	std::string _telescopeName, _observer, _objectName;
	std::string _origin, _originComment;
	std::vector<std::string> _history;
	std::vector<Dimension> _extraDimensions;
	std::map<std::string, std::string> _extraStringKeywords;
	std::map<std::string, double> _extraNumKeywords;
	double _timeDirectionStart, _timeDirectionInc;
	
	void julianDateToYMD(double jd, int &year, int &month, int &day) const;
	void writeHeaders(fitsfile*& fptr, const std::string& filename) const;
	void writeHeaders(fitsfile*& fptr, const std::string& filename, const std::vector<Dimension>& extraDimensions) const;
	void writeImage(fitsfile* fptr, const std::string& filename, const double* image, long* currentPixel) const;
	void writeImage(fitsfile* fptr, const std::string& filename, const float* image, long* currentPixel) const;
	template<typename NumType>
	void writeImage(fitsfile* fptr, const std::string& filename, const NumType* image, long* currentPixel) const;
	
	std::string _multiFilename;
	fitsfile *_multiFPtr;
	std::vector<long> _currentPixel;
};

#endif