Jump to content
  • Advertisement
Sign in to follow this  
havresylt

cube-ray intersection

This topic is 3408 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hi, I'm new here, and I have a question. I'm currently doing a ray tracing engine where I want to test for intersections between rays and cubes. A ray is represented by 2 vectors, a point p, and a (normalized) direction d. A cube is represented by a center vector, (normalized and orthogonal) up, right and out vectors, and a side length. Here is how I thought when implementing the function: A cube is determined by 6 inequalities, where a point x is inside the cube if all three are satisfied simultaneously: -side/2 <= (x - center) dot (right) <= side/2 -side/2 <= (x - center) dot (up) <= side/2 -side/2 <= (x - center) dot (out) <= side/2 Now for a ray: x = p + a*d where 'a' is a scalar. If this is inserted above: -side/2 <= (p + a*d - center) dot (right) <= side/2 -side/2 <= (p + a*d - center) dot (up) <= side/2 -side/2 <= (p + a*d - center) dot (out) <= side/2 Now we want to determine if it is possible to find an 'a' that satisfies all of them. Rearranging: (center - p) dot (right)-side/2 <= a*(d) dot (right) <= (center - p) dot (right) + side/2 (center - p) dot (up)-side/2 <= a*(d) dot (up) <= (center - p) dot (up) + side/2 (center - p) dot (out)-side/2 <= a*(d) dot (out) <= (center - p) dot (out) + side/2 Making sure the '(d) dot (X)' factors are non-zero, we can divide by them: [(center - p).dot(right)-side/2]/d.dot(right) <= a <= [(center - p).dot(right) + side/2]/d.dot(right) [(center - p).dot(up)-side/2]/d.dot(up) <= a <= [(center - p).dot(up) + side/2]/d.dot(up) [(center - p).dot(out)-side/2]/d.dot(out) <= a <= [(center - p).dot(out) + side/2]/d.dot(out) Now, finding that largest value of the left-hand sides, and the smallest one of the right-hand sides, we can determine a single interval in which 'a' has to lie: max(LHS) <= a <= min(RHS) Now there are a couple of cases. If max(LHS) > min(RHS) there can't possibly be any intersections. If both of them are negative, the entire cube lies in the wrong direction. Otherwise, we want to find the smallest positive value of 'a'. ------ To me, this seems to be a valid way. It doesn't work however. I get all kinds of strange results. Sometimes I get intersections, and sometimes not, depending on the direction of right and up (which shouldn't have an effect). I'm pretty sure I've managed to implement it correctly, so I would be very grateful for some input of what could be my mistake. Perhaps, there might be a smarter way to do this altogether (eg. by making use of 6 quads). I still would like to get what I'm doing wrong.

Share this post


Link to post
Share on other sites
Advertisement
Quote:
Original post by havresylt
Making sure the '(d) dot (X)' factors are non-zero, we can divide by them:


I think this is where you are going wrong.

The thing is that (d) dot (X) can be negative, and if you multiply/divide your equation by a negative number, you'd have to flip the <= into a >=

Share this post


Link to post
Share on other sites
I personally prefer the following way (using column vectors here):

(1) The ray is given as
r(t) := r0 + t * d
in the global co-ordinate system.

(2) The cube's local co-ordinate frame (w.r.t. the global co-ordinate system) is given as
M

Inside this frame the cube has the origin 0, and the normals of the 6 planes of the cube are
(1,0,0)t, (-1,0,0)t, (0,1,0)t, (0,-1,0)t, (0,0,1)t, (0,0,-1)t

The cube's half edge length should be h.

(3) Then the ray can be transformed into the cube's local frame
r(t)' := M-1 * r(t) = r0' + t * d'

Using this ray has the following advantages: A translation, rotation, and scaling of the cube is considered; the planes of the cube are the principal unit vectors and hence reduce the computational effort (see below).

(4) Compute the travel distance on the ray's track for a side of the cube. E.g. for the right side which is given by the plane equation
sR := ( p - h * x ) . x = 0
where p denotes an arbitray point and x the unit vector in x direction.

Setting in the local ray gives us
( r0' + tR * d' - h * x ) . x = 0
( r0' + tR * d' ) . x = h
r0,x' + tR * dx' = h
due to the nature of x. This can simply be solved for tR
tR = ( h - r0,x' ) / dx'
which must be protected against division by zero, of course. If dx' is zero then the ray is parallel to the plane.

(5) Assume that a hit was detected in (4). Setting in tR into the local ray gives us a hit point in the local frame
r(tR)'
Actually we need not compute the x component of the above since it will be h. But the y component and the z component need to be inside [-h,+h] to hit the cube's side (else the infine plane in which the side lies is hit).

(6) Perform (4) and (5) for all sides.


Please check the math twice, I'm not sure I've done all well here.

[Edited by - haegarr on January 22, 2009 3:46:52 PM]

Share this post


Link to post
Share on other sites
Quote:
Original post by quasar3d
The thing is that (d) dot (X) can be negative, and if you multiply/divide your equation by a negative number, you'd have to flip the <= into a >=


Thanks! This a very good point! :)

I've fixed it by swapping the right- and left-hand sides in the current inequality. Unfortunately it still doesn't make any difference though.


Quote:
Original post by haegarr
I personally prefer the following way (using column vectors here):

...


Thank you for your long and thorough answer!

Well, to me this is the same thing. First you transform the vector into the system of the cube (subtracting by the center of the cube, and then taking projections on its orthogonal axes). These are the quantities I called '(x-center).dot(X)'; where X={right, up, out}.

Then, adopting your notation, one sets x = r0 + t*d and solves for t. In my case, I'm combining all the inequalities and try to solve them simultaneously, which seems OK, but it could be that I'm doing something wrong there.

Any other ideas?

Share this post


Link to post
Share on other sites
Quote:
Original post by havresylt
Quote:
Original post by quasar3d
The thing is that (d) dot (X) can be negative, and if you multiply/divide your equation by a negative number, you'd have to flip the <= into a >=


Thanks! This a very good point! :)

I've fixed it by swapping the right- and left-hand sides in the current inequality. Unfortunately it still doesn't make any difference though.


You should only swap it if the result of d dot X is negative, right?

Share this post


Link to post
Share on other sites
Quote:
Original post by quasar3d
Quote:
Original post by havresylt
Quote:
Original post by quasar3d
The thing is that (d) dot (X) can be negative, and if you multiply/divide your equation by a negative number, you'd have to flip the <= into a >=


Thanks! This a very good point! :)

I've fixed it by swapping the right- and left-hand sides in the current inequality. Unfortunately it still doesn't make any difference though.


You should only swap it if the result of d dot X is negative, right?


yes of course. sorry if I wasn't clear about that.

I still haven't got any further. I might try copying haegarr's solution literally (although it's the same as mine imo).

Share this post


Link to post
Share on other sites
You could of course copy his solution, but I'm actually quite curious as well why your solution won't work, so maybe you could post the actual code?

Share this post


Link to post
Share on other sites
Here's my code, hope you find an error :P


float Cube::intersect(const Ray& r) const{

Vec3f diff = center - r.getOrigin();
float max = MAX_FLT;
float min = -max;
float DOT, A, B;

if(DOT = r.getDirection().dot(right) != 0 ){
A = (diff.dot(right) + side/2)/DOT;
B = (diff.dot(right) - side/2)/DOT;
if(DOT <0){
float tmp = A;
A = B;
B = tmp;
}
max = A < max ? A : max;
min = B > min ? B : min;
}

if(DOT = r.getDirection().dot(up) != 0 ){
A = (diff.dot(up) + side/2)/DOT;
B = (diff.dot(up) - side/2)/DOT;
if(DOT <0){
float tmp = A;
A = B;
B = tmp;
}
max = A < max ? A : max;
min = B > min ? B : min;
}

if(DOT = r.getDirection().dot(out) != 0 ){
A = (diff.dot(out) + side/2)/DOT;
B = (diff.dot(out) - side/2)/DOT;
if(DOT <0){
float tmp = A;
A = B;
B = tmp;
}
max = A < max ? A : max;
min = B > min ? B : min;
}
if(min < max){
if(min > 0)
return min;
else if(max > 0)
return max;
}

return -1.0f; //miss
}



Share this post


Link to post
Share on other sites
Quote:
Original post by havresylt
I might try copying haegarr's solution literally (although it's the same as mine imo).

They are similar but not exactly the same. 1st my approach doesn't suffer from any inequality reversal, because it computes actual hit points and checks whether they are in the valid range. Your approach works on the track of the ray instead. And 2nd is the implementation of my approach very straight-forward if the matrix inversion and vectors transforms are out-sourced to a math library (what they should).

The simplicity of testing the ray against axis parallel planes is also well known in the domain of bounding volumes, namely Axis-aligned bounding boxes (AABBs).

Your geometry (i.e. a cube) is simple here, and hence any overhead is visible on the intersection test. But the trick to transform the ray into the object local space is often worth a look. Think e.g. of an ellipsoid with an arbitrary position and orientation. Using the local space trick, the ellipsoid can be virtually converted to a unit sphere around 0 so that the intersection test itself becomes trivial.

Share this post


Link to post
Share on other sites
Quote:
Original post by havresylt
Here's my code, hope you find an error :P

Maybe I'm wrong, but I have a problem with your max/min computation. IMHO you compute a minimum in the max variable and a maximum in the min variable.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

Participate in the game development conversation and more when you create an account on GameDev.net!

Sign me up!