Ray-Quadric Intersection

Represent the ray as:

x = a + t b

where a, b and x are 3-vectors

Represent the quadric as:

x' A x + x' B + C = 0

where A is a symmetric 3x3 matrix, B a 3x1 vector, and C a scalar. Note that x' A x is the dot product of x and A x.

Substituting in gives the quadratic equation:

              2
        l2 * t  + l1 * t + l0 = 0

where:

	l2 = b' A b

        l1 = 2 b' A a + b' B

        l0 = a' A a + b' B + C

If the discriminant of this equation is less than 0, the line does not intersect the quadric. Solving the equation and substituting the values of t into the ray equation will give you the points.

int Torus::line_intersect(float ax, float ay, float az, 
			  float bx, float by, float bz,
			  int * num_intersections,
			  float * intersections,
			  float * shade) const
//              
// Intersect ray
//          [ax]     [bx]              
//	    [ay] + t [by]
//	    [az]     [bz]
// with torus at origin, major radius R, minor radius r
//
{

  // This struct is syntactic sugar so that a.b,
  // (the dot product) looks just right (:-)
  struct { float a,b; } a;
  a.b = ax*bx + ay*by + az*bz;
  a.a = ax*ax + ay*ay + az*az;

  // Set up quartic in t:
  //
  //  4     3     2
  // t + A t + B t + C t + D = 0
  //
  float R2 = R*R;
  float K = a.a - r*r - R2;
  float A = 4*a.b;
  float B = 2*(2*a.b*a.b + K + 2*R2*bz*bz);
  float C = 4*(K*a.b + 2*R2*az*bz);
  float D = K*K + 4*R2*(az*az - r*r);

  // Solve quartic...
  double roots[4];
  int nroots = solve_quartic(A,B,C,D,roots,SOLVE_ALL);

  *num_intersections = 0;
  while(nroots--)
    {
      float t = roots[nroots];
      float x = ax + t*bx;
      float y = ay + t*by;
      float l = R*(M_PI/2 - atan2(y,x));
      if (l <= vlength && l >= 0)
        intersections[(*num_intersections)++] = t;
    }
}