Поиск точки в 3D-поли в питоне

Я проверил версию QHull (сверху) и решение для линейного программирования (например, см. этот вопрос). Пока что использование QHull кажется лучшим выбором, хотя я могу упустить некоторые оптимизации с scipy.spatial LP.

import numpy
import numpy.random
from numpy import zeros, ones, arange, asarray, concatenate
from scipy.optimize import linprog

from scipy.spatial import ConvexHull

def pnt_in_cvex_hull_1(hull, pnt):
    '''
    Checks if `pnt` is inside the convex hull.
    `hull` -- a QHull ConvexHull object
    `pnt` -- point array of shape (3,)
    '''
    new_hull = ConvexHull(concatenate((hull.points, [pnt])))
    if numpy.array_equal(new_hull.vertices, hull.vertices): 
        return True
    return False


def pnt_in_cvex_hull_2(hull_points, pnt):
    '''
    Given a set of points that defines a convex hull, uses simplex LP to determine
    whether point lies within hull.
    `hull_points` -- (N, 3) array of points defining the hull
    `pnt` -- point array of shape (3,)
    '''
    N = hull_points.shape[0]
    c = ones(N)
    A_eq = concatenate((hull_points, ones((N,1))), 1).T   # rows are x, y, z, 1
    b_eq = concatenate((pnt, (1,)))
    result = linprog(c, A_eq=A_eq, b_eq=b_eq)
    if result.success and c.dot(result.x) == 1.:
        return True
    return False


points = numpy.random.rand(8, 3)
hull = ConvexHull(points, incremental=True)
hull_points = hull.points[hull.vertices, :]
new_points = 1. * numpy.random.rand(1000, 3)

где

%%time
in_hull_1 = asarray([pnt_in_cvex_hull_1(hull, pnt) for pnt in new_points], dtype=bool)

производит:

CPU times: user 268 ms, sys: 4 ms, total: 272 ms
Wall time: 268 ms

и

%%time
in_hull_2 = asarray([pnt_in_cvex_hull_2(hull_points, pnt) for pnt in new_points], dtype=bool)

производит

CPU times: user 3.83 s, sys: 16 ms, total: 3.85 s
Wall time: 3.85 s

Аналогичный вопрос был задан здесь, но с упором на эффективность.

Подход scipy.spatial.ConvexHull, предложенный здесь @Brian и @fatalaccidents работает, но работает очень медленно, если вам нужно проверить более одной точки.

Что ж, наиболее эффективное решение также исходит из scipy.spatial, но использует тесселяцию Delaunay:

from scipy.spatial import Delaunay

Delaunay(poly).find_simplex(point) >= 0  # True if point lies within poly

Это работает, потому что -1 возвращает .find_simplex(point), если точка не находится ни в одном из симплексов (т.е. вне триангуляции). (Примечание: он работает в N измерениях, а не только в 2/3D.)

Сравнение производительности

Сначала за один балл:

import numpy
from scipy.spatial import ConvexHull, Delaunay

def in_poly_hull_single(poly, point):
    hull = ConvexHull(poly)
    new_hull = ConvexHull(np.concatenate((poly, [point])))
    return np.array_equal(new_hull.vertices, hull.vertices)

poly = np.random.rand(65, 3)
point = np.random.rand(3)

%timeit in_poly_hull_single(poly, point)
%timeit Delaunay(poly).find_simplex(point) >= 0

Результат:

2.63 ms ± 280 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1.49 ms ± 153 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Таким образом, подход Delaunay быстрее. Но это зависит от размера полигона! Я обнаружил, что для многоугольника, состоящего более чем из ~65 точек, подход Delaunay становится все более медленным, в то время как подход ConvexHull остается почти постоянным по скорости.

Для нескольких точек:

def in_poly_hull_multi(poly, points):
    hull = ConvexHull(poly)
    res = []
    for p in points:
        new_hull = ConvexHull(np.concatenate((poly, [p])))
        res.append(np.array_equal(new_hull.vertices, hull.vertices))
    return res

points = np.random.rand(10000, 3)

%timeit in_poly_hull_multi(poly, points)
%timeit Delaunay(poly).find_simplex(points) >= 0

Результат:

155 ms ± 9.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1.81 ms ± 106 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Таким образом, Delaunay дает экстремальное увеличение скорости; не говоря уже о том, как долго нам придется ждать 10 000 баллов или более. В таком случае размер полигона больше не имеет большого значения.

Таким образом, Delaunay не только намного быстрее, но и очень лаконичен в коде.

Спасибо всем, кто прокомментировал. Для тех, кто ищет ответ на этот вопрос, я нашел тот, который работает для некоторых случаев (но не для сложных случаев).

Я использую scipy.spatial.ConvexHull, как предложил шонгололо, но с небольшим поворотом. Я создаю трехмерную выпуклую оболочку облака точек, а затем добавляю точку, которую я проверяю, в «новое» облако точек и создаю новую трехмерную выпуклую оболочку. Если они идентичны, то я предполагаю, что он должен быть внутри выпуклой оболочки. Я все равно был бы признателен, если бы у кого-то был более надежный способ сделать это, так как я считаю это немного хакерским. Код будет выглядеть примерно так:

from scipy.spatial import ConvexHull

def pnt_in_pointcloud(points, new_pt):
    hull = ConvexHull(points)
    new_pts = points + new_pt
    new_hull = ConvexHull(new_pts)
    if hull == new_hull: 
        return True
    else:
        return False

Надеюсь, это поможет кому-то в будущем искать ответ! Спасибо!