计算单位圆内两个多边形的交点
Calculating the Intersection of two polygons in a unit circle
给定一个单位圆和两个内接多边形，我该如何计算这两个多边形的并集交集 (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
三种方法
- 身材匀称
- 在 Python 而不是 Java 中重做 RaffleBuffle 方法（检查 his/her post 的描述）
- Monte Carlo模拟
身材匀称
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:
'''
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:
# 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
Monte Carlo模拟
- 在点的 min/max 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:
'''
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
假设多边形是非自相交的，即围绕圆的点的顺序是单调的，那么我相信有一种相对简单的方法可以确定 IoU 值，而不需要一般的形状包。
假设每个多边形的点围绕圆顺时针排列。我们可以通过增加 x 轴的角度 w.r.t 进行排序来确保这一点，或者如果我们发现带符号的区域是负数，则反转点。
将两个多边形的点合并到一个列表中，L，跟踪每个点属于哪个多边形。我们还需要能够为每个点确定原始多边形中的上一个点和下一个点。
按增加 x 轴的角度 w.r.t 排序 L。
如果输入多边形相交，则在迭代 L 时从一个多边形到另一个多边形的过渡次数将大于两次。
遍历 L。如果遇到属于不同多边形的连续点，则第一个点及其下一个点与第二个点及其前一个点之间的线的交点将属于两个多边形之间的交点。
将步骤 4 中识别的每个点添加到新的多边形中，I。 I的点会依次遇到
每个多边形的面积之和等于它们的并集加上交集的面积，因为这将被计算两次。
因此 IoU 的值将由 I 的面积除以两个多边形的面积之和减去 [=17= 的面积得到].
唯一需要的几何是使用 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);
}
给定一个单位圆和两个内接多边形，我该如何计算这两个多边形的并集交集 (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
三种方法
- 身材匀称
- 在 Python 而不是 Java 中重做 RaffleBuffle 方法（检查 his/her post 的描述）
- Monte Carlo模拟
身材匀称
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:
'''
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:
# 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
Monte Carlo模拟
- 在点的 min/max 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:
'''
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
假设多边形是非自相交的，即围绕圆的点的顺序是单调的，那么我相信有一种相对简单的方法可以确定 IoU 值，而不需要一般的形状包。
假设每个多边形的点围绕圆顺时针排列。我们可以通过增加 x 轴的角度 w.r.t 进行排序来确保这一点，或者如果我们发现带符号的区域是负数，则反转点。
将两个多边形的点合并到一个列表中，
L，跟踪每个点属于哪个多边形。我们还需要能够为每个点确定原始多边形中的上一个点和下一个点。按增加 x 轴的角度 w.r.t 排序
L。如果输入多边形相交，则在迭代
L时从一个多边形到另一个多边形的过渡次数将大于两次。遍历
L。如果遇到属于不同多边形的连续点，则第一个点及其下一个点与第二个点及其前一个点之间的线的交点将属于两个多边形之间的交点。将步骤 4 中识别的每个点添加到新的多边形中，
I。I的点会依次遇到每个多边形的面积之和等于它们的并集加上交集的面积，因为这将被计算两次。
因此
IoU的值将由I的面积除以两个多边形的面积之和减去 [=17= 的面积得到].
唯一需要的几何是使用 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);
}