Optimize triangle intersections (~12% speedup)

- Precompute AB, AC and AB^AC in members of Triangle, to avoid
  calculating them each time in intersectTriangle.

- Reorder exit checks a bit, to enable quicker culling of
  triangles we don't intersect.

- Precompute 1/Dcr to avoid divisions.

src/mesh.cpp

@@ -104,10 +104,8 @@ bool Mesh::intersectTriangle(const Ray& ray, const Triangle& T, IntersectionInfo
 	// backface culling?
 	if (backfaceCulling && dot(ray.dir, T.gnormal) > 0) return false;
 	const Vector& A = vertices[T.v[0]];
-	const Vector& B = vertices[T.v[1]];
-	const Vector& C = vertices[T.v[2]];
 
-	if (T.intersect(ray, A, B, C, info.dist, lambda2, lambda3)) {
+	if (T.intersectFast(ray, A, info.dist, lambda2, lambda3)) {
 		info.geom = this;
 		info.ip = ray.start + ray.dir * info.dist;
 		if (faceted || normals.empty()) {
@@ -279,6 +277,9 @@ void Mesh::prepareTriangles()
 		Vector AB = B - A;
 		Vector AC = C - A;
 		t.gnormal = AB ^ AC;
+		t.AB = AB;
+		t.AC = AC;
+		t.ABcrossAC = t.gnormal;
 		t.gnormal.normalize();
 
 		if (!uvs.empty() && !normals.empty()) {

src/triangle.cpp

@@ -39,15 +39,49 @@ bool Triangle::intersect(Ray ray, const Vector& A, const Vector& B, const Vector
 	double Dcr = det(AB, AC, D);
 
 	if (fabs(Dcr) < 1e-12) return false;
+	double rDcr = 1 / Dcr;
 
 	Vector H = ray.start - A;
 
-	double lambda2 = det(H, AC, D) / Dcr;
-	double lambda3 = det(AB, H, D) / Dcr;
-	double gamma   = det(AB, AC, H) / Dcr;
+	double gamma   = det(AB, AC, H) * rDcr;
+	if (gamma < 0 || gamma > minDist) return false;
+
+	double lambda2 = det(H, AC, D) * rDcr;
+	if (lambda2 < 0 || lambda2 > 1) return false;
+
+	double lambda3 = det(AB, H, D) * rDcr;
+
+	if (lambda3 < 0 || lambda3 > 1) return false;
+
+	double lambda1 = 1 - (lambda2 + lambda3);
+	if (lambda1 < 0) return false;
+
+	minDist = gamma;
+	l2 = lambda2;
+	l3 = lambda3;
+	return true;
+}
+
+bool Triangle::intersectFast(Ray ray, const Vector& A, double& minDist, double& l2, double& l3) const
+{
+	Vector D = -ray.dir;
+
+	double Dcr = dot(ABcrossAC, D);
+
+	if (fabs(Dcr) < 1e-12) return false;
+	double rDcr = 1 / Dcr;
+
+	Vector H = ray.start - A;
 
+	double gamma   = dot(ABcrossAC, H) * rDcr;
 	if (gamma < 0 || gamma > minDist) return false;
-	if (lambda2 < 0 || lambda2 > 1 || lambda3 < 0 || lambda3 > 1) return false;
+
+	double lambda2 = det(H, AC, D) * rDcr;
+	if (lambda2 < 0 || lambda2 > 1) return false;
+
+	double lambda3 = det(AB, H, D) * rDcr;
+
+	if (lambda3 < 0 || lambda3 > 1) return false;
 
 	double lambda1 = 1 - (lambda2 + lambda3);
 	if (lambda1 < 0) return false;

src/triangle.h

@@ -32,7 +32,9 @@ struct Triangle {
 	int t[3]; //!< holds indices to the three texture coordinates of the triangle (indexes in the `uvs' array)
 	Vector gnormal; //!< The geometric normal of the mesh (AB ^ AC, normalized)
 	Vector dNdx, dNdy; //!< tangent and binormal vectors for this triangle
+	Vector AB, AC, ABcrossAC;
 
 	static bool intersect(Ray ray, const Vector& A, const Vector& B, const Vector& C, double& minDist,
 						  double& l2, double& l3);
+	bool intersectFast(Ray ray, const Vector& A, double& minDist, double& l2, double& l3) const;
 };