计算单位圆中两个多边形的交点

Question

给定一个单位圆和两个内接多边形，如何计算两个多边形的并集交点（IoU）

假设我有一个关于单位圆的多边形坐标的Numpy数组，如下所示：

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),
    ]

预期IoU输出应为：0.124

Answer 1

假设多边形是非自交的，即圆周围点的顺序是单调的，那么我相信有一种相对简单的方法来确定IoU值，而不需要一般的形状包

1. 假设每个多边形的点围绕圆顺时针排列。我们可以通过增加x轴的角度w.r.t或反转点（如果我们发现有符号区域为负数）来进行排序来确保这一点

2. 将两个多边形中的点合并到一个列表中，L，跟踪每个点属于哪个多边形。我们还需要能够为每个点确定原始多边形中的上一点和下一点

3. 通过增加x轴的w.r.t角度对L进行排序

4. 如果输入多边形相交，则在通过L迭代时从一个多边形到另一个多边形的过渡次数将大于两个

5. 在L中迭代。如果遇到属于不同多边形的连续点，则第一个点及其下一个点与第二个点及其上一个点之间的直线交点将属于两个多边形之间的交点

6. 将步骤4中标识的每个点添加到新多边形I。将按顺序遇到I点

7. 每个多边形的面积之和等于它们的并集加上交点的面积，因为这将被计算两次

8. 因此IoU的值将由I的面积除以两个多边形的面积之和减去I的面积得出

所需的唯一几何图形是使用Shoelace formula计算简单多边形的面积，并根据步骤5的要求确定两条线段之间的交点

下面是一些Java代码（Ideone）来说明-您可能可以在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);

输出：

0.12403616470027268

和支持代码：

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);    
}

Answer 2

三种方法

• 匀称
• 用Python而不是Java重做RaffleBuffle方法（查看他/她的帖子以了解描述）
• 蒙特卡罗模拟

匀称

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

Python中的RaffleBuffle方法

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))
    

测试

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

蒙特卡罗模拟

• 在点的最小/最大x和y范围内生成随机点
• 计算任一多边形中的点数（即并集）
• 计算两个多边形中的点数（即交点）
• 交叉点的点数与并集的点数之比就是答案

代码

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))  

测试

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