From 4485a86ff847bf6db41e451387bb12de804d145f Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 6 Jul 2017 11:10:10 +0200 Subject: [PATCH] Integrate code --- ChangeLog | 8 ++++++++ NEWS | 204 +++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- include/GSkyMap.hpp | 6 ++++-- include/GSkyRegions.hpp | 4 ++-- inst/cta/src/GCTACubeBackground.cpp | 12 +++++++----- inst/cta/src/GCTACubeEdisp.cpp | 31 +++++++++++++++++-------------- inst/cta/src/GCTACubeExposure.cpp | 31 +++++++++++++++++-------------- inst/cta/src/GCTACubePsf.cpp | 31 +++++++++++++++++-------------- inst/cta/src/GCTAEventList.cpp | 8 +++++--- pyext/GSkyMap.i | 3 ++- pyext/GSkyRegions.i | 5 +++-- src/sky/GSkyMap.cpp | 177 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------- src/sky/GSkyRegions.cpp | 162 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------- 13 files changed, 305 insertions(+), 377 deletions(-) diff --git a/ChangeLog b/ChangeLog index fdbba6f..cf1810a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2017-07-06 + + * Version 1.4.0 released + ========================= + + Add GSkyRegions::overlaps() and GSkyMap::overlaps() classes (#1980) + + 2017-06-06 * Version 1.3.0 released diff --git a/NEWS b/NEWS index 19e3dc9..096ee0a 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ -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&, std::vector&) -- 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 diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index e4bcb4c..6ac187e 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -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; diff --git a/include/GSkyRegions.hpp b/include/GSkyRegions.hpp index 406fb38..f582392 100644 --- a/include/GSkyRegions.hpp +++ b/include/GSkyRegions.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/inst/cta/src/GCTACubeBackground.cpp b/inst/cta/src/GCTACubeBackground.cpp index 9c02e46..711eab5 100644 --- a/inst/cta/src/GCTACubeBackground.cpp +++ b/inst/cta/src/GCTACubeBackground.cpp @@ -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& 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& 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()); diff --git a/inst/cta/src/GCTACubeEdisp.cpp b/inst/cta/src/GCTACubeEdisp.cpp index f748546..1e9b5c2 100644 --- a/inst/cta/src/GCTACubeEdisp.cpp +++ b/inst/cta/src/GCTACubeEdisp.cpp @@ -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& 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 @@ -915,6 +905,19 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation& 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 { diff --git a/inst/cta/src/GCTACubeExposure.cpp b/inst/cta/src/GCTACubeExposure.cpp index f17892e..37a3ca8 100644 --- a/inst/cta/src/GCTACubeExposure.cpp +++ b/inst/cta/src/GCTACubeExposure.cpp @@ -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& 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 @@ -693,6 +683,19 @@ void GCTACubeExposure::fill_cube(const GCTAObservation& 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 { diff --git a/inst/cta/src/GCTACubePsf.cpp b/inst/cta/src/GCTACubePsf.cpp index f26767b..b939381 100644 --- a/inst/cta/src/GCTACubePsf.cpp +++ b/inst/cta/src/GCTACubePsf.cpp @@ -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& 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 (obs.response()); @@ -824,6 +814,19 @@ void GCTACubePsf::fill_cube(const GCTAObservation& 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 { diff --git a/inst/cta/src/GCTAEventList.cpp b/inst/cta/src/GCTAEventList.cpp index 779e400..4a15c2f 100644 --- a/inst/cta/src/GCTAEventList.cpp +++ b/inst/cta/src/GCTAEventList.cpp @@ -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(m_events).swap(m_events); std::vector(m_columns).swap(m_columns); - + // Return return; } diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index 329f0a2..34c18a1 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/pyext/GSkyRegions.i b/pyext/GSkyRegions.i index efbb521..e101435 100644 --- a/pyext/GSkyRegions.i +++ b/pyext/GSkyRegions.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 4862e2e..d7cc156 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -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& 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& 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(®ion); + 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& dir1, const GSkyDir& 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. diff --git a/src/sky/GSkyRegions.cpp b/src/sky/GSkyRegions.cpp index d42355e..f47da12 100644 --- a/src/sky/GSkyRegions.cpp +++ b/src/sky/GSkyRegions.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& 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& 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& 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; } - -- libgit2 0.21.4