Find the intersection of a line with a plane
Finding the intersection of an infinite ray with a plane in 3D is an important topic in collision detection.
You are encouraged to solve this task according to the task description, using any language you may know.
- 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
<lang 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)))</lang>
- Output:
The ray intersects the plane at (0, -5, 5)
Action!
<lang 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</lang>
- 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
<lang 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;</lang>
- Output:
( 0.000,-5.000, 5.000)
APL
<lang 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
</lang>
- Output:
I 0 ¯5 5
AutoHotkey
<lang AutoHotkey>/*
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
}</lang> Examples:<lang AutoHotkey>; 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</lang>
- 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. <lang C>
- 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; } </lang> 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#
<lang csharp>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);
}
}
}</lang>
- Output:
The ray intersects the plane at (0.00, -5.00, 5.00)
C++
<lang cpp>#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; }</lang>
- Output:
The ray intersects the plane at (0, -5, 5)
D
<lang D>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);
}</lang>
- Output:
The ray intersects the plane at (0.000000,-5.000000,5.000000)
F#
<lang fsharp>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</lang>
- Output:
The ray intersects the plane at (0.00, -5.00, 5.00)
Factor
<lang factor>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 .</lang>
- Output:
The ray intersects the plane at { 0 -5 5 }
FreeBASIC
<lang 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</lang>
- Output:
line intersects the plane at (0, -5, 5)
Go
<lang go>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)
}</lang>
- Output:
The ray intersects the plane at (0, -5, 5)
Groovy
<lang groovy>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")
}
}</lang>
- Output:
The ray intersects the plane at (0.0, -5.0, 5.0)
Haskell
Note that V3 is implemented similarly in the external library linear. <lang Haskell>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</lang>
- Output:
The ray intersects the plane at (0.0, -5.0, 5.0)
J
Solution: <lang j>mp=: +/ .* NB. matrix product p=: mp&{: %~ -~&{. mp {:@] NB. solve intersectLinePlane=: [ +/@:* 1 , p NB. substitute</lang> Example Usage: <lang j> 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</lang>
Java
<lang Java>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);
}
}</lang>
- Output:
The ray intersects the plane at (0.000000, -5.000000, 5.000000)
jq
Works with gojq, the Go implementation of jq
In the following, a 3d vector is represented by a JSON array: [x, y, z] <lang jq># 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)</lang>
- Output:
The ray intersects the plane at: [0,-5,5]
Julia
<lang julia>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 $ψ")</lang>
- Output:
Intersection at [0.0, -5.0, 5.0]
Kotlin
<lang scala>// 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")
}</lang>
- Output:
The ray intersects the plane at (0.0, -5.0, 5.0)
Lua
<lang 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))</lang>
- Output:
The ray intersects the plane at (0, -5, 5)
Maple
<lang 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);</lang>
- Output:
[["name of the object",px],["form of the object",point3d],["coordinates of the point",[0,-5,5]]]
Mathematica / Wolfram Language
<lang Mathematica>RegionIntersection[InfiniteLine[{0, 0, 10}, {0, -1, -1}], InfinitePlane[{0, 0, 5}, {{0, 1, 0}, {1, 0, 0}}]]</lang>
- Output:
Point[{0, -5, 5}]
MATLAB
<lang MATLAB>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;</lang>
- Output:
<lang MATLAB>>> intersectPoint([0 -1 -1], [0 0 10], [0 0 1], [0 0 5])
ans =
0 -5 5
</lang>
Modula-2
<lang modula2>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.</lang>
Nim
<lang 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</lang>
- Output:
Intersection at (x: 0.0, y: -5.0, z: 5.0)
Perl
<lang perl>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"; </lang>
- 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
<lang 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]
).
</lang>
- Output:
The ray intersects the plane at (0.000000, -5.000000, 5.000000)
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]
R
<lang R>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)
}</lang>
- Output:
<lang R>>>intersect_point(c(0, -1, -1), c(0, 0, 10), c(0, 0, 1), c(0, 0, 5)) [1] 0 -5 5</lang>
Racket
<lang racket>#lang racket
- 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!
Raku
(formerly Perl 6)
<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 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) ) );</lang>
- Output:
Intersection at point: (0 -5 5)
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 within it. <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 */ 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 </lang>
- 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)
Ruby
<lang ruby>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()</lang>
- Output:
The ray intersects the plane at Vector[0.0, -5.0, 5.0]
Rust
<lang Rust>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 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));
} </lang>
Scala
<lang 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")
}</lang>
- Output:
See it in running in your browser by ScalaFiddle (JavaScript).
Sidef
<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]
Visual Basic .NET
<lang vbnet>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</lang>
- Output:
The ray intersects the plane at (0.00, -5.00, 5.00)
Wren
<lang ecmascript>class Vector3D {
construct new(x, y, z) {
_x = x
_y = y
_z = z
}
x { _x }
y { _y }
z { _z }
+(v) { Vector3D.new(_x + v.x, _y + v.y, _z + v.z) }
-(v) { Vector3D.new(_x - v.x, _y - v.y, _z - v.z) }
*(s) { Vector3D.new(s * _x, s * _y, s * _z) }
dot(v) { _x * v.x + _y * v.y + _z * v.z }
toString { "(%(_x), %(_y), %(_z))" }
}
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 = Vector3D.new(0, -1, -1) var rp = Vector3D.new(0, 0, 10) var pn = Vector3D.new(0, 0, 1) var pp = Vector3D.new(0, 0, 5) var ip = intersectPoint.call(rv, rp, pn, pp) System.print("The ray intersects the plane at %(ip).")</lang>
- Output:
The ray intersects the plane at (0, -5, 5).
zkl
<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 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);
}</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