Find the intersection of a line with a plane

Find the intersection of a line with a plane is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Finding the intersection of an infinite ray with a plane in 3D is an important topic in collision detection.

Task

Find the point of intersection for the infinite ray with direction (0,-1,-1) passing through position (0, 0, 10) with the infinite plane with a normal vector of (0, 0, 1) and which passes through [0, 0, 5].

Perl 6

Works with: Rakudo version 2016.11

Translation of: Python

<lang perl6>class Line {

   has $.P0; # point
   has $.u⃗;  # ray

} class Plane {

   has $.V0; # point
   has $.n⃗;  # normal

}

sub infix:<∙> ( @a, @b where +@a == +@b ) { [+] @a «*» @b } # dot product

sub line-plane-intersection ($𝑳, $𝑷) {

   my $cos = $𝑷.n⃗ ∙ $𝑳.u⃗; # cosine between normal & ray
   return 'Vectors are orthoganol; no intersection or line within plane'
     if $cos == 0;
   my $𝑊 = $𝑳.P0 «-» $𝑷.V0;      # difference between P0 and V0
   my $S𝐼 = -($𝑷.n⃗ ∙ $𝑊) / $cos;  # line segment where it intersets the plane
   $𝑊 «+» $S𝐼 «*» $𝑳.u⃗ «+» $𝑷.V0; # point where line intersects the plane
}

say 'Intersection at point: ', line-plane-intersection(

    Line.new( :P0(0,0,10), :u⃗(0,-1,-1) ),
   Plane.new( :V0(0,0, 5), :n⃗(0, 0, 1) )
 );</lang>

Output:

Intersection at point: (0 -5 5)

Python

Based on the approach at geomalgorithms.com

<lang python>#!/bin/python from __future__ import print_function import numpy as np

def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):

ndotu = planeNormal.dot(rayDirection) if abs(ndotu) < epsilon: raise RuntimeError("no intersection or line is within plane")

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

if __name__=="__main__": #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

Psi = LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint) print ("intersection at", Psi)</lang>

Output:

intersection at [ 0 -5  5]

Racket

Translation of: Sidef

<lang racket>#lang racket

Translation of: Sidef

vectors are represented by lists

(struct Line (P0 u⃗))

(struct Plane (V0 n⃗))

(define (· a b) (apply + (map * a b)))

(define (line-plane-intersection L P)

 (match-define (cons (Line P0 u⃗) (Plane V0 n⃗)) (cons L P))  
 (define cos (· n⃗ u⃗))
 (when (zero? cos) (error "vectors are orthoganal"))
 (define W (map - P0 V0))
 (define *SI (let ((SI (- (/ (· n⃗ W) cos)))) (λ (n) (* SI n))))
 (map + W (map *SI u⃗) V0))

(module+ test

 (require rackunit)
 (check-equal?
  (line-plane-intersection (Line '(0 0 10) '(0 -1 -1))
                           (Plane '(0 0 5) '(0 0 1)))
  '(0 -5 5)))</lang>

Output:

No output -- all tests passed!

REXX

version 1

This program does NOT handle the case when the line is parallel to or within the plane. <lang rexx>/* REXX */ Parse Value '0 0 1' With n.1 n.2 n.3 /* Normal Vector of the plane */ Parse Value '0 0 5' With p.1 p.2 p.3 /* Point in the plane */ Parse Value '0 0 10' With a.1 a.2 a.3 /* Point of the line */ Parse Value '0 -1 -1' With v.1 v.2 v.3 /* Vector of the line */

a=n.1 b=n.2 c=n.3 d=n.1*p.1+n.2*p.2+n.3*p.3 /* Parameter form of the plane */ Say a'*x +' b'*y +' c'*z =' d

t=(d-(a*a.1+b*a.2+c*a.3))/(a*v.1+b*v.2+c*v.3)

x=a.1+t*v.1 y=a.2+t*v.2 z=a.3+t*v.3

Say 'Intersection: P('||x','y','z')'</lang>

Output:

0*x + 0*y + 1*z = 5
Intersection: P(0,-5,5)

version 2

handle the case that the line is parallel to the plane or lies withi it.
Beautification of plane and line definitions is yet tbd. <lang rexx>/*REXX*/ Parse Value '1 2 3' With n.1 n.2 n.3 Parse Value '3 3 3' With p.1 p.2 p.3 Parse Value '0 2 4' With a.1 a.2 a.3 Parse Value '3 2 1' With v.1 v.2 v.3

a=n.1 b=n.2 c=n.3 d=n.1*p.1+n.2*p.2+n.3*p.3 /* Parameter form of the plane */ pd= If a<>0 Then pd=a'*x' If b<>0 Then pd=pd '+' b'*y' If c<>0 Then pd=pd '+' c'*z' Say 'Planes definition:' pd '=' d

ip=0 Do i=1 To 3

 ip=ip+n.i*v.i
 dd=n.1*a.1+n.2*a.2+n.3*a.3
 End

If ip=0 Then Do

 If dd=d Then
   Say 'Line is part of the plane'
 Else
   Say 'Line is parallel to the plane'
 Exit
 End

t=(d-(a*a.1+b*a.2+c*a.3))/(a*v.1+b*v.2+c*v.3)

x=a.1+t*v.1 y=a.2+t*v.2 z=a.3+t*v.3

ld='x='a.1; If v.1<>0 Then ld=ld'+t*'v.1 ld=ld '; y='a.2; If v.2<>0 Then ld=ld'+t*'v.2 ld=ld '; z='a.3; If v.3<>0 Then ld=ld'+t*'v.3 Say 'Line definition:' ld

Say 'Intersection: P('||x','y','z')'</lang>

Output:

Plane's definition: 1*x + 2*y + 3*z = 18
Line definition: x=0+t*3 ; y=2+t*2 ; z=4+t*1
Intersection: P(0.6,2.4,4.2)

Sidef

Translation of: Perl 6

<lang ruby>struct Line {

   P0,       # point
   u⃗,        # ray

}

struct Plane {

   V0,       # point
   n⃗,        # normal

}

func dot_prod(a, b) { a »*« b -> sum }

func line_plane_intersection(𝑳, 𝑷) {

   var cos = dot_prod(𝑷.n⃗, 𝑳.u⃗) ->
    || return 'Vectors are orthogonal'
   var 𝑊 = (𝑳.P0 »-« 𝑷.V0)
   var S𝐼 = -(dot_prod(𝑷.n⃗, 𝑊) / cos)
   𝑊 »+« (𝑳.u⃗ »*» S𝐼) »+« 𝑷.V0

}

say ('Intersection at point: ', line_plane_intersection(

        Line(P0: [0,0,10], u⃗: [0,-1,-1]),
       Plane(V0: [0,0, 5], n⃗: [0, 0, 1]),

))</lang>

Output:

Intersection at point: [0, -5, 5]

zkl

Translation of: Perl 6

Translation of: Python

<lang zkl>class Line { fcn init(pxyz, ray_xyz) { var pt=pxyz, ray=ray_xyz; } } class Plane{ fcn init(pxyz, normal_xyz){ var pt=pxyz, normal=normal_xyz; } }

fcn dotP(a,b){ a.zipWith('*,b).sum(0.0); } # dot product --> x fcn linePlaneIntersection(line,plane){

  cos:=dotP(plane.normal,line.ray); # cosine between normal & ray
  _assert_((not cos.closeTo(0,1e-6)),
     "Vectors are orthoganol; no intersection or line within plane");
  w:=line.pt.zipWith('-,plane.pt); # difference between P0 and V0
  si:=-dotP(plane.normal,w)/cos;   # line segment where it intersets the plane
     # point where line intersects the plane:
  //w.zipWith('+,line.ray.apply('*,si)).zipWith('+,plane.pt);  // or
  w.zipWith('wrap(w,r,pt){ w + r*si + pt },line.ray,plane.pt);

}</lang> <lang zkl>println("Intersection at point: ", linePlaneIntersection(

  Line( T(0.0, 0.0, 10.0), T(0.0, -1.0, -1.0) ),
  Plane(T(0.0, 0.0,  5.0), T(0.0,  0.0,  1.0) ))

);</lang>

Output:

Intersection at point: L(0,-5,5)