GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcs.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcs.hpp - Abstract world coordinate system base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2019 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GWcs.hpp
23  * @brief Abstract world coordinate system base class definition
24  * @author Juergen Knoedlseder
25  */
26 
27 #ifndef GWCS_HPP
28 #define GWCS_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <vector>
32 #include <string>
33 #include <map>
34 #include <iostream>
35 #include "GSkyProjection.hpp"
36 #include "GLog.hpp"
37 #include "GFitsHDU.hpp"
38 #include "GSkyDir.hpp"
39 #include "GSkyPixel.hpp"
40 
41 #ifdef _OPENMP
42 #include <omp.h>
43 #endif
44 
45 
46 /***********************************************************************//**
47  * @class GWcs
48  *
49  * @brief Abstract world coordinate system base class
50  ***************************************************************************/
51 class GWcs : public GSkyProjection {
52 
53 public:
54  // Constructors and destructors
55  GWcs(void);
56  GWcs(const std::string& coords,
57  const double& crval1, const double& crval2,
58  const double& crpix1, const double& crpix2,
59  const double& cdelt1, const double& cdelt2);
60  explicit GWcs(const GFitsHDU& hdu);
61  GWcs(const GWcs& wcs);
62  virtual ~GWcs(void);
63 
64  // Operators
65  virtual GWcs& operator=(const GWcs& wcs);
66 
67  // Pure virtual base class methods
68  virtual void clear(void) = 0;
69  virtual GWcs* clone(void) const = 0;
70  virtual std::string classname(void) const = 0;
71  virtual std::string code(void) const = 0;
72  virtual std::string name(void) const = 0;
73  virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
74 
75  // Implemented virtual base class methods
76  virtual int size(void) const;
77  virtual void read(const GFitsHDU& hdu);
78  virtual void write(GFitsHDU& hdu) const;
79  virtual double solidangle(const GSkyPixel& pixel) const;
80  virtual GSkyDir pix2dir(const GSkyPixel& pixel) const;
81  virtual GSkyPixel dir2pix(const GSkyDir& dir) const;
82 
83  // Other methods
84  void set(const std::string& coords,
85  const double& crval1, const double& crval2,
86  const double& crpix1, const double& crpix2,
87  const double& cdelt1, const double& cdelt2);
88  double crval(const int& inx) const;
89  double crpix(const int& inx) const;
90  double cdelt(const int& inx) const;
91 
92 private:
93  // Static constants (set in GWcs.cpp)
94  static const int PVN = 32;
95  static const double UNDEFINED;
96 
97 protected:
98  // Protected methods
99  void init_members(void);
100  void copy_members(const GWcs& wcs);
101  void free_members(void);
102  void set_members(const std::string& coords,
103  const double& crval1, const double& crval2,
104  const double& crpix1, const double& crpix2,
105  const double& cdelt1, const double& cdelt2);
106  virtual bool compare(const GSkyProjection& proj) const;
107  bool undefined(const double& value) const;
108 
109  // Methods adapted from wcslib::wcs.c
110  void wcs_ini(int naxis);
111  void wcs_set(void) const;
112  void wcs_set_radesys(void) const;
113  void wcs_set_ctype(void) const;
114  void wcs_p2s(int ncoord, int nelem, const double* pixcrd, double* imgcrd,
115  double* phi, double* theta, double* world, int* stat) const;
116  void wcs_s2p(int ncoord, int nelem, const double* world,
117  double* phi, double* theta, double* imgcrd,
118  double* pixcrd, int* stat) const;
119  std::string wcs_print(const GChatter& chatter = NORMAL) const;
120  std::string wcs_print_value(const double& value) const;
121 
122  // Methods adapted from wcslib::cel.c
123  void cel_ini(void) const;
124  void cel_set(void) const;
125  void cel_x2s(int nx, int ny, int sxy, int sll,
126  const double* x, const double* y,
127  double* phi, double* theta,
128  double* lng, double* lat, int* stat) const;
129  void cel_s2x(int nlng, int nlat, int sll, int sxy,
130  const double* lng, const double* lat,
131  double* phi, double* theta,
132  double* x, double* y, int* stat) const;
133 
134  // Methods adapted from wcslib::sph.c
135  void sph_x2s(int nphi, int ntheta, int spt, int sll,
136  const double* phi, const double* theta,
137  double* lng, double* lat) const;
138  void sph_s2x(int nlng, int nlat, int sll, int spt,
139  const double* lng, const double* lat,
140  double* phi, double* theta) const;
141 
142  // Methods adapted from wcslib::spc.c
143  void spc_ini(void);
144 
145  // Methods adapted from wcslib::lin.c
146  void lin_ini(int naxis);
147  void lin_set(void) const;
148  void lin_p2x(int ncoord, int nelem, const double* pixcrd, double* imgcrd) const;
149  void lin_x2p(int ncoord, int nelem, const double* imgcrd, double* pixcrd) const;
150  void lin_matinv(const std::vector<double>& mat, std::vector<double>& inv) const;
151 
152  // Methods adapted from wcslib::prj.c
153  void prj_ini(void) const;
154  void prj_off(const double& phi0, const double& theta0) const;
155  int prj_bchk(const double& tol,
156  const int& nphi,
157  const int& ntheta,
158  const int& spt,
159  double* phi,
160  double* theta,
161  int* stat) const;
162  virtual void prj_set(void) const = 0;
163  virtual void prj_x2s(int nx, int ny, int sxy, int spt,
164  const double* x, const double* y,
165  double* phi, double* theta, int* stat) const = 0;
166  virtual void prj_s2x(int nphi, int ntheta, int spt, int sxy,
167  const double* phi, const double* theta,
168  double* x, double* y, int* stat) const = 0;
169 
170  // World Coordinate System parameters
171  mutable bool m_wcsset; //!< WCS information is set
172  int m_naxis; //!< Number of axes
173  std::vector<double> m_crval; //!< CRVALia keyvalues for each coord axis
174  std::vector<std::string> m_cunit; //!< CUNITia keyvalues for each coord axis
175  mutable std::vector<std::string> m_ctype; //!< CTYPEia keyvalues for each coord axis
176  mutable std::vector<std::string> m_ctype_c; //!< CTYPEia comments for each coord axis
177  mutable double m_lonpole; //!< LONPOLEa keyvalue
178  mutable double m_latpole; //!< LATPOLEa keyvalue
179  double m_restfrq; //!< RESTFRQa keyvalue
180  double m_restwav; //!< RESTWAVa keyvalue
181  mutable std::string m_radesys; //!< RADESYS keyvalue
182  mutable double m_equinox; //!< EQUINOX keyvalue
183  std::vector<double> m_cd; //!< CDi_ja linear transformation matrix
184  std::vector<double> m_crota; //!< CROTAia keyvalues for each coord axis
185  int m_lng; //!< Longitude axis
186  int m_lat; //!< Latitude axis
187  int m_spec; //!< Spectral axis
188 
189  // Linear transformation parameters
190  mutable bool m_linset; //!< Linear transformation is set
191  mutable bool m_unity; //!< Signals unity PC matrix
192  std::vector<double> m_crpix; //!< CRPIXja keyvalues for each pixel axis
193  std::vector<double> m_pc; //!< PCi_ja linear transformation matrix
194  std::vector<double> m_cdelt; //!< CDELTia keyvalues for each coord axis
195  mutable std::vector<double> m_piximg; //!< Pixel to image transformation matrix
196  mutable std::vector<double> m_imgpix; //!< Image to pixel transformation matrix
197 
198  // Celestial transformation parameters
199  mutable bool m_celset; //!< Celestial transformation is set
200  mutable bool m_offset; //!< Force (x,y) = (0,0) at (phi_0,theta_0)
201  mutable double m_phi0; //!< Native azimuth angle of fiducial point
202  mutable double m_theta0; //!< Native zenith angle of fiducial point
203  mutable double m_ref[4]; //!< Celestial coordinates of fiducial
204  // point and native coordinates of celestial
205  // pole
206  mutable double m_euler[5];//!< Euler angles and functions thereof
207  mutable int m_latpreq; //!< LATPOLEa requirement
208  mutable bool m_isolat; //!< True if |latitude| is preserved
209 
210  // Projection parameters
211  mutable bool m_prjset; //!< Projection is set
212  mutable double m_r0; //!< Radius of the generating sphere
213  mutable double m_pv[PVN]; //!< Projection parameters
214  mutable bool m_bounds; //!< Enable strict bounds checking
215  mutable double m_x0; //!< Fiducial x offset
216  mutable double m_y0; //!< Fiducial y offset
217  mutable std::vector<double> m_w; //!< Intermediate values
218 
219  // Spectral transformation parameters
220  //struct spcprm spc
221 
222  // Pre-computation cache
223  mutable bool m_has_pix2dir_cache; //!< Has valid pix2dir cache value
224  mutable bool m_has_dir2pix_cache; //!< Has valid dir2pix cache value
225  mutable GSkyDir m_last_pix2dir_dir; //!< Last sky direction for pix2dir
226  mutable GSkyDir m_last_dir2pix_dir; //!< Last sky direction for dir2pix
227  mutable GSkyPixel m_last_pix2dir_pix; //!< Last pixel for pix2dir
228  mutable GSkyPixel m_last_dir2pix_pix; //!< Last pixel for dir2pix
229 
230  // Thread locking
231  void init_lock(const int& lock_id=0) const;
232  void set_lock(const int& lock_id=0) const;
233  void unset_lock(const int& lock_id=0) const;
234 
235  #ifdef _OPENMP
236  mutable std::map<int, omp_lock_t> m_locks;
237  #endif
238 };
239 
240 
241 /***********************************************************************//**
242  * @brief Return dimension of projection
243  *
244  * @return Dimension of projection.
245  *
246  * Returns the dimension of the projection.
247  ***************************************************************************/
248 inline
249 int GWcs::size(void) const
250 {
251  return 2;
252 }
253 
254 
255 /***********************************************************************//**
256  * @brief Check if value is undefined
257  *
258  * @param[in] value Value to check
259  * @return True if @p value is undefined, false otherwise
260  ***************************************************************************/
261 inline
262 bool GWcs::undefined(const double& value) const
263 {
264  return (value==UNDEFINED);
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Initializes an OpenMP lock with a specific name
270  *
271  * @param[in] lock_id Identifier of the lock that should be initialized
272  ***************************************************************************/
273 inline
274 void GWcs::init_lock(const int& lock_id) const
275 {
276  #ifdef _OPENMP
277  if (m_locks.count(lock_id) == 0) {
278  m_locks[lock_id] = omp_lock_t();
279  omp_init_lock(&m_locks[lock_id]);
280  }
281  #endif
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Sets an OpenMP lock with a specific name
287  *
288  * @param[in] lock_id Identifier of the lock that should be set
289  ***************************************************************************/
290 inline
291 void GWcs::set_lock(const int& lock_id) const
292 {
293  #ifdef _OPENMP
294  omp_set_lock(&m_locks[lock_id]);
295  #endif
296 }
297 
298 
299 /***********************************************************************//**
300  * @brief Releases a previously set OpenMP lock
301  *
302  * @param[in] lock_id Identifier of the lock to be released
303  ***************************************************************************/
304 inline
305 void GWcs::unset_lock(const int& lock_id) const
306 {
307  #ifdef _OPENMP
308  omp_unset_lock(&m_locks[lock_id]);
309  #endif
310 }
311 
312 #endif /* GWCS_HPP */
bool m_linset
Linear transformation is set.
Definition: GWcs.hpp:190
void wcs_p2s(int ncoord, int nelem, const double *pixcrd, double *imgcrd, double *phi, double *theta, double *world, int *stat) const
Pixel-to-world transformation.
Definition: GWcs.cpp:1491
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
GSkyPixel m_last_pix2dir_pix
Last pixel for pix2dir.
Definition: GWcs.hpp:227
void cel_ini(void) const
Initialise celestial transformation parameters.
Definition: GWcs.cpp:1826
std::vector< double > m_imgpix
Image to pixel transformation matrix.
Definition: GWcs.hpp:196
void set_members(const std::string &coords, const double &crval1, const double &crval2, const double &crpix1, const double &crpix2, const double &cdelt1, const double &cdelt2)
Set World Coordinate System parameters.
Definition: GWcs.cpp:1077
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of sky map pixel.
Definition: GWcs.cpp:682
void prj_off(const double &phi0, const double &theta0) const
Compute fiducial offset to force (x,y) = (0,0) at (phi0,theta0)
Definition: GWcs.cpp:3154
int m_naxis
Number of axes.
Definition: GWcs.hpp:172
void prj_ini(void) const
Initialise projection parameters.
Definition: GWcs.cpp:3124
std::vector< double > m_crpix
CRPIXja keyvalues for each pixel axis.
Definition: GWcs.hpp:192
void set(const std::string &coords, const double &crval1, const double &crval2, const double &crpix1, const double &crpix2, const double &cdelt1, const double &cdelt2)
Set World Coordinate System parameters.
Definition: GWcs.cpp:849
Sky direction class interface definition.
std::vector< std::string > m_ctype
CTYPEia keyvalues for each coord axis.
Definition: GWcs.hpp:175
void lin_p2x(int ncoord, int nelem, const double *pixcrd, double *imgcrd) const
Pixel-to-world linear transformation.
Definition: GWcs.cpp:2832
std::vector< double > m_piximg
Pixel to image transformation matrix.
Definition: GWcs.hpp:195
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
virtual void prj_set(void) const =0
bool undefined(const double &value) const
Check if value is undefined.
Definition: GWcs.hpp:262
int m_lng
Longitude axis.
Definition: GWcs.hpp:185
bool m_celset
Celestial transformation is set.
Definition: GWcs.hpp:199
int m_latpreq
LATPOLEa requirement.
Definition: GWcs.hpp:207
void wcs_s2p(int ncoord, int nelem, const double *world, double *phi, double *theta, double *imgcrd, double *pixcrd, int *stat) const
World-to-pixel transformation.
Definition: GWcs.cpp:1557
int m_spec
Spectral axis.
Definition: GWcs.hpp:187
void cel_x2s(int nx, int ny, int sxy, int sll, const double *x, const double *y, double *phi, double *theta, double *lng, double *lat, int *stat) const
Pixel-to-world celestial transformation.
Definition: GWcs.cpp:2151
void spc_ini(void)
Initialise spectral transformation parameters.
Definition: GWcs.cpp:2694
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:877
double m_ref[4]
Celestial coordinates of fiducial.
Definition: GWcs.hpp:203
void unset_lock(const int &lock_id=0) const
Releases a previously set OpenMP lock.
Definition: GWcs.hpp:305
Abstract sky projection base class definition.
virtual double solidangle(const GSkyPixel &pixel) const
Returns solid angle of pixel in units of steradians.
Definition: GWcs.cpp:477
GSkyPixel m_last_dir2pix_pix
Last pixel for dir2pix.
Definition: GWcs.hpp:228
std::vector< double > m_cdelt
CDELTia keyvalues for each coord axis.
Definition: GWcs.hpp:194
void sph_s2x(int nlng, int nlat, int sll, int spt, const double *lng, const double *lat, double *phi, double *theta) const
Rotation in the pixel-to-world direction.
Definition: GWcs.cpp:2489
double m_equinox
EQUINOX keyvalue.
Definition: GWcs.hpp:182
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
void cel_set(void) const
Setup of celestial transformation.
Definition: GWcs.cpp:1881
virtual std::string code(void) const =0
void wcs_ini(int naxis)
Initialise World Coordinate System.
Definition: GWcs.cpp:1157
Abstract FITS extension base class definition.
std::vector< double > m_crval
CRVALia keyvalues for each coord axis.
Definition: GWcs.hpp:173
Sky map pixel class.
Definition: GSkyPixel.hpp:74
double m_latpole
LATPOLEa keyvalue.
Definition: GWcs.hpp:178
bool m_unity
Signals unity PC matrix.
Definition: GWcs.hpp:191
void init_lock(const int &lock_id=0) const
Initializes an OpenMP lock with a specific name.
Definition: GWcs.hpp:274
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:900
std::vector< double > m_crota
CROTAia keyvalues for each coord axis.
Definition: GWcs.hpp:184
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
virtual void prj_s2x(int nphi, int ntheta, int spt, int sxy, const double *phi, const double *theta, double *x, double *y, int *stat) const =0
double m_lonpole
LONPOLEa keyvalue.
Definition: GWcs.hpp:177
double m_restfrq
RESTFRQa keyvalue.
Definition: GWcs.hpp:179
std::vector< std::string > m_ctype_c
CTYPEia comments for each coord axis.
Definition: GWcs.hpp:176
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:923
virtual GWcs & operator=(const GWcs &wcs)
Assignment operator.
Definition: GWcs.cpp:180
void wcs_set_ctype(void) const
Set CTYPEa keywords.
Definition: GWcs.cpp:1378
virtual void write(GFitsHDU &hdu) const
Write WCS definition into FITS HDU header.
Definition: GWcs.cpp:332
virtual GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel of sky direction.
Definition: GWcs.cpp:766
std::vector< double > m_cd
CDi_ja linear transformation matrix.
Definition: GWcs.hpp:183
virtual void read(const GFitsHDU &hdu)
Read WCS definition from FITS header.
Definition: GWcs.cpp:220
double m_pv[PVN]
Projection parameters.
Definition: GWcs.hpp:213
void sph_x2s(int nphi, int ntheta, int spt, int sll, const double *phi, const double *theta, double *lng, double *lat) const
Rotation in the pixel-to-world direction.
Definition: GWcs.cpp:2253
virtual bool compare(const GSkyProjection &proj) const
Compares sky projection.
Definition: GWcs.cpp:1116
double m_x0
Fiducial x offset.
Definition: GWcs.hpp:215
GChatter
Definition: GTypemaps.hpp:33
GSkyDir m_last_pix2dir_dir
Last sky direction for pix2dir.
Definition: GWcs.hpp:225
double m_restwav
RESTWAVa keyvalue.
Definition: GWcs.hpp:180
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcs.cpp:1602
GSkyDir m_last_dir2pix_dir
Last sky direction for dir2pix.
Definition: GWcs.hpp:226
double m_theta0
Native zenith angle of fiducial point.
Definition: GWcs.hpp:202
void lin_ini(int naxis)
Initialise linear transformation parameters.
Definition: GWcs.cpp:2712
std::string wcs_print_value(const double &value) const
Helper function for value printing.
Definition: GWcs.cpp:1797
void wcs_set_radesys(void) const
Set radesys and equinox members.
Definition: GWcs.cpp:1315
bool m_isolat
True if |latitude| is preserved.
Definition: GWcs.hpp:208
virtual ~GWcs(void)
Destructor.
Definition: GWcs.cpp:158
double m_phi0
Native azimuth angle of fiducial point.
Definition: GWcs.hpp:201
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
bool m_has_pix2dir_cache
Has valid pix2dir cache value.
Definition: GWcs.hpp:223
bool m_offset
Force (x,y) = (0,0) at (phi_0,theta_0)
Definition: GWcs.hpp:200
std::string m_radesys
RADESYS keyvalue.
Definition: GWcs.hpp:181
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
int prj_bchk(const double &tol, const int &nphi, const int &ntheta, const int &spt, double *phi, double *theta, int *stat) const
Performs bounds checking on native spherical coordinates.
Definition: GWcs.cpp:3203
Information logger class definition.
GWcs(void)
Void constructor.
Definition: GWcs.cpp:80
double m_y0
Fiducial y offset.
Definition: GWcs.hpp:216
void copy_members(const GWcs &wcs)
Copy class members.
Definition: GWcs.cpp:979
void set_lock(const int &lock_id=0) const
Sets an OpenMP lock with a specific name.
Definition: GWcs.hpp:291
void lin_matinv(const std::vector< double > &mat, std::vector< double > &inv) const
Invert linear transformation matrix.
Definition: GWcs.cpp:2979
virtual GWcs * clone(void) const =0
Clones object.
void cel_s2x(int nlng, int nlat, int sll, int sxy, const double *lng, const double *lat, double *phi, double *theta, double *x, double *y, int *stat) const
World-to-pixel celestial transformation.
Definition: GWcs.cpp:2197
std::vector< double > m_pc
PCi_ja linear transformation matrix.
Definition: GWcs.hpp:193
bool m_prjset
Projection is set.
Definition: GWcs.hpp:211
virtual void clear(void)=0
Clear object.
Abstract sky projection base class.
void lin_set(void) const
Initialise linear transformation parameters.
Definition: GWcs.cpp:2762
static const double UNDEFINED
Definition: GWcs.hpp:95
int m_lat
Latitude axis.
Definition: GWcs.hpp:186
bool m_wcsset
WCS information is set.
Definition: GWcs.hpp:171
Sky map pixel class definition.
double m_euler[5]
Euler angles and functions thereof.
Definition: GWcs.hpp:206
static const int PVN
Definition: GWcs.hpp:94
std::vector< std::string > m_cunit
CUNITia keyvalues for each coord axis.
Definition: GWcs.hpp:174
Sky direction class.
Definition: GSkyDir.hpp:62
virtual void prj_x2s(int nx, int ny, int sxy, int spt, const double *x, const double *y, double *phi, double *theta, int *stat) const =0
virtual std::string name(void) const =0
virtual int size(void) const
Return dimension of projection.
Definition: GWcs.hpp:249
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
void lin_x2p(int ncoord, int nelem, const double *imgcrd, double *pixcrd) const
World-to-pixel linear transformation.
Definition: GWcs.cpp:2908
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
bool m_bounds
Enable strict bounds checking.
Definition: GWcs.hpp:214
bool m_has_dir2pix_cache
Has valid dir2pix cache value.
Definition: GWcs.hpp:224
virtual std::string classname(void) const =0
Return class name.