From 961a8778a36b764b2a2384f0b66b19b3ee35cb00 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 8 Feb 2015 20:52:32 +0100 Subject: [PATCH] Add scaling operators to GSkymap class. --- ChangeLog | 3 ++- NEWS | 5 ++++- include/GSkymap.hpp | 2 ++ pyext/GSkymap.i | 11 +++++++++++ src/sky/GSkymap.cpp | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---- test/test_GSky.cpp | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++--- test/test_GSky.py | 46 ++++++++++++++++++++++++++++++---------------- 7 files changed, 162 insertions(+), 25 deletions(-) diff --git a/ChangeLog b/ChangeLog index 2897739..db1d53c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2015-02-05 Juergen Knoedlseder +2015-02-08 Juergen Knoedlseder * Version 1.0.0 released ======================== @@ -13,6 +13,7 @@ Add GBilinear class Add WCS projections (MOL) Add GXmlElement::has_attribute() method + Add GSkymap scaling operators (#1296) 2015-01-07 Juergen Knoedlseder diff --git a/NEWS b/NEWS index de4b0ed..74ce9af 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.0.0 Author(s) : Juergen Knoedlseder -5 February 2015 +8 February 2015 1. Introduction @@ -81,6 +81,9 @@ None The following WCS projections have been added: * GWcsMOL (Mollweide's projection) +Scaling operators operator*=(double) and operator/=(double) have been added +to GSkymap to allow scaling of all pixel values by a given value. + 13. Support module ------------------ diff --git a/include/GSkymap.hpp b/include/GSkymap.hpp index 6df77ff..b560ec3 100644 --- a/include/GSkymap.hpp +++ b/include/GSkymap.hpp @@ -107,7 +107,9 @@ public: GSkymap& operator+=(const GSkymap& map); GSkymap& operator-=(const GSkymap& map); GSkymap& operator*=(const GSkymap& map); + GSkymap& operator*=(const double& factor); GSkymap& operator/=(const GSkymap& map); + GSkymap& operator/=(const double& factor); double& operator()(const int& index, const int& map = 0); const double& operator()(const int& index, const int& map = 0) const; double& operator()(const GSkyPixel& pixel, const int& map = 0); diff --git a/pyext/GSkymap.i b/pyext/GSkymap.i index 5715214..1be3dd0 100644 --- a/pyext/GSkymap.i +++ b/pyext/GSkymap.i @@ -131,6 +131,7 @@ public: GSkymap& operator+=(const GSkymap& map); GSkymap& operator-=(const GSkymap& map); GSkymap& operator*=(const GSkymap& map); + GSkymap& operator*=(const double& factor); // Methods void clear(void); @@ -211,4 +212,14 @@ public: self->operator/=(map); return (*self); } + // Python 2.x operator/= + GSkymap __idiv__(const double& factor) { + self->operator/=(factor); + return (*self); + } + // Python 3.x operator/= + GSkymap __itruediv__(const double& factor) { + self->operator/=(factor); + return (*self); + } }; diff --git a/src/sky/GSkymap.cpp b/src/sky/GSkymap.cpp index 560515e..c7aa8c4 100644 --- a/src/sky/GSkymap.cpp +++ b/src/sky/GSkymap.cpp @@ -357,7 +357,7 @@ GSkymap& GSkymap::operator+=(const GSkymap& map) } // endfor: looped over all pixels - // Return present sky map + // Return this object return *this; } @@ -408,7 +408,7 @@ GSkymap& GSkymap::operator-=(const GSkymap& map) } // endfor: looped over all pixels - // Return present sky map + // Return this object return *this; } @@ -458,7 +458,31 @@ GSkymap& GSkymap::operator*=(const GSkymap& map) } // endfor: looped over all pixels - // Return present sky map + // Return this object + return *this; +} + + +/***********************************************************************//** + * @brief Sky map scaling operator + * + * @param[in] factor Scale factor. + * @return Sky map. + * + * Multiplies all pixels of the sky map by the given scale @p factor. + ***************************************************************************/ +GSkymap& GSkymap::operator*=(const double& factor) +{ + // Compute total number of pixels + int n = npix() * nmaps(); + + // Loop over all pixels + double* pixel = m_pixels; + for (int i = 0; i < n; ++i) { + *pixel++ *= factor; + } + + // Return this object return *this; } @@ -517,7 +541,41 @@ GSkymap& GSkymap::operator/=(const GSkymap& map) } // endfor: looped over all pixels - // Return present sky map + // Return this object + return *this; +} + + +/***********************************************************************//** + * @brief Sky map division operator + * + * @param[in] factor Scale factor. + * @return Sky map. + * + * @exception GException::invalid_argument + * Division by zero error. + * + * Divides all pixels of the sky map by the given @p factor. + ***************************************************************************/ +#define G_OP_UNARY_DIV2 "GSkymap::operator/=(double&)" +GSkymap& GSkymap::operator/=(const double& factor) +{ + // Check for division by zero + if (factor == 0.0) { + std::string msg = "Trying to divide sky map pixels by zero."; + throw GException::invalid_argument(G_OP_UNARY_DIV2, msg); + } + + // Compute total number of pixels + int n = npix() * nmaps(); + + // Loop over all pixels + double* pixel = m_pixels; + for (int i = 0; i < n; ++i) { + *pixel++ /= factor; + } + + // Return this object return *this; } diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 06d9bed..f7b13e5 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -851,7 +851,7 @@ void TestGSky::test_GSkymap(void) map_dst += map_src; // Check total in destination map. Note that the total in the destination - // map should be 100 times the total in the source maps as the operator + // map should be 100 times the total in the source map as the operator // is expected to perform strict bi-linear interpolation double total_dst = 0.0; for (int pix = 0; pix < map_dst.npix(); ++pix) { @@ -860,7 +860,7 @@ void TestGSky::test_GSkymap(void) } } total_dst /= 100.0; - test_value(total_dst, total_src, 1.0e-3, "Test operator+"); + test_value(total_dst, total_src, 1.0e-3, "Test operator+(GSkymap)"); // Subtract pixels from destination map map_dst -= map_src; @@ -873,7 +873,55 @@ void TestGSky::test_GSkymap(void) total_dst += map_dst(pix,k); } } - test_value(total_dst, 0.0, 1.0e-3, "Test operator+"); + test_value(total_dst, 0.0, 1.0e-3, "Test operator-(GSkymap)"); + + // Check map multiplication + GSkymap test_map = map_src; + test_map *= map_src; + double total_test = 0.0; + double total_ref = 0.0; + for (int pix = 0; pix < test_map.npix(); ++pix) { + for (int k = 0; k < test_map.nmaps(); ++k) { + total_test += test_map(pix,k); + total_ref += map_src(pix,k) * map_src(pix,k); + } + } + test_value(total_test, total_ref, 1.0e-3, "Test operator*(GSkymap)"); + + // Check map division + test_map = map_src; + test_map /= map_src; + total_test = 0.0; + total_ref = 0.0; + for (int pix = 0; pix < test_map.npix(); ++pix) { + for (int k = 0; k < test_map.nmaps(); ++k) { + total_test += test_map(pix,k); + total_ref += map_src(pix,k) / map_src(pix,k); + } + } + test_value(total_test, total_ref, 1.0e-3, "Test operator/(GSkymap)"); + + // Check map scaling + test_map = map_src; + test_map *= 3.3; + total_test = 0.0; + for (int pix = 0; pix < test_map.npix(); ++pix) { + for (int k = 0; k < test_map.nmaps(); ++k) { + total_test += test_map(pix,k); + } + } + test_value(total_test, total_src*3.3, 1.0e-3, "Test operator*(double)"); + + // Check map division + test_map = map_src; + test_map /= 3.3; + total_test = 0.0; + for (int pix = 0; pix < test_map.npix(); ++pix) { + for (int k = 0; k < test_map.nmaps(); ++k) { + total_test += test_map(pix,k); + } + } + test_value(total_test, total_src/3.3, 1.0e-3, "Test operator/(double)"); // Save maps map_src.save("test_map_src.fits", true); diff --git a/test/test_GSky.py b/test/test_GSky.py index b574a44..e45d491 100755 --- a/test/test_GSky.py +++ b/test/test_GSky.py @@ -332,31 +332,45 @@ class Test(GPythonTestSuite): # Addition operator map += map_b - self.test_value(map[0], 2.0); - self.test_value(map[1], 4.0); - self.test_value(map[2], 6.0); - self.test_value(map[3], 8.0); + self.test_value(map[0], 2.0) + self.test_value(map[1], 4.0) + self.test_value(map[2], 6.0) + self.test_value(map[3], 8.0) # Multiplication operator map *= map_b - self.test_value(map[0], 2.0); - self.test_value(map[1], 8.0); - self.test_value(map[2], 18.0); - self.test_value(map[3], 32.0); + self.test_value(map[0], 2.0) + self.test_value(map[1], 8.0) + self.test_value(map[2], 18.0) + self.test_value(map[3], 32.0) # Subtraction operator map -= map_b - self.test_value(map[0], 1.0); - self.test_value(map[1], 6.0); - self.test_value(map[2], 15.0); - self.test_value(map[3], 28.0); + self.test_value(map[0], 1.0) + self.test_value(map[1], 6.0) + self.test_value(map[2], 15.0) + self.test_value(map[3], 28.0) # Division operator map /= map_b - self.test_value(map[0], 1.0); - self.test_value(map[1], 3.0); - self.test_value(map[2], 5.0); - self.test_value(map[3], 7.0); + self.test_value(map[0], 1.0) + self.test_value(map[1], 3.0) + self.test_value(map[2], 5.0) + self.test_value(map[3], 7.0) + + # Scaling operator + map *= 2.0 + self.test_value(map[0], 2.0) + self.test_value(map[1], 6.0) + self.test_value(map[2], 10.0) + self.test_value(map[3], 14.0) + + # Division operator + map /= 2.0 + self.test_value(map[0], 1.0) + self.test_value(map[1], 3.0) + self.test_value(map[2], 5.0) + self.test_value(map[3], 7.0) # Return return -- libgit2 0.21.4