How to Intersect Two Polygons

A simple algorithm for polygon intersection

I understand the original poster was looking for a simple solution, but unfortunately there really is no simple solution.

Nevertheless, I've recently created an open-source freeware clipping library (written in Delphi, C++ and C#) which clips all kinds of polygons (including self-intersecting ones). This library is pretty simple to use: http://sourceforge.net/projects/polyclipping/ .

How to intersect two polygons?

What I think you should do

Do not attempt to do this yourself if you can possibly help it. Instead, use one of the many available polygon intersection algorithms that already exist.

I was strongly considering the following codebase on the strength of their demonstration code and the fact that they mentioned their handling of most/all of the weird cases. You would need to donate an amount (of you/your company's choice) if you use it commercially, but it's worth it to get a robust version of this kind of code.

http://www.cs.man.ac.uk/~toby/gpc/

What I actually did was to use a polygon-intersection algorithm that is part of the Java2D libraries. You can possibly find something similar in MS's own C# libraries to use.

There are other options out there as well; look for "polygon clipper" or "polygon clipping", since the same basic algorithms that handle polygon intersection also tend to be usable for the general clipping cases.

Once you actually have a polygon clipping library, you just need to subtract polygon B from polygon A to get your first piece of output, and intersect polygons A and B to get your second piece of output.

How to roll your own, for the hopelessly masochistic

When I was considering rolling my own, I found the Weiler-Atherton algorithm to have the most potential for general polygon-cutting. I used the following as a reference:

http://cs1.bradley.edu/public/jcm/weileratherton.html

http://en.wikipedia.org/wiki/Weiler-Atherton

The details, as they say, are too dense to include here, but I have no doubt that you'll be able to find references on Weiler-Atherton for years to come. Essentially, you split all the points into those that are entering the final polygon or exiting the final polygon, then you form a graph out of all the points, and then walk the graph in the appropriate directions in order to extract all the polygon pieces you want. By changing the way you define and treat the "entering" and "exiting" polygons, you can achieve several possible polygon intersections (AND, OR, XOR, etc.).

It's actually fairly implementable, but like with any computational geometry code, the devil is in the degeneracies.

How to intersect more than 2 polygons using Turf.js?

turf.intersect() can only intersect 2 polygons.

If you want to intersect multiple polygons you can intersect each polygon with each other polygon with turf.intersect() and then combine the results using turf.combine().

Here is some example code:

const polygonA = ...;
const polygonB = ...;
const polygonC = ...;

const allIntersections = {
  type: 'FeatureCollections',
  features: [
    turf.intersect(polygonA, polygonB),
    turf.intersect(polygonA, polygonC),
    turf.intersect(polygonB, polygonC),
  ],
};

const combinedIntersection = turf.combine(allIntersections);

Calculating the Intersection of two polygons in a unit circle

Three Approaches

Shapely
Redo RaffleBuffle approach in Python rather than Java (check his/her post for description)
Monte Carlo Simulation

Shapely

from shapely.geometry import Polygon

p1 = [  (-0.708, 0.707),
        (0.309, -0.951),
        (0.587, -0.809)]  

p2 = [  (1, 0),
        (0, 1),
        (-1, 0),
        (0, -1),
        (0.708, -0.708)]

# Create polygons from coordinates
poly1 = Polygon(p1)
poly2 = Polygon(p2)

# ratio of intersection area to union area
print(poly1.intersection(poly2).area/poly1.union(poly2).area)
# Output: 0.12403616470027264

RaffleBuffle Approach in Python

import math
import numpy as np

class Point:
    def __init__(self, x, y, id = None):
        self.x = x
        self.y = y
        self.id = id
        self.prev = None
        self.next = None
        
    def __repr__(self):
        result = f"{self.x} {self.y} {self.id}"
        result += f" Prev: {self.prev.x} {self.prev.y} {self.prev.id}" if self.prev else ""
        result += f" Next: {self.next.x} {self.next.y} {self.next.id}" if self.next else ""
        return result
        
class Poly:
    def __init__(self, pts):
        self.pts = [p for p in pts]
        
        ' Sort polynomial coordinates based upon angle and radius in clockwize direction'
        # Origin is the centroid of points
        origin = Point(*[sum(pt.x for pt in self.pts)/len(self.pts), sum(pt.y for pt in self.pts)/len(self.pts)])

        # Sort points by incresing angle around centroid based upon angle and distance to centroid
        self.pts.sort(key=lambda p: clockwiseangle_and_distance(p, origin))
        
    def __add__(self, other):
        ' Overload for adding two polynomials '
        return Poly(self.pts + other.pts)
        
    def assign_chain(self):
        ' Assign prev and next '
        n = len(self.pts)
        for i in range(n):
            self.pts[i].next = self.pts[ (i + 1) % n]
            self.pts[i].prev = self.pts[(i -1) % n]
        return self
            
    def area(self):
        '''
            area enclosed by polynomial coordinates '
        
            Source: https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates
        '''
        x = np.array([p.x for p in self.pts])
        y = np.array([p.y for p in self.pts])
        return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
        
    def intersection(self):
        ' Intersection of coordinates with different ids '
        pts = self.pts
        
        n = len(pts)

        intersections = []
        for i in range(n):
            j = (i -1) % n
            if pts[j].id != pts[i].id:
                get_j = [pts[j], pts[j].next]
                get_i = [pts[i].prev, pts[i]]
                pt_intersect = line_intersect(get_j, get_i)
                if pt_intersect:
                    intersections.append(pt_intersect)
                
        return intersections
        
            
    def __str__(self):
        return '\n'.join(str(pt) for pt in self.pts)
    
    def __repr__(self):
        return '\n'.join(str(pt) for pt in self.pts)
    
def clockwiseangle_and_distance(point, origin):
    # Source: https://stackoverflow.com/questions/41855695/sorting-list-of-two-dimensional-coordinates-by-clockwise-angle-using-python
    # Vector between point and the origin: v = p - o
    vector = [point.x-origin.x, point.y-origin.y]
    refvec = [0, 1]
    
    # Length of vector: ||v||
    lenvector = math.hypot(vector[0], vector[1])
    # If length is zero there is no angle
    if lenvector == 0:
        return -math.pi, 0
    # Normalize vector: v/||v||
    normalized = [vector[0]/lenvector, vector[1]/lenvector]
    dotprod  = normalized[0]*refvec[0] + normalized[1]*refvec[1]     # x1*x2 + y1*y2
    diffprod = refvec[1]*normalized[0] - refvec[0]*normalized[1]     # x1*y2 - y1*x2
    angle = math.atan2(diffprod, dotprod)
    # Negative angles represent counter-clockwise angles so we need to subtract them 
    # from 2*pi (360 degrees)
    if angle < 0:
        return 2*math.pi+angle, lenvector
    # I return first the angle because that's the primary sorting criterium
    # but if two vectors have the same angle then the shorter distance should come first.
    return angle, lenvector

def line_intersect(segment1, segment2):
    """ returns a (x, y) tuple or None if there is no intersection 
    
        segment1 and segment2 are two line segements
        
        specified by their starting/ending points
        
        Source: https://rosettacode.org/wiki/Find_the_intersection_of_two_lines#Python
        
    """
    Ax1, Ay1 = segment1[0].x, segment1[0].y     # starting point in line segment 1
    Ax2, Ay2 = segment1[1].x, segment1[1].y     # ending point in line segment 1
    Bx1, By1 = segment2[0].x, segment2[0].y     # starting point in line segment 2
    Bx2, By2 = segment2[1].x, segment2[1].y     # ending point in line segment 2
    
    d = (By2 - By1) * (Ax2 - Ax1) - (Bx2 - Bx1) * (Ay2 - Ay1)
    
    if d:
        uA = ((Bx2 - Bx1) * (Ay1 - By1) - (By2 - By1) * (Ax1 - Bx1)) / d
        uB = ((Ax2 - Ax1) * (Ay1 - By1) - (Ay2 - Ay1) * (Ax1 - Bx1)) / d
    else:
        return
    if not(0 <= uA <= 1 and 0 <= uB <= 1):
        return
    x = Ax1 + uA * (Ax2 - Ax1)
    y = Ay1 + uA * (Ay2 - Ay1)
 
    return Point(x, y, None)

    
def polygon_iou(coords1, coords2):
    '''
        Calculates IoU of two 2D polygons based upon coordinates
    '''
    # Make polynomials ordered clockwise and assign ID (0 and 1)
    poly1 = Poly(Point(*p, 0) for p in coords1)   # counter clockwise with ID 0
    poly2 = Poly(Point(*p, 1) for p in coords2)   # counter clockwise with ID 1
    
    # Assign previous and next coordinates in polynomial chain
    poly1.assign_chain()
    poly2.assign_chain()
    
    # Polygon areas
    area1 = poly1.area()
    area2 = poly2.area()
            
    # Combine both polygons into one counter clockwise
    poly = poly1 + poly2
    
    # Get interesections
    intersections = poly.intersection()
    
    # IoU based upon intersection and sum of areas
    if intersections:
        area_intersection = Poly(intersections).area()
        result = area_intersection/(area1 + area2 - area_intersection)
    else:
        result = 0.0
        
    return result

print(polygon_iou(p1, p2))

Test

p1 = [  (-0.708, 0.707),
        (0.309, -0.951),
        (0.587, -0.809)]  

p2 = [  (1, 0),
        (0, 1),
        (-1, 0),
        (0, -1),
        (0.708, -0.708)]

print(polygon_iou(p1, p2))
# Output: 0.12403616470027277

Monte Carlo Simulation

Generate random points in the min/max x and y range of the points
Count the number of points in either polygon (i.e. union)
Count the number of points in both polygon (i.e. intersection)
Ratio of the number of points in the intersection to the number in the union is the answer

Code

import math
from random import uniform

def ray_tracing_method(x, y, poly):
    '''
       Determines if point x, y is inside polygon poly

       Source: "What's the fastest way of checking if a point is inside a polygon in python"
             at URL: https://stackoverflow.com/questions/36399381/whats-the-fastest-way-of-checking-if-a-point-is-inside-a-polygon-in-python

    '''
    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

def intersection_union(p1, p2, num_throws = 1000000):
    '''
        Computes the intersection untion for polygons p1, p2
        (assuming that p1, and p2 are inside at polygon of radius 1)
    '''
    # Range of values
    p = p1 + p2
    xmin = min(x[0] for x in p)
    xmax = max(x[0] for x in p)
    ymin = min(x[1] for x in p)
    ymax = max(x[1] for x in p)

    # Init counts
    in_union = 0
    in_intersect = 0
    throws = 0

    while (throws < num_throws):
        # Choose random x, y position in rectangle
        x_pos = uniform (xmin, xmax)
        y_pos = uniform (ymin, ymax)

        # Only test points inside unit circle
        # Check if points are inside p1 & p2
        in_p1 = ray_tracing_method(x_pos, y_pos, p1)
        in_p2 = ray_tracing_method(x_pos, y_pos, p2)

        if in_p1 or in_p2:
            in_union += 1             # in union

        if in_p1 and in_p2:
            in_intersect += 1         # in intersetion

        throws += 1   
        
    return in_intersect/in_union
        
print(intersection_union(p1, p2))