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