Sign in to follow this  
Endar

Woo's method

Recommended Posts

Endar    668
Taking a look at Woo's method from here for testing to see if a ray intersects with an AABB.
Quote:
From the six planes, reject those with back-facing normals From the three planes remaining, compute the t distance on the ray-plane test Select the farthest of the three distances Test if that point is inside the box's boundaries
And this is what woo's method consists of:
Quote:
  • Six dot products to check the first step
  • Three point-ray tests
  • A few comparisons for the last step
Now, I know that the dot product can show that two vectors are perpendicular, but remind me again, what does the dot product do in this case that will show that three of the planes are back-facing compared to the side of the aabb that the ray is facing?

Share this post


Link to post
Share on other sites
jyk    2094
As you note, the dot product of two perpendicular vectors is zero. If the dot product is not zero, the sign tells you whether the vectors are pointing generally in the same direction, or generally in the opposite direction (very generally). Drawing some pictures will make this clear. In the ray-AABB test you're describing, presumably the planes for which the dot product between the ray direction and the plane normal is positive are discarded, as the ray will intersect the corresponding parallel plane first (unless the ray is facing away from the box, in which case the value of t will be negative and the test will return false).

Share this post


Link to post
Share on other sites
Endar    668
So, if we consider 2d vectors and a unit circle, vectors that are in the same quadrant, going out from the origin, therefore going roughly in the same direction have positive dot products. Vectors that are in opposite quadrants have negative dot products. Right?


[Edited by - Endar on November 26, 2005 3:17:10 AM]

Share this post


Link to post
Share on other sites
haegarr    7372
The comparison with the "aligned" unit circle and quadrants isn't good. Look here:

Please consider the formula behind dot product:
v1 dot v2 = ||v1|| * ||v2|| * cos<v1,v2>

As you could see, there is the cosine of the angle between the both vectors. The cosine is +1 for 0 degree, and it is axis symmetrical. So the cosine is positive (and hence the entire dot product) as long as the angle between the two vectors is less than 90 degree. Due to the symmetry, if fixing v1, v2 could be in the range of (+90,-90) degree around v1! That makes a total of 180 degree, and hence the dot product is positive if v2 is in the same half sphere (or half circle in 2D) where v1 defines the one pole.

Similarly, the dot product will become negative if v2 is in the other half sphere (or half circle) where v1 defines the one pole.

EDIT: Look at the fine ASCII art of jyk below.

Share this post


Link to post
Share on other sites
jyk    2094
Here's a little diagram to illustrate:
 A__  |  __B
|\ | /|
\ | /
\ | /
\|/
------------>C
|
|
|
|
|
dot(A,C) is negative, and dot(B,C) is positive. In general, any vector that points to the left of the vertical line (or plane in 3d) has a negative dot product with C, and any vector that points to the right has a positive dot product.

Share this post


Link to post
Share on other sites
Endar    668
Presumably it's possible that eventually, by chance, a ray will be perpendicular to some faces. In that case, since we're picking 3, will it really matter which one I pick?

Eg. I have a ray in the direction of (0,0,-1), now, won't this be perpendicular to both (1,0,0) and (-1,0,0)? So, will it matter which I pick?

Share this post


Link to post
Share on other sites
haegarr    7372
If I understand you correctly, you ask which planes should remain after Woo's first step if the ray is also (partly) axis aligned, so that it is parallel to a couple of the planes of the BB?

If you mean that, and you pick any of the planes, then you have to compute the travel distance in the second step. But you could not due to the parallelism. So, without having it programmed yet, I assume you have no choice other than dropping all of those belonging planes. IMHO in such case you have no senseful possibility to pick 3 planes.

Quote:
Original post by Endar
Eg. I have a ray in the direction of (0,0,-1), now, won't this be perpendicular to both (1,0,0) and (-1,0,0)?

Definitely, yes.

Share this post


Link to post
Share on other sites
Endar    668
Quote:
Original post by haegarr
If you mean that, and you pick any of the planes, then you have to compute the travel distance in the second step. But you could not due to the parallelism. So, without having it programmed yet, I assume you have no choice other than dropping all of those belonging planes. IMHO in such case you have no senseful possibility to pick 3 planes.


Okay, I'm a little lost. Are you saying that I should forget about any planes that they ray is perpendicular to? So, in my little example, I would end up with only two planes to do the rest of the test against? The ones to which the ray is not perpendicular. Right?

Share this post


Link to post
Share on other sites
haegarr    7372
The problem will come with the second step, where you have to compute the hit point on the plane, correct?

If the ray is parallel to the plane, then 2 possible solutions may occur:

(1) The origin of the ray is outside the plane. In such case there is no hit point at all.

(2) The origin of the ray is inside the plane. In such case you have infinitely many hit points.

(EDIT: Try it: Compute the t parameter by the Ray-Plane test. You will encounter a division by zero due to the (n dot dir) term in the denominator.)

With both solutions you have computational problems since they simply makes no sense for the algo, as far as I see. So yes, I suggest to drop those planes, too, even if that means that you may have less than 3 planes.

BTW: The worst case would be 1 remaining plane.

Share this post


Link to post
Share on other sites
haegarr    7372
Appendix: The algo reduces the count of candidate planes to 1, and at last tests the hit on that plane against the volume. So I don't see a problem if the count is restricted formerly to 1 in the worst case. Even in the cases of parallel rays the concluding test seems me to work well.

Attention: You should not drop planes that are perpendicular to the ray, but planes that are _parallel_! The _normals_ of parallel planes are perpendicular to the ray! Please take care of this.

[Edited by - haegarr on November 26, 2005 4:26:29 AM]

Share this post


Link to post
Share on other sites
Endar    668
Quote:
Original post by haegarr
Attention: You should not drop planes that are perpendicular to the ray, but planes that are _parallel_! The _normals_ of parallel planes are perpendicular to the ray! Please take care of this.


Yeah, that's what I meant. I wrote planes, but was thinking about their normals.

Share this post


Link to post
Share on other sites
haegarr    7372
Quote:
Original post by Endar
Yeah, that's what I meant. I wrote planes, but was thinking about their normals.

Okay okay, I would be sure only you're not dropping the good planes ;-)

Share this post


Link to post
Share on other sites
Endar    668
Okay, check my common sense.

I have a aabb from (-1,-1,-1) to (2,1,1). I have a ray starting at (3,0,0) going in the direction (-1,0,0).

Now, correct me if I'm wrong, but, shouldn't the intersection point be (2,0,0), not "No intersection found"?

Share this post


Link to post
Share on other sites
haegarr    7372
I assume that the normals of the BB's planes are defined to point to the outside. (If they are defined to point inside the BB, then in the following the statements using the result of the dot product has to be negated.)

Then the first step of the algo defines the left plane (w/ the normal (-1,0,0) to be dropped since it is backfaced. You may detect this by the result of the dot product of the ray and the normal being positive.

Furthurmore the front and back planes w/ the normals (0,0,+/-1) as well as the top and bottom plane w/ normals (0,+/-1,0) are parallel to the ray. IMHO they have to be dropped.

The remaining plane is the right one w/ normal (1,0,0) that points against the given ray. So its dot product will be negative and hence the plane survives the first step of the algo.

Next, the second step of the algo could be done since at least one plane has still survived (IMHO at everytime at least one plane will be there). Computing the hit point will result in (2,0,0).

Proove: The ray equation is:
r(t) := (3,0,0)^T + t * (-1,0,0)^T
A valid plane equation of the remaining plane is (notice that infinitely many valid equations exist)
( p - (2,1,1)^T ) dot (1,0,0)^T = 0
so that the 4 plane parameters are computed to be (simply by elimiating the brackets above)
A := 1, B := 0, C := 0, D := -2
Putting into the equation for t (see the site you've cited in the orig. post)
t := - ( (1,0,0)^T dot (3,0,0)^T -2 ) / ( (1,0,0)^T dot (-1,0,0)^T ) = - (3-2) / (-1) = 1
Putting this into the equation for the ray yields in
r(t=1) = (3,0,0)^T + 1 * (-1,0,0)^T = (2,0,0)^T
q.e.d.

I don't see why "No intersection found" should come out. IMHO all works fine.

Share this post


Link to post
Share on other sites
Endar    668
It looks like my function that finds the intersection point returns (-2,0,0) instead of (2,0,0).

Edit:: From some preliminary (2) tests, it looks like the components should just be multiplied by -1 to get the correct answer. Should I just do this? Or do you think that this is indicative of a possibly bigger problem that would cause more normal cases to fail?

[Edited by - Endar on November 26, 2005 11:20:22 PM]

Share this post


Link to post
Share on other sites
haegarr    7372
You should never do a justification only to achieve the expected result. Please ask "why (or why not) does it come out". If you have a reasonable answer found, and the work-around fits into the answer, ok. But if you hide a bug, you don't know at what next opportunity it will have another drawback.

(Ok, I earn my money with programming, and so I may have another relation to that stuff. However, I write unit tests for my classes where ever possible. That costs time, of course, but gives (a) the safety that all goes the right way, and (b) allows quick prooves in cases where changes was made to the code. And (c), sometimes it let me think of cases that where not in my mind at the time of writing the code.)

You may take into account to post the belonging code for inspection by another pair of eyes.

EDIT:
Quote:
Original post by Endar
From some preliminary (2) tests, it looks like the components should just be multiplied by -1 to get the correct answer.

What components do you mean? Inverting the ray should not change anything, since the hit point will be the same.

If you use
r'(t) := -r(t) = (-3,0,0)^T + t * (1,0,0)^T
you get
t := - ( (1,0,0)^T dot (-3,0,0)^T -2 ) / ( (1,0,0)^T dot (1,0,0)^T ) = - (-3-2) / (1) = 5
and finally
r'(t=-5) = (-3,0,0)^T + 5 * (1,0,0)^T = (2,0,0)^T
This should work since it is just another representation of the same ray. However, it will fail in general, since negating the ray's origin will generally not yield in just another representation of the same ray.

Inverting the normal of the plane would be the other possibility, but then the formula of t as is given would no longer be correct. (However, the result would be (4,0,0) if I'm right, what still differs from (-2,0,0).)

[Edited by - haegarr on November 27, 2005 4:40:13 AM]

Share this post


Link to post
Share on other sites
Trienco    2555
Ok, this is just me being lazy, but the method as such should work fine for all bounding boxes. Seeing how we are using aabb all this can be made a lot simpler.

The normal for the "right" plane is always (1,0,0) (or -1,0,0 I usually like them pointing inward for bounding volumes). A dot product with this will ALWAYS return the direction vectors x.

Bottom line: scrap the dot products, all you want to know is if your ray is basically going left, right or doesn't change in x direction.

ie. you get six tests like this:
float t[3]={0}, *ptr=t;
if (RayDir.x<0) *(ptr++)=(bbMax.x-RayP.x)/RayDir.x;


meaning, if your ray is moving left, divice the distance between the start of the ray and the right plane of the box by the "distance" the ray travels along x per "unit" of movement to get the right t. Here the t's are stored in an array. They will be positive and not more than three.

float max=(t[0]>t[1]) ? (t[0]>t[2]?t[0]:t[2]) : (t[1]>t[2]?t[1]:t[2]);


Yes, ugly and I'm sure the stl has some max function for three operands (or at least vectors). Once you have the max t the rest is just getting the point of intersection and checking if it is between all the min and max x/y/z coords.

In other words, for AABB you shouldn't even start throwing the dot products and normals around. Save them for later when you want to do the same test for arbitrarily oriented boxes. That way it should be easier to visualize.

Or start with micro optimizing. I don't like the look of 7 if's and up to 3 / in a function that's potentially called a lot. If's make baby caches cry. On the other hand I would trust the compiler until some profiler tells me that this very function is making my app crawl. And then I would look for ways to call it less often ,-)

Share this post


Link to post
Share on other sites
haegarr    7372
Quote:
Original post by Trienco
The normal for the "right" plane is always (1,0,0) (or -1,0,0 I usually like them pointing inward for bounding volumes). A dot product with this will ALWAYS return the direction vectors x.

Bottom line: scrap the dot products, all you want to know is if your ray is basically going left, right or doesn't change in x direction.

Right you are: This is a suitable (and "should be done") optimization for AABB. I suggest to do such things after having the things working in general. IMHO the general solution introduces the best understanding of the principles.

Quote:
Original post by Trienco
ie. you get six tests like this:
float t[3]={0}, *ptr=t;
if (RayDir.x<0) *(ptr++)=(bbMax.x-RayP.x)/RayDir.x;

I don't see the need for 6 tests. Why to distinct between x>0 and x<0 in this example? That should not be necessary.

Share this post


Link to post
Share on other sites
Trienco    2555
Quote:
Original post by haegarr
I don't see the need for 6 tests. Why to distinct between x>0 and x<0 in this example? That should not be necessary.


For a second you had me wondering the same thing. Then I remembered that <0 uses a different plane and hence point than >0 and nestling the <0 >0 test inside yet another !=0 test didn't seem too useful.

However, it came to me, that checking for a line segment is much more likely than checking for an infinite ray. So it might be easier to do

float t, max=0;
if (rayDir.x<0) {
t=bbMax.x-rayP.x;
if (t>rayDir.x) return 0;
t /= rayDir.x;
if (t>max) max=t;
}


Of course rayDir would then have to be the vector to the line segments end point. Anyway, if t>1 then the plane is out of reach. As we use the max distance anyway it doesn't matter if others would be closer and we can still abort the test right there. But I would probably run some tests if the even less cache friendly version is doing any good.

Share this post


Link to post
Share on other sites
haegarr    7372
Ah, you use the sign of the x co-ordinate (as an example) to drop the backfacing planes. _That_ is the primary sense behind the "if" clause. That is truly a nice thing, too ;-)

For the moment that leads to 1 to 3 candidate planes as the full Woo algo also does. However, inside the clause one furthur investigates only the nearest plane (w.r.t. x). This is possible because the other 2 planes (if not parallel) will be investigated by the other "if" clauses, since their respective clauses will pick them accordingly (i.e. they are frontfaced w.r.t. the y and z co-ordinates of the ray). Due to the maximum finding on-the-fly, the farthest hit distance will survive.

That's really tricky, and I like it ;-) So you're right in the sense that 6 "if" clauses (*) with 1 plane investigation inside each one is the optimum one could reach. Otherwise one would get only 3 "if" clauses, but each one with 3 plane investigations inside.
(*) from which only at most 3 will be entered due to nearly pairwise mutual exclusive conditions

Quote:
Original post by Trienco
and nestling the <0 >0 test inside yet another !=0 test didn't seem too useful.

I agree definitely.

Quote:
Original post by Trienco
However, it came to me, that checking for a line segment is much more likely than checking for an infinite ray.

Hmm. That method requires the ray's track to be of a suitable length. As fas as I see is both the size of the BB as well as the origin of the ray unrestricted, and hence a suitable length unknown. Computing it would introduce possibly more effort than could be saved. Choosing an arbitrary but great length means to very seldom get the case of early dropping. No, I'm not convinced that this is better than to use the infinite ray.

[Edited by - haegarr on November 27, 2005 6:30:37 AM]

Share this post


Link to post
Share on other sites
Trienco    2555
Quote:
Original post by haegarr
No, I'm not convinced that this is better than to use the infinite ray.


Not in a situation, where you usually have a very long segment anyway. For picking it might not be useful (unless the world is really huge compared to the view distance and absolutely no other method to limit the number of box tests is used). For collisions maybe, as the movement is usually very small per frame. But as I said, I would test it with a profiler and real world data sets to see if it hurts more than it helps.

Though against static objects I would think about calculating once when I would hit some object and only update when my direction or velocity changes. Else all the repeated small step tests are a bit like what I will just call the "are we there yet"-syndrome.

Well, everything at its own time. Back to quaternions and trying to build a "better" binary md5 bastard format.

Share this post


Link to post
Share on other sites
Endar    668
Quote:
Original post by haegarr
You may take into account to post the belonging code for inspection by another pair of eyes.

......

What components do you mean? Inverting the ray should not change anything, since the hit point will be the same.


1. Here's my code for inspection.
2. The components I was talking about were the components of the vector that represent the intersection point.


#define PLANE_EPSILON 0.0001f



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

class CPlane
{
public:

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.


/**
* Default constructor
*/

CPlane()
: m_normal(0,0,0), m_distance(0)
{ }

/**
* 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
}

/**
* Assignment operator
* \param o A const reference to a CPlane object.
* \return A reference to 'this' CPlane object.
*/

CPlane& operator=(const CPlane& p)
{
if( this != &p ){
m_normal = p.m_normal;
m_distance = p.m_distance;
}

return *this;
}

/**
* 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) const
{
/*
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) const
{
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) const
{
float t;

/* 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

// get the distance from the origin point of the ray along direction 'd' to 'this' plane
if( !distance(o,dir,t) )
return false; // no intersection point

// calcuate the intersection point by starting with the origin of the ray and
// scaling the ray's direction vector by the distance between the start point and
// the plane.
intersection = o + (dir * t);

return true;
}


/**
* Returns 'true' if there is an intersection between the ray and 'this' plane.
* \param o The origin point of the vector.
* \param d The direction of the vector (normalized)
* \param dist The distance from the origin point to 'this' plane
* \return 'true' if the distance is calculated, else 'false
*/

bool distance(const util::CVector3d& o, const util::CVector3d& d, float& dist) const
{
// 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
// this finds the distance between the origin point and the plane in the direction
// of 'd'
dist = (-(m_normal.dotProduct(o) + m_distance)) / m_normal.dotProduct(dir);


return true;
}

};







/**
* Default constructor
* NOTE:: min is the bottom left corner toward the negative z-axis
* NOTE:: max is the top right corner toward the positive z-axis
*/

aabbox::aabbox()
: min(util::CVector3d(-1,-1,-1)), max(util::CVector3d(1,1,1))
{

}

/**
* Constructor
* NOTE:: min is the bottom left corner toward the negative z-axis
* NOTE:: max is the top right corner toward the positive z-axis
* \param nmin The new min value of the AABB
* \param nmax The new max value of the AABB
*/

aabbox::aabbox(const util::CVector3d& nmin, const util::CVector3d& nmax )
: min(nmin), max(nmax)
{

}

/**
* Copy constructor
* NOTE:: min is the bottom left corner toward the negative z-axis
* NOTE:: max is the top right corner toward the positive z-axis
*/

aabbox::aabbox(const aabbox&a)
: min(a.min), max(a.max)
{

}

/**
* Overloaded assignment operator.
* \param a A const reference to an aabbox object
* \return A reference to 'this' aabbox object
*/

aabbox& aabbox::operator=(const aabbox& a)
{
if( &a != this ){
min = a.min;
max = a.max;
}

return *this;
}

/**
* Add a point to the inside of the AABB and increase
* the AABB max and min points if necessary to ensure that
* point <= max && point >= min
*/

void aabbox::addPoint(const util::CVector3d& p)
{
if( p.x > max.x ) max.x = p.x;
if( p.y > max.y ) max.y = p.y;
if( p.z > max.z ) max.z = p.z;

if( p.x < min.x ) min.x = p.x;
if( p.y < min.y ) min.y = p.y;
if( p.z < min.z ) min.z = p.z;

}

/**
* Check if a point is within the bounds of the aabbox.
* \param p A const reference to a CVector3d object
* \return 'true' if the point is within the AABB, else 'false'
*/

bool aabbox::checkPoint(const util::CVector3d& p) const
{
// if point 'p' is within the AABB bounds
if( p.x >= min.x && p.x <= max.x )
if( p.y >= min.y && p.y <= max.y )
if( p.z >= min.z && p.z <= max.z )
return true;

return false;
}

/**
* Reset the bounding box.
* \param initValue a reference to a CVector3d object that is the new min and max.
*/

void aabbox::reset(const util::CVector3d& initValue)
{
min = initValue;
max = initValue;
}

/**
* Return if the aabbox and the passed parameter aabbox intersect.
* \param b Another aabbox to test for intersection with 'this' aabbox
* \return 'true' if the two boxes intersect, else 'false'
*/

bool aabbox::intersection(const aabbox& b) const
{
if( (this->min.x > b.max.x) || (b.min.x > this->max.x) )
return false;

if( (this->min.y > b.max.y) || (b.min.y > this->max.y) )
return false;

if( (this->min.z > b.max.z) || (b.min.z > this->max.z) )
return false;

return true; // the boxes intersect
}

/**
* Return true if there is an intersection between the passed ray
* and 'this' aabbox.
* This creates 6 planes, one for each side of the aabb and then uses
* Woo's method to find if there is an intersection with the ray.
* http://gamealgorithms.tar.hu/ch22lev1sec2.html
*
* \param o A const reference to a CVector3d object that is the origin point of the ray.
* \param d A const reference to a CVector3d object that is the direction of the ray (normalized vector)
* \param intersection A reference to a CVector3d object in which the intersection point will be stored
* if there is a valid point.
* \return 'true' if there is an intersection, else 'false'
*/

bool aabbox::intersection(const util::CVector3d& o, const util::CVector3d& d, util::CVector3d& intersection) const
{
util::CPlane sides[6]; // 6 planes making up each side of the aabb
util::CPlane can[3]; // the 3 candidate planes
util::CVector3d ntemp; // temp vector for holding the aabb normals

// top of the box
ntemp = util::CVector3d(0,1,0);
sides[0] = util::CPlane( ntemp, ntemp.dotProduct(max) );

// front of the box
ntemp = util::CVector3d(0,0,1);
sides[1] = util::CPlane( ntemp, ntemp.dotProduct(max) );

// right of the box
ntemp = util::CVector3d(1,0,0);
sides[2] = util::CPlane( ntemp, ntemp.dotProduct(max) );

// bottom of the box
ntemp = util::CVector3d(0,-1,0);
sides[3] = util::CPlane( ntemp, ntemp.dotProduct(min) );

// left of the box
ntemp = util::CVector3d(-1,0,0);
sides[4] = util::CPlane( ntemp, ntemp.dotProduct(min) );

// back of the box
ntemp = util::CVector3d(0,0,-1);
sides[5] = util::CPlane( ntemp, ntemp.dotProduct(min) );


// select the three different planes
int count=0; // the number of planes to test against
for(int i=0, j=0; (i < 6) && (j < 3); i++){
// If the dot product is less than 0, this indicates that
// the vectors are going in roughly opposite directions, which
// is what we want, because with woo's method, the back-facing faces
// are discarded (faces with normals in roughly the same direction
// as the ray).
// This also ignores planes if they are perpendicular to the ray.
if( d.dotProduct(sides[i].m_normal) < 0 ){
can[j] = sides[i];
j++;
count++;
}
}

// By this point we have at most 3 planes to test against.

// Woo's method consists of choosing the 3 best planes
// (or however many we can get, if less than 3),
// and we find the intersection distances between the ray and
// the planes. The largest of these distances corresponds to a
// potential hit. If a potential hit is found, the intersection
// point is computed. If it is within the AABB, then it is a
// real hit.
float largest = 0.0f;
float dist = 0.0f;
int hit = 0;
for(i=0; i < count; i++){
// get the distance to this plane from the ray origin along the ray
if( can[i].distance(o,d,dist) ){
if( dist > largest ){
largest = dist; // new largest
hit = i;
}
}
}

// By this point we should have a single plane which is the best potential
// hit, according to Woo's method. Now, we calculate the intersection point,
// and if it is within the aabbox, then we have an actual hit, and the coords
// of the hit.

util::CVector3d hit_point;

// get the intersection point between the ray and the specified plane
if( can[hit].intersectionPoint( o, d, hit_point ) ){
// if the intersection point is within the AABB
if( checkPoint( hit_point ) ){
intersection = hit_point; // set the intersection point (a real hit)
return true;
}
}

return false;
}




int main()
{

util::aabbox box;
box.addPoint(util::CVector3d(-1,-1,-1) );
box.addPoint(util::CVector3d(2,1,1) );

// intersection right side, with ray. Intersection should be (2,0,0)
util::CVector3d o(3,0,0);
util::CVector3d d(-1,0,0);

util::CVector3d p;
if( box.intersection(o,d,p) ){
printf("\n Right Intersection found: %.3f %.3f %.3f\n", p.x, p.y, p.z);
}
else
printf("\n Right Nothing found \n");


// intersection left side, with ray. Intersection should be (-1,0.5,0)
o = util::CVector3d(-2,0.5f,0);
d = util::CVector3d(1,0,0);

if( box.intersection(o,d,p) ){
printf("\n Left Intersection found: %.3f %.3f %.3f\n", p.x, p.y, p.z);
}
else
printf("\n Left Nothing found \n");

return 0;
}



It prints the results:


Right Nothing found

Left Intersection found: 1.000 0.500 0.000

Share this post


Link to post
Share on other sites
haegarr    7372
Here comes my first inspection results: class CPlane

IMHO several of your boolean conditions are malformed. Those of normalization may possibly play no role (dependent on the usage) but there are also essentials.
Look at the lines prepended by /* haegarr */ and compare them with your previous code.

bool pointOnPlane(const util::CVector3d& p) const
{
/*
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);
/* haegarr */ return ( dot >= m_distance - PLANE_EPSILON && dot <= m_distance + PLANE_EPSILON );
}

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

float len = dir.getLengthSquared(); // who wants to do an unecessary sqrt?
/* haegarr */ 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);
/* haegarr */ if( dot >= 0.0f - PLANE_EPSILON && dot <= 0.0f + PLANE_EPSILON )
return false; // no intersection

return true; // there is an intersection
/* haegarr */ // even better would be:
/* haegarr */ // return dot < 0.0f - PLANE_EPSILON || dot > 0.0f + PLANE_EPSILON;
}

bool intersectionPoint(const util::CVector3d& o, const util::CVector3d& d, util::CVector3d& intersection) const
{
float t;

/* calculate the intersection point */

util::CVector3d dir(d);

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

// get the distance from the origin point of the ray along direction 'd' to 'this' plane
if( !distance(o,dir,t) )
return false; // no intersection point

// calcuate the intersection point by starting with the origin of the ray and
// scaling the ray's direction vector by the distance between the start point and
// the plane.
intersection = o + (dir * t);

return true;
}

bool distance(const util::CVector3d& o, const util::CVector3d& d, float& dist) const
{
// 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?
/* haegarr */ 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
// this finds the distance between the origin point and the plane in the direction
// of 'd'
dist = (-(m_normal.dotProduct(o) + m_distance)) / m_normal.dotProduct(dir);


return true;
}



E.g. pointOnPlane: If dot-product minus distance nearly 0, then return true. Can be transformed into dot-product nearly distance, what means dot-product at least distance minus tolerance _and_ dot-product at most distance plus tolerance.

Exactly the inverse at the condition of normalization: if the vector's (squared) length is _not_ nearly 1, then ... I've applied the deMorgan rule here.

BTW: Notice that there is at least one comment that tolds that a vector is already normalized, but the routine also works on unnormalized vectors (due to self normalizing). You should hold comments and implementation synchrone.

I look furthur on the code, but I expect the corrections above already changing the overall result.

Share this post


Link to post
Share on other sites
haegarr    7372
Here my 2nd thought (I could not do all at once, sorry):

IMHO the computation of the distance in CPlane isn't correct. Please check my explanation:

The plane is given by the Hess' Normal Form (don't know whether this is the correct term in english)
( r - p ) dot n = r dot n - p dot n = 0
where r is an arbitrary point, p is a point on the plane, and n is the plane's normal.
From this the distance D could be computed as
D := p dot n
so that
r dot n = D
is equivalent.
Now, the arbitrary point is set to the ray
r(t) := o + t * dir
so that
( o + t * dir) dot n = o dot n + t * dir dot n = D
t * dir dot n = D - o dot n
and finally
t = ( D - o dot n ) / ( dir dot n )
If you please compare this with your implementation
dist = (-(m_normal.dotProduct(o) + m_distance)) / m_normal.dotProduct(dir);
then you see a difference at the sign of D (your "m_distance")!

Looking at the site you've cited as the source of the formula, you'll see the plane equation:
A'x + B'y + C'z + D' = 0
Comparing to what I've written above, you see that the D' of this formula is the "inverse distance"
D' := - p dot n
(something that could often be found in literature, but one must be aware of this!

So, please replace your dist formula by
dist = ( m_distance - m_normal.dotProduct(o) ) / m_normal.dotProduct(dir);

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