spacetelescope / spherical_geometry

A Python package for handling spherical polygons that represent arbitrary regions of the sky
http://spherical-geometry.readthedocs.io/
61 stars 31 forks source link

SphericalPolygon.intersection sometimes returns the complement of the intersection #278

Open Hyfro-Cyrup opened 3 months ago

Hyfro-Cyrup commented 3 months ago

I have noticed that in some cases, calculating the intersection of two simple polygons can result in their complement. I have not delved deeply into the graph module to find out why this happens or produce other counterexamples, but here is a very simple one:

from spherical_geometry.polygon import SphericalPolygon

p1 = SphericalPolygon.from_cone(90, 0, 100)
p2 = SphericalPolygon.from_cone(270, 0, 100)

print(p1.contains_lonlat(0, 0)) # true
print(p2.contains_lonlat(0, 0)) # true
print(p1.intersection(p2).contains_lonlat(0, 0)) # erroneously false

I defined two caps on either side of the sphere whose intersection should be a thin strip around the prime meridian, or at very least empty if I've defined my polygons wrong.

Not sure why this is happening. I'll try and read through the graph module when I have time unless someone has an idea first. This may be related to #125; they mentioned something similar happening, but the consensus there seemed to be about self-intersection.

version 1.3.1

Hyfro-Cyrup commented 3 months ago

I believe I have narrowed the issue down to this part of the intersection code:

    poly = self._trace()
    # If multiple polygons, the inside point can only be in one
    if len(poly._polygons) == 1 and not self._contains_inside_point(poly):
        poly = poly.invert_polygon()

The polygon created from the trace doesn't have the right interior, so (I believe) it tries to check for that and correct it using invert_polygon. The issue gets created when invert_polygon does not invert the polygon because of an issue with _get_new_outside: it assumes that points antipodal to the polygon's defining points lie outside the polygon. In my example with a circle of radius >=90, this is clearly false.

_get_new_inside appears to have a similar issue, in that it doesn't take into account which side of the polygon is actually inside, though I don't think that's affecting the intersection bug. Here's a snippet that shows both issues

from spherical_geometry.polygon import SphericalPolygon

p1 = SphericalPolygon.from_cone(90, 0, 100)

print(p1.contains_point(p1.polygons[0]._find_new_inside())) # erroneously false
print(p1.contains_point(p1.polygons[0]._find_new_outside())) # erroneously true