Skrzyżowanie linii 3D z płaszczyzną

Oto przykład Pythona, który znajduje przecięcie linii i płaszczyzny.

Gdzie płaszczyzna może być punktem i normalnym lub wektorem 4d (postać normalna), W poniższych przykładach (podano kod dla obu).

Zauważ również, że ta funkcja oblicza wartość reprezentującą, gdzie punkt znajduje się w linii (o nazwie fac w kodzie poniżej). Możesz też chcieć to zwrócić, ponieważ wartości od 0 do 1 przecinają segment linii - co może być przydatne dla dzwoniący.

Inne szczegóły odnotowane w kodzie-komentarze.

Uwaga: Ten przykład używa czystych funkcji, bez żadnych zależności - aby ułatwić przenoszenie się do innych języków. Z typem danych Vector i przeciążeniem operatora, może być bardziej zwięzły(w przykładzie poniżej).

# intersection function
def isect_line_plane_v3(p0, p1, p_co, p_no, epsilon=1e-6):
    """
    p0, p1: define the line
    p_co, p_no: define the plane:
        p_co is a point on the plane (plane coordinate).
        p_no is a normal vector defining the plane direction;
             (does not need to be normalized).

    return a Vector or None (when the intersection can't be found).
    """

    u = sub_v3v3(p1, p0)
    dot = dot_v3v3(p_no, u)

    if abs(dot) > epsilon:
        # the factor of the point between p0 -> p1 (0 - 1)
        # if 'fac' is between (0 - 1) the point intersects with the segment.
        # otherwise:
        #  < 0.0: behind p0.
        #  > 1.0: infront of p1.
        w = sub_v3v3(p0, p_co)
        fac = -dot_v3v3(p_no, w) / dot
        u = mul_v3_fl(u, fac)
        return add_v3v3(p0, u)
    else:
        # The segment is parallel to plane
        return None

# ----------------------
# generic math functions

def add_v3v3(v0, v1):
    return (
        v0[0] + v1[0],
        v0[1] + v1[1],
        v0[2] + v1[2],
        )


def sub_v3v3(v0, v1):
    return (
        v0[0] - v1[0],
        v0[1] - v1[1],
        v0[2] - v1[2],
        )


def dot_v3v3(v0, v1):
    return (
        (v0[0] * v1[0]) +
        (v0[1] * v1[1]) +
        (v0[2] * v1[2])
        )


def len_squared_v3(v0):
    return dot_v3v3(v0, v0)


def mul_v3_fl(v0, f):
    return (
        v0[0] * f,
        v0[1] * f,
        v0[2] * f,
        )

Jeśli płaszczyzna jest zdefiniowana jako wektor 4d (postać normalna), musimy znaleźć punkt na płaszczyźnie, a następnie obliczyć przecięcie jak wcześniej (patrz p_co zadanie).

def isect_line_plane_v3_4d(p0, p1, plane, epsilon=1e-6):    
    u = sub_v3v3(p1, p0)
    dot = dot_v3v3(plane, u)

    if abs(dot) > epsilon:
        # calculate a point on the plane
        # (divide can be omitted for unit hessian-normal form).
        p_co = mul_v3_fl(plane, -plane[3] / len_squared_v3(plane))

        w = sub_v3v3(p0, p_co)
        fac = -dot_v3v3(plane, w) / dot
        u = mul_v3_fl(u, fac)
        return add_v3v3(p0, u)
    else:
        return None

Dla dalszych odniesień, zostało to zaczerpnięte z blendera i zaadaptowane do Pythona. isect_line_plane_v3() w math_geom.c

Dla jasności, oto wersje używające mathutils API (które można modyfikować dla innych bibliotek matematycznych z przeciążeniem operatora).

# point-normal plane
def isect_line_plane_v3(p0, p1, p_co, p_no, epsilon=1e-6):
    u = p1 - p0
    dot = p_no * u
    if abs(dot) > epsilon:
        w = p0 - p_co
        fac = -(plane * w) / dot
        return p0 + (u * fac)
    else:
        return None


# normal-form plane
def isect_line_plane_v3_4d(p0, p1, plane, epsilon=1e-6):    
    u = p1 - p0
    dot = plane.xyz * u
    if abs(dot) > epsilon:
        p_co = plane.xyz * (-plane[3] / plane.xyz.length_squared)

        w = p0 - p_co
        fac = -(plane * w) / dot
        return p0 + (u * fac)
    else:
        return None

Używanie numpy i Pythona:

#Based on http://geomalgorithms.com/a05-_intersect-1.html
from __future__ import print_function
import numpy as np

epsilon=1e-6

#Define plane
planeNormal = np.array([0, 0, 1])
planePoint = np.array([0, 0, 5]) #Any point on the plane

#Define ray
rayDirection = np.array([0, -1, -1])
rayPoint = np.array([0, 0, 10]) #Any point along the ray

ndotu = planeNormal.dot(rayDirection) 

if abs(ndotu) < epsilon:
    print ("no intersection or line is within plane")

w = rayPoint - planePoint
si = -planeNormal.dot(w) / ndotu
Psi = w + si * rayDirection + planePoint

print ("intersection at", Psi)

Na podstawie tego kodu Matlab (minus sprawdzanie przecięcia), w Pythonie

# n: normal vector of the Plane 
# V0: any point that belongs to the Plane 
# P0: end point 1 of the segment P0P1
# P1:  end point 2 of the segment P0P1
n = np.array([1., 1., 1.])
V0 = np.array([1., 1., -5.])
P0 = np.array([-5., 1., -1.])
P1 = np.array([1., 2., 3.])

w = P0 - V0;
u = P1-P0;
N = -np.dot(n,w);
D = np.dot(n,u)
sI = N / D
I = P0+ sI*u
print I

Wynik

[-3.90909091  1.18181818 -0.27272727]

Sprawdziłem to graficznie wygląda na to, że działa,

Uważam, że jest to bardziej solidna implementacja łącza współdzielonego przed

To pytanie jest stare, ale ponieważ jest takie o wiele wygodniejsze rozwiązanie, pomyślałem, że może komuś pomóc.

Przecięcia płaszczyzn i linii są dość eleganckie, gdy są wyrażone w współrzędnych jednorodnych, ale załóżmy, że chcemy tylko rozwiązania:

Istnieje wektor 4x1 p, który opisuje płaszczyznę taką, że p^t * x =0 dla dowolnego punktu jednorodnego na płaszczyźnie. Następnie Oblicz współrzędne pluckera dla linii L = ab^T-ba^t gdzie a = {point_1; 1}, b = {point_2; 1}, oba 4x1 na linia

Oblicz: x = L*p = {x0,x1,x2,x3}

X_intersect=({x0,x1,x2}/x3) gdzie jeśli x3 jest zerem, nie ma przecięcia w sensie euklidesowym.

Jeśli masz dwa punkty p i q, które definiują linię I płaszczyznę w ogólnej postaci kartezjańskiej ax + przez + cz + d = 0, możesz użyć metody parametrycznej.

Jeśli potrzebujesz tego do celów kodowania, oto fragment kodu javascript:

/**
* findLinePlaneIntersectionCoords (to avoid requiring unnecessary instantiation)
* Given points p with px py pz and q that define a line, and the plane
* of formula ax+by+cz+d = 0, returns the intersection point or null if none.
*/
function findLinePlaneIntersectionCoords(px, py, pz, qx, qy, qz, a, b, c, d) {
    var tDenom = a*(qx-px) + b*(qy-py) + c*(qz-pz);
    if (tDenom == 0) return null;

    var t = - ( a*px + b*py + c*pz + d ) / tDenom;

    return {
        x: (px+t*(qx-px)),
        y: (py+t*(qy-py)),
        z: (pz+t*(qz-pz))
    };
}

// Example (plane at y = 10  and perpendicular line from the origin)
console.log(JSON.stringify(findLinePlaneIntersectionCoords(0,0,0,0,1,0,0,1,0,-10)));

// Example (no intersection, plane and line are parallel)
console.log(JSON.stringify(findLinePlaneIntersectionCoords(0,0,0,0,0,1,0,1,0,-10)));