Shapely Algorithm – What Algorithm Does Shapely Use to Check if Two Polygons Intersect

algorithmintersectpolygonpythonshapely

What is the algorithm that Shapely used to check if two polygons intersect?

from shapely.geometry import Polygon
p1 = Polygon([(0,0), (3,0), (3,1), (1,1), (1,2), (3,2), (3,3), (0,3)])
p2 = Polygon([(4,0), (5,0), (5,1.5), (2,1.5), (2,1.2), (4,1.2)])
print(p1.intersects(p2))

I had read the source code of Shapely but didn't find the implementation.

Best Answer

A quick search of the Shapely code base leads to impl.py which:

"""Implementation of the intermediary layer between Shapely and GEOS This is layer number 2 from the list below.

  1. geometric objects: the Python OO API.
  2. implementation map: an abstraction that permits different backends.
  3. backend: callable objects that take Shapely geometric objects as arguments and, with GEOS as a backend, translate them to C data structures.
  4. GEOS library: algorithms implemented in C++. Shapely 1.2 includes a GEOS backend and it is the default. """

So the next place to look is GEOS and Geometry.ccp and it's intersection method. This is a wrapper for HeuristicOverlay which calls OverlayNGRobust::overlay which is where we start to see references to JTS, where the API docs state:

The obvious naive algorithm for intersection detection (comparing every segment with every other) has unacceptably slow performance. There is a large literature of faster algorithms for intersection detection. Unfortunately, many of them involve substantial code complexity. JTS tries to balance code simplicity with performance gains. It uses some simple techniques to produce substantial performance gains for common types of input data.

After some more poking around in the JTS code we eventually come to RelateComputer which calculates the IntersectionMatrix between two geometries (which includes intersection).

Any deeper than that and you will need to wait for Dr JTS to answer.

Related Question