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