The following is a compilation of notes on computing the intersection of a ray with a torus, taken from newsgroup articles and Eric Haine's Ray Tracing Newsletter.
For one thing, you can refer to Graphics Gems II for the article on "Intersecting a Ray with an Elliptical Torus", written by Joseph M. Cychosz. Note the following errata:

Not included in the gem "Intersecting a Ray with an Elliptical Torus" is the following code which compensates for the order of the roots being transposed. It was not included since it would make the gem dependent on the behaviour of the root solver. Insert after the call to SolveQuartic (middle of p. 579):

    /* SolveQuartic returns root pairs in reversed order. */
    if  ( *nhits > 1 )
        m = rhits[0]; u = rhits[1]; rhits[0] = u; rhits[1] = m;
    if  ( *nhits > 3 )
        m = rhits[2]; u = rhits[3]; rhits[2] = u; rhits[3] = m;
Also, note that the torus intersector may not handle degenerate surfaces (e.g. where the major radius < minor radius, or major radius < 0) as desired.

Quartic Roots, and "Intro to RT" Errata, by Larry Gritz (vax5.cit.cornell.edu!cfwy)

I was trying to find the solution to a quartic equation to solve a ray-torus intersection test. I've gotten lots of replies, but generally fall into one of two categories: either solve the quartic equation (I forget the reference now, but I'll send you either a reference or the formula if you want [It's in the _CRC Standard Math Tables_ book - Eric Haines]), or by some iterative method. Everybody says (and I have confirmed this experimentally) that solving the equation is very numerically unstable. I have chosen to use the Laguerre method from Press, et. al. _Numerical Recipes in C_, which is slow, but seems to work, and finds all roots without needing to specify brackets for the roots. (An advantage since although I can bracket all possible real roots with the bounding box test that I already do, I'm not really sure how many roots lie within the brackets.)

What actually turns out to be a bigger problem is that I got the quartic coefficients from the SIGGRAPH '88 Ray Tracing Course Notes (on page 3-14 of Pat Hanrahan's section).

[Larry and I thrashed this out over a number of passes (boy, I wish I had access to Mathematica or somesuch), and came out with the following corrected equation set for those on page 93 of An Introduction to Ray Tracing:

        a4 & a3 - Pat's are OK.
        a2 = 2(x1^2+y1^2+z1^2)((x0^2+y0^2+z0^2)-(a^2+b^2))
           + 4 * (x0x1+y0y1+z0z1)^2 + 4a^2z1^2
        a1 = 4 * (x0x1+y0y1+z0z1)((x0^2+y0^2+z0^2)-(a^2+b^2))
           + 8a^2 * z0 * z1
        a0 = ((x0^2+y0^2+z0^2)-(a^2+b^2))^2 - 4a^2(b^2-z0^2)

Eric Haines


Reply from Bob Webber (webber@fender.rutgers.edu):

For planar curves we have J. Dennis Lawrence's A Catalog of Special Plane Curves (Dover 1972) to satisfy those times when one wakes up in the middle of the night, racking one's mind trying to remember the equation for the hippopede. However, for 3-d, the best I have seen is Pat Hanrahan's A Survey of Ray-Surface Intersection Algorithms that appears in Andrew Glassner's An Introduction to Ray Tracing (Academic Press 1989). There we find, among other things, the equation for a torus as:

   (x**2 + y**2 + z**2 - (a**2 + b**2))**2 = 4 a**2 (b**2 - z**2)

This describes a torus centered at the origin defined by a circle of radius b being swept along a center defined by a circle of radius a. It is derived from considering the intersection of a plane with a torus that yields the two circles:

   ((x - a)**2 + z**2 - b**2)((x + a)**2 + z**2 - b**2) = 0

[if you are unfamiliar with this construction, it is worthwhile pausing here and savouring how this equation actually works -- sometimes the equations are prettier than the pictures] and then spinning this intersection by replacing x**2 with x**2 + y**2 (after some algebraic simplification, which converted the above to:

     (x**2 + z**2 - (a**2 + b**2))**2 = 4 a**2 (b**2 - z**2)

The section includes a reference to an unpublished 1982 note by Hanrahan entitled: Tori: algebraic definitions and display algorithms. The general scheme for a number of variations on the torus is to start with a quartic curve of the form f(x**2,y)=0 and then substitute x**2+y**2 for x**2 and z for y.


Correct Roots for Torus Intersection, by Haakan "Zap" Andersson (zap@lage.lysator.liu.se)

When using the torii-intersector you gave me [I sent the code from Graphics Gems II, Joseph Cychosz's torus intersector, I believe. - EAH], we have the trouble of getting the wrong side when rmajor < rminor and above all when rmajor < 0, which are both quite valid and well defined torii shapes.

Ok, here's what you do:

Check if rmajor > rminor.

If so, hit is OK as is.

If not, calculate the distance between the hit you got and torus center. We will be needing the distance SQUARED, so no need for a square root. Put the distance squared in variable SqrDist.

If rmajor < rminor && rmajor > 0
   if SqrDist < (rminor*rminor - rmajor*rmajor)
      then discard the hit, it's on the wrong surface.
If rmajor < rminor && rmajor < 0
   if SqrDist > (rminor*rminor - rmajor*rmajor)
      then discard the hit, it's on the wrong surface.
Voila!
Joe Cychosz (3ksnn64@ecn.purdue.edu) responds:

I did test the case where the torus was degenerate, however I did not consider the application sense of it not really being a wanted surface. I will study this further and incorporate your suggestion.

Also, not include in the GEM is the following code which compensates for the order of the roots being transposed. I didn't include it since it would make the GEM dependent on the behaviour of the root solver.

Insert after the call to SolveQuartic:

/*      SolveQuartic returns root pairs in reversed order.              */
        if  ( *nhits > 1 )
            m = rhits[0]; u = rhits[1]; rhits[0] = u; rhits[1] = m;
        if  ( *nhits > 3 )
            m = rhits[2]; u = rhits[3]; rhits[2] = u; rhits[3] = m;

From hanwen@stack.urc.tue.nl (Han-Wen Nienhuys) / Mar 28 1994

maybe you're interested in this: (or should I publish it in Raytracing news?)

Newsgroups: comp.graphics
Subject: Re: Symbolic Algebra check on Torus/Ray intersect needed.
Robin J Green <R.J.Green-SE1@cs.bham.ac.uk> wrote:
Does anyone have access to a symbolic algebra system at all? I need to check the LONG computation for a ray-torus intersection.
x = (xo + t * xd)
y = (yo + t * yd)
z = (zo + t * zd)

into the implicit equation for a torus (a=major rad, b=minor rad):

(x^2+y^2+z^2-(a^2+b^2))-4a^2(b^2-z^2) = 0

You must have forgotten a ^2 here.

I actually always use the equation

(x^2 + y^2 + z^2 + R^2 - r^2)^2 - 4R^2(x^2+y^2) = 0

For a torus lying in the XY plane, major radius R and minor radius r

If you look closely, you'll notice that the left term of the left hand side looks similar to the equation of a sphere. Filling X = P + t*D (X is a function of t, D and P are vectors, D = (xd,yd,zd) and P = (xo,yo,zo)) in the left side amounts to filling in the ray equation X = P + D*t in the equation of a sphere with radius sqrt(R^2 - r^2). Let me rewrite the equation to

((X,X) + R^2 - r^2)^2 - Stuff = 0

!!!!(let (a,b) denote the dot product between vectors a and b)

first we find (X,X):

(X,X) = (P+t*D, P+t*d) = (P,P) + 2*t*(D,P) + t^2*(D,D)

(we've just used the linearity of the inner-product... sounds good, huh? :) so let

b0 = (P,P) + R^2 - r^2
b1 = 2*(D,P)
b2 = (D,D)

then the left term in the torusequation equals

(b0 + b1*t + b2*t^2)^2 

which is

(1)  (b0^2 + 2*b0*b1*t + (b1^2+ 2*b2*b0)*t^2 + 2*b1*b2*t^3 + b2^2*t^4)

now we have the left part, the right part ("Stuff") works analogously:

let Q=(x,y,0)= Q0 + t*Qd

with Q0 = (x0, y0, 0) and Qd = (xd,yd, 0), then the right part is

4*R^2*(Q,Q)

which is

c0 + c1*t + c2 * t^2

with

c0 =4*R^2* (Q0,Q0), c1 = 8*R^2*(Q0,Qd) and c2 = 4*R^2*(Qd,Qd)

if you substract the c0, c1, c2 from the coefficients of t^? in eq. (1), you get the desired result.

a0 = b0^2 - c0
a1 = 2*b0*b1 - c1
a2 = (b1^2+ 2*b2*b0) - c2
a3 = 2*b1*b2
a4 = b2^2
and rearrange the resulting formula in terms of a quartic in t. Some partial derivitives in x,y and z for normal calculation would be useful too.

very unsmart. Check out the Rayshade code for a better way: to find the normal to a torus in the XY plane, you project the intersection points onto the xy plane, and then you radially procject it onto the heartline (what are they called in english :) ? the difference between the intersection point and the projection is the normal.

algorithmically: let X=(x,y,z) be the intersection point, then

p := (x,y,0);
l := veclen(p);
p := Majrad * p / l;
n := X - p;
normalise(n);

computes the normal n. For actual implementation, see the Rayce raytracer, wuarchive, /graphics/graphics/ray/rayce

BTW, I am typing this without the sources or a piece of paper, so this might still contain a few errors.
a4 = (xd^2+yd^2+zd^2)^2
a3 = 4(xo*xd+yo*yd+zo*zd)*(xd^2+yd^2+zd^2)
a2 = 2*(xd^2+yd^2+zd^2)*((xo^2+yo^2+zo^2)-(a^2+b^2))+
                                      4*(xo*xd+yo*yd+zo*zd)-4*a^2*zd^2
a1 = 4*(xo*xd+yo*yd+zo*zd)*((xo^2+yo^2+zo^2)-(a^2+b^2))+8*a^2*zo*zd
a0 = ((xo^2+yo^2+zo^2)-(a^2+b^2))^2-4a^2(b^2-zo^2)

I must be stupid. Why does anybody use those hairy formulas? The symbolic way to get those formulas is a lot easier to read, you make less mistakes, and is likely to be more efficient.