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

三种方法

匀称
在 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

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

问题描述

解决方案

推荐阅读