# How is polygon intersection implemented in JTS/Shapely?

I am using Shapely, which from the questions

uses the JTS as a backend. However, none of the answers to these questions specify which algorithm is used in JTS to compute the intersection of two polygons. I've checked the question Are there any references which describe the algorithms used in JTS? on the JTS website, but with no luck.

More specifically, if I write

``````from shapely.geometry import Polygon
import math

sqrt2 = math.sqrt(2)

# this is a square
A = Polygon([(1,1),(1,-1),(-1,-1),(-1,1)])

# this is the square A rotated by 45 degrees counter-clockwise
B = Polygon([(0,sqrt2),(sqrt2,0),(0,-sqrt2),(-sqrt2,0)])

# intersection of A and B, which is an octagon
C = A.intersection(B)
``````

which algorithm is used by JTS to compute the polygon `C`? Is it the Sutherland-Hodgman algorithm?

• Your rings do not close, so they are invalid. Sep 15, 2021 at 15:35
• Not sure what you mean... Sep 15, 2021 at 15:56
• I'm not clear on how this question differs from gis.stackexchange.com/questions/379207/…? Sep 15, 2021 at 16:07
• @IanTurton the answer to that question is concerned with the `intersects` method and not the `intersection` method, which I am asking about. I've checked the code related to dr_jts's comments, but it still does not answer my question. Sep 15, 2021 at 16:30
• From the JTS faq - E. Chan, J. Ng. A General and Efficient Implementation of Geometric Operators and Predicates; Proceedings of the 5th International Symposium on Advances in Spatial Databases, 1997. seems to be the answer you are looking for. Check citeseer for a copy Sep 16, 2021 at 13:14

From @user30184's comment, I checked the code and first found this:

``````  public static Geometry overlay(Geometry geom0, Geometry geom1,
int opCode, PrecisionModel pm)
{
OverlayNG ov = new OverlayNG(geom0, geom1, pm, opCode);
Geometry geomOv = ov.getResult();
return geomOv;
}
``````

I then checked the `getResult()` method and found this:

``````    else {
// handle case where both inputs are formed of edges (Lines and Polygons)
result = computeEdgeOverlay();
}
``````

Checking the `computeEdgeOverlay` function led me to this:

``````    return extractResult(opCode, graph);
``````

I then checked the `extractResult` function and found this:

``````  private Geometry extractResult(int opCode, OverlayGraph graph) {
boolean isAllowMixedIntResult = ! isStrictMode;

//--- Build polygons
List<OverlayEdge> resultAreaEdges = graph.getResultAreaEdges();
PolygonBuilder polyBuilder = new PolygonBuilder(resultAreaEdges, geomFact);
List<Polygon> resultPolyList = polyBuilder.getPolygons();
boolean hasResultAreaComponents = resultPolyList.size() > 0;
``````

Finally, I checked the `getResultAreaEdges` method and found the following comment after it:

``````  /**
* Traverse the star of DirectedEdges, linking the included edges together.
* To link two dirEdges, the <code>next</code> pointer for an incoming dirEdge
* is set to the next outgoing edge.
* <p>
* DirEdges are only linked if:
* <ul>
* <li>they belong to an area (i.e. they have sides)
* <li>they are marked as being in the result
* </ul>
* <p>
* Edges are linked in CCW order (the order they are stored).
* This means that rings have their face on the Right
* (in other words,
* the topological location of the face is given by the RHS label of the DirectedEdge)
* <p>
* PRECONDITION: No pair of dirEdges are both marked as being in the result
*/
``````

This comment includes the words `CCW order`, `RHS label`, and `in the result`, which highly suggests that the algorithm that is being used can be found in the paper

"An Edge Labeling Approach to Concave Polygon Clipping" by Klamer Schutte (1995)

since these terms explicitly appear in this paper. This is also the second paper mentioned in the JTS FAQ question Are there any references which describe the algorithms used in JTS?. However, there is no way to say for sure that this is true, and I believe that special cases are handled using the algorithms found in the paper

"A General and Efficient Implementation of Geometric Operators and Predicates" by Edward P. F. Chan and Jimmy N. H. Ng (1997)

which is also mentioned in the JTS FAQ.

The overlay algorithm used in JTS and GEOS is my own design. I developed it before I saw either of those papers referenced above. But the techniques are similar (there's a fairly small design space for geometry overlay).

It's worth noting that although the concept of computing the noded and labelled arrangement of the input geometries is relatively common, the approaches used to provide efficient computation and geometrical robustness are important as well, and are often glossed over in academic implementations.

• Thanks a lot for this answer. Is there a paper/reference I could read to understand the overlay algorithm you designed a bit better? I am not sure which of the papers in my answer are relevant. Apr 15, 2022 at 12:03
• There is no paper specifically on the JTS approach - the source code is the reference! Apr 16, 2022 at 2:05