Adjust sky rectangle code (#3248)

Jürgen Knödlseder
1 parent bed0350a
Showing 5 changed files with 99 additions and 83 deletions
include/GSkyRegionRect.hpp
src/sky/GSkyRegionCircle.cpp
src/sky/GSkyRegionRect.cpp
src/sky/GSkyRegions.cpp
test/test_GSky.cpp
@@ -30,10 +30,8 @@
 /* __ Includes ___________________________________________________________ */
 #include <vector>
 #include <string>
-//#include "GMath.hpp"
 #include "GSkyDir.hpp"
 #include "GSkyRegion.hpp"
-//#include "GTools.hpp"
 /* __ Forward declarations _______________________________________________ */
@@ -28,9 +28,10 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
+#include "GTools.hpp"
 #include "GSkyRegionCircle.hpp"
 #include "GSkyRegionRect.hpp"
-#include "GTools.hpp"
+#include "GSkyRegionMap.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_RADIUS                          "GSkyRegionCircle::radius(double&)"
@@ -70,10 +71,10 @@ GSkyRegionCircle::GSkyRegionCircle(void) : GSkyRegion()
  * @brief Direction constructor
  *
  * @param[in] centre Centre sky direction.
- * @param[in] radius Region radius [deg].
+ * @param[in] radius Region radius (degrees).
  ***************************************************************************/
-GSkyRegionCircle::GSkyRegionCircle(const GSkyDir& centre, const double& radius) :
-                  GSkyRegion()
+GSkyRegionCircle::GSkyRegionCircle(const GSkyDir& centre,
+                                   const double& radius) : GSkyRegion()
 {
     // Initialise members
     init_members();
@@ -93,11 +94,12 @@ GSkyRegionCircle::GSkyRegionCircle(const GSkyDir&amp; centre, const double&amp; radius) 
 /***********************************************************************//**
  * @brief Direction constructor
  *
- * @param[in] ra Right Ascension of region centre [deg].
- * @param[in] dec Declination of region centre [deg].
- * @param[in] radius Region radius [deg].
+ * @param[in] ra Right Ascension of region centre (degrees).
+ * @param[in] dec Declination of region centre (degrees).
+ * @param[in] radius Region radius (degrees).
  ***************************************************************************/
-GSkyRegionCircle::GSkyRegionCircle(const double& ra, const double& dec,
+GSkyRegionCircle::GSkyRegionCircle(const double& ra,
+                                   const double& dec,
                                    const double& radius) : GSkyRegion()
 {
     // Initialise members
@@ -243,7 +245,7 @@ GSkyRegionCircle* GSkyRegionCircle::clone(void) const
 /***********************************************************************//**
  * @brief Set radius of circular region
  *
- * @param[in] radius Radius [deg].
+ * @param[in] radius Radius (degrees).
  *
  * @exception GException::invalid_argument
  *            Radius value is less than 0.
@@ -256,10 +258,9 @@ void GSkyRegionCircle::radius(const double&amp; radius)
     // Check if radius is valid
     if (radius < 0.0) {
         std::string msg =
-            "A radius of "+gammalib::str(radius)+" degrees has been"
-            " specified for a circular sky region.\n"
-            "The radius of a circular sky region can't be less than 0"
-            " degrees.";
+            "A negative radius of "+gammalib::str(radius)+" degrees has been "
+            "specified for a circular sky region. Please specify a "
+            "non-negative radius.";
         throw GException::invalid_argument(G_RADIUS, msg);
     }
@@ -296,8 +297,7 @@ void GSkyRegionCircle::read(const std::string&amp; line)
     if (region_def.find("circle") == std::string::npos) {
         std::string msg =
             "Unable to find the key word \"circle\" in provided string"
-            " \""+line+"\".\n"
-            "The \"circle\" key word is mandatory.";
+            " \""+line+"\". The \"circle\" key word is mandatory.";
         throw GException::invalid_value(G_READ, msg);
     }
@@ -314,8 +314,8 @@ void GSkyRegionCircle::read(const std::string&amp; line)
     std::vector<std::string> values = gammalib::split(circlestring,",");
     if (values.size() != 3) {
         std::string msg =
-            "Invalid number of "+gammalib::str(values.size())+" arguments"
-            " after the \"circle\" key word in provided string \""+line+"\".\n"
+            "Invalid number of "+gammalib::str(values.size())+" arguments "
+            "after the \"circle\" key word in provided string \""+line+"\". "
             "Exactly 3 arguments are expected.";
         throw GException::invalid_value(G_READ, msg);
     }
@@ -343,10 +343,9 @@ void GSkyRegionCircle::read(const std::string&amp; line)
     }
     else {
         std::string msg =
-            "Unsupported coordinate system \""+system+"\" in provided string"
-            " \""+line+"\".\n"
-            "Only the following coordinate systems are supported: "
-            "\"fk5\", \"icrs\" and \"galactic\".";
+            "Unsupported coordinate system \""+system+"\" in provided string "
+            "\""+line+"\". Only the following coordinate systems are "
+            "supported: \"fk5\", \"icrs\" and \"galactic\".";
         throw GException::invalid_value(G_READ, msg);
     }
@@ -479,11 +478,11 @@ std::string GSkyRegionCircle::print(const GChatter&amp; chatter) const
         // Append sky circle information
         result.append("\n"+gammalib::parformat("Right Ascension of centre"));
-        result.append(gammalib::str(m_centre.ra_deg())+" deg");
+        result.append(gammalib::str(ra())+" deg");
         result.append("\n"+gammalib::parformat("Declination of centre"));
-        result.append(gammalib::str(m_centre.dec_deg())+" deg");
+        result.append(gammalib::str(dec())+" deg");
         result.append("\n"+gammalib::parformat("Radius"));
-        result.append(gammalib::str(m_radius)+" deg");
+        result.append(gammalib::str(radius())+" deg");
     } // endif: chatter was not silent
@@ -539,9 +538,6 @@ bool GSkyRegionCircle::contains(const GSkyRegion&amp; reg) const
         // Calculate angular distance between the centres
         double ang_dist = m_centre.dist_deg(regcirc->centre());
-        // Calculate angular distance between the centres
-        double ang_dist = m_centre.dist_deg(regcirc->centre());
-
         // Check if the region is contained in this
         if ((ang_dist + regcirc->radius()) <= m_radius) {
             fully_inside = true;
@@ -557,16 +553,16 @@ bool GSkyRegionCircle::contains(const GSkyRegion&amp; reg) const
               dynamic_cast<const GSkyRegionRect*>(&reg);
         // Check containment of all edges
-        fully_inside = contains(regrect->get_corner(0)) &
-                       contains(regrect->get_corner(1)) &
-                       contains(regrect->get_corner(2)) &
+        fully_inside = contains(regrect->get_corner(0)) &&
+                       contains(regrect->get_corner(1)) &&
+                       contains(regrect->get_corner(2)) &&
                        contains(regrect->get_corner(3));
     }
     // ... otherwise throw an exception
     else {
         throw GException::feature_not_implemented(G_CONTAINS,
-              "Cannot compare to region type \""+reg.type()+"\" yet");
+              "Cannot compare to region type \""+reg.type()+"\" yet.");
     }
     // Return value
@@ -581,13 +577,16 @@ bool GSkyRegionCircle::contains(const GSkyRegion&amp; reg) const
  *
  * @exception GException::feature_not_implemented
  *            Regions differ in type.
+ *
+ * @todo Implement checks for rectangles and maps
  ***************************************************************************/
 bool GSkyRegionCircle::overlaps(const GSkyRegion& reg) const
 {
     // Initialise return value
     bool overlap = false;
-    // If other region is circle use a simple way to calculate
+    // If region is circle then check overlap by comparing the centre
+    // distances
     if (reg.type() == "Circle") {
         // Create circular region from reg
@@ -602,12 +601,36 @@ bool GSkyRegionCircle::overlaps(const GSkyRegion&amp; reg) const
             overlap = true;
         }
-    }
+    } // endif: region was circle
+
+    // ... otherwise if region is rectangle then check overlap
+    else if (reg.type() == "Rect") {
+
+        // Create rectangular region from reg
+        const GSkyRegionRect* regrect =
+              dynamic_cast<const GSkyRegionRect*>(&reg);
+
+        // Check overlap with circle
+        overlap = regrect->overlaps(*this);
+
+    } // endif: region was rectangle
+
+    // ... otherwise if region is map then check overlap with map
+    else if (reg.type() == "Map") {
+
+        // Create map from reg
+        const GSkyRegionMap* regmap =
+              dynamic_cast<const GSkyRegionMap*>(&reg);
+
+        // Check overlap with circle
+        overlap = regmap->overlaps(*this);
+
+    } // endif: region was map
     // ... otherwise throw an exception
     else {
-        throw GException::feature_not_implemented(G_OVERLAPS,
-              "Cannot compare two different region types yet");
+        throw GException::feature_not_implemented(G_CONTAINS,
+              "Cannot compare to region type \""+reg.type()+"\" yet.");
     }
     // Return value
@@ -28,7 +28,6 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
-//#include "GMath.hpp"
 #include "GTools.hpp"
 #include "GSkyRegionCircle.hpp"
 #include "GSkyRegionRect.hpp"
@@ -598,8 +597,7 @@ bool GSkyRegionRect::contains(const GSkyRegion&amp; reg) const
     // ... otherwise throw an exception
     else {
         throw GException::feature_not_implemented(G_CONTAINS,
-              "Cannot compare rectangular region with other than rectangle or"
-              " circle yet.");
+              "Cannot compare to region type \""+reg.type()+"\" yet.");
     }
     // Return value
@@ -615,14 +613,15 @@ bool GSkyRegionRect::contains(const GSkyRegion&amp; reg) const
  * @exception GException::feature_not_implemented
  *            Regions differ in type.
  *
- * @todo: - Improve implementation for rectangle-rectangle
+ * @todo Improve implementation for rectangle-rectangle
+ * @todo Implement check for map
  ***************************************************************************/
 bool GSkyRegionRect::overlaps(const GSkyRegion& reg) const
 {
     // Initialise return value
-    bool is_overlapping = false;
+    bool overlap = false;
-    // If other region is circle use a simple way to calculate
+    // If region is circle then check overlap between circle and rectangle
     if (reg.type() == "Circle") {
         // Create circular region from reg
@@ -641,33 +640,46 @@ bool GSkyRegionRect::overlaps(const GSkyRegion&amp; reg) const
             // Check extension of circle along right ascension axis
             if (std::abs(local_centre.ra_deg()) <= (0.5 * m_width + circ_rad)) {
-                // Circle is fully contained in this rectangle
-                is_overlapping = true;
+                // Circle is overlapping with rectangle
+                overlap = true;
             }
         }
-    } // Region was of type "Circle"
+    } // endif: region was of type "Circle"
+    // ... otherwise if region is a rectangle then check overlap between
+    // two rectangles
     else if (reg.type() == "Rect") {
         // Create rectangular region from reg
         const GSkyRegionRect* regrect =
               dynamic_cast<const GSkyRegionRect*>(&reg);
-        // Dirty kludge: compare vs map
+        // Dirty kludge: compare with map
         GSkyRegionMap regmap = GSkyRegionMap(regrect);
-        is_overlapping = regmap.overlaps(*this);
+        overlap = regmap.overlaps(*this);
-    } // Region was of type "Rect"
+    } // endif: region was of type "Rect"
+
+    // ... otherwise if region is map then check overlap with map
+    else if (reg.type() == "Map") {
+
+        // Create map from reg
+        const GSkyRegionMap* regmap =
+              dynamic_cast<const GSkyRegionMap*>(&reg);
+
+        // Check overlap with circle
+        overlap = regmap->overlaps(*this);
+
+    } // endif: region was map
     // ... otherwise throw an exception
     else {
-        throw GException::feature_not_implemented(G_OVERLAPS,
-              "Cannot compare rectangular region with other than rectangle or"
-              " circle yet.");
+        throw GException::feature_not_implemented(G_CONTAINS,
+              "Cannot compare to region type \""+reg.type()+"\" yet.");
     }
     // Return value
-    return is_overlapping;
+    return overlap;
 }
 /***************************************************************************
  *                 GSkyRegions.cpp - Sky region container class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2018 by Pierrick Martin                             *
+ *  copyright (C) 2013-2020 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -536,7 +536,7 @@ void GSkyRegions::load(const GFilename&amp; filename)
                 coordsys = "icrs";
             }
-            // If region is a circle
+            // If region is a circle then read circle
             if (std::string::npos != fileline.find("circle")) {
                 // Create instance of GSkyRegion object
@@ -558,32 +558,14 @@ void GSkyRegions::load(const GFilename&amp; filename)
 					region.read(newfileline);
 					append(region);
 				}
-			}
-            // If region is a rectangle
-            if (std::string::npos != fileline.find("box")) {
+			} // endif: region was circle
-                // Create instance of GSkyRegion object
-                GSkyRegionRect region;
-
-                // If coordinate system and region defined on the same line
-                if ((std::string::npos != fileline.find("fk5")) ||
-                    (std::string::npos != fileline.find("icrs")) ||
-                    (std::string::npos != fileline.find("galactic"))) {
-                    region.read(fileline);
-                    append(region);
-                }
-
-                // else, prepend the coordinate system
-                else {
-                coordsys = "icrs";
-            }
-
-            // If region is a circle
-            if (std::string::npos != fileline.find("circle")) {
+            // ... otherwise if region is a box then read rectangle
+            else if (std::string::npos != fileline.find("box")) {
                 // Create instance of GSkyRegion object
-                GSkyRegionCircle region;
+                GSkyRegionRect region;
                 // If coordinate system and region defined on the same line
                 if ((std::string::npos != fileline.find("fk5")) ||
@@ -601,8 +583,10 @@ void GSkyRegions::load(const GFilename&amp; filename)
                     region.read(newfileline);
                     append(region);
                 }
-            }
-        }
+
+            } // endif: region was box
+
+        } // endwhile: looped over file
         // Close file
         ds9file.close();
@@ -1913,7 +1913,7 @@ void TestGSky::test_GSkyRegionRect(void)
     test_try("Test posang assignment after");
     try {
         GSkyRegionRect rect(refdir_radeczerozero, ref_width, ref_height, ref_pa_deg);
-        rect.posang_deg(-1.0);
+        rect.posang(-1.0);
         test_try_success();
     }
     catch (GException::invalid_argument &e) {
@@ -1932,7 +1932,7 @@ void TestGSky::test_GSkyRegionRect(void)
     double check_height = rect.width();
     test_value(ref_height, check_height, 1.0e-10, "Test height assignment");
-    double check_pa_deg = rect.posang_deg();
+    double check_pa_deg = rect.posang();
     test_value(ref_pa_deg, check_pa_deg, 1.0e-10, "Test posang assignment");
     // Check solid angle assignment
@@ -1941,7 +1941,7 @@ void TestGSky::test_GSkyRegionRect(void)
     test_value(solidangle_check, solidangle, 1.0e-10, "Test solid angle");
     // Initialise objects
-    GSkyDir refdir_radeczerozero = GSkyDir();
+    refdir_radeczerozero = GSkyDir();
     GSkyDir test_dir_0 = GSkyDir();
     GSkyDir test_dir_1 = GSkyDir();
     GSkyDir test_dir_2 = GSkyDir();
@@ -2050,7 +2050,6 @@ void TestGSky::test_GSkyRegionRect(void)
     test_assert(refregion.overlaps(refregion_raoffset),"test3 for overlap");
     test_assert(!refregion.overlaps(refregion_decpole),"test4 for overlap");
-
     // Exit test
     return;
 }