GammaLib  1.7.0.dev
 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-2018 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 "GSkyRegionMap.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_CONTAINS "GSkyRegionMap::contains(GSkyRegion&)"
41 
42 /* __ Macros _____________________________________________________________ */
43 
44 /* __ Coding definitions _________________________________________________ */
45 
46 /* __ Debug definitions __________________________________________________ */
47 
48 /* __ Prototype __________________________________________________________ */
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  ***************************************************************************/
61 {
62  // Initialise members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief FITS filename constructor
72  *
73  * @param[in] filename FITS file.
74  *
75  * Constructs region from a FITS file.
76  ***************************************************************************/
78 {
79  // Initialise members
80  init_members();
81 
82  // Load FITS file
83  load(filename);
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief Sky map constructor
92  *
93  * @param[in] map Sky map.
94  ***************************************************************************/
96 {
97  // Initialise members
98  init_members();
99 
100  // Set sky map
101  this->map(map);
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Region conversion constructor
110  *
111  * @param[in] region Pointer to sky region.
112  *
113  * Constructs a sky region map from any region type.
114  ***************************************************************************/
116 {
117  // Initialise members
118  init_members();
119 
120  // Get dynamic casts
121  const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(region);
122  const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(region);
123 
124  // If region is of type GSkyRegionMap then copy class members
125  if (map != NULL) {
126  copy_members(*map);
127  }
128 
129  // ... otherwise, if region is of type GSkyRegionCircle then set region
130  // map from region circle
131  else if (circle != NULL) {
132  set_region_circle(circle);
133  }
134 
135  // Return
136  return;
137 }
138 
139 
140 /***********************************************************************//**
141  * @brief Copy constructor
142  *
143  * @param[in] region Sky region map.
144  ***************************************************************************/
146 {
147  // Initialise members
148  init_members();
149 
150  // Copy members
151  copy_members(region);
152 
153  // Return
154  return;
155 }
156 
157 
158 /***********************************************************************//**
159  * @brief Destructor
160  ***************************************************************************/
162 {
163  // Free members
164  free_members();
165 
166  // Return
167  return;
168 }
169 
170 
171 /*==========================================================================
172  = =
173  = Operators =
174  = =
175  ==========================================================================*/
176 
177 /***********************************************************************//**
178  * @brief Assignment operator
179  *
180  * @param[in] region Sky region map.
181  * @return Sky region map.
182  ***************************************************************************/
184 {
185  // Execute only if object is not identical
186  if (this != &region) {
187 
188  // Free members
189  free_members();
190 
191  // Initialise private members
192  init_members();
193 
194  // Copy members
195  copy_members(region);
196 
197  } // endif: object was not identical
198 
199  // Return this object
200  return *this;
201 }
202 
203 
204 /*==========================================================================
205  = =
206  = Public methods =
207  = =
208  ==========================================================================*/
209 
210 /***********************************************************************//**
211  * @brief Clear instance
212  ***************************************************************************/
214 {
215  // Free members
216  free_members();
217  this->GSkyRegion::free_members();
218 
219  // Initialise private members
220  this->GSkyRegion::init_members();
221  init_members();
222 
223  // Return
224  return;
225 }
226 
227 
228 /***********************************************************************//**
229  * @brief Clone sky region map
230  *
231  * @return Deep copy to sky region map.
232  ***************************************************************************/
234 {
235  // Clone sky map region
236  return new GSkyRegionMap(*this);
237 }
238 
239 
240 /***********************************************************************//**
241  * @brief Set sky map
242  *
243  * @param[in] map Sky map
244  *
245  * Sets the sky map of a region map.
246  ***************************************************************************/
247 void GSkyRegionMap::map(const GSkyMap& map)
248 {
249  // Set map object
250  m_map = map;
251 
252  // Set non-zero pixel indices
254 
255  // Compute solid angle
257 
258  // Return
259  return;
260 }
261 
262 
263 /***********************************************************************//**
264  * @brief Load region map from FITS file
265  *
266  * @param[in] filename FITS file name.
267  ***************************************************************************/
268 void GSkyRegionMap::load(const GFilename& filename)
269 {
270  // Free members
271  free_members();
272 
273  // Initialise private members
274  init_members();
275 
276  // Load map
277  m_map.load(filename);
278 
279  // Set non-zero indices
281 
282  // Compute solid angle
284 
285  // Return
286  return;
287 }
288 
289 
290 /***********************************************************************//**
291  * @brief Create region from a string in DS9 format
292  *
293  * @param[in] line String describing region in DS9 region format.
294  *
295  * Empty implementation of a virtual function because map is not supported
296  * as DS9 region.
297  *
298  * @todo Translate any DS9 region into a map.
299  ***************************************************************************/
300 void GSkyRegionMap::read(const std::string& line)
301 {
302  // Do nothing because map is not supported as DS9 region
303  return;
304 }
305 
306 
307 /***********************************************************************//**
308  * @brief Write string describing region in DS9 region format
309  *
310  * @return DS9 region string
311  *
312  * Returns an empty string.
313  *
314  * @todo Translate region map into a DS9 polygon.
315  ***************************************************************************/
316 std::string GSkyRegionMap::write(void) const
317 {
318  // Allocate string
319  std::string result;
320 
321  // Return string
322  return result;
323 }
324 
325 
326 /***********************************************************************//**
327  * @brief Print region description
328  *
329  * @param[in] chatter Chattiness.
330  * @return String containing region information.
331  ***************************************************************************/
332 std::string GSkyRegionMap::print(const GChatter& chatter) const
333 {
334  // Initialise result string
335  std::string result;;
336 
337  // Continue only if chatter is not silent
338  if (chatter != SILENT) {
339 
340  // Append string
341  result.append("=== GSkyRegionMap ===");
342  result.append("\n(");
343  result.append(m_map.print());
344  result.append(")");
345 
346  } // endif: chatter was not silent
347 
348  // Return result
349  return result;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Check if a given direction is contained in this region
355  *
356  * @param[in] dir Sky direction.
357  * @return True if sky direction is in region, false otherwise.
358  *
359  * Checks if sky direction is contained in the region map.
360  ***************************************************************************/
361 bool GSkyRegionMap::contains(const GSkyDir& dir) const
362 {
363  // Initialize
364  bool contain = false;
365 
366  // If direction is within map boundaries
367  if (m_map.contains(dir)) {
368 
369  // Convert sky direction into pixel index
370  int i = m_map.dir2inx(dir);
371 
372  // Set containment flag
373  contain = (m_map(i,0) != 0.0);
374 
375  } // endif: map contains direction
376 
377  // Return value
378  return contain;
379 }
380 
381 
382 /***********************************************************************//**
383  * @brief Checks if a given region is fully contained within this region
384  *
385  * @param[in] reg Sky region.
386  * @return True if the region is contained within this region, false
387  * otherwise.
388  *
389  * @exception GException::feature_not_implemented
390  * Specified sky region type not supported by method.
391  ***************************************************************************/
392 bool GSkyRegionMap::contains(const GSkyRegion& reg) const
393 {
394  // Initialise return value
395  bool inside = true;
396 
397  // Case of a region map
398  if (reg.type() == "Map") {
399 
400  // Cast to map type
401  const GSkyRegionMap* inreg = dynamic_cast<const GSkyRegionMap*>(&reg);
402 
403  // Retrieve map data
404  const GSkyMap& regmap = inreg->map();
405  const std::vector<int>& regnzvec = inreg->nonzero_indices();
406 
407  // Loop over input map non-zero pixels and check if they are all in
408  for (int i = 0; i < regnzvec.size(); ++i) {
409  if (!contains(regmap.inx2dir(regnzvec[i]))) {
410  inside = false;
411  break;
412  }
413  }
414 
415  } // endcase: region map
416 
417  // Case of a circular region
418  else if (reg.type() == "Circle") {
419 
420  // Create circular region from reg
421  const GSkyRegionCircle* inreg = dynamic_cast<const GSkyRegionCircle*>(&reg);
422 
423  // Retrieve circle data
424  const GSkyDir& centre = inreg->centre();
425  double rad = inreg->radius();
426 
427  // Test a number of points along the circle
428  const int numtest = 360;
429  double rotangle = 360.0 / numtest;
430  GSkyDir testdir;
431  for (int i = 0; i < numtest; ++i) {
432  testdir = centre;
433  testdir.rotate_deg(i*rotangle, rad);
434  if (!contains(testdir)) {
435  inside = false;
436  break;
437  }
438  }
439 
440  } // endcase: circular region
441 
442  // Other cases not implemented
443  else {
445  "Method only implemented for cicular and map regions.");
446  }
447 
448  // Return result
449  return inside;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Checks if a given region is overlapping with this region
455  *
456  * @param[in] reg Sky region.
457  * @return True if the region is overlapping, false otherwise.
458  ***************************************************************************/
459 bool GSkyRegionMap::overlaps(const GSkyRegion& reg) const
460 {
461  // Initialise return value
462  bool overlap = false;
463 
464  // Loop over non-zero pixels until direction is contained into region
465  for (int i = 0; i < m_nonzero_indices.size(); ++i) {
467  if (reg.contains(dir)) {
468  overlap = true;
469  break;
470  }
471  }
472 
473  // Return
474  return overlap;
475 }
476 
477 
478 /*==========================================================================
479  = =
480  = Private methods =
481  = =
482  ==========================================================================*/
483 
484 /***********************************************************************//**
485  * @brief Initialise class members
486  ***************************************************************************/
488 {
489  // Initialie base class members
490  m_type = "Map";
491 
492  // Initialise members
493  m_map.clear();
494  m_nonzero_indices.clear();
495 
496  //Return
497  return;
498 }
499 
500 
501 /***********************************************************************//**
502  * @brief Copy class members
503  *
504  * @param[in] map Sky region map.
505  ***************************************************************************/
507 {
508  // Copy members
509  m_map = map.m_map;
511 
512  // Compute solid angle
514 
515  // Return
516  return;
517 }
518 
519 
520 /***********************************************************************//**
521  * @brief Delete class members
522  ***************************************************************************/
524 {
525  // No memory to be freed, return
526  return;
527 }
528 
529 
530 /***********************************************************************//**
531  * @brief Compute solid angle
532  ***************************************************************************/
534 {
535  // Initialise solid angle
536  m_solid = 0.0;
537 
538  // Loop over all non-zero pixels of map and add up solid angles
539  for (int i = 0; i < m_nonzero_indices.size(); ++i) {
541  }
542 
543  // Return
544  return;
545 }
546 
547 
548 /***********************************************************************//**
549  * @brief Create an array of non-zero pixel indices
550  ***************************************************************************/
552 {
553  // Clear zero pixel array
554  m_nonzero_indices.clear();
555 
556  // Loop over map pixels and find non-zero pixels
557  for (int i = 0; i < m_map.npix(); ++i) {
558  if (m_map(i,0) != 0.0) {
559  m_nonzero_indices.push_back(i);
560  }
561  }
562 
563  // Return
564  return;
565 }
566 
567 
568 /***********************************************************************//**
569  * @brief Set region map from region circle
570  *
571  * @param[in] circle Region circle
572  *
573  * Sets a region map from a region circle. The method allocates a sky map
574  * with a default resolution of 1/20 of the circle radius. In any case, the
575  * sky map resolution is comprised within [0.0002, 0.1] degrees. The default
576  * size of the sky map is 2.5 times the circle radius, and in any case, at
577  * least 10 sky map pixels will be allocated. The sky map will be in
578  * celestial coordinates and TAN projection.
579  ***************************************************************************/
581 {
582  // Set sky map resolution from circle radius and constrain the
583  // resolution to [0.0002, 0.1] deg
584  double resolution = 0.05 * circle->radius();
585  if (resolution < 0.0002) {
586  resolution = 0.0002;
587  }
588  else if (resolution > 0.1) {
589  resolution = 0.1;
590  }
591 
592  // Set number of sky map pixels
593  int nbins = int(2.5 * circle->radius() / resolution);
594  if (nbins < 10) {
595  nbins = 10;
596  }
597 
598  // Set sky map
599  GSkyMap map("TAN", "CEL", circle->ra(), circle->dec(),
601  nbins, nbins);
602 
603  // Set sky map pixels
604  for (int i = 0; i < map.npix(); ++i) {
605 
606  // Compute distance to circle centre
607  double distance = circle->centre().dist_deg(map.inx2dir(i));
608 
609  // If distance is smaller or equal to radius then set map
610  // pixel to 1
611  if (distance <= circle->radius()) {
612  map(i) = 1.0;
613  }
614 
615  } // endfor: looped over sky map pixels
616 
617  // Set sky map as region map
618  this->map(map);
619 
620  // Return
621  return;
622 }
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:280
double m_solid
Solid angle subtended by the region (sr)
Definition: GSkyRegion.hpp:96
void copy_members(const GSkyRegionMap &region)
Copy class members.
Sky direction class interface definition.
const GSkyMap & map(void) const
Return sky map.
void free_members(void)
Delete class members.
Definition: GSkyRegion.cpp:168
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.
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:997
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:2557
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:1321
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:321
void init_members(void)
Initialise class members.
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:1879
#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:1090
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.
const GSkyDir & centre(void) const
Return circular region centre.
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.
std::string print(const GChatter &chatter=NORMAL) const
Print region description.
void load(const GFilename &filename)
Load region map from FITS file.
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1770
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2281
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
Sky direction class.
Definition: GSkyDir.hpp:62
int dir2inx(const GSkyDir &dir) const
Returns pixel index for a given sky direction.
Definition: GSkyMap.cpp:1435
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.