Предисловия
Ваш вопрос является специализацией " Является ли точка внутри правильного многоугольника? " обычными, я имею в виду не самопересекающимися или многоугольниками, как это часто можно найти в ГИС-системе.,И в продолжение, потому что вы спрашиваете о треугольнике, я буду считать, что многоугольник выпуклый.
Что интересно, так это то, что перекрестное произведение являетсяключ, чтобы решить это.Когда вы имеете дело с векторами в 2D плоскости, перекрестные произведения ортогональны этой плоскости.Полезная информация для извлечения: указывает ли она вверх или вниз?
Будут возникать арифметические ошибки с плавающей точкой, и они становятся критическими, когда перекрестный продукт получаетсяблизко к нулю, но не равно, тогда у него будет знак вместо нуля.
Чтобы проверить, находится ли ваша точка внутри многоугольника, он просто сводится к проверке, все ли перекрестные произведения между ребром иточки имеют одинаковые знаки, например: 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
Класс инициализируется со списком 2D точек (по умолчанию он ваш).
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, лучшее, что вы можете сделать, это настроить 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)
Мы видим, что некоторые точки, расположенные очень близко к краю, но внутри или снаружи многоугольника, классифицируются как "накрай».Это из-за критерия эпсилон-шара.Вы также можете заметить, что точки распределены неравномерно (независимо от того, использовал ли я linspace), поскольку невозможно выразить 10 как целую степень 2.
Заключение
Приведенное выше решение является обобщением вашей проблемы, выполненной в O(n).Это может показаться излишним, но оно общее (работает для любого правильного многоугольника) и всеобъемлющее (оно опирается на хорошо известную геометрическую концепцию).
На самом делеалгоритм просто проверяет, находится ли точка на одной стороне всех ребер многоугольника при прохождении по пути.Если это так, то он приходит к выводу, что точка находится внутри многоугольника, то есть!
Приведенное выше решение отклоняется от курса из-за арифметической ошибки с плавающей точкой, поскольку оно основано на вычислении с плавающей точкой (см. Точку 10).К счастью, используя тест с эпсилон-шариком, мы можем смягчить его.
Reference
Если вы хотите глубже понять арифметику конечной точности, я бы посоветовал вам прочитать превосходную книгу: Точность и устойчивость численных алгоритмов, Дж. Хайм .
Бонус
Ниже приведено сравнение всех ответов с пробным набором данных:
Мы можем дать некоторое представление о различного рода «ошибках», подчеркнутыхэта мягкая проверка:
- Точки
11 и 12 ошибочно классифицированы как @MagedSaeed из-за недостатка проекта; @MahmoudElshahat используют равенство с плавающей запятой в своей логике,вот почему он возвращает разные результаты для эквивалентного класса точек (например: точки 0 до 2, точки многоугольника); @SamMason использует простейшую логику с проверкой эпсилон-шара.Кажется, что он имеет наибольшую частоту ошибок, потому что он не различает внутри и на границе (мы можем игнорировать все точки, где точный ответ 0 и его функция вернула 1).ИМО, это лучший из представленных алгоритмов if-then, который работает в O(1).Кроме того, его можно легко обновить, чтобы учесть логику 3 состояний; - Точка
10 была специально разработана для вызова логики на границе, точка явно находится за пределами многоугольника от расстояния с величиной вокруг машины.точность.Вот где арифметическая ошибка с плавающей точкой становится существенной, а логика более нечеткой: насколько далеко от края допустимы и поддаются обработке?