Przecięcie linii nD z wypukłym kadłubem w Pythonie

Question

Stworzyłem wypukły kadłub za pomocą scipy.spatial.ConvexHull. Muszę obliczyć punkt przecięcia między wypukłym kadłubem i promieniem, zaczynając od 0 i w kierunku jakiegoś innego określonego punktu. Wiadomo, że wypukły kadłub zawiera 0, więc skrzyżowanie powinno być zagwarantowane. Rozmiar problemu może wynosić od 2 do 5. Próbowałem wyszukiwać google, ale nie znalazłem odpowiedzi. Mam nadzieję, że jest to powszechny problem ze znanymi rozwiązaniami w geometrii obliczeniowej. Dziękuję Ci.

Answer 1

Według qhull.org, punkty x powierzchni kadłuba wypukłego zweryfikować V.x+b=0, gdzie V i b są podane przez hull.equations. Jeśli U jest wektorem promienia rozpoczynającego się w O, równanie on ray wynosi x=αU, α>0. więc przecięcie promienia i aspektu to x = αU = -b/(V.U) U. Unikalny punkt przecięcia z kadłubem odpowiadają minimum wartości dodatnie α:

następny kod dać:

from pylab import * 
from scipy.spatial import ConvexHull 

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

Jest to czysty roztwór numpy więc szybko. Exemple do 1 miliona punktów w [-1,1]^3 kostki:

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

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

Answer 2

Jak wspomniał Ante w komentarzach, musisz znaleźć najbliższy punkt przecięcia wszystkich linii/płaszczyzn/hiper-płaszczyzn w kadłubie.

Aby znaleźć punkt przecięcia promienia z hiperpłaszczyzną, należy wykonać rzut punktowy znormalizowanego promienia z normalną hiperpłaszczyzną, która powie, jak daleko w kierunku hiperpłaszczyzny normalnej poruszasz się dla każdej odległości jednostki wzdłuż promienia .

Jeśli produkt punktowy jest ujemny, oznacza to, że hiperpłaszczyzna znajduje się w przeciwnym kierunku promienia, jeśli zero, oznacza to, że promień jest równoległy do ​​niego i nie przecina się.

Po uzyskaniu dodatniego produktu punktowego można ustalić, w jakiej odległości znajduje się hiperpłaszczyzna w kierunku promienia, dzieląc odległość płaszczyzny w kierunku płaszczyzny normalnej przez produkt punktowy. Na przykład, jeśli płaszczyzna znajduje się w odległości 3 jednostek, a iloczyn jest równy 0,5, to otrzymasz tylko 0,5 jednostki bliżej każdej jednostki poruszającej się wzdłuż promienia, więc hiperpłaszczyzna jest 3/0,5 = 6 jednostek dalej w kierunku promienia .

Po obliczeniu odległości dla wszystkich hiperpłaszczyzn i znalezieniu najbliższego, punkt przecięcia jest po prostu promieniem pomnożonym przez najbliższą odległość.

Oto rozwiązanie Pythona (normalizacja funkcji wynosi here)

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) 

kod testowy w 2D (rozwiązanie uogólnia większych wymiarów)

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 

Wydajność:

[ 0.66666667 -0.66666667]