PDA

View Full Version : Ray-plane collision

User137
16-03-2008, 02:55 PM
Anyone have source code or tips how to make function that finds out intersection point with line and plane in 3D space? We know line's start point and direction vectors and also position and normal vector for plane.

This is what i have come up with but it doesn't work when i move either line or plane from 0,0,0 point.

function Dot(const a,b: P3f): single;
begin
result:=a.x*b.x + a.y*b.y + a.z*b.z;
end;

function vDec(const a,b: P3f): P3f;
begin
result.x:=a.x-b.x; result.y:=a.y-b.y; result.z:=a.z-b.z;
end;

procedure rayPlaneCollision(const rayOrigin, rayDirection,
planeOrigin, planeNormal: P3f; intersection: pP3f);
var u,d: single;
begin
d:=dot(planeNormal,rayDirection);
if d=0 then exit; // No intersection
u:=dot(planeNormal,vDec(planeOrigin,rayOrigin)) /d;
intersection.x:=rayOrigin.x+u*rayDirection.x;
intersection.y:=rayOrigin.y+u*rayDirection.y;
intersection.z:=rayOrigin.z+u*rayDirection.z;
end;

Update: It does work when planeOrigin and rayorigin are at exact same (anywhere) point, but that is cheating because that is also the same as intersection point.

cronodragon
16-03-2008, 09:04 PM
I have the Segment-Plane collision function which I translated from C to Pascal, from the book Real-Time Collision Detection. You might be able to emulate the ray using a very large segment:

//------------------------------------------------------------------------------
function ch_IntersecaSegmentoPlano&#40;
const S1, S2&#58; gr_rVector3D; // Puntos del segmento
const P1&#58; gr_rPlano3D;
var T&#58; Single;
var P&#58; gr_rVector3D
&#41;&#58; Boolean;

var
AB&#58; gr_rVector3D;

begin
Result &#58;= True;
AB &#58;= gr_Resta&#40;S2, S1&#41;;
T &#58;= &#40;P1.PP-gr_ProductoPunto&#40;P1.N, S1&#41;&#41;/gr_ProductoPunto&#40;P1.N, AB&#41;;

if &#40;T >= 0&#41; and &#40;T <= 1&#41; then
begin
P &#58;= gr_Suma&#40;S1, gr_Multiplica&#40;AB, T&#41;&#41;;
Exit;
end;

Result &#58;= False;
end;

Sorry for the spanish, but it should be easy to understand (if you don't I can translate it). I highly recomend that book that explains very well most of the collision routines.

User137
16-03-2008, 11:21 PM
It looks exactly same function as i have, just more wrapped in functions and added range check at the end (0..1 is the range between line start and end). I'm wondering if i have a typo somewhere... I did it according to more than 1 math site, just couldn't find any source code.

Btw in your function, if i assume right, you may get same division by zero which is handled in my code.

cronodragon
17-03-2008, 01:53 AM
You are right! The author has a note about that division by zero problem. I just translated the content of the book, but haven't deeply tested it. There is a entire chapter at the end of the book dealing with those errors. It uses integer arithmetic, this is my new translation:

function ch_IntersecaSegmentoPlano(
const S1, S2: gr_rVector3D; // Puntos del segmento
const P1: gr_rPlano3D;
var T: Single;
var P: gr_rVector3D
): Boolean;

var
AB: gr_rVector3D;
TN, TD, K: Int64;
V: gr_rVector3D;

begin
Result := False;

{ versi??n de la p?°gina 176:
AB := gr_Resta(S2, S1);
T := (P1.PP-gr_ProductoPunto(P1.N, S1))/gr_ProductoPunto(P1.N, AB);

if (T >= 0) and (T <= 1) then
begin
P := gr_Suma(S1, gr_Multiplica(AB, T));
Exit;
end;

Result := False;
}
AB := gr_Resta(S2, S1);
TN := Trunc(P1.PP-gr_ProductoPunto(P1.N, S1));
TD := Trunc(gr_ProductoPunto(P1.N, AB));

if TD = 0 then Exit;

if TD < 0 then
begin
TN := -TN;
TD := -TD;
end;

if (TN < 0) or (TN > TD) then Exit;

V := gr_CreaVector3D;
K := TD-1;

if TD > 0 then K := -K;
if P1.N.X > 0 then V.X := +K else
if P1.N.X < 0 then V.X := -K;
if P1.N.Y > 0 then V.Y := +K else
if P1.N.Y < 0 then V.Y := -K;
if P1.N.Z > 0 then V.Z := +K else
if P1.N.Z < 0 then V.Z := -K;

P := gr_Suma(S1, gr_Divide(gr_Suma(gr_Multiplica(AB, TN), V), TD));
Result := True;
end;

Please let me know if it works for you, I don't want to have that bug in my code neither.

paul_nicholls
17-03-2008, 02:52 AM
For your interest, I found a raytracing site (http://www.devmaster.net/articles/raytracing_series/part1.php) that includes ray-plane intersections if this helps :-)

I have made a Delphi version of the raytracer so I know the code works.

cheers,
Paul

User137
17-03-2008, 08:45 AM
You are right! The author has a note about that division by zero problem. I just translated the content of the book, but haven't deeply tested it. There is a entire chapter at the end of the book dealing with those errors. It uses integer arithmetic, this is my new translation:

function ch_IntersecaSegmentoPlano(
...
end;

Please let me know if it works for you, I don't want to have that bug in my code neither.

I have few questions to make it work. What does gr_rPlano3D type contain? How do i use its .PP variable? And why int64 instead of floating point single?
Edit: Oh, think it has 2 vectors in it, N for normal and PP for position :)

int PlanePrim&#58;&#58;Intersect&#40; Ray& a_Ray, float& a_Dist &#41;
&#123;
float d = DOT&#40; m_Plane.N, a_Ray.GetDirection&#40;&#41; &#41;;
if &#40;d != 0&#41;
&#123;
float dist = -&#40;DOT&#40; m_Plane.N, a_Ray.GetOrigin&#40;&#41; &#41; + m_Plane.D&#41; / d;
if &#40;dist > 0&#41;
&#123;
if &#40;dist < a_Dist&#41;
&#123;
a_Dist = dist;
return HIT;
&#125;
&#125;
&#125;
return MISS;
&#125;
Will let know how it works out.

And finally, of course the error can be in my test app, but i refuse to believe it :P

User137
17-03-2008, 09:07 AM
Noo, there is now variable in both functions that i don't understand. P1.PP in first code and m_Plane.D in C++ source. They must be some single floating point value from plane but what is it?

Maybe means this from plane equation:

function GetD(const p,n: P3f): single;
begin
//result:=-(a*x+b*y+c*z);
result:=-(n.x*p.x + n.y*p.y + n.z*p.z);
end;

Getting some results finally with the C++ function :) Not all perfect yet but now there is adjustments to do on test application...

Final edit: It all works now, all versions of now 3 different ones :D Of course it was faulty test application to cause. Properly rendered plane again after drawing its normal that didn't appear where it should've been...

Thanks all.

This being the fastest:
function rayPlaneCollision(const rayOrigin, rayDirection,
planeOrigin, planeNormal: P3f; intersection: pP3f): single;
var d,numer,denom: single;
begin
result:=-1;
d:=planeNormal.x*planeOrigin.x+planeNormal.y*plane Origin.y+planeNormal.z*planeOrigin.z;
numer:=-planeNormal.x*rayOrigin.x-planeNormal.y*rayOrigin.y-planeNormal.z*rayOrigin.z+d;
denom:=planeNormal.x*rayDirection.x+planeNormal.y* rayDirection.y+planeNormal.z*rayDirection.z;
if denom=0 then exit;
result:=numer/denom;
if intersection<>nil then begin
intersection.x:=rayOrigin.x+result*rayDirection.x;
intersection.y:=rayOrigin.y+result*rayDirection.y;
intersection.z:=rayOrigin.z+result*rayDirection.z;
end;

cronodragon
17-03-2008, 03:08 PM
Noo, there is now variable in both functions that i don't understand. P1.PP in first code and m_Plane.D in C++ source. They must be some single floating point value from plane but what is it?

PP stands for "producto punto" in spanish, that means "dot product". For your C++ code the D must stand for dot product. The plane is defined with just the Normal of the plane, which is a 3D vector, and the dot product.

The author is using integer arithmetic, that's why he uses those Int64.

Anyway, it's great it works for you! :D