Sergio Hernández Cadena / gammalib

Commit 48b8040c7d7aa913677cb8628775b24706177c48

Authored by Jürgen Knödlseder
1 parent e4d91b6e

Avoid pixel validity check in GSkyMap::flux() methods (#1785)

Inline Side-by-side
Showing 1 changed file with 30 additions and 10 deletions
src/sky/GSkyMap.cpp
  View file @48b8040
... ... @@ -71,6 +71,8 @@
71 71 #define G_DIR2PIX "GSkyMap::dir2pix(GSkyDir&)"
72 72 #define G_FLUX1 "GSkyMap::flux(int&)"
73 73 #define G_FLUX2 "GSkyMap::flux(GSkyPixel&)"
  74 +#define G_FLUX3 "GSkyMap::flux(GSkyRegion&, int&)"
  75 +#define G_FLUX4 "GSkyMap::flux(GSkyRegions&, int&)"
74 76 #define G_SOLIDANGLE1 "GSkyMap::solidangle(int&)"
75 77 #define G_SOLIDANGLE2 "GSkyMap::solidangle(GSkyPixel&)"
76 78 #define G_OVERLAPS "GSkyMap::overlaps(GSkyRegion&)"
... ... @@ -1730,6 +1732,9 @@ double GSkyMap::flux(const GSkyPixel& pixel, const int& map) const
1730 1732 * @param[in] map Map index [0,...,nmaps()-1].
1731 1733 * @return Flux in sky region.
1732 1734 *
  1735 + * @exception GException::out_of_range
  1736 + * Map index out of valid range
  1737 + *
1733 1738 * The method returns the total flux in all sky map pixels that fall within
1734 1739 * the specified region. Containment is tested using the
1735 1740 * GSkyRegion::contains() method. Sky map values are multiplied by the solid
... ... @@ -1737,14 +1742,20 @@ double GSkyMap::flux(const GSkyPixel& pixel, const int& map) const
1737 1742 ***************************************************************************/
1738 1743 double GSkyMap::flux(const GSkyRegion& region, const int& map) const
1739 1744 {
  1745 + // Throw an error if map index is not in valid range
  1746 + if (map < 0 || map >= m_num_maps) {
  1747 + throw GException::out_of_range(G_FLUX3, "Map index",
  1748 + map, m_num_maps);
  1749 + }
  1750 +
1740 1751 // Initialise flux
1741 1752 double flux = 0.0;
1742 1753  
1743 1754 // Loop over all sky pixels and sum flux of all pixels that fall within
1744 1755 // the sky region
1745   - for (int i = 0; i < m_num_pixels; ++i) {
  1756 + for (int i = 0, ipix = m_num_pixels * map; i < m_num_pixels; ++i, ++ipix) {
1746 1757 if (region.contains(inx2dir(i))) {
1747   - flux += (*this)(i,map) * solidangle(i);
  1758 + flux += m_pixels(ipix) * solidangle(i);
1748 1759 }
1749 1760 }
1750 1761  
... ... @@ -1760,6 +1771,9 @@ double GSkyMap::flux(const GSkyRegion&amp; region, const int&amp; map) const
1760 1771 * @param[in] map Map index [0,...,nmaps()-1].
1761 1772 * @return Flux in sky region.
1762 1773 *
  1774 + * @exception GException::out_of_range
  1775 + * Map index out of valid range
  1776 + *
1763 1777 * The method returns the total flux in all sky map pixels that fall within
1764 1778 * the specified regions. Containment is tested using the
1765 1779 * GSkyRegions::contains() method. Sky map values are multiplied by the solid
... ... @@ -1767,14 +1781,20 @@ double GSkyMap::flux(const GSkyRegion&amp; region, const int&amp; map) const
1767 1781 ***************************************************************************/
1768 1782 double GSkyMap::flux(const GSkyRegions& regions, const int& map) const
1769 1783 {
  1784 + // Throw an error if map index is not in valid range
  1785 + if (map < 0 || map >= m_num_maps) {
  1786 + throw GException::out_of_range(G_FLUX4, "Map index",
  1787 + map, m_num_maps);
  1788 + }
  1789 +
1770 1790 // Initialise flux
1771 1791 double flux = 0.0;
1772 1792  
1773 1793 // Loop over all sky pixels and sum flux of all pixels that fall within
1774 1794 // the sky region
1775   - for (int i = 0; i < m_num_pixels; ++i) {
  1795 + for (int i = 0, ipix = m_num_pixels * map; i < m_num_pixels; ++i, ++ipix) {
1776 1796 if (regions.contains(inx2dir(i))) {
1777   - flux += (*this)(i,map) * solidangle(i);
  1797 + flux += m_pixels(ipix) * solidangle(i);
1778 1798 }
1779 1799 }
1780 1800  
... ... @@ -1797,15 +1817,15 @@ GNdarray GSkyMap::flux(void) const
1797 1817 GNdarray flux_array = GNdarray(m_num_maps);
1798 1818  
1799 1819 // Loop over maps and store sums in output array
1800   - for (int i = 0; i < m_num_maps; ++i) {
  1820 + for (int i = 0, ipix = 0; i < m_num_maps; ++i) {
1801 1821  
1802 1822 // Initialise flux for this map
1803 1823 double flux = 0.0;
1804 1824  
1805 1825 // Compute map flux
1806   - for (int k = 0; k < m_num_pixels; ++k) {
  1826 + for (int k = 0; k < m_num_pixels; ++k, ++ipix) {
1807 1827 double solidangle = this->solidangle(k);
1808   - flux += (*this)(k,i) * solidangle;
  1828 + flux += m_pixels(ipix) * solidangle;
1809 1829 }
1810 1830  
1811 1831 // Store map flux
... ... @@ -1924,13 +1944,13 @@ double GSkyMap::solidangle(const GSkyRegion&amp; region) const
1924 1944  
1925 1945  
1926 1946 /***********************************************************************//**
1927   - * @brief Returns flux within sky regions
  1947 + * @brief Returns solid angle within sky regions
1928 1948 *
1929   - * @param[in] region Sky region.
  1949 + * @param[in] regions Sky regions.
1930 1950 * @return Solid angle within sky region.
1931 1951 *
1932 1952 * The method returns the total solid angle of all sky map pixels that fall
1933   - * within the specified region. Containment is tested using the
  1953 + * within the specified regions. Containment is tested using the
1934 1954 * GSkyRegions::contains() method.
1935 1955 ***************************************************************************/
1936 1956 double GSkyMap::solidangle(const GSkyRegions& regions) const
... ...