Archived

This topic is now archived and is closed to further replies.

davidofqls

Calculating ray-triangle intercection

Recommended Posts

davidofqls    122
I''m looking for an algorithm to find the intersection of a ray with a triangle, that will tell at what point they intersect if they do. Does anyone know of such an algorithm? David

Share this post


Link to post
Share on other sites
Eric Lengyel    3461
You want to look up "barycentric coordinates".

For a thorough discussion of the ray-triangle intersection problem, see Section 5.2.1 of Mathematics for 3D Game Programming and Computer Graphics, Section Edition:

http://www.terathon.com/books/mathgames2.html

[edited by - Eric Lengyel on January 23, 2004 1:23:00 AM]

Share this post


Link to post
Share on other sites
digisoap    241
Stop spamming the boards w/ a link to your book. If you know the answer - which I would hope you do, since you wrote a book about it - then provide the answer. People are not going to pay $50.00+ if they need an answer to a single question.

quote:
Original post by Eric Lengyel
You want to look up "barycentric coordinates".

For a thorough discussion of the ray-triangle intersection problem, see Section 5.2.1 of Mathematics for 3D Game Programming and Computer Graphics, Section Edition :

http://www.terathon.com/books/mathgames2.html

[edited by - Eric Lengyel on January 23, 2004 1:23:00 AM]



-Nick "digisoap" Robalik
Web & Print Design, 2D & 3D Illustration and Animation, Game Design
http://www.digital-soapbox.com
nick@digital-soabox.com

Share this post


Link to post
Share on other sites
I have to agree with the mod. If you were to actually help people, and just have a link to your book in your sig, eventually people would notice that this really nice, helpful, knowlegable guy happens to have a book out. I''d be more inclined to look at a book of someone who''s made themselves respected in the forums.

Given the number of free triangle-ray intersection algorithms on the net, just pointing someone to your book is extremely cheesy, spammy, and makes you out as self-serving "anything to make a buck" kinda guy rather than helpful.

For the OP: Try this link.

Share this post


Link to post
Share on other sites
Rocket05    152
Simple:

Given the two points of the line, L1 and L2, check if the line intersects the plane of the triangle

let N be the normal of the triangle, and * represent a dot product, D be the distance part of the plane equation, in the form of Ax + By + Cz = D


float d1,d2;
d1 = N*L1 - D;
d2 = N*L2 - D;

if((d1 < 0.0f && d2 >= 0.0f) || (d1 >= 0.0f && d2 < 0.0f))
{
// the line intersects the plane
// check how much, and get an intersection point
Vector3 linedirection = (L2-L1);
float linelength = linedirection.length();
float percent = d1/linelength;

// percent now contains a value between 0 and 1
Vector3 intersection = L1 + (linedirection*percent);
// we now have an intersection point, now see if its in
// the triangle
// remember that for any point inside a convex polygon,
// the angle between the point and all verticies will always
// have the sum of 360 degrees (2PI for radians).

// assuming the points of the triangle are labled v1,v2,v3
float angle = 0.0f;
Vector3 direction1 = v1 - intersection;
Vector3 direction2 = v2 - intersection;
Vector3 direction3 = v3 - intersection;
direction1.normalize();
direction2.normalize();
direction3.normalize();

angle += direction1 * direction2;
angle += direction2 * direction3;
angle += direction3 * direction1;

angle -= 2*PI;
if(fabsf(angle) < (EPSILON))
{
// line intersects the triangle!
}
}


i hate to give code w/o an explanation, but i had to write this quick off the top of my head. if you need more explanation, just say so and i''ll try to explain it better

and to others, yes i know there is probably a better/faster/more optimized way to do this, but this was off the top of my head. if you have a better version, please post it.

Share this post


Link to post
Share on other sites
You give the code, I'll explain...

The Plane (Triangle explained later)

Given a plane in 3D space and a point that lies on it, there is exactly one vector that's perpendicular to it. Thus, a plane can be described with a point and a vector. The perpendicular vector is given the name "normal." Because the normal must be perpendicular to the plane, the dot product btw a vector on the plane and the normal must be 0, because the dot product of two vectors has the identity:

A dot B = |A||B|cos(t)
where t is the angle btw them, and |A| is the magnitude of A.
As you can see, the COS of a 90 degree (or pi/2 radian) angle is 0.

Knowing this, one can build the plane equation as the dot of the normal N = (a, b, c), a point know to be on the plane p0 = (x0, y0, z0), and another point on the plane p1 = (x1, y1, z1). For p1 and p0 to be on the plane defined by N, their difference, which returns a vector from one to the other, must dot with N to be 0, thus

N dot (p1-p0) = 0

The vector from the subtraction points to p1 from p0, but wheter the direction is parallel or antiparallel to this is trivial. Now for discussion the dot operator will be changed to the * operator, but on scalars it'll mean multiplication. The above equation is now expanded to:

a(x1-x0) + b(y1-y0) + c(z1-z0) = 0

and the properties of the dot product say that it's distributive, so

a*x1 + b*y1 + c*z1 - a*x0 - b*y0 - c*z0 = 0

To form the infamous plane equation we call -N*p0 = d, and make p1 = p. Now

a*x + b*y + c*z + d = 0

This is the equation for a plane. Inserting a points x, y, and z components will make this equation true if the point is on a plane. Or, better yet, if a point is on a plane, inserting the point and representing the equation in vector notation, this must be true:

N*p + d = 0

The Line

This is the way I learned it, so we're going to represent the ray as a parametric line. Here is the equation for one:

p = (x,y,z)
p(t) = p0 + v*t (i'll explain)

Here, p is is a point defined by the component-wise parametric functions x(t), y(t), and z(t), but because we are smart, we can do the problem in vector form. In the above equation, p0 is the starting point (and has nothing to do with the plane derivation's p0), and v is a vector that is conincident with the line. When t varies and v*t is added to p0, a point in the diretion of v, offset by p0, will be returned, pending t is not 0 (leaving the point to be just p0, no direction I mean). So, how would we define this v and p0? Well, we could do it with two points. This comes in handy in 3D graphics, b/c geometric info (ie a triangle) is defined by points, and a 2 points can define a line. The points are thusly called p0, and p1. Now how do we build v, you ask? Well, you must ask yourself another question: how do you want to vary t? Well, I'll answer it for you. You can vary t any way you'd like, what I mean is when p(t) goes from p0 to p1, would you like t to be 1? |v|? or what? Well, it doesn't matter for us, so I say 1. That means, p(1) has to equal p1, and t varies from [0...1] as p varies from [p0...p1], get it? Well...

p(1) = p1 = p0 + v*1
v = p1 - p0

Now v points toward p1 and starts at p0, and has a magnitude of the distance btw p1 and p0. Now we have the line. The points making up the line can define the edge of a triange, say, and the line will be the triangle's edge.

The Intersection

You already have an equation for a point on a line (p(t)), and the plane equation, so now all you have to do is find for what value of t, p(t) satisfies the plane equation. Can you say substitution?

N*p(t) + d = 0
N*(p0 + v*t) + d = 0
N*p0 + N*v*t + d = 0

t = (-N*p0 - d)/(N*v)

Inserting t back into the p(t) will give you the point!

p(t) = p0 - v*(N*p0 + d)/(N*v)

Therefore, given two points that define a line, and a plane, you can find the intersection point on the line, if there is one.

The Trinagle

A triange lies on only one plane right? So how would we find that plane? We'll use the cross product. The cross product returns a vector perpendicular to the two vectors it's done one, thus forming an orthogonal basis. In our case, this vector will be our normal, N. Now, the vectors on the plane that we're crossing are the v1-v0 and v2-v0, made from the triangle's three vertices. Since we know all three have to be on this plane, we have to dot N with one of them to get the d component of the plane. Thus the plane equation components can be formed like so

N = (v1-v0) cross (v2-v0)
d = N dot v0

I chose to dot it with v0 because it seems logical

My own implementation with overloaded operators (ie the intersection can be coded in mathmatical notation):

Vector.h

#pragma once

//#include "Matrix.cpp"


namespace Vector
{
//using namespace Matrix;


template<class T>
class Vector3
{
public:
//default constructor, does nothing

Vector3();
//init w/ 3 floats e.g.

Vector3(T, T, T);
//copy constructor

Vector3(const Vector3 &);
//indexing operators

T& operator[](unsigned int i);
const T& operator[](unsigned int i) const;
//assignment

Vector3& operator=(const Vector3&);
//v plus-equals u

Vector3& operator+=(const Vector3&);
//v minus-equals u

Vector3& operator-=(const Vector3&);
//v cross-equals u

Vector3& operator^=(const Vector3&);
//v times-equals matrix

//Vector3& operator*=(const Matrix::Matrix3<T>&);

//v times-equals scalar

Vector3& operator*=(T);
//v over-equlas scalar

Vector3& operator/=(T);
//v plus u

Vector3 operator+(const Vector3&) const;
//v minus u

Vector3 operator-(const Vector3&) const;
//v dot u

T operator*(const Vector3&) const;
//v cross u

Vector3 operator^(const Vector3&) const;
//v times matrix

//Vector3 operator*(const Matrix3<T>&) const;

//v over scalar

Vector3 operator/(T) const;
//negated v

Vector3 operator-() const;
//length of v

T magnitude() const;
//normalised v, |v| = 1

Vector3 getNormalise() const;
void normalise();
//v times scalar

Vector3 operator*(T s) const;
//scalar times v

friend Vector3 operator*(T s, const Vector3 &v)
{
return Vector3(s*v.x, s*v.y, s*v.z);
}
//matrix times this

/*friend Vector3 operator*(const Matrix3<T> &m, const Vector3 &v)
{
return Vector3<T>(
v[0]*m[0] + v[1]*m[1] + v[2]*m[2],
v[0]*m[3] + v[1]*m[4] + v[2]*m[5],
v[0]*m[6] + v[1]*m[7] + v[2]*m[8]);
}*/

private:
T x, y, z;
};
}//namespace Vector



Vector.cpp

#pragma once

#include <assert.h>
//#include "Math.h"

#include "Vector.h"

namespace Math
{
template<class T>
inline T Sqrt(T x)
{
return (T)sqrt((float)x);
}
}

namespace Vector
{
//using namespace Matrix;


template<class T>
Vector3<T>::Vector3()
{
}
template<class T>
Vector3<T>::Vector3(T X, T Y, T Z)
: x(X), y(Y), z(Z)
{
}
template<class T>
Vector3<T>::Vector3(const Vector3 &v)
: x(v.x), y(v.y), z(v.z)
{
}
template<class T>
inline T& Vector3<T>::operator[](unsigned int i)
{
assert(i < 3);
return *(&x+i);
}
template<class T>
inline const T& Vector3<T>::operator[](unsigned int i) const
{
assert(i < 3);
return *(&x+i);
}
template<class T>
inline Vector3<T>& Vector3<T>::operator=(const Vector3 &v)
{
x = v.x;
y = v.y;
z = v.z;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator+=(const Vector3 &v)
{
(*this) = (*this)+v;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator-=(const Vector3 &v)
{
(*this) = (*this)-v;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator^=(const Vector3 &v)
{
(*this) = (*this)^v;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator*=(const Matrix3<T> &m)
{
(*this) = (*this)*m;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator*=(T s)
{
(this*) = (this*)*s;
return *this;
}
template<class T>
inline Vector3<T>& Vector3<T>::operator/=(T s)
{
if(s != 0)
s = 1/s;
(*this) *= s;
return *this;
}
template<class T>
inline Vector3<T> Vector3<T>::operator+(const Vector3 &v) const
{
return Vector3(x+v.x, y+v.y, z+v.z);
}
template<class T>
inline Vector3<T> Vector3<T>::operator-(const Vector3 &v) const
{
return Vector3(x-v.x, y-v.y, z-v.z);
}
template<class T>
inline T Vector3<T>::operator*(const Vector3 &v) const
{
return x*v.x + y*v.y + z*v.z;
}
template<class T>
inline Vector3<T> Vector3<T>::operator^(const Vector3 &v) const
{
return Vector3(
y*v.z - z*v.y,
z*v.x - x*v.z,
x*v.y - y*v.x);
}
/*template<class T>
inline Vector3<T> Vector3<T>::operator*(const Matrix3<T> &m) const
{
return Vector3(
x*m[0] + y*m[3] + z*m[6],
x*m[1] + y*m[4] + z*m[7],
x*m[2] + y*m[5] + z*m[8]);
}*/

template<class T>
inline Vector3<T> Vector3<T>::operator/(T s) const
{
if(s != 0)
s = 1/s;
return Vector3(x*s, y*s, z*s);
}
template<class T>
inline Vector3<T> Vector3<T>::operator-() const
{
return Vector3(-x, -y, -z);
}
template<class T>
inline T Vector3<T>::magnitude() const
{
return Math::Sqrt<T>(x*x + y*y + z*z);
}
template<class T>
inline Vector3<T> Vector3<T>::getNormalise() const
{
T m = magnitude();
return (*this)/m;
}
template<class T>
inline void Vector3<T>::normalise()
{
T m = magnitude();
(*this) /= m;
}
template<class T>
inline Vector3<T> Vector3<T>::operator*(T s) const
{
return Vector3(x*s, y*s, z*s);
}
}//namespace Vector



Plane.h

#pragma once

#include "Vector.cpp"

namespace Plane
{
using namespace Vector;

template<class T>
class Plane3
{
public:
Plane3();
Plane3(T, T, T, T);
Plane3(const Vector3<T>&, const Vector3<T>&, const Vector3<T>&);
Plane3(const Plane3&);
Vector3<T> lineIntersect(const Vector3<T>&, const Vector3<T>&);
public:
Vector3<T> n;
T d;
};
}//namespace Plane



Plane.cpp

#pragma once

#include "Plane.h"

namespace Plane
{
using namespace Vector;

template<class T>
Plane3<T>::Plane3()
{
}
template<class T>
Plane3<T>::Plane3(T x, T y, T z, T D)
: n(Vector3<T>(x, y, z)), d(D)
{
}
template<class T>
Plane3<T>::Plane3(const Vector3<T> &p0, const Vector3<T> &p1,
const Vector3<T> &p2)
: n((p0-p2)^(p1-p2)), d(-n*p2)
{
}
template<class T>
Plane3<T>::Plane3(const Plane3& p)
: n(p.n), d(p.d)
{
}
template<class T>
inline Vector3<T> Plane3<T>::lineIntersect(const Vector3<T> &p0,
const Vector3<T> &p1)
{
Vector3<T> v = p1-p0;
Vector3<T> p = p0-v*(n*p0+d)/(n*v);
return p;
}
}//namespace Plane



Main.cpp

#include <stdio.h>
#include "Vector.cpp"
//#include "Matrix.cpp"

#include "Plane.cpp"
//#include "Math.h"

//#include "List.cpp"


using namespace Vector;
using namespace Matrix;
using namespace Plane;

void PrintVec3f(Vector3<float>&);
void PrintMat3f(Matrix3<float>&);



void main(void)
{

Vector3<float> x(1.0f, 0.f, 0.0f), y(0.0f, 1.0f, 0.0f), a(1.0f, 1.0f, 0.0f);
Plane3<float> xyPlane(x, y, a);
Vector3<float> p0(0.0f, 0.0f, 10.0f), p1(0.0f, 0.0f, -100.0f);
Vector3<float> p2 = xyPlane.lineIntersect(p0, p1);
PrintVec3f(p2);
}

void PrintVec3f(Vector3<float> &v)
{
printf("Vec3f = (%f, %f, %f) \n", v[0], v[1], v[2]);
}
void PrintMat3f(Matrix3<float> &m)
{
printf("Mat3f = \n");
printf("{%f, %f, %f} \n", m[0], m[1], m[2]);
printf("{%f, %f, %f} \n", m[3], m[4], m[5]);
printf("{%f, %f, %f} \n", m[6], m[7], m[8]);
}


--The point? You can do math-like expressions as shown above:

template
inline Vector3 Plane3::lineIntersect(const Vector3 &p0,
const Vector3 &p1)
{
Vector3 v = p1-p0;
Vector3 p = p0-v*(n*p0+d)/(n*v);
return p;
}


--Should compile in VC++ 6.0

EDIT: I meant cos t, not sin t in the dot identity



[edited by - temp_ie_cant_thinkof_name on March 15, 2004 2:09:00 PM]

Share this post


Link to post
Share on other sites
McDusty    122
Check out ''Fast minimum storage ray/triangle intersection'' - Mol97. (c++ code included)

http://www.graphics.cornell.edu/pubs/1997/MT97.html

This method gives you the barycentric coordinates and the distance along the ray that the triangle intersection occurs.

hitPoint = ray.origin + t * ray.direction

where t is the distance along the ray

Share this post


Link to post
Share on other sites