What is the fastest way to find the point of intersection between a ray and a polygon?

Pretty much as the question asks. Answers preferably in pseudo code and referenced. The correct answer should value speed over simplicity.

From stackoverflow

• See Intersections of Rays, Segments, Planes and Triangles in 3D. You can find ways to triangulate polygons.

  If you really need ray/polygon intersection, it's on 16.9 of Real-Time Rendering (13.8 for 2nd ed).

  We first compute the intersection between the ray and [the plane of the ploygon] pie_p, which is easily done by replacing x by the ray.

   n_p DOT (o + td) + d_p = 0 <=> t = (-d_p - n_p DOT o) / (n_p DOT d)
  

  If the denominator |n_p DOT d| < epsilon, where epsilon is very small number, then the ray is considered parallel to the polygon plane and no intersection occurs. Otherwise, the intersection point, p, of the ray and the polygon plane is computed: p = o + td. Thereafter, the problem of deciding whether p is inside the polygon is reduced from three to two dimentions...

  See the book for more details.

  Chris Lloyd : This was more of a general reference question, but its true that soft surfer has the most up to date algorithm.

• Whole book chapters have been devoted to this particular requirement - it's too long to describe a suitable algorithm here. I'd suggest reading any number of reference works in computer graphics, in particular:

  • Introduction to Ray Tracing, ed. Andrew S. Glassner, ISBN 0122861604
  MrZebra : Good book, I wrote my own ray tracer following it.

  Alnitak : Me too :) It doesn't do polygons yet though!

• For ray-triangle/square here's a very good explanation
  http://www.blackpawn.com/texts/pointinpoly/default.html

• struct point
  {
      float x
      float y
      float z
  }
  
  struct ray
  {
      point R1
      point R2
  }
  
  struct polygon
  {
      point P[]
      int count
  }
  
  float dotProduct(point A, point B)
  {
      return A.x*B.x + A.y*B.y + A.z*B.z
  }
  
  point crossProduct(point A, point B)
  {
      return point(A.y*B.z-A.z*B.y, A.z*B.x-A.x*B.z, A.x*B.y-A.y*B.x)
  }
  
  point vectorSub(point A, point B)
  {
      return point(A.x-B.x, A.y-B.y, A.z-B.z) 
  }
  
  point scalarMult(float a, Point B)
  {
      return point(a*B.x, a*B.y, a*B.z)
  }
  
  bool findIntersection(ray Ray, polygon Poly, point& Answer)
  {
      point plane_normal = crossProduct(vectorSub(Poly.P[1], Poly.P[0]), vectorSub(Poly.P[2], Poly.P[0]))
  
      float denominator = dotProduct(vectorSub(Ray.R2, Poly.P[0]), plane_normal)
  
      if (denominator == 0) { return FALSE } // ray is parallel to the polygon
  
      float ray_scalar = dotProduct(vectorSub(Poly.P[0], Ray.R1), plane_normal)
  
      Answer = vectorAdd(Ray.R1, scalarMult(ray_scalar, Ray.R2))
  
      // verify that the point falls inside the polygon
  
      point test_line = vectorSub(Answer, Poly.P[0])
      point test_axis = crossProduct(plane_normal, test_line)
  
      bool point_is_inside = FALSE
  
      point test_point = vectorSub(Poly.P[1], Answer)
      bool prev_point_ahead = (dotProduct(test_line, test_point) > 0)
      bool prev_point_above = (dotProduct(test_axis, test_point) > 0)
  
      bool this_point_ahead
      bool this_point_above
  
      int index = 2;
      while (index < Poly.count)
      {
          test_point = vectorSub(Poly.P[index], Answer)
          this_point_ahead = (dotProduct(test_line, test_point) > 0)
  
          if (prev_point_ahead OR this_point_ahead)
          {
              this_point_above = (dotProduct(test_axis, test_point) > 0)
  
              if (prev_point_above XOR this_point_above)
              {
                  point_is_inside = !point_is_inside
              }
          }
  
          prev_point_ahead = this_point_ahead
          prev_point_above = this_point_above
          index++
      }
  
      return point_is_inside
  }