Find the intersection of a line with a plane

From Rosetta Code
Task
Find the intersection of a line with a plane
You are encouraged to solve this task according to the task description, using any language you may know.

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

11l

F intersection_point(ray_direction, ray_point, plane_normal, plane_point)
   R ray_point - ray_direction * dot(ray_point - plane_point, plane_normal) / dot(ray_direction, plane_normal)

print(‘The ray intersects the plane at ’intersection_point((0.0, -1.0, -1.0), (0.0, 0.0, 10.0), (0.0, 0.0, 1.0), (0.0, 0.0, 5.0)))
Output:
The ray intersects the plane at (0, -5, 5)

Action!

INCLUDE "D2:REAL.ACT" ;from the Action! Tool Kit

DEFINE REALPTR="CARD"
TYPE VectorR=[REALPTR x,y,z]

PROC PrintVector(VectorR POINTER v)
  Print("(") PrintR(v.x)
  Print(",") PrintR(v.y)
  Print(",") PrintR(v.z)
  Print(")")
RETURN

PROC Vector(REAL POINTER vx,vy,vz VectorR POINTER v)
  v.x=vx v.y=vy v.z=vz
RETURN

PROC VectorSub(VectorR POINTER a,b,res)
  RealSub(a.x,b.x,res.x)
  RealSub(a.y,b.y,res.y)
  RealSub(a.z,b.z,res.z)
RETURN

PROC VectorDot(VectorR POINTER a,b REAL POINTER res)
  REAL tmp1,tmp2

  RealMult(a.x,b.x,res)
  RealMult(a.y,b.y,tmp1)
  RealAdd(res,tmp1,tmp2)
  RealMult(a.z,b.z,tmp1)
  RealAdd(tmp1,tmp2,res)
RETURN

PROC VectorMul(VectorR POINTER a REAL POINTER b VectorR POINTER res)
  RealMult(a.x,b,res.x)
  RealMult(a.y,b,res.y)
  RealMult(a.z,b,res.z)
RETURN

BYTE FUNC IsZero(REAL POINTER a)
  CHAR ARRAY s(10)

  StrR(a,s)
  IF s(0)=1 AND s(1)='0 THEN
    RETURN (1)
  FI
RETURN (0)

BYTE FUNC Intersection(VectorR POINTER
  rayVector,rayPoint,planeNormal,planePoint,result)
  
  REAL tmpx,tmpy,tmpz,prod1,prod2,prod3
  VectorR tmp

  Vector(tmpx,tmpy,tmpz,tmp)

  VectorSub(rayPoint,planePoint,tmp)
  VectorDot(tmp,planeNormal,prod1)
  VectorDot(rayVector,planeNormal,prod2)

  IF IsZero(prod2) THEN
    RETURN (1)
  FI

  RealDiv(prod1,prod2,prod3)
  VectorMul(rayVector,prod3,tmp)
  VectorSub(rayPoint,tmp,result)
RETURN (0)

PROC Test(VectorR POINTER rayVector,rayPoint,planeNormal,planePoint)
  BYTE res
  REAL px,py,pz
  VectorR p

  Vector(px,py,pz,p)
  res=Intersection(rayVector,rayPoint,planeNormal,planePoint,p)

  Print("Ray vector: ")
  PrintVector(rayVector) PutE()
  Print("Ray point: ")
  PrintVector(rayPoint) PutE()
  Print("Plane normal: ")
  PrintVector(planeNormal) PutE()
  Print("Plane point: ")
  PrintVector(planePoint) PutE()

  IF res=0 THEN
    Print("Intersection point: ")
    PrintVector(p) PutE()
  ELSEIF res=1 THEN
    PrintE("There is no intersection")
  FI
  PutE()
RETURN

PROC Main()
  REAL r0,r1,r5,r10,rm1    
  VectorR rayVector,rayPoint,planeNormal,planePoint

  Put(125) PutE() ;clear screen

  ValR("0",r0) ValR("1",r1) ValR("5",r5)
  ValR("10",r10) ValR("-1",rm1)

  Vector(r0,rm1,rm1,rayVector)
  Vector(r0,r0,r10,rayPoint)
  Vector(r0,r0,r1,planeNormal)
  Vector(r0,r0,r5,planePoint)
  Test(rayVector,rayPoint,planeNormal,planePoint)

  Vector(r1,r1,r0,rayVector)
  Vector(r1,r1,r0,rayPoint)
  Vector(r0,r0,r1,planeNormal)
  Vector(r5,r1,r0,planePoint)
  Test(rayVector,rayPoint,planeNormal,planePoint)
RETURN
Output:

Screenshot from Atari 8-bit computer

Ray vector: (0,-1,-1)
Ray point: (0,0,10)
Plane normal: (0,0,1)
Plane point: (0,0,5)
Intersection point: (0,-5,5)

Ray vector: (1,1,0)
Ray point: (1,1,0)
Plane normal: (0,0,1)
Plane point: (5,1,0)
There is no intersection

Ada

with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;

procedure Intersection is

   type Real is new Long_Float;

   package Real_Arrays is
      new Ada.Numerics.Generic_Real_Arrays (Real => Real);
   use Real_Arrays;

   package Real_IO is
      new Ada.Text_IO.Float_IO (Num => Real);

   subtype Vector_3D is Real_Vector (1 .. 3);

   function Line_Plane_Intersection (Line_Vector  : in Vector_3D;
                                     Line_Point   : in Vector_3D;
                                     Plane_Normal : in Vector_3D;
                                     Plane_Point  : in Vector_3D)
                                    return Vector_3D
   is
      Diff  : constant Vector_3D := Line_Point - Plane_Point;
      Denom : constant Real      := Line_Vector * Plane_Normal;
   begin
      if Denom = 0.0 then
         raise Constraint_Error with "Line does not intersect plane";
      end if;
      declare
         Scale : constant Real :=
           -Real'(Diff * Plane_Normal) / Denom;
         Point : constant Vector_3D :=
           Diff + Plane_Point + Scale * Line_Vector;
      begin
         return Point;
      end;
   end Line_Plane_Intersection;

   procedure Put (V : in Vector_3D) is
      use Ada.Text_IO, Real_IO;
   begin
      Put ("(");
      Put (V (1));  Put (",");
      Put (V (2));  Put (",");
      Put (V (3));  Put (")");
   end Put;

begin
   Real_IO.Default_Exp := 0;
   Real_IO.Default_Aft := 3;

   Put (Line_Plane_Intersection (Line_Vector  => (0.0, -1.0, -1.0),
                                 Line_Point   => (0.0,  0.0, 10.0),
                                 Plane_Normal => (0.0,  0.0,  1.0),
                                 Plane_Point  => (0.0,  0.0,  5.0)));
end Intersection;
Output:
( 0.000,-5.000, 5.000)

APL

⍝ Find the intersection of a line with a plane
⍝ The intersection I belongs to a line defined by point L and vector V, translates to: 
⍝ A real parameter t exists, that satisfies I = L + tV 
⍝ I belongs to the plan defined by point P and normal vector N. This means that any two points of the plane make a vector 
⍝ normal to vector N. As I and P belong to the plane, the vector IP is normal to N.
⍝ This translates to: The scalar product IP.N = 0.
⍝ (P - I).N = 0 <=> (P - L - tV).N = 0
⍝ Using distributivity, then associativity, the following equations are established:
⍝ (P - L - tV).N = (P - L).N - (tV).N = (P - L).N - t(V.N) = 0
⍝ Which allows to resolve t: t = ((P - L).N) ÷ (V.N)
⍝ In APL, A.B is coded +/A x B
  V  0 ¯1 ¯1
  L  0 0 10
  N  0 0 1
  P  0 0 5 
  dot  { +/  ×  }
  t  ((P - L) dot N) ÷ V dot N
  I  L + t × V
Output:
  I
0 ¯5 5

Arturo

Translation of: Nim
define :vector [x, y, z][]

addv: function [v1 :vector, v2 :vector]->
    to :vector @[v1\x+v2\x, v1\y+v2\y, v1\z+v2\z]

subv: function [v1 :vector, v2 :vector]->
    to :vector @[v1\x-v2\x, v1\y-v2\y, v1\z-v2\z]

mulv: function [v1 :vector, v2 :vector :floating][
    if? is? :vector v2
        -> return sum @[v1\x*v2\x v1\y*v2\y v1\z*v2\z]
    else 
        -> return to :vector @[v1\x*v2, v1\y*v2, v1\z*v2]
]

intersect: function [lV, lP, pV, pP][
    tdenom: mulv pV lV
    if zero? tdenom -> return to :vector @[∞, ∞, ]
    t: (mulv pV subv pP lP) / tdenom
    return addv mulv lV t lP
]

coords: intersect to :vector @[0.0, neg 1.0, neg 1.0]
                  to :vector @[0.0, 0.0, 10.0]
                  to :vector @[0.0, 0.0, 1.0]
                  to :vector @[0.0, 0.0, 5.0]

print ["Intersection at:" coords]
Output:
Intersection at: [x:0.0 y:-5.0 z:5.0]

AutoHotkey

/*
; https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form
l	= line vector
lo	= point on the line
n	= plane normal vector
Po	= point on the plane

if (l . n) = 0						; line and plane are parallel.
if (Po - lo) . n = 0					; line is contained in the plane.

(P - Po) . n = 0					; vector equation of plane.
P = lo + l * d						; vector equation of line.
((lo + l * d) - Po) . n = 0				; Substitute line into plane equation.
(l . n) * d + (lo - Po) . n = 0				; Expanding. 
d = ((Po - lo) . n) / (l . n)				; solving for d.
P = lo + l * ((Po - lo) . n) / (l . n)			; solving P.
*/

intersectPoint(l, lo, n, Pn ){
	if (Vector_dot(Vector_sub(Pn, lo), n) = 0)	; line is contained in the plane
		return [1]
	if (Vector_dot(l, n) = 0)			; line and plane are parallel
		return [0]
    diff  := Vector_Sub(Pn, lo)				; (Po - lo)
    prod1 := Vector_Dot(diff, n)			; ((Po - lo) . n)
    prod2 := Vector_Dot(l, n)				; (l . n)
    d := prod1 / prod2					; d = ((Po - lo) . n) / (l . n)
    return Vector_Add(lo, Vector_Mul(l, d))		; P = lo + l * d
}

Vector_Add(v, w){
    return [v.1+w.1, v.2+w.2, v.3+w.3]
}
Vector_Sub(v, w){
    return [v.1-w.1, v.2-w.2, v.3-w.3]
}
Vector_Mul(v, s){
    return [s*v.1, s*v.2, s*v.3]
}
Vector_Dot(v, w){
    return v.1*w.1 + v.2*w.2 + v.3*w.3
}
Examples:
; task
l1	:= [0, -1, -1]
lo1	:= [0, 0, 10]
n1	:= [0, 0, 1]
Po1	:= [0, 0, 5]

; line on plane
l2	:= [1, 1, 0]
lo2	:= [1, 1, 0]
n2	:= [0, 0, 1]
Po2	:= [5, 1, 0]

; line parallel to plane
l3	:= [1, 1, 0]
lo3	:= [1, 1, 1]
n3	:= [0, 0, 1]
Po3	:= [5, 1, 0]

output := ""
loop 3
{
	result := ""
	ip := intersectPoint(l%A_Index%, lo%A_Index%, n%A_Index%, Po%A_Index%)
	for i, v in ip
		result .= v ", "
	output .= Trim(result, ", ") "`n"
}
MsgBox % output
return
Output:
0.000000, -5.000000, 5.000000
1	; line on plane
0	; line parallel to plane

C

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

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

C#

using System;

namespace FindIntersection {
    class Vector3D {
        private double x, y, z;

        public Vector3D(double x, double y, double z) {
            this.x = x;
            this.y = y;
            this.z = z;
        }

        public static Vector3D operator +(Vector3D lhs, Vector3D rhs) {
            return new Vector3D(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
        }

        public static Vector3D operator -(Vector3D lhs, Vector3D rhs) {
            return new Vector3D(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
        }

        public static Vector3D operator *(Vector3D lhs, double rhs) {
            return new Vector3D(lhs.x * rhs, lhs.y * rhs, lhs.z * rhs);
        }

        public double Dot(Vector3D rhs) {
            return x * rhs.x + y * rhs.y + z * rhs.z;
        }

        public override string ToString() {
            return string.Format("({0:F}, {1:F}, {2:F})", x, y, z);
        }
    }

    class Program {
        static Vector3D IntersectPoint(Vector3D rayVector, Vector3D rayPoint, Vector3D planeNormal, Vector3D planePoint) {
            var diff = rayPoint - planePoint;
            var prod1 = diff.Dot(planeNormal);
            var prod2 = rayVector.Dot(planeNormal);
            var prod3 = prod1 / prod2;
            return rayPoint - rayVector * prod3;
        }

        static void Main(string[] args) {
            var rv = new Vector3D(0.0, -1.0, -1.0);
            var rp = new Vector3D(0.0, 0.0, 10.0);
            var pn = new Vector3D(0.0, 0.0, 1.0);
            var pp = new Vector3D(0.0, 0.0, 5.0);
            var ip = IntersectPoint(rv, rp, pn, pp);
            Console.WriteLine("The ray intersects the plane at {0}", ip);
        }
    }
}
Output:
The ray intersects the plane at (0.00, -5.00, 5.00)

C++

Translation of: Java
#include <iostream>
#include <sstream>

class Vector3D {
public:
	Vector3D(double x, double y, double z) {
		this->x = x;
		this->y = y;
		this->z = z;
	}

	double dot(const Vector3D& rhs) const {
		return x * rhs.x + y * rhs.y + z * rhs.z;
	}

	Vector3D operator-(const Vector3D& rhs) const {
		return Vector3D(x - rhs.x, y - rhs.y, z - rhs.z);
	}

	Vector3D operator*(double rhs) const {
		return Vector3D(rhs*x, rhs*y, rhs*z);
	}

	friend std::ostream& operator<<(std::ostream&, const Vector3D&);

private:
	double x, y, z;
};

std::ostream & operator<<(std::ostream & os, const Vector3D &f) {
	std::stringstream ss;
	ss << "(" << f.x << ", " << f.y << ", " << f.z << ")";
	return os << ss.str();
}

Vector3D intersectPoint(Vector3D rayVector, Vector3D rayPoint, Vector3D planeNormal, Vector3D planePoint) {
	Vector3D diff = rayPoint - planePoint;
	double prod1 = diff.dot(planeNormal);
	double prod2 = rayVector.dot(planeNormal);
	double prod3 = prod1 / prod2;
	return rayPoint - rayVector * prod3;
}

int main() {
	Vector3D rv = Vector3D(0.0, -1.0, -1.0);
	Vector3D rp = Vector3D(0.0, 0.0, 10.0);
	Vector3D pn = Vector3D(0.0, 0.0, 1.0);
	Vector3D pp = Vector3D(0.0, 0.0, 5.0);
	Vector3D ip = intersectPoint(rv, rp, pn, pp);

	std::cout << "The ray intersects the plane at " << ip << std::endl;

	return 0;
}
Output:
The ray intersects the plane at (0, -5, 5)

D

Translation of: Kotlin
import std.stdio;

struct Vector3D {
    private real x;
    private real y;
    private real z;

    this(real x, real y, real z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    auto opBinary(string op)(Vector3D rhs) const {
        static if (op == "+" || op == "-") {
            mixin("return Vector3D(x" ~ op ~ "rhs.x, y" ~ op ~ "rhs.y, z" ~ op ~ "rhs.z);");
        }
    }

    auto opBinary(string op : "*")(real s) const {
        return Vector3D(s*x, s*y, s*z);
    }

    auto dot(Vector3D rhs) const {
        return x*rhs.x + y*rhs.y + z*rhs.z;
    }

    void toString(scope void delegate(const(char)[]) sink) const {
        import std.format;

        sink("(");
        formattedWrite!"%f"(sink, x);
        sink(",");
        formattedWrite!"%f"(sink, y);
        sink(",");
        formattedWrite!"%f"(sink, z);
        sink(")");
    }
}

auto intersectPoint(Vector3D rayVector, Vector3D rayPoint, Vector3D planeNormal, Vector3D planePoint) {
    auto diff = rayPoint - planePoint;
    auto prod1 = diff.dot(planeNormal);
    auto prod2 = rayVector.dot(planeNormal);
    auto prod3 = prod1 / prod2;
    return rayPoint - rayVector * prod3;
}

void main() {
    auto rv = Vector3D(0.0, -1.0, -1.0);
    auto rp = Vector3D(0.0,  0.0, 10.0);
    auto pn = Vector3D(0.0,  0.0,  1.0);
    auto pp = Vector3D(0.0,  0.0,  5.0);
    auto ip = intersectPoint(rv, rp, pn, pp);
    writeln("The ray intersects the plane at ", ip);
}
Output:
The ray intersects the plane at (0.000000,-5.000000,5.000000)

EasyLang

Translation of: Lua
proc minus . l[] r[] res[] .
   len res[] 3
   for i to 3
      res[i] = l[i] - r[i]
   .
.
func dot l[] r[] .
   for i to 3
      res += l[i] * r[i]
   .
   return res
.
proc scale f . l[] .
   for i to 3
      l[i] = l[i] * f
   .
.
proc inter_point rv[] rp[] pn[] pp[] . res[] .
   minus rp[] pp[] dif[]
   prd1 = dot dif[] pn[]
   prd2 = dot rv[] pn[]
   scale (prd1 / prd2) rv[]
   minus rp[] rv[] res[]
.
rv[] = [ 0.0 -1.0 -1.0 ]
rp[] = [ 0.0 0.0 10.0 ]
pn[] = [ 0.0 0.0 1.0 ]
pp[] = [ 0.0 0.0 5.0 ]
inter_point rv[] rp[] pn[] pp[] res[]
print res[]

Evaldraw

Translation of: C

Makes use of the intersectionPoint function to intersect 9 lines with 1 moving plane in a realtime demo.

Grid of 3x3 3d points intersecting a 3D plane
Shows 3x3 grid of lines intersecting a plane. Gridlines drawn between intersection points. Intersection "time" value projected from 3d intersection point to 2d screen rendering.
struct vec{x,y,z;};
enum{GRIDRES=3} // Keep a NxN grid of intersection results.
static vec intersections[GRIDRES][GRIDRES];
static vec ipos = {0,5,-15};
static vec ileft = {-1,0,0};
static vec iup   = {0,-1,0};
static vec ifor  = {0,0,1};
()
{
  cls(0); clz(1e32);
  
  setcam( ipos.x, ipos.y, ipos.z,
    ileft.x, ileft.y, ileft.z, // flip right basis to left
    iup.x, iup.y, iup.z, // flip down basis to up
    ifor.x, ifor.y, ifor.z);
  
  t=klock(0);
  vec planePoint = {0,5,0};   // Plane Position
  vec pN = {cos(t),1,sin(t)}; // PlaneNormal, un-normalized
  normalize(pN);

  for(x=0; x<GRIDRES; x++)
  for(z=0; z<GRIDRES; z++)
  {
     scale = 4.5; halfgrid = scale*(GRIDRES-1)/2;
     vec lineVector = {0,1,0}; // Direction of line
     vec linePoint ={-halfgrid+scale*x, 5, -halfgrid+scale*z};
     
     if (vecdot( lineVector, pN ) == 0 )
     {
       moveto(0,0); printf("Line and Plane dont intersect.");
     } else {
       vec isect;
       isect_time = intersectionPoint(lineVector, linePoint, pN, planePoint, isect);
       intersections[x][z] = isect; // Store for drawing grid
       
       //setcol(255,255,0);  drawsph(isect.x, isect.y, isect.z, .1);
       setcol(255,0,0); line(linePoint, isect);
       unproject(isect);
       setfont(8,12); setcol(255,255,255); printf("t=%2.1f", isect_time);
     }
  }
  // drawgridPlane
  setcol(255,0,255);
  for(i=0; i<GRIDRES; i++)
  for(j=0; j<GRIDRES; j++) {
    vec p00 = intersections[i][j];
    vec p10 = intersections[(i+1)%GRIDRES][j];
    vec p01 = intersections[i][(j+1)%GRIDRES]; // oob wraps to 0 anyhow
    line(p00,p10);
    line(p00,p01);
  }
  setcol(192,192,192); moveto(0,0); printf("Line vs Plane intersection");
}
intersectionPoint(vec lineVector, vec linePoint, vec planeNormal, vec planePoint, vec isect){
   vec diff; vecsub(diff,linePoint,planePoint);
   vec pd; vecadd(pd, diff,planePoint);
   t = -vecdot(diff,planeNormal) / vecdot(lineVector,planeNormal);
   vec scaledVec; vecscalar(scaledVec, lineVector, t);
   vecadd(isect, pd, scaledVec);
   return t;
}
line(vec a, vec b) { moveto(a.x,a.y,a.z); lineto(b.x,b.y,b.z); }
// -------------------------------------- VECTOR MATH
vecScalar( vec out, vec a, s ) {
   out.x = a.x * s;
   out.y = a.y * s;
   out.z = a.z * s;
}
vecAdd( vec out, vec a, vec b) {
   out.x = a.x + b.x;
   out.y = a.y + b.y;
   out.z = a.z + b.z;   
}
vecAdd( vec out, vec b) {
   out.x += b.x;
   out.y += b.y;
   out.z += b.z;   
}
vecSub( vec out, vec a, vec b) {
   out.x = a.x - b.x;
   out.y = a.y - b.y;
   out.z = a.z - b.z;   
}
vecCross( vec out, vec a, vec b) {
   out.x = a.y*b.z - a.z*b.y;
   out.y = a.z*b.x - a.x*b.z;
   out.z = a.x*b.y - a.y*b.x;
}
vecDot( vec a, vec b) {
   return a.x*b.x + a.y*b.y + a.z*b.z;
}
length( vec v ) {
   return sqrt( vecdot(v,v) );
}
normalize( vec v ) {
   len = length(v);
   if ( len ) { v.x /= len; v.y /= len; v.z /= len; }
}
unproject(vec pt) { // unproject a 3D screenpoint
  vec from_eye; vecsub(from_eye, pt, ipos);
  nx = vecdot(from_eye, ileft);
  ny = vecdot(from_eye, iup);
  nz = vecdot(from_eye, ifor);
  if (nz <= 0.5) return; // behind eye
  f = xres/2/nz; // 90 degree projection
  moveto(nx*f + xres/2, ny*f + yres/2 );
}

F#

Translation of: C#
open System

type Vector(x : double, y : double, z : double) =
    member this.x = x
    member this.y = y
    member this.z = z
    static member (-) (lhs : Vector, rhs : Vector) =
        Vector(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z)
    static member (*) (lhs : Vector, rhs : double) =
        Vector(lhs.x * rhs, lhs.y * rhs, lhs.z * rhs)
    override this.ToString() =
        String.Format("({0:F}, {1:F}, {2:F})", x, y, z)

let Dot (lhs:Vector) (rhs:Vector) =
    lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z

let IntersectPoint rayVector rayPoint planeNormal planePoint =
    let diff = rayPoint - planePoint
    let prod1 = Dot diff planeNormal
    let prod2 = Dot rayVector planeNormal
    let prod3 = prod1 / prod2
    rayPoint - rayVector * prod3

[<EntryPoint>]
let main _ = 
    let rv = Vector(0.0, -1.0, -1.0)
    let rp = Vector(0.0, 0.0, 10.0)
    let pn = Vector(0.0, 0.0, 1.0)
    let pp = Vector(0.0, 0.0, 5.0)
    let ip = IntersectPoint rv rp pn pp
    Console.WriteLine("The ray intersects the plane at {0}", ip)

    0 // return an integer exit code
Output:
The ray intersects the plane at (0.00, -5.00, 5.00)

Factor

Translation of: 11l
USING: io locals math.vectors prettyprint ;

:: intersection-point ( rdir rpt pnorm ppt -- loc )
    rpt rdir pnorm rpt ppt v- v. v*n rdir pnorm v. v/n v- ;

"The ray intersects the plane at " write
{ 0 -1 -1 } { 0 0 10 } { 0 0 1 } { 0 0 5 } intersection-point .
Output:
The ray intersects the plane at { 0 -5 5 }

FreeBASIC

' version 11-07-2018
' compile with: fbc -s console

Type vector3d
    Dim As Double x, y ,z
    Declare Constructor ()
    Declare Constructor (ByVal x As Double, ByVal y As Double, ByVal z As Double)
End Type

Constructor vector3d()
    This.x = 0
    This.y = 0
    This.z = 0
End Constructor

Constructor vector3d(ByVal x As Double, ByVal y As Double, ByVal z As Double)
    This.x = x
    This.y = y
    This.z = z
End Constructor

Operator + (lhs As vector3d, rhs As vector3d) As vector3d
    Return Type(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z)
End Operator

Operator - (lhs As vector3d, rhs As vector3d) As vector3d
    Return Type(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z)
End Operator

Operator * (lhs As vector3d, d As Double) As vector3d 
    Return Type(lhs.x * d, lhs.y * d, lhs.z * d)
End Operator

Function dot(lhs As vector3d, rhs As vector3d) As Double
    Return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z
End Function

Function tostring(vec As vector3d) As String
    Return "(" + Str(vec.x) + ", " + Str(vec.y) + ", " + Str(vec.z) + ")"
End Function

Function intersectpoint(rayvector As vector3d, raypoint As vector3d, _
                    planenormal As vector3d, planepoint As vector3d) As vector3d

    Dim As vector3d diff = raypoint - planepoint
    Dim As Double prod1 = dot(diff, planenormal)
    Dim As double prod2 = dot(rayvector, planenormal)
    Return raypoint - rayvector * (prod1 / prod2)

End Function

' ------=< MAIN >=------

Dim As vector3d rv = Type(0, -1, -1)
Dim As vector3d rp = Type(0,  0, 10)
Dim As vector3d pn = Type(0,  0,  1)
Dim As vector3d pp = Type(0,  0,  5)
Dim As vector3d ip = intersectpoint(rv, rp, pn, pp)

print
Print "line intersects the plane at "; tostring(ip)

' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
line intersects the plane at (0, -5, 5)

Go

Translation of: Kotlin
package main

import "fmt"

type Vector3D struct{ x, y, z float64 }

func (v *Vector3D) Add(w *Vector3D) *Vector3D {
    return &Vector3D{v.x + w.x, v.y + w.y, v.z + w.z}
}

func (v *Vector3D) Sub(w *Vector3D) *Vector3D {
    return &Vector3D{v.x - w.x, v.y - w.y, v.z - w.z}
}

func (v *Vector3D) Mul(s float64) *Vector3D {
    return &Vector3D{s * v.x, s * v.y, s * v.z}
}

func (v *Vector3D) Dot(w *Vector3D) float64 {
    return v.x*w.x + v.y*w.y + v.z*w.z
}

func (v *Vector3D) String() string {
    return fmt.Sprintf("(%v, %v, %v)", v.x, v.y, v.z)
}

func intersectPoint(rayVector, rayPoint, planeNormal, planePoint *Vector3D) *Vector3D {
    diff := rayPoint.Sub(planePoint)
    prod1 := diff.Dot(planeNormal)
    prod2 := rayVector.Dot(planeNormal)
    prod3 := prod1 / prod2
    return rayPoint.Sub(rayVector.Mul(prod3))
}

func main() {
    rv := &Vector3D{0.0, -1.0, -1.0}
    rp := &Vector3D{0.0, 0.0, 10.0}
    pn := &Vector3D{0.0, 0.0, 1.0}
    pp := &Vector3D{0.0, 0.0, 5.0}
    ip := intersectPoint(rv, rp, pn, pp)
    fmt.Println("The ray intersects the plane at", ip)
}
Output:
The ray intersects the plane at (0, -5, 5)

Groovy

Translation of: Java
class LinePlaneIntersection {
    private static class Vector3D {
        private double x, y, z

        Vector3D(double x, double y, double z) {
            this.x = x
            this.y = y
            this.z = z
        }

        Vector3D plus(Vector3D v) {
            return new Vector3D(x + v.x, y + v.y, z + v.z)
        }

        Vector3D minus(Vector3D v) {
            return new Vector3D(x - v.x, y - v.y, z - v.z)
        }

        Vector3D multiply(double s) {
            return new Vector3D(s * x, s * y, s * z)
        }

        double dot(Vector3D v) {
            return x * v.x + y * v.y + z * v.z
        }

        @Override
        String toString() {
            return "($x, $y, $z)"
        }
    }

    private static Vector3D intersectPoint(Vector3D rayVector, Vector3D rayPoint, Vector3D planeNormal, Vector3D planePoint) {
        Vector3D diff = rayPoint - planePoint
        double prod1 = diff.dot(planeNormal)
        double prod2 = rayVector.dot(planeNormal)
        double prod3 = prod1 / prod2
        return rayPoint - rayVector * prod3
    }

    static void main(String[] args) {
        Vector3D rv = new Vector3D(0.0, -1.0, -1.0)
        Vector3D rp = new Vector3D(0.0, 0.0, 10.0)
        Vector3D pn = new Vector3D(0.0, 0.0, 1.0)
        Vector3D pp = new Vector3D(0.0, 0.0, 5.0)
        Vector3D ip = intersectPoint(rv, rp, pn, pp)
        println("The ray intersects the plane at $ip")
    }
}
Output:
The ray intersects the plane at (0.0, -5.0, 5.0)

Haskell

Translation of: Kotlin

Note that V3 is implemented similarly in the external library linear.

import Control.Applicative (liftA2)
import Text.Printf (printf)

data V3 a = V3 a a a
    deriving Show

instance Functor V3 where
    fmap f (V3 a b c) = V3 (f a) (f b) (f c)

instance Applicative V3 where
    pure a = V3 a a a
    V3 a b c <*> V3 d e f = V3 (a d) (b e) (c f)

instance Num a => Num (V3 a) where
    (+) = liftA2 (+)
    (-) = liftA2 (-)
    (*) = liftA2 (*)
    negate = fmap negate
    abs = fmap abs
    signum = fmap signum
    fromInteger = pure . fromInteger

dot ::Num a => V3 a -> V3 a -> a
dot a b = x + y + z
  where
    V3 x y z = a * b

intersect :: Fractional a => V3 a -> V3 a -> V3 a -> V3 a -> V3 a
intersect rayVector rayPoint planeNormal planePoint =
    rayPoint - rayVector * pure prod3
  where
    diff = rayPoint - planePoint
    prod1 = diff `dot` planeNormal
    prod2 = rayVector `dot` planeNormal
    prod3 = prod1 / prod2

main = printf "The ray intersects the plane at (%f, %f, %f)\n" x y z
  where
    V3 x y z = intersect rv rp pn pp :: V3 Double
    rv = V3 0 (-1) (-1)
    rp = V3 0 0 10
    pn = V3 0 0 1
    pp = V3 0 0 5
Output:
The ray intersects the plane at (0.0, -5.0, 5.0)

J

Solution:

mp=: +/ .*                          NB. matrix product
p=: mp&{: %~ -~&{. mp {:@]          NB. solve
intersectLinePlane=: [ +/@:* 1 , p  NB. substitute

Example Usage:

   Line=: 0 0 10 ,: 0 _1 _1   NB. Point, Ray
   Plane=: 0 0 5 ,: 0 0 1     NB. Point, Normal
   Line intersectLinePlane Plane
0 _5 5

Java

Translation of: Kotlin
public class LinePlaneIntersection {
    private static class Vector3D {
        private double x, y, z;

        Vector3D(double x, double y, double z) {
            this.x = x;
            this.y = y;
            this.z = z;
        }

        Vector3D plus(Vector3D v) {
            return new Vector3D(x + v.x, y + v.y, z + v.z);
        }

        Vector3D minus(Vector3D v) {
            return new Vector3D(x - v.x, y - v.y, z - v.z);
        }

        Vector3D times(double s) {
            return new Vector3D(s * x, s * y, s * z);
        }

        double dot(Vector3D v) {
            return x * v.x + y * v.y + z * v.z;
        }

        @Override
        public String toString() {
            return String.format("(%f, %f, %f)", x, y, z);
        }
    }

    private static Vector3D intersectPoint(Vector3D rayVector, Vector3D rayPoint, Vector3D planeNormal, Vector3D planePoint) {
        Vector3D diff = rayPoint.minus(planePoint);
        double prod1 = diff.dot(planeNormal);
        double prod2 = rayVector.dot(planeNormal);
        double prod3 = prod1 / prod2;
        return rayPoint.minus(rayVector.times(prod3));
    }

    public static void main(String[] args) {
        Vector3D rv = new Vector3D(0.0, -1.0, -1.0);
        Vector3D rp = new Vector3D(0.0, 0.0, 10.0);
        Vector3D pn = new Vector3D(0.0, 0.0, 1.0);
        Vector3D pp = new Vector3D(0.0, 0.0, 5.0);
        Vector3D ip = intersectPoint(rv, rp, pn, pp);
        System.out.println("The ray intersects the plane at " + ip);
    }
}
Output:
The ray intersects the plane at (0.000000, -5.000000, 5.000000)

jq

Works with: jq

Works with gojq, the Go implementation of jq

In the following, a 3d vector is represented by a JSON array: [x, y, z]

# add as many as you please
def addVector:
  transpose | add;

# . - y
def minusVector(y):
  [.[0] - y[0], .[1] - y[1], .[2] - y[2]];

# scalar multiplication: . * s
def multVector(s):
    map(. * s);

def dot(y):
  .[0] * y[0] + .[1] * y[1] + .[2] * y[2];

def intersectPoint($rayVector; $rayPoint; $planeNormal; $planePoint):
  ($rayPoint | minusVector($planePoint)) as $diff
  | ($diff|dot($planeNormal)) as $prod1
  | ($rayVector|dot($planeNormal)) as $prod2
  | $rayPoint | minusVector($rayVector | multVector(($prod1 / $prod2) )) ;

 
def rv : [0, -1, -1];
def rp : [0,  0, 10];
def pn : [0,  0,  1];
def pp : [0,  0,  5];

"The ray intersects the plane at:",
intersectPoint(rv; rp; pn; pp)
Output:
The ray intersects the plane at:
[0,-5,5]


Julia

Works with: Julia version 0.6
Translation of: Python
function lineplanecollision(planenorm::Vector, planepnt::Vector, raydir::Vector, raypnt::Vector)
    ndotu = dot(planenorm, raydir)
    if ndotu  0 error("no intersection or line is within plane") end

    w  = raypnt - planepnt
    si = -dot(planenorm, w) / ndotu
    ψ  = w .+ si .* raydir .+ planepnt
    return ψ
end

# Define plane
planenorm = Float64[0, 0, 1]
planepnt  = Float64[0, 0, 5]

# Define ray
raydir = Float64[0, -1, -1]
raypnt = Float64[0,  0, 10]

ψ = lineplanecollision(planenorm, planepnt, raydir, raypnt)
println("Intersection at ")
Output:
Intersection at [0.0, -5.0, 5.0]

Kotlin

// version 1.1.51

class Vector3D(val x: Double, val y: Double, val z: Double) {

    operator fun plus(v: Vector3D) = Vector3D(x + v.x, y + v.y, z + v.z)

    operator fun minus(v: Vector3D) = Vector3D(x - v.x, y - v.y, z - v.z)

    operator fun times(s: Double) = Vector3D(s * x, s * y, s * z)

    infix fun dot(v: Vector3D) = x * v.x + y * v.y + z * v.z

    override fun toString() = "($x, $y, $z)"
}

fun intersectPoint(
    rayVector: Vector3D,
    rayPoint: Vector3D,
    planeNormal: Vector3D,
    planePoint: Vector3D
): Vector3D {
    val diff  = rayPoint - planePoint
    val prod1 = diff dot planeNormal
    val prod2 = rayVector dot planeNormal
    val prod3 = prod1 / prod2
    return rayPoint - rayVector * prod3
}

fun main(args: Array<String>) {
    val rv = Vector3D(0.0, -1.0, -1.0)
    val rp = Vector3D(0.0,  0.0, 10.0)
    val pn = Vector3D(0.0,  0.0,  1.0)
    val pp = Vector3D(0.0,  0.0,  5.0)
    val ip = intersectPoint(rv, rp, pn, pp)
    println("The ray intersects the plane at $ip")
}
Output:
The ray intersects the plane at (0.0, -5.0, 5.0)

Lua

function make(xval, yval, zval)
    return {x=xval, y=yval, z=zval}
end

function plus(lhs, rhs)
    return make(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z)
end

function minus(lhs, rhs)
    return make(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z)
end

function times(lhs, scale)
    return make(scale * lhs.x, scale * lhs.y, scale * lhs.z)
end

function dot(lhs, rhs)
    return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z
end

function tostr(val)
    return "(" .. val.x .. ", " .. val.y .. ", " .. val.z .. ")"
end

function intersectPoint(rayVector, rayPoint, planeNormal, planePoint)
    diff = minus(rayPoint, planePoint)
    prod1 = dot(diff, planeNormal)
    prod2 = dot(rayVector, planeNormal)
    prod3 = prod1 / prod2
    return minus(rayPoint, times(rayVector, prod3))
end

rv = make(0.0, -1.0, -1.0)
rp = make(0.0, 0.0, 10.0)
pn = make(0.0, 0.0, 1.0)
pp = make(0.0, 0.0, 5.0)
ip = intersectPoint(rv, rp, pn, pp)
print("The ray intersects the plane at " .. tostr(ip))
Output:
The ray intersects the plane at (0, -5, 5)

Maple

geom3d:-plane(P, [geom3d:-point(p1,0,0,5), [0,0,1]]);
geom3d:-line(L, [geom3d:-point(p2,0,0,10), [0,-1,-1]]);
geom3d:-intersection(px,L,P);
geom3d:-detail(px);
Output:
[["name of the object",px],["form of the object",point3d],["coordinates of the point",[0,-5,5]]]

Mathematica / Wolfram Language

RegionIntersection[InfiniteLine[{0, 0, 10}, {0, -1, -1}], InfinitePlane[{0, 0, 5}, {{0, 1, 0}, {1, 0, 0}}]]
Output:
Point[{0, -5, 5}]


MATLAB

Translation of: Kotlin
function point = intersectPoint(rayVector, rayPoint, planeNormal, planePoint)

pdiff = rayPoint - planePoint;
prod1 = dot(pdiff, planeNormal);
prod2 = dot(rayVector, planeNormal);
prod3 = prod1 / prod2;

point = rayPoint - rayVector * prod3;
Output:
>> intersectPoint([0 -1 -1], [0 0 10], [0 0 1], [0 0 5])

ans =

     0    -5     5

Modula-2

MODULE LinePlane;
FROM RealStr IMPORT RealToStr;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;

TYPE
    Vector3D = RECORD
        x,y,z : REAL;
    END;

PROCEDURE Minus(lhs,rhs : Vector3D) : Vector3D;
VAR out : Vector3D;
BEGIN
    RETURN Vector3D{lhs.x-rhs.x, lhs.y-rhs.y, lhs.z-rhs.z};
END Minus;

PROCEDURE Times(a : Vector3D; s : REAL) : Vector3D;
BEGIN
    RETURN Vector3D{a.x*s, a.y*s, a.z*s};
END Times;

PROCEDURE Dot(lhs,rhs : Vector3D) : REAL;
BEGIN
    RETURN lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
END Dot;

PROCEDURE ToString(self : Vector3D);
VAR buf : ARRAY[0..63] OF CHAR;
BEGIN
    WriteString("(");
    RealToStr(self.x,buf);
    WriteString(buf);
    WriteString(", ");
    RealToStr(self.y,buf);
    WriteString(buf);
    WriteString(", ");
    RealToStr(self.z,buf);
    WriteString(buf);
    WriteString(")");
END ToString;

PROCEDURE IntersectPoint(rayVector,rayPoint,planeNormal,planePoint : Vector3D) : Vector3D;
VAR
    diff : Vector3D;
    prod1,prod2,prod3 : REAL;
BEGIN
    diff := Minus(rayPoint,planePoint);
    prod1 := Dot(diff, planeNormal);
    prod2 := Dot(rayVector, planeNormal);
    prod3 := prod1 / prod2;
    RETURN Minus(rayPoint, Times(rayVector, prod3));
END IntersectPoint;

VAR ip : Vector3D;
BEGIN
    ip := IntersectPoint(Vector3D{0.0,-1.0,-1.0},Vector3D{0.0,0.0,10.0},Vector3D{0.0,0.0,1.0},Vector3D{0.0,0.0,5.0});

    WriteString("The ray intersects the plane at ");
    ToString(ip);
    WriteLn;

    ReadChar;
END LinePlane.

Nim

type Vector = tuple[x, y, z: float]


func `+`(v1, v2: Vector): Vector =
  ## Add two vectors.
  (v1.x + v2.x, v1.y + v2.y, v1.z + v2.z)

func `-`(v1, v2: Vector): Vector =
  ## Subtract a vector to another one.
  (v1.x - v2.x, v1.y - v2.y, v1.z - v2.z)

func `*`(v1, v2: Vector): float =
  ## Compute the dot product of two vectors.
  v1.x * v2.x + v1.y * v2.y + v1.z * v2.z

func `*`(v: Vector; k: float): Vector =
  ## Multiply a vector by a scalar.
  (v.x * k, v.y * k, v.z * k)


func intersection(lineVector, linePoint, planeVector, planePoint: Vector): Vector =
  ## Return the coordinates of the intersection of two vectors.

  let tnum = planeVector * (planePoint - linePoint)
  let tdenom = planeVector * lineVector
  if tdenom == 0: return (Inf, Inf, Inf)  # No intersection.
  let t = tnum / tdenom
  result = lineVector * t + linePoint

let coords = intersection(lineVector = (0.0, -1.0, -1.0),
                          linePoint = (0.0, 0.0, 10.0),
                          planeVector = (0.0, 0.0, 1.0),
                          planePoint = (0.0, 0.0, 5.0))
echo "Intersection at ", coords
Output:
Intersection at (x: 0.0, y: -5.0, z: 5.0)

Perl

Translation of: Raku
package  Line; sub new { my ($c, $a) = @_; my $self = { P0 => $a->{P0}, u => $a->{u} } } # point / ray
package Plane; sub new { my ($c, $a) = @_; my $self = { V0 => $a->{V0}, n => $a->{n} } } # point / normal

package main;

sub dot { my $p; $p    += $_[0][$_] * $_[1][$_] for 0..@{$_[0]}-1; $p } # dot product
sub vd  { my @v; $v[$_] = $_[0][$_] - $_[1][$_] for 0..@{$_[0]}-1; @v } # vector difference
sub va  { my @v; $v[$_] = $_[0][$_] + $_[1][$_] for 0..@{$_[0]}-1; @v } # vector addition
sub vp  { my @v; $v[$_] = $_[0][$_] * $_[1][$_] for 0..@{$_[0]}-1; @v } # vector product

sub line_plane_intersection {
    my($L, $P) = @_;

    my $cos = dot($L->{u}, $P->{n});     # cosine between normal & ray
    return 'Vectors are orthogonol; no intersection or line within plane' if $cos == 0;

    my @W = vd($L->{P0},$P->{V0});       # difference between P0 and V0
    my $SI = -dot($P->{n}, \@W) / $cos;  # line segment where it intersets the plane

    my @a = vp($L->{u},[($SI)x3]);
    my @b = va($P->{V0},\@a);
    va(\@W,\@b);
}

my $L =  Line->new({ P0=>[0,0,10], u=>[0,-1,-1]});
my $P = Plane->new({ V0=>[0,0,5 ], n=>[0, 0, 1]});
print 'Intersection at point: ', join(' ', line_plane_intersection($L, $P)) . "\n";
Output:
Intersection at point: 0 -5 5

Phix

with javascript_semantics
function dot(sequence a, b) return sum(sq_mul(a,b)) end function
 
function intersection_point(sequence line_vector,line_point,plane_normal,plane_point)
    atom a = dot(line_vector,plane_normal)
    if a=0 then return "no intersection" end if
    sequence diff = sq_sub(line_point,plane_point)
    return sq_add(sq_add(diff,plane_point),sq_mul(-dot(diff,plane_normal)/a,line_vector))
end function
 
?intersection_point({0,-1,-1},{0,0,10},{0,0,1},{0,0,5})
?intersection_point({3,2,1},{0,2,4},{1,2,3},{3,3,3})
?intersection_point({1,1,0},{0,0,1},{0,0,3},{0,0,0}) -- (parallel to plane)
?intersection_point({1,1,0},{1,1,0},{0,0,3},{0,0,0}) -- (line within plane)
Output:
{0,-5,5}
{0.6,2.4,4.2}
"no intersection"
"no intersection"

Picat

Translation of: Java
Works with: Picat
plus(U, V) = {U[1] + V[1], U[2] + V[2], U[3] + V[3]}.

minus(U, V) = {U[1] - V[1], U[2] - V[2], U[3] - V[3]}.

times(U, S) = {U[1] * S, U[2] * S, U[3] * S}.

dot(U, V) = U[1] * V[1] + U[2] * V[2] + U[3] * V[3].

intersect_point(RayVector, RayPoint, PlaneNormal, PlanePoint) = IntersectPoint =>
    Diff = minus(RayPoint, PlanePoint),
    Prod1 = dot(Diff, PlaneNormal),
    Prod2 = dot(RayVector, PlaneNormal),
    Prod3 = Prod1 / Prod2,
    IntersectPoint = minus(RayPoint, times(RayVector, Prod3)).

main =>
    RayVector = {0.0, -1.0, -1.0},
    RayPoint = {0.0, 0.0, 10.0},
    PlaneNormal = {0.0, 0.0, 1.0},
    PlanePoint = {0.0, 0.0, 5.0},
    IntersectPoint = intersect_point(RayVector, RayPoint, PlaneNormal, PlanePoint),
    printf("The ray intersects the plane at (%f, %f, %f)\n",
        IntersectPoint[1],
        IntersectPoint[2],
        IntersectPoint[3]
    ).
Output:
The ray intersects the plane at (0.000000, -5.000000, 5.000000)


Prolog

Translation of: Picat
Works with: GNU Prolog
Works with: SWI Prolog
:- initialization(main).

vector_plus(U, V, W) :-
    U = p(X1, Y1, Z1),
    V = p(X2, Y2, Z2),
    X3 is X1 + X2,
    Y3 is Y1 + Y2,
    Z3 is Z1 + Z2,
    W = p(X3, Y3, Z3).

vector_minus(U, V, W) :-
    U = p(X1, Y1, Z1),
    V = p(X2, Y2, Z2),
    X3 is X1 - X2,
    Y3 is Y1 - Y2,
    Z3 is Z1 - Z2,
    W = p(X3, Y3, Z3).

vector_times(U, S, V) :-
    U = p(X1, Y1, Z1),
    X2 is X1 * S,
    Y2 is Y1 * S,
    Z2 is Z1 * S,
    V = p(X2, Y2, Z2).

vector_dot(U, V, S) :-
    U = p(X1, Y1, Z1),
    V = p(X2, Y2, Z2),
    S is X1 * X2 + Y1 * Y2 + Z1 * Z2.

intersect_point(RayVector, RayPoint, PlaneNormal, PlanePoint, IntersectPoint) :-
    vector_minus(RayPoint, PlanePoint, Diff),
    vector_dot(Diff, PlaneNormal, Prod1),
    vector_dot(RayVector, PlaneNormal, Prod2),
    Prod3 is Prod1 / Prod2,
    vector_times(RayVector, Prod3, Times),
    vector_minus(RayPoint, Times, IntersectPoint).

main :-
    RayVector = p(0.0, -1.0, -1.0),
    RayPoint = p(0.0, 0.0, 10.0),
    PlaneNormal = p(0.0, 0.0, 1.0),
    PlanePoint = p(0.0, 0.0, 5.0),
    intersect_point(RayVector, RayPoint, PlaneNormal, PlanePoint, p(X, Y, Z)),
    format("The ray intersects the plane at (~f, ~f, ~f)\n", [X, Y, Z]).
Output:

The ray intersects the plane at (0.000000, -5.000000, 5.000000)

Python

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]

R

Translation of: MATLAB
intersect_point <- function(ray_vec, ray_point, plane_normal, plane_point) {

  pdiff <- ray_point - plane_point
  prod1 <- pdiff %*% plane_normal
  prod2 <- ray_vec %*% plane_normal
  prod3 <- prod1 / prod2
  point <- ray_point - ray_vec * as.numeric(prod3)

  return(point)
}
Output:
>>intersect_point(c(0, -1, -1), c(0, 0, 10), c(0, 0, 1), c(0, 0, 5))
[1]  0 -5  5

Racket

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!

Raku

(formerly Perl 6)

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)

With a geometric algebra library

See task geometric algebra

use Clifford:ver<6.2.1>;

# We pick a (non-degenerate) projective basis and
# we define the dual and meet operators.
my $I = [∧] my ($i, $j, $k, $l) = @e;
sub prefix:<∗>($M) { $M/$I }
sub infix:<∨>($A, $B) { ∗((∗$B)∧(∗$A)) }

my $direction = -$j - $k;

# Homogeneous coordinates of (X, Y, Z) are (X, Y, Z, 1)
my $point = 10*$k + $l;

# A projective line is a bivector
my $line = $direction$point;

# A projective plane is a trivector
my $plane = (5*$k + $l) ∧ ($k*-$i$j$k);

# The intersection is the meet
my $m = $line$plane;

# Affine coordinates of (X, Y, Z, W) are (X/W, Y/W, Z/W)
say $m/($m·$l) X· ($i, $j, $k);
Output:
(0 -5 5)

REXX

version 1

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

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)

RPL

≪ → rd rp pn pp 
  ≪ rd rp pp - pn DOT * rd pn DOT /
≫ ≫ 'INTLP' STO
[ 0 -1 -1 ] [ 0 0 0 ] [ 0 0 1 ] [ 0 0 5 ] INTLP
Output:
1: [ 0 -5 -5 ]

Ruby

Translation of: C#
require "matrix"

def intersectPoint(rayVector, rayPoint, planeNormal, planePoint)
    diff = rayPoint - planePoint
    prod1 = diff.dot planeNormal
    prod2 = rayVector.dot planeNormal
    prod3 = prod1 / prod2
    return rayPoint - rayVector * prod3
end

def main
    rv = Vector[0.0, -1.0, -1.0]
    rp = Vector[0.0, 0.0, 10.0]
    pn = Vector[0.0, 0.0, 1.0]
    pp = Vector[0.0, 0.0, 5.0]
    ip = intersectPoint(rv, rp, pn, pp)
    puts "The ray intersects the plane at %s" % [ip]
end

main()
Output:
The ray intersects the plane at Vector[0.0, -5.0, 5.0]

Rust

Translation of: Kotlin
use std::ops::{Add, Div, Mul, Sub};

#[derive(Copy, Clone, Debug, PartialEq)]
struct V3<T> {
    x: T,
    y: T,
    z: T,
}

impl<T> V3<T> {
    fn new(x: T, y: T, z: T) -> Self {
        V3 { x, y, z }
    }
}

fn zip_with<F, T, U>(f: F, a: V3<T>, b: V3<T>) -> V3<U>
where
    F: Fn(T, T) -> U,
{
    V3 {
        x: f(a.x, b.x),
        y: f(a.y, b.y),
        z: f(a.z, b.z),
    }
}

impl<T> Add for V3<T>
where
    T: Add<Output = T>,
{
    type Output = Self;

    fn add(self, other: Self) -> Self {
        zip_with(<T>::add, self, other)
    }
}

impl<T> Sub for V3<T>
where
    T: Sub<Output = T>,
{
    type Output = Self;

    fn sub(self, other: Self) -> Self {
        zip_with(<T>::sub, self, other)
    }
}

impl<T> Mul for V3<T>
where
    T: Mul<Output = T>,
{
    type Output = Self;

    fn mul(self, other: Self) -> Self {
        zip_with(<T>::mul, self, other)
    }
}

impl<T> V3<T>
where
    T: Mul<Output = T> + Add<Output = T>,
{
    fn dot(self, other: Self) -> T {
        let V3 { x, y, z } = self * other;
        x + y + z
    }
}

impl<T> V3<T>
where
    T: Mul<Output = T> + Copy,
{
    fn scale(self, scalar: T) -> Self {
        self * V3 {
            x: scalar,
            y: scalar,
            z: scalar,
        }
    }
}

fn intersect<T>(
    ray_vector: V3<T>,
    ray_point: V3<T>,
    plane_normal: V3<T>,
    plane_point: V3<T>,
) -> V3<T>
where
    T: Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + Copy,
{
    let diff = ray_point - plane_point;
    let prod1 = diff.dot(plane_normal);
    let prod2 = ray_vector.dot(plane_normal);
    let prod3 = prod1 / prod2;
    ray_point - ray_vector.scale(prod3)
}

fn main() {
    let rv = V3::new(0.0, -1.0, -1.0);
    let rp = V3::new(0.0, 0.0, 10.0);
    let pn = V3::new(0.0, 0.0, 1.0);
    let pp = V3::new(0.0, 0.0, 5.0);
    println!("{:?}", intersect(rv, rp, pn, pp));
}

Scala

object LinePLaneIntersection extends App {
  val (rv, rp, pn, pp) =
    (Vector3D(0.0, -1.0, -1.0), Vector3D(0.0, 0.0, 10.0), Vector3D(0.0, 0.0, 1.0), Vector3D(0.0, 0.0, 5.0))
  val ip = intersectPoint(rv, rp, pn, pp)

  def intersectPoint(rayVector: Vector3D,
                     rayPoint: Vector3D,
                     planeNormal: Vector3D,
                     planePoint: Vector3D): Vector3D = {
    val diff = rayPoint - planePoint
    val prod1 = diff dot planeNormal
    val prod2 = rayVector dot planeNormal
    val prod3 = prod1 / prod2
    rayPoint - rayVector * prod3
  }

  case class Vector3D(x: Double, y: Double, z: Double) {
    def +(v: Vector3D) = Vector3D(x + v.x, y + v.y, z + v.z)
    def -(v: Vector3D) = Vector3D(x - v.x, y - v.y, z - v.z)
    def *(s: Double) = Vector3D(s * x, s * y, s * z)
    def dot(v: Vector3D): Double = x * v.x + y * v.y + z * v.z
    override def toString = s"($x, $y, $z)"
  }
  
  println(s"The ray intersects the plane at $ip")
}
Output:
See it in running in your browser by ScalaFiddle (JavaScript).

Sidef

Translation of: Raku
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]

Visual Basic .NET

Translation of: C#
Module Module1

    Class Vector3D
        Private ReadOnly x As Double
        Private ReadOnly y As Double
        Private ReadOnly z As Double

        Sub New(nx As Double, ny As Double, nz As Double)
            x = nx
            y = ny
            z = nz
        End Sub

        Public Function Dot(rhs As Vector3D) As Double
            Return x * rhs.x + y * rhs.y + z * rhs.z
        End Function

        Public Shared Operator +(ByVal a As Vector3D, ByVal b As Vector3D) As Vector3D
            Return New Vector3D(a.x + b.x, a.y + b.y, a.z + b.z)
        End Operator

        Public Shared Operator -(ByVal a As Vector3D, ByVal b As Vector3D) As Vector3D
            Return New Vector3D(a.x - b.x, a.y - b.y, a.z - b.z)
        End Operator

        Public Shared Operator *(ByVal a As Vector3D, ByVal b As Double) As Vector3D
            Return New Vector3D(a.x * b, a.y * b, a.z * b)
        End Operator

        Public Overrides Function ToString() As String
            Return String.Format("({0:F}, {1:F}, {2:F})", x, y, z)
        End Function
    End Class

    Function IntersectPoint(rayVector As Vector3D, rayPoint As Vector3D, planeNormal As Vector3D, planePoint As Vector3D) As Vector3D
        Dim diff = rayPoint - planePoint
        Dim prod1 = diff.Dot(planeNormal)
        Dim prod2 = rayVector.Dot(planeNormal)
        Dim prod3 = prod1 / prod2
        Return rayPoint - rayVector * prod3
    End Function

    Sub Main()
        Dim rv = New Vector3D(0.0, -1.0, -1.0)
        Dim rp = New Vector3D(0.0, 0.0, 10.0)
        Dim pn = New Vector3D(0.0, 0.0, 1.0)
        Dim pp = New Vector3D(0.0, 0.0, 5.0)
        Dim ip = IntersectPoint(rv, rp, pn, pp)
        Console.WriteLine("The ray intersects the plane at {0}", ip)
    End Sub

End Module
Output:
The ray intersects the plane at (0.00, -5.00, 5.00)

Wren

Translation of: Kotlin
Library: Wren-vector
import "./vector" for Vector3

var intersectPoint = Fn.new { |rayVector, rayPoint, planeNormal, planePoint|
    var diff  = rayPoint - planePoint
    var prod1 = diff.dot(planeNormal)
    var prod2 = rayVector.dot(planeNormal)
    var prod3 = prod1 / prod2
    return rayPoint - rayVector*prod3
}

var rv = Vector3.new(0, -1, -1)
var rp = Vector3.new(0,  0, 10)
var pn = Vector3.new(0,  0,  1)
var pp = Vector3.new(0,  0,  5)
var ip = intersectPoint.call(rv, rp, pn, pp)
System.print("The ray intersects the plane at %(ip).")
Output:
The ray intersects the plane at (0, -5, 5).

XPL0

Translation of: Wren
include xpllib;

func real IntersectPoint; real RayVector, RayPoint, PlaneNormal, PlanePoint;
real Diff(3), Prod1, Prod2, Prod3, Prod(3);
[VSub(Diff, RayPoint, PlanePoint);
Prod1:= VDot(Diff, PlaneNormal);
Prod2:= VDot(RayVector, PlaneNormal);
Prod3:= Prod1 / Prod2;
return VSub(Diff, RayPoint, VMul(Prod, RayVector, Prod3));
];

real RV, RP, PN, PP, IP;
[RV:= [0., -1., -1.];
 RP:= [0.,  0., 10.];
 PN:= [0.,  0.,  1.];
 PP:= [0.,  0.,  5.];
 IP:= IntersectPoint(RV, RP, PN, PP);
Print("The ray intersects the plane at %1.1f, %1.1f, %1.1f\n", IP(0), IP(1), IP(2));
]
Output:
The ray intersects the plane at 0.0, -5.0, 5.0

zkl

Translation of: Raku
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 orthogonol; 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