GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GWcsSFL.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GWcsSFL.cpp - Sanson-Flamsteed (SFL) projection class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-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 GWcsSFL.cpp
23  * @brief Sanson-Flamsteed (SFL) 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 "GWcsSFL.hpp"
34 #include "GWcsRegistry.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_PRJ_X2S "GWcsSFL::prj_x2s(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_sfl_registry(&g_wcs_sfl_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 GWcsSFL::GWcsSFL(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  // Initialise class members
92  init_members();
93 
94  // Setup WCS derived parameters
95  wcs_set();
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Copy constructor
104  *
105  * @param[in] wcs World Coordinate System.
106  ***************************************************************************/
107 GWcsSFL::GWcsSFL(const GWcsSFL& wcs) : GWcs(wcs)
108 {
109  // Initialise class members for clean destruction
110  init_members();
111 
112  // Copy members
113  copy_members(wcs);
114 
115  // Return
116  return;
117 }
118 
119 
120 /***********************************************************************//**
121  * @brief Destructor
122  ***************************************************************************/
124 {
125  // Free members
126  free_members();
127 
128  // Return
129  return;
130 }
131 
132 
133 /*==========================================================================
134  = =
135  = Operators =
136  = =
137  ==========================================================================*/
138 
139 /***********************************************************************//**
140  * @brief Assignment operator
141  *
142  * @param[in] wcs World Coordinate System.
143  * @return World Coordinate System.
144  ***************************************************************************/
146 {
147  // Execute only if object is not identical
148  if (this != &wcs) {
149 
150  // Copy base class members
151  this->GWcs::operator=(wcs);
152 
153  // Free members
154  free_members();
155 
156  // Initialise private members for clean destruction
157  init_members();
158 
159  // Copy members
160  copy_members(wcs);
161 
162  } // endif: object was not identical
163 
164  // Return this object
165  return *this;
166 }
167 
168 
169 /*==========================================================================
170  = =
171  = Public methods =
172  = =
173  ==========================================================================*/
174 
175 /***********************************************************************//**
176  * @brief Clear Sanson-Flamsteed projection
177  *
178  * Resets the Sanson-Flamsteed projection to an clean initial state.
179  ***************************************************************************/
180 void GWcsSFL::clear(void)
181 {
182  // Free class members (base and derived classes, derived class first)
183  free_members();
184  this->GWcs::free_members();
186 
187  // Initialise members
189  this->GWcs::init_members();
190  init_members();
191 
192  // Return
193  return;
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Clone Sanson-Flamsteed projection
199  *
200  * @return Pointer to deep copy of Sanson-Flamsteed projection.
201  ***************************************************************************/
203 {
204  return new GWcsSFL(*this);
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Print Sanson-Flamsteed projection information
210  *
211  * @param[in] chatter Chattiness.
212  * @return String containing Sanson-Flamsteed projection information.
213  ***************************************************************************/
214 std::string GWcsSFL::print(const GChatter& chatter) const
215 {
216  // Initialise result string
217  std::string result;
218 
219  // Continue only if chatter is not silent
220  if (chatter != SILENT) {
221 
222  // Append header
223  result.append("=== GWcsSFL ===");
224 
225  // Append information
226  result.append(wcs_print(chatter));
227 
228  } // endif: chatter was not silent
229 
230  // Return result
231  return result;
232 }
233 
234 
235 /*==========================================================================
236  = =
237  = Private methods =
238  = =
239  ==========================================================================*/
240 
241 /***********************************************************************//**
242  * @brief Initialise class members
243  ***************************************************************************/
245 {
246  // Return
247  return;
248 }
249 
250 
251 /***********************************************************************//**
252  * @brief Copy class members
253  *
254  * @param[in] wcs World Coordinate System.
255  ***************************************************************************/
257 {
258  // Return
259  return;
260 }
261 
262 
263 /***********************************************************************//**
264  * @brief Delete class members
265  ***************************************************************************/
267 {
268  // Return
269  return;
270 }
271 
272 
273 /***********************************************************************//**
274  * @brief Setup of projection
275  *
276  * This method sets up the projection information. The method has been
277  * adapted from the wcslib function prj.c::sflset.
278  *
279  * Given and/or returned:
280  * m_r0 Reset to 180/pi if 0.
281  * m_phi0 Reset to 0.0 if undefined.
282  * m_theta0 Reset to 0.0 if undefined.
283  *
284  * Returned:
285  * m_x0 Fiducial offset in x.
286  * m_y0 Fiducial offset in y.
287  * m_w[0] r0*(pi/180)
288  * m_w[1] (180/pi)/r0
289  ***************************************************************************/
290 void GWcsSFL::prj_set(void) const
291 {
292  // Signal that projection has been set (needs to be done before calling
293  // the prj_off() method to avoid an endless loop)
294  m_prjset = true;
295 
296  // Initialise projection parameters
297  m_w.assign(2, 0.0);
298 
299  // Precompute
300  if (m_r0 == 0.0) {
302  m_w[0] = 1.0;
303  m_w[1] = 1.0;
304  }
305  else {
306  m_w[0] = m_r0 * gammalib::deg2rad;
307  m_w[1] = 1.0 / m_w[0];
308  }
309 
310  // Compute fiducial offset
311  prj_off(0.0, 0.0);
312 
313  // Return
314  return;
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Pixel-to-spherical deprojection
320  *
321  * @param[in] nx X vector length.
322  * @param[in] ny Y vector length (0=no replication).
323  * @param[in] sxy Input vector step.
324  * @param[in] spt Output vector step.
325  * @param[in] x Vector of projected x coordinates.
326  * @param[in] y Vector of projected y coordinates.
327  * @param[out] phi Longitude of the projected point in native spherical
328  * coordinates [deg].
329  * @param[out] theta Latitude of the projected point in native spherical
330  * coordinates [deg].
331  * @param[out] stat Status return value for each vector element (always 0)
332  *
333  * Deproject pixel (x,y) coordinates in the plane of projection to native
334  * spherical coordinates (phi,theta).
335  *
336  * This method has been adapted from the wcslib function prj.c::sflx2s().
337  * The interface follows very closely that of wcslib. In contrast to the
338  * wcslib routine, however, the method assumes that the projection has been
339  * setup previously (as this will be done by the constructor).
340  ***************************************************************************/
341 void GWcsSFL::prj_x2s(int nx, int ny, int sxy, int spt,
342  const double* x, const double* y,
343  double* phi, double* theta, int* stat) const
344 {
345  // Initialize projection if required
346  if (!m_prjset) {
347  prj_set();
348  }
349 
350  // Set value replication length mx,my
351  int mx;
352  int my;
353  if (ny > 0) {
354  mx = nx;
355  my = ny;
356  }
357  else {
358  mx = 1;
359  my = 1;
360  ny = nx;
361  }
362 
363  // Initialise status code and statistics
364  int status = 0;
365  int n_invalid = 0;
366 
367 
368  // Do x dependence
369  const double* xp = x;
370  int rowoff = 0;
371  int rowlen = nx * spt;
372  for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
373  double s = m_w[1] * (*xp + m_x0);
374  double* phip = phi + rowoff;
375  for (int iy = 0; iy < my; ++iy, phip += rowlen) {
376  *phip = s;
377  }
378  }
379 
380  // Do y dependence
381  const double* yp = y;
382  double* phip = phi;
383  double* thetap = theta;
384  int* statp = stat;
385  for (int iy = 0; iy < ny; ++iy, yp += sxy) {
386  double yj = *yp + m_y0;
387  double s = std::cos(yj/m_r0);
388 
389  // Initialise status
390  int istat = 0;
391 
392  // Invert s
393  if (s == 0.0) {
394  istat = 1;
395  status = 3;
396  n_invalid++;
397  #if defined(G_DEBUG_PRJ)
398  std::cout << "prj_x2s(Phi)..:";
399  std::cout << " nx=" << nx;
400  std::cout << " ny=" << ny;
401  std::cout << " ix=" << ix;
402  std::cout << " iy=" << iy;
403  std::cout << " yp=" << *yp;
404  std::cout << " yj=" << yj;
405  std::cout << " phip=" << *phip;
406  std::cout << " s=" << s;
407  std::cout << std::endl;
408  #endif
409  }
410  else {
411  s = 1.0 / s;
412  }
413 
414  // ...
415  double t = m_w[1] * yj;
416 
417  // ...
418  for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
419  *phip *= s;
420  *thetap = t;
421  *(statp++) = istat;
422  }
423 
424  } // endfor: y dependence
425 
426  // Handle status code
427  if (status == 3) {
428  throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
429  }
430 
431  // Do boundary checking
432  status = prj_bchk(1.0e-12, nx, my, spt, phi, theta, stat);
433 
434  // Handle status code
435  if (status != 0) {
436  std::string msg = "One or more of the (x, y) coordinates were invalid "
437  "for the "+name()+" ("+code()+") projection.";
439  }
440 
441  // Return
442  return;
443 }
444 
445 
446 /***********************************************************************//**
447  * @brief Generic spherical-to-pixel projection
448  *
449  * @param[in] nphi Longitude vector length.
450  * @param[in] ntheta Latitude vector length (0=no replication).
451  * @param[in] spt Input vector step.
452  * @param[in] sxy Output vector step.
453  * @param[in] phi Longitude vector of the projected point in native spherical
454  * coordinates [deg].
455  * @param[in] theta Latitude vector of the projected point in native spherical
456  * coordinates [deg].
457  * @param[out] x Vector of projected x coordinates.
458  * @param[out] y Vector of projected y coordinates.
459  * @param[out] stat Status return value for each vector element (always 0)
460  *
461  * Project native spherical coordinates (phi,theta) to pixel (x,y)
462  * coordinates in the plane of projection.
463  *
464  * This method has been adapted from the wcslib function prj.c::sfls2x().
465  * The interface follows very closely that of wcslib. In contrast to the
466  * wcslib routine, however, the method assumes that the projection has been
467  * setup previously (as this will be done by the constructor).
468  ***************************************************************************/
469 void GWcsSFL::prj_s2x(int nphi, int ntheta, int spt, int sxy,
470  const double* phi, const double* theta,
471  double* x, double* y, int* stat) const
472 {
473  // Initialize projection if required
474  if (!m_prjset) {
475  prj_set();
476  }
477 
478  // Set value replication length mphi,mtheta
479  int mphi;
480  int mtheta;
481  if (ntheta > 0) {
482  mphi = nphi;
483  mtheta = ntheta;
484  }
485  else {
486  mphi = 1;
487  mtheta = 1;
488  ntheta = nphi;
489  }
490 
491  // Do phi dependence
492  const double* phip = phi;
493  int rowoff = 0;
494  int rowlen = nphi * sxy;
495  for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
496  double xi = m_w[0] * (*phip);
497  double* xp = x + rowoff;
498  for (int itheta = 0; itheta < mtheta; ++itheta, xp += rowlen) {
499  *xp = xi;
500  }
501  }
502 
503  // Do theta dependence
504  const double* thetap = theta;
505  double* xp = x;
506  double* yp = y;
507  int* statp = stat;
508  for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
509  double xi = gammalib::cosd(*thetap);
510  double eta = m_w[0] * (*thetap) - m_y0;
511  for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
512  *xp = xi * (*xp) - m_x0;
513  *yp = eta;
514  *(statp++) = 0;
515  }
516  }
517 
518  // Return
519  return;
520 }
virtual std::string name(void) const
Return projection name.
Definition: GWcsSFL.hpp:115
void init_members(void)
Initialise class members.
Definition: GWcs.cpp:949
Interface definition for the WCS registry class.
void init_members(void)
Initialise class members.
Definition: GWcsSFL.cpp:244
Sanson-Flamsteed (SFL) projection class definition.
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
GWcsSFL(void)
Void constructor.
Definition: GWcsSFL.cpp:64
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
virtual ~GWcsSFL(void)
Destructor.
Definition: GWcsSFL.cpp:123
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: GWcsSFL.cpp:341
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: GWcsSFL.cpp:469
void free_members(void)
Delete class members.
Definition: GWcs.cpp:1053
Sanson-Flamsteed (SFL) projection class definition.
Definition: GWcsSFL.hpp:43
World Coordinate Projection registry class interface definition.
const double deg2rad
Definition: GMath.hpp:43
void wcs_set(void) const
Setup of World Coordinate System.
Definition: GWcs.cpp:1249
virtual GWcsSFL * clone(void) const
Clone Sanson-Flamsteed projection.
Definition: GWcsSFL.cpp:202
virtual GWcs & operator=(const GWcs &wcs)
Assignment operator.
Definition: GWcs.cpp:180
virtual std::string print(const GChatter &chatter=NORMAL) const
Print Sanson-Flamsteed projection information.
Definition: GWcsSFL.cpp:214
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:1602
virtual void clear(void)
Clear Sanson-Flamsteed projection.
Definition: GWcsSFL.cpp:180
GWcsSFL & operator=(const GWcsSFL &wcs)
Assignment operator.
Definition: GWcsSFL.cpp:145
void free_members(void)
Delete class members.
std::vector< double > m_w
Intermediate values.
Definition: GWcs.hpp:217
const GWcsSFL g_wcs_sfl_seed
Definition: GWcsSFL.cpp:51
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
double cosd(const double &angle)
Compute cosine of angle in degrees.
Definition: GMath.cpp:134
void prj_set(void) const
Setup of projection.
Definition: GWcsSFL.cpp:290
#define G_PRJ_X2S
Definition: GWcsSFL.cpp:37
virtual std::string code(void) const
Return projection code.
Definition: GWcsSFL.hpp:101
Exception handler interface definition.
bool m_prjset
Projection is set.
Definition: GWcs.hpp:211
void init_members(void)
Initialise class members.
void copy_members(const GWcsSFL &wcs)
Copy class members.
Definition: GWcsSFL.cpp:256
void free_members(void)
Delete class members.
Definition: GWcsSFL.cpp:266
const double rad2deg
Definition: GMath.hpp:44
double m_r0
Radius of the generating sphere.
Definition: GWcs.hpp:212
Mathematical function definitions.