+// http://www.cs.lth.se/home/Tomas_Akenine_Moller/raytri/raytri.c
+// http://little3d.free.fr/ressources/jgt%20Fast,%20Minumum%20Storage%20Ray-Triangle%20Intersection.htm
+// http://www.acm.org/jgt/papers/MollerTrumbore97/
+
+/* Ray-Triangle Intersection Test Routines */
+/* Different optimizations of my and Ben Trumbore's */
+/* code from journals of graphics tools (JGT) */
+/* http://www.acm.org/jgt/ */
+/* by Tomas Moller, May 2000 */
+
+/* code rewritten to do tests on the sign of the determinant */
+/* the division is at the end in the code */
+// cosmetics change by H.J :
+// make u & v locals since we don't use them, use sg functions
+static bool intersect_triangle(const double orig[3], const double dir[3],
+ const double vert0[3], const double vert1[3], const double vert2[3],
+ double *t)
+{
+ double u, v;
+ double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+
+ const SGDfloat eps = 1e-4;
+
+ /* find vectors for two edges sharing vert0 */
+ sgdSubVec3(edge1, vert1, vert0);
+ sgdSubVec3(edge2, vert2, vert0);
+
+ /* begin calculating determinant - also used to calculate U parameter */
+ sgdVectorProductVec3(pvec, dir, edge2);
+
+ /* if determinant is near zero, ray lies in plane of triangle */
+ double det = sgdScalarProductVec3(edge1, pvec);
+
+ if (det > eps)
+ {
+ /* calculate distance from vert0 to ray origin */
+ sgdSubVec3(tvec, orig, vert0);
+
+ /* calculate U parameter and test bounds */
+ u = sgdScalarProductVec3(tvec, pvec);
+ if (u < 0.0 || u > det)
+ return false;
+
+ /* prepare to test V parameter */
+ sgdVectorProductVec3(qvec, tvec, edge1);
+
+ /* calculate V parameter and test bounds */
+ v = sgdScalarProductVec3(dir, qvec);
+ if (v < 0.0 || u + v > det)
+ return false;
+
+ }
+ else if(det < -eps)
+ {
+ /* calculate distance from vert0 to ray origin */
+ sgdSubVec3(tvec, orig, vert0);
+
+ /* calculate U parameter and test bounds */
+ u = sgdScalarProductVec3(tvec, pvec);
+ if (u > 0.0 || u < det)
+ return false;
+
+ /* prepare to test V parameter */
+ sgdVectorProductVec3(qvec, tvec, edge1);
+
+ /* calculate V parameter and test bounds */
+ v = sgdScalarProductVec3(dir, qvec) ;
+ if (v > 0.0 || u + v < det)
+ return false;
+ }
+ else return false; /* ray is parallell to the plane of the triangle */
+
+ /* calculate t, ray intersects triangle */
+ *t = sgdScalarProductVec3(edge2, qvec) / det;
+
+ return true;
+}
+
+