1

How can I determine if a polygon is self-touching(!) like the blue one, but not self-intersecting like the red one

enter image description here

Vertices for the two polygons shown in the images is here:

from shapely.geometry import Polygon
from shapely.affinity import translate
import pylab

pylab.figure()

points = [[0, 0], [0, 100], [50, 80], [0, 60], [50, 30]]
polygon = Polygon(points)
bounds = polygon.bounds

x,y = polygon.exterior.xy
pylab.plot(x,y, 'b')
pylab.plot(0, 60, 'ko')

points = [[0, 0], [0, 100], [50, 80], [-10, 60], [50, 30]]
polygon = Polygon(points)
polygon = translate(polygon, xoff=70)

x,y = polygon.exterior.xy
pylab.plot(x,y, 'r')

pylab.axis('equal')
pylab.grid(True)
pylab.show()
nickponline
  • 25,354
  • 32
  • 99
  • 167

1 Answers1

1

An important remark is that the intersection point on the self touching polygon intersects with polygon segments three times, whereas on the strictly self intersecting polygon, intersection points intersect only twice with the polygon segments.

According to this (hopefully, it generalizes well), we can build a method that get intersection points and then checks how many times this point intersects with the polygon segments:

  • A strictly self intersecting polygon will only have at least one intersection point that intersects twice with the polygon
  • A strictly self touching polygon will have all its intersection points intersecting three times.

To speed up the computation, we can use SRTree (see this post). I also use the polygon to get consecutive segments and points with this post.

from shapely.geometry import Polygon, Point, LineString
from shapely.affinity import translate
from shapely import STRtree


def get_polygon_points(polygon: Polygon) -> list[Point]:
    x, y = polygon.exterior.coords.xy
    return [Point(xx, yy) for xx, yy in zip(x, y)]

def get_polygon_segments(polygon: Polygon) -> list[LineString]:
    points = get_polygon_points(polygon)
    return [
        LineString([[p1.x, p1.y], [p2.x, p2.y]]) for p1, p2 in zip(points[:-1], points[1:])
    ]

def get_intersecting_points(polygon: Polygon) -> list[Point]:
    segments = get_polygon_segments(polygon)
    tree = STRtree(segments)
    intersection_points = set()
    for i, seg in enumerate(segments):
        intersecting_segments_idx = tree.query(seg)
        for idx in intersecting_segments_idx:
            # the condition checks that intersecting segments are not consecutive segments
            if idx != i and idx != (i + 1) % len(segments) and idx != (i - 1) % len(segments):
                intersection = seg.intersection(segments[idx])
                intersection_points.add((intersection.x, intersection.y))

    return [Point(*p) for p in intersection_points]


def self_touching_polygon(polygon: Polygon) -> bool:
    intersection_points = get_intersecting_points(polygon)
    segments = get_polygon_segments(polygon)
    tree = STRtree(segments)

    return all(len(tree.query(p)) == 3 for p in intersection_points)

def strictly_self_intersecting_polygon(polygon: Polygon) -> bool:
    return len(get_intersecting_points(polygon)) > 0 and not self_touching_polygon(polygon)

blue_coords = [[0, 0], [0, 100], [50, 80], [0, 60], [50, 30]]
blue_polygon = Polygon(blue_coords)

red_coords = [[0, 0], [0, 100], [50, 80], [-10, 60], [50, 30]]
red_polygon = Polygon(red_coords)

print(self_touching_polygon(red_polygon), strictly_self_intersecting_polygon(red_polygon))
print(self_touching_polygon(blue_polygon), strictly_self_intersecting_polygon(blue_polygon))
ftorre
  • 501
  • 1
  • 9