The DDA Algorithm, explained interactively
I've written a number of voxel raytracers, and all of them (all the good ones, at least) use the Digital Differential Analyzer Algorithm for raycasting. I've only ever "implemented" this algorithm once, by copying a reference somewhere. Since then, I've used and re-used that code. But I have something to admit: I've never really understood how it works. I've seen it explained well, explained poorly, explained on video, and still, I just don't get it. It's one of those geometric algorithms that's hard to wrap your head around, at least for me. But recently, I ran into some issues with a raytracer, and I decided to try to fully understand it. And I made this blog post, for others like me.
All code in this article is editable. Examples will update as you modify them.
DDA, in 2D, is an algorithm that iterates over all of the grid squares intersected by a ray. For those who haven't seen the algorithm, it looks like this:
If you understand the above, have a nice day :) If not, this article has an explanation.
3Blue1Brown teaches us that a good way to to explain something is to derive it from scratch. So how would you derive this line drawing algorithm? First, we've got our arguments: the Ray Origin ro
, and Ray Direction rd
. RD has to be normalized. In this article, we calculate it using an angle, like this:
Let's plot those on a grid.
The first grid space we need to fill in is pretty simple–it's the grid space that ro
is in. We can floor ro to get the exact integer space.
Not super exciting, but an important step.
Now, we have to figure out what the next grid space is. See if you can guess the answer: move the point and the arrow above and guess where the next grid space would be.
At first glance, you might think it has something to do with rdAngle
, like if the angle is facing up, then the grid space above the dot should be filled in, if the angle is facing to the right, the one to the right should be filled in, etc.
But if you , you'll notice that this quickly falls apart.
We need to take into account the distance from the origin point to the different grid squares. If you've spent some time recently with geometry or a graphing calculator, you might already have figured out that this is actually a line intersection problem. If you , you can see that the ro->rd vector intersects with the lines x = 1 and y = 1. (In this article, we are assuming the center of the grid is the origin.) Let's put those intersections on the screen:
The yellow dot is the intersection between ro->rd and y = 1, and the green dot is the intersection between ro->rd and x = 1.
Now, let's call the distance between ro and the yellow dot lengthY.
Similarly, we can call the distance between ro and green dot lengthX.
Now see that when lengthX is less than lengthY, the raycasting line intersects x = 1 before it intersects y = 1, so the next grid space is the one to the right. Similarly, when lengthY is less than lengthX, the raycasting line intersects y = 1 before it intersects x = 1, so the next grid space is the one above.
Great! We've taken our first steps. Before we move to three squares (can you imagine?!), let's talk about that intersection code.
Importantly, we don't actually need to calculate the intersection points. We only need lengthX and lengthY, which are both distances. For those, we can consult Pythagoras. If we imagine the yellow line above as the diagonal of a right triangle, we can calculate its length using its legs, drawn below as pink and blue.
Also note that we're going to be using ro = {0,0} for a little while, for simplicity's sake :')
The length of the pink line is just 1.
We can calculate the color of the blue line with some basic algebra. We want to find the point where the yellow line (presented in Slope-intercept form mx+b) intersects y = 1.
Using the pythagorean theorem, we can then calculate the length of the line.
Similarly,
We can verify that these values are correct by multiplying them by rd and drawing the resulting vector on the screen:
And now all we have to do is use the comparison trick we had earlier.
Great! Now let's get to space three. If we move to the right, we want to start measuring the distance to the next line, x = 2. The distance to this line is double the distance to x = 1. On the other hand, if we move up, we want to measure the distance to y = 2, which is just double the distance to y = 1. Let's try doubling our distances depending on the branch.
At each step, we want to increase the distToX if we move right, and we want to increase the distToY if we move up. And we want to increase yDist by the distance between rows along rd, and xDist by the distance between columns along rd. The distance between columns is just the distance from {0,0} to x = 1, which is the value we calculated earlier. Same for the distance between rows.
We still have to deal with cases where ro != {0,0}. The initial distToX & distToY change depending on where the ro is inside the grid square. Let's call that value roFract, since it represents the fractional portion of ro.
When y = 0, the intial distToY is distBetweenColumns, as shown above. When y = 0.5, the intial distToY is half distBetweenColumns. In the general case, to calculate the initial distToY, we can just multiply the distBetweenColumns by 1 - roFract.y. We can do the same for distToX.
We're almost done–we just have to make sure that we're calculating the right distance when rd.x < 0 or rd.y < 0. When rd.x < 0, we're measuring the distance to x = 0 instead of x = 1, so distToX = roFract.x * distBetweenColumns: we don't have to subtract roFract from one. Similarly, when rd.y < 0, distToY = roFract.y * distBetweenRows.
Also, if rd.y < 0, instead of increasing gridSpace.y, we need to decrease gridSpace.y. We generalize for both cases by simply adding sign(rd.y) instead of increasing by the constant 1. We do the same for rd.x and gridSpace.x.
And that's it! Some modifications are used in the real world:
- distBetweenRows & distBetweenColumns are combined into a single variable, sometimes called rayUnitStepSize.
- sign(rd.x) & sign(rd.y) are combined into a single variable, sometimes called step.
- To avoid the comparison in the loop, you can use a boolean vector.