Find the intersection of a line with a plane

From Rosetta Code
Revision as of 18:08, 25 December 2016 by rosettacode>Craigd (→‎{{header|zkl}}: rewrite a line for clarity)
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[1]

<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]

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)

References