在 python 中找到直线和圆的交点的最有效方法是什么？
What is most efficient way to find the intersection of a line and a circle in python?
我有一个由很多点组成的多边形。我想找到多边形和圆的交点。提供[x0,y0]的圆心和r0的半径，我写了一个粗略的函数来简单地求解圆和直线的二次方程。但是将多边形的每条线段逐条求交点效率如何呢？ 有没有更高效的方法？
我知道 sympy 已经提供了获取不同几何图形交集的功能。而且如果我想处理很多多边形，那么与通过我自己的函数计算它相比，像 sympy 这样的外部库的效率如何？
def LineIntersectCircle(p,lsp,lep):
# p is the circle parameter, lsp and lep is the two end of the line
x0,y0,r0 = p
x1,y1 = lsp
x2,y2 = esp
if x1 == x2:
if abs(r0) >= abs(x1 - x0):
p1 = x1, y0 - sqrt(r0**2 - (x1-x0)**2)
p2 = x1, y0 + sqrt(r0**2 - (x1-x0)**2)
inp = [p1,p2]
# select the points lie on the line segment
inp = [p for p in inp if p[1]>=min(y1,y2) and p[1]<=max(y1,y2)]
else:
inp = []
else:
k = (y1 - y2)/(x1 - x2)
b0 = y1 - k*x1
a = k**2 + 1
b = 2*k*(b0 - y0) - 2*x0
c = (b0 - y0)**2 + x0**2 - r0**2
delta = b**2 - 4*a*c
if delta >= 0:
p1x = (-b - sqrt(delta))/(2*a)
p2x = (-b + sqrt(delta))/(2*a)
p1y = k*x1 + b0
p2y = k*x2 + b0
inp = [[p1x,p1y],[p2x,p2y]]
# select the points lie on the line segment
inp = [p for p in inp if p[0]>=min(x1,x2) and p[0]<=max(x1,x2)]
else:
inp = []
return inp
我认为您用于查找两个交点坐标的公式无法进一步优化。唯一的改进（这在数值上很重要）是区分两种情况：|x_2-x_1| >= |y_2-y_1| 和 |x_2-x1| < |y_2-y1| 以便数量 k 始终在 -1 和 1 之间（在您的计算中您可以得到如果 |x_2-x_1| 非常小，则数值错误非常高）。您可以交换 x-s 和 y-s 以将一种情况减少到另一种情况。
您还可以进行初步检查：如果两个端点都在圆的内部，则没有交点。通过计算从点到圆心的平方距离，这变成了一个不使用平方根函数的简单公式。另一个检查："whether the line is outside the circle" 已经在您的代码中实现并且对应于 delta < 0。如果您有很多小段，这两个检查在大多数情况下应该给出一个快捷答案（没有交集）。
一个低成本的空间分区可能是把飞机分成9块
这是一张蹩脚的图表。想象一下，如果你愿意的话，这些线刚好接触到圆圈。
| |
__|_|__
__|O|__
| |
| |
我们感兴趣的区域中有 8 个在圆圈周围。中间的正方形对于便宜的测试没有太大用处，但是你可以在圆圈内放置一个 r/sqrt(2) 的正方形，这样它的角就可以接触到圆圈。
让我们标记区域
A |B| C
__|_|__
D_|O|_E
| |
F |G| H
中间的 r/sqrt(2) 的正方形我们称之为 J
我们将调用图中所示中心方块中不在 J、Z
中的点集
用字母代码标记多边形的每个顶点。
现在我们可以很快看到
AA => Outside
AB => Outside
AC => Outside
...
AJ => Intersects
BJ => Intersects
...
JJ => Inside
这可以变成查找table
因此，根据您的数据集，您可能已经为自己节省了大量工作。但是，任何具有 Z 端点的内容都需要进行测试。
我想您的问题可能是理论上如何以最快的方式执行此操作。但是如果你想快速做到这一点，你真的应该使用用 C/C++.
编写的东西
我已经习惯了 Shapely，所以我将提供一个示例来说明如何使用此库执行此操作。 python 有很多几何库。我会在这个答案的末尾列出它们。
from shapely.geometry import LineString
from shapely.geometry import Point
p = Point(5,5)
c = p.buffer(3).boundary
l = LineString([(0,0), (10, 10)])
i = c.intersection(l)
print i.geoms[0].coords[0]
(2.8786796564403576, 2.8786796564403576)
print i.geoms[1].coords[0]
(7.121320343559642, 7.121320343559642)
Shapely 中有点违反直觉的是，圆是带缓冲区的点的边界。这就是为什么我 p.buffer(3).boundry
此外，交集 i 是一个几何形状的列表，在这种情况下它们都是点，这就是为什么我必须从 i.geoms[]
中得到它们
有 another Whosebug question 为感兴趣的人详细介绍了这些库。
编辑因为评论：
Shapely 基于 GEOS (trac.osgeo.org/geos)，它是用 C++ 构建的，比您在 python 中原生编写的任何东西都要快得多。 SymPy 似乎是基于 mpmath (mpmath.org)，它似乎也是 python，但似乎集成了很多非常复杂的数学。自己实施可能需要大量工作，而且可能不如 GEOS C++ 实施快。
这是一个计算圆与由两个 (x, y) 点定义的直线或线段的交点的解决方案：
def circle_line_segment_intersection(circle_center, circle_radius, pt1, pt2, full_line=True, tangent_tol=1e-9):
""" Find the points at which a circle intersects a line-segment. This can happen at 0, 1, or 2 points.
:param circle_center: The (x, y) location of the circle center
:param circle_radius: The radius of the circle
:param pt1: The (x, y) location of the first point of the segment
:param pt2: The (x, y) location of the second point of the segment
:param full_line: True to find intersections along full line - not just in the segment. False will just return intersections within the segment.
:param tangent_tol: Numerical tolerance at which we decide the intersections are close enough to consider it a tangent
:return Sequence[Tuple[float, float]]: A list of length 0, 1, or 2, where each element is a point at which the circle intercepts a line segment.
Note: We follow: http://mathworld.wolfram.com/Circle-LineIntersection.html
"""
(p1x, p1y), (p2x, p2y), (cx, cy) = pt1, pt2, circle_center
(x1, y1), (x2, y2) = (p1x - cx, p1y - cy), (p2x - cx, p2y - cy)
dx, dy = (x2 - x1), (y2 - y1)
dr = (dx ** 2 + dy ** 2)**.5
big_d = x1 * y2 - x2 * y1
discriminant = circle_radius ** 2 * dr ** 2 - big_d ** 2
if discriminant < 0: # No intersection between circle and line
return []
else: # There may be 0, 1, or 2 intersections with the segment
intersections = [
(cx + (big_d * dy + sign * (-1 if dy < 0 else 1) * dx * discriminant**.5) / dr ** 2,
cy + (-big_d * dx + sign * abs(dy) * discriminant**.5) / dr ** 2)
for sign in ((1, -1) if dy < 0 else (-1, 1))] # This makes sure the order along the segment is correct
if not full_line: # If only considering the segment, filter out intersections that do not fall within the segment
fraction_along_segment = [(xi - p1x) / dx if abs(dx) > abs(dy) else (yi - p1y) / dy for xi, yi in intersections]
intersections = [pt for pt, frac in zip(intersections, fraction_along_segment) if 0 <= frac <= 1]
if len(intersections) == 2 and abs(discriminant) <= tangent_tol: # If line is tangent to circle, return just one point (as both intersections have same location)
return [intersections[0]]
else:
return intersections
我有一个由很多点组成的多边形。我想找到多边形和圆的交点。提供[x0,y0]的圆心和r0的半径，我写了一个粗略的函数来简单地求解圆和直线的二次方程。但是将多边形的每条线段逐条求交点效率如何呢？ 有没有更高效的方法？
我知道 sympy 已经提供了获取不同几何图形交集的功能。而且如果我想处理很多多边形，那么与通过我自己的函数计算它相比，像 sympy 这样的外部库的效率如何？
def LineIntersectCircle(p,lsp,lep):
# p is the circle parameter, lsp and lep is the two end of the line
x0,y0,r0 = p
x1,y1 = lsp
x2,y2 = esp
if x1 == x2:
if abs(r0) >= abs(x1 - x0):
p1 = x1, y0 - sqrt(r0**2 - (x1-x0)**2)
p2 = x1, y0 + sqrt(r0**2 - (x1-x0)**2)
inp = [p1,p2]
# select the points lie on the line segment
inp = [p for p in inp if p[1]>=min(y1,y2) and p[1]<=max(y1,y2)]
else:
inp = []
else:
k = (y1 - y2)/(x1 - x2)
b0 = y1 - k*x1
a = k**2 + 1
b = 2*k*(b0 - y0) - 2*x0
c = (b0 - y0)**2 + x0**2 - r0**2
delta = b**2 - 4*a*c
if delta >= 0:
p1x = (-b - sqrt(delta))/(2*a)
p2x = (-b + sqrt(delta))/(2*a)
p1y = k*x1 + b0
p2y = k*x2 + b0
inp = [[p1x,p1y],[p2x,p2y]]
# select the points lie on the line segment
inp = [p for p in inp if p[0]>=min(x1,x2) and p[0]<=max(x1,x2)]
else:
inp = []
return inp
我认为您用于查找两个交点坐标的公式无法进一步优化。唯一的改进（这在数值上很重要）是区分两种情况：|x_2-x_1| >= |y_2-y_1| 和 |x_2-x1| < |y_2-y1| 以便数量 k 始终在 -1 和 1 之间（在您的计算中您可以得到如果 |x_2-x_1| 非常小，则数值错误非常高）。您可以交换 x-s 和 y-s 以将一种情况减少到另一种情况。
您还可以进行初步检查：如果两个端点都在圆的内部，则没有交点。通过计算从点到圆心的平方距离，这变成了一个不使用平方根函数的简单公式。另一个检查："whether the line is outside the circle" 已经在您的代码中实现并且对应于 delta < 0。如果您有很多小段，这两个检查在大多数情况下应该给出一个快捷答案（没有交集）。
一个低成本的空间分区可能是把飞机分成9块
这是一张蹩脚的图表。想象一下，如果你愿意的话，这些线刚好接触到圆圈。
| | __|_|__ __|O|__ | | | |
我们感兴趣的区域中有 8 个在圆圈周围。中间的正方形对于便宜的测试没有太大用处，但是你可以在圆圈内放置一个 r/sqrt(2) 的正方形，这样它的角就可以接触到圆圈。
让我们标记区域
A |B| C __|_|__ D_|O|_E | | F |G| H
中间的 r/sqrt(2) 的正方形我们称之为 J
我们将调用图中所示中心方块中不在 J、Z
用字母代码标记多边形的每个顶点。
现在我们可以很快看到
AA => Outside AB => Outside AC => Outside ... AJ => Intersects BJ => Intersects ... JJ => Inside
这可以变成查找table
因此，根据您的数据集，您可能已经为自己节省了大量工作。但是，任何具有 Z 端点的内容都需要进行测试。
我想您的问题可能是理论上如何以最快的方式执行此操作。但是如果你想快速做到这一点，你真的应该使用用 C/C++.
编写的东西我已经习惯了 Shapely，所以我将提供一个示例来说明如何使用此库执行此操作。 python 有很多几何库。我会在这个答案的末尾列出它们。
from shapely.geometry import LineString
from shapely.geometry import Point
p = Point(5,5)
c = p.buffer(3).boundary
l = LineString([(0,0), (10, 10)])
i = c.intersection(l)
print i.geoms[0].coords[0]
(2.8786796564403576, 2.8786796564403576)
print i.geoms[1].coords[0]
(7.121320343559642, 7.121320343559642)
Shapely 中有点违反直觉的是，圆是带缓冲区的点的边界。这就是为什么我 p.buffer(3).boundry
此外，交集 i 是一个几何形状的列表，在这种情况下它们都是点，这就是为什么我必须从 i.geoms[]
有 another Whosebug question 为感兴趣的人详细介绍了这些库。
编辑因为评论：
Shapely 基于 GEOS (trac.osgeo.org/geos)，它是用 C++ 构建的，比您在 python 中原生编写的任何东西都要快得多。 SymPy 似乎是基于 mpmath (mpmath.org)，它似乎也是 python，但似乎集成了很多非常复杂的数学。自己实施可能需要大量工作，而且可能不如 GEOS C++ 实施快。
这是一个计算圆与由两个 (x, y) 点定义的直线或线段的交点的解决方案：
def circle_line_segment_intersection(circle_center, circle_radius, pt1, pt2, full_line=True, tangent_tol=1e-9):
""" Find the points at which a circle intersects a line-segment. This can happen at 0, 1, or 2 points.
:param circle_center: The (x, y) location of the circle center
:param circle_radius: The radius of the circle
:param pt1: The (x, y) location of the first point of the segment
:param pt2: The (x, y) location of the second point of the segment
:param full_line: True to find intersections along full line - not just in the segment. False will just return intersections within the segment.
:param tangent_tol: Numerical tolerance at which we decide the intersections are close enough to consider it a tangent
:return Sequence[Tuple[float, float]]: A list of length 0, 1, or 2, where each element is a point at which the circle intercepts a line segment.
Note: We follow: http://mathworld.wolfram.com/Circle-LineIntersection.html
"""
(p1x, p1y), (p2x, p2y), (cx, cy) = pt1, pt2, circle_center
(x1, y1), (x2, y2) = (p1x - cx, p1y - cy), (p2x - cx, p2y - cy)
dx, dy = (x2 - x1), (y2 - y1)
dr = (dx ** 2 + dy ** 2)**.5
big_d = x1 * y2 - x2 * y1
discriminant = circle_radius ** 2 * dr ** 2 - big_d ** 2
if discriminant < 0: # No intersection between circle and line
return []
else: # There may be 0, 1, or 2 intersections with the segment
intersections = [
(cx + (big_d * dy + sign * (-1 if dy < 0 else 1) * dx * discriminant**.5) / dr ** 2,
cy + (-big_d * dx + sign * abs(dy) * discriminant**.5) / dr ** 2)
for sign in ((1, -1) if dy < 0 else (-1, 1))] # This makes sure the order along the segment is correct
if not full_line: # If only considering the segment, filter out intersections that do not fall within the segment
fraction_along_segment = [(xi - p1x) / dx if abs(dx) > abs(dy) else (yi - p1y) / dy for xi, yi in intersections]
intersections = [pt for pt, frac in zip(intersections, fraction_along_segment) if 0 <= frac <= 1]
if len(intersections) == 2 and abs(discriminant) <= tangent_tol: # If line is tangent to circle, return just one point (as both intersections have same location)
return [intersections[0]]
else:
return intersections