GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSkyRegionRectangle.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSkyRegionRectangle.cpp - Rectangular sky region class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2019-2021 by Andreas Specovius *
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 GSkyRegionRectangle.hpp
23  * @brief Rectangular sky region implementation
24  * @author Andreas Specovius
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GSkyRegionCircle.hpp"
33 #include "GSkyRegionRectangle.hpp"
34 #include "GSkyRegionMap.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_WIDTH "GSkyRegionRectangle::width(double&)"
38 #define G_HEIGHT "GSkyRegionRectangle::height(double&)"
39 #define G_READ "GSkyRegionRectangle::read(std::string&)"
40 #define G_CONTAINS "GSkyRegionRectangle::contains(GSkyRegion&)"
41 #define G_OVERLAPS "GSkyRegionRectangle::overlaps(GSkyRegion&)"
42 #define G_CORNER "GSkyRegionRectangle::corner(int&)"
43 
44 /* __ Macros _____________________________________________________________ */
45 
46 /* __ Coding definitions _________________________________________________ */
47 
48 /* __ Debug definitions __________________________________________________ */
49 
50 /* __ Prototype __________________________________________________________ */
51 
52 
53 /*==========================================================================
54  = =
55  = Constructors/destructors =
56  = =
57  ==========================================================================*/
58 
59 /***********************************************************************//**
60  * @brief Void constructor
61  ***************************************************************************/
63 {
64  // Initialise members
65  init_members();
66 
67  // Return
68  return;
69 }
70 
71 
72 /***********************************************************************//**
73  * @brief Sky direction constructor
74  *
75  * @param[in] centre Centre of rectangle.
76  * @param[in] width Region width (degrees).
77  * @param[in] height Region height (degrees).
78  * @param[in] posang Position angle (degrees).
79  ***************************************************************************/
81  const double& width,
82  const double& height,
83  const double& posang) : GSkyRegion()
84 {
85  // Initialise members
86  init_members();
87 
88  // Set members
89  this->centre(centre);
90  this->width(width);
91  this->height(height);
92  this->posang(posang);
93 
94  // Compute solid angle
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Right Ascension and Declination constructor
104  *
105  * @param[in] ra Right Ascension of region centre (degrees).
106  * @param[in] dec Declination of region centre (degrees).
107  * @param[in] width Region width (degrees).
108  * @param[in] height Region height (degrees).
109  * @param[in] posang Position angle (degrees).
110  ***************************************************************************/
112  const double& dec,
113  const double& width,
114  const double& height,
115  const double& posang) : GSkyRegion()
116 {
117  // Initialise members
118  init_members();
119 
120  // Set members
121  this->centre(ra, dec);
122  this->width(width);
123  this->height(height);
124  this->posang(posang);
125 
126  // Compute solid angle
128 
129  // Return
130  return;
131 }
132 
133 
134 /***********************************************************************//**
135  * @brief String constructor
136  *
137  * @param[in] line DS9 region file line.
138  *
139  * Constructs region from a DS9 region file line.
140  ***************************************************************************/
142 {
143  // Initialise members
144  init_members();
145 
146  // Read information from DS9 region file line
147  read(line);
148 
149  // Return
150  return;
151 
152 }
153 
154 
155 /***********************************************************************//**
156  * @brief Copy constructor
157  *
158  * @param[in] region Rectangular sky region.
159  ***************************************************************************/
161  GSkyRegion(region)
162 {
163  // Initialise members
164  init_members();
165 
166  // Copy members
167  copy_members(region);
168 
169  // Return
170  return;
171 }
172 
173 
174 /***********************************************************************//**
175  * @brief Destructor
176  ***************************************************************************/
178 {
179  // Free members
180  free_members();
181 
182  // Return
183  return;
184 }
185 
186 
187 /*==========================================================================
188  = =
189  = Operators =
190  = =
191  ==========================================================================*/
192 
193 /***********************************************************************//**
194  * @brief Assignment operator
195  *
196  * @param[in] region Rectangular sky region.
197  * @return Rectangular sky region.
198  ***************************************************************************/
200 {
201  // Execute only if object is not identical
202  if (this != &region) {
203 
204  // Copy base class members
205  this->GSkyRegion::operator=(region);
206 
207  // Free members
208  free_members();
209 
210  // Initialise private members for clean destruction
211  init_members();
212 
213  // Copy members
214  copy_members(region);
215 
216  } // endif: object was not identical
217 
218  // Return this object
219  return *this;
220 }
221 
222 
223 /*==========================================================================
224  = =
225  = Public methods =
226  = =
227  ==========================================================================*/
228 
229 /***********************************************************************//**
230  * @brief Clear instance
231  ***************************************************************************/
233 {
234  // Free members
235  free_members();
236  this->GSkyRegion::free_members();
237 
238  // Initialise private members
239  this->GSkyRegion::init_members();
240  init_members();
241 
242  // Return
243  return;
244 }
245 
246 
247 /***********************************************************************//**
248  * @brief Clone rectangular sky region
249  *
250  * @return Deep copy to rectangular sky region.
251  ***************************************************************************/
253 {
254  // Clone rectangular sky region
255  return new GSkyRegionRectangle(*this);
256 }
257 
258 
259 /***********************************************************************//**
260  * @brief Set width of rectangular region
261  *
262  * @param[in] width Rectangle width (degrees).
263  *
264  * @exception GException::invalid_argument
265  * Width value is less than 0.
266  *
267  * Sets the @p width of the rectangular sky region. Only non-negative widths
268  * are allowed.
269  ***************************************************************************/
270 void GSkyRegionRectangle::width(const double& width)
271 {
272  // Throw an exception if the width is less than zero
273  if (width < 0.0) {
274  std::string msg =
275  "A negative width of "+gammalib::str(width)+" degrees has been "
276  "specified for a rectangular sky region. Please specify a "
277  "non-negative width.";
279  }
280 
281  // Set radius value
282  m_width = width;
283 
284  // Compute the solid angle
286 
287  // Return
288  return;
289 }
290 
291 
292 /***********************************************************************//**
293  * @brief Set height of rectangular region
294  *
295  * @param[in] height Rectangle height (degrees).
296  *
297  * @exception GException::invalid_argument
298  * Height value is less than 0.
299  *
300  * Sets the @p height of the rectangular sky region. Only non-negative
301  * heights are allowed.
302  *
303  * Note that for a position angle of zero, the height axis is pointing to
304  * celestial north.
305  ***************************************************************************/
306 void GSkyRegionRectangle::height(const double& height)
307 {
308  if (height < 0.0) {
309  std::string msg =
310  "A negative height of "+gammalib::str(height)+" degrees has been "
311  "specified for a rectangular sky region. Please specify a "
312  "non-negative height.";
314  }
315 
316  // Set radius value
317  m_height = height;
318 
319  // Compute the solid angle
321 
322  // Return
323  return;
324 }
325 
326 
327 /***********************************************************************//**
328  * @brief Read region from DS9 string
329  *
330  * @param[in] line String in DS9 format.
331  *
332  * @exception GException::invalid_value
333  * Invalid value found in DS9 format string.
334  ***************************************************************************/
335 void GSkyRegionRectangle::read(const std::string& line)
336 {
337  // Clear the current instance
338  clear();
339 
340  // Split the string into 2 parts seperated by #
341  std::vector<std::string> substrings = gammalib::split(line,"#");
342  std::string region_def = (!substrings.empty()) ? substrings[0] : "";
343  std::string comment = (substrings.size() > 1) ? substrings[1] : "";
344 
345  // Finding the box
346  if (region_def.find("box") == std::string::npos) {
347  std::string msg =
348  "Unable to find the key word \"box\" in provided string "
349  "\""+line+"\". The \"box\" key word is mandatory.";
350  throw GException::invalid_value(G_READ, msg);
351  }
352 
353  // Get the the coordinate system of the values
354  std::string system = gammalib::split(region_def, ";")[0];
355 
356  // Get the substring of the important values
357  unsigned pos = region_def.find("box(");
358  unsigned end = region_def.find(")");
359  std::string regionstring = region_def.substr(pos+4, end);
360  regionstring.erase(regionstring.find(")"), 1);
361 
362  // Get the values of the region x,y,and radius
363  std::vector<std::string> values = gammalib::split(regionstring,",");
364  if (values.size() != 5) {
365  std::string msg =
366  "Invalid number of "+gammalib::str(values.size())+" arguments "
367  "after the \"box\" key word in provided string \""+line+"\". "
368  "Exactly 5 arguments are expected.";
369  throw GException::invalid_value(G_READ, msg);
370  }
371  double x = gammalib::todouble(values[0]);
372  double y = gammalib::todouble(values[1]);
373  double width = gammalib::todouble(values[2]);
374  double height = gammalib::todouble(values[3]);
375  double posang = gammalib::todouble(values[4]);
376 
377  // Get radius units
378  if (gammalib::contains(values[2], "'")) {
379  width /= 60.0;
380  }
381  else if (gammalib::contains(values[2], "\"")) {
382  width /= 3600.0;
383  }
384  if (gammalib::contains(values[3], "'")) {
385  height /= 60.0;
386  }
387  else if (gammalib::contains(values[3], "\"")) {
388  height /= 3600.0;
389  }
390  if (gammalib::contains(values[4], "'")) {
391  posang /= 60.0;
392  }
393  else if (gammalib::contains(values[4], "\"")) {
394  posang /= 3600.0;
395  }
396 
397  // Initialise centre direction
398  GSkyDir centre = GSkyDir();
399 
400  // Set the values in the correct system
401  if (system == "fk5" || system == "icrs") {
402  centre.radec_deg(x,y);
403  }
404  else if (system == "galactic") {
405  centre.lb_deg(x,y);
406  }
407  else {
408  std::string msg =
409  "Unsupported coordinate system \""+system+"\" in provided string "
410  "\""+line+"\". Only the following coordinate systems are "
411  "supported: \"fk5\", \"icrs\" and \"galactic\".";
412  throw GException::invalid_value(G_READ, msg);
413  }
414 
415  // Set members
416  this->centre(centre);
417  this->width(width);
418  this->height(height);
419  this->posang(posang);
420 
421  // Compute solid angle
423 
424  // Check if there is a given name for the region and set it
425  std::vector<std::string>comments = gammalib::split(comment, " ");
426  for (int i = 0; i < comments.size(); ++i) {
427  if (gammalib::contains(comments[i], "text")) {
428  std::vector<std::string> attributes = gammalib::split(comments[i], "=");
429  if (attributes.size() < 2) {
430  std::string msg =
431  "Invalid character sequence encountered in provided "
432  "string \""+line+"\". An attribute of the type "
433  "\"text=Name\" is expected.";
434  throw GException::invalid_value(G_READ, msg);
435  }
436  m_name = attributes[1];
437  }
438  }
439 
440  // Return
441  return;
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Write region into a string
447  *
448  * @return String to be written in a DS9 region file
449  *
450  * Writes a DS9 region into a string. The region name is only written if it
451  * is defined.
452  ***************************************************************************/
453 std::string GSkyRegionRectangle::write(void) const
454 {
455  // Allocate string
456  std::string result;
457 
458  // Set string
459  result.append("fk5;box(");
460  result.append(gammalib::str(m_centre.ra_deg()));
461  result.append(",");
462  result.append(gammalib::str(m_centre.dec_deg()));
463  result.append(",");
464  result.append(gammalib::str(m_width));
465  result.append(",");
466  result.append(gammalib::str(m_height));
467  result.append(",");
468  result.append(gammalib::str(m_posang));
469  result.append(")");
470 
471  // Optionally add region name
472  if (m_name.length() > 0) {
473  result.append(" # text=");
474  result.append(m_name);
475  }
476 
477  // Return string
478  return result;
479 }
480 
481 
482 /***********************************************************************//**
483  * @brief Checks if sky direction lies within region
484  *
485  * @param[in] dir Sky direction.
486  * @return True if sky direction is contained in region.
487  *
488  * A sky direction lies within a region when its distance to the region
489  * centre is not larger than the region extension in both axes directions.
490  ***************************************************************************/
492 {
493  // Transform sky direction to local coordinate system
494  GSkyPixel local = dir_to_local(dir);
495 
496  // Check containment using local coordinate system
497  bool contains = this->contains(local);
498 
499  // Return containment flag
500  return contains;
501 }
502 
503 
504 /***********************************************************************//**
505  * @brief Checks if region is fully contained within this region
506  *
507  * @param[in] reg Sky region.
508  * @return True if region is contained in region.
509  *
510  * @exception GException::feature_not_implemented
511  * Not all region types supported currently.
512  ***************************************************************************/
514 {
515  // Initialise return value
516  bool contains = false;
517 
518  // If the other region is a circle then check whether it is fully
519  // contained within the rectangle
520  if (reg.type() == "Circle") {
521 
522  // Create circular region from reg
523  const GSkyRegionCircle* circle =
524  static_cast<const GSkyRegionCircle*>(&reg);
525 
526  // Transform circle center to local coordinate system
527  GSkyPixel local = dir_to_local(circle->centre());
528  double x_deg = local.x() * gammalib::rad2deg;
529  double y_deg = local.y() * gammalib::rad2deg;
530 
531  // Get circle radius in degrees
532  double radius = circle->radius();
533 
534  // Check whether circle is contained within rectangle
535  if ((std::abs(y_deg) + radius) <= 0.5 * m_height) {
536  if ((std::abs(x_deg) + radius) <= 0.5 * m_width) {
537  contains = true;
538  }
539  }
540 
541  } // endif: other region was of type "Circle"
542 
543  // ... otherwise, if the other region is a rectangle then check whether
544  // all corners fall within the rectangle
545  else if (reg.type() == "Rectangle") {
546 
547  // Create rectangular region from reg
548  const GSkyRegionRectangle* rect =
549  static_cast<const GSkyRegionRectangle*>(&reg);
550 
551  // Initialise containment flag to true
552  contains = true;
553 
554  // Loop over the four corners of the rectangle :regrect:
555  for (int i = 0; i < 4; ++i) {
556 
557  // Get sky direction of current corner
558  GSkyDir corner = rect->corner(i);
559 
560  // If corner is not contained then set containment flag to false
561  // and break
562  if (!this->contains(corner)) {
563  contains = false;
564  break;
565  }
566 
567  } // endfor: looped over the four corners
568 
569  } // endif: other region was a rectangle
570 
571  // ... otherwise, if the other region is a map then check whether it
572  // is contained within the rectangle
573  else if (reg.type() == "Map") {
574 
575  // Create map from reg
576  const GSkyRegionMap* map = static_cast<const GSkyRegionMap*>(&reg);
577 
578  // Get non-zero indices
579  const std::vector<int>& indices = map->nonzero_indices();
580 
581  // Initialise containment flag to true
582  contains = true;
583 
584  // Loop over all indices
585  for (int i = 0; i < indices.size(); ++i) {
586 
587  // Get sky direction of map pixel
588  GSkyDir dir = map->map().inx2dir(indices[i]);
589 
590  // If pixel is not contained then set containment flag to false
591  // and break
592  if (!this->contains(dir)) {
593  contains = false;
594  break;
595  }
596 
597  } // endfor: looped over non-zero map indices
598 
599  } // endif: other region was of type "Map"
600 
601  // ... otherwise throw an exception
602  else {
604  "Cannot compare to region type \""+reg.type()+"\" yet.");
605  }
606 
607  // Return containment flag
608  return contains;
609 }
610 
611 
612 /***********************************************************************//**
613  * @brief Checks if region is overlapping with this region
614  *
615  * @param[in] reg Sky region.
616  * @return True if region overlaps with region.
617  *
618  * @exception GException::feature_not_implemented
619  * Regions differ in type.
620  *
621  * Checks whether a region overlaps with the rectangle.
622  *
623  * @todo Improve implementation for rectangle-rectangle
624  ***************************************************************************/
626 {
627  // Initialise return value
628  bool overlap = false;
629 
630  // If the region is a circle then check overlap between circle and
631  // rectangle using the alogrithm given under #292 at
632  // https://stackoverflow.com/questions/401847/circle-rectangle-collision-detection-intersection
633  if (reg.type() == "Circle") {
634 
635  // Create circular region from reg
636  const GSkyRegionCircle* circle =
637  static_cast<const GSkyRegionCircle*>(&reg);
638 
639  // Transform circle center to local coordinate system
640  GSkyPixel local = dir_to_local(circle->centre());
641  double x_deg = std::abs(local.x() * gammalib::rad2deg);
642  double y_deg = std::abs(local.y() * gammalib::rad2deg);
643 
644  // Get circle radius in degrees
645  double radius = circle->radius();
646 
647  // Compute half width and height
648  double half_width = 0.5 * width();
649  double half_height = 0.5 * height();
650 
651  // Assess overlap
652  if (x_deg > (half_width + radius)) {
653  overlap = false;
654  }
655  else if (y_deg > (half_height + radius)) {
656  overlap = false;
657  }
658  else if (x_deg <= half_width) {
659  overlap = true;
660  }
661  else if (y_deg <= half_height) {
662  overlap = true;
663  }
664  else {
665  double dx = x_deg - half_width;
666  double dy = y_deg - half_height;
667  double d2 = dx * dx + dy * dy;
668  if (d2 <= radius*radius) {
669  overlap = true;
670  }
671  }
672 
673  } // endif: region was of type "Circle"
674 
675  // ... otherwise if region is a rectangle then check overlap between
676  // two rectangles
677  else if (reg.type() == "Rectangle") {
678 
679  // Create rectangular region from reg
680  const GSkyRegionRectangle* rect =
681  static_cast<const GSkyRegionRectangle*>(&reg);
682 
683  // Dirty kludge: compare with map
684  GSkyRegionMap map = GSkyRegionMap(rect);
685  overlap = map.overlaps(*this);
686 
687  } // endif: region was a rectangle
688 
689  // ... otherwise if region is map then check overlap with map
690  else if (reg.type() == "Map") {
691 
692  // Create map from reg
693  const GSkyRegionMap* map = static_cast<const GSkyRegionMap*>(&reg);
694 
695  // Check overlap with circle
696  overlap = map->overlaps(*this);
697 
698  } // endif: region was map
699 
700  // ... otherwise throw an exception
701  else {
703  "Cannot compare to region type \""+reg.type()+"\" yet.");
704  }
705 
706  // Return value
707  return overlap;
708 }
709 
710 
711 /***********************************************************************//**
712  * @brief Return sky direction of one corner of the rectangle
713  *
714  * @param[in] index Corner index [0,3]
715  * @return Sky direction of corner.
716  *
717  * @exception GException::out_of_range
718  * Corner index is not in [0,3].
719  *
720  * Returns the sky direction of one corner of the rectangle. The corner
721  * @p index is counted counterclockwise from the positive Right Ascension
722  * and Declination directions:
723  *
724  * * `index = 0` : \f$(\alpha + \frac{\rm width}{2}, \delta + \frac{\rm height}{2})\f$
725  * * `index = 1` : \f$(\alpha + \frac{\rm width}{2}, \delta - \frac{\rm height}{2})\f$
726  * * `index = 2` : \f$(\alpha - \frac{\rm width}{2}, \delta + \frac{\rm height}{2})\f$
727  * * `index = 3` : \f$(\alpha - \frac{\rm width}{2}, \delta - \frac{\rm height}{2})\f$
728  *
729  * where \f$\alpha\f$ and \f$\delta\f$ are the Right Ascension and
730  * Declination of the rectangle centre, and \f${\rm width}\f$ and
731  * \f${\rm height}\f$ are the width and height of the rectangle.
732  ***************************************************************************/
733 GSkyDir GSkyRegionRectangle::corner(const int& index) const
734 {
735  // Assert index is in [0,3]
736  if ((index < 0) || (index > 3)) {
737  throw GException::out_of_range(G_CORNER, "Corner index", index, 3);
738  }
739 
740  // Compute the offset to the corners in degrees
741  double dx = 0.5 * m_width * (1 - 2*int(index >= 2));
742  double dy = 0.5 * m_height * (1 - 2*int((index==1) || (index==2)));
743 
744  // Convert offset into radians
745  dx *= gammalib::deg2rad;
746  dy *= gammalib::deg2rad;
747 
748  // Define corner in local coordinates
749  GSkyPixel local(dx,dy);
750 
751  // Transform to global sky direction
752  GSkyDir corner = local_to_dir(local);
753 
754  // Return
755  return corner;
756 }
757 
758 
759 /***********************************************************************//**
760  * @brief Print rectangular region
761  *
762  * @param[in] chatter Chattiness
763  * @return String containing region information.
764  ***************************************************************************/
765 std::string GSkyRegionRectangle::print(const GChatter& chatter) const
766 {
767  // Initialise result string
768  std::string result;;
769 
770  // Continue only if chatter is not silent
771  if (chatter != SILENT) {
772 
773  // Append header
774  result.append("=== GSkyRegionRectangle ===");
775 
776  // Append sky region information
777  result.append("\n"+gammalib::parformat("Right Ascension of centre"));
778  result.append(gammalib::str(ra())+" deg");
779  result.append("\n"+gammalib::parformat("Declination of centre"));
780  result.append(gammalib::str(dec())+" deg");
781  result.append("\n"+gammalib::parformat("Width"));
782  result.append(gammalib::str(width())+" deg");
783  result.append("\n"+gammalib::parformat("Height"));
784  result.append(gammalib::str(height())+" deg");
785  result.append("\n"+gammalib::parformat("PA"));
786  result.append(gammalib::str(posang())+" deg");
787 
788  } // endif: chatter was not silent
789 
790  // Return result
791  return result;
792 }
793 
794 
795 /*==========================================================================
796  = =
797  = Private methods =
798  = =
799  ==========================================================================*/
800 
801 /***********************************************************************//**
802  * @brief Initialise class members
803  ***************************************************************************/
805 {
806  // Initialise members
807  m_type = "Rectangle";
808  m_centre.clear();
809  m_width = 0.0;
810  m_height = 0.0;
811  m_posang = 0.0;
812 
813  //Return
814  return;
815 }
816 
817 
818 /***********************************************************************//**
819  * @brief Copy class members
820  *
821  * @param[in] region Rectangular sky region.
822  ***************************************************************************/
824 {
825  // Copy attributes
826  m_centre = region.m_centre;
827  m_width = region.m_width;
828  m_height = region.m_height;
829  m_posang = region.m_posang;
830 
831  // Return
832  return;
833 }
834 
835 
836 /***********************************************************************//**
837  * @brief Delete class members
838  ***************************************************************************/
840 {
841  // Return
842  return;
843 }
844 
845 
846 /***********************************************************************//**
847  * @brief Compute solid angle of rectangle
848  *
849  * This method computes the solid angle of the rectangle in steradians and
850  * stores the result in the `m_solid` member.
851  *
852  * The solid angle is computed using
853  *
854  * \f[
855  * \Omega = {\rm width} \times {\rm height}
856  * \f]
857  *
858  * where \f${\rm width}\f$ and \f${\rm height}\f$ are the width and the
859  * height of the rectangle expressed in radians.
860  ***************************************************************************/
862 {
863  // Compute solid angle
865 
866  // Return
867  return;
868 }
869 
870 
871 /***********************************************************************//**
872  * @brief Checks if local direction lies within region
873  *
874  * @param[in] local Sky direction in local coordinate system.
875  *
876  * A local direction lies within a region when its distance to the region
877  * centre is not larger than the region extension in both axes directions.
878  ***************************************************************************/
879 bool GSkyRegionRectangle::contains(const GSkyPixel& local) const
880 {
881  // Initialise return value
882  bool contains = false;
883 
884  // Get local y coordinate in degrees
885  double y = local.y() * gammalib::rad2deg;
886 
887  // Continue if the y coordinate is within half of the height
888  if ((std::abs(y) - 1.0e-10) <= 0.5 * m_height) {
889 
890  // Get local x coordinate in degrees
891  double x = local.x() * gammalib::rad2deg;
892 
893  // If the x coordinate is within half of the width the local
894  // coordinate is contained within the rectangle
895  if ((std::abs(x) - 1.0e-10) <= 0.5 * m_width) {
896  contains = true;
897  }
898 
899  } // endif: y coordinate is within half of the height
900 
901  // Return containment
902  return contains;
903 }
904 
905 
906 /***********************************************************************//**
907  * @brief Transform sky direction to local rectangle coordinates
908  *
909  * @param[in] dir Sky direction.
910  * @return Local cartesian region object coordinates.
911  *
912  * Transform the sky direction @p dir to the local cartesian coordinate
913  * system. The local coordinates \f$(x,y)\f$ are defined by
914  *
915  * \f[
916  * x = d \sin(\phi)
917  * \f]
918  *
919  * and
920  *
921  * \f[
922  * y = d \cos(\phi)
923  * \f]
924  *
925  * where \f$d\f$ is the angular distance between the sky direction @p dir
926  * and the centre of the rectangle in radians, and \f$\phi\f$ is the
927  * position angle of the sky direction @p dir with respect to the position
928  * angle of the rectangle, computed counter clockwise.
929  ***************************************************************************/
931 {
932  // Get distance and polar angle relative to rectangle centre
933  double dist = m_centre.dist(dir);
934  double pa = m_centre.posang(dir) - m_posang * gammalib::deg2rad;
935 
936  // Compute corresponding x and y coordinate from polar coords
937  double x = dist * std::sin(pa);
938  double y = dist * std::cos(pa);
939 
940  // Create local coordinates
941  GSkyPixel local(x,y);
942 
943  // Return local coordinates
944  return local;
945 }
946 
947 
948 /***********************************************************************//**
949  * @brief Transform local rectangle coordinates to sky direction
950  *
951  * @param[in] local Local rectangle coordinates.
952  * @return Sky direction.
953  *
954  * Transform the local coordinates @p local to a sky direction. See
955  * GSkyRegionRectangle::dir_to_local for the relation between local rectangle
956  * coordinates and the sky direction.
957  ***************************************************************************/
959 {
960  // Get local coordinates (radians)
961  double x = local.x();
962  double y = local.y();
963 
964  // Compute dist and position angle in polar coords
965  double dist = std::sqrt(x*x + y*y);
966  double pa = gammalib::pihalf - std::atan2(y, x) +
968 
969  // Create output global sky direction
970  GSkyDir dir(m_centre);
971  dir.rotate(pa, dist);
972 
973  // Return sky direction
974  return dir;
975 }
bool contains(const std::string &str, const std::string &substring)
Checks if a substring is in a string.
Definition: GTools.cpp:1342
void compute_solid_angle(void)
Compute solid angle of rectangle.
GSkyPixel dir_to_local(const GSkyDir &dir) const
Transform sky direction to local rectangle coordinates.
double m_solid
Solid angle subtended by the region (sr)
Definition: GSkyRegion.hpp:96
#define G_CORNER
Rectangular sky region class interface definition.
void free_members(void)
Delete class members.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
void x(const double &x)
Set x value of sky map pixel.
Definition: GSkyPixel.hpp:204
std::string print(const GChatter &chatter=NORMAL) const
Print rectangular region.
GSkyRegionRectangle * clone(void) const
Clone rectangular sky region.
#define G_READ
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
GSkyRegionRectangle & operator=(const GSkyRegionRectangle &region)
Assignment operator.
double dec(void) const
Return rectangular region centre Declination.
void free_members(void)
Delete class members.
Definition: GSkyRegion.cpp:168
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
void init_members(void)
Initialise class members.
Definition: GSkyRegion.cpp:136
const std::string & type(void) const
Return region type.
Definition: GSkyRegion.hpp:137
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:983
GSkyDir local_to_dir(const GSkyPixel &local) const
Transform local rectangle coordinates to sky direction.
Gammalib tools definition.
void y(const double &y)
Set y value of sky map pixel.
Definition: GSkyPixel.hpp:221
Interface for the circular sky region class.
const double & posang(void) const
Return region position angle (in degrees)
Interface for the rectangular sky region class.
std::string m_type
Type of the region (circle, rectangle,...)
Definition: GSkyRegion.hpp:94
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:1118
#define G_HEIGHT
std::string write(void) const
Write region into a string.
const double & height(void) const
Return region height extension (in degrees)
GSkyDir m_centre
Centre or reference point of the region.
const double & radius(void) const
Return circular region radius (in degrees)
void map(const GSkyMap &map)
Set sky map.
Sky map pixel class.
Definition: GSkyPixel.hpp:74
const double & width(void) const
Return region width extension (in degrees)
#define G_CONTAINS
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
const std::vector< int > & nonzero_indices(void) const
Get non-zero index vector.
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
double m_width
Width of the region (degrees)
std::string m_name
Name of the region.
Definition: GSkyRegion.hpp:95
bool overlaps(const GSkyRegion &reg) const
Checks if a given region is overlapping with this region.
bool contains(const GSkyDir &dir) const
Checks if sky direction lies within region.
void rotate(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:378
GSkyRegionRectangle(void)
Void constructor.
virtual GSkyRegion & operator=(const GSkyRegion &region)
Assignment operator.
Definition: GSkyRegion.cpp:106
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
GChatter
Definition: GTypemaps.hpp:33
Interface for a sky region map.
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
void copy_members(const GSkyRegionRectangle &region)
Copy class members.
bool overlaps(const GSkyRegion &reg) const
Checks if region is overlapping with this region.
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1151
void clear(void)
Clear instance.
const GSkyDir & centre(void) const
Return circular region centre.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:164
double ra(void) const
Return rectangular region centre Right Ascension.
Sky region map class interface definition.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:278
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
GSkyDir corner(const int &index) const
Return sky direction of one corner of the rectangle.
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
double m_height
Height of the region (degrees)
const GSkyDir & centre(void) const
Return rectangular region centre.
#define G_WIDTH
double m_posang
Position angle, counterclockwise from North (degrees)
const double rad2deg
Definition: GMath.hpp:44
void init_members(void)
Initialise class members.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
Sky direction class.
Definition: GSkyDir.hpp:62
void read(const std::string &line)
Read region from DS9 string.
Circular sky region class interface definition.
virtual ~GSkyRegionRectangle(void)
Destructor.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:926
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489