Sign in to follow this  
Endar

Ray and plane intersection

Recommended Posts

Endar    668
I've decided that it's pretty much time to write this code which I have been putting off for a while: AABB code for ray intersection. I was looking at Woo's method for finding the intersection between a ray and an AABB, so, I'll need a plane object and some functions to find the intersection between a ray and a plane. Here's what I've got so far:
/**
 * A class to represent a plane.
 */
class CPlane
{
private:

	util::CVector3d m_normal;		///< The normal of the plane
	float m_shiftConstant;		///< The plane-shift constant, the amount that the plane
					///< is shifted along the axes. Is calculated by the dot product
					///< between the normal and a point on the plane.

public:

	/**
	 * Constructor
	 * \param n A const reference to a CVector3d object that is the normal of the plane
	 * \param s The plane-shift constant (the dot product with the normal and a point on the plane
	 */
	CPlane( const util::CVector3d& n, float s )
		:	m_normal(n), m_shiftConstant(s)
	{
		m_normal.normalize();	// make sure that the vector is normalized
	}

	/**
	 * Constructor
	 * \param n A const reference to a CVector3d object that is the normal of the plane
	 * \param p A point on the plane
	 */
	CPlane( const util::CVector3d& n, const util::CVector3d& p )
		:	m_normal(n)
	{
		m_normal.normalize();		// make sure that the vector is normalized
		m_shiftConstant = n.dotProduct(p);	// calc the plane-shift constant
	}

	/**
	 * Calculate if a point is on the plane.
	 * \param p A const reference to a CVector3d object that represents a point
	 * \return A boolean value that represents whether the passed point is on the plane
	 */
	bool pointOnPlane(const util::CVector3d& p)
	{
		/*
			Because the plane shift constant is the dot product between the normal of the plane,
			to find out if the point is on the plane all we need to do is calc the dot product
			of the normal and the passed point and see if that equals the shift constant we have
		 */
		return ( m_normal.dotProduct(p) == m_shiftConstant );
	}


	/**
	 * Find if the passed ray intersects with this plane.
	 * \param o The origin point of the vector.
	 * \param d The direction of the vector (normalized)
	 * \return 'true' if the plane and ray intersect, else 'false'
	 */
	bool intersect(const util::CVector3d& o, const util::CVector3d& d)
	{
		util::CVector3d dir(d);

		if( dir.getLengthSquared() != 1.0f )	// who wants to do an unecessary sqrt?
			dir.normalize();					// make sure the vector is normalized

		// if direction of the ray is perpendicular to the plane normal, then the plane
		// and the direction are parallel, so there won't even be an intersection,
		// or, the ray is part of the plane (completely on it)
		if( m_normal.dotProduct(dir) == 0.0f )
			return false;		// no intersection
	}

	/**
	 * Return the point of intersection between the ray and plane.
	 * \param o The origin point of the vector.
	 * \param d The direction of the vector (normalized)
	 * \return The 3d coords of where the intersection point on the plane if there is an intersection.
	 * Else it just returns (0,0,0)
	 */
	util::CVector3d intersectionPoint(const util::CVector3d& o, const util::CVector3d& d)
	{
		// check to see if there is an intersection point
		if( !intersect(o,d) )
			return CVector3d(0,0,0);

		/* calculate the intersection point */

		util::CVector3d dir(d);
		
		if( dir.getLengthSquared() != 1.0f )	// who wants to do an unecessary sqrt?
			dir.normalize();		// make sure the vector is normalized


		// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rayplane_intersection.htm
		float t = (-(m_normal.dotProduct(o) + m_shiftConstant)) / m_normal.dotProduct(dir);

		// calcuate the intersection point
		util::CVector3d i;
		i.x = m_normal.x + (dir.x * t);
		i.y = m_normal.y + (dir.y * t);
		i.z = m_normal.z + (dir.z * t);

		return i;	// return the point of intersection
	}

};




Can anyone see anything wrong with it? (I haven't found anything, but, better to be safe)

Share this post


Link to post
Share on other sites
jyk    2094
Here's your code with some observations/suggestions added:


class CPlane
{
private:

util::CVector3d m_normal;
// 'constant' or 'distance' are more common names for this field,
// and might be easier to read and type
float m_shiftConstant;

public:

CPlane( const util::CVector3d& n, float s )
: m_normal(n), m_shiftConstant(s)
{
m_normal.normalize();

// If you're going to normalize, you have to normalize the distance also:

// float invLength = 1.0f / m_normal.length();
// m_normal *= invLength;
// m_shiftConstant *= invLength;
}

CPlane( const util::CVector3d& n, const util::CVector3d& p )
: m_normal(n)
{
m_normal.normalize();
m_shiftConstant = n.dotProduct(p);
}

bool pointOnPlane(const util::CVector3d& p)
{
// This will hardly ever, if ever, return true due to floating-
// point roundoff. For this reason, point-on-plane or point-plane
// classification tests are usually performed with an epsilon.
return ( m_normal.dotProduct(p) == m_shiftConstant );
}

bool intersect(const util::CVector3d& o, const util::CVector3d& d)
{
util::CVector3d dir(d);

// Same problem here - the squared length will almost never be exactly
// 1, even if the vector is normalized. Also, for this particular test
// the line direction probably doesn't even need to be normalized.
if( dir.getLengthSquared() != 1.0f )
dir.normalize();

// Same problem - dot product will almost never be exactly 0.
if( m_normal.dotProduct(dir) == 0.0f )
return false;
}

util::CVector3d intersectionPoint(const util::CVector3d& o, const util::CVector3d& d)
{
// This probably isn't the best way to indicate no intersection. For one thing,
// what if there is an intersection and it just happens to be (0,0,0)? Unlikely,
// but not impossible. In any case, better to make the function return a boolean
// result, and make the intersection point a non-const reference argument.
if( !intersect(o,d) )
return CVector3d(0,0,0);

util::CVector3d dir(d);

// Same problem with the exact float compare, and again normalization is
// not required.
if( dir.getLengthSquared() != 1.0f )
dir.normalize();

float t = (-(m_normal.dotProduct(o) + m_shiftConstant)) / m_normal.dotProduct(dir);

// This is picky, but since 'i' is usually used as an index it can be
// confusing to see it used for other types of variables. 'intersection'
// might be better.
util::CVector3d i;

// Here's a place where overloaded operators for the vector class would be nice...
i.x = m_normal.x + (dir.x * t);
i.y = m_normal.y + (dir.y * t);
i.z = m_normal.z + (dir.z * t);

return i;
}

};



I'll also add that ray-AABB intersection code is usually written to take advantage of the parallel axis-aligned planes. But, the concepts in the above code are still applicable.

Share this post


Link to post
Share on other sites
Endar    668

bool pointOnPlane(const util::CVector3d& p)
{
// This will hardly ever, if ever, return true due to floating-
// point roundoff. For this reason, point-on-plane or point-plane
// classification tests are usually performed with an epsilon.
return ( m_normal.dotProduct(p) == m_shiftConstant );
}


<s>First of all, I have to read up on point on plane tests,</s> and second, what size epsilon would take into account a small amount of error, but still be small enough not to invalidate the test if the vector was wrong?

Edit::
So the actual test I'm doing is fine, except for the fact that I don't acount for float rounding-errors?

Edit2::

CPlane( const util::CVector3d& n, float s )
: m_normal(n), m_shiftConstant(s)
{
m_normal.normalize();

// If you're going to normalize, you have to normalize the distance also:

// float invLength = 1.0f / m_normal.length();
// m_normal *= invLength;
// m_shiftConstant *= invLength;
}


Normalizing the distance? What is that, and why should I normalize the distance as well? The normal is normalized because it is a direction, but why the distance? Wouldn't that completely change the position of the plane?

[Edited by - Endar on November 24, 2005 3:02:29 AM]

Share this post


Link to post
Share on other sites
Endar    668
Okay, here it is again with some changes, most of the ones that you suggested jyk. And, I'm always stuck when its possible that I'll need to return an error message when any value that could be returned could be valid. I did change this beacuse it was really the only way, but it seems a little inelegant. But hey, maybe it's just me. [smile]


#define PLANE_EPSILON 0.0001f



/**
* A class to represent a plane.
*/

class CPlane
{
private:

util::CVector3d m_normal; ///< The normal of the plane
float m_distance; ///< The plane-shift constant, the distance of the plane
///< from the origin. Is calculated by the dot product
///< between the normal and a point on the plane.

public:

/**
* Constructor
* \param n A const reference to a CVector3d object that is the normal of the plane
* \param s The plane-shift constant (the dot product with the normal and a point on the plane
*/

CPlane( const util::CVector3d& n, float s )
: m_normal(n), m_distance(s)
{
m_normal.normalize(); // make sure that the vector is normalized
}

/**
* Constructor
* \param n A const reference to a CVector3d object that is the normal of the plane
* \param p A point on the plane
*/

CPlane( const util::CVector3d& n, const util::CVector3d& p )
: m_normal(n)
{
m_normal.normalize(); // make sure that the vector is normalized
m_distance = n.dotProduct(p); // calc the plane-shift constant
}

/**
* Calculate if a point is on the plane.
* \param p A const reference to a CVector3d object that represents a point
* \return A boolean value that represents whether the passed point is on the plane
*/

bool pointOnPlane(const util::CVector3d& p)
{
/*
Because the plane shift constant is the dot product between the normal of the plane,
to find out if the point is on the plane all we need to do is calc the dot product
of the normal and the passed point and see if that equals the shift constant we have
*/

float dot = m_normal.dotProduct(p);
return ( dot >= m_distance - PLANE_EPSILON || dot <= m_distance + PLANE_EPSILON );
}


/**
* Find if the passed ray intersects with this plane.
* \param o The origin point of the vector.
* \param d The direction of the vector (normalized)
* \return 'true' if the plane and ray intersect, else 'false'
*/

bool intersect(const util::CVector3d& o, const util::CVector3d& d)
{
util::CVector3d dir(d);

float len = dir.getLengthSquared(); // who wants to do an unecessary sqrt?
if( !(len >= 1.0f - PLANE_EPSILON || len <= 1.0f + PLANE_EPSILON) )
dir.normalize(); // make sure the vector is normalized

// if direction of the ray is perpendicular to the plane normal, then the plane
// and the direction are parallel, so there won't even be an intersection,
// or, the ray is part of the plane (completely on it)
float dot = m_normal.dotProduct(dir);
if( !(dot >= 0.0f - PLANE_EPSILON || dot <= 0.0f + PLANE_EPSILON) )
return false; // no intersection

return true; // there is an intersection
}

/**
* Return the point of intersection between the ray and plane.
* \param o The origin point of the vector.
* \param d The direction of the vector (normalized)
* \param intersection A reference to a vector in which the
* intersection point, if there is one, will be stored.
* \return 'true' if there is an intersection point, else 'false'
*/

bool intersectionPoint(const util::CVector3d& o, const util::CVector3d& d, util::CVector3d& intersection)
{
// check to see if there is an intersection point
if( !intersect(o,d) )
return false;

/* calculate the intersection point */

util::CVector3d dir(d);

float len = dir.getLengthSquared(); // who wants to do an unecessary sqrt?
if( !(len >= 1.0f - PLANE_EPSILON || len <= 1.0f + PLANE_EPSILON) )
dir.normalize(); // make sure the vector is normalized


// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rayplane_intersection.htm
float t = (-(m_normal.dotProduct(o) + m_distance)) / m_normal.dotProduct(dir);

// calcuate the intersection point
intersection = m_normal + (dir * t);

return true;
}

};



The epsilon value probably isn't small enough, and, would it be okay to use the same epsilon value all over? Or would I have to use a one when calculating the length and a different one when calculating the dot product?

Share this post


Link to post
Share on other sites
timw    598
your epslon is small enough. another suggestion I'd use fabs() to avoid the two checks as that is simply a bit set and should be faster.

Tim

Share this post


Link to post
Share on other sites
timw    598
if( !(dot >= 0.0f - PLANE_EPSILON || dot <= 0.0f + PLANE_EPSILON) )
return false; // no intersection
//to
if( fabs(dot) > PLANE_EPSILON)

also it's probably better to just assume it's normalized rather then make the check. I did that all over the place in my ray tracer for fear of it, but realized later I just had redundant code everywhere. normalize it durring it's creation and you know it's ok, as long as you dont modify it. that's my opinion anyway.

on second thought, maybe fabs() wont speed it up lol. it seems it's probably inlined as a and 0x7fffffff.

Share this post


Link to post
Share on other sites
Endar    668
Oh, okay. But that just for the dot product, right? I mean, you can't do something like that when you're checking the length can you? Because in that case it's 1.0f and not 0.0f.

Share this post


Link to post
Share on other sites
timw    598
Quote:

you can't do something like that when you're checking the length can you?

na you couldn't. I'm thinking fabs isn't the best hack lol. but definitly I'd say, dont check for unit length all the time, it really doesnt matter. just enforce the caller to provide unit length. even if it's not unit length, the calculation will still be correct(unless its crazy large magnitude, in which case it might effect precision) but that's just a dumb caller. I don't think it's woruth the 3 multiplies and two additions, especailly cus stuff like this gets called so much.

happy holiday!

Tim

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this