   cube-ray intersection
  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/2Now for a ray:x = p + a*dwhere '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/2Now 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/2Making 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.   
Quote:Original post by havresyltMaking 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 >=;  I personally prefer the following way (using column vectors here):(1) The ray is given asr(t) := r0 + t * din the global co-ordinate system.(2) The cube's local co-ordinate frame (w.r.t. the global co-ordinate system) is given asMInside 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)tThe cube's half edge length should be h.(3) Then the ray can be transformed into the cube's local framer(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 equationsR := ( p - h * x ) . x = 0where 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 = hr0,x' + tR * dx' = hdue to the nature of x. This can simply be solved for tRtR = ( 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 framer(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];Quote:Original post by quasar3dThe 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 haegarrI 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?;Quote:Original post by havresyltQuote:Original post by quasar3dThe 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?;Quote:Original post by quasar3dQuote:Original post by havresyltQuote:Original post by quasar3dThe 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).   ;  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?   ;  Here's my code, hope you find an error :Pfloat 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};Quote:Original post by havresyltI 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.   ;Quote:Original post by havresyltHere's my code, hope you find an error :PMaybe 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.   