Francois Brun / gammalib

Commit 7fc53dcf70be8eff9e65685b061a46dac88c1a1e

Authored by knodlseder
1 parent 7f701615

implement dummy methods (to be overloaded by approriate derived class methods if applicable)

Inline Side-by-side
Showing 2 changed files with 285 additions and 40 deletions
include/GWcs.hpp
  View file @7fc53dc
... ... @@ -55,18 +55,18 @@ public:
55 55  
56 56 // Pure virtual methods (not implemented)
57 57 virtual std::string type(void) const = 0;
58   - virtual GSkyDir pix2dir(const int& pix) = 0;
59   - virtual GSkyDir xy2dir(const GSkyPixel& pix) = 0;
60   - virtual int dir2pix(GSkyDir dir) const = 0;
61   - virtual GSkyPixel dir2xy(GSkyDir dir) const = 0;
62   - virtual double omega(const int& pix) const = 0;
63   - virtual double omega(const GSkyPixel& pix) const = 0;
64 58  
65   - // Virtual methods (implemented)
  59 + // Virtual methods
66 60 virtual void read(const GFitsHDU* hdu);
67 61 virtual void write(GFitsHDU* hdu);
68 62 virtual std::string coordsys(void) const;
69 63 virtual void coordsys(const std::string& coordsys);
  64 + virtual double omega(const int& pix) const;
  65 + virtual double omega(const GSkyPixel& pix) const;
  66 + virtual GSkyDir pix2dir(const int& pix);
  67 + virtual int dir2pix(GSkyDir dir) const;
  68 + virtual GSkyDir xy2dir(const GSkyPixel& pix);
  69 + virtual GSkyPixel dir2xy(GSkyDir dir) const;
70 70  
71 71 private:
72 72 // Private methods
... ... @@ -76,13 +76,17 @@ private:
76 76 virtual GWcs* clone(void) const = 0;
77 77  
78 78 protected:
  79 + // Typdefs
  80 + typedef void (GWcs::*_wcspf)(GVector* coord) const;
  81 +
79 82 // Protected methods
80 83 GVector wcs_dir2native(GSkyDir dir) const;
81 84 GSkyDir wcs_native2dir(GVector native) const;
  85 + void wcs_init(const double& theta0);
82 86 GVector wcs_getpole(const double& theta0);
83 87 GMatrix wcs_get_rot(void);
84 88 void dump_wcs(std::ostream& os) const;
85   -
  89 +
86 90 // Protected data area (astr structure)
87 91 int m_coordsys; //!< 0=celestial, 1=galactic
88 92 int m_reverse; //!< Reverse axes (1=true)
... ... @@ -92,7 +96,7 @@ protected:
92 96 GVector m_npole; //!< Native coordinate of North Pole
93 97 GMatrix m_cd; //!< Astrometry parameters (2x2 matrix, deg/pixel)
94 98 GVector m_pv2; //!< Projection parameters (up to 21)
95   -
  99 +
96 100 // Derived parameters
97 101 double m_theta0; //!< Native latitude of the fiducial point
98 102 GVector m_refval; //!< Ordered value of reference point
... ... @@ -101,6 +105,8 @@ protected:
101 105 GVector m_native_pole; //!< Coordinates of native pole
102 106 GMatrix m_rot; //!< Rotation matrix
103 107 GMatrix m_trot; //!< Transpose of rotation matrix
  108 + _wcspf m_std2nat; //!< Standard to native projection function
  109 + _wcspf m_nat2std; //!< Native to standard projection function
104 110 };
105 111  
106 112 #endif /* GWCS_HPP */
... ...
src/sky/GWcs.cpp
  View file @7fc53dc
... ... @@ -9,7 +9,6 @@
9 9 * the Free Software Foundation; either version 2 of the License, or *
10 10 * (at your option) any later version. *
11 11 * *
12   - * ----------------------------------------------------------------------- *
13 12 ***************************************************************************/
14 13 /**
15 14 * @file GWcs.cpp
... ... @@ -29,12 +28,18 @@
29 28 #define G_WCS "GWcs::GWcs(std::string,double,double,double,double,double,"\
30 29 "double)"
31 30 #define G_COORDSYS_SET "GWcs::coordsys(std::string)"
  31 +#define G_OMEGA1 "GWcs::omega(int)"
  32 +#define G_OMEGA2 "GWcs::omega(GSkyPixel)"
  33 +#define G_XY2DIR "GWcs::xy2dir(GSkyPixel)"
  34 +#define G_DIR2XY "GWcs::dir2xy(GSkyDir)"
32 35  
33 36 /* __ Macros _____________________________________________________________ */
34 37  
35 38 /* __ Coding definitions _________________________________________________ */
36 39  
37 40 /* __ Debug definitions __________________________________________________ */
  41 +//#define G_DIR2XY_DEBUG // Debug dir2xy
  42 +//#define G_XY2DIR_DEBUG // Debug xy2dir
38 43  
39 44 /* __ Local prototypes ___________________________________________________ */
40 45  
... ... @@ -254,7 +259,7 @@ std::string GWcs::coordsys(void) const
254 259 *
255 260 * @param[in] coordsys Coordinate system (EQU/CEL/E/C or GAL/G)
256 261 *
257   - * @exception GException::wcs_bad_coords
  262 + * @exception GException::wcs_bad_coords
258 263 * Invalid coordsys parameter.
259 264 *
260 265 * Set coordinate system from std::string. Each of the following is
... ... @@ -280,6 +285,209 @@ void GWcs::coordsys(const std::string&amp; coordsys)
280 285 }
281 286  
282 287  
  288 +/***********************************************************************//**
  289 + * @brief Returns solid angle of pixel (DUMMY METHOD)
  290 + *
  291 + * @param[in] pix Pixel number (0,1,...,m_num_pixels).
  292 + *
  293 + * @exception GException::wcs
  294 + * Method not defined for projection.
  295 + *
  296 + * This method should be overloaded by the appropriate method of the derived
  297 + * class. It thus should never get called.
  298 + ***************************************************************************/
  299 +double GWcs::omega(const int& pix) const
  300 +{
  301 + // Set error message
  302 + std::string message = "GWcs method not defined for " +
  303 + this->type() + " projection.";
  304 +
  305 + // Throw error
  306 + throw GException::wcs(G_OMEGA1, message);
  307 +
  308 + // Define constant direction
  309 + const double omega = 0.0;
  310 +
  311 + // Return solid angle
  312 + return omega;
  313 +}
  314 +
  315 +
  316 +/***********************************************************************//**
  317 + * @brief Returns solid angle of pixel (DUMMY METHOD)
  318 + *
  319 + * @param[in] pix Sky pixel.
  320 + *
  321 + * @exception GException::wcs
  322 + * Method not defined for projection.
  323 + *
  324 + * This method should be overloaded by the appropriate method of the derived
  325 + * class. It thus should never get called.
  326 + ***************************************************************************/
  327 +double GWcs::omega(const GSkyPixel& pix) const
  328 +{
  329 + // Set error message
  330 + std::string message = "GWcs method not defined for " +
  331 + this->type() + " projection.";
  332 +
  333 + // Throw error
  334 + throw GException::wcs(G_OMEGA2, message);
  335 +
  336 + // Define constant direction
  337 + const double omega = 0.0;
  338 +
  339 + // Return solid angle
  340 + return omega;
  341 +}
  342 +
  343 +
  344 +/***********************************************************************//**
  345 + * @brief Returns sky direction of pixel (DUMMY METHOD)
  346 + *
  347 + * @param[in] pix Pixel number (0,1,...,m_num_pixels).
  348 + *
  349 + * @exception GException::wcs
  350 + * Method not defined for projection.
  351 + *
  352 + * This method should be overloaded by the appropriate method of the derived
  353 + * class. It thus should never get called.
  354 + ***************************************************************************/
  355 +GSkyDir GWcs::pix2dir(const int& pix)
  356 +{
  357 + // Set error message
  358 + std::string message = "GWcs method not defined for " +
  359 + this->type() + " projection.";
  360 +
  361 + // Throw error
  362 + throw GException::wcs(G_XY2DIR, message);
  363 +
  364 + // Define constant direction
  365 + const GSkyDir dir;
  366 +
  367 + // Return
  368 + return dir;
  369 +}
  370 +
  371 +
  372 +/***********************************************************************//**
  373 + * @brief Returns pixel for sky direction (DUMMY METHOD)
  374 + *
  375 + * @param[in] dir Sky direction.
  376 + *
  377 + * @exception GException::wcs
  378 + * Method not defined for projection.
  379 + *
  380 + * This method should be overloaded by the appropriate method of the derived
  381 + * class. It thus should never get called.
  382 + ***************************************************************************/
  383 +int GWcs::dir2pix(GSkyDir dir) const
  384 +{
  385 + // Set error message
  386 + std::string message = "GWcs method not defined for " +
  387 + this->type() + " projection.";
  388 +
  389 + // Throw error
  390 + throw GException::wcs(G_DIR2XY, message);
  391 +
  392 + // Return
  393 + return 0;
  394 +}
  395 +
  396 +
  397 +/***********************************************************************//**
  398 + * @brief Returns sky direction of pixel
  399 + *
  400 + * @param[in] pix Sky pixel.
  401 + *
  402 + * @exception GException::wcs_no_proj_fct
  403 + * No WCS projection function has been declared.
  404 + *
  405 + * Note that pixel indices start from 0.
  406 + ***************************************************************************/
  407 +GSkyDir GWcs::xy2dir(const GSkyPixel& pix)
  408 +{
  409 + // Set constants
  410 + const int __permutation[2] = {1,0};
  411 +
  412 + // Throw an error if no projection function exists
  413 + if (this->m_std2nat == NULL)
  414 + throw GException::wcs_no_proj_fct(G_XY2DIR);
  415 +
  416 + // Determine offset
  417 + GVector offset = GVector(pix.x(), pix.y());
  418 +
  419 + // Determine xy
  420 + GVector xy = m_cd * (offset - m_refpix);
  421 +
  422 + // Swap result if coordinates are reversed
  423 + if (m_reverse)
  424 + xy = perm(xy, __permutation);
  425 +
  426 + // Perform CAR map projection
  427 + GVector native = xy;
  428 + (this->*m_nat2std)(&native);
  429 +
  430 + // Get sky direction for native coordinates
  431 + GSkyDir dir = wcs_native2dir(native);
  432 +
  433 + // Debug: Dump transformation steps
  434 + #if defined(G_XY2DIR_DEBUG)
  435 + std::cout << "xy2dir: pixel=" << offset << " xy=" << xy
  436 + << " native=" << native << " dir=" << dir << std::endl;
  437 + #endif
  438 +
  439 + // Return
  440 + return dir;
  441 +}
  442 +
  443 +
  444 +/***********************************************************************//**
  445 + * @brief Returns pixel of sky direction
  446 + *
  447 + * @param[in] dir Sky direction.
  448 + *
  449 + * @exception GException::wcs_no_proj_fct
  450 + * No WCS projection function has been declared.
  451 + *
  452 + * Note that pixel indices start from 0.
  453 + ***************************************************************************/
  454 +GSkyPixel GWcs::dir2xy(GSkyDir dir) const
  455 +{
  456 + // Set constants
  457 + const int __permutation[2] = {1,0};
  458 +
  459 + // Throw an error if no projection function exists
  460 + if (this->m_std2nat == NULL)
  461 + throw GException::wcs_no_proj_fct(G_DIR2XY);
  462 +
  463 + // Get native coordinates for sky direction in radians
  464 + GVector native = wcs_dir2native(dir);
  465 +
  466 + // Perform CAR map projection
  467 + GVector xy = native;
  468 + (this->*m_std2nat)(&xy);
  469 +
  470 + // Swap result if coordinates are reversed
  471 + if (m_reverse)
  472 + xy = perm(xy, __permutation);
  473 +
  474 + // Determine pixel offset
  475 + GVector offset = m_invcd * xy + m_refpix;
  476 +
  477 + // Set sky pixel
  478 + GSkyPixel pixel(offset(0), offset(1));
  479 +
  480 + // Debug: Dump transformation steps
  481 + #if defined(G_DIR2XY_DEBUG)
  482 + std::cout << "dir2xy: dir=" << dir << " native=" << native
  483 + << " xy=" << xy << " pixel=" << offset << std::endl;
  484 + #endif
  485 +
  486 + // Return sky pixel
  487 + return pixel;
  488 +}
  489 +
  490 +
283 491 /*==========================================================================
284 492 = =
285 493 = GWcs protected methods =
... ... @@ -304,17 +512,17 @@ GVector GWcs::wcs_dir2native(GSkyDir dir) const
304 512 phi1 = dir.l();
305 513 theta1 = dir.b();
306 514 }
307   -
  515 +
308 516 // Define right hand side vector for matrix equation
309 517 double ct = cos(theta1);
310 518 double x0 = ct * cos(phi1);
311 519 double x1 = ct * sin(phi1);
312 520 double x2 = sin(theta1);
313 521 GVector x(x0,x1,x2);
314   -
  522 +
315 523 // Find solution
316 524 GVector b = m_rot * x;
317   -
  525 +
318 526 // Account for possible roundoff in 3rd element
319 527 if (b(2) < -1.0) b(2) = -1.0;
320 528 if (b(2) > +1.0) b(2) = +1.0;
... ... @@ -323,7 +531,7 @@ GVector GWcs::wcs_dir2native(GSkyDir dir) const
323 531 double theta = asin(b(2));
324 532 double phi = atan2(b(1),b(0));
325 533 GVector native(phi, theta);
326   -
  534 +
327 535 // Return native coordinates
328 536 return native;
329 537 }
... ... @@ -364,34 +572,66 @@ GSkyDir GWcs::wcs_native2dir(GVector native) const
364 572 dir.radec(phi, theta);
365 573 else
366 574 dir.lb(phi, theta);
367   -
  575 +
368 576 // Return sky direction
369 577 return dir;
370 578 }
371 579  
372 580  
373 581 /***********************************************************************//**
  582 + * @brief Initialises derived WCS parameters
  583 + *
  584 + * @param[in] theta0 Native latitude of the fiducial point.
  585 + ***************************************************************************/
  586 +void GWcs::wcs_init(const double& theta0)
  587 +{
  588 + // Set native latitude of the fiducial point
  589 + m_theta0 = theta0;
  590 +
  591 + // Compute the coordinates of native pole for non-polar projection
  592 + m_native_pole = wcs_getpole(m_theta0);
  593 +
  594 + // Get matrix for rotation between sky and native coordinates
  595 + m_rot = wcs_get_rot();
  596 +
  597 + // Get transpose matrix for inverse rotation
  598 + m_trot = transpose(m_rot);
  599 +
  600 + // Return
  601 + return;
  602 +}
  603 +
  604 +
  605 +/***********************************************************************//**
374 606 * @brief Compute the coordinates of native pole for non-polar projection
375 607 *
376 608 * @param[in] theta0 Native latitude of the fiducial point.
377 609 *
378   - * Returns the coordinates of the native pole for non-polar projections.
379   - * The result is in spherical coordinates and in units of radians.
  610 + * For non-polar (cylindrical or conic) projections, the native pole is
  611 + * not at the reference point, and this method is used to determine the
  612 + * position of the native pole. See section 2.4 of the paper
  613 + * "Representation of Celestial Coordinates in FITS" by
  614 + * Calabretta & Greisen (2002) A&A, 395, 1077, also available at
  615 + * http://www.aoc.nrao.edu/~egreisen.
  616 + * This method returns the coordinates of the native pole for non-polar
  617 + * projections in spherical coordinates. Units are radians.
  618 + * The code has largely be inspired by the IDL routine wcs_getpole.pro
  619 + * (update of May 1998).
380 620 ***************************************************************************/
381 621 GVector GWcs::wcs_getpole(const double& theta0)
382 622 {
383 623 // Declare result
384 624 double alpha_p;
385 625 double delta_p;
386   -
  626 +
387 627 // Compute crval in radians
388 628 double alpha_0 = m_refval(0) * deg2rad;
389 629 double delta_0 = m_refval(1) * deg2rad;
390   -
  630 +
391 631 // Get longitude and latitude of pole
392 632 double lonpole = m_npole(0);
393 633 double latpole = m_npole(1);
394   -
  634 +
395 635 // If theta0=90 then the coordinate of the native pole are given by crval
396 636 if (theta0 == 90.0) {
397 637 alpha_p = alpha_0;
... ... @@ -407,7 +647,7 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
407 647 double sd = sin(delta_0);
408 648 double cd = sin(delta_0);
409 649 double tand = sin(delta_0);
410   -
  650 +
411 651 // Case A: theta0 = 0
412 652 if (theta0 == 0.0) {
413 653 if (delta_0 == 0.0 && lonpole == 90.0)
... ... @@ -423,7 +663,7 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
423 663 else
424 664 alpha_p = alpha_0 - atan2(sp/cd, -tan(delta_p)*tand);
425 665 }
426   -
  666 +
427 667 // Case B: theta0 != 0
428 668 else {
429 669 double ctheta = cos(theta0*deg2rad);
... ... @@ -467,9 +707,9 @@ GVector GWcs::wcs_getpole(const double&amp; theta0)
467 707 } // endelse: sdelt != 1
468 708 } // endelse: cd != 0
469 709 } // endelse: term2 != 0
470   - } // endelse: Case B
  710 + } // endelse: Case B
471 711 } // endelse: theta0 != 90
472   -
  712 +
473 713 // Set result vector
474 714 GVector pole(alpha_p,delta_p);
475 715  
... ... @@ -487,7 +727,7 @@ GMatrix GWcs::wcs_get_rot(void)
487 727 {
488 728 // Allocate rotation matrix
489 729 GMatrix r(3,3);
490   -
  730 +
491 731 // Compute useful quantities relating to reference angles
492 732 double sp = sin(m_npole(0)*deg2rad);
493 733 double cp = cos(m_npole(0)*deg2rad);
... ... @@ -495,7 +735,7 @@ GMatrix GWcs::wcs_get_rot(void)
495 735 double ca = cos(m_native_pole(0));
496 736 double sd = sin(m_native_pole(1));
497 737 double cd = cos(m_native_pole(1));
498   -
  738 +
499 739 // Compute rotation matrix
500 740 r(0,0) = -sa*sp - ca*cp*sd;
501 741 r(0,1) = ca*sp - sa*cp*sd;
... ... @@ -506,7 +746,7 @@ GMatrix GWcs::wcs_get_rot(void)
506 746 r(2,0) = ca*cd;
507 747 r(2,1) = sa*cd;
508 748 r(2,2) = sd;
509   -
  749 +
510 750 // Return rotation matrix
511 751 return r;
512 752 }
... ... @@ -535,9 +775,7 @@ void GWcs::dump_wcs(std::ostream&amp; os) const
535 775 << m_invcd(1,1) << std::endl;
536 776 os << " Coordinate of native pole .: " << m_native_pole*rad2deg
537 777 << " deg" << std::endl;
538   -// os << m_rot << std::endl;
539   -// os << m_trot << std::endl;
540   -
  778 +
541 779 // Return
542 780 return;
543 781 }
... ... @@ -563,7 +801,7 @@ void GWcs::init_members(void)
563 801 m_npole = GVector(180.0, 0.0);
564 802 m_cd = GMatrix(2,2);
565 803 m_pv2 = GVector(21);
566   -
  804 +
567 805 // Initialise derived members
568 806 m_theta0 = 0.0;
569 807 m_refval = GVector(2);
... ... @@ -573,6 +811,10 @@ void GWcs::init_members(void)
573 811 m_rot = GMatrix(3,3);
574 812 m_trot = GMatrix(3,3);
575 813  
  814 + // Initialise projection function pointers
  815 + m_std2nat = NULL;
  816 + m_nat2std = NULL;
  817 +
576 818 // Return
577 819 return;
578 820 }
... ... @@ -594,7 +836,7 @@ void GWcs::copy_members(const GWcs&amp; wcs)
594 836 m_npole = wcs.m_npole;
595 837 m_cd = wcs.m_cd;
596 838 m_pv2 = wcs.m_pv2;
597   -
  839 +
598 840 // Copy derived attributes
599 841 m_theta0 = wcs.m_theta0;
600 842 m_refval = wcs.m_refval;
... ... @@ -604,6 +846,10 @@ void GWcs::copy_members(const GWcs&amp; wcs)
604 846 m_rot = wcs.m_rot;
605 847 m_trot = wcs.m_trot;
606 848  
  849 + // Copy projection function pointers
  850 + m_std2nat = wcs.m_std2nat;
  851 + m_nat2std = wcs.m_nat2std;
  852 +
607 853 // Return
608 854 return;
609 855 }
... ... @@ -624,10 +870,3 @@ void GWcs::free_members(void)
624 870 = GWcs friends =
625 871 = =
626 872 ==========================================================================*/
627   -
628   -
629   -/*==========================================================================
630   - = =
631   - = Other functions used by GWcs =
632   - = =
633   - ==========================================================================*/
... ...