Given a unit circle and two polygons inscribed, how can I go about calculating the Intersection over Union(IoU) of the two polygons?
Assume I have a Numpy array of the polygon coordinates with respect to the unit circle as:
Polygon 1: [
(-0.708, 0.707),
(0.309, -0.951),
(0.587, -0.809),
]
Polygon 2: [
(1, 0),
(0, 1),
(-1, 0),
(0, -1),
(0.708, -0.708),
]
The expected IoU output should be: 0.124
Three Approaches
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
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))
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)]
intersection_union(p1, p2)
# Out: 0.12418051627698147
Assuming polygons are non self-intersecting, i.e. the ordering of the points around the circle is monotonic, then I believe there is a relatively simple way to determine the IoU value without requiring a general shape package.
Assume the points of each polygon are ordered clockwise around the circle. We can ensure this by sorting by increasing angle w.r.t the x-axis or reversing the points if we find the signed area is negative.
Combine points from both polygons into a single list, L, keeping track of which polygon each point belongs to. We also need to be able to determine the previous and next point in the original polygon for each point.
Sort L by increasing angle w.r.t the x-axis.
If the input polygons intersect the number of transitions from one polygon to the other while iterating through L will be greater than two.
Iterate through L. If successive points are encountered belonging to different polygons then the intersection of the lines between the 1st point and its next point and the 2nd point and its previous point will belong to the intersection between the two polygons.
Add each point identified in step 4 to a new polygon, I. Points of I will be encountered in order.
The sum of the areas of each polygon will be equal to their union plus the area of the intersection, since this will be counted twice.
The value of IoU will therefore be given by the area of I divided by the sum of the areas of the two polygons minus the area of I.
The only geometry required is to calculate the area of a simple polygon using the Shoelace formula, and to determine the point of intersection between two line segments, required by step 5.
Here's some Java code (Ideone) to illustrate - you can probably make it a lot more compact in Python.
double[][] coords = {{-0.708, 0.707, 0.309, -0.951, 0.587, -0.809},
{1, 0, 0, 1, -1, 0, 0, -1, 0.708, -0.708}};
double areaSum = 0;
List<CPoint> pts = new ArrayList<>();
for(int p=0; p<coords.length; p++)
{
List<CPoint> poly = new ArrayList<>();
for(int j=0; j<coords[p].length; j+=2)
{
poly.add(new CPoint(p, coords[p][j], coords[p][j+1]));
}
double area = area(poly);
if(area < 0)
{
area = -area;
Collections.reverse(poly);
}
areaSum += area;
pts.addAll(poly);
int n = poly.size();
for(int i=0, j=n-1; i<n; j=i++)
{
poly.get(i).prev = poly.get(j);
poly.get(j).next = poly.get(i);
}
}
pts.sort((a, b) -> Double.compare(a.theta, b.theta));
List<Point2D> intersections = new ArrayList<>();
int n = pts.size();
for(int i=0, j=n-1; i<n; j=i++)
{
if(pts.get(j).id != pts.get(i).id)
{
intersections.add(intersect(pts.get(j), pts.get(j).next, pts.get(i).prev, pts.get(i)));
}
}
double areaInt = area(intersections);
double iou = areaInt/(areaSum - areaInt);
System.out.println(iou);
Output:
0.12403616470027268
And supporting code:
static class CPoint extends Point2D.Double
{
int id;
double theta;
CPoint prev, next;
public CPoint(int id, double x, double y)
{
super(x, y);
this.id = id;
theta = Math.atan2(y, x);
if(theta < 0) theta = 2*Math.PI + theta;
}
}
static double area(List<? extends Point2D> poly)
{
double area = 0;
for(int i=0, j=poly.size()-1; i<poly.size(); j=i++)
area += (poly.get(j).getX() * poly.get(i).getY()) - (poly.get(i).getX() * poly.get(j).getY());
return Math.abs(area)/2;
}
// https://rosettacode.org/wiki/Find_the_intersection_of_two_lines#Java
static Point2D intersect(Point2D p1, Point2D p2, Point2D p3, Point2D p4)
{
double a1 = p2.getY() - p1.getY();
double b1 = p1.getX() - p2.getX();
double c1 = a1 * p1.getX() + b1 * p1.getY();
double a2 = p4.getY() - p3.getY();
double b2 = p3.getX() - p4.getX();
double c2 = a2 * p3.getX() + b2 * p3.getY();
double delta = a1 * b2 - a2 * b1;
return new Point2D.Double((b2 * c1 - b1 * c2) / delta, (a1 * c2 - a2 * c1) / delta);
}