Need ray-sphere intersection code

[MattDiamond]

This is going to sound incredibly lazy, but I really need code to compute the intersection of a ray and sphere. I can convert whatever structures it uses to my own, never fear.

There are a billion descriptions of how to do this on the 'net. Half of them describe an algebraic solution (solve a quadratic equation) which I suspect isn't as good for me as the geometric approach which uses vectors, dot products and so on. Most of the references have no code; some seem to disagree on whether they need to check for certain special cases explicitly or not. I found a function somewhere but my testing indicates there may be a problem.

I could write my own given a little more time and a little more sleep but I'm out of both of those. :-( So, please, debugged code or *robust* pseudo-code, anyone?

Or, I could post what I'm using and someone can tell me whether it looks reasonable or not.

Thanks! Just a few more hurdles and I'll finally be able to post an alpha version of Asteroid Rally for people to try out...

[Jammer]

Check out the E3Ray3D_IntersectSphere() function in the Quesa library. It's in the E3Math.c file. Quesa is here. (http://quesa.designcommunity.com/index.html)

[DoG]

I'll give it a shot.

Basically, the dot product and cross product approach is all you need.

The center of the sphere is at C, its radius is r, your ray starts at V and has direction D


1. Map the (C-V) vector onto the D vector. (project (C-V) on D)=CV'
2. Find the vector perpendicular to D that points to C: ((C-V) - CV')=R
3. if the magnitude of R is less than the radius of the sphere, you have an intersection. if (abs(R) <= r) then intersection

I cannot remember how to find the exact points of intersection, but this should get you started. The point (C+R) should mark the middle of the 2 intersections, and I guess there is some trig to figure out the rest.

[MattDiamond]

Amazing! Both your answers are better than anything I found with Google. I've had pretty good luck getting hints via Google for most topics, but for some reason I came up short on this one. Too much of a hurry, maybe.

I'll poke around with these (and any other answers I get) tonight, and hopefully this will get me past my hump! Fingers crossed...

Thanks!

[MattDiamond]

Originally posted by DoooG

I cannot remember how to find the exact points of intersection, but this should get you started. The point (C+R) should mark the middle of the 2 intersections, and I guess there is some trig to figure out the rest.

Yes, that much I had figured out. One method is that the distance d from (C+R) to the intersection points can be computed from the Pythagoras theorem:
d = sqrt(r^2 - len(C+R)^2)
Maybe there are ways I could avoid taking the square root but for the number of collisions I expect each frame (usually 0, sometimes 1 or 2) it's not a problem.

It's weird, I had the worst time with this problem. I knew the pieces but I could not for the life of me put together the whole solution. "Mental block" is the only way to describe it. Now I can feel the block falling away. I may or may not use the Quesa code, but either way I feel much better now that I understand how to do it myself. Thanks, Dooog and Jammer...

[codemattic]

Glad you got a handle on it. But let me pile on with one more link (even though its a little late) - <http://wild-magic.com/>. Eberly's library has tons of collision/intersection code. If you find it useful make sure to buy his book!

cheers

[MattDiamond]

The wild-magic license is clear-cut: I can use the code, period. I may yet do that. But using theirs would require more work. The Quesa code is very close to what I need, but their license doesn't seem to talk about the terms under which I can use that one function, whether as is, or with modifications. They seem more concerned with the library as a whole...

Does anyone know whether I can crib this one function from Quesa, and under what terms? I want to do the right thing. I will be changing the data structures and adding a case for the end of my line segment, but the logic would be mostly theirs.

[Johan]

Here is the code I use for sphere intersection, maybe it helps.
float ray_sphereIntersection(sphere_t* sphere, ray_t* ray) {
float b = vdot(ray->d, vsub(ray->o, sphere->p));
float c = vsquaredist(ray->o, sphere->p) - (sphere->r * sphere->r);
float d = b * b - c;

if (d < 0.0f) {
return -1.0f;
}
return -b - sqrt(d);
}

.johan

[Jammer]

Originally posted by MattDiamond
I will be changing the data structures and adding a case for the end of my line segment, but the logic would be mostly theirs.

Actually, "the logic" is not theirs as it comes from 'Real Time Rendering', section 10.3.2., as they note at the top of the function. Frankly, as long as you are making changes to the code, and you will have to to get your structures to fit, then I think you're fairly safe since the whole library is LGPL'd anyway.

[MacFiend]

magicsofware.com .... gee they aint lying :p
/me has been touched by God :wow:

[codemattic]

Originally posted by MacFiend
magicsofware.com .... gee they aint lying :p
yeah - and his game physics book/cd/library comes out in December!

Codemattic (who is saving his pennies)

[lpetrich]

OK, I'll give the complete formulas here:

I'll Use DoooG's notation:

sphere with center C and radius R
line with point V and direction D

The vector from the sphere's center to the line point S = V - C

The distance along the line is

L = - (S.D)/(D.D) +/- sqrt((R^2-(S.S))/(D.D) + ((S.D)/(D.D))^2)

and the intersection point's position is V + L*D

The square-root term can be expressed in a way that may be more convenient for calculation, if the line's direction is close to its point's direction to the circle's center.

Let p = (S x D)^2 -- the absolute square of the cross product of S and D. Then

L = ( - (S.D) +/- sqrt((D.D)*R^2 - P) )/(D.D)

The discriminant (the term in the square root) can be used as a test for intersection.

[lpetrich]

I'm guessing that the sphere is being used as a bounding box. Other shapes have been used as bounding boxes, notably vertical cylinders and rectangular solids.

The circle part of a cylinder can be handled with the same mathematics as with the sphere, except that the cross-product is a scalar, not a vector quantity:

a x b = a[1]*b[2] - a[2]*b[1]

An object with a bounding box can have its bounding box follow its rotations by making its bounding box also be rotated (oriented bounding box). One has to calculate a collision with each of its six bounding planes, thoough the calculation is somewhat simplified by those planes coming in parallel pairs.

One can simplify a bit further by constructing an axis-aligned one; find the minimum and maximum coordinate values of all the vertices of an oriented bounding box and use those to construct a new bounding box with bounding planes along the coordinate axes.

[MattDiamond]

Actually, the sphere is the very thing I'm intersecting with, it's not a bounding box. Aren't simplifying assumptions grand? :)

In the end I got the Quesa/Real Time Rendering version of the algorithm working very quickly. Thanks to everyone's help I've finally posted a playable version of Asteroid Rally (see thread in the uDG 2003 forum for info.)

Thanks, all!

[DoG]

As it happens, I am working on the subject, and have just finished a tidbit about intersection detection.

http://web.axelero.hu/g5d6/yag/collision_ray-tri_ray-sphere.pdf

It has ray/triangle and ray/sphere collision details, outlining the steps to be taken, on a whole 3 pages. No code, but implementing it should be easy. I don't know if the paper contains the fastest methods, but I think they are logically correct.

If there is interest, I am planning on implementing ray/cylinder, sphere/triangle, cylinder/sphere, cylinder/triangle, ray/box, box/cylinder, sphere/box, and I would document it too, in a similar format, some graphs thrown in when I get the chance to make them.

Any feedback welcome, cheers

[MattDiamond]

Is the purpose to present easily understandable/implementable methods, rather than speed? Seems worthwhile to me.

Note that there may be a cheaper way to test whether a point V is in a triangle (Px,Py,Pz). From memory, I think you can add up the angles that the point V makes with each pair of triangle vertices. I.e.
acos( PxV dot VPy ) + acos( PyV dot VPz ) + acos ( PzV dot VPx )

If the point is in the triangle, the angles add up to 360 degrees (2*pi). Cost is 3 acos's, plus some dot products.

First, I don't know if that's a cheaper solution overall than yours, and second, if the point "in" the triangle isn't exactly coplanar with the triangle vertices you won't get exactly 360 degrees. Nonetheless, I thought I'd mention this... (Hmm, maybe there's an even faster way by checking the dot products directly rather than taking the acos. In effect, you just want to know that V is "in between" every pair of triangle legs.)

Well I guess someone's always going to know a faster algorithm. I just googled and found one that uses cross products and supposedly beats the pants off the 360-degree test:
http://www.blackpawn.com/texts/pointinpoly/default.html
Not sure what you want to do about this though. To not be seen as trying to compete with the magicsoftware's and Graphic Gem's of the world, maybe you just need to make your purpose clear that you are not necessarily going after speed. At the same time, people will want to know that your solutions aren't hideously slow, yes? In which case you may want to glance at Fundamentals of Comp Graphics, or Real Time Rendering, or Graphics Gems, just to make sure there isn't some widely known solution that you hadn't considered.

But if the point of the exercise is that you feel like working this stuff out for yourself, then ignore the above ramblings.

[DoG]

Originally posted by MattDiamond
Is the purpose to present easily understandable/implementable methods, rather than speed? Seems worthwhile to me.

Note that there may be a cheaper way to test whether a point V is in a triangle (Px,Py,Pz). From memory, I think you can add up the angles that the point V makes with each pair of triangle vertices. I.e.
acos( PxV dot VPy ) + acos( PyV dot VPz ) + acos ( PzV dot VPx )

If the point is in the triangle, the angles add up to 360 degrees (2*pi). Cost is 3 acos's, plus some dot products.



Doesnt work that way, because you have to take the unit length vectors for the acos() to work. Thats where the normalization comes from, though it is only a single division in the end, not 6 per acos(). The 3 sqrt() and 4 acos() are the expensive parts. But as said in the paper, this is the very last step and infrequently executed, so i dont know how much optimizing it brings.

But if the point of the exercise is that you feel like working this stuff out for yourself, then ignore the above ramblings.

Yea, I am implementing all this, and I have to write it down, wont remember it in a few months otherwise. I am writing things for me to understand, basically, but of course it would be nice if its useful to others. Simple stuff first, and who is ready can dig (and understand) the advanced stuff.

PS: I looked at that other test, and it seems like he is just shifting the sqrt to another part to the SameSide(p,a, b,c) because he has to use vector projection against the normal to find out what is on the same side. Though it seems like no more acos(). I will look into this, too, seems worthwile to do.

[lpetrich]

First, I wish to apologize to MattDiamond for misunderstanding what he had wanted to do. In a bboard dedicated to game creation, I'd expected the discussion to be about simplified ways of handling 3D models, and ball-shaped objects are usually not the most common objects in such games.

As to the triangle problem, there are several ways to do it.

One way is to see whether a point is on the same side of all the sides. Here is a sidedness test:

Triangle has vertices v1, v2, and v3

Sidedness relative to v1-v2, v2-v3, and v3-v1 sides:

S1 = (v - v1)x(v2 - v1)
S2 = (v - v2)x(v3 - v2)
S3 = (v - v3)x(v1 - v3)

where x is the 2D cross product.

For v to be inside the triangle, S1, S2, and S3 must have the same sign. Which sign will depend on whether one had done a clockwise or counterclockwise choice of vertices v1, v2, and v3.

The area of this triangle is (1/2)(v2 - v1)x(v3 - v1) with the sign depending on the vertex direction.

The above can easily be generalized to convex polygons with any number of vertices; one can even test for the convexity at a vertex by finding the area of the triangle with it and its two neighbors.

[Jesse]

The most robust and efficient way I've found of doing ray-triangle intersections is to use Plucker coordinates. I haven't done any tests, but I assume it's faster than the other methods 'cause it doesn't involve any sqrts or trig functions. It also makes for nice, clean code. Best of all, it's very robust in that there are no special cases, and rays don't fall 'in the cracks' between tris.

On another topic, I'm also trying to come up with some code for ray-cylinder intersection. I gather that there are some quadratic equations involved, but I've never solved these in code and I'm not sure how to do it. I'm working through Dave Eberly's code, but haven't quite deciphered it yet. So I too would be interested in any info on ray-cylinder intersection algorithms.

[DoG]

Originally posted by Jesse

On another topic, I'm also trying to come up with some code for ray-cylinder intersection. I gather that there are some quadratic equations involved, but I've never solved these in code and I'm not sure how to do it. I'm working through Dave Eberly's code, but haven't quite deciphered it yet. So I too would be interested in any info on ray-cylinder intersection algorithms.

I was thinking about this. my approach would be to project the cylinder and ray into a plane for which the normal is the center line of the cylinder.

Then the solution is reduced to finding a circle with a line intersection, which is really simple.

The same works for two cylinders, in which case we'd get a circle and a box.

[Jesse]

I was thinking about this. my approach would be to project the cylinder and ray into a plane for which the normal is the center line of the cylinder.

Then the solution is reduced to finding a circle with a line intersection, which is really simple.

The same works for two cylinders, in which case we'd get a circle and a box.


This would only give you true-false, right, not the points of intersection in 3D? I'm hoping to get the points of intersection also...

Eberly's approach first transforms the ray into local space, in which the cylinder axis is the z axis. This simplifies things somewhat, but it's still pretty complicated.

[DoG]

Originally posted by Jesse
This would only give you true-false, right, not the points of intersection in 3D? I'm hoping to get the points of intersection also...

Eberly's approach first transforms the ray into local space, in which the cylinder axis is the z axis. This simplifies things somewhat, but it's still pretty complicated.

No, it would actually give the intersection coordinates, and from what you say, Eberly's way is the exact same thing I suggested.

You know the length of the ray/line segment in local coordinates, and the intersection points. You can get the ratio of intersection to segment length, which you can apply to the original distance vector of the ray to find the intersection points in 3D space.

[lvnguyen]

OK, I'll give the complete formulas here:

I'll Use DoooG's notation:

sphere with center C and radius R
line with point V and direction D

The vector from the sphere's center to the line point S = V - C

The distance along the line is

L = - (S.D)/(D.D) +/- sqrt((R^2-(S.S))/(D.D) + ((S.D)/(D.D))^2)

and the intersection point's position is V + L*D

The square-root term can be expressed in a way that may be more convenient for calculation, if the line's direction is close to its point's direction to the circle's center.

Let p = (S x D)^2 -- the absolute square of the cross product of S and D. Then

L = ( - (S.D) +/- sqrt((D.D)*R^2 - P) )/(D.D)

The discriminant (the term in the square root) can be used as a test for intersection.

Does anyone know if this derivation is from a paper or text? I want to be able to reference it