Пересечение линии nD с выпуклой оболочкой в ​​Python

Question

Я создал выпуклую оболочку, используя scipy.spatial.ConvexHull. Мне нужно вычислить точку пересечения между выпуклой оболочкой и лучом, начиная с 0 и в направлении некоторой другой определенной точки. Известно, что выпуклая оболочка содержит 0, поэтому пересечение должно быть гарантировано. Размер проблемы может варьироваться от 2 до 5. Я пробовал поиск в Google, но не нашел ответа. Я надеюсь, что это общая проблема с известными решениями в вычислительной геометрии. Спасибо.

Answer 1

Согласно qhull.org, пункты x грани выпуклой оболочки проверить V.x+b=0, где V а также b даны hull.equations, (. обозначает точечное произведение здесь. V нормальный вектор длины один.)

Если V - нормаль, b - смещение, а x - точка внутри выпуклой оболочки, то Vx+b <0.

Если U является вектором луча, начинающегося в O, уравнение луча x=αU, α>0, поэтому пересечение луча является x = αU = -b/(V.U) U, Уникальная точка пересечения с корпусом соответствует минимуму положительных значений α:

Следующий код дать это:

import numpy as np
from scipy.spatial import ConvexHull

def hit(U,hull):
    eq=hull.equations.T
    V,b=eq[:-1],eq[-1]
    alpha=-b/np.dot(V,U)
    return np.min(alpha[alpha>0])*U

Это чистое решение, так что это быстро. Пример для 1 миллиона очков в [-1,1]^3 куб:

In [13]: points=2*np.random.rand(1e6,3)-1;hull=ConvexHull(points)

In [14]: %timeit x=hit(np.ones(3),hull) 
#array([ 0.98388702,  0.98388702,  0.98388702])
10000 loops, best of 3: 30 µs per loop

Answer 2

Как упомянуто Анте в комментариях, вам нужно найти ближайшее пересечение всех линий / плоскостей / гиперплоскостей в корпусе.

Чтобы найти пересечение луча с гиперплоскостью, произведите точечное произведение нормализованного луча на нормаль гиперплоскости, которая скажет вам, как далеко в направлении нормали гиперплоскости вы продвигаетесь на каждое единичное расстояние вдоль луча.

Если скалярное произведение отрицательно, это означает, что гиперплоскость находится в направлении, противоположном лучу, если ноль, то это означает, что луч параллелен ему и не будет пересекаться.

Получив положительное точечное произведение, вы можете определить, насколько далеко гиперплоскость находится в направлении луча, поделив расстояние плоскости в направлении нормали плоскости к точечному произведению. Например, если плоскость находится на расстоянии 3 единицы, а точечное произведение равно 0,5, то вы приближаетесь только на 0,5 единицы к каждой единице, которую вы перемещаете вдоль луча, поэтому гиперплоскость находится на расстоянии 3 / 0,5 = 6 единиц в направлении луча.,

Как только вы вычислили это расстояние для всех гиперплоскостей и нашли ближайшую, точка пересечения - это просто луч, умноженный на ближайшее расстояние.

Вот решение в Python (функция нормализации отсюда):

def normalize(v):
    norm = np.linalg.norm(v)
    if norm == 0: 
       return v
    return v / norm

def find_hull_intersection(hull, ray_point):
    # normalise ray_point
    unit_ray = normalize(ray_point)
    # find the closest line/plane/hyperplane in the hull:
    closest_plane = None
    closest_plane_distance = 0
    for plane in hull.equations:
        normal = plane[:-1]
        distance = plane[-1]
        # if plane passes through the origin then return the origin
        if distance == 0:
            return np.multiply(ray_point, 0) # return n-dimensional zero vector 
        # if distance is negative then flip the sign of both the
        # normal and the distance:       
        if distance < 0:
            np.multiply(normal, -1);
            distance = distance * -1
        # find out how much we move along the plane normal for
        # every unit distance along the ray normal:
        dot_product = np.dot(normal, unit_ray)
        # check the dot product is positive, if not then the
        # plane is in the opposite direction to the rayL
        if dot_product > 0:  
            # calculate the distance of the plane
            # along the ray normal:          
            ray_distance = distance / dot_product
            # is this the closest so far:
            if closest_plane is None or ray_distance < closest_plane_distance:
                closest_plane = plane
                closest_plane_distance = ray_distance
    # was there no valid plane? (should never happen):
    if closest_plane is None:
        return None
    # return the point along the unit_ray of the closest plane,
    # which will be the intersection point
    return np.multiply(unit_ray, closest_plane_distance)

Тестовый код в 2D (решение обобщается на более высокие измерения):

from scipy.spatial import ConvexHull
import numpy as np

points = np.array([[-2, -2], [2, 0], [-1, 2]])
h = ConvexHull(points)
closest_point = find_hull_intersection(h, [1, -1])
print closest_point

выход:

[ 0.66666667 -0.66666667]