GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSkyRegionMap.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSkyRegionMap.cpp - Sky region map class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2020 by Pierrick Martin *
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 GSkyRegionMap.cpp
23  * @brief Sky region map implementation
24  * @author Pierrick Martin
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GSkyDir.hpp"
34 #include "GSkyMap.hpp"
35 #include "GSkyRegion.hpp"
36 #include "GSkyRegionCircle.hpp"
37 #include "GSkyRegionRectangle.hpp"
38 #include "GSkyRegionMap.hpp"
39 
40 /* __ Method name definitions ____________________________________________ */
41 #define G_CONTAINS "GSkyRegionMap::contains(GSkyRegion&)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 /* __ Prototype __________________________________________________________ */
50 
51 
52 /*==========================================================================
53  = =
54  = Constructors/destructors =
55  = =
56  ==========================================================================*/
57 
58 /***********************************************************************//**
59  * @brief Void constructor
60  ***************************************************************************/
62 {
63  // Initialise members
64  init_members();
65 
66  // Return
67  return;
68 }
69 
70 
71 /***********************************************************************//**
72  * @brief FITS filename constructor
73  *
74  * @param[in] filename FITS file.
75  *
76  * Constructs region from a FITS file.
77  ***************************************************************************/
79 {
80  // Initialise members
81  init_members();
82 
83  // Load FITS file
84  load(filename);
85 
86  // Return
87  return;
88 }
89 
90 
91 /***********************************************************************//**
92  * @brief Sky map constructor
93  *
94  * @param[in] map Sky map.
95  ***************************************************************************/
97 {
98  // Initialise members
99  init_members();
100 
101  // Set sky map
102  this->map(map);
103 
104  // Return
105  return;
106 }
107 
108 
109 /***********************************************************************//**
110  * @brief Region conversion constructor
111  *
112  * @param[in] region Pointer to sky region.
113  *
114  * Constructs a sky region map from any region type.
115  ***************************************************************************/
117 {
118  // Initialise members
119  init_members();
120 
121  // Get dynamic casts
122  const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(region);
123  const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(region);
124  const GSkyRegionRectangle* rect = dynamic_cast<const GSkyRegionRectangle*>(region);
125 
126  // If region is of type GSkyRegionMap then copy class members
127  if (map != NULL) {
128  copy_members(*map);
129  }
130 
131  // ... otherwise, if region is of type GSkyRegionCircle then set region
132  // map from region circle
133  else if (circle != NULL) {
134  set_region_circle(circle);
135  }
136 
137  // If region is of type GSkyRegionRectangle than set region map fron
138  // region rect
139  else if (rect != NULL) {
140  set_region_rectangle(rect);
141  }
142 
143  // Return
144  return;
145 }
146 
147 
148 /***********************************************************************//**
149  * @brief Copy constructor
150  *
151  * @param[in] region Sky region map.
152  ***************************************************************************/
154 {
155  // Initialise members
156  init_members();
157 
158  // Copy members
159  copy_members(region);
160 
161  // Return
162  return;
163 }
164 
165 
166 /***********************************************************************//**
167  * @brief Destructor
168  ***************************************************************************/
170 {
171  // Free members
172  free_members();
173 
174  // Return
175  return;
176 }
177 
178 
179 /*==========================================================================
180  = =
181  = Operators =
182  = =
183  ==========================================================================*/
184 
185 /***********************************************************************//**
186  * @brief Assignment operator
187  *
188  * @param[in] region Sky region map.
189  * @return Sky region map.
190  ***************************************************************************/
192 {
193  // Execute only if object is not identical
194  if (this != &region) {
195 
196  // Free members
197  free_members();
198 
199  // Initialise private members
200  init_members();
201 
202  // Copy members
203  copy_members(region);
204 
205  } // endif: object was not identical
206 
207  // Return this object
208  return *this;
209 }
210 
211 
212 /*==========================================================================
213  = =
214  = Public methods =
215  = =
216  ==========================================================================*/
217 
218 /***********************************************************************//**
219  * @brief Clear instance
220  ***************************************************************************/
222 {
223  // Free members
224  free_members();
225  this->GSkyRegion::free_members();
226 
227  // Initialise private members
228  this->GSkyRegion::init_members();
229  init_members();
230 
231  // Return
232  return;
233 }
234 
235 
236 /***********************************************************************//**
237  * @brief Clone sky region map
238  *
239  * @return Deep copy to sky region map.
240  ***************************************************************************/
242 {
243  // Clone sky map region
244  return new GSkyRegionMap(*this);
245 }
246 
247 
248 /***********************************************************************//**
249  * @brief Set sky map
250  *
251  * @param[in] map Sky map
252  *
253  * Sets the sky map of a region map.
254  ***************************************************************************/
255 void GSkyRegionMap::map(const GSkyMap& map)
256 {
257  // Set map object
258  m_map = map;
259 
260  // Set non-zero pixel indices
262 
263  // Compute solid angle
265 
266  // Return
267  return;
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Load region map from FITS file
273  *
274  * @param[in] filename FITS file name.
275  ***************************************************************************/
276 void GSkyRegionMap::load(const GFilename& filename)
277 {
278  // Free members
279  free_members();
280 
281  // Initialise private members
282  init_members();
283 
284  // Load map
285  m_map.load(filename);
286 
287  // Set non-zero indices
289 
290  // Compute solid angle
292 
293  // Return
294  return;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Create region from a string in DS9 format
300  *
301  * @param[in] line String describing region in DS9 region format.
302  *
303  * Empty implementation of a virtual function because map is not supported
304  * as DS9 region.
305  *
306  * @todo Translate any DS9 region into a map.
307  ***************************************************************************/
308 void GSkyRegionMap::read(const std::string& line)
309 {
310  // Do nothing because map is not supported as DS9 region
311  return;
312 }
313 
314 
315 /***********************************************************************//**
316  * @brief Write string describing region in DS9 region format
317  *
318  * @return DS9 region string
319  *
320  * Returns an empty string.
321  *
322  * @todo Translate region map into a DS9 polygon.
323  ***************************************************************************/
324 std::string GSkyRegionMap::write(void) const
325 {
326  // Allocate string
327  std::string result;
328 
329  // Return string
330  return result;
331 }
332 
333 
334 /***********************************************************************//**
335  * @brief Print region description
336  *
337  * @param[in] chatter Chattiness.
338  * @return String containing region information.
339  ***************************************************************************/
340 std::string GSkyRegionMap::print(const GChatter& chatter) const
341 {
342  // Initialise result string
343  std::string result;;
344 
345  // Continue only if chatter is not silent
346  if (chatter != SILENT) {
347 
348  // Append string
349  result.append("=== GSkyRegionMap ===");
350  result.append("\n(");
351  result.append(m_map.print());
352  result.append(")");
353 
354  } // endif: chatter was not silent
355 
356  // Return result
357  return result;
358 }
359 
360 
361 /***********************************************************************//**
362  * @brief Check if a given direction is contained in this region
363  *
364  * @param[in] dir Sky direction.
365  * @return True if sky direction is in region, false otherwise.
366  *
367  * Checks if sky direction is contained in the region map.
368  ***************************************************************************/
369 bool GSkyRegionMap::contains(const GSkyDir& dir) const
370 {
371  // Initialize
372  bool contain = false;
373 
374  // If direction is within map boundaries
375  if (m_map.contains(dir)) {
376 
377  // Convert sky direction into pixel index
378  int i = m_map.dir2inx(dir);
379 
380  // Set containment flag
381  contain = (m_map(i,0) != 0.0);
382 
383  } // endif: map contains direction
384 
385  // Return value
386  return contain;
387 }
388 
389 
390 /***********************************************************************//**
391  * @brief Checks if a given region is fully contained within this region
392  *
393  * @param[in] reg Sky region.
394  * @return True if the region is contained within this region, false
395  * otherwise.
396  *
397  * @exception GException::feature_not_implemented
398  * Specified sky region type not supported by method.
399  ***************************************************************************/
400 bool GSkyRegionMap::contains(const GSkyRegion& reg) const
401 {
402  // Initialise return value
403  bool inside = true;
404 
405  // Case of a region map
406  if (reg.type() == "Map") {
407 
408  // Cast to map type
409  const GSkyRegionMap* inreg = static_cast<const GSkyRegionMap*>(&reg);
410 
411  // Retrieve map data
412  const GSkyMap& regmap = inreg->map();
413  const std::vector<int>& regnzvec = inreg->nonzero_indices();
414 
415  // Loop over input map non-zero pixels and check if they are all in
416  for (int i = 0; i < regnzvec.size(); ++i) {
417  if (!contains(regmap.inx2dir(regnzvec[i]))) {
418  inside = false;
419  break;
420  }
421  }
422 
423  } // endcase: region map
424 
425  // Case of a circular region
426  else if (reg.type() == "Circle") {
427 
428  // Create circular region from reg
429  const GSkyRegionCircle* inreg =
430  static_cast<const GSkyRegionCircle*>(&reg);
431 
432  // Retrieve circle data
433  const GSkyDir& centre = inreg->centre();
434  double rad = inreg->radius();
435 
436  // Test a number of points along the circle
437  const int numtest = 360;
438  double rotangle = 360.0 / numtest;
439  GSkyDir testdir;
440  for (int i = 0; i < numtest; ++i) {
441  testdir = centre;
442  testdir.rotate_deg(i*rotangle, rad);
443  if (!contains(testdir)) {
444  inside = false;
445  break;
446  }
447  }
448 
449  } // endcase: circular region
450 
451  // Case of a rectangular region
452  else if (reg.type() == "Rectangle") {
453 
454  // Create rectangular region from reg
455  const GSkyRegionRectangle* rect =
456  static_cast<const GSkyRegionRectangle*>(&reg);
457 
458  // Create region map from rectangle
459  GSkyRegionMap rect_map(rect);
460 
461  // Retrieve map data
462  const GSkyMap& regmap = rect_map.map();
463  const std::vector<int>& regnzvec = rect_map.nonzero_indices();
464 
465  // Loop over input map non-zero pixels and check if they are all in
466  for (int i = 0; i < regnzvec.size(); ++i) {
467  if (!contains(regmap.inx2dir(regnzvec[i]))) {
468  inside = false;
469  break;
470  }
471  }
472  } // endcase: rectangular region
473 
474  // Other cases not implemented
475  else {
477  "Method only implemented for cicular, rectangular and map regions.");
478  }
479 
480  // Return result
481  return inside;
482 }
483 
484 
485 /***********************************************************************//**
486  * @brief Checks if a given region is overlapping with this region
487  *
488  * @param[in] reg Sky region.
489  * @return True if the region is overlapping, false otherwise.
490  ***************************************************************************/
491 bool GSkyRegionMap::overlaps(const GSkyRegion& reg) const
492 {
493  // Initialise return value
494  bool overlap = false;
495 
496  // Loop over non-zero pixels until direction is contained into region
497  for (int i = 0; i < m_nonzero_indices.size(); ++i) {
499  if (reg.contains(dir)) {
500  overlap = true;
501  break;
502  }
503  }
504 
505  // Return
506  return overlap;
507 }
508 
509 
510 /*==========================================================================
511  = =
512  = Private methods =
513  = =
514  ==========================================================================*/
515 
516 /***********************************************************************//**
517  * @brief Initialise class members
518  ***************************************************************************/
520 {
521  // Initialie base class members
522  m_type = "Map";
523 
524  // Initialise members
525  m_map.clear();
526  m_nonzero_indices.clear();
527 
528  //Return
529  return;
530 }
531 
532 
533 /***********************************************************************//**
534  * @brief Copy class members
535  *
536  * @param[in] map Sky region map.
537  ***************************************************************************/
539 {
540  // Copy members
541  m_map = map.m_map;
543 
544  // Compute solid angle
546 
547  // Return
548  return;
549 }
550 
551 
552 /***********************************************************************//**
553  * @brief Delete class members
554  ***************************************************************************/
556 {
557  // No memory to be freed, return
558  return;
559 }
560 
561 
562 /***********************************************************************//**
563  * @brief Compute solid angle
564  ***************************************************************************/
566 {
567  // Initialise solid angle
568  m_solid = 0.0;
569 
570  // Loop over all non-zero pixels of map and add up solid angles
571  for (int i = 0; i < m_nonzero_indices.size(); ++i) {
573  }
574 
575  // Return
576  return;
577 }
578 
579 
580 /***********************************************************************//**
581  * @brief Create an array of non-zero pixel indices
582  ***************************************************************************/
584 {
585  // Clear zero pixel array
586  m_nonzero_indices.clear();
587 
588  // Loop over map pixels and find non-zero pixels
589  for (int i = 0; i < m_map.npix(); ++i) {
590  if (m_map(i,0) != 0.0) {
591  m_nonzero_indices.push_back(i);
592  }
593  }
594 
595  // Return
596  return;
597 }
598 
599 
600 /***********************************************************************//**
601  * @brief Set region map from region circle
602  *
603  * @param[in] circle Region circle
604  *
605  * Sets a region map from a region circle. The method allocates a sky map
606  * with a default resolution of 1/20 of the circle radius. In any case, the
607  * sky map resolution is comprised within [0.0002, 0.1] degrees. The default
608  * size of the sky map is 2.5 times the circle radius, and in any case, at
609  * least 10 sky map pixels will be allocated. The sky map will be in
610  * celestial coordinates and TAN projection.
611  ***************************************************************************/
613 {
614  // Set sky map resolution from circle radius and constrain the
615  // resolution to [0.0002, 0.1] deg
616  double resolution = 0.05 * circle->radius();
617  if (resolution < 0.0002) {
618  resolution = 0.0002;
619  }
620  else if (resolution > 0.1) {
621  resolution = 0.1;
622  }
623 
624  // Set number of sky map pixels
625  int nbins = int(2.5 * circle->radius() / resolution);
626  if (nbins < 10) {
627  nbins = 10;
628  }
629 
630  // Set sky map
631  GSkyMap map("TAN", "CEL", circle->ra(), circle->dec(),
633  nbins, nbins);
634 
635  // Set sky map pixels
636  for (int i = 0; i < map.npix(); ++i) {
637 
638  // Compute distance to circle centre
639  double distance = circle->centre().dist_deg(map.inx2dir(i));
640 
641  // If distance is smaller or equal to radius then set map
642  // pixel to 1
643  if (distance <= circle->radius()) {
644  map(i) = 1.0;
645  }
646 
647  } // endfor: looped over sky map pixels
648 
649  // Set sky map as region map
650  this->map(map);
651 
652  // Return
653  return;
654 }
655 
656 
657 /***********************************************************************//**
658  * @brief Set region map from region rectangle
659  *
660  * @param[in] rect Region rectangle
661  *
662  * Sets a region map from a region rectangle.
663  *
664  * A skymap is allocated with a resolution goal of 0.01 degrees (same for
665  * both axes). The particular resolution of either axes depends on the width
666  * over height of the RoI of the rotated rectangle. The map extent is chosen
667  * dynamically to match the RoI of the rotated rectangle. Three sparse pixels
668  * are added at either side of the map. The sky map is created in celestial
669  * coordinates and TAN projection.
670  ***************************************************************************/
672 {
673  // Define desired resolution (in degrees)
674  const double resolution_goal = 0.01;
675 
676  // Try to estimate extent of the (rotated) rectangle:
677  // Loop over the rectangle corners to get the extension
678  double dxmax = 0.0;
679  double dymax = 0.0;
680  for (int icorner = 0; icorner < 4; ++icorner) {
681 
682  // Get skydir of current corner
683  GSkyDir corner = rect->corner(icorner);
684 
685  // Compute polar props in global coords relative to rectangle centre
686  double dist = rect->centre().dist(corner);
687  double pa = rect->centre().posang(corner);
688 
689  // Compute required map extension along RA/Dec for this corner
690  double dx = dist * std::sin(pa);
691  double dy = dist * std::cos(pa);
692 
693  // Select the largest extension
694  if (dy > dymax) {
695  dymax = dy;
696  }
697  if (dx > dxmax) {
698  dxmax = dx;
699  }
700 
701  } // endfor: looped over all corners
702 
703  // Compute width of map RoI in degrees
704  double map_width = 2.0 * dxmax * gammalib::rad2deg;
705  double map_height = 2.0 * dymax * gammalib::rad2deg;
706 
707  // Compute num of bins in x direction, round
708  int nbins_x = int(map_width / resolution_goal + 0.5);
709  int nbins_y = int(map_height / resolution_goal + 0.5);
710 
711  // Ensure a minimum num of bins
712  if (nbins_x < 10) {
713  nbins_x = 10;
714  }
715  if (nbins_y < 10) {
716  nbins_y = 10;
717  }
718 
719  // Compute actual resolution using nbins
720  double resolution_x = map_width / nbins_x;
721  double resolution_y = map_height / nbins_y;
722 
723  // Add leakage pixels
724  nbins_x += 6;
725  nbins_y += 6;
726 
727  // Create fine sky map
728  GSkyMap map_fine("TAN", "CEL", rect->ra(), rect->dec(),
729  -resolution_x, resolution_y,
730  nbins_x, nbins_y);
731 
732  // Set sky map pixels that are contained in rectangle to 1
733  for (int ipix = 0; ipix < map_fine.npix(); ++ipix) {
734  if (rect->contains(map_fine.inx2dir(ipix))) {
735  map_fine(ipix) = 1.0;
736  }
737  }
738 
739  // Set sky map as region map
740  this->map(map_fine);
741 
742  // Return
743  return;
744 }
void set_nonzero_indices(void)
Create an array of non-zero pixel indices.
GSkyRegionMap(void)
Void constructor.
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
double m_solid
Solid angle subtended by the region (sr)
Definition: GSkyRegion.hpp:96
Rectangular sky region class interface definition.
void copy_members(const GSkyRegionMap &region)
Copy class members.
Sky direction class interface definition.
double dec(void) const
Return rectangular region centre Declination.
const GSkyMap & map(void) const
Return sky map.
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 compute_solid_angle(void)
Compute solid angle.
void init_members(void)
Initialise class members.
Definition: GSkyRegion.cpp:136
const std::string & type(void) const
Return region type.
Definition: GSkyRegion.hpp:137
bool contains(const GSkyDir &dir) const
Check if a given direction is contained in this region.
Gammalib tools definition.
Interface for the circular sky region class.
Interface for the rectangular sky region class.
Sky map class definition.
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
const double & radius(void) const
Return circular region radius (in degrees)
void map(const GSkyMap &map)
Set sky map.
GSkyRegionMap & operator=(const GSkyRegionMap &region)
Assignment operator.
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2787
void clear(void)
Clear instance.
const std::vector< int > & nonzero_indices(void) const
Get non-zero index vector.
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
GSkyRegionMap * clone(void) const
Clone sky region map.
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1331
virtual bool contains(const GSkyDir &dir) const =0
bool overlaps(const GSkyRegion &reg) const
Checks if a given region is overlapping with this region.
Filename class.
Definition: GFilename.hpp:62
virtual ~GSkyRegionMap(void)
Destructor.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:424
bool contains(const GSkyDir &dir) const
Checks if sky direction lies within region.
void init_members(void)
Initialise class members.
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:2023
#define G_CONTAINS
double resolution(const GModelSpatial *model)
Determine resolution of spatial model.
GChatter
Definition: GTypemaps.hpp:33
Interface for a sky region map.
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
void free_members(void)
Delete class members.
double dec(void) const
Return circular region centre Declination.
void read(const std::string &line)
Create region from a string in DS9 format.
std::vector< int > m_nonzero_indices
Vector of non-zero pixel indices.
void set_region_circle(const GSkyRegionCircle *circle)
Set region map from region circle.
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1151
const GSkyDir & centre(void) const
Return circular region centre.
double ra(void) const
Return rectangular region centre Right Ascension.
std::string write(void) const
Write string describing region in DS9 region format.
Sky region map class interface definition.
GSkyMap m_map
The region map.
void set_region_rectangle(const GSkyRegionRectangle *rect)
Set region map from region rectangle.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
std::string print(const GChatter &chatter=NORMAL) const
Print region description.
void load(const GFilename &filename)
Load region map from FITS file.
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.
const GSkyDir & centre(void) const
Return rectangular region centre.
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1858
const double rad2deg
Definition: GMath.hpp:44
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2511
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
Sky direction class.
Definition: GSkyDir.hpp:62
int dir2inx(const GSkyDir &dir) const
Returns pixel index for a given sky direction.
Definition: GSkyMap.cpp:1445
Circular sky region class interface definition.
Abstract sky region base class interface definition.
Mathematical function definitions.
double ra(void) const
Return circular region centre Right Ascension.