Eirik Hatlen / gammalib

Commit ca6d6d37239a755eb77055fd1ba8cd78173555d3

Authored by Jürgen Knödlseder
1 parent 1c7cb600

Adjust sky region classes (#3248) 

Renamed the GSkyRegionRect::get_corner() method to GSkyRegionRect::corner() and make local coordinate methods provide. Use GSkyPixel to carry local coordinates in units of radians. Extent the GSkyRegionRect::contains() and GSkyRegionRect::overlaps() method to cover all other possible sky region types.

Adjust GSkyRegionCircle and GSkyRegionMap to the new interface and make sure that their contains() and overlaps() method deal with all other regions.

Adjust the GSkyRegionRect unit test accordingly.
Inline Side-by-side
Showing 6 changed files with 318 additions and 229 deletions
include/GSkyRegionRect.hpp
  View file @ca6d6d3
... ... @@ -34,6 +34,7 @@
34 34 #include "GSkyRegion.hpp"
35 35  
36 36 /* __ Forward declarations _______________________________________________ */
  37 +class GSkyPixel;
37 38  
38 39  
39 40 /***********************************************************************//**
... ... @@ -97,26 +98,26 @@ public:
97 98 void read(const std::string& line);
98 99 std::string write(void) const;
99 100 bool contains(const GSkyDir& dir) const;
100   - bool contains_local(const GSkyDir& locdir) const;
101 101 bool contains(const GSkyRegion& reg) const;
102 102 bool overlaps(const GSkyRegion& reg) const;
103   - GSkyDir transform_to_local(const GSkyDir& skydir) const;
104   - GSkyDir transform_to_global(const GSkyDir& locdir) const;
105   - GSkyDir get_corner(const int& index) const;
  103 + GSkyDir corner(const int& index) const;
106 104 std::string print(const GChatter& chatter = NORMAL) const;
107 105  
108 106 protected:
109 107 // Protected methods
110   - void init_members(void);
111   - void copy_members(const GSkyRegionRect& region);
112   - void free_members(void);
113   - void compute_solid_angle(void);
  108 + void init_members(void);
  109 + void copy_members(const GSkyRegionRect& region);
  110 + void free_members(void);
  111 + void compute_solid_angle(void);
  112 + bool contains(const GSkyPixel& local) const;
  113 + GSkyPixel dir_to_local(const GSkyDir& dir) const;
  114 + GSkyDir local_to_dir(const GSkyPixel& local) const;
114 115  
115 116 // Protected members
116   - GSkyDir m_centre; //!< Centre or reference point of the region
117   - double m_width; //!< Width of the region (degrees)
118   - double m_height; //!< Height of the region (degrees)
119   - double m_posang; //!< Position angle, counterclockwise from North (degrees)
  117 + GSkyDir m_centre; //!< Centre or reference point of the region
  118 + double m_width; //!< Width of the region (degrees)
  119 + double m_height; //!< Height of the region (degrees)
  120 + double m_posang; //!< Position angle, counterclockwise from North (degrees)
120 121 };
121 122  
122 123  
... ...
pyext/GSkyRegionRect.i
  View file @ca6d6d3
... ... @@ -70,12 +70,9 @@ public:
70 70 void read(const std::string& line);
71 71 std::string write(void) const;
72 72 bool contains(const GSkyDir& dir) const;
73   - bool contains_local(const GSkyDir& locdir) const;
74 73 bool contains(const GSkyRegion& reg) const;
75 74 bool overlaps(const GSkyRegion& reg) const;
76   - GSkyDir transform_to_local(const GSkyDir& skydir) const;
77   - GSkyDir transform_to_global(const GSkyDir& locdir) const;
78   - GSkyDir get_corner(const int& index) const;
  75 + GSkyDir corner(const int& index) const;
79 76 };
80 77  
81 78  
... ...
src/sky/GSkyRegionCircle.cpp
  View file @ca6d6d3
... ... @@ -526,47 +526,78 @@ bool GSkyRegionCircle::contains(const GSkyDir&amp; dir) const
526 526 bool GSkyRegionCircle::contains(const GSkyRegion& reg) const
527 527 {
528 528 // Initialise return value
529   - bool fully_inside = false;
  529 + bool contains = false;
530 530  
531   - // If other region is circle use a simple way to calculate
  531 + // If the other region is a circle then check whether it is fully
  532 + // container in the circle
532 533 if (reg.type() == "Circle") {
533 534  
534 535 // Create circular region from reg
535   - const GSkyRegionCircle* regcirc =
536   - dynamic_cast<const GSkyRegionCircle*>(&reg);
  536 + const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(&reg);
537 537  
538 538 // Calculate angular distance between the centres
539   - double ang_dist = m_centre.dist_deg(regcirc->centre());
  539 + double ang_dist = m_centre.dist_deg(circle->centre());
540 540  
541 541 // Check if the region is contained in this
542   - if ((ang_dist + regcirc->radius()) <= m_radius) {
543   - fully_inside = true;
  542 + if ((ang_dist + circle->radius()) <= m_radius) {
  543 + contains = true;
544 544 }
545 545  
546   - }
  546 + } // endif: other region was of type "Circle"
547 547  
548   - // If other region is rectangle check corners
  548 + // ... otherwise, if the other region is a rectangle then check whether
  549 + // all four corners are contained within the circle
549 550 else if (reg.type() == "Rect") {
550 551  
551 552 // Create rectangular region from reg
552   - const GSkyRegionRect* regrect =
553   - dynamic_cast<const GSkyRegionRect*>(&reg);
554   -
555   - // Check containment of all edges
556   - fully_inside = contains(regrect->get_corner(0)) &&
557   - contains(regrect->get_corner(1)) &&
558   - contains(regrect->get_corner(2)) &&
559   - contains(regrect->get_corner(3));
  553 + const GSkyRegionRect* rect = dynamic_cast<const GSkyRegionRect*>(&reg);
  554 +
  555 + // Check containment of all corners
  556 + contains = this->contains(rect->corner(0)) &&
  557 + this->contains(rect->corner(1)) &&
  558 + this->contains(rect->corner(2)) &&
  559 + this->contains(rect->corner(3));
  560 +
560 561 }
561 562  
  563 + // ... otherwise, if the other region is a map then check whether it
  564 + // is contained within the rectangle
  565 + else if (reg.type() == "Map") {
  566 +
  567 + // Create map from reg
  568 + const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(&reg);
  569 +
  570 + // Get non-zero indices
  571 + std::vector<int> indices = map->nonzero_indices();
  572 +
  573 + // Initialise containment flag to true
  574 + contains = true;
  575 +
  576 + // Loop over all indices
  577 + for (int i = 0; i < indices.size(); ++i) {
  578 +
  579 + // Get sky direction of map pixel
  580 + GSkyDir dir = map->map().inx2dir(i);
  581 +
  582 + // If pixel is not contained then set containment flag to false
  583 + // and break
  584 + if (!this->contains(dir)) {
  585 + contains = false;
  586 + break;
  587 + }
  588 +
  589 + } // endfor: looped over non-zero map indices
  590 +
  591 + } // endif: other region was of type "Map"
  592 +
562 593 // ... otherwise throw an exception
563 594 else {
564 595 throw GException::feature_not_implemented(G_CONTAINS,
565 596 "Cannot compare to region type \""+reg.type()+"\" yet.");
566 597 }
567 598  
568   - // Return value
569   - return fully_inside;
  599 + // Return containment flag
  600 + return contains;
570 601 }
571 602  
572 603  
... ... @@ -590,14 +621,13 @@ bool GSkyRegionCircle::overlaps(const GSkyRegion&amp; reg) const
590 621 if (reg.type() == "Circle") {
591 622  
592 623 // Create circular region from reg
593   - const GSkyRegionCircle* regcirc =
594   - dynamic_cast<const GSkyRegionCircle*>(&reg);
  624 + const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(&reg);
595 625  
596 626 // Calculate angular distance between the centres
597   - double ang_dist = m_centre.dist_deg(regcirc->centre());
  627 + double ang_dist = m_centre.dist_deg(circle->centre());
598 628  
599 629 // Check if the distance is smaller than the sum of both radii
600   - if (ang_dist <= (m_radius + regcirc->radius())) {
  630 + if (ang_dist <= (m_radius + circle->radius())) {
601 631 overlap = true;
602 632 }
603 633  
... ... @@ -607,11 +637,10 @@ bool GSkyRegionCircle::overlaps(const GSkyRegion&amp; reg) const
607 637 else if (reg.type() == "Rect") {
608 638  
609 639 // Create rectangular region from reg
610   - const GSkyRegionRect* regrect =
611   - dynamic_cast<const GSkyRegionRect*>(&reg);
  640 + const GSkyRegionRect* rect = dynamic_cast<const GSkyRegionRect*>(&reg);
612 641  
613 642 // Check overlap with circle
614   - overlap = regrect->overlaps(*this);
  643 + overlap = rect->overlaps(*this);
615 644  
616 645 } // endif: region was rectangle
617 646  
... ... @@ -619,11 +648,10 @@ bool GSkyRegionCircle::overlaps(const GSkyRegion&amp; reg) const
619 648 else if (reg.type() == "Map") {
620 649  
621 650 // Create map from reg
622   - const GSkyRegionMap* regmap =
623   - dynamic_cast<const GSkyRegionMap*>(&reg);
  651 + const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(&reg);
624 652  
625 653 // Check overlap with circle
626   - overlap = regmap->overlaps(*this);
  654 + overlap = map->overlaps(*this);
627 655  
628 656 } // endif: region was map
629 657  
... ... @@ -638,7 +666,6 @@ bool GSkyRegionCircle::overlaps(const GSkyRegion&amp; reg) const
638 666 }
639 667  
640 668  
641   -
642 669 /*==========================================================================
643 670 = =
644 671 = Private methods =
... ...
src/sky/GSkyRegionMap.cpp
  View file @ca6d6d3
... ... @@ -677,7 +677,7 @@ void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect)
677 677 for (int icorner = 0; icorner < 4; ++icorner) {
678 678  
679 679 // Get skydir of current corner
680   - GSkyDir corner = rect->get_corner(icorner);
  680 + GSkyDir corner = rect->corner(icorner);
681 681  
682 682 // Compute polar props in global coords relative to rectangle centre
683 683 double dist = rect->centre().dist(corner);
... ...
src/sky/GSkyRegionRect.cpp
  View file @ca6d6d3
... ... @@ -39,7 +39,7 @@
39 39 #define G_READ "GSkyRegionRect::read(std::string&)"
40 40 #define G_CONTAINS "GSkyRegionRect::contains(GSkyRegion&)"
41 41 #define G_OVERLAPS "GSkyRegionRect::overlaps(GSkyRegion&)"
42   -#define G_GET_CORNER "GSkyRegionRect::get_corner(int&)"
  42 +#define G_CORNER "GSkyRegionRect::corner(int&)"
43 43  
44 44 /* __ Macros _____________________________________________________________ */
45 45  
... ... @@ -102,8 +102,8 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir&amp; centre,
102 102 /***********************************************************************//**
103 103 * @brief Right Ascension and Declination constructor
104 104 *
105   - * @param[in] ra Right Ascension of region centre [deg].
106   - * @param[in] dec Declination of region centre [deg].
  105 + * @param[in] ra Right Ascension of region centre (degrees).
  106 + * @param[in] dec Declination of region centre (degrees).
107 107 * @param[in] width Region width (degrees).
108 108 * @param[in] height Region height (degrees).
109 109 * @param[in] posang Position angle (degrees).
... ... @@ -258,7 +258,7 @@ GSkyRegionRect* GSkyRegionRect::clone(void) const
258 258 /***********************************************************************//**
259 259 * @brief Set width of rectangular region
260 260 *
261   - * @param[in] width Width (degrees).
  261 + * @param[in] width Rectangle width (degrees).
262 262 *
263 263 * @exception GException::invalid_argument
264 264 * Width value is less than 0.
... ... @@ -291,7 +291,7 @@ void GSkyRegionRect::width(const double&amp; width)
291 291 /***********************************************************************//**
292 292 * @brief Set height of rectangular region
293 293 *
294   - * @param[in] height Height (degrees).
  294 + * @param[in] height Rectangle height (degrees).
295 295 *
296 296 * @exception GException::invalid_argument
297 297 * Height value is less than 0.
... ... @@ -488,43 +488,14 @@ std::string GSkyRegionRect::write(void) const
488 488 ***************************************************************************/
489 489 bool GSkyRegionRect::contains(const GSkyDir& dir) const
490 490 {
491   - // Compute sky direction in local coordinate system
492   - GSkyDir locdir = transform_to_local(dir);
  491 + // Transform sky direction to local coordinate system
  492 + GSkyPixel local = dir_to_local(dir);
493 493  
494 494 // Check containment using local coordinate system
495   - bool dir_is_in = contains_local(locdir);
  495 + bool contains = this->contains(local);
496 496  
497   - // Return bool
498   - return dir_is_in;
499   -}
500   -
501   -
502   -/***********************************************************************//**
503   - * @brief Checks if local direction lies within region
504   - *
505   - * @param[in] locdir Sky direction in local coordinate system.
506   - *
507   - * A local direction lies within a region when its distance to the region
508   - * centre is not larger than the region extension in both axes directions.
509   - ***************************************************************************/
510   -bool GSkyRegionRect::contains_local(const GSkyDir& locdir) const
511   -{
512   - // Initialise return value
513   - bool dir_is_in = false;
514   -
515   - // Check containment: declination axis
516   - if ((std::abs(locdir.dec_deg()) - 1.0e-10) <= 0.5 * m_height) {
517   -
518   - // Check containment: right ascension axis
519   - if ((std::abs(locdir.ra_deg()) - 1.0e-10) <= 0.5 * m_width) {
520   -
521   - // Sky direction is inside the rectangle
522   - dir_is_in = true;
523   - }
524   - }
525   -
526   - // Return bool
527   - return dir_is_in;
  497 + // Return containment flag
  498 + return contains;
528 499 }
529 500  
530 501  
... ... @@ -539,60 +510,88 @@ bool GSkyRegionRect::contains_local(const GSkyDir&amp; locdir) const
539 510 bool GSkyRegionRect::contains(const GSkyRegion& reg) const
540 511 {
541 512 // Initialise return value
542   - bool is_fully_inside = false;
  513 + bool contains = false;
543 514  
544   - // If other region is circle use a simple way to calculate
  515 + // If the other region is a circle then check whether it is fully
  516 + // contained within the rectangle
545 517 if (reg.type() == "Circle") {
546 518  
547 519 // Create circular region from reg
548   - const GSkyRegionCircle* regcirc =
549   - dynamic_cast<const GSkyRegionCircle*>(&reg);
  520 + const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(&reg);
550 521  
551 522 // Transform circle center to local coordinate system
552   - GSkyDir local_centre = transform_to_local(regcirc->centre());
553   -
554   - // Get circle radius
555   - const double circ_rad = regcirc->radius();
556   -
557   - // Check extension of circle along declination axis
558   - if ((std::abs(local_centre.dec_deg()) + circ_rad) <= 0.5 * m_height) {
  523 + GSkyPixel local = dir_to_local(circle->centre());
  524 + double x_deg = local.x() * gammalib::rad2deg;
  525 + double y_deg = local.y() * gammalib::rad2deg;
559 526  
560   - // Check extension of circle along right ascension axis
561   - if ((std::abs(local_centre.ra_deg()) + circ_rad) <= 0.5 * m_width) {
  527 + // Get circle radius in degrees
  528 + double radius = circle->radius();
562 529  
563   - // Circle is fully contained in this rectangle
564   - is_fully_inside = true;
  530 + // Check whether circle is contained within rectangle
  531 + if ((std::abs(y_deg) + radius) <= 0.5 * m_height) {
  532 + if ((std::abs(x_deg) + radius) <= 0.5 * m_width) {
  533 + contains = true;
565 534 }
566 535 }
567   - } // Region was of type "Circle"
568 536  
569   - // If other region is rectangle use a simple way to calculate
  537 + } // endif: other region was of type "Circle"
  538 +
  539 + // ... otherwise, if the other region is a rectangle then check whether
  540 + // all corners fall within the rectangle
570 541 else if (reg.type() == "Rect") {
571 542  
572 543 // Create rectangular region from reg
573   - const GSkyRegionRect* regrect =
574   - dynamic_cast<const GSkyRegionRect*>(&reg);
  544 + const GSkyRegionRect* rect = dynamic_cast<const GSkyRegionRect*>(&reg);
  545 +
  546 + // Initialise containment flag to true
  547 + contains = true;
575 548  
576 549 // Loop over the four corners of the rectangle :regrect:
577   - for (int icorner = 0; icorner < 4; ++icorner) {
  550 + for (int i = 0; i < 4; ++i) {
578 551  
579   - // Get skydir of current corner
580   - GSkyDir corner = regrect->get_corner(icorner);
  552 + // Get sky direction of current corner
  553 + GSkyDir corner = rect->corner(i);
581 554  
582   - // Check corner containment
583   - if (!contains(corner)) {
  555 + // If corner is not contained then set containment flag to false
  556 + // and break
  557 + if (!this->contains(corner)) {
  558 + contains = false;
584 559 break;
585 560 }
586   - // Did we arrive at the last corner?
587   - else if (icorner == 3) {
588 561  
589   - // As we reached this point, the three prior corners were inside
590   - // and the last corner now is also inside:
591   - // Rectangle is fully inside
592   - is_fully_inside = true;
  562 + } // endfor: looped over the four corners
  563 +
  564 + } // endif: other region was of type "Rect"
  565 +
  566 + // ... otherwise, if the other region is a map then check whether it
  567 + // is contained within the rectangle
  568 + else if (reg.type() == "Map") {
  569 +
  570 + // Create map from reg
  571 + const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(&reg);
  572 +
  573 + // Get non-zero indices
  574 + std::vector<int> indices = map->nonzero_indices();
  575 +
  576 + // Initialise containment flag to true
  577 + contains = true;
  578 +
  579 + // Loop over all indices
  580 + for (int i = 0; i < indices.size(); ++i) {
  581 +
  582 + // Get sky direction of map pixel
  583 + GSkyDir dir = map->map().inx2dir(i);
  584 +
  585 + // If pixel is not contained then set containment flag to false
  586 + // and break
  587 + if (!this->contains(dir)) {
  588 + contains = false;
  589 + break;
593 590 }
594   - } // Looped over the four corners
595   - } // Region was of type "Rect"
  591 +
  592 + } // endfor: looped over non-zero map indices
  593 +
  594 + } // endif: other region was of type "Map"
596 595  
597 596 // ... otherwise throw an exception
598 597 else {
... ... @@ -600,8 +599,8 @@ bool GSkyRegionRect::contains(const GSkyRegion&amp; reg) const
600 599 "Cannot compare to region type \""+reg.type()+"\" yet.");
601 600 }
602 601  
603   - // Return value
604   - return is_fully_inside;
  602 + // Return containment flag
  603 + return contains;
605 604 }
606 605  
607 606  
... ... @@ -614,36 +613,34 @@ bool GSkyRegionRect::contains(const GSkyRegion&amp; reg) const
614 613 * Regions differ in type.
615 614 *
616 615 * @todo Improve implementation for rectangle-rectangle
617   - * @todo Implement check for map
618 616 ***************************************************************************/
619 617 bool GSkyRegionRect::overlaps(const GSkyRegion& reg) const
620 618 {
621 619 // Initialise return value
622 620 bool overlap = false;
623 621  
624   - // If region is circle then check overlap between circle and rectangle
  622 + // If the region is a circle then check overlap between circle and
  623 + // rectangle
625 624 if (reg.type() == "Circle") {
626 625  
627 626 // Create circular region from reg
628   - const GSkyRegionCircle* regcirc =
629   - dynamic_cast<const GSkyRegionCircle*>(&reg);
  627 + const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(&reg);
630 628  
631 629 // Transform circle center to local coordinate system
632   - GSkyDir local_centre = transform_to_local(regcirc->centre());
633   -
634   - // Get circle radius
635   - const double circ_rad = regcirc->radius();
  630 + GSkyPixel local = dir_to_local(circle->centre());
  631 + double x_deg = local.x() * gammalib::rad2deg;
  632 + double y_deg = local.y() * gammalib::rad2deg;
636 633  
637   - // Check extension of circle along declination axis
638   - if (std::abs(local_centre.dec_deg()) <= (0.5 * m_height + circ_rad)) {
  634 + // Get circle radius in degrees
  635 + double radius = circle->radius();
639 636  
640   - // Check extension of circle along right ascension axis
641   - if (std::abs(local_centre.ra_deg()) <= (0.5 * m_width + circ_rad)) {
642   -
643   - // Circle is overlapping with rectangle
  637 + // Check whether circle overlaps within rectangle
  638 + if (std::abs(y_deg) <= (0.5 * m_height + radius)) {
  639 + if (std::abs(x_deg) <= (0.5 * m_width + radius)) {
644 640 overlap = true;
645 641 }
646 642 }
  643 +
647 644 } // endif: region was of type "Circle"
648 645  
649 646 // ... otherwise if region is a rectangle then check overlap between
... ... @@ -651,12 +648,11 @@ bool GSkyRegionRect::overlaps(const GSkyRegion&amp; reg) const
651 648 else if (reg.type() == "Rect") {
652 649  
653 650 // Create rectangular region from reg
654   - const GSkyRegionRect* regrect =
655   - dynamic_cast<const GSkyRegionRect*>(&reg);
  651 + const GSkyRegionRect* rect = dynamic_cast<const GSkyRegionRect*>(&reg);
656 652  
657 653 // Dirty kludge: compare with map
658   - GSkyRegionMap regmap = GSkyRegionMap(regrect);
659   - overlap = regmap.overlaps(*this);
  654 + GSkyRegionMap map = GSkyRegionMap(rect);
  655 + overlap = map.overlaps(*this);
660 656  
661 657 } // endif: region was of type "Rect"
662 658  
... ... @@ -664,11 +660,10 @@ bool GSkyRegionRect::overlaps(const GSkyRegion&amp; reg) const
664 660 else if (reg.type() == "Map") {
665 661  
666 662 // Create map from reg
667   - const GSkyRegionMap* regmap =
668   - dynamic_cast<const GSkyRegionMap*>(&reg);
  663 + const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(&reg);
669 664  
670 665 // Check overlap with circle
671   - overlap = regmap->overlaps(*this);
  666 + overlap = map->overlaps(*this);
672 667  
673 668 } // endif: region was map
674 669  
... ... @@ -684,99 +679,47 @@ bool GSkyRegionRect::overlaps(const GSkyRegion&amp; reg) const
684 679  
685 680  
686 681 /***********************************************************************//**
687   - * @brief Transform a sky direction to the local cartesian coordinate system.
688   - *
689   - * @param[in] skydir Sky direction in global coordinate system.
690   - * @return Sky direction in local cartesian region object coordinates.
691   - *
692   - * Transform the sky direction @p skydir to the local cartesian coordinate
693   - * system. The origin of the local coordinate system is fixed to the center
694   - * of the rectangle and aligned in +ra (width) and +dec (height). Please
695   - * note, that a GSkyDir object is returned, even if the local coordsys is
696   - * is cartesian.
697   - ***************************************************************************/
698   -GSkyDir GSkyRegionRect::transform_to_local(const GSkyDir& skydir) const
699   -{
700   - // Get distance and polar angle relative to rectangle centre
701   - double dist = m_centre.dist(skydir);
702   - double pa = m_centre.posang(skydir) - m_posang * gammalib::deg2rad;
703   -
704   - // Compute corresponding x and y coordinate from polar coords
705   - double x = dist * std::sin(pa);
706   - double y = dist * std::cos(pa);
707   -
708   - // Create output sky direction
709   - GSkyDir transformed;
710   - transformed.radec(x, y);
711   -
712   - // Return
713   - return transformed;
714   -}
715   -
716   -
717   -/***********************************************************************//**
718   - * @brief Transform a local coordinate to the global coordinate system.
719   - *
720   - * @param[in] locdir Direction in local coordinate system.
721   - * @return Sky direction in global coordinate system.
722   - *
723   - * Transform the local direction @p locdir to the global coordinate system.
724   - * The origin of the local coordinate system is fixed to the center of the
725   - * rectangle and aligned in +ra (width) and +dec (height).
726   - ***************************************************************************/
727   -GSkyDir GSkyRegionRect::transform_to_global(const GSkyDir& locdir) const
728   -{
729   - // // Initialise centre in local coordinates (maybe make class member?)
730   - GSkyDir centre_local;
731   - centre_local.radec(0,0);
732   -
733   - // Access x and y in local coordsys
734   - double x = locdir.ra();
735   - double y = locdir.dec();
736   -
737   - // Compute dist and position angle in polar coords
738   - double dist = std::sqrt(x*x + y*y);
739   - double pa = gammalib::pihalf - std::atan2(y, x) +
740   - m_posang * gammalib::deg2rad;
741   -
742   - // Create output global sky direction
743   - GSkyDir transformed(m_centre);
744   - transformed.rotate(pa, dist);
745   -
746   - // Return
747   - return transformed;
748   -}
749   -
750   -
751   -/***********************************************************************//**
752   - * @brief Get the sky direction of a corner of the rectangle, there are four.
  682 + * @brief Return sky direction of one corner of the rectangle
753 683 *
754 684 * @param[in] index Corner index [0,3]
755   - * @return Sky direction of corner in global coordinate system.
  685 + * @return Sky direction of corner.
756 686 *
757 687 * @exception GException::out_of_range
758 688 * Corner index is not in [0,3].
759 689 *
760   - * The index counts CCW starting from +RA,+DEC: +R+D, +R-D, -R-D, -R+D.
  690 + * Returns the sky direction of one corner of the rectangle. The corner
  691 + * @p index is counted counterclockwise from the positive Right Ascension
  692 + * and Declination directions:
  693 + *
  694 + * * `index = 0` : \f$(\alpha + \frac{\rm width}{2}, \delta + \frac{\rm height}{2})\f$
  695 + * * `index = 1` : \f$(\alpha + \frac{\rm width}{2}, \delta - \frac{\rm height}{2})\f$
  696 + * * `index = 2` : \f$(\alpha - \frac{\rm width}{2}, \delta + \frac{\rm height}{2})\f$
  697 + * * `index = 3` : \f$(\alpha - \frac{\rm width}{2}, \delta - \frac{\rm height}{2})\f$
  698 + *
  699 + * where \f$\alpha\f$ and \f$\delta\f$ are the Right Ascension and
  700 + * Declination of the rectangle centre, and \f${\rm width}\f$ and
  701 + * \f${\rm height}\f$ are the width and height of the rectangle.
761 702 ***************************************************************************/
762   -GSkyDir GSkyRegionRect::get_corner(const int& index) const
  703 +GSkyDir GSkyRegionRect::corner(const int& index) const
763 704 {
764 705 // Assert index is in [0,3]
765 706 if ((index < 0) || (index > 3)) {
766   - throw GException::out_of_range(G_GET_CORNER, "Corner index", index, 3);
  707 + throw GException::out_of_range(G_CORNER, "Corner index", index, 3);
767 708 }
768 709  
769   - // Compute the offset to the corners
770   - // Index counts CCW starting from +RA,+DEC: +R+D, +R-D, -R-D, -R+D
  710 + // Compute the offset to the corners in degrees
771 711 double dx = 0.5 * m_width * (1 - 2*int(index >= 2));
772 712 double dy = 0.5 * m_height * (1 - 2*int((index==1) || (index==2)));
773 713  
774   - // Define local corner sky direction
775   - GSkyDir corner;
776   - corner.radec_deg(dx,dy);
  714 + // Convert offset into radians
  715 + dx *= gammalib::deg2rad;
  716 + dy *= gammalib::deg2rad;
  717 +
  718 + // Define corner in local coordinates
  719 + GSkyPixel local(dx,dy);
777 720  
778 721 // Transform to global sky direction
779   - corner = transform_to_global(corner);
  722 + GSkyDir corner = local_to_dir(local);
780 723  
781 724 // Return
782 725 return corner;
... ... @@ -871,7 +814,19 @@ void GSkyRegionRect::free_members(void)
871 814  
872 815  
873 816 /***********************************************************************//**
874   - * @brief Compute solid angle (sr)
  817 + * @brief Compute solid angle of rectangle
  818 + *
  819 + * This method computes the solid angle of the rectangle in steradians and
  820 + * stores the result in the `m_solid` member.
  821 + *
  822 + * The solid angle is computed using
  823 + *
  824 + * \f[
  825 + * \Omega = {\rm width} \times {\rm height}
  826 + * \f]
  827 + *
  828 + * where \f${\rm width}\f$ and \f${\rm height}\f$ are the width and the
  829 + * height of the rectangle expressed in radians.
875 830 ***************************************************************************/
876 831 void GSkyRegionRect::compute_solid_angle(void)
877 832 {
... ... @@ -881,3 +836,110 @@ void GSkyRegionRect::compute_solid_angle(void)
881 836 // Return
882 837 return;
883 838 }
  839 +
  840 +
  841 +/***********************************************************************//**
  842 + * @brief Checks if local direction lies within region
  843 + *
  844 + * @param[in] locdir Sky direction in local coordinate system.
  845 + *
  846 + * A local direction lies within a region when its distance to the region
  847 + * centre is not larger than the region extension in both axes directions.
  848 + ***************************************************************************/
  849 +bool GSkyRegionRect::contains(const GSkyPixel& local) const
  850 +{
  851 + // Initialise return value
  852 + bool contains = false;
  853 +
  854 + // Get local y coordinate in degrees
  855 + double y = local.y() * gammalib::rad2deg;
  856 +
  857 + // Continue if the y coordinate is within half of the height
  858 + if ((std::abs(y) - 1.0e-10) <= 0.5 * m_height) {
  859 +
  860 + // Get local x coordinate in degrees
  861 + double x = local.x() * gammalib::rad2deg;
  862 +
  863 + // If the x coordinate is within half of the width the local
  864 + // coordinate is contained within the rectangle
  865 + if ((std::abs(x) - 1.0e-10) <= 0.5 * m_width) {
  866 + contains = true;
  867 + }
  868 +
  869 + } // endif: y coordinate is within half of the height
  870 +
  871 + // Return containment
  872 + return contains;
  873 +}
  874 +
  875 +
  876 +/***********************************************************************//**
  877 + * @brief Transform sky direction to local rectangle coordinates
  878 + *
  879 + * @param[in] dir Sky direction.
  880 + * @return Local cartesian region object coordinates.
  881 + *
  882 + * Transform the sky direction @p dir to the local cartesian coordinate
  883 + * system. The local coordinates \f$(x,y)\f$ are defined by
  884 + *
  885 + * \f[
  886 + * x = d \sin(\phi)
  887 + * \f]
  888 + *
  889 + * and
  890 + *
  891 + * \f[
  892 + * y = d \cos(\phi)
  893 + * \f]
  894 + *
  895 + * where \f$d\f$ is the angular distance between the sky direction @p dir
  896 + * and the centre of the rectangle in radians, and \f$\phi\f$ is the
  897 + * position angle of the sky direction @p dir with respect to the position
  898 + * angle of the rectangle, computed counter clockwise.
  899 + ***************************************************************************/
  900 +GSkyPixel GSkyRegionRect::dir_to_local(const GSkyDir& dir) const
  901 +{
  902 + // Get distance and polar angle relative to rectangle centre
  903 + double dist = m_centre.dist(dir);
  904 + double pa = m_centre.posang(dir) - m_posang * gammalib::deg2rad;
  905 +
  906 + // Compute corresponding x and y coordinate from polar coords
  907 + double x = dist * std::sin(pa);
  908 + double y = dist * std::cos(pa);
  909 +
  910 + // Create local coordinates
  911 + GSkyPixel local(x,y);
  912 +
  913 + // Return local coordinates
  914 + return local;
  915 +}
  916 +
  917 +
  918 +/***********************************************************************//**
  919 + * @brief Transform local rectangle coordinates to sky direction
  920 + *
  921 + * @param[in] local Local rectangle coordinates.
  922 + * @return Sky direction.
  923 + *
  924 + * Transform the local coordinates @p local to a sky direction. See
  925 + * GSkyRegionRect::dir_to_local for the relation between local rectangle
  926 + * coordinates and the sky direction.
  927 + ***************************************************************************/
  928 +GSkyDir GSkyRegionRect::local_to_dir(const GSkyPixel& local) const
  929 +{
  930 + // Get local coordinates (radians)
  931 + double x = local.x();
  932 + double y = local.y();
  933 +
  934 + // Compute dist and position angle in polar coords
  935 + double dist = std::sqrt(x*x + y*y);
  936 + double pa = gammalib::pihalf - std::atan2(y, x) +
  937 + m_posang * gammalib::deg2rad;
  938 +
  939 + // Create output global sky direction
  940 + GSkyDir dir(m_centre);
  941 + dir.rotate(pa, dist);
  942 +
  943 + // Return sky direction
  944 + return dir;
  945 +}
... ...
test/test_GSky.cpp
  View file @ca6d6d3
... ... @@ -1983,35 +1983,37 @@ void TestGSky::test_GSkyRegionRect(void)
1983 1983 test_assert(!refregion.contains(test_dir_2),"test for containment");
1984 1984  
1985 1985 // Corners
1986   - test_assert(refregion.contains(refregion.get_corner(0)),"test for corner containment");
1987   - test_assert(refregion.contains(refregion.get_corner(1)),"test for corner containment");
1988   - test_assert(refregion.contains(refregion.get_corner(2)),"test for corner containment");
1989   - test_assert(refregion.contains(refregion.get_corner(3)),"test for corner containment");
  1986 + test_assert(refregion.contains(refregion.corner(0)),"test for corner containment");
  1987 + test_assert(refregion.contains(refregion.corner(1)),"test for corner containment");
  1988 + test_assert(refregion.contains(refregion.corner(2)),"test for corner containment");
  1989 + test_assert(refregion.contains(refregion.corner(3)),"test for corner containment");
1990 1990  
1991 1991  
1992 1992 // Check GSkyDir containment for rotated rectangle
1993 1993 GSkyRegionRect refregion_rot(refdir_radeczerozero, 10, 10, 45);
1994 1994  
1995 1995 // Check corner containment
1996   - test_assert(refregion_rot.contains(refregion_rot.get_corner(0)),"test for corner containment");
1997   - test_assert(refregion_rot.contains(refregion_rot.get_corner(1)),"test for corner containment");
1998   - test_assert(refregion_rot.contains(refregion_rot.get_corner(2)),"test for corner containment");
1999   - test_assert(refregion_rot.contains(refregion_rot.get_corner(3)),"test for corner containment");
  1996 + test_assert(refregion_rot.contains(refregion_rot.corner(0)),"test for corner containment");
  1997 + test_assert(refregion_rot.contains(refregion_rot.corner(1)),"test for corner containment");
  1998 + test_assert(refregion_rot.contains(refregion_rot.corner(2)),"test for corner containment");
  1999 + test_assert(refregion_rot.contains(refregion_rot.corner(3)),"test for corner containment");
2000 2000  
2001 2001 // Check GSkyDir containment
2002   - test_dir_0.radec_deg(refregion_rot.get_corner(0).ra_deg()+1, 0);
2003   - test_dir_1.radec_deg(refregion_rot.get_corner(0).ra_deg() , 0);
2004   - test_dir_2.radec_deg(refregion_rot.get_corner(0).ra_deg()-1, 0);
  2002 + test_dir_0.radec_deg(refregion_rot.corner(0).ra_deg()+1, 0);
  2003 + test_dir_1.radec_deg(refregion_rot.corner(0).ra_deg() , 0);
  2004 + test_dir_2.radec_deg(refregion_rot.corner(0).ra_deg()-1, 0);
2005 2005 test_assert(!refregion_rot.contains(test_dir_0),"test for containment");
2006 2006 test_assert(refregion_rot.contains(test_dir_1),"test for containment");
2007 2007 test_assert(refregion_rot.contains(test_dir_2),"test for containment");
2008 2008  
2009 2009 // Check local<->global coordinate transformations
  2010 + /*
2010 2011 test_dir_0.radec_deg(50,50);
2011 2012 test_dir_1 = refregion_rot.transform_to_local(test_dir_0);
2012 2013 test_dir_2 = refregion_rot.transform_to_global(test_dir_1);
2013 2014 test_value(test_dir_0.ra(), test_dir_2.ra(), 1.0e-10, "Test coord trafo");
2014 2015 test_value(test_dir_0.dec(), test_dir_2.dec(), 1.0e-10, "Test coord trafo");
  2016 + */
2015 2017  
2016 2018 // Check other rects containment
2017 2019 GSkyDir refdir_raoffset = GSkyDir();
... ...