GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcsARC.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcsARC.cpp - Zenithal/azimuthal equidistant (ARC) projection class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2018 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 GWcsARC.cpp
23  * @brief Zenithal/azimuthal equidistant (ARC) projection class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GMath.hpp"
33 #include "GTools.hpp"
34 #include "GWcsARC.hpp"
35 #include "GWcsRegistry.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 #define G_PRJ_X2S "GWcsARC::prj_x2s(int, int, int, int, double*, double*,"\
39  " double*, double*, int*)"
40 #define G_PRJ_S2X "GWcsARC::prj_s2x(int, int, int, int, double*, double*,"\
41  " double*, double*, int*)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 /* __ Local prototypes ___________________________________________________ */
50 
51 /* __ Constants __________________________________________________________ */
52 
53 /* __ Globals ____________________________________________________________ */
55 const GWcsRegistry g_wcs_arc_registry(&g_wcs_arc_seed);
56 
57 
58 /*==========================================================================
59  = =
60  = Constructors/destructors =
61  = =
62  ==========================================================================*/
63 
64 /***********************************************************************//**
65  * @brief Void constructor
66  ***************************************************************************/
68 {
69  // Initialise class members
70  init_members();
71 
72  // Return
73  return;
74 }
75 
76 
77 /***********************************************************************//**
78  * @brief Constructor
79  *
80  * @param[in] coords Coordinate system.
81  * @param[in] crval1 X value of reference pixel.
82  * @param[in] crval2 Y value of reference pixel.
83  * @param[in] crpix1 X index of reference pixel (starting from 1).
84  * @param[in] crpix2 Y index of reference pixel (starting from 1).
85  * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
86  * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
87  ***************************************************************************/
88 GWcsARC::GWcsARC(const std::string& coords,
89  const double& crval1, const double& crval2,
90  const double& crpix1, const double& crpix2,
91  const double& cdelt1, const double& cdelt2) :
92  GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
93 
94 {
95  // Initialise class members
96  init_members();
97 
98  // Setup WCS derived parameters
99  wcs_set();
100 
101  // Return
102  return;
103 }
104 
105 
106 /***********************************************************************//**
107  * @brief Copy constructor
108  *
109  * @param[in] wcs World Coordinate System.
110  ***************************************************************************/
111 GWcsARC::GWcsARC(const GWcsARC& wcs) : GWcs(wcs)
112 {
113  // Initialise class members for clean destruction
114  init_members();
115 
116  // Copy members
117  copy_members(wcs);
118 
119  // Return
120  return;
121 }
122 
123 
124 /***********************************************************************//**
125  * @brief Destructor
126  ***************************************************************************/
128 {
129  // Free members
130  free_members();
131 
132  // Return
133  return;
134 }
135 
136 
137 /*==========================================================================
138  = =
139  = Operators =
140  = =
141  ==========================================================================*/
142 
143 /***********************************************************************//**
144  * @brief Assignment operator
145  *
146  * @param[in] wcs World Coordinate System.
147  * @return World Coordinate System.
148  ***************************************************************************/
150 {
151  // Execute only if object is not identical
152  if (this != &wcs) {
153 
154  // Copy base class members
155  this->GWcs::operator=(wcs);
156 
157  // Free members
158  free_members();
159 
160  // Initialise private members for clean destruction
161  init_members();
162 
163  // Copy members
164  copy_members(wcs);
165 
166  } // endif: object was not identical
167 
168  // Return this object
169  return *this;
170 }
171 
172 
173 /*==========================================================================
174  = =
175  = Public methods =
176  = =
177  ==========================================================================*/
178 
179 /***********************************************************************//**
180  * @brief Clear instance
181  *
182  * This method properly resets the object to an initial state.
183  ***************************************************************************/
184 void GWcsARC::clear(void)
185 {
186  // Free class members (base and derived classes, derived class first)
187  free_members();
188  this->GWcs::free_members();
190 
191  // Initialise members
193  this->GWcs::init_members();
194  init_members();
195 
196  // Return
197  return;
198 }
199 
200 
201 /***********************************************************************//**
202  * @brief Clone instance
203  *
204  * @return Pointer to deep copy of Zenithal/azimuthal equidistant projection.
205  ***************************************************************************/
207 {
208  return new GWcsARC(*this);
209 }
210 
211 
212 
213 /***********************************************************************//**
214  * @brief Print WCS information
215  *
216  * @param[in] chatter Chattiness.
217  * @return String containing WCS information.
218  ***************************************************************************/
219 std::string GWcsARC::print(const GChatter& chatter) const
220 {
221  // Initialise result string
222  std::string result;
223 
224  // Continue only if chatter is not silent
225  if (chatter != SILENT) {
226 
227  // Append header
228  result.append("=== GWcsARC ===");
229 
230  // Append information
231  result.append(wcs_print(chatter));
232 
233  } // endif: chatter was not silent
234 
235  // Return result
236  return result;
237 }
238 
239 
240 /*==========================================================================
241  = =
242  = Private methods =
243  = =
244  ==========================================================================*/
245 
246 /***********************************************************************//**
247  * @brief Initialise class members
248  ***************************************************************************/
250 {
251  // Return
252  return;
253 }
254 
255 
256 /***********************************************************************//**
257  * @brief Copy class members
258  *
259  * @param[in] wcs World Coordinate System.
260  ***************************************************************************/
262 {
263  // Return
264  return;
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Delete class members
270  ***************************************************************************/
272 {
273  // Return
274  return;
275 }
276 
277 
278 /***********************************************************************//**
279  * @brief Setup of projection
280  *
281  * This method sets up the projection information. The method has been
282  * adapted from the wcslib function prj.c::arcset.
283  *
284  * Given and/or returned:
285  * m_r0 Reset to 180/pi if 0.
286  * m_phi0 Reset to 0.0 if undefined.
287  * m_theta0 Reset to 90.0 if undefined.
288  *
289  * Returned:
290  * m_x0 Fiducial offset in x.
291  * m_y0 Fiducial offset in y.
292  * m_w[0] r0*(pi/180)
293  * m_w[1] (180/pi)/r0
294  ***************************************************************************/
295 void GWcsARC::prj_set(void) const
296 {
297  // Signal that projection has been set (needs to be done before calling
298  // the prj_off() method to avoid an endless loop)
299  m_prjset = true;
300 
301  // Initialise projection parameters
302  m_w.assign(2, 0.0);
303 
304  // Set undefined parameters
305  if (m_r0 == 0.0) {
307  m_w[0] = 1.0;
308  m_w[1] = 1.0;
309  }
310  else {
311  m_w[0] = m_r0 * gammalib::deg2rad;
312  m_w[1] = 1.0/m_w[0];
313  }
314 
315  // Compute fiducial offset
316  prj_off(0.0, 0.0);
317 
318  // Return
319  return;
320 }
321 
322 
323 /***********************************************************************//**
324  * @brief Pixel-to-spherical deprojection
325  *
326  * @param[in] nx X vector length.
327  * @param[in] ny Y vector length (0=no replication).
328  * @param[in] sxy Input vector step.
329  * @param[in] spt Output vector step.
330  * @param[in] x Vector of projected x coordinates.
331  * @param[in] y Vector of projected y coordinates.
332  * @param[out] phi Longitude of the projected point in native spherical
333  * coordinates [deg].
334  * @param[out] theta Latitude of the projected point in native spherical
335  * coordinates [deg].
336  * @param[out] stat Status return value for each vector element (always 0)
337  *
338  * @exception GException::wcs_invalid_x_y
339  * One or more of the (x,y) coordinates were invalid, as indicated
340  * by the stat vector.
341  *
342  * Deproject pixel (x,y) coordinates in the plane of projection to native
343  * spherical coordinates (phi,theta).
344  *
345  * This method has been adapted from the wcslib function prj.c::aitx2s().
346  * The interface follows very closely that of wcslib. In contrast to the
347  * wcslib routine, however, the method assumes that the projection has been
348  * setup previously (as this will be done by the constructor).
349  ***************************************************************************/
350 void GWcsARC::prj_x2s(int nx, int ny, int sxy, int spt,
351  const double* x, const double* y,
352  double* phi, double* theta, int* stat) const
353 {
354  // Initialize projection if required
355  if (!m_prjset) {
356  prj_set();
357  }
358 
359  // Set value replication length mx,my
360  int mx;
361  int my;
362  if (ny > 0) {
363  mx = nx;
364  my = ny;
365  }
366  else {
367  mx = 1;
368  my = 1;
369  ny = nx;
370  }
371 
372  // Do x dependence
373  const double* xp = x;
374  int rowoff = 0;
375  int rowlen = nx * spt;
376  for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
377  double xj = *xp + m_x0;
378  double* phip = phi + rowoff;
379  for (int iy = 0; iy < my; ++iy, phip += rowlen) {
380  *phip = xj;
381  }
382  }
383 
384  // Do y dependence
385  const double* yp = y;
386  double* phip = phi;
387  double* thetap = theta;
388  int* statp = stat;
389  for (int iy = 0; iy < ny; ++iy, yp += sxy) {
390  double yj = *yp + m_y0;
391  double yj2 = yj * yj;
392  for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
393  double xj = *phip;
394  double r = std::sqrt(xj*xj + yj2);
395  if (r == 0.0) {
396  *phip = 0.0;
397  *thetap = 90.0;
398  }
399  else {
400  *phip = gammalib::atan2d(xj, -yj);
401  *thetap = 90.0 - r * m_w[1];
402  }
403  *(statp++) = 0;
404  }
405  }
406 
407  // Do boundary checking
408  int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat);
409 
410  // Handle status code
411  if (status != 0) {
412  std::string msg = "One or more of the (x, y) coordinates were invalid "
413  "for the "+name()+" ("+code()+") projection.";
415  }
416 
417  // Return
418  return;
419 }
420 
421 
422 /***********************************************************************//**
423  * @brief Generic spherical-to-pixel projection
424  *
425  * @param[in] nphi Longitude vector length.
426  * @param[in] ntheta Latitude vector length (0=no replication).
427  * @param[in] spt Input vector step.
428  * @param[in] sxy Output vector step.
429  * @param[in] phi Longitude vector of the projected point in native spherical
430  * coordinates [deg].
431  * @param[in] theta Latitude vector of the projected point in native spherical
432  * coordinates [deg].
433  * @param[out] x Vector of projected x coordinates.
434  * @param[out] y Vector of projected y coordinates.
435  * @param[out] stat Status return value for each vector element (always 0)
436  *
437  * @exception GException::wcs_invalid_phi_theta
438  * One or more of the (phi,theta) coordinates were invalid, as
439  * indicated by the stat vector.
440  *
441  * Project native spherical coordinates (phi,theta) to pixel (x,y)
442  * coordinates in the plane of projection.
443  *
444  * This method has been adapted from the wcslib function prj.c::aits2x().
445  * The interface follows very closely that of wcslib. In contrast to the
446  * wcslib routine, however, the method assumes that the projection has been
447  * setup previously (as this will be done by the constructor).
448  ***************************************************************************/
449 void GWcsARC::prj_s2x(int nphi, int ntheta, int spt, int sxy,
450  const double* phi, const double* theta,
451  double* x, double* y, int* stat) const
452 {
453  // Initialize projection if required
454  if (!m_prjset) {
455  prj_set();
456  }
457 
458  // Set value replication length mphi,mtheta
459  int mphi;
460  int mtheta;
461  if (ntheta > 0) {
462  mphi = nphi;
463  mtheta = ntheta;
464  }
465  else {
466  mphi = 1;
467  mtheta = 1;
468  ntheta = nphi;
469  }
470 
471  // Initialise status code and statistics
472  //int status = 0;
473  //int n_invalid = 0;
474 
475  // Do phi dependence
476  const double* phip = phi;
477  int rowoff = 0;
478  int rowlen = nphi * sxy;
479  for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
480  double sinphi;
481  double cosphi;
482  gammalib::sincosd(*phip, &sinphi, &cosphi);
483  double* xp = x + rowoff;
484  double* yp = y + rowoff;
485  for (int itheta = 0; itheta < mtheta; ++itheta) {
486  *xp = sinphi;
487  *yp = cosphi;
488  xp += rowlen;
489  yp += rowlen;
490  }
491  }
492 
493  // Do theta dependence
494  const double* thetap = theta;
495  double* xp = x;
496  double* yp = y;
497  int* statp = stat;
498  for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
499 
500  // Compute r
501  double r = m_w[0] * (90.0 - *thetap);
502 
503  // Do phi dependence
504  for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
505  *xp = r * (*xp) - m_x0;
506  *yp = -r * (*yp) - m_y0;
507  *(statp++) = 0;
508  } // endfor: phi
509 
510  } // endfor: theta
511 
512  // Handle status code
513  /*
514  if (status == 4) {
515  throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
516  }
517  */
518 
519  // Return
520  return;
521 }
GWcsARC(void)
Void constructor.
Definition: GWcsARC.cpp:67
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
Interface definition for the WCS registry class.
Zenithal/azimuthal equidistant (ARC) projection class definition.
Definition: GWcsARC.hpp:43
void prj_set(void) const
Setup of projection.
Definition: GWcsARC.cpp:295
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
virtual std::string name(void) const
Return projection name.
Definition: GWcsARC.hpp:115
GWcsARC & operator=(const GWcsARC &wcs)
Assignment operator.
Definition: GWcsARC.cpp:149
virtual std::string code(void) const
Return projection code.
Definition: GWcsARC.hpp:101
Zenithal/azimuthal equidistant (ARC) projection class definition.
virtual ~GWcsARC(void)
Destructor.
Definition: GWcsARC.cpp:127
Gammalib tools definition.
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
virtual void clear(void)
Clear instance.
Definition: GWcsARC.cpp:184
void prj_x2s(int nx, int ny, int sxy, int spt, const double *x, const double *y, double *phi, double *theta, int *stat) const
Pixel-to-spherical deprojection.
Definition: GWcsARC.cpp:350
World Coordinate Projection registry class interface definition.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
virtual GWcsARC * clone(void) const
Clone instance.
Definition: GWcsARC.cpp:206
void sincosd(const double &angle, double *s, double *c)
Compute sine and cosine of angle in degrees.
Definition: GMath.cpp:342
virtual GWcs & operator=(const GWcs &wcs)
Assignment operator.
Definition: GWcs.cpp:180
double m_x0
Fiducial x offset.
Definition: GWcs.hpp:215
GChatter
Definition: GTypemaps.hpp:33
virtual std::string print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcsARC.cpp:219
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcs.cpp:1602
void prj_s2x(int nphi, int ntheta, int spt, int sxy, const double *phi, const double *theta, double *x, double *y, int *stat) const
Generic spherical-to-pixel projection.
Definition: GWcsARC.cpp:449
void free_members(void)
Delete class members.
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
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
double m_y0
Fiducial y offset.
Definition: GWcs.hpp:216
double atan2d(const double &y, const double &x)
Compute arc tangens in degrees.
Definition: GMath.cpp:308
Exception handler interface definition.
bool m_prjset
Projection is set.
Definition: GWcs.hpp:211
void init_members(void)
Initialise class members.
void init_members(void)
Initialise class members.
Definition: GWcsARC.cpp:249
const GWcsARC g_wcs_arc_seed
Definition: GWcsARC.cpp:54
const double rad2deg
Definition: GMath.hpp:44
void free_members(void)
Delete class members.
Definition: GWcsARC.cpp:271
#define G_PRJ_X2S
Definition: GWcsARC.cpp:38
void copy_members(const GWcsARC &wcs)
Copy class members.
Definition: GWcsARC.cpp:261
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
Mathematical function definitions.