Intersection Ligne-Plan 3D

Voici un exemple Python qui trouve l'intersection d'une ligne et d'un plan.

Où le plan peut être soit un point et une normale, soit un vecteur 4d (forme normale), Dans les exemples ci-dessous (le code pour les deux est fourni).

Notez également que cette fonction calcule une valeur représentant où le point est sur la ligne, (appelée fac dans le code ci-dessous). Vous voudrez peut-être renvoyer cela aussi, car les valeurs de 0 à 1 croisent le segment de ligne - ce qui peut être utile pour le appelant.

Autres détails notés dans le code-commentaires.

Note: Cet exemple utilise des fonctions pures, sans aucune dépendance - pour faciliter le passage à d'autres langues. Avec un type de données Vector et une surcharge d'opérateur, il peut être plus concis (inclus dans l'exemple ci-dessous).

# 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,
        )

Si le plan est défini comme un vecteur 4d (Forme Normale) , nous devons trouver un point sur le plan, puis calculer l'intersection comme précédemment (voir p_co affectation).

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

Pour plus de référence, cela a été pris de Blender et adapté à Python. isect_line_plane_v3() dans math_geom.c

Pour plus de clarté, voici des versions utilisant l'APImathutils (qui peut être modifiée pour d'autres bibliothèques mathématiques avec une surcharge d'opérateur) .

# 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

En utilisant numpy et python:

#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)

Basé sur CE code Matlab (moins les contrôles d'intersection), en Python

# 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

Résultat

[-3.90909091  1.18181818 -0.27272727]

Je l'ai vérifié graphiquement ça semble marcher,

Je crois que c'est une implémentation plus robuste du lien partagé avant

Cette question Est ancienne mais comme il existe une solution beaucoup plus pratique, j'ai pensé que cela pourrait aider quelqu'un.

Les intersections de plan et de ligne sont assez élégantes lorsqu'elles sont exprimées en coordonnées homogènes, mais supposons que vous voulez juste la solution:

Il existe un vecteur 4x1 p qui décrit le plan tel que p^T * x = 0 pour tout point homogène sur le plan. Ensuite, calculez les coordonnées de plucker pour la ligne L = ab^t-ba^t où a = {point_1; 1}, b = {point_2; 1}, les deux 4x1 sur la ligne

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

X_intersect=({x0, x1, x2}/x3) où si x3 est zéro, il n'y a pas d'intersection au sens euclidien.

Si vous avez deux points p et q qui définissent une ligne, et un plan sous la forme cartésienne générale ax + by + cz + d = 0, Vous pouvez utiliser la méthode paramétrique.

Si vous en avez besoin à des fins de codage, voici un extrait 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)));