Find the intersection of a line with a plane: Difference between revisions
Walterpachl (talk | contribs) (add REXX) |
Walterpachl (talk | contribs) m (→{{header|REXX}}: disclaimer for line being parallel to or within the plane) |
||
Line 103: | Line 103: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
This program does NOT handle the case when the line is parallel to or within the plane. |
|||
<lang rexx>/* REXX */ |
<lang rexx>/* REXX */ |
||
Parse Value '0 0 1' With n.1 n.2 n.3 /* Normal Vector of the plane */ |
Parse Value '0 0 1' With n.1 n.2 n.3 /* Normal Vector of the plane */ |
||
Line 124: | Line 125: | ||
{{out}} |
{{out}} |
||
<pre>0*x + 0*y + 1*z = 5 |
<pre>0*x + 0*y + 1*z = 5 |
||
Intersection: P(0,-5,5)</pre> |
Intersection: P(0,-5,5)</pre> |
||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
Revision as of 21:09, 3 July 2017
Finding the intersection of an infinite ray with a plane in 3D is an important topic in collision detection.
- Task
Find the point of intersection for the infinite ray with direction (0,-1,-1) passing through position (0, 0, 10) with the infinite plane with a normal vector of (0, 0, 1) and which passes through [0, 0, 5].
Perl 6
<lang perl6>class Line {
has $.P0; # point has $.u⃗; # ray
} class Plane {
has $.V0; # point has $.n⃗; # normal
}
sub infix:<∙> ( @a, @b where +@a == +@b ) { [+] @a «*» @b } # dot product
sub line-plane-intersection ($𝑳, $𝑷) {
my $cos = $𝑷.n⃗ ∙ $𝑳.u⃗; # cosine between normal & ray return 'Vectors are orthoganol; no intersection or line within plane' if $cos == 0; my $𝑊 = $𝑳.P0 «-» $𝑷.V0; # difference between P0 and V0 my $S𝐼 = -($𝑷.n⃗ ∙ $𝑊) / $cos; # line segment where it intersets the plane $𝑊 «+» $S𝐼 «*» $𝑳.u⃗ «+» $𝑷.V0; # point where line intersects the plane }
say 'Intersection at point: ', line-plane-intersection(
Line.new( :P0(0,0,10), :u⃗(0,-1,-1) ), Plane.new( :V0(0,0, 5), :n⃗(0, 0, 1) ) );</lang>
- Output:
Intersection at point: (0 -5 5)
Python
Based on the approach at geomalgorithms.com[1]
<lang python>#!/bin/python from __future__ import print_function import numpy as np
def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
ndotu = planeNormal.dot(rayDirection) if abs(ndotu) < epsilon: raise RuntimeError("no intersection or line is within plane")
w = rayPoint - planePoint si = -planeNormal.dot(w) / ndotu Psi = w + si * rayDirection + planePoint return Psi
if __name__=="__main__":
#Define plane
planeNormal = np.array([0, 0, 1])
planePoint = np.array([0, 0, 5]) #Any point on the plane
#Define ray rayDirection = np.array([0, -1, -1]) rayPoint = np.array([0, 0, 10]) #Any point along the ray
Psi = LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint) print ("intersection at", Psi)</lang>
- Output:
intersection at [ 0 -5 5]
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!
REXX
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)
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]
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 orthoganol; no intersection or line within plane"); w:=line.pt.zipWith('-,plane.pt); # difference between P0 and V0 si:=-dotP(plane.normal,w)/cos; # line segment where it intersets the plane # point where line intersects the plane: //w.zipWith('+,line.ray.apply('*,si)).zipWith('+,plane.pt); // or w.zipWith('wrap(w,r,pt){ w + r*si + pt },line.ray,plane.pt);
}</lang> <lang zkl>println("Intersection at point: ", linePlaneIntersection(
Line( T(0.0, 0.0, 10.0), T(0.0, -1.0, -1.0) ), Plane(T(0.0, 0.0, 5.0), T(0.0, 0.0, 1.0) ))
);</lang>
- Output:
Intersection at point: L(0,-5,5)