GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcsMOL.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcsMOL.cpp - Mollweide's projection class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-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 GWcsMOL.cpp
23  * @brief Mollweide's 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 "GWcsMOL.hpp"
35 #include "GWcsRegistry.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 #define G_PRJ_X2S "GWcsMOL::prj_x2s(int, int, int, int, double*, double*,"\
39  " double*, double*, int*)"
40 #define G_PRJ_S2X "GWcsMOL::prj_s2x(int, int, int, int, double*, double*,"\
41  " double*, double*, int*)"
42 
43 /* __ Macros _____________________________________________________________ */
44 #define copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
45 
46 /* __ Coding definitions _________________________________________________ */
47 
48 /* __ Debug definitions __________________________________________________ */
49 //#define G_DEBUG_PRJ //!< Debug GWcsMOL::prj_x2s
50 
51 /* __ Local prototypes ___________________________________________________ */
52 
53 /* __ Constants __________________________________________________________ */
54 
55 /* __ Globals ____________________________________________________________ */
57 const GWcsRegistry g_wcs_mol_registry(&g_wcs_mol_seed);
58 
59 
60 /*==========================================================================
61  = =
62  = Constructors/destructors =
63  = =
64  ==========================================================================*/
65 
66 /***********************************************************************//**
67  * @brief Void constructor
68  ***************************************************************************/
70 {
71  // Initialise class members
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief Constructor
81  *
82  * @param[in] coords Coordinate system.
83  * @param[in] crval1 X value of reference pixel.
84  * @param[in] crval2 Y value of reference pixel.
85  * @param[in] crpix1 X index of reference pixel (starting from 1).
86  * @param[in] crpix2 Y index of reference pixel (starting from 1).
87  * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
88  * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
89  ***************************************************************************/
90 GWcsMOL::GWcsMOL(const std::string& coords,
91  const double& crval1, const double& crval2,
92  const double& crpix1, const double& crpix2,
93  const double& cdelt1, const double& cdelt2) :
94  GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
95 
96 {
97  // Initialise class members
98  init_members();
99 
100  // Setup WCS derived parameters
101  wcs_set();
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Copy constructor
110  *
111  * @param[in] wcs World Coordinate System.
112  ***************************************************************************/
113 GWcsMOL::GWcsMOL(const GWcsMOL& wcs) : GWcs(wcs)
114 {
115  // Initialise class members for clean destruction
116  init_members();
117 
118  // Copy members
119  copy_members(wcs);
120 
121  // Return
122  return;
123 }
124 
125 
126 /***********************************************************************//**
127  * @brief Destructor
128  ***************************************************************************/
130 {
131  // Free members
132  free_members();
133 
134  // Return
135  return;
136 }
137 
138 
139 /*==========================================================================
140  = =
141  = Operators =
142  = =
143  ==========================================================================*/
144 
145 /***********************************************************************//**
146  * @brief Assignment operator
147  *
148  * @param[in] wcs World Coordinate System.
149  * @return World Coordinate System.
150  ***************************************************************************/
152 {
153  // Execute only if object is not identical
154  if (this != &wcs) {
155 
156  // Copy base class members
157  this->GWcs::operator=(wcs);
158 
159  // Free members
160  free_members();
161 
162  // Initialise private members for clean destruction
163  init_members();
164 
165  // Copy members
166  copy_members(wcs);
167 
168  } // endif: object was not identical
169 
170  // Return this object
171  return *this;
172 }
173 
174 
175 /*==========================================================================
176  = =
177  = Public methods =
178  = =
179  ==========================================================================*/
180 
181 /***********************************************************************//**
182  * @brief Clear instance
183  *
184  * This method properly resets the object to an initial state.
185  ***************************************************************************/
186 void GWcsMOL::clear(void)
187 {
188  // Free class members (base and derived classes, derived class first)
189  free_members();
190  this->GWcs::free_members();
192 
193  // Initialise members
195  this->GWcs::init_members();
196  init_members();
197 
198  // Return
199  return;
200 }
201 
202 
203 /***********************************************************************//**
204  * @brief Clone instance
205  *
206  * @return Pointer to deep copy of Mollweide projection.
207  ***************************************************************************/
209 {
210  return new GWcsMOL(*this);
211 }
212 
213 
214 
215 /***********************************************************************//**
216  * @brief Print WCS information
217  *
218  * @param[in] chatter Chattiness (defaults to NORMAL).
219  * @return String containing WCS information.
220  ***************************************************************************/
221 std::string GWcsMOL::print(const GChatter& chatter) const
222 {
223  // Initialise result string
224  std::string result;
225 
226  // Continue only if chatter is not silent
227  if (chatter != SILENT) {
228 
229  // Append header
230  result.append("=== GWcsMOL ===");
231 
232  // Append information
233  result.append(wcs_print(chatter));
234 
235  } // endif: chatter was not silent
236 
237  // Return result
238  return result;
239 }
240 
241 
242 /*==========================================================================
243  = =
244  = Private methods =
245  = =
246  ==========================================================================*/
247 
248 /***********************************************************************//**
249  * @brief Initialise class members
250  ***************************************************************************/
252 {
253  // Return
254  return;
255 }
256 
257 
258 /***********************************************************************//**
259  * @brief Copy class members
260  *
261  * @param[in] wcs World Coordinate System.
262  ***************************************************************************/
264 {
265  // Return
266  return;
267 }
268 
269 
270 /***********************************************************************//**
271  * @brief Delete class members
272  ***************************************************************************/
274 {
275  // Return
276  return;
277 }
278 
279 
280 /***********************************************************************//**
281  * @brief Setup of projection
282  *
283  * This method sets up the projection information. The method has been
284  * adapted from the wcslib function prj.c::molset (version 4.20).
285  *
286  * Given and/or returned:
287  * m_r0 Reset to 180/pi if 0.
288  * m_phi0 Reset to 0.0 if undefined.
289  * m_theta0 Reset to 0.0 if undefined.
290  *
291  * Returned:
292  * m_x0 Fiducial offset in x.
293  * m_y0 Fiducial offset in y.
294  * m_w[0] sqrt(2)*r0
295  * m_w[1] sqrt(2)*r0/90
296  * m_w[2] 1/(sqrt(2)*r0)
297  * m_w[3] 90/r0
298  * m_w[4] 2/pi
299  ***************************************************************************/
300 void GWcsMOL::prj_set(void) const
301 {
302  // Signal that projection has been set (needs to be done before calling
303  // the prj_off() method to avoid an endless loop)
304  m_prjset = true;
305 
306  // Initialise projection parameters
307  m_w.assign(5, 0.0);
308 
309  // Set undefined parameters
310  if (m_r0 == 0.0) {
312  }
313 
314  // Precompute
315  m_w[0] = gammalib::sqrt_two * m_r0;
316  m_w[1] = m_w[0]/90.0;
317  m_w[2] = 1.0/m_w[0];
318  m_w[3] = 90.0/m_r0;
319  m_w[4] = 2.0/gammalib::pi;
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
343  *
344  * @exception GException::wcs_invalid_x_y
345  * One or more of the (x,y) coordinates were invalid, as indicated
346  * by the stat vector.
347  *
348  * Deproject pixel (x,y) coordinates in the plane of projection to native
349  * spherical coordinates (phi,theta).
350  *
351  * This method has been adapted from the wcslib function prj.c::molx2s()
352  * (version 4.20). The interface follows very closely that of wcslib. In
353  * contrast to the wcslib routine, however, the method assumes that the
354  * projection has been setup previously (as this will be done by the
355  * constructor).
356  ***************************************************************************/
357 void GWcsMOL::prj_x2s(int nx, int ny, int sxy, int spt,
358  const double* x, const double* y,
359  double* phi, double* theta, int* stat) const
360 {
361  // Set tolerance
362  const double tol = 1.0e-12;
363 
364  // Initialize projection if required
365  if (!m_prjset) {
366  prj_set();
367  }
368 
369  // Set value replication length mx,my
370  int mx;
371  int my;
372  if (ny > 0) {
373  mx = nx;
374  my = ny;
375  }
376  else {
377  mx = 1;
378  my = 1;
379  ny = nx;
380  }
381 
382  // Initialise status code and statistics
383  int status = 0;
384  int n_invalid = 0;
385 
386  // Do x dependence
387  const double* xp = x;
388  int rowoff = 0;
389  int rowlen = nx * spt;
390  for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
391  double xj = *xp + m_x0;
392  double s = m_w[3] * xj;
393  double t = std::abs(xj) - tol;
394  double* phip = phi + rowoff;
395  double* thetap = theta + rowoff;
396  for (int iy = 0; iy < my; ++iy, phip += rowlen, thetap += rowlen) {
397  *phip = s;
398  *thetap = t;
399  }
400  }
401 
402  // Do y dependence
403  const double* yp = y;
404  double* phip = phi;
405  double* thetap = theta;
406  int* statp = stat;
407  for (int iy = 0; iy < ny; ++iy, yp += sxy) {
408 
409  // Initialise status
410  int istat = 0;
411 
412  // Compute some stuff
413  double yj = *yp + m_y0;
414  double y0 = yj / m_r0;
415  double r = 2.0 - y0 * y0;
416  double s = 0.0;
417 
418  // Check limits
419  if (r <= tol) {
420  if (r < -tol) {
421  istat = 1;
422  status = 3;
423  n_invalid++;
424  #if defined(G_DEBUG_PRJ)
425  std::cout << "prj_x2s(Phi)..:";
426  std::cout << " nx=" << nx;
427  std::cout << " ny=" << ny;
428  std::cout << " iy=" << iy;
429  std::cout << " yp=" << *yp;
430  std::cout << " yj=" << yj;
431  std::cout << " y0=" << y0;
432  std::cout << " phip=" << *phip;
433  std::cout << " r=" << r;
434  std::cout << " s=" << s;
435  std::cout << std::endl;
436  #endif
437  }
438  else {
439  istat = -1; // OK if fabs(x) < tol whence phi = 0.0
440  }
441  r = 0.0;
442  s = 0.0;
443  }
444  else {
445  r = std::sqrt(r);
446  s = 1.0 / r;
447  }
448 
449  // Compute z and t
450  double z = yj * m_w[2];
451  if (std::abs(z) > 1.0) {
452  if (std::abs(z) > 1.0+tol) {
453  istat = 1;
454  status = 3;
455  n_invalid++;
456  #if defined(G_DEBUG_PRJ)
457  std::cout << "prj_x2s(Phi)..:";
458  std::cout << " nx=" << nx;
459  std::cout << " ny=" << ny;
460  std::cout << " iy=" << iy;
461  std::cout << " yp=" << *yp;
462  std::cout << " yj=" << yj;
463  std::cout << " y0=" << y0;
464  std::cout << " phip=" << *phip;
465  std::cout << " r=" << r;
466  std::cout << " s=" << s;
467  std::cout << " z=" << z;
468  std::cout << std::endl;
469  #endif
470  z = 0.0;
471  }
472  else {
473  z = copysign(1.0, z) + y0 * r / gammalib::pi;
474  }
475  }
476  else {
477  z = std::asin(z) * m_w[4] + y0 * r / gammalib::pi;
478  }
479  if (std::abs(z) > 1.0) {
480  if (std::abs(z) > 1.0+tol) {
481  istat = 1;
482  status = 3;
483  n_invalid++;
484  #if defined(G_DEBUG_PRJ)
485  std::cout << "prj_x2s(Phi)..:";
486  std::cout << " nx=" << nx;
487  std::cout << " ny=" << ny;
488  std::cout << " iy=" << iy;
489  std::cout << " yp=" << *yp;
490  std::cout << " yj=" << yj;
491  std::cout << " y0=" << y0;
492  std::cout << " phip=" << *phip;
493  std::cout << " r=" << r;
494  std::cout << " s=" << s;
495  std::cout << " z=" << z;
496  std::cout << std::endl;
497  #endif
498  z = 0.0;
499  }
500  else {
501  z = copysign(1.0, z);
502  }
503  }
504  double t = std::asin(z) * gammalib::rad2deg;
505 
506  // Do ...
507  for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt, statp += spt) {
508 
509  // Set status flag
510  if (istat < 0 && *thetap >= 0.0) {
511  *statp = 1;
512  status = 3;
513  n_invalid++;
514  }
515  else {
516  *statp = 0;
517  }
518 
519  // Set values
520  *phip *= s;
521  *thetap = t;
522 
523  }
524 
525  } // endfor: did y dependence
526 
527  // Handle status code
528  if (status == 3) {
529  throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
530  }
531 
532  // Do boundary checking
533  status = prj_bchk(1.0e-11, nx, my, spt, phi, theta, stat);
534 
535  // Handle status code
536  if (status != 0) {
537  std::string msg = "One or more of the (x, y) coordinates were invalid "
538  "for the "+name()+" ("+code()+") projection.";
540  }
541 
542  // Return
543  return;
544 }
545 
546 
547 /***********************************************************************//**
548  * @brief Generic spherical-to-pixel projection
549  *
550  * @param[in] nphi Longitude vector length.
551  * @param[in] ntheta Latitude vector length (0=no replication).
552  * @param[in] spt Input vector step.
553  * @param[in] sxy Output vector step.
554  * @param[in] phi Longitude vector of the projected point in native spherical
555  * coordinates [deg].
556  * @param[in] theta Latitude vector of the projected point in native spherical
557  * coordinates [deg].
558  * @param[out] x Vector of projected x coordinates.
559  * @param[out] y Vector of projected y coordinates.
560  * @param[out] stat Status return value for each vector element (always 0)
561  *
562  * @exception GException::wcs_invalid_phi_theta
563  * One or more of the (phi,theta) coordinates were invalid, as
564  * indicated by the stat vector.
565  *
566  * Project native spherical coordinates (phi,theta) to pixel (x,y)
567  * coordinates in the plane of projection.
568  *
569  * This method has been adapted from the wcslib function prj.c::mols2x()
570  * (version 4.20). The interface follows very closely that of wcslib. In
571  * contrast to the wcslib routine, however, the method assumes that the
572  * projection has been setup previously (as this will be done by the
573  * constructor).
574  ***************************************************************************/
575 void GWcsMOL::prj_s2x(int nphi, int ntheta, int spt, int sxy,
576  const double* phi, const double* theta,
577  double* x, double* y, int* stat) const
578 {
579  // Set tolerance
580  const double tol = 1.0e-13;
581 
582  // Initialize projection if required
583  if (!m_prjset) {
584  prj_set();
585  }
586 
587  // Set value replication length mphi,mtheta
588  int mphi;
589  int mtheta;
590  if (ntheta > 0) {
591  mphi = nphi;
592  mtheta = ntheta;
593  }
594  else {
595  mphi = 1;
596  mtheta = 1;
597  ntheta = nphi;
598  }
599 
600  // Initialise status code and statistics
601  int status = 0;
602  int n_invalid = 0;
603 
604  // Do phi dependence
605  const double* phip = phi;
606  int rowoff = 0;
607  int rowlen = nphi * sxy;
608  for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
609  double xi = m_w[1] * (*phip);
610  double* xp = x + rowoff;
611  for (int itheta = 0; itheta < mtheta; ++itheta) {
612  *xp = xi;
613  xp += rowlen;
614  }
615  }
616 
617  // Do theta dependence
618  const double* thetap = theta;
619  double* xp = x;
620  double* yp = y;
621  int* statp = stat;
622  for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
623 
624  // Compute xi and eta
625  double xi = 0.0;
626  double eta = 0.0;
627  if (std::abs(*thetap) == 90.0) {
628  xi = 0.0;
629  eta = copysign(m_w[0], *thetap);
630  }
631  else if (*thetap == 0.0) {
632  xi = 1.0;
633  eta = 0.0;
634  }
635  else {
636  double u = gammalib::pi * std::sin(*thetap * gammalib::deg2rad);
637  double v0 = -gammalib::pi;
638  double v1 = gammalib::pi;
639  double v = u;
640  for (int k = 0; k < 100; ++k) {
641  double resid = (v - u) + std::sin(v);
642  if (resid < 0.0) {
643  if (resid > -tol) {
644  break;
645  }
646  v0 = v;
647  }
648  else {
649  if (resid < tol) {
650  break;
651  }
652  v1 = v;
653  }
654  v = 0.5 * (v0 + v1);
655  }
656  double gamma = 0.5 * v;
657  xi = std::cos(gamma);
658  eta = std::sin(gamma) * m_w[0];
659  }
660  eta -= m_y0;
661 
662  // Do ...
663  for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy, statp += sxy) {
664  *xp = xi * (*xp) - m_x0;
665  *yp = eta;
666  *statp = 0;
667  }
668 
669  } // endfor: theta dependence
670 
671  // Handle status code
672  if (status == 4) {
674  }
675 
676  // Return
677  return;
678 }
Mollweide&#39;s projection class definition.
Definition: GWcsMOL.hpp:43
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
Interface definition for the WCS registry class.
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: GWcsMOL.cpp:357
virtual GWcsMOL * clone(void) const
Clone instance.
Definition: GWcsMOL.cpp:208
virtual std::string name(void) const
Return projection name.
Definition: GWcsMOL.hpp:115
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
const double pi
Definition: GMath.hpp:35
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
void copy_members(const GWcsMOL &wcs)
Copy class members.
Definition: GWcsMOL.cpp:263
Mollweide&#39;s projection class definition.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
void init_members(void)
Initialise class members.
Definition: GWcsMOL.cpp:251
#define G_PRJ_S2X
Definition: GWcsMOL.cpp:40
virtual void clear(void)
Clear instance.
Definition: GWcsMOL.cpp:186
virtual std::string code(void) const
Return projection code.
Definition: GWcsMOL.hpp:101
#define G_PRJ_X2S
Definition: GWcsMOL.cpp:38
Gammalib tools definition.
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
const double sqrt_two
Definition: GMath.hpp:52
virtual std::string print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcsMOL.cpp:221
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
#define copysign(X, Y)
Definition: GWcsMOL.cpp:44
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
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
GWcsMOL & operator=(const GWcsMOL &wcs)
Assignment operator.
Definition: GWcsMOL.cpp:151
std::string wcs_print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GWcs.cpp:1602
void prj_set(void) const
Setup of projection.
Definition: GWcsMOL.cpp:300
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
void free_members(void)
Delete class members.
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
void free_members(void)
Delete class members.
Definition: GWcsMOL.cpp:273
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
const GWcsMOL g_wcs_mol_seed
Definition: GWcsMOL.cpp:56
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
bool m_prjset
Projection is set.
Definition: GWcs.hpp:211
void init_members(void)
Initialise class members.
GWcsMOL(void)
Void constructor.
Definition: GWcsMOL.cpp:69
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: GWcsMOL.cpp:575
const double rad2deg
Definition: GMath.hpp:44
virtual ~GWcsMOL(void)
Destructor.
Definition: GWcsMOL.cpp:129
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
Mathematical function definitions.