Bresenham(ish) line, skipping, raycasting.

Haven't been up to much lately - just trying to have some holidays. I have been tentatively poking at some ideas for raycasting. Not that I have any particular idea in mind but it might be something to try on the parallella.

Part of that is to step through the cells looking for ray hits, which is where bresenham's line drawing algorithm comes in. I'm working toward a heirarchical data structure with wider branching than a typical quad/oct tree and so because of that I need to be able to step in arbitrary amounts and not just per-cell. Fairly cheaply.

I don't really like the description on wikipedia and although there are other explanations I ended up deriving the maths myself which allows me to handle arbitrary steps. Unfortunately this does need a division but I can't see how that can be avoided in the general case.

If one takes a simple expression for a line:

  y = mx + c
  m = dy / dx
    = (y1 - y0) / (x1 - x0)

And implements it directly using real arithmetic assuming x increments by a fixed amount:

 dx = x1-x0;
 dy = y1-y0;
 for (x=0; x<dx; x+=step) {
   x = x0 + x;
   y = round((real)x * dy / dx) + y0;
   plot(x, y);
 }

Since it steps incrementally the multiply/division can be replaced by addition:

  incy = (real)step * dy / dx;
  ry = 0;
  for (x=0; x<dx; x+=step) {
    x = x0 + x;
    y = round(ry) + y0;
    plot(x, y);
    ry += incy;
  }

But ... floating point aren't 'real' numbers they are just approximations, so this leads to rounding errors not to mention the rounding and type conversions needed.

Converting the equation to using rational numbers solves this problem.

  a = dy * step / dx;        // whole part of dy / dx
  b = dy * step - a * dx;    // remainder
  ry = 0;                    // whole part of y
  fy = 0;                    // fractional part of y
  for (x=0; x<dx; x+=step) {
    x = x0 + x;
    y = ry + y0;
    plot(x, y);
    ry += a;                 // increment by whole part
    fy += b;                 // tally remainder
    if (fy > dx) {           // is it a whole fraction yet?
      fy -= dx;              // carry it over to the whole part
      ry += 1;
    }
  }

Where y = ry + fy / dx (whole part plus fractional part).

For rendering into the centre of pixels one shifts everything by dx/2 but this cannot be accurately represented using integers. So instead the upper test on the fractional overflow is shifted by dx/2 by multiplying both sides by 2 which leads to:;

    if (fy * 2 > dx) {
      fy -= dx;
      ry += 1;
    }

This suits my requirements because step can be changed at any point during the algorithm (together with a newly calculated a and b) to step over any arbitrary number of locations. It only needs one division per step size.

The black outline is calculated using floating point/rounded, and the filled-colour dots are calculated using the final algorithm. The green are calculated by stepping 8 pixels at a time and the blue are calculated by starting at each green location with the same (y = ry + fy/dx) and running the equation with single stepping values for 3 pixels.

(i did this kinda shit on amiga coding 20 years ago; it's just been so long and outside of what i normally hack on that i always need to do a refresher when I revisit things like this).

Unfortunately for raycasting things are a bit more involved so I will have to see whether I can be bothered getting any further than this.

Update: So yeah after all that I realised it isn't much use for raycasting anyway since you start with floats from the view transform - converting this to integer coordinates will already be an approximation. And since hardware float multiply is just as fast as float add may as well just use the first equation with floats and pre-calculate dy/dx so any rounding errors will be consistent and not accumulate.

So here's a plot of a 2d ray through a 16-way (4x4) tree. Cyan X is when it hits a vertical cell wall, purple X is when it hits a horizontal cell wall, and black circle is when it hits a non-transparent leaf cell. The black boxes are the nodes and leaf elements in the tree, and the pink boxes are those which are touched by the ray. I just filled the tree with some random circles.

The code turned out to be pretty simple and consists of just two steps:

First the location is converted to an integer (using truncation) and that is used to lookup into the tree. Because each level of the tree is 4x4, 2 bits are used for both x and y at each level so the lookup at each level is a simple bit shift and mask. If the given cell is in a non-transparent leaf node it is a hit.

Then it calculates the amount it has to step in X to get beyond the current cell bounds: this will simply be MIN(cell.x + cell.width - loc.x, (cell.y + cell.height - loc.y) * slopex). The cell size is calculated implicitly based on the node depth.

The 2d case isn't all that interesting to me but it's easier to visualise/test.

Hmm, on the other hand this still has some pesky edge cases and rounding issues to deal with so perhaps an integer approach is still useful.