Find the intersection of a line with a plane

From Rosetta Code
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].

C[edit]

Straightforward application of the intersection formula, prints usage on incorrect invocation.

 
/*Abhishek Ghosh, 12th October 2017*/
 
#include<stdio.h>
 
typedef struct{
double x,y,z;
}vector;
 
vector addVectors(vector a,vector b){
return (vector){a.x+b.x,a.y+b.y,a.z+b.z};
}
 
vector subVectors(vector a,vector b){
return (vector){a.x-b.x,a.y-b.y,a.z-b.z};
}
 
double dotProduct(vector a,vector b){
return a.x*b.x + a.y*b.y + a.z*b.z;
}
 
vector scaleVector(double l,vector a){
return (vector){l*a.x,l*a.y,l*a.z};
}
 
vector intersectionPoint(vector lineVector, vector linePoint, vector planeNormal, vector planePoint){
vector diff = subVectors(linePoint,planePoint);
 
return addVectors(addVectors(diff,planePoint),scaleVector(-dotProduct(diff,planeNormal)/dotProduct(lineVector,planeNormal),lineVector));
}
 
int main(int argC,char* argV[])
{
vector lV,lP,pN,pP,iP;
 
if(argC!=5)
printf("Usage : %s <line direction, point on line, normal to plane and point on plane given as (x,y,z) tuples separated by space>");
else{
sscanf(argV[1],"(%lf,%lf,%lf)",&lV.x,&lV.y,&lV.z);
sscanf(argV[3],"(%lf,%lf,%lf)",&pN.x,&pN.y,&pN.z);
 
if(dotProduct(lV,pN)==0)
printf("Line and Plane do not intersect, either parallel or line is on the plane");
else{
sscanf(argV[2],"(%lf,%lf,%lf)",&lP.x,&lP.y,&lP.z);
sscanf(argV[4],"(%lf,%lf,%lf)",&pP.x,&pP.y,&pP.z);
 
iP = intersectionPoint(lV,lP,pN,pP);
 
printf("Intersection point is (%lf,%lf,%lf)",iP.x,iP.y,iP.z);
}
}
 
return 0;
}
 

Invocation and output:

C:\rosettaCode>linePlane.exe (0,-1,-1) (0,0,10) (0,0,1) (0,0,5)
Intersection point is (0.000000,-5.000000,5.000000)

Perl 6[edit]

Works with: Rakudo version 2016.11
Translation of: Python
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 orthogonal; 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 intersects 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) )
);
Output:
Intersection at point: (0 -5 5)

Python[edit]

Based on the approach at geomalgorithms.com[1]

#!/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)
Output:
intersection at [ 0 -5  5]

Racket[edit]

Translation of: Sidef
#lang racket
;; {{trans|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)))
Output:

No output -- all tests passed!

REXX[edit]

version 1[edit]

This program does NOT handle the case when the line is parallel to or within the plane.

/* 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')'
Output:
0*x + 0*y + 1*z = 5
Intersection: P(0,-5,5)

version 2[edit]

handle the case that the line is parallel to the plane or lies within it.

/*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 */
Select
When a=0 Then
pd=''
When a=1 Then
pd='x'
When a=-1 Then
pd='-x'
Otherwise
pd=a'*x'
End
pd=pd
yy=mk2('y',b)
Select
When left(yy,1)='-' Then
pd=pd '-' substr(yy,2)
When left(yy,1)='0' Then
Nop
Otherwise
pd=pd '+' yy
End
zz=mk2('z',c)
Select
When left(zz,1)='-' Then
pd=pd '-' substr(zz,2)
When left(zz,1)='0' Then
Nop
Otherwise
pd=pd '+' zz
End
pd=pd '=' d
 
Say 'Plane definition:' pd
 
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=mk('x',a.1,v.1) ';' mk('y',a.2,v.2) ';' mk('z',a.3,v.3)
Say 'Line definition:' ld
 
Say 'Intersection: P('||x','y','z')'
Exit
 
Mk: Procedure
/*---------------------------------------------------------------------
* build part of line definition
*--------------------------------------------------------------------*/

Parse Arg v,aa,vv
If aa<>0 Then
res=v'='aa
Else
res=v'='
Select
When vv=0 Then
res=res||'0'
When vv=-1 Then
res=res||'-t'
When vv<0 Then
res=res||vv'*t'
Otherwise Do
If res=v'=' Then Do
If vv=1 Then
res=res||'t'
Else
res=res||vv'*t'
End
Else Do
If vv=1 Then
res=res||'+t'
Else
res=res||'+'vv'*t'
End
End
End
Return res
 
mk2: Procedure
/*---------------------------------------------------------------------
* build part of plane definition
*--------------------------------------------------------------------*/

Parse Arg v,u
Select
When u=0 Then
res=''
When u=1 Then
res=v
When u=-1 Then
res='-'v
When u<0 Then
res=u'*'v
Otherwise Do
If pd<>'' Then
res='+'u'*'v
Else
res=u'*'v
End
End
Return res
Output:
Plane definition: x+2*y+3*z=18
Line definition: x=3*t ; y=2+2*t ; z=4+t
Intersection: P(0.6,2.4,4.2)

Sidef[edit]

Translation of: Perl 6
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]),
))
Output:
Intersection at point: [0, -5, 5]

zkl[edit]

Translation of: Perl 6
Translation of: Python
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);
}
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) ))
);
Output:
Intersection at point: L(0,-5,5)

References[edit]

  1. http://geomalgorithms.com/a05-_intersect-1.html