7

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?

10
  • Your rings do not close, so they are invalid.
    – Vince
    Commented Sep 15, 2021 at 15:35
  • Not sure what you mean...
    – mhdadk
    Commented Sep 15, 2021 at 15:56
  • I'm not clear on how this question differs from gis.stackexchange.com/questions/379207/…?
    – Ian Turton
    Commented 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.
    – mhdadk
    Commented Sep 15, 2021 at 16:30
  • 1
    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
    – Ian Turton
    Commented Sep 16, 2021 at 13:14

2 Answers 2

5

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
   */
  public void linkResultDirectedEdges()

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.

5

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.

2
  • 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.
    – mhdadk
    Commented Apr 15, 2022 at 12:03
  • 1
    There is no paper specifically on the JTS approach - the source code is the reference!
    – dr_jts
    Commented Apr 16, 2022 at 2:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.