What is the use of epsilon in the Möller–Trumbore intersection algorithm?

The Möller–Trumbore intersection algorithm compares the determinant with an epsilon. What is this epsilon and how is its value chosen?

The epsilon is there to counter floating-point accuracy problems. Note that the code checks that the algorithm is between $-\epsilon$ and $\epsilon$; that is, close to zero. If the determinant is zero, then the ray lies in the plane of the triangle and can be safely said to not intersect it (as the alternative, an edge-on intersection, should be caught by neighbouring triangles in the mesh).

The Cramer's rule computation of the intersection point then divides by the determinant, and here's where the floating-point accuracy comes in. If the determinant is merely very close to zero, not exactly equal to it, then its reciprocal will be very large. We multiply the other determinants by this reciprocal, so small floating-point errors in their values will be blown up into significance and make the difference between an intersection and a non-intersection.

I haven't tested it, but I suspect that leaving out the epsilon will cause intersections on the silhouette of a mesh (i.e. with triangles close to parallel with the ray direction) to become very nondeterministic. This would result in nearby rays either missing or hitting the triangles unpredictably and giving you speckling.

The value of the epsilon depends on the floating-point precision you're using. It should be small enough that only very grazing intersections are culled, but large enough that it avoids the speckling problem. It can be determined by experiment, or by just eyeballing it. The implementation you linked to uses $\epsilon = 10^{-7}$, which is a reasonable choice for single-precision arithmetic; you may want to use $10^{-9}$ or $10^{-10}$ for double-precision computations.