GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcsAZP.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcsAZP.cpp - Zenithal/azimuthal perspective (AZP) projection class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-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 GWcsAZP.cpp
23  * @brief Zenithal/azimuthal perspective (AZP) 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 "GWcsAZP.hpp"
34 #include "GWcsRegistry.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_PRJ_SET "GWcsAZP::prj_set()"
38 #define G_PRJ_X2S "GWcsAZP::prj_x2s(int, int, int, int, double*, double*,"\
39  " double*, double*, int*)"
40 #define G_PRJ_S2X "GWcsAZP::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_azp_registry(&g_wcs_azp_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 GWcsAZP::GWcsAZP(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 GWcsAZP::GWcsAZP(const GWcsAZP& 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 GWcsAZP::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 World Coordinate System.
205  ***************************************************************************/
207 {
208  return new GWcsAZP(*this);
209 }
210 
211 
212 
213 /***********************************************************************//**
214  * @brief Print WCS information
215  *
216  * @param[in] chatter Chattiness (defaults to NORMAL).
217  * @return String containing WCS information.
218  ***************************************************************************/
219 std::string GWcsAZP::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("=== GWcsAZP ===");
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::azpset.
283  *
284  * @exception GException::wcs_invalid_parameter
285  * PV(1) or PV(2) are invalid.
286  *
287  * Given:
288  * m_pv[1] Distance parameter, mu in units of r0.
289  * m_pv[2] Tilt angle, gamma in degrees.
290  *
291  * Given and/or returned:
292  * m_r0 Reset to 180/pi if 0.
293  * m_phi0 Reset to 0.0 if undefined.
294  * m_theta0 Reset to 90.0 if undefined.
295  *
296  * Returned:
297  * m_x0 Fiducial offset in x.
298  * m_y0 Fiducial offset in y.
299  * m_w[0] r0*(mu+1)
300  * m_w[1] tan(gamma)
301  * m_w[2] sec(gamma)
302  * m_w[3] cos(gamma)
303  * m_w[4] sin(gamma)
304  * m_w[5] asin(-1/mu) for |mu| >= 1, -90 otherwise
305  * m_w[6] mu*cos(gamma)
306  * m_w[7] 1 if |mu*cos(gamma)| < 1, 0 otherwise
307  ***************************************************************************/
308 void GWcsAZP::prj_set(void) const
309 {
310  // Signal that projection has been set (needs to be done before calling
311  // the prj_off() method to avoid an endless loop)
312  m_prjset = true;
313 
314  // Initialise projection parameters
315  m_w.assign(8, 0.0);
316 
317  // Set undefined parameters
318  if (undefined(m_pv[1])) m_pv[1] = 0.0;
319  if (undefined(m_pv[2])) m_pv[2] = 0.0;
320  if (m_r0 == 0.0) m_r0 = gammalib::rad2deg;
321 
322  // Precompute
323  m_w[0] = m_r0*(m_pv[1] + 1.0);
324  if (m_w[0] == 0.0) {
325  std::string message = "PV(1)=-1 encountered.";
327  }
328  m_w[3] = gammalib::cosd(m_pv[2]);
329  if (m_w[3] == 0.0) {
330  std::string message = "cos(PV(2))=0 encountered.";
332  }
333  m_w[2] = 1.0/m_w[3];
334  m_w[4] = gammalib::sind(m_pv[2]);
335  m_w[1] = m_w[4] / m_w[3];
336  if (std::abs(m_pv[1]) > 1.0) {
337  m_w[5] = gammalib::asind(-1.0/m_pv[1]);
338  }
339  else {
340  m_w[5] = -90.0;
341  }
342  m_w[6] = m_pv[1] * m_w[3];
343  m_w[7] = (std::abs(m_w[6]) < 1.0) ? 1.0 : 0.0;
344 
345  // Compute fiducial offset
346  prj_off(0.0, 90.0);
347 
348  // Return
349  return;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Pixel-to-spherical deprojection
355  *
356  * @param[in] nx X vector length.
357  * @param[in] ny Y vector length (0=no replication).
358  * @param[in] sxy Input vector step.
359  * @param[in] spt Output vector step.
360  * @param[in] x Vector of projected x coordinates.
361  * @param[in] y Vector of projected y coordinates.
362  * @param[out] phi Longitude of the projected point in native spherical
363  * coordinates [deg].
364  * @param[out] theta Latitude of the projected point in native spherical
365  * coordinates [deg].
366  * @param[out] stat Status return value for each vector element (always 0)
367  *
368  * @exception GException::wcs_invalid_x_y
369  * One or more of the (x,y) coordinates were invalid, as indicated
370  * by the stat vector.
371  *
372  * Deproject pixel (x,y) coordinates in the plane of projection to native
373  * spherical coordinates (phi,theta).
374  *
375  * This method has been adapted from the wcslib function prj.c::azpx2s().
376  * The interface follows very closely that of wcslib. In contrast to the
377  * wcslib routine, however, the method assumes that the projection has been
378  * setup previously (as this will be done by the constructor).
379  ***************************************************************************/
380 void GWcsAZP::prj_x2s(int nx, int ny, int sxy, int spt,
381  const double* x, const double* y,
382  double* phi, double* theta, int* stat) const
383 {
384  // Set tolerance
385  const double tol = 1.0e-13;
386 
387  // Initialize projection if required
388  if (!m_prjset) {
389  prj_set();
390  }
391 
392  // Set value replication length mx,my
393  int mx;
394  int my;
395  if (ny > 0) {
396  mx = nx;
397  my = ny;
398  }
399  else {
400  mx = 1;
401  my = 1;
402  ny = nx;
403  }
404 
405  // Initialise status code and statistics
406  int status = 0;
407  int n_invalid = 0;
408 
409  // Do x dependence
410  const double* xp = x;
411  int rowoff = 0;
412  int rowlen = nx * spt;
413  for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
414  double xj = *xp + m_x0;
415  double* phip = phi + rowoff;
416  for (int iy = 0; iy < my; ++iy, phip += rowlen) {
417  *phip = xj;
418  }
419  }
420 
421  // Do y dependence
422  const double* yp = y;
423  double* phip = phi;
424  double* thetap = theta;
425  int* statp = stat;
426  for (int iy = 0; iy < ny; ++iy, yp += sxy) {
427  double yj = *yp + m_y0;
428  double yc = yj * m_w[3];
429  double yc2 = yc*yc;
430  double q = m_w[0] + yj*m_w[4];
431  for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
432  double xj = *phip;
433  double r = sqrt(xj*xj + yc2);
434  if (r == 0.0) {
435  *phip = 0.0;
436  *thetap = 90.0;
437  *(statp++) = 0;
438 
439  }
440  else {
441  *phip = gammalib::atan2d(xj, -yc);
442  double s = r / q;
443  double t = s * m_pv[1]/std::sqrt(s*s + 1.0);
444  s = gammalib::atan2d(1.0, s);
445  if (std::abs(t) > 1.0) {
446  if (std::abs(t) > 1.0+tol) {
447  *thetap = 0.0;
448  *(statp++) = 1;
449  status = 3;
450  n_invalid++;
451  continue;
452  }
453  t = (t < 0) ? -90.0 : 90.0;
454  }
455  else {
456  t = gammalib::asind(t);
457  }
458  double a = s - t;
459  double b = s + t + 180.0;
460  if (a > 90.0) a -= 360.0;
461  if (b > 90.0) b -= 360.0;
462  *thetap = (a > b) ? a : b;
463  *(statp++) = 0;
464  }
465  }
466  }
467 
468  // Handle status code
469  if (status == 3) {
470  throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
471  }
472 
473  // Return
474  return;
475 }
476 
477 
478 /***********************************************************************//**
479  * @brief Generic spherical-to-pixel projection
480  *
481  * @param[in] nphi Longitude vector length.
482  * @param[in] ntheta Latitude vector length (0=no replication).
483  * @param[in] spt Input vector step.
484  * @param[in] sxy Output vector step.
485  * @param[in] phi Longitude vector of the projected point in native spherical
486  * coordinates [deg].
487  * @param[in] theta Latitude vector of the projected point in native spherical
488  * coordinates [deg].
489  * @param[out] x Vector of projected x coordinates.
490  * @param[out] y Vector of projected y coordinates.
491  * @param[out] stat Status return value for each vector element (always 0)
492  *
493  * @exception GException::wcs_invalid_phi_theta
494  * One or more of the (phi,theta) coordinates were invalid, as
495  * indicated by the stat vector.
496  *
497  * Project native spherical coordinates (phi,theta) to pixel (x,y)
498  * coordinates in the plane of projection.
499  *
500  * This method has been adapted from the wcslib function prj.c::azps2x().
501  * The interface follows very closely that of wcslib. In contrast to the
502  * wcslib routine, however, the method assumes that the projection has been
503  * setup previously (as this will be done by the constructor).
504  ***************************************************************************/
505 void GWcsAZP::prj_s2x(int nphi, int ntheta, int spt, int sxy,
506  const double* phi, const double* theta,
507  double* x, double* y, int* stat) const
508 {
509  // Initialize projection if required
510  if (!m_prjset) {
511  prj_set();
512  }
513 
514  // Set value replication length mphi,mtheta
515  int mphi;
516  int mtheta;
517  if (ntheta > 0) {
518  mphi = nphi;
519  mtheta = ntheta;
520  }
521  else {
522  mphi = 1;
523  mtheta = 1;
524  ntheta = nphi;
525  }
526 
527  // Initialise status code and statistics
528  int status = 0;
529  int n_invalid = 0;
530 
531  // Do phi dependence
532  const double* phip = phi;
533  int rowoff = 0;
534  int rowlen = nphi * sxy;
535  for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
536  double sinphi;
537  double cosphi;
538  gammalib::sincosd(*phip, &sinphi, &cosphi);
539  double* xp = x + rowoff;
540  double* yp = y + rowoff;
541  for (int itheta = 0; itheta < mtheta; ++itheta) {
542  *xp = sinphi;
543  *yp = cosphi;
544  xp += rowlen;
545  yp += rowlen;
546  }
547  }
548 
549  // Do theta dependence
550  const double* thetap = theta;
551  double* xp = x;
552  double* yp = y;
553  int* statp = stat;
554  for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
555 
556  // Compute sin(theta) and cos(theta)
557  double sinthe;
558  double costhe;
559  gammalib::sincosd(*thetap, &sinthe, &costhe);
560 
561  // ...
562  for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
563 
564  double s = m_w[1]*(*yp);
565  double t = (m_pv[1] + sinthe) + costhe*s;
566 
567  // ...
568  if (t == 0.0) {
569  *xp = 0.0;
570  *yp = 0.0;
571  *(statp++) = 1;
572  status = 4;
573  n_invalid++;
574  }
575 
576  // ...
577  else {
578  double r = m_w[0]*costhe/t;
579 
580  // Bounds checking
581  int istat = 0;
582  if (m_bounds) {
583 
584  // Is there overlap?
585  if (*thetap < m_w[5]) {
586  istat = 1;
587  status = 4;
588  n_invalid++;
589  }
590 
591  // Is there divergence?
592  else if (m_w[7] > 0.0) {
593  t = m_pv[1] / sqrt(1.0 + s*s);
594  if (std::abs(t) <= 1.0) {
595  s = gammalib::atand(-s);
596  t = gammalib::asind(t);
597  double a = s - t;
598  double b = s + t + 180.0;
599  if (a > 90.0) a -= 360.0;
600  if (b > 90.0) b -= 360.0;
601  if (*thetap < ((a > b) ? a : b)) {
602  istat = 1;
603  status = 4;
604  n_invalid++;
605  }
606  }
607  }
608  }
609 
610  // Set value
611  *xp = r*(*xp) - m_x0;
612  *yp = -r*(*yp) * m_w[2] - m_y0;
613  *(statp++) = istat;
614  } // endelse
615  } // endfor: phi
616  } // endfor: theta
617 
618  // Handle status code
619  if (status == 4) {
621  }
622 
623  // Return
624  return;
625 }
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
virtual ~GWcsAZP(void)
Destructor.
Definition: GWcsAZP.cpp:127
Interface definition for the WCS registry class.
Zenithal/azimuthal perspective (AZP) projection class definition.
Definition: GWcsAZP.hpp:43
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
#define G_PRJ_X2S
Definition: GWcsAZP.cpp:38
#define G_PRJ_S2X
Definition: GWcsAZP.cpp:40
GWcsAZP(void)
Void constructor.
Definition: GWcsAZP.cpp:67
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: GWcsAZP.cpp:505
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
bool undefined(const double &value) const
Check if value is undefined.
Definition: GWcs.hpp:262
const GWcsAZP g_wcs_azp_seed
Definition: GWcsAZP.cpp:54
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
World Coordinate Projection registry class interface definition.
double atand(const double &value)
Compute arc tangens in degrees.
Definition: GMath.cpp:282
#define G_PRJ_SET
Definition: GWcsAZP.cpp:37
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
virtual std::string print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcsAZP.cpp:219
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
GWcsAZP & operator=(const GWcsAZP &wcs)
Assignment operator.
Definition: GWcsAZP.cpp:149
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
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: GWcsAZP.cpp:380
double m_pv[PVN]
Projection parameters.
Definition: GWcs.hpp:213
double m_x0
Fiducial x offset.
Definition: GWcs.hpp:215
GChatter
Definition: GTypemaps.hpp:33
void prj_set(void) const
Setup of projection.
Definition: GWcsAZP.cpp:308
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcs.cpp:1602
void free_members(void)
Delete class members.
void free_members(void)
Delete class members.
Definition: GWcsAZP.cpp:271
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
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
double cosd(const double &angle)
Compute cosine of angle in degrees.
Definition: GMath.cpp:134
virtual void clear(void)
Clear instance.
Definition: GWcsAZP.cpp:184
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: GWcsAZP.cpp:249
const double rad2deg
Definition: GMath.hpp:44
virtual GWcsAZP * clone(void) const
Clone instance.
Definition: GWcsAZP.cpp:206
Zenithal/azimuthal perspective (AZP) projection class definition.
double sind(const double &angle)
Compute sine of angle in degrees.
Definition: GMath.cpp:163
void copy_members(const GWcsAZP &wcs)
Copy class members.
Definition: GWcsAZP.cpp:261
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
Mathematical function definitions.
bool m_bounds
Enable strict bounds checking.
Definition: GWcs.hpp:214