Copyright (c) 2009, Jesus Mena
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Jamie Heather (view profile)
Nice vectorised implementation with handy options for choosing between ray/line/segments and handling numerical precision issues - thanks!
mirshad (view profile)
the code is ok...
i have a mistake.thank u
mirshad (view profile)
magnitude==andaze
i changed it formerly.
mirshad (view profile)
tnx.this code, give me a incorrect response for this inputs;
LINE = [0 0 1 1 1 2 ];
tri = [0 0 5;5 0 0;0 5 0;0 0 5];
magnitude=sqrt((LINE(4)-LINE(1))^2+(LINE(5)-LINE(2))^2+(LINE(6)-LINE(3))^2);
direction=[LINE(4)-LINE(1)/andaze LINE(5)-LINE(2)/andaze LINE(6)-LINE(3)/andaze];
direction=[LINE(4)-LINE(1)/andaze LINE(5)-LINE(2)/andaze LINE(6)-LINE(3)/andaze];
rayTriangleIntersection(LINE(1:3),direction,tri(1,:),tri(2,:),tri(3,:))
Kapil Dev (view profile)
Thank you very much. Really very useful piece of code.
dberm22 (view profile)
Works great, but too slow for what I needed. Mostly not your fault, but I needed to call this function more than a million times, and the built-in matlab routines you called were weighing it down.
I took your code and sped it up. for my test with 90,000 calls, it takes your function 16s, and mine only takes 1.5. That's a speedup of 10x! Here's what I have:
function [flag, u, v, t] = fastRayTriangleIntersection (o, d, p0, p1, p2)
% Ray/triangle intersection using the algorithm proposed by Möller and Trumbore (1997).
%
% Input:
% o : origin.
% d : direction.
% p0, p1, p2: vertices of the triangle.
% Output:
% flag: (0) Reject, (1) Intersect.
% u,v: barycentric coordinates.
% t: distance from the ray origin.
% Author:
% Originally written by Jesus Mena, edited by David Berman (dberm22@gmail.com)
epsilon = 0.00001;
e1 = p1-p0;
e2 = p2-p0;
q = [d(2)*e2(3)-d(3)*e2(2), d(3)*e2(1)-d(1)*e2(3), d(1)*e2(2)-d(2)*e2(1)]; %cross product
a = e1(1)*q(1)+e1(2)*q(2)+e1(3)*q(3); % determinant of the matrix M
if (a>-epsilon && a<epsilon)
% the vector is parallel to the plane (the intersection is at infinity)
flag=0;
u=0;
v=0;
t=0;
return;
end;
f = 1/a;
s = o-p0;
u = f*(s(1)*q(1)+s(2)*q(2)+s(3)*q(3));
if (u<0.0)
% the intersection is outside of the triangle
flag=0;
u=0;
v=0;
t=0;
return;
end;
r = [s(2)*e1(3)-s(3)*e1(2), s(3)*e1(1)-s(1)*e1(3), s(1)*e1(2)-s(2)*e1(1)];
v = f*(d(1)*r(1)+d(2)*r(2)+d(3)*r(3));
if (v<0.0 || u+v>1.0)
% the intersection is outside of the triangle
flag=0;
u=0;
v=0;
t=0;
return;
end;
t = f*(e2(1)*r(1)+e2(2)*r(2)+e2(3)*r(3)); % verified!
flag = 1;
return;
end
Thanks for the basis!
dberm22 (view profile)
Works great, but too slow for what I needed. Mostly not your fault, but I needed to call this function more than a million times, and the built-in matlab routines you called were weighing it down.
I took your code and sped it up. for my test with 90,000 calls, it takes your function 16s, and mine only takes 2. That's a speedup of 8x! Here's what I have:
function [flag, u, v, t] = fastRayTriangleIntersection (o, d, p0, p1, p2)
% Ray/triangle intersection using the algorithm proposed by Möller and Trumbore (1997).
%
% Input:
% IMPORTANT: ALL INPUTS NEED TO BE 3 ELEMENT ROW VECTORS
% o : origin.
% d : direction.
% p0, p1, p2: vertices of the triangle.
% Output:
% flag: (0) Reject, (1) Intersect.
% u,v: barycentric coordinates.
% t: distance from the ray origin.
% Author:
% Originally written by Jesus Mena, edited by David Berman (dberm22@gmail.com)
epsilon = 0.00001;
e1 = p1-p0;
e2 = p2-p0;
q = [d(2)*e2(3)-d(3)*e2(2), d(3)*e2(1)-d(1)*e2(3), d(1)*e2(2)-d(2)*e2(1)]; %cross product
a = e1*q'; % determinant of the matrix M
if (a>-epsilon && a<epsilon)
% the vector is parallel to the plane (the intersection is at infinity)
flag=0;
u=0;
v=0;
t=0;
return;
end;
f = 1/a;
s = o-p0;
u = f*s*q';
if (u<0.0)
% the intersection is outside of the triangle
flag=0;
u=0;
v=0;
t=0;
return;
end;
r = [s(2)*e1(3)-s(3)*e1(2), s(3)*e1(1)-s(1)*e1(3), s(1)*e1(2)-s(2)*e1(1)];
v = f*d*r';
if (v<0.0 || u+v>1.0)
% the intersection is outside of the triangle
flag=0;
u=0;
v=0;
t=0;
return;
end;
t = f*e2*r'; % verified!
flag = 1;
return;
end
Only caveat is that all inputs need to be 3 element ROW vectors, and there is absolutely no error checking. If you can live with that, then this will do it!
Thanks for the basis!
Alexander (view profile)
Thank you! I found this to be very useful! Thanks also to Wouter for the suggestion of using "sum" and an own "cross" function. I can't confirm factor 10, but factor 3 isn't so bad as well!
zf (view profile)
HARI KRISHNAN (view profile)
Alec (view profile)
Sunghun Jung (view profile)
Thank you vey much for this useful file.
Btw, could anyone give me a tip how to calculate the direction vector if I have two 3 dimensional nodes?
Wouter (view profile)
Excellent, it helped me a lot.
You can speed this function up considerably by not using the Matlab functions cross and dot, these are unnecessarily slow. Make your own cross product function, and use "sum(a.*b)" instead of "dot(a,b), and it'll will be at least a factor 10 faster.
Christian (view profile)