Integrate code

Jürgen Knödlseder
1 parent 1d494a96
Showing 13 changed files with 305 additions and 377 deletions
ChangeLog
NEWS
include/GSkyMap.hpp
include/GSkyRegions.hpp
inst/cta/src/GCTACubeBackground.cpp
inst/cta/src/GCTACubeEdisp.cpp
inst/cta/src/GCTACubeExposure.cpp
inst/cta/src/GCTACubePsf.cpp
inst/cta/src/GCTAEventList.cpp
pyext/GSkyMap.i
pyext/GSkyRegions.i
src/sky/GSkyMap.cpp
src/sky/GSkyRegions.cpp
+2017-07-06
+
+        * Version 1.4.0 released
+        =========================
+
+        Add GSkyRegions::overlaps() and GSkyMap::overlaps() classes (#1980)
+
+
 2017-06-06
  
         * Version 1.3.0 released
-New Features and Important Changes in GammaLib 1.3.0
+New Features and Important Changes in GammaLib 1.4.0
  
-6 June 2017
+6 July 2017
  
  
 1. Introduction
@@ -12,15 +12,7 @@ made since the last release of GammaLib.
 2. Public interface modifications
 ---------------------------------
 The following classes have been added:
-- GPhases
-- GWcsGLS
-- GWcsSFL
-- GModelSpectralSmoothBrokenPlaw
-- GModelSpatialRadialProfile
-- GModelSpatialRadialProfileGauss
-- GModelSpatialRadialProfileDMBurkert
-- GModelSpatialRadialProfileDMEinasto
-- GModelSpatialRadialProfileDMZhao
+- n.a.
  
 The following classes have been removed:
 - n.a.
@@ -29,48 +21,14 @@ The following classes have been renamed:
 - n.a.
  
 The following methods have been added:
-- GEbounds::set(GEnergies&)
-- GEbounds::emin(int&, GEnergy&)
-- GEbounds::emax(int&, GEnergy&)
-- GEnergies::set(GEbounds&)
-- GFilename::type()
-- GMatrix::operator+(double&)
-- GMatrix::operator-(double&)
-- GMatrix::operator+=(double&)
-- GMatrix::operator-=(double&)
-- GMatrixSymmetric::operator+(double&)
-- GMatrixSymmetric::operator-(double&)
-- GMatrixSymmetric::operator+=(double&)
-- GMatrixSymmetric::operator-=(double&)
-- GMatrixSparse::operator+(double&)
-- GMatrixSparse::operator-(double&)
-- GMatrixSparse::operator+=(double&)
-- GMatrixSparse::operator-=(double&)
-- GCTAResponseIrf::rspname()
-- GCTARoi::contains(GCTAInstDir&)
-- GCTAEdisp2D::fetch()
-- GCTAAeff::ebounds()
-- GCTAAeff2D::ebounds()
-- GCTAAeffArf::ebounds()
-- GCTAAeffPerfTable::ebounds()
-- GCTAEventList::append_column(GFitsTableCol&)
-- GCTAEventList::has_phase(bool&)
-- GCTAEventList::has_phase()
-- GCTAEventList::has_detxy(bool&)
-- GCTAEventList::has_detxy()
-- GCTAEventList::has_mc_id(bool&)
-- GCTAEventList::has_mc_id()
-- GCTAEventList::set_mc_id_names(std::vector<int>&, std::vector<std::string>&)
-- GCTAEventAtom::mc_id(int&)
-- GCTAEventAtom::mc_id()
+- GSkyRegion::overlaps(GSkyRegions&)
+- GSkyMap::overlaps(GSkyRegion&)
  
 The following methods have been removed:
-- GFitsTable::size()
+- n.a.
  
 The following methods have been renamed:
-- GFitsImage::size() => GFitsImage::npix()
-- GFitsHDU::size() => GFitsHDU::cards()
-- GFitsTableCol::length() => GFitsTableCol::nrows()
+- n.a.
  
 The arguments for the following methods have been changed:
 - n.a.
@@ -81,19 +39,7 @@ The return value of the following methods has been changed:
  
 3. Configuration
 -----------------
-The Python unit tests now raise exceptions instead of exiting with a non-
-zero return code so that the user is not thrown out of Python after calling
-gammalib.test() (#1983).
-
-If the code originates from the Git repository the configuration now switches
-to developer mode which implies that the Python bindings are always built.
-This solves the problem that occured for developers that could not get rid
-of existing Python wrappers, even if they modified the Python interface
-(#1985).
-
-The shared libraries needed by the Python wrappers are now copied into a
-build/gammalib directory for unit testing instead of creating a symbolic
-link since on shared VM file systems symbolic links may not work (#1987).
+None
  
  
 4. Application module
@@ -108,57 +54,17 @@ None
  
 6. FITS I/O Interface
 ---------------------
-FITS image type conversion constructors have been added to convert a FITS
-image of any type into a specific storage class (#1986).
-
-A map between header card keywords and header card pointers was added to
-GFitsHeader to speed up the access of FITS header cards (#2120).
-
-The GFitsImage::size() method has been renamed to GFitsImage::npix() since
-the size() method overloaded the GFitsHDU::size() method (#1858).
-
-The GFitsTable::size() method has been removed since the meaning of the
-method was ambiguous. It returned the number of table columns, but the 
-GFitsTable::ncols() does exactly the same thing, hence the method was
-redundant (#1857).
-
-The GFitsHDU::size() method has been renamed to GFitsHDU::cards() since the
-meaning of the size() method was not clear (#1859).
-
-The GFitsTableCol::length() method was renamed to GFitsTableCol::nrows()
-to clarify the meaning of thie method (#1860).
+None
  
  
 7. Linear algebra module
 ------------------------
-Add operators for addition and subtraction of scalars to all matrix classes
-(#1990).
+None
  
  
 8. Model module
 ---------------
-The GModelSpectralNodes class now throws an exception if the energy or
-intensity of a node is set to a non-positive value (#990).
-
-The GModelSpectralSmoothBrokenPlaw class that implements a smoothly broken
-power law spectrum has been added (#1948).
-
-The scale of the composite spatial models can now be fitted by adding the
-free_scale="1" attribute to the XML file (#1992).
-
-Add GModelSpatialRadialProfile class to handle arbitrary radial profiles.
-The method deals with pre-computation and normalization of the profiles,
-and implements the eval() and mc() methods. The following derived classes
-have been implemented:
-- GModelSpatialRadialProfileGauss: radial Gaussian profile
-- GModelSpatialRadialProfileDMBurkert: Burkert Dark Matter profile
-- GModelSpatialRadialProfileDMEinasto: Einasto Dark Matter profile
-- GModelSpatialRadialProfileDMZhao: Zhao Dark matter profile
-Please note that these models are not deeply tested so far and the code
-and interface is still subject to change (#1520).
-
-The names of spatial model parameters are now always set to RA/DEC even if
-the input parameters were GLON/GLAT (#2113).
+None
  
  
 9. Numerics module
@@ -168,39 +74,27 @@ None
  
 10. Observation module
 ----------------------
-A set() method has been added to GEbounds and GEnergies to allow converting
-energy boundaries into an energy container and vice versa. In addition,
-conversion constructors have been added (#1977).
-
-The GEbounds::set_lin() and GEbounds::set_log() method now throw exceptions
-if the minimum energy is larger than the maximum energy. In addition, the
-GEbounds::set_log() method throws an exception if either the minimum or
-maximum energy are zero. Furthermore, the methods clear the energy boundaries
-if the requested number of intervals is non-zero (#1989). 
-
-The GPhases class has been added to handle phase selections of events (#1996).
+None
  
  
 11. Optimization module
 -----------------------
-The Levenberg-Marquardt optimiser convergence criterion is now scaled with the
-step size so that a small step size does no longer lead to a stop of the
-iterative process (#2022).
+None
  
  
 12. Sky module
 --------------
-The GWcsGLS and GWcsSFL classes have been added to support the GLS and SFL
-World Coordinate System projections (#1851).
+A GSkyRegions::overlaps(GSkyRegions&) method has been added to check the
+overlap of any region in a given container with any region in another
+container (#1980).
  
-The GWcs base class now also supports CD keywords instead of CDELT keywords
-(#1951).
+A GSkyMap::overlaps(GSkyRegions&) method has been added to check the overlap
+of a sky region with a sky map (#1980).
  
  
 13. Support module
 ------------------
-A GFilename::type() method has been added that returns the type of the file
-based on the file type extension (e.g. ".fits") (#1863).
+None
  
  
 14. XML module
@@ -235,62 +129,10 @@ None
  
 20. CTA interface
 -----------------
-Correct the containment radius computation in GCTAPsf2D (#1876).
-
-The GCTAResponseIrf::rspname() method has been added that returns the name
-of the response function (attribute m_rspname).
-
-Do no longer interpret a CTA calibration database as a file path in the
-GCTAResponseIrf::read() method (#1969).
-
-Add GCTARoi::contains() method that checks whether an instrument direction
-is contained within the region of interest.
-
-Lazy loading of the energy dispersion in GCTAEdisp2D has been implemented by
-adding a fetch() method to the class. Now the energy dispersion is only loaded
-when needed, for example for function evaluation or Monte Carlo sampling (#1940).
-
-The GCTAEdisp2D::compute_ebounds_src() method now loops over the migration to
-enhance the precision in the true energy boundaries that are returned by this
-method. In addition, the validity range of the table is now stored in members
-of the class, and a non-zero energy dispersion is only returned for valid
-elements of the response table.
-
-An GCTAAeff::ebounds() method has been added to the abstract interface for
-effective area responses to return the valid energy range of the effective
-area (#1989).
-
-The GCTAModelAeffBackground::mc() could get stuck in the rejection loop if the
-effective area returned is always zero. A counter has been added that leads to
-the exit of the loop in case that 1000 subsequent zero values are returned
-(#1988).
-
-A GCTAEventList::append_column() method has been added to allow appending of
-FITS columns to an event list. The method checks whether the length of the
-FITS column is compatible with the number of events that are present in the
-event list (#1991).
-
-Add GCTAEventList::has_phase() and GCTAEventList::has_detxy() getter and setter
-methods to manipulate the event list information flags (#1982).
-
-Add GCTAEventList::has_mc_id() and GCTAEventAtom::mc_id() getter and setter methods
-to set the Monte Carlo identifier of a given event. This allows to attribute a
-Monte Carlo simulated event to a specific model (#1483).
-
-Add GCTAEventList::set_mc_id_names() method to store Monte Carlo metadata into an
-event list and to write this information into the FITS header. Furthermore, the
-information is read back when loading an event list (#1483).
-
-The loading and saving of CTA response cubes was put into an OpenMP critial zone to
-avoid FITS file pointer exceptions (#2021).
-
-The CTA background template HDUCLAS2 was changed from "BGD" to "BKG". Both strings
-are supported by GammaLib.
-
-The maximum radius of the King PSF has been increased from 0.7 to 2 degrees to cope
-with some extreme cases. Now also the maximum PSF radius is computed for a given
-energy so that if the PSF is smaller than the maximum radius of 2 degrees a smaller
-region of interest will be used for the computation (#2112).
+The GCTACubeExposure::fill_cube(), GCTACubePsf::fill_cube(),
+GCTACubeEdisp::fill_cube(), and GCTACubeBackground::fill() methods now
+exclude any observations from the filling for which the RoI does not
+overlap with the cube definition (#1980).
  
  
 21. COMPTEL interface
@@ -35,7 +35,6 @@
 #include "GSkyPixel.hpp"
 #include "GSkyProjection.hpp"
 #include "GBilinear.hpp"
-#include "GSkyRegionCircle.hpp"
  
 /* __ Forward declarations _______________________________________________ */
 class GFilename;
@@ -47,6 +46,8 @@ class GFitsImage;
 class GFitsImageDouble;
 class GMatrix;
 class GVector;
+class GSkyRegion;
+class GSkyRegionCircle;
  
  
 /***********************************************************************//**
@@ -159,9 +160,9 @@ public:
     double                  flux(const GSkyPixel& pixel, const int& map = 0) const;
     double                  solidangle(const int& index) const;
     double                  solidangle(const GSkyPixel& pixel) const;
-    bool                    overlaps(const GSkyRegionCircle& reg) const;
     bool                    contains(const GSkyDir& dir) const;
     bool                    contains(const GSkyPixel& pixel) const;
+    bool                    overlaps(const GSkyRegion& region) const;
     const GSkyProjection*   projection(void) const;
     void                    projection(const GSkyProjection& proj);
     const double*           pixels(void) const;
@@ -196,6 +197,7 @@ private:
                                  const GSkyDir& dir3) const;
     double            solidangle(const GSkyDir& dir1, const GSkyDir& dir2,
                                  const GSkyDir& dir3, const GSkyDir& dir4) const;
+    bool              overlaps_circle(const GSkyRegionCircle& region) const;
     bool              is_healpix(const GFitsHDU& hdu) const;
     bool              is_wcs(const GFitsHDU& hdu) const;
  
 /***************************************************************************
  *                  GSkyRegions.hpp - Sky region container class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2016 by Pierrick Martin                             *
+ *  copyright (C) 2013-2017 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -83,7 +83,7 @@ public:
     void              reserve(const int& num);
     void              extend(const GSkyRegions& regions);
     bool              contains(const GSkyDir& dir) const;
-    bool              overlaps(const GSkyRegion& reg) const;
+    bool              overlaps(const GSkyRegion& region) const;
     bool              overlaps(const GSkyRegions& regions) const;
     void              load(const GFilename& filename);
     void              save(const GFilename& filename) const;
@@ -32,15 +32,16 @@
 #include "GException.hpp"
 #include "GMath.hpp"
 #include "GRan.hpp"
+#include "GLog.hpp"
 #include "GFits.hpp"
 #include "GFitsBinTable.hpp"
 #include "GObservations.hpp"
+#include "GSkyRegionCircle.hpp"
 #include "GCTAEventCube.hpp"
 #include "GCTAObservation.hpp"
 #include "GCTARoi.hpp"
 #include "GCTAInstDir.hpp"
 #include "GCTACubeBackground.hpp"
-#include "GLog.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_READ                             "GCTACubeBackground::read(GFits&)"
@@ -405,16 +406,17 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
         GSkyRegionCircle obs_reg(roi.centre().dir(), roi.radius());
         if (!m_cube.overlaps(obs_reg)) {
             if (log != NULL) {
-                *log << "Skipping ";
+                *log << "Skipping unbinned ";
                 *log << cta->instrument();
                 *log << " observation ";
                 *log << "\"" << cta->name() << "\"";
-                *log << " (id=" << cta->id() << ") as it lies outside the cube";
+                *log << " (id=" << cta->id() << ") since it does not overlap ";
+                *log << "with the background cube.";
                 *log << std::endl;
             }
             continue;
         } // endif: m_cube overlaps observation
-        
+
         // Announce observation usage
         if (log != NULL) {
             *log << "Including ";
@@ -423,7 +425,7 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
             *log << "\" (id=" << cta->id() << ")";
             *log << " in background cube computation." << std::endl;
         }
-        
+
         // Set GTI of actual observations as the GTI of the event cube
         eventcube.gti(cta->gti());
  
@@ -35,6 +35,7 @@
 #include "GFitsImage.hpp"
 #include "GFitsTable.hpp"
 #include "GObservations.hpp"
+#include "GSkyRegionCircle.hpp"
 #include "GCTACubeEdisp.hpp"
 #include "GCTAObservation.hpp"
 #include "GCTAResponseIrf.hpp"
@@ -866,20 +867,9 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
                               "to specify an RoI for this observation.";
             throw GException::invalid_value(G_FILL_CUBE, msg);
         }
-        
-        // Make sure the observation falls within the bounds of the cube
-        GSkyRegionCircle obs_reg(roi.centre().dir(), roi.radius());
-        if (!m_cube.overlaps(obs_reg)) {
-            if (log != NULL) {
-                *log << "Skipping unbinned ";
-                *log << obs.instrument();
-                *log << " observation ";
-                *log << "\"" << obs.name() << "\"";
-                *log << " (id=" << obs.id() << ") which does not overlap Edisp cube.";
-                *log << std::endl;
-            }
-            return;
-        }
+
+        // Convert RoI into a circular region for overlap checking
+        GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
  
         // Extract response from observation
         const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
@@ -915,6 +905,19 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
             }
         }
  
+        // Skip observation if observation is outside the bounds of the cube
+        else if (!m_cube.overlaps(roi_reg)) {
+            if (log != NULL) {
+                *log << "Skipping unbinned ";
+                *log << obs.instrument();
+                *log << " observation ";
+                *log << "\"" << obs.name() << "\"";
+                *log << " (id=" << obs.id() << ") since it does not overlap ";
+                *log << "with the energy dispersion cube.";
+                *log << std::endl;
+            }
+        }
+
         // ... otherwise continue
         else {
  
@@ -32,6 +32,7 @@
 #include "GMath.hpp"
 #include "GLog.hpp"
 #include "GObservations.hpp"
+#include "GSkyRegionCircle.hpp"
 #include "GCTACubeExposure.hpp"
 #include "GCTAObservation.hpp"
 #include "GCTAResponseIrf.hpp"
@@ -655,20 +656,9 @@ void GCTACubeExposure::fill_cube(const GCTAObservation&amp; obs, GLog* log)
                               "to specify an RoI for this observation.";
             throw GException::invalid_value(G_FILL_CUBE, msg);
         }
-        
-        // Check that the observation overlaps with the cube
-        GSkyRegionCircle obs_reg(roi.centre().dir(), roi.radius());
-        if (!m_cube.overlaps(obs_reg)) {
-            if (log != NULL) {
-                *log << "Skipping unbinned ";
-                *log << obs.instrument();
-                *log << " observation ";
-                *log << "\"" << obs.name() << "\"";
-                *log << " (id=" << obs.id() << ") ";
-                *log << "which lies outside exposure cube";
-                *log << std::endl;
-            }
-        }
+
+        // Convert RoI into a circular region for overlap checking
+        GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
  
         // Extract response from observation
         const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
@@ -693,6 +683,19 @@ void GCTACubeExposure::fill_cube(const GCTAObservation&amp; obs, GLog* log)
             }
         }
  
+        // Skip observation if observation is outside the bounds of the cube
+        else if (!m_cube.overlaps(roi_reg)) {
+            if (log != NULL) {
+                *log << "Skipping unbinned ";
+                *log << obs.instrument();
+                *log << " observation ";
+                *log << "\"" << obs.name() << "\"";
+                *log << " (id=" << obs.id() << ") since it does not overlap ";
+                *log << "with the exposure cube.";
+                *log << std::endl;
+            }
+        }
+
         // ... otherwise continue
         else {
  
@@ -32,6 +32,7 @@
 #include "GMath.hpp"
 #include "GLog.hpp"
 #include "GObservations.hpp"
+#include "GSkyRegionCircle.hpp"
 #include "GCTACubePsf.hpp"
 #include "GCTAObservation.hpp"
 #include "GCTAResponseIrf.hpp"
@@ -787,20 +788,9 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
             throw GException::invalid_value(G_FILL_CUBE, msg);
         }
  
-        // Make sure the observation falls within the bounds of the cube
-        GSkyRegionCircle obs_reg(roi.centre().dir(), roi.radius());
-        if (!m_cube.overlaps(obs_reg)) {
-            if (log != NULL) {
-                *log << "Skipping unbinned ";
-                *log << obs.instrument();
-                *log << " observation ";
-                *log << "\"" << obs.name() << "\"";
-                *log << " (id=" << obs.id() << ") which does not overlap PSF cube.";
-                *log << std::endl;
-            }
-            return;
-        }
-        
+        // Convert RoI into a circular region for overlap checking
+        GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
+
         // Extract response from observation
         const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
                                      (obs.response());
@@ -824,6 +814,19 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
             }
         }
  
+        // Skip observation if observation is outside the bounds of the cube
+        else if (!m_cube.overlaps(roi_reg)) {
+            if (log != NULL) {
+                *log << "Skipping unbinned ";
+                *log << obs.instrument();
+                *log << " observation ";
+                *log << "\"" << obs.name() << "\"";
+                *log << " (id=" << obs.id() << ") since it does not overlap ";
+                *log << "with the point spread function cube.";
+                *log << std::endl;
+            }
+        }
+
         // ... otherwise continue
         else {
  
@@ -797,11 +797,13 @@ void GCTAEventList::dispose(void) const
     // Clear events and additional columns
     m_events.clear();
     m_columns.clear();
-    
-    // Now actually remove the memory reserved for these objects
+
+    // Now actually remove the memory reserved for these objects. See
+    // https://prateekvjoshi.com/2013/10/20/c-vector-memory-release/
+    // for an explanation of this code.
     std::vector<GCTAEventAtom>(m_events).swap(m_events);
     std::vector<GFitsTableCol*>(m_columns).swap(m_columns);
-    
+
     // Return
     return;
 }
 /***************************************************************************
  *                         GSkyMap.i - Sky map class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -171,6 +171,7 @@ public:
     double                  solidangle(const GSkyPixel& pixel) const;
     bool                    contains(const GSkyDir& dir) const;
     bool                    contains(const GSkyPixel& pixel) const;
+    bool                    overlaps(const GSkyRegion& region) const;
     const GSkyProjection*   projection(void) const;
     void                    projection(const GSkyProjection& proj);
     const double*           pixels(void) const;
 /***************************************************************************
  *                  GSkyRegions.i - Sky region container class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2016 by Pierrick Martin                             *
+ *  copyright (C) 2013-2017 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -56,7 +56,8 @@ public:
     void              reserve(const int& num);
     void              extend(const GSkyRegions& regions);
     bool              contains(const GSkyDir& dir) const;
-    bool              overlaps(const GSkyRegion& reg) const;
+    bool              overlaps(const GSkyRegion& region) const;
+    bool              overlaps(const GSkyRegions& regions) const;
     void              load(const GFilename& filename);
     void              save(const GFilename& filename) const;
     const GFilename&  filename(void) const;
@@ -44,6 +44,8 @@
 #include "GMatrix.hpp"
 #include "GVector.hpp"
 #include "GVOClient.hpp"
+#include "GSkyRegion.hpp"
+#include "GSkyRegionCircle.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_CONSTRUCT_HPX                "GSkyMap::GSkyMap(std::string&, int&,"\
@@ -68,6 +70,7 @@
 #define G_FLUX2                                   "GSkyMap::flux(GSkyPixel&)"
 #define G_SOLIDANGLE1                             "GSkyMap::solidangle(int&)"
 #define G_SOLIDANGLE2                       "GSkyMap::solidangle(GSkyPixel&)"
+#define G_OVERLAPS                           "GSkyMap::overlaps(GSkyRegion&)"
 #define G_EXTRACT                              "GSkyMap::extract(int&, int&)"
 #define G_READ                                     "GSkyMap::read(GFitsHDU&)"
 #define G_SET_WCS     "GSkyMap::set_wcs(std::string&, std::string&, double&,"\
@@ -1716,72 +1719,11 @@ void GSkyMap::projection(const GSkyProjection&amp; proj)
 }
  
  
-/***********************************************************************
- * @brief Checks whether or not a circular region overlaps with this map
- *
- * @param[in] reg Circular region
- * @return True or False
- *
- * The check is done by first testing whether the central pointing
- * position of the observation falls into the sky map. If this is false, 
- * then pixel positions that border the sky map are tested for whether or
- * not they fall into the observation region. The positions tested can be
- * visualized as follows, where '*' marks the positions tested.
- *
- *     *   *   *   *   *   *
- *       +---+---+---+---+
- *     * |0,2|1,2|2,2|3,2| *
- *       +---+---+---+---+
- *     * |0,1|1,1|2,1|3,1| *
- *       +---+---+---+---+
- *     * |0,0|1,0|2,0|3,0| *
- *       +---+---+---+---+
- *     *   *   *   *   *   *
- *
- ***************************************************************************/
-bool GSkyMap::overlaps(const GSkyRegionCircle& reg) const
-{
-    // Check if the center of the region is inside the map
-    if (contains(reg.centre())) {
-        return true;
-    } else {
-        // Loop through each of the outer bins and check if the map
-        // overlaps with them (within some delta)
-        GSkyPixel pix1(0.0,0.0);
-        GSkyPixel pix2(0.0,0.0);
-        
-        // Check the bottom and top bins
-        pix1.y(-1);
-        pix2.y(ny());
-        for (int xbin=-1; xbin <= nx(); xbin+=1.0) {
-            pix1.x(xbin);
-            pix2.x(xbin);
-            if ((reg.contains(pix2dir(pix1))) || (reg.contains(pix2dir(pix2)))) {
-                return true;
-            }
-        } // endfor: loop on xbin
-        
-        // Check the left and right bins
-        // Check the bottom and top bins
-        pix1.x(-1);
-        pix2.x(nx());
-        for (double ybin=0; ybin < ny(); ybin+=1.0) {
-            pix1.y(ybin);
-            pix2.y(ybin);
-            if ((reg.contains(pix2dir(pix1))) || (reg.contains(pix2dir(pix2)))) {
-                return true;
-            }
-        } // endfor: loop on ybin
-
-    }
-    return false;
-}
-
-
 /***********************************************************************//**
- * @brief Verifies if sky direction falls in map
+ * @brief Checks if sky direction falls in map
  *
  * @param[in] dir Sky direction.
+ * @return True if sky direction falls into map, false otherwise.
  *
  * This method checks if the specified sky direction falls within the pixels
  * covered by the skymap. The method uses the dir2xy method to convert the
@@ -1837,6 +1779,38 @@ bool GSkyMap::contains(const GSkyPixel&amp; pixel) const
  
  
 /***********************************************************************//**
+ * @brief Checks whether a region overlaps with this map
+ *
+ * @param[in] region Region
+ * @return True if region overlaps with map, false otherwise
+ *
+ * @exception GException::feature_not_implemented
+ *            Region is not a circular sky region
+ ***************************************************************************/
+bool GSkyMap::overlaps(const GSkyRegion& region) const
+{
+    // Initialise overlap
+    bool overlap = false;
+
+    // If the region is a circular region than test on overlap
+    if (region.type() == "Circle") {
+        const GSkyRegionCircle* circle =
+              static_cast<const GSkyRegionCircle*>(&region);
+        overlap = overlaps_circle(*circle);
+    }
+
+    // ... otherwise signal that method is not implemented
+    else {
+        std::string msg = "Method not implemented for non-circular regions.";
+        throw GException::feature_not_implemented(G_OVERLAPS, msg);
+    }
+
+    // Return overlap flag
+    return overlap;
+}
+
+
+/***********************************************************************//**
  * @brief Extract maps into a new sky map object
  *
  * @param[in] map First map to extract
@@ -3033,6 +3007,83 @@ double GSkyMap::solidangle(const GSkyDir&amp; dir1, const GSkyDir&amp; dir2,
  
  
 /***********************************************************************//**
+ * @brief Checks whether a circular region overlaps with this map
+ *
+ * @param[in] region Circular region
+ * @return True if circular region overlaps with map, false otherwise
+ *
+ * The check is done by first testing whether the central pointing
+ * position of the observation falls into the sky map. If this is false, 
+ * then pixel positions that border the sky map are tested for whether or
+ * not they fall into the observation region. The positions tested can be
+ * visualized as follows, where '*' marks the positions tested.
+ *
+ *     *   *   *   *   *   *
+ *       +---+---+---+---+
+ *     * |0,2|1,2|2,2|3,2| *
+ *       +---+---+---+---+
+ *     * |0,1|1,1|2,1|3,1| *
+ *       +---+---+---+---+
+ *     * |0,0|1,0|2,0|3,0| *
+ *       +---+---+---+---+
+ *     *   *   *   *   *   *
+ *
+ ***************************************************************************/
+bool GSkyMap::overlaps_circle(const GSkyRegionCircle& region) const
+{
+    // Initialise overlap with true
+    bool overlap = false;
+
+    // If the centre of the region is inside the map then signal overlap ...
+    if (contains(region.centre())) {
+        overlap = true;
+    }
+
+    // ... otherwise loop through each of the outer bins and check if the
+    // region overlaps with them (within some delta)
+    else {
+
+        // Initialise test pixels
+        GSkyPixel pix1(0.0,0.0);
+        GSkyPixel pix2(0.0,0.0);
+
+        // Check the bottom and top bins
+        pix1.y(-1.0);
+        pix2.y(double(ny()));
+        for (double xbin = -1.0; xbin <= nx(); xbin += 1.0) {
+            pix1.x(xbin);
+            pix2.x(xbin);
+            if ((region.contains(pix2dir(pix1))) ||
+                (region.contains(pix2dir(pix2)))) {
+                overlap = true;
+                break;
+            }
+        } // endfor: loop on xbin
+
+        // If no overlap has been found then check now the left and right
+        // bins
+        if (!overlap) {
+            pix1.x(-1.0);
+            pix2.x(double(nx()));
+            for (double ybin = 0.0; ybin < ny(); ybin += 1.0) {
+                pix1.y(ybin);
+                pix2.y(ybin);
+                if ((region.contains(pix2dir(pix1))) ||
+                    (region.contains(pix2dir(pix2)))) {
+                    overlap = true;
+                    break;
+                }
+            } // endfor: loop on ybin
+        }
+
+    } // endelse: check border bins
+
+    // Return overlap flag
+    return overlap;
+}
+
+
+/***********************************************************************//**
  * @brief Check if HDU contains HEALPix data
  *
  * @param[in] hdu FITS Header Data Unit.
 /***************************************************************************
  *                 GSkyRegions.cpp - Sky region container class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2016 by Pierrick Martin                             *
+ *  copyright (C) 2013-2017 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -408,6 +408,90 @@ void GSkyRegions::extend(const GSkyRegions&amp; regions)
  
  
 /***********************************************************************//**
+ * @brief Check if direction is contained in one of the regions
+ *
+ * @param[in] dir Sky direction.
+ * @return True if the sky direction is container in one of the regions,
+ *         falls otherwise.
+ *
+ * Checks if a sky direction is contained in one of the regions of the
+ * container.
+ ***************************************************************************/
+bool GSkyRegions::contains(const GSkyDir& dir) const
+{
+	 // Initialise return value
+	 bool overlap = false;
+	 
+	 // Loop over regions
+	 for (int i = 0; i < size(); ++i) {
+		 overlap = m_regions[i]->contains(dir);
+		 if (overlap) {
+            break;
+        }
+	 }
+	 
+	 // Return result
+	 return overlap;
+}
+
+ 
+/***********************************************************************//**
+ * @brief Check if region overlaps one of the regions
+ *
+ * @param[in] region Sky region.
+ * @return True if @p region overlaps with one of the regions in the region
+ *         container, false otherwise
+ *
+ * Tells if region overlaps one of the regions
+ ***************************************************************************/
+bool GSkyRegions::overlaps(const GSkyRegion& region) const
+{
+    // Initialise return value
+    bool overlap = false;
+    
+    // Loop over regions
+    for (int i = 0; i < size(); ++i) {
+        overlap = m_regions[i]->overlaps(region);
+        if (overlap) {
+            break;
+        }
+    }
+    
+    // Return result
+    return overlap;
+}
+
+
+/***********************************************************************//**
+ * @brief Check if any of the regions in two containers overlap
+ *
+ * @param[in] regions Sky region container.
+ * @return True if one of the regions in the @p regions container overlaps
+ *         with one of the regions in the container, false otherwise
+ *
+ * Tells if all regions in  overlaps one of the regions. Note, this method
+ * returns true if ANY of the regions in the two containers overlap with 
+ * each other
+ ***************************************************************************/
+bool GSkyRegions::overlaps(const GSkyRegions& regions) const
+{
+    // Initialize return value
+    bool overlap = false;
+    
+    // Loop over each region in the container
+    for (int i = 0; i < size(); ++i) {
+        overlap = regions.overlaps(*m_regions[i]);
+        if (overlap) {
+            break;
+        }
+    }
+    
+    // Return result
+    return overlap;
+}
+
+
+/***********************************************************************//**
  * @brief Load regions from DS9 region file
  *
  * @param[in] filename DS9 region filename.
@@ -545,7 +629,7 @@ void GSkyRegions::save(const GFilename&amp; filename) const
 /***********************************************************************//**
  * @brief Print regions
  *
- * @param[in] chatter Chattiness
+ * @param[in] chatter Chattiness.
  * @return String containing region container information.
  *
  * Prints all regions into a string.
@@ -577,79 +661,6 @@ std::string GSkyRegions::print(const GChatter&amp; chatter) const
 }
  
  
-/***********************************************************************//**
- * @brief Tells if direction is contained in one of the regions
- *
- * @param[in] dir A sky direction
- * @return True or False
- *
- * Tells if direction is contained in one of the regions
- ***************************************************************************/
- bool GSkyRegions::contains(const GSkyDir& dir) const
- {
-	 // Initialise return value
-	 bool dir_is_in = false;
-	 
-	 // Loop over regions
-	 for (int i = 0; i < size(); ++i) {
-		 dir_is_in=m_regions[i]->contains(dir);
-		 if (dir_is_in) break;
-	 }
-	 
-	 // Return result
-	 return dir_is_in;
- }
-
- 
-/***********************************************************************//**
- * @brief Tells if region overlaps one of the regions
- *
- * @param[in] reg A sky region
- * @return True or False
- *
- * Tells if region overlaps one of the regions
- ***************************************************************************/
-bool GSkyRegions::overlaps(const GSkyRegion& reg) const
-{
-    // Initialise return value
-    bool reg_is_in = false;
-    
-    // Loop over regions
-    for (int i = 0; i < size(); ++i) {
-        reg_is_in=m_regions[i]->overlaps(reg);
-        if (reg_is_in) break;
-    }
-    
-    // Return result
-    return reg_is_in;
-}
-
-
-/***********************************************************************
- * @brief Tells if any of the regions in two containers overlap
- *
- * @param[in] regions GSkyRegions object containing a list of sky regions
- * @return True or False
- *
- * Tells if all regions in  overlaps one of the regions. Note, this method
- * returns true if ANY of the regions in the two containers overlap with 
- * each other
- ***************************************************************************/
-bool GSkyRegions::overlaps(const GSkyRegions& regions) const
-{
-    // Initialize return value
-    bool reg_is_in = false;
-    
-    // Loop over each region in the container
-    for (int i=0; i < size(); i++) {
-        reg_is_in = regions.overlaps( *m_regions[i] );
-        if (reg_is_in) break;
-    }
-    
-    // Return result
-    return reg_is_in;
-}
-
 /*==========================================================================
  =                                                                         =
  =                              Private methods                            =
@@ -712,4 +723,3 @@ void GSkyRegions::free_members(void)
     // Return
     return;
 }
-