I settled on the half-planes solution. This algorithm defines a plane at each edge of the triangle by crossing the triangle's normal and the edge direction. To check if a triangle is inside the triangle, you simply check to see if a point (on the plane of the triangle) is on the same side of all of those planes.

This costs a little extra storage space, but it's a very fast test.

So. How does the code pan out? Well, pointInsideTriangle is simple enough:

distanceToPlane :: Shape -> Vector -> Float

distanceToPlane (Plane !norm !dist) !pos = pos `dot3` norm + dist

distanceToPlane _ _ = undefined

pointInsideTriangle :: Triangle -> Position -> Bool

pointInsideTriangle !tri !point = foldr (&&) True $ map (\pln -> (distanceToPlane pln point) >= 0) (halfPlanes tri)

This is best read right-to-left. The map applies a lambda function to each of the half planes. This function simply returns a bool, telling us if a point is on the correct side of a triangle. This yields a list of bools. We need to turn that into a single result; this is done using a fold to reduce many values to a single Bool. I've deliberately used a right-fold so that the fold lazily consumes only the necessary parts of the list. It can early-out.

So there we have it. One line of Haskell (with a couple of helper functions) to do the test.

The rest of the Haskell code for a ray-triangle intersection is pretty straightforward:

intersectRayTriangle :: Ray -> Object -> Triangle -> Bool -> Maybe (Float, Triangle)

intersectRayTriangle !ray !obj !triangle !doubleSided

| doubleSided == False && (direction ray) `dot3` (normal $ plane triangle) > 0 = Nothing

| otherwise = case shapeClosestIntersect (plane triangle) ray obj of

Nothing -> Nothing

Just (dist', _) -> if pointInsideTriangle triangle (pointAlongRay ray dist')

then Just (dist', triangle)

else Nothing

Here we have a simple guard condition to permit backface culling. If the triangle is not double-sided, we do a backface culling test and possibly reject it. If not, we fall through to the default case and intersect the ray against the plane, and then test that intersection point for containment in the triangle.

Intersecting against a list is an interesting case. A simple version would just map against the list, and then attempt to find the closest intersection out of the resulting list. However, if we use tail recursion, we can add a parameter to hold our current "state" and therefore maintain a current-closest intersection. This eliminates a second search-the-list step, and also permits an early-reject optimisation:

intersectRayTriangleList :: [Triangle] -> Int -> Maybe (Float, Int) -> Ray -> Object -> Maybe (Float, Int)

intersectRayTriangleList !(x:xs) !index !currentResult !currentRay !obj = intersectRayTriangleList xs (index + 1) newResult newRay obj

where

(newRay, newResult) = case intersectRayTriangle currentRay obj x False of

Nothing -> (currentRay, currentResult)

Just (dist, _) -> (shortenRay currentRay dist, Just (dist, index))

intersectRayTriangleList [] _ !currentResult _ _ = currentResult

The pattern of "searching a list of x and return the closest" occurs frequently in raytracing. This is a potential candidate to be factored out into a re-usable function, or possibly even a monad.

The resulting code is also quite efficient. I'm now running at less than a minute for my test scene of the Cornell Box with 8x8 distributed raytracing. This is 30% quicker than my previously MÃ¶ller-Trumbore ray-triangle intersection test.

Half-planes are interesting!

ReplyDeleteI think you should implement pointInsideTriangle with "all". And I'd probably do intersectRayTriangle with the Maybe monad, using "guard".

One more thing... doubleSided == False? Dude!

Using 'all' saved 4 seconds :-)

ReplyDeleteThanks for the pointer, I guess I should get to know some of these standard functions a little better!

Could this be the same Tom Hammersley from wolvo? Haskell... Surely not.. ;-)

ReplyDeleteI guess I'm a little out-of-date with this stuff but why not barrycentric coords?

Yep, it's me alright.. the very same.

ReplyDeleteBarycentric is good and all, but the half-planes were a nicer (more natural) match to FP programming.

FWIW I use barycentric for the interpolation of normals etc :-)