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