ICFP 2008 #3 - The End, but with a bit of ray tracing
Real life has sort of overtaken things and it's clear that I won't be able to get as far as I'd like in the ICFP contest so I've decided to bow out. Given that, I figure that I might as well save some time for those still in the running and impart some of my ray tracing knowledge. So I'm giving away two little routines that ought to be helpful to those still in the running. I've translated them to Python as it's a nice form of executable pseudocode. But first, here's a tiny vector class. It's obviously very incomplete, but it's just enough for the two functions I'll give below:
from math import * class vector: def __init__( self, x, y ): self.x = x self.y = y def __add__( self, other ): return vector( self.x + other.x, self.y + other.y ) def __sub__( self, other ): return vector( self.x - other.x, self.y - other.y ) def __mul__( self, other ): return vector( self.x * other, self.y * other ) def dot( self, other ): return self.x * other.x + self.y * other.y
Now, on to the fun stuff. First, you'll undoubtedly want to be able to answer the question: does a line from here to there cross an obstacle. Here's some code that will answer that:
def intersects_circle( start, end, center, radius ): direction = end - start offset = center - start offsetsq = offset.dot( offset ) radiussq = radius * radius if offsetsq < radiussq: return True midchord = offset.dot( direction ) if midchord < 0.0: return False invdirectionsq = 1.0 / direction.dot( direction ) midchord = midchord * invdirectionsq midchordsq = midchord * midchord halfsq = ( radiussq - offsetsq ) * invdirectionsq + midchordsq if halfsq < 0.0: return False hit = midchord - sqrt( halfsq ) return hit <= 1.0
The first test there is a quick check to see whether the start point is already within the radius of the obstacle. If so, then obviously part of the line segment is in the obstacle. The second test looks to see whether the line's direction even heads in any way towards the obstacle. And the third test determines if any part of the line from the start point towards the end point and on out to infinity touches the circle. If it does, the final test in the last statement is whether it reaches the closer edge of the circle before the end point. If you want to know exactly where it first enters the circle, then that's just start + direction * hit.
If you've got a lot of obstacles to test those rays against, you can build a tree out of them and use that to reduce the search time to be logarithmic. One of the simpler ways to do that is with a bounding volume heirarchy. This is a tree where the obstacles will be at the leaves and each interior node will contain the data for a shape (a bounding volume) that completely encloses its all of its subtrees. To test a ray against this tree you start at the root of the tree and see if your ray intersects the root's bounding volume. If so, you recurse to its children and do the same test against them and so forth. If you've built a good BVH, you can quickly cull away large groups of obstacles that you obviously can't hit.
Building the BVH is a bit tricky. There's an old approach (the MacDonald-Booth algorithm) which can work here. Take the first obstacle and make it the root of the tree. Then for each new obstacle that you find, begin at the root of the tree and expand the bounding volume of that node to include the new obstacle. Examine its children and figure out which child's bounding volume would have the smallest perimeter if you expanded it to add the new obstacle. Then descend to that child, expand its bounding volume and check its children. Eventually, you'll reach a leaf, so then you replace that with a new interior node that holds both the original obstacle at that leaf and the new one. This approach can generate highly inferior trees if you're not careful (especially for something like this the trees may end up fairly linear), but it does have the advantage of being fast and incremental.
But what do you want to use as the bounding volume? Circles can be tricky to expand, but doable. They also tend be very loose in terms of space. Boxes are dead simple deal with and are reasonably tight. If you do decide to use boxes, here's how to check whether any part of a line segment goes through a box defined by two opposing corners:
def intersects_box( start, end, corner1, corner2 ): direction = end - start if direction.x == 0: direction.x = 1e-9 if direction.y == 0: direction.y = 1e-9 inverse = vector( 1.0 / direction.x, 1.0 / direction.y ) xmin = ( corner1.x - start.x ) * inverse.x xmax = ( corner2.x - start.x ) * inverse.x ymin = ( corner1.y - start.y ) * inverse.y ymax = ( corner2.y - start.y ) * inverse.y if xmin > xmax: xmin, xmax = xmax, xmin if ymin > ymax: ymin, ymax = ymax, ymin near = max( 0.0, xmin, ymin ) far = min( 1.0, xmax, ymax ) return near < far
Note that for both these functions, there are some simple precomputations that you can do to speed them up if you'll be calling them a lot on the same line segment, box or circle. You'd probably want to cache the inverse direction in the box test, for example. I leave that to you. Still, as for the algorithms themselves, these are pretty much the state of the art in ray tracing so you'll find them tough to beat.
Hope that saves some of you some time! Good luck!