Intersezione della linea nD con lo scafo convesso in Python

Question

Ho creato uno scafo convesso usando scipy.spatial.ConvexHull. Ho bisogno di calcolare il punto di intersezione tra lo scafo convesso e un raggio, iniziando da 0 e nella direzione di qualche altro punto definito. Lo scafo convesso è noto per contenere 0, quindi l'intersezione dovrebbe essere garantita. La dimensione del problema può variare tra 2 e 5. Ho provato alcune ricerche su google ma non ho trovato una risposta. Spero che questo sia un problema comune con soluzioni note nella geometria computazionale. Grazie.

Answer 1

Secondo qhull.org, i punti x di un lato dello scafo convesso verificare V.x+b=0, dove V e sono forniti da hull.equations. Se U è un vettore del raggio che inizia con O, l'equazione di raggio è x=αU, α>0. quindi l'intersezione di raggio di una faccetta è x = αU = -b/(V.U) U. Il punto di intersezione unico con lo scafo corrisponde al minimo dei valori positivi di α:

il codice successivo dare:

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 

È una soluzione NumPy puro così è veloce. Un esempio per 1 milione di punti nel [-1,1]^3 cubo:

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

Come menzionato da Ante nei commenti, è necessario trovare l'intersezione più vicina di tutte le linee/piani/iperpiani nello scafo.

Per trovare l'intersezione del raggio con l'iperpiano, fare un prodotto punto del raggio normalizzato con l'iperpiano normale, che indicherà quanto lontano nella direzione dell'ipercapiente normale si sposta per ogni distanza unità lungo il raggio .

Se il prodotto punto è negativo, significa che l'iperpiano si trova nella direzione opposta del raggio, se zero significa che il raggio è parallelo ad esso e non si intersecherà.

Una volta ottenuto un prodotto con punto positivo, è possibile calcolare quanto lontano l'iperpiano si trova nella direzione del raggio, dividendo la distanza del piano nella direzione del piano normale per il prodotto punto. Ad esempio se il piano è a 3 unità di distanza, e il prodotto punto è 0,5, allora ottieni solo 0,5 unità più vicine per ogni unità che muovi lungo il raggio, quindi l'iperpiano è 3/0,5 = 6 unità di distanza nella direzione del raggio .

Dopo aver calcolato questa distanza per tutti gli iperplan e trovato quello più vicino, il punto di intersezione è solo il raggio moltiplicato per la distanza più vicina.

Ecco una soluzione in Python (normalizzare funzione è da 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) 

Codice di prova in 2D (la soluzione generalizza di dimensioni superiori):

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 

uscita:

[ 0.66666667 -0.66666667]