Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix map pixel solid angle computation #2276

Merged
merged 4 commits into from Jul 9, 2019
Merged

Conversation

@adonath
Copy link
Member

@adonath adonath commented Jul 5, 2019

The computation of the solid angle in WCSGeom.solid_angle() was based on the approximation of rectangular pixels. This lead to problems for WCS projections, where strong distortions of the pixels occurs, such as an all sky AIT projection or local projections at high latitudes (see #2150 (comment), TODO comment in the code, and privately reported by @lmohrmann).

This PR changes the current implementation to a method that splits the pixel into two spherical triangles and uses the spherical excess theorem to compute the area of the triangles. A short comparison between the old and new method is presented in https://github.com/adonath/notebooks-public/blob/master/solid_angle.ipynb.

From a quick test it seems the new method is ~4 times slower. For now I would propose to change to the new slower, but more correct method and later think about optimisations, or alternative faster methods.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 5, 2019

@adonath - Thanks!

I wanted to fix this, but didn't get to it.

My idea was to use this:
https://github.com/astropy/regions/blob/ffdf4d4bf3eadb9498a2cc4d62b359933db63df5/regions/_utils/wcs_helpers.py#L13

I think it's very similar to what you do, except that it uses the local planar approximation within a given pixel. I think that will give pretty much the same precision to what you do, for some projections and positions it will be better, and for some a tiny bit worse. Neither approximation is 100% correct for very large pixels, you assume pixel edges are great circles, which isn't the case. I would guess the planar approximation is usually a bit better for typical pixel sizes that we have, which are often 0.01 deg or 0.1 deg or 1 deg, but not larger.

I think performance doesn't matter much, but still I would suggest that you put the planar approximation or just assume your two triangles are the same size, which removes some code and should speed up by factor ~ 2.

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jul 7, 2019

Thanks @adonath. This looks good.
You might also test in your notebook a TAN projection at the northern pole.

An option to have a faster method on average could be to compute the solid angle size analytically for the projection when it is feasible, CAR obviously and probably TAN as well.

@cdeil cdeil added the bug label Jul 8, 2019
@cdeil cdeil added this to the 0.13 milestone Jul 8, 2019
@cdeil cdeil added this to To do in gammapy.maps via automation Jul 8, 2019
@cdeil cdeil changed the title Improve precision of solid angle computation Fix map pixel solid angle computation Jul 8, 2019
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

@adonath - I've labeled this "bug" and changed the title to "Fix ...".

The previous algorithm really was incorrect, assuming the angle was always 90 deg, where really it's often not, and combined with the fact that we haven't adjusted the background norm so far, the background estimates so far were incorrect whenever there was a non-trivial WCS.

@adonath
Copy link
Member Author

@adonath adonath commented Jul 8, 2019

Thanks for the comments @cdeil and @registerrier! I've updated the notebook (https://github.com/adonath/notebooks-public/blob/master/solid_angle.ipynb) with a TAN projection at the pole and also added the simplified method @cdeil proposed (assuming a planar triangle and symmetry of the pixel). It seems the simple method is precise enough, 2 times faster and less code. So I'll change to the simplified method, adapt the tests and merge.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

@adonath - Great, thanks!

I'm not sure what you mean by "planar triangle".

I would think that "planar" means using the differential WCS scale like here, not calling a spherical distance function any more.

You still call and use the spherical distance of the pixel edges, no?
That's fine with me, just don't call it "planar approximation" in the comments or something like that.

(I'm not really 100% sure how the spherical excess and the differential local WCS relate. I think they are different algorithms, but I'm not really sure if they don't work out to be basically the same.)

@adonath
Copy link
Member Author

@adonath adonath commented Jul 8, 2019

@cdeil The question is whether you compute the area of the triangle assuming great circles or straight lines as sides. In one case you just use area = a * b * sin(C), in the other you need some spherical excess theorem. That's what I meant by "planar triangle".

However the symmetric pixel assumption, where the area is computed from one triangle only, leads to asymmetries in the lat direction. While the effect is minimal and probably negligible
in practice, it leads to somehow un-expected results for the unit tests we have. We compute solid angles for 3x3 pixel images at lon=0 and those come out not symmetric in lat direction, which I don't like, because the result is un-intuitive and must be explained and I'm sure it will lead to user complains in future (in fact we already had the issue in #2035). I think I'm leaning towards the 2 triangle solution again, which assures the symmetry in lat direction...

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

Well, as you know, I'm -1 on the two-triangle solution. Like any other method it's an approximation, and in practice these 0.01 % differences in pixel solid angle are really irrelevant, we should go for the least amount and simplest code possible.

If really ever higher-precision is needed, then they scalable and simple way to do it is to downsample the WCS e.g. split each pixel into 4 or 100, compute their areas and sum. That's what we do e.g. in regions, or also in model evaluation, so doing it again is better than something special and new, like two triangles.

@adonath
Copy link
Member Author

@adonath adonath commented Jul 8, 2019

@cdeil In terms of absolute precision I completely agree, but still my point is, that in unit-tests as well as for users the following result is completely un-intuitive and must be explained:

from gammapy.maps import WcsGeom
from astropy import units as u
geom = WcsGeom.create(skydir=(0, 0), coordsys="GAL", npix=(3, 3), binsz=1 * u.deg)
print(geom.solid_angle().value)
[[0.00030461 0.00030461 0.00030461]
 [0.00030461 0.00030461 0.00030461]
 [0.00030451 0.00030451 0.00030451]]

While the 2 triangle solution gives a symmetric array, which is what you naively expect:

[[0.00030456 0.00030456 0.00030456]
 [0.00030461 0.00030461 0.00030461]
 [0.00030456 0.00030456 0.00030456]]

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jul 8, 2019

I also prefer the 2 triangles solution since it is required to solve an issue that triggered a former update of the code PR #2035.

The actual difference in complexity is really tiny. Both solutions require complex code anyway.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

I don't think the unit tests should drive our algorithms, and having 9 pixels in the reference values that basically all have the same solid angle isn't useful either, and users don't look at our unit tests anyways. If you say in the docstring of solid_angle that it's an approximation that's fine. And if you really care about providing better accuracy, offer a subpixel argument like https://astropy-regions.readthedocs.io/en/latest/masks.html or how I think we already have it in other places in Gammapy. (I don't think that's needed, but that would be a good way to introduce better accuracy method if needed)

But @adonath - I see you're fond of the 2-triangle method. OK to merge this in if you think it's better.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

An option to have a faster method on average could be to compute the solid angle size analytically for the projection when it is feasible, CAR obviously and probably TAN as well.

Just in case it isn't clear: I'm strongly -1 to add such code. Solid angle is computed once in setup, not in the likelihood loop. Performance doesn't matter. I strongly prefer to have as little and as simple code in Gammapy as possible, only add special code where really needed.

@registerrier
Copy link
Contributor

@registerrier registerrier commented Jul 8, 2019

Well, if performance does not matter, precision does. I find very bothering to have a solid_angle method that return something without any indication on the actual precision achieved and to expect a user to know how to oversample its WCS to reach a decent precision.

Code simplicity cannot always be the main driver.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

@registerrier - Then putting an oversample=2 that split each pixel in 4 parts is better. It gives more precise results than two triangles, and gives the user a clear and simple option to get more precision if needed.

@cdeil
Copy link
Member

@cdeil cdeil commented Jul 8, 2019

We already have such an oversample factor for map-based computations to get higher precision here: https://docs.gammapy.org/0.12/api/gammapy.cube.PSFKernel.html#gammapy.cube.PSFKernel.from_gauss

@adonath adonath force-pushed the improve_solid_angle branch from f241d59 to 2a72df0 Jul 8, 2019
@adonath
Copy link
Member Author

@adonath adonath commented Jul 9, 2019

@cdeil I kept the 2 triangle version, because of the symmetry argument. I think for these kind of algorithms (or similar one like bounding box computations etc.) it easily happens, that you get shifts by half a pixel or one pixel and I want this to be covered by simple unit-tests, which are easy to understand.

I think it's clear that we keep the current implementation until we find that we need a more precise method (which might never happen...). If the need for more precise solid angle comes up, we'll re-discuss our options and decide between oversampling vs. analytical solutions, based on the use-case we found.

@adonath adonath merged commit b12ee74 into gammapy:master Jul 9, 2019
8 of 9 checks passed
gammapy.maps automation moved this from To do to Done Jul 9, 2019
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 9, 2019

@adonath - Thanks for fixing this bug!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
gammapy.maps
  
Done
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants