Trouble with ray-triangle intersection optimization

Dear Group

As a personal project of mine, I wanted to integrate the M�ller-Trumbore (see white paper at algorithm for ray-triangle intersection into Radiance. The reason why I chose this algorithm is because it was simple to implement and for front facing, early division can be avoided. I have read the white paper and the math all checks out. The algorithm has been out for some time and has been implemented in numerous other raytracers. I have not heard any problems with the algorithm but in my own implementation neither the algorithm for front or back facing seem to work properly. I have check the code, which was adapted from the source code given in the paper, but it still doesn�t work like it should. I have thrown my hands up because I do not know what I have done wrong. I have considered the Segura-Feito method (implemented in o_mesh.c) but the problem here is that this method isn�t ray-triangle intersection but segment-triangle intersection. The only way that I know of how to calculate the ray segment length is find the ray length via ray-plane intersection. Unfortunately, this requires division which negates the primary advantage of the Segura-Feito method. Here is the modified version of o_face.c that I have been working on for M�ller-Trumbore. Any help�. Please?



#ifndef lint
static const char RCSid[] = "$Id: o_face.c,v 2.4 2003/03/11 17:08:55 greg Exp $";
* o_face.c - compute ray intersection with faces.
#include "copyright.h"
#include "ray.h"
#include "face.h"
#include "standard.h"
#include "object.h"
#include "fvect.h"
o_face(o, r) /* compute intersection with polygonal face */
register RAY *r;
  double rdot; /* direction . normal */
  double t; /* distance to intersection */
  FVECT pisect; /* intersection point */
  register FACE *f; /* face record */
  register int i;

  f = getface(o);
   * First, we find the distance to the plane containing the
   * face. If this distance is less than zero or greater
   * than a previous intersection, we return. Otherwise,
   * we determine whether in fact the ray intersects the
   * face. The ray intersects the face if the
   * point of intersection with the plane of the face
   * is inside the face.
  /* compute dist. to plane */
  if(f->nv == 3)
      if (intersect_triangle(r,f) == 0)

      rdot = -DOT(r->rdir, f->norm);
      if (rdot <= FTINY && rdot >= -FTINY) /* ray parallels plane */
        t = FHUGE;
        t = (DOT(r->rorg, f->norm) - f->offset) / rdot;
      if (t <= FTINY || t >= r->rot) /* not good enough */
      /* compute intersection */
      for (i = 0; i < 3; i++)
        pisect[i] = r->rorg[i] + r->rdir[i]*t;
      if (!inface(pisect, f)) /* ray intersects face? */
      r->rot = t;
      VCOPY(r->rop, pisect);
      r->rod = rdot;
  r->ro = o;
  VCOPY(r->ron, f->norm);
  r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
  r->uv[0] = r->uv[1] = 0.0;
  r->rox = NULL;
  return(1); /* hit */
/* The following function is adapted from the M�ller-Trumbore algorithm ***/
intersect_triangle(r, f)
FACE *f;
register RAY *r;
  double *u, *v, *t;
  double det,inv_det, d;
  int i;
  double rdot;
  FVECT tvec, pvec, qvec, edge1, edge2;

  VSUB(edge1, VERTEX(f,1), VERTEX(f,0));
  VSUB(edge2, VERTEX(f,2), VERTEX(f,0));

  /* begin calculating determinant - also used to calculate U parameter */
  VCROSS(pvec, r->rdir, edge2);
  /* if determinant is near zero, ray lies in plane of triangle */
  det = DOT(edge1, pvec);
  if (det > -FTINY && det < FTINY)
    return 0;
  inv_det = 1.0 / det;
  /* calculate distance from vert0 to ray origin */
  VSUB(tvec, r->rorg, VERTEX(f,0));
  /* calculate U parameter and test bounds */
  *u = DOT(tvec, pvec)* inv_det;
  if (*u < 0.0 || *u > 1.0)
    return 0;
  /* prepare to test V parameter */
  VCROSS(qvec, tvec, edge1);
  /* calculate V parameter and test bounds */
  *v = DOT(r->rdir, qvec) * inv_det;
  if (*v < 0.0 || *u + *v > 1.0)
    return 0;
  /* calculate t, scale parameters, ray intersects triangle */

  *t = DOT(edge2, qvec) * inv_det;

  if (*t <= FTINY || *t >= r->rot) /* not good enough */
  r->rot = *t;
  rdot = -DOT(r->rdir, f->norm);
  r->rod = rdot;
  for (i = 0; i < 3; i++)
    r->rop[i] = r->rorg[i] + r->rdir[i]* r->rot;
  return 1;


MSN Movies - Trailers, showtimes, DVD's, and the latest news from Hollywood!