Walk a line between two points in a 3D voxel space visiting all cells

Question 1

Came up with this, or see: http://jsfiddle.net/wivlaro/mkaWf/6/

function visitAll(gx0, gy0, gz0, gx1, gy1, gz1, visitor) {

    var gx0idx = Math.floor(gx0);
    var gy0idx = Math.floor(gy0);
    var gz0idx = Math.floor(gz0);

    var gx1idx = Math.floor(gx1);
    var gy1idx = Math.floor(gy1);
    var gz1idx = Math.floor(gz1);

    var sx = gx1idx > gx0idx ? 1 : gx1idx < gx0idx ? -1 : 0;
    var sy = gy1idx > gy0idx ? 1 : gy1idx < gy0idx ? -1 : 0;
    var sz = gz1idx > gz0idx ? 1 : gz1idx < gz0idx ? -1 : 0;

    var gx = gx0idx;
    var gy = gy0idx;
    var gz = gz0idx;

    //Planes for each axis that we will next cross
    var gxp = gx0idx + (gx1idx > gx0idx ? 1 : 0);
    var gyp = gy0idx + (gy1idx > gy0idx ? 1 : 0);
    var gzp = gz0idx + (gz1idx > gz0idx ? 1 : 0);

    //Only used for multiplying up the error margins
    var vx = gx1 === gx0 ? 1 : gx1 - gx0;
    var vy = gy1 === gy0 ? 1 : gy1 - gy0;
    var vz = gz1 === gz0 ? 1 : gz1 - gz0;

    //Error is normalized to vx * vy * vz so we only have to multiply up
    var vxvy = vx * vy;
    var vxvz = vx * vz;
    var vyvz = vy * vz;

    //Error from the next plane accumulators, scaled up by vx*vy*vz
    // gx0 + vx * rx === gxp
    // vx * rx === gxp - gx0
    // rx === (gxp - gx0) / vx
    var errx = (gxp - gx0) * vyvz;
    var erry = (gyp - gy0) * vxvz;
    var errz = (gzp - gz0) * vxvy;

    var derrx = sx * vyvz;
    var derry = sy * vxvz;
    var derrz = sz * vxvy;

    do {
        visitor(gx, gy, gz);

        if (gx === gx1idx && gy === gy1idx && gz === gz1idx) break;

        //Which plane do we cross first?
        var xr = Math.abs(errx);
        var yr = Math.abs(erry);
        var zr = Math.abs(errz);

        if (sx !== 0 && (sy === 0 || xr < yr) && (sz === 0 || xr < zr)) {
            gx += sx;
            errx += derrx;
        }
        else if (sy !== 0 && (sz === 0 || yr < zr)) {
            gy += sy;
            erry += derry;
        }
        else if (sz !== 0) {
            gz += sz;
            errz += derrz;
        }

    } while (true);
}

Question 2

As far as I remember the original Bresenham algorithm assumes that movement along diagonals is allowed, in your case is makes sense to disallow it.

But the main idea is the same - for every voxel answer the question "what's next?"

Every voxel has 6 faces each leading to a different neighbour. Just check centre of which voxel is closer to the line than others. That's the next voxel.

Note: this assumes that voxel has the same size along every axis, if that's not the case, you should calculate modified distance (every component should be divided by voxel size along corresponding axis)

Question 3

I think 3d Bresenham is the way to go, just tweaked a bit. As a first pass at the problem, proceed as Bresenham, but be suspicious when you're about to take a step, or you've just taken a step, as these are the places where the line could pass through extra cells.

For simplicity, let's assume that z is dominant, meaning that z increments every step. The 3d Bresenham question is: "when do we increment/decrement in x or y?" The answer is when accumulated error in x reaches .5, or when the error in y does, or both.

For your case, I think you need to have a secondary threshold that uses slopeY = deltaY/deltaZto decide if the line is about to cross into a neighboring cell. If stepZ is the change in z along the line for each pixel, then a test like error > .5 - slopeY/stepZ should tell you to get cells on both sides of the line in y. A similar test tells you if you have to get the extra cell in x. If you have to get the extra cell in both x AND y, then you have to get the cell diagonal to the Bresenham cell as well.

If you've detected that you added a cell in y before the increment, you won't add a cell after. If you haven't added a y cell before, you will have to after, unless you happened to pass through a cell corner. How you handle that depends on your use-case.

These are my thoughts on the issue, I haven't tested anything, but something like it should work.

Question 4

Here is a public link to a recent port of my voxel ray from C++ into javascript:

https://github.com/jeremykentbgross/EmpathicCivGameEngine/blob/master/engine/public/scripts/Ray2D.js

Note: the port is currently in 2D on a quadtree (instead of 3D on an octtree), but only because one dimension is commented out for my 2D javascript engine. It works fine in my 3D C++ engine (where I ported it from) so if you uncomment the Z axis lines it will work. The file also has a lot of inline comments on how the math works.

You should also reference the RayTracer2D.js (in the same directory) which uses the ray to find all intersected objects and their intersection points in the order they are hit.

For reference the quad tree structure it is tracing through is also in the same folder: QuadTree.js

Note that you could also ray trace lower LOD's simply by limiting how deep you traverse into the tree during the trace.

Hope that helps.

Question 5

https://code.activestate.com/recipes/578112-bresenhams-line-algorithm-in-n-dimensions/

Here is a numpy implementation for N-D bresenham line drawing in case someone came across this thread from googling 'bresenham 3d python'.