GammaLib

Back in 2009, I asked the question already to Mark Calabretta (Mister WCS). He was not aware of some software doing this job:

managed to get a general conversion. Do you know about a collection
of formula that can be used to derive the solid angles of individual
image pixels or are there some subroutines in the wcslib to do this
task (or are there some subroutines available for idl)?

Dear Jürgen,

It’s difficult in general to evaluate pixel areas analytically - a
square pixel in the plane of projection may become severely distorted
when projected back onto the sphere. Cylindrical projections and equal
area projections are really the only easy cases, numerical methods must
be used for the others. However, I’m not aware of any software that
can do this.

Regards,
Mark Calabretta

I put the priority of this issue to High as we should really not forget about this. Maybe we need a maths geek who works out the formula, or who can propose a solid numerical method.

https://github.com/astrofrog/python-reprojection/blob/master/reproject/overlapArea.c
is now under a BSD license, i.e. I think we can include a copy in gammalib.

Tom and I are currently improving it and once it’s in good shape (no global variables, good test coverage) in a week or two we can put a copy in gammalib:
https://github.com/astrofrog/python-reprojection/issues

The core functionality implemented there is flux-conserving reprojection, not only computation of pixel solid angles.
I think we want flux-conserving reprojection methods in gammalib as well, e.g. it could come in handy to stack runs into survey maps or to reproject background models in the FOV system onto sky maps.

Jürgen, do you agree?
If yes, can you please propose a place and outline what API you think would be nice (class method? function?)

Ellis, could you please implement the method to compute the solid angle of a given pixel described here?
http://mail.scipy.org/pipermail/astropy/2013-December/002940.html

I replaced the existing solid angle script the following:

{
    GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
    GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
    GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
    GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));

    GVector vec1 = dir1.celvector();
    GVector vec2 = dir2.celvector();
    GVector vec3 = dir3.celvector();
    GVector vec4 = dir4.celvector();

    double angle1 = std::acos(cross(vec2, (cross(vec1, vec2))) * cross(vec2, (cross(vec3, vec2))));
    double angle2 = std::acos(cross(vec3, (cross(vec2, vec3))) * cross(vec3, (cross(vec4, vec3))));
    double angle3 = std::acos(cross(vec4, (cross(vec3, vec4))) * cross(vec4, (cross(vec1, vec4))));
    double angle4 = std::acos(cross(vec1, (cross(vec4, vec1))) * cross(vec1, (cross(vec2, vec1))));

    double solidangle = (angle1 + angle2 + angle3 + angle4) - (2 * gammalib::pi);

    return solidangle;
}

Which gives clearly incorrect results: I would expect the map to have pixel values of approx 1 but they are around 1e-4 (and quite variable) - see attached ds9 screenshot - the test script is here: https://gist.github.com/ellisowen/8473696

I shall start de-bugging, but wondered if you could see anything obviously wrong with the method/script?

Also, would it be best to include this as an alternative to the original solidangle calculation, or as a replacement?

(Branch is here: https://github.com/ellisowen/gammalib/compare/issue_1017)

I got the following result with this code:

    //
    double a;
    double b;
    double c;

    // Angle 1
    a = dir1.dist(dir3);
    b = dir4.dist(dir3);
    c = dir4.dist(dir1);
    double angle1 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));

    // Angle 2
    a = dir4.dist(dir2);
    b = dir2.dist(dir3);
    c = dir3.dist(dir4);
    double angle2 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));

    // Angle 3
    a = dir1.dist(dir3);
    b = dir1.dist(dir2);
    c = dir2.dist(dir3);
    double angle3 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));

    // Angle 4
    a = dir4.dist(dir2);
    b = dir4.dist(dir1);
    c = dir1.dist(dir2);
    double angle4 = std::acos((std::cos(a) - std::cos(b)*std::cos(c))/(std::sin(b)*std::sin(c)));

Looks reasonable. The code uses spherical triangles to compute the inner angles of the Polygon. I guess you code attempts to do the same thing, but maybe some stuff got mixed up.

My code is not yet optimized as I repeat distance computations. But I just wanted to see whether this works.

Just think: can still be optimized a little more as there are repeated cos and sin calls.

Okay, I think I also have figured out the problem with my method. I avoid calling cos and sin so much, so maybe it would need less optimizing?

I’ll push it to https://github.com/ellisowen/gammalib/compare/issue_1017 when I’ve done a couple of checks!

(Changed % done by accident so just returned it to 50%)

Pushed to https://github.com/ellisowen/gammalib/compare/issue_1017 - please review.

Attached check projection ds9 screenshot (it’s not very clear, but values seem reasonable).

Sorry to say that there is still some problem.

I tried your code with a pixel size of (0.5 x 0.5) instead of (1.0 x 1.0), and when you integrate over the map you get a total solid angle that is 4 times too large. With the trigonometric formula the correct values are obtained.

For the (1.0 x 1.0) pixel size I still get an error of a few percent when being close to the pole, so the formula seem only to be approximate.

Can you point to where you took the cross-product / dot-product formula from? Could it be that they only work for cartesian coordinates?

I modified the code so that both methods can be toggled using a compiler flag. This allows you to investigate this further. Pushed into issue_1017_jk.

I optimized the code a little more for the trigonometric version, and also take care of the possibility that a pixel touches the pole, reducing the polygon to a triangle. I tested the code with the script attached.

Below also a number of all sky maps produced by the script for various projections. They seem to be all okay (at least visually). The script now also computes the total solid angle of the map, and at least non of the projections led to a NaN. AIT and CAR give 4 pi as expected. Eventually, some more testing would be nice.

Owen Ellis wrote:

The factor of 4 issue you mentioned I think is because I just noticed I hadn’t adapted the term ((gammalib::pi * (1.0/180.0)) in the solid angle formula to account for different pixel sizes.

However, I’d tend to agree with Christoph - the formula in my code clearly isn’t quite right (even with the above problem corrected), and it’s not yet clear to me how to fix the other issues. So I think it would perhaps better to just use your working code?

For the moment I left all options in the code, can be selected through #define. In the long run we may remove the obsolete code, but we better keep for the moment the possibility to play with the options.

I now merged all the stuff in devel.