GammaLib

Currently, GSkyMap does not consistently handle the upper edge boundary. Here’s a quick example script demonstrating the issue. It creates a 10x10 GSkyMap and tries to access the value at the top-center border of the map:

#include <GammaLib.hpp>
#include <iostream>

int main()
{
    // Create a test skymap
    GSkyMap skymap("CAR","CEL",0, 0, -0.5, 0.5, 10, 10);
    skymap = 3.14159;

    // Create a test pixel on the top-center border of the map
    GSkyPixel pix(5.5, 9.5);

    std::cout << "Skymap Info:\n" << skymap.print() << std::endl;
    std::cout << "Pixel Info:\n" << pix.print() << std::endl;

    std::cout << "\nPixel is contained: ";
    if (skymap.contains(pix)) {
        std::cout << "yes\n\n";
    } else {
        std::cout << "no\n\n";
    }

    // Access the map value by pixel
    std::cout << "Trying to access map value by pixel...\n";
    try {
        double val(skymap(pix, 0));
        std::printf("   succeeded (val = %f)\n", val);
    } catch (...) {
        std::cout << "   failed\n";
    }

    // Access the pixel by index
    int indx = skymap.pix2inx(pix);
    std::printf("Trying to access map value by index (index = %d)...\n", indx);
    try {
        double val(skymap(indx, 0));
        std::printf("   succeeded (val = %f)\n", val);
    } catch (...) {
        std::cout << "   failed\n";
    }

    return 0;
}

which yields the following output:

Pixel is contained: yes

Trying to access map value by pixel...
   succeeded (val = 0.000000)
Trying to access map value by index (index = 106)...
   failed

It appears that the border counts as being contained in the map, however accessing directly via the pixel does not return a value from in the map and accessing by index actually throws an exception. Here’s a breakdown of my understanding of what’s going on:

• skymap.contains(GSkyPixel&) tests that x-pixel value is within [-0.5, m_num_x-0.5] with the upper bound inclusive (same for the y-pixel value), so pixel values that fall on any of the borders are considered as being inside the map.
• skymap(int&) simply returns the pixel value at indx, or throws an exception if indx is out of bounds.
• skymap(GSkyPixel&) first checks that the pixel is contained (which in the above example it is) then converts the pixel into an index and returns the pixel value at that index. However, in the above example the index is out-of-bounds, the pixel queried is beyond the end of the array and the returned value will be unpredictable.

The simplest change would be to make the upper bound on the skymap in x & y non-inclusive, which would prevent 'contains’ from returning true resulting in both issues involving the contains() and skymap(GSkyPixel&) methods to be resolved. Alternatively, a check could be done for the special case when the pixel falls on the upper boundary and be handled accordingly by pix2inx.