前言
您的问题是“一个点在正多边形内吗？”的专业化，我的意思是不是我们经常在 GIS 系统中发现的自相交或多多边形。并且由于您要求三角形，因此我假设多边形是凸的。
有趣的是，叉积是解决它的关键。当您在 2D 平面中处理向量时，叉积与该平面正交。要提取的有用信息是：它是向上还是向下？
Float Arithmetic 会发生错误，当叉积接近于零但不等于时，它变得很关键，然后它将有符号而不是 null。
要检查您的点是否在多边形内，只需检查边缘和点之间的所有叉积是否具有相同的符号，例如：sign(h x w) = -1。
同理，要检查多边形是否凸，就要检查所有连续边的叉积是否具有相同的符号，例如：sign(u x v) = -1。
实施
让我们构建一个小类来检查一个点是否在正凸多边形的内部（在边缘或外部）：
import numpy as np
class cpoly:
def __init__(self, points=[[0,0], [0,1], [1,0]], assert_convexity=True):
"""
Initialize 2D Polygon with a sequence of 2D points
"""
self._points = np.array(points)
assert self.p.shape[0] >= 3
assert self.p.shape[1] == 2
assert self.is_convex or not(assert_convexity)
@property
def n(self):
return self.p.shape[0]
@property
def p(self):
return self._points
@property
def is_convex(self):
"""
Check convexity of the polygon (operational for a non intersecting polygon)
"""
return self.contains()
def contains(self, p=None, debug=False, atol=2e-16):
"""
Check if a 2D convex polygon contains a point (also used to assess convexity)
Returns:
-1: Point is oustide the polygon
0: Point is close to polygon edge (epsilon ball)
+1: Point is inside the polygon
"""
s = None
c = False
n = self.n
for k in range(n):
# Vector Differences:
d1 = self.p[(k+1)%n,:] - self.p[k%n,:]
if p:
d2 = p - self.p[k%n,:]
else:
d2 = self.p[(k+2)%n,:] - self.p[(k+1)%n,:]
# Cross Product:
z = np.cross(d1, d2)
if np.allclose(z, 0, atol=atol):
s_ = 0
c = True
else:
s_ = np.sign(z)
# Debug Helper:
if debug:
print("k = %d, d1 = %s, d2 = %s, z = %.32f, s = %d" % (k, d1, d2, z, s_))
# Check if cross product sign change (excluded null, when point is colinear with the segment)
if s and (s_ != s) and not(s_ == 0):
# Nota: Integer are exact if float representable, therefore comparizons are correct
return -1
s = s_
if c:
return 0
else:
return 1
def plot(self, axe=None):
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
if not(axe):
fig, axe = plt.subplots()
axe.plot(self.p[:,0], self.p[:,1], 'x', markersize=10, label='Points $p_i$')
axe.add_patch(Polygon(self.p, alpha=0.4, label='Area'))
axe.set_xlabel("$x$")
axe.set_ylabel("$y$")
axe.set_title("Polygon")
axe.set_aspect('equal')
axe.legend(bbox_to_anchor=(1,1), loc='upper left')
axe.grid()
return axe.get_figure(), axe
这个类被初始化为一个二维点列表（默认是你的）。
p = cpoly()
在我的设置中，浮点精度大约是：
e = np.finfo(np.double).eps # 2.220446049250313e-16
我们创建一个试验数据集用于测试目的：
p = cpoly()
r = [
[0,0], [0,1], [1,0], # Polygon vertices
[0,0.5], [-e,0.6], [e,0.4], [0.1, 0.1], [1,1],
[0.5+e,0.5], [0.3-e,0.7], [0.7+e/10,0.3],
[0, 1.2], [1.2, 0.], # Those points make your logic fails
[0.2,0.8], [0.1,0.9],
[0.8+10*e,0.2], [0.9+10*e,0.1]
]
如果我们调整您的函数以获得合规的输出：
def isInsideTriangle(x,y):
if x == 0 or y == 0 or y == 1-x:
return 0
elif x > 1 or y > 1 or x < 0 or y < 0 or y > 1-x:
return -1
else:
return 1
然后我们检查试验点，看看我们的两个函数是否表现良好：
_, axe = p.plot()
cols = {-1: 'red', 0:'orange', 1:'green'}
for x in r:
q1 = p.contains(x)
q2 = isInsideTriangle(*x)
print(q1==q2)
axe.plot(*x, 'o', markersize=4, color=cols[q1])
使用此设置，所有点都被正确分类。但是你可以看到你的算法是有缺陷的。主要是以下几行：
if x == 0 or y == 0 or y == 1-x:
未能拒绝[0, 1.2] 和[1.2, 0.]。
靠近边缘的点和浮点运算错误
即使对于像[0.2,0.8] 这样的点正好在边缘上，浮点错误也会导致点错误分类。以下几点将使叉积不完全等于零。详情见下文：
p.contains([0.2,0.8], debug=True) # True
# k = 0, d1 = [0 1], d2 = [0.2 0.8], z = -0.20000000000000001110223024625157, s = -1
# k = 1, d1 = [ 1 -1], d2 = [ 0.2 -0.2], z = 0.00000000000000005551115123125783, s = 0
# k = 2, d1 = [-1 0], d2 = [-0.8 0.8], z = -0.80000000000000004440892098500626, s = -1
这就是为什么我们必须添加一个半径为atol 的球来检查给定公差是否实际上为零：
if np.allclose(z, 0, atol=atol):
s_ = 0
# ...
else:
s_ = np.sign(z)
这意味着我们必须接受距离边缘足够近的点（从两侧）被认为包含在多边形中。这是 Float Arithmetic 所固有的，您能做的最好的事情就是将 atol 调整为您的应用程序可接受的值。或者你可以找出另一个不存在这个问题的逻辑或数据模型。
缩放到精确比例
如果我们放大到精确比例以查看靠近边缘的情况，我们得到：
n = 40
lim = np.array([0.5,0.5]) + [-n*e/2, +n*e/2]
x = np.linspace(lim[0], lim[1], 30)
X, Y = np.meshgrid(x, x)
x, y = X.reshape(-1), Y.reshape(-1)
_, axe = p.plot()
axe.set_xlim(lim)
axe.set_ylim(lim)
for r in zip(x,y):
q = p.contains(r)
axe.plot(*r, 'o', color=cols[q], markersize=2)
我们看到一些非常靠近边缘但在多边形内部或外部的点被归类为“在边缘上”。这是因为 epsilon 球准则。您还可以观察到点的间距不相等（无论我是否使用linspace），因为不可能将10 表示为2 的整数幂。
结论
上面的解决方案是您的问题的概括，在O(n) 中执行。它可能看起来有点矫枉过正，但它是通用（适用于任何正多边形）和全面（它依赖于众所周知的几何概念）。
实际上，该算法只是检查点在通过路径时是否位于多边形所有边的同一侧。如果是，那么它就断定该点在多边形内部，即！
上述解决方案当然会受到浮点算术错误的影响，因为它依赖于浮点计算（见点10）。幸运的是，使用 epsilon 球测试，我们可以缓解它。
参考
如果你想更深入地理解有限精度算术，我建议你阅读这本优秀的书：Accuracy and Stability of Numerical Algorithms, J. Higham。
奖金
下面是所有答案与试验数据集的比较：
我们可以就软检查强调的不同类型的“错误”提供一些背景信息：
- 由于设计缺陷，
11 和 12 点被 @MagedSaeed 错误分类；
-
@MahmoudElshahat 在其逻辑中使用直接浮点相等，这就是为什么它为等价类的点返回不同的结果（例如：点 0 到 2，多边形点）；
-
@SamMason 使用最简单的逻辑和 epsilon 球测试。它似乎有最大的错误率，因为它不区分内部和边界上（我们可以忽略确切答案为0 并且其函数返回1 的所有点）。 IMO，这是在O(1) 中运行的最佳if-then 算法。此外，它可以轻松更新以考虑 3 种状态逻辑；
- 点
10 是故意设计来挑战边界处的逻辑，从距离上看，该点明显位于多边形之外，其数量级与机器精度相当。这是浮点算术误差变得显着并且逻辑更加模糊的地方：离边缘多远是可以接受和易于处理的？