#define EPSILON 0.000001

int
intersect_triangle_barycentric(
   double orig[3], double dir[3],
   double vert0[3], double vert1[3], double vert2[3],
   double planeq[4], int i1, int i2,
   double *t, double *alpha, double *beta)
{
   double	dot, dot2;
   double	point[2];
   double	u0, v0, u1, v1, u2, v2;

   /* is ray parallel to plane? */
   dot = dir[0] * planeq[0] + dir[1] * planeq[1] + dir[2] * planeq[2];
   if (dot < EPSILON && dot > -EPSILON)		/* or use culling check */
      return 0;

   /* find distance to plane and intersection point */
   dot2 = orig[0] * planeq[0] +
      orig[1] * planeq[1] + orig[2] * planeq[2];
   *t = -(planeq[3] + dot2 ) / dot;
   point[0] = orig[i1] + dir[i1] * *t;
   point[1] = orig[i2] + dir[i2] * *t;

   /* begin barycentric intersection algorithm */
   u0 = point[0] - vert0[i1];
   v0 = point[1] - vert0[i2];
   u1 = vert1[i1] - vert0[i1];
   u2 = vert2[i1] - vert0[i1];
   v1 = vert1[i2] - vert0[i2];
   v2 = vert2[i2] - vert0[i2];

   /* calculate and compare barycentric coordinates */
   if (u1 == 0) {		/* uncommon case */
      *beta = u0 / u2;
      if (*beta < 0 || *beta > 1)
	       return 0;
      *alpha = (v0 - *beta * v2) / v1;
   }
   else {			/* common case, used for this analysis */
      *beta = (v0 * u1 - u0 * v1) / (v2 * u1 - u2 * v1);
      if (*beta < 0 || *beta > 1)
	       return 0;
      *alpha = (u0 - *beta * u2) / u1;
   }

   if (*alpha < 0 || (*alpha + *beta) > 1.0)
      return 0;

   return 1;
}