Voxel raycasting algorithm not working correctly

Question

I am using Unity with a compute shader to render to a texture. So far I have been checking points along the rays in units of 1 just for testing, so I know it works. Now I am simply trying to write a function that increments the ray to stop at the next voxel face, basically. And I think I understand the algorithm but my implementation is not working correctly. I included pictures of the problem below.

EDIT: Got it working, it is just really slow at 1080p. I can't really do more than a view distance of 100 without it being unplayable. If I used something like OpenGL or Vulcan, would it be significantly faster? I thought this would be fast enough because while it is Unity it's basically directX11 when it's all on the compute shader. I was really hoping this would work. And advice on how to speed it up would be much appreciated.

I later optimized the crap out of it and it is a tiny bit faster but still slow so I think it has something to do with Unity.

Here is the full compute shader:

#pragma kernel CSMain

RWTexture2D<float4> Result; // the actual array of pixels the player sees
float width; // in pixels
float height;

StructuredBuffer<int> voxelMaterials; // for now just getting a flat voxel array
int voxelBufferRowSize;
StructuredBuffer<float3> rayDirections; // I'm now actually using it as points instead of directions
float maxRayDistance;

float3 playerCameraPosition; // relative to the voxelData, ie the first voxel's bottom, back, left corner position, no negative coordinates
float3 playerWorldForward;
float3 playerWorldRight;
float3 playerWorldUp;

[numthreads(8, 8, 1)]
void CSMain(uint3 id : SV_DispatchThreadID)
{
    Result[id.xy] = float4(0, 0, 0, 0); // setting the pixel to black by default
    float3 pointHolder = playerCameraPosition; // initializing the first point to the player's position
    float3 p = rayDirections[id.x + (id.y * width)]; // vector transformation getting the world space directions of the rays relative to the player
    float3 u1 = p.x * playerWorldRight;
    float3 u2 = p.y * playerWorldUp;
    float3 u3 = p.z * playerWorldForward;
    float3 direction = u1 + u2 + u3; // the direction to that point

    float distanceTraveled = 0;
    while (distanceTraveled < maxRayDistance) 
    {           
        // finding the distances to the next voxel on all axises
        float3 distancesXYZ = { 1000, 1000, 1000 };
        if (direction.x > 0) {
            distancesXYZ.x = (ceil(pointHolder.x) - pointHolder.x) / direction.x;
            if (distancesXYZ.x == 0) {
                distancesXYZ.x = 1 / direction.x;
            }
        }
        else if (direction.x < 0) {
            distancesXYZ.x = (floor(pointHolder.x) - pointHolder.x) / direction.x;
            if (distancesXYZ.x == 0) {
                distancesXYZ.x = 1 / abs(direction.x);
            }
        }
        if (direction.y > 0) {
            distancesXYZ.y = (ceil(pointHolder.y) - pointHolder.y) / direction.y;
            if (distancesXYZ.y == 0) {
                distancesXYZ.y = 1 / direction.y;
            }
        }
        else if (direction.y < 0) {
            distancesXYZ.y = (floor(pointHolder.y) - pointHolder.y) / direction.y;
            if (distancesXYZ.y == 0) {
                distancesXYZ.y = 1 / abs(direction.y);
            }
        }
        if (direction.z > 0) {
            distancesXYZ.z = (ceil(pointHolder.z) - pointHolder.z) / direction.z;
            if (distancesXYZ.z == 0) {
                distancesXYZ.z = 1 / direction.z;
            }
        }
        else if (direction.z < 0) {
            distancesXYZ.z = (floor(pointHolder.z) - pointHolder.z) / direction.z;
            if (distancesXYZ.z == 0) {
                distancesXYZ.z = 1 / abs(direction.z);
            }
        }

        int face = 0; // 1 = x, 2 = y, 3 = z
        // finding smallest distance along the direction to the next voxel
        float smallestDistance = 1000;
        if (distancesXYZ.x < smallestDistance) {
            smallestDistance = distancesXYZ.x;
            face = 1;
        }
        if (distancesXYZ.y < smallestDistance) {
            smallestDistance = distancesXYZ.y;
            face = 2;
        }
        if (distancesXYZ.z < smallestDistance) {
            smallestDistance = distancesXYZ.z;
            face = 3;
        }
        if (face == 0) {
            break;
        }

        pointHolder += direction * smallestDistance;
        distanceTraveled += smallestDistance;

        // convert the point into a voxel index and check if a voxel exists there

        int3 voxelIndexXYZ = { -1,-1,-1 }; // the integer coordinates within the buffer

        if (face == 1) {
            if (direction.x >= 0) {
                voxelIndexXYZ.x = floor(pointHolder.x);
            }
            else {
                if (voxelIndexXYZ.x == 0) {
                    break; // breaking and leaving the pixel black because this is out of bounds of the voxel buffer
                }
                voxelIndexXYZ.x = ceil(pointHolder.x - 1);
            }
            voxelIndexXYZ.y = floor(pointHolder.y);
            voxelIndexXYZ.z = floor(pointHolder.z);
        }
        else if (face == 2) {
            if (direction.y >= 0) {
                voxelIndexXYZ.y = floor(pointHolder.y);
            }
            else {
                if (voxelIndexXYZ.y == 0) {
                    break; // breaking and leaving the pixel black because this is out of bounds of the voxel buffer
                }
                voxelIndexXYZ.y = ceil(pointHolder.y - 1);
            }
            voxelIndexXYZ.x = floor(pointHolder.x);
            voxelIndexXYZ.z = floor(pointHolder.z);
        }
        else if (face == 3) {
            if (direction.z >= 0) {
                voxelIndexXYZ.z = floor(pointHolder.z);
            }
            else {
                if (voxelIndexXYZ.z == 0) {
                    break; // breaking and leaving the pixel black because this is out of bounds of the voxel buffer
                }
                voxelIndexXYZ.z = ceil(pointHolder.z - 1);
            }
            voxelIndexXYZ.y = floor(pointHolder.y);
            voxelIndexXYZ.x = floor(pointHolder.x);
        }
        else {
            break;
        }

        //check if voxelIndexXYZ is within bounds of the voxel buffer
        if (voxelIndexXYZ.x < voxelBufferRowSize && voxelIndexXYZ.x >= 0 && 
            voxelIndexXYZ.y < voxelBufferRowSize && voxelIndexXYZ.y >= 0 && 
            voxelIndexXYZ.z < voxelBufferRowSize && voxelIndexXYZ.z >= 0)
        {
            int voxelIndex = voxelIndexXYZ.x + (voxelIndexXYZ.z * voxelBufferRowSize) + (voxelIndexXYZ.y * (voxelBufferRowSize * voxelBufferRowSize)); // the voxel index in the flat array

            if (voxelMaterials[voxelIndex] == 1) { // if the voxel has a material ID of 1, using integers to represent materials, 0 is empty
                Result[id.xy] = float4((distanceTraveled / maxRayDistance) * 2, (float)voxelIndex / (voxelBufferRowSize * voxelBufferRowSize * voxelBufferRowSize), pointHolder.z, 0); // giving it a crazy color
                break;
            }
        }
    }
}

That code produces this:

and in case you were curious this is what it looks like marching the ray by 1 unit. It was taken at a different time without the middle pole of voxels and a smaller view distance which changes the color.

The bottom is two layers deep. Marching the ray by one actually looks cool except it violently wiggles around when the player moves.

Answer 1

(Apologies, I'm not familiar with Unity/HLSL at all, so my syntax may be wonky.)

Let's restart, with

float4  eye;     // Player camera position, .w = 0
float4  proj0;   // Corner of the projection plane, .w = 0
float4  projX;   // Projection plane X axis, .w = 0
float4  projY;   // Projection plane Y axis, .w = 0

so that picture plane pixel uint3 id is at proj0 + id.x*projX + id.y*projY in the voxel coordinates. (We'll use the w component for the distance, saving a number of operations; that's why these are four-component vectors.)

The ray starting point position and unit direction vector unitdir are then

float4  position = proj0 + id.x*projX + id.y*projY;
float4  unitdir = normalize(start - eye);

We can save some operations if we set

position.w = 0.0f;
unitdir.w = 1.0f;

here, but that's the sort of optimization you'll need to do yourself here.

Next, we need the distances to the very first faces the ray intersects (as measured along the ray itself):

float xstart, ystart, zstart;

if (unitdir.x > 0) {
    xstart = (1.0f + floor(unitdir.x) - unitdir.x) / unitdir.x;
} else
if (unitdir.x < 0) {
    xstart = (floor(unitdir.x) - unitdir.x) / unitdir.x;
} else {
    xstart = +INF;
}

if (unitdir.y > 0) {
    ystart = (1.0f + floor(unitdir.y) - unitdir.y) / unitdir.y;
} else
if (unitdir.y < 0) {
    ystart = (floor(unitdir.y) - unitdir.y) / unitdir.y;
} else {
    ystart = +INF;
}

if (unitdir.y > 0) {
    zstart = (1.0f + floor(unitdir.z) - unitdir.z) / unitdir.z;
} else
if (unitdir.y < 0) {
    zstart = (floor(unitdir.z) - unitdir.z) / unitdir.z;
} else {
    zstart = +INF;
}

The part in parentheses is the distance along the axis, so dividing by the same axis unit vector component we get the length along the ray.

Note that xstart, ystart, and zstart are all positive.

Again, the above can be computed much more efficiently. Now, we initialize the first positions where the ray intersects the first face of each type,

float4  xnext = position + xstart * unitdir;
float4  ynext = position + ystart * unitdir;
float4  znext = position + zstart * unitdir;
xnext.w = xstart;  // This can be omitted, since position.w = 0 and unitdir.w = 1
ynext.w = ystart;  // This can be omitted, since position.w = 0 and unitdir.w = 1
znext.w = zstart;  // This can be omitted, since position.w = 0 and unitdir.w = 1

and when we "use" one of those vectors, we add the corresponding delta,

float4  xdelta = (1.0f / unitdir.x) * unitdir;
float4  ydelta = (1.0f / unitdir.y) * unitdir;
float4  zdelta = (1.0f / unitdir.z) * unitdir;
xdelta.w = 1.0f / unitdir.x;  // This can be omitted, since unitdir.w = 1
ydelta.w = 1.0f / unitdir.y;  // This can be omitted, since unitdir.w = 1
zdelta.w = 1.0f / unitdir.z;  // This can be omitted, since unitdir.w = 1

The actual ray traversal loop is then

int4   cell;  /* Cell integer coordinates */
float4 cellf; /* Fractional cell coordinates */
uint   face;  /* Face, edge, or vertex */

while (1) {
    if (position.w >= maxRayDistance) {
        face = 0;
        break;
    }

    if (xnext.w < ynext.w && xnext.w < znext.w) {
        position = xnext;
        xnext += xdelta;
        face = 1;  // X face
    } else
    if (ynext.w < xnext.w && ynext.w < znext.w) {
        position = ynext;
        ynext += ydelta;
        face = 2;  // Y face
    } else
    if (znext.w < xnext.w && znext.w < ynext.w) {
        position = znext;
        znext += zdelta;
        face = 4;  // Z face
    } else
    if (xnext.w == ynext.w && xnext.w < znext.w) {
        position = xnext;
        xnext += xdelta;
        ynext += ydelta;
        face = 3;  // XY edge
    } else
    if (xnext.w == znext.w && xnext.w < ynext.w) {
        position = xnext;
        xnext += xdelta;
        znext += zdelta;
        face = 5;  // XZ edge
    } else
    if (ynext.w == znext.w && ynext.w < xnext.w) {
        position = ynext;
        ynext += ydelta;
        znext += zdelta;
        face = 6;  // YZ edge
    } else
    if (xnext.w == ynext.w && xnext.w == znext.w) {
        position = xnext;
        xnext += xdelta;
        ynext += ydelta;
        znext += zdelta;
        face = 7;  // XYZ vertex
    } else {
        // .w are all +INF
        face = 0;
        break;
    }

    cellf = modf(position, cell);

    // We intersected with 'face' at cell
    //     cell.x, cell.y, cell.z
    // with fractional coordinates
    //     cellf.x, cellf.y, cellf.z
    // The exact coordinates are
    //     position.x = cell.x + cellf.x
    //     position.y = cell.y + cellf.y
    //     position.z = cell.z + cellf.z
    // but note that cellf can be negative.

    // TODO: Examine the voxel buffer.
    // Note: position.w is the distance from the projection plane
    //       to the face intersection point.  This can be useful
    //       for depth buffers.
}

It is up to you if you use face at all, but all eight cases do need to be handled. (It is unlikely the eighth case ever occurs, but it's good to be thorough, methinks.)

Answer 2

Here is a full C99/C11 implementation of a simple voxel renderer.

First, let's define vectors.h, for float3, float4, and int4 support:

// SPDX-License-Identifier: CC0-1.0
//
#ifndef  VECTORS_H
#define  VECTORS_H
#include <math.h>

/*
 *  float3 support
*/

typedef struct {
    float   x;
    float   y;
    float   z;
} float3;

static inline float3  Float3(const float x, const float y, const float z)
{
    const float3  result = { x, y, z };
    return result;
}

static inline float3  float3_sub3(const float3 a, const float3 b)
{
    const float3  result = { a.x - b.x, a.y - b.y, a.z - b.z };
    return result;
}

static inline float3  float3_add3(const float3 a, const float3 b)
{
    const float3  result = { a.x + b.x, a.y + b.y, a.z + b.z };
    return result;
}

static inline float  float3_length(const float3 a)
{
    return sqrtf(a.x*a.x + a.y*a.y + a.z*a.z);
}

static inline float  float3_dot(const float3 a, const float3 b)
{
    return a.x*b.x + a.y*b.y + a.z*b.z;
}

static inline float3  float3_cross(const float3 a, const float3 b)
{
    const float3  result = { a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x };
    return result;
}

static inline float3  float3_mul1(const float3 a, const float b)
{
    const float3  result = { a.x*b, a.y*b, a.z*b };
    return result;
}

static inline float3  float3_scale_to_length(const float3 a, const float b)
{
    const float n = float3_length(a);
    const float3  result = { a.x*b/n, a.y*b/n, a.z*b/n };
    return result;
}

/*
 *  float4 support
*/

typedef struct {
    float   x;
    float   y;
    float   z;
    float   w;
} float4;

static inline float4  Float4(const float x, const float y, const float z, const float w)
{
    const float4  result = { x, y, z, w };
    return result;
}

/* float4_length(a) == || a ||, Euclidean length of vector a */
static inline float  float4_length(const float4 a)
{
    return sqrtf(a.x*a.x + a.y*a.y + a.z*a.z + a.w*a.w);
}

/* float4_normalize(a) == a / || a || */
static inline float4  float4_normalize(const float4 a)
{
    const float   n = float4_length(a);
    const float4  result = { a.x/n, a.y/n, a.z/n, a.w/n };
    return result;
}

/* float4_sign(a): Component-wise sign: -1.0, 0.0, or +1.0 */
static inline float4  float4_sign(const float4 a)
{
    const float4  result = { (a.x < 0.0f) ? -1.0f : (a.x > 0.0f) ? +1.0f : 0.0f,
                             (a.y < 0.0f) ? -1.0f : (a.y > 0.0f) ? +1.0f : 0.0f,
                             (a.z < 0.0f) ? -1.0f : (a.z > 0.0f) ? +1.0f : 0.0f,
                             (a.w < 0.0f) ? -1.0f : (a.w > 0.0f) ? +1.0f : 0.0f };
    return result;
}

/* float4_max(a, b): Component-wise maximum */
static inline float4  float4_max4(const float4 a, const float4 b)
{
    const float4  result = { (a.x >= b.x) ? a.x : b.x,
                             (a.y >= b.y) ? a.y : b.y,
                             (a.z >= b.z) ? a.z : b.z,
                             (a.w >= b.w) ? a.w : b.w };
    return result;
}

/* float4_min(a, b): Component-wise minimum */
static inline float4  float4_min4(const float4 a, const float4 b)
{
    const float4  result = { (a.x <= b.x) ? a.x : b.x,
                             (a.y <= b.y) ? a.y : b.y,
                             (a.z <= b.z) ? a.z : b.z,
                             (a.w <= b.w) ? a.w : b.w };
    return result;
}

/* float4_add(a, b) == a + b */
static inline float4  float4_add4(const float4 a, const float4 b)
{
    const float4  result = { a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w };
    return result;
}

/* float4_sub4(a, b) == a - b */
static inline float4  float4_sub4(const float4 a, const float4 b)
{
    const float4  result = { a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w };
    return result;
}

/* float4_floor4(a): Round each component towards negative infinity */
static inline float4  float4_floor4(const float4 a)
{
    const float4  result = { floorf(a.x), floorf(a.y), floorf(a.z), floorf(a.w) };
    return result;
}

/* float4_mul1(a, b) = { a.x*b, a.y*b, a.z*b, a.w*b } */
static inline float4  float4_mul1(const float4 a, const float b)
{
    const float4  result = { a.x*b, a.y*b, a.z*b, a.w*b };
    return result;
}

/* float4_div1(a, b) = { a.x/b, a.y/b, a.z/b, a.w/b } */
static inline float4  float4_div1(const float4 a, const float b)
{
    const float4  result = { a.x/b, a.y/b, a.z/b, a.w/b };
    return result;
}

/* float4_div4(a, b) = { a.x/b.x, a.y/b.y, a.z/b.x, a.w/b.w } */
static inline float4  float4_div4(const float4 a, const float4 b)
{
    const float4  result = { a.x/b.x, a.y/b.y, a.z/b.z, a.w/b.w };
    return result;
}

/*
 * int4 support
*/

typedef struct {
    int     x;
    int     y;
    int     z;
    int     w;
} int4;

static inline int4  Int4(const int x, const int y, const int z, const int w)
{
    const int4  result = { x, y, z, w };
    return result;
}

static inline int4  int4_add4(const int4 a, const int4 b)
{
    const int4  result = { a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w };
    return result;
}

/* float4_int4(a): Component-wise cast to int. */
static inline int4  float4_int4(const float4 a)
{
    const int4  result = { (int)floorf(a.x), (int)floorf(a.y), (int)floorf(a.z), (int)floorf(a.w) };
    return result;
}

#endif /* VECTORS_H */

Here is the example renderer, render.c:

// SPDX-License-Identifier: CC0-1.0
//
// Compile using e.g.
//      gcc -DSPHERE -Wall -Wextra -O2 render.c -lm -o render
// or a simpler test with
//      gcc -Wall -Wextra -O2 render.c -lm -o render


#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include "vectors.h"


float4           voxel_rgba[256][8];    /* Look-up table for voxel faces; all components [0..1]. x is red, y is green, z is blue, w is opacity (0=transparent, 1=opaque). */
int4             voxel_size;
unsigned char   *voxel_cell = NULL;
size_t           voxel_xstride = 0;     /* Typically 1 */
size_t           voxel_ystride = 0;     /* Typically voxel_size.x */
size_t           voxel_zstride = 0;     /* Typically voxel_size.x * voxel_size.y */


/* Trace one voxel ray starting at (projection plane) pos, with eye/camera at eye.
 * Return the color and distance { .x=red, .y=green, .z=blue, .w=distance }
*/
float4 voxel_ray(float4 eye, float4 pos, float maxdist)
{
    /* On input, eye and pos are really 3-component vectors; we need the fourth one to be zero,
       so that it won't affect the unit direction length vector below. */
    eye.w = 0.0f;
    pos.w = 0.0f;

    /* Ray unit direction vector */
    float4  dir = float4_normalize(float4_sub4(pos, eye));
    /* We rely on the .w component to track length, so set that one now. */
    dir.w = 1.0f;

    float4  posf = float4_sub4(pos, float4_floor4(pos));
    /* Note: 0 <= posf.x < 1,
             0 <= posf.y < 1,
             0 <= posf.z < 1. */

    /* Find first intersections with a voxel cell wall (*next),
       and the delta to the consecutive following intersections */
    float4  xnext,  ynext,  znext;
    float4  xdelta, ydelta, zdelta;

    if (dir.x > 0.0f) {
        xnext  = float4_add4(pos, float4_mul1(dir, (1.0f - posf.x) / dir.x));
        xdelta = float4_div1(dir, dir.x);
    } else
    if (dir.x < 0.0f) {
        xnext  = float4_add4(pos, float4_mul1(dir, -posf.x / dir.x));
        xdelta = float4_div1(dir, -dir.x);
    } else {
        xnext  = Float4(0.0f, 0.0f, 0.0f, maxdist);
        xdelta = Float4(0.0f, 0.0f, 0.0f, 0.0f);
    }

    if (dir.y > 0.0f) {
        ynext  = float4_add4(pos, float4_mul1(dir, (1.0f - posf.y) / dir.y));
        ydelta = float4_div1(dir, dir.y);
    } else
    if (dir.y < 0.0f) {
        ynext  = float4_add4(pos, float4_mul1(dir, -posf.y / dir.y));
        ydelta = float4_div1(dir, -dir.y);
    } else {
        ynext  = Float4(0.0f, 0.0f, 0.0f, maxdist);
        ydelta = Float4(0.0f, 0.0f, 0.0f, 0.0f);
    }

    if (dir.z > 0.0f) {
        znext  = float4_add4(pos, float4_mul1(dir, (1.0f - posf.z) / dir.z));
        zdelta = float4_div1(dir, dir.z);
    } else
    if (dir.z < 0.0f) {
        znext  = float4_add4(pos, float4_mul1(dir, -posf.z / dir.z));
        zdelta = float4_div1(dir, -dir.z);
    } else {
        znext  = Float4(0.0f, 0.0f, 0.0f, maxdist);
        zdelta = Float4(0.0f, 0.0f, 0.0f, 0.0f);
    }

    float4  color = { 0.0f, 0.0f, 0.0f, 0.0f };  /* Transparent! */

    while (1) {
        unsigned char  intersection = 0;    /* 1:x, 2:y, 4:z */

        if (pos.w >= maxdist) {
            pos.w = maxdist;
            break;
        }
        if (color.w >= 1.0f)
            break;

        /* Pick the closest next step first. */
        pos = xnext;
        if (pos.w > ynext.w) {
            pos = ynext;
        }
        if (pos.w > znext.w) {
            pos = znext;
        }

        /* Update intersection and prepare for the next step. */
        if (pos.w >= xnext.w) {
            intersection |= 1;
            xnext = float4_add4(xnext, xdelta);
        }
        if (pos.w >= ynext.w) {
            intersection |= 2;
            ynext = float4_add4(ynext, ydelta);
        }
        if (pos.w >= znext.w) {
            intersection |= 4;
            znext = float4_add4(znext, zdelta);
        }

        /* If pos.w == INF, we have intersection = 0. */
        if (!intersection) {
            pos.w = maxdist;
            break;
        }

        /* Position within the wraparound voxel space. */
        float4  temp = float4_floor4(pos);
        int4    posi = float4_int4(temp);
        /* We could use the fractional positive sub-voxel coordinates posf,
              posf = float4_sub4(pos, temp);
           where 0 <= posf.x < 1, 0 <= posf.y < 1, 0 <= posf.z < 1
           and if (intersection & 1), posf.x = 0 (except for rounding errors),
               if (intersection & 2), posf.y = 0 (except for rounding errors),
               if (intersection & 4), posf.z = 0 (except for rounding errors),
           for interpolation etc.
        */

        /* Adjust cell coordinates so that each cell always defines an outer wall. */
        if ((intersection & 1) && (dir.x < 0.0f)) --posi.x;
        if ((intersection & 2) && (dir.y < 0.0f)) --posi.y;
        if ((intersection & 4) && (dir.z < 0.0f)) --posi.z;

        /* Ensure posi is within the positive voxel space. */
        posi.x = posi.x % voxel_size.x; if (posi.x < 0) posi.x += voxel_size.x;
        posi.y = posi.y % voxel_size.y; if (posi.y < 0) posi.y += voxel_size.y;
        posi.z = posi.z % voxel_size.z; if (posi.z < 0) posi.z += voxel_size.z;

        /* Look up the voxel cell properties for this intersection. */
        float4  c = voxel_rgba[voxel_cell[ (size_t)posi.x * voxel_xstride
                                         + (size_t)posi.y * voxel_ystride
                                         + (size_t)posi.z * voxel_zstride ]][ intersection ];

        if (c.w >= 1.0f) {
            /* Opaque; good, ray ends here. Blend 'c' behind 'color'. */
            color = float4_add4(color, float4_mul1(c, 1.0f - color.w));
            break;
        } else
        if (c.w > 0.0f) {
            /* Blend color 'color' *behind* color 'c'. */
            color = float4_add4(color, float4_mul1(c, 1.0f - color.w));
        }
    }
    color.w = pos.w;
    return color;
}

void renderPPM(FILE *outppm, FILE *outpgm, int width, int height, const float3 eye, const float3 forward, const float3 right, const float maxdist)
{
    /* Assume 'right' is perpendicular to 'forward'.  'up' is perpendicular to both, with length (height/width) times that of 'right'. */
    const float3  up = float3_scale_to_length(float3_cross(forward, right), float3_length(right) * (float)height / (float)width);

    /* Image plane corner, rowstart = eye + forward - right + up */
    float3  rowstart = float3_add3(float3_sub3(float3_add3(eye, forward), right), up);

    /* Delta vectors per pixel for the image plane */
    const float3  dx = float3_mul1(right, 2.0f / (float)width);
    const float3  dy = float3_mul1(up,   -2.0f / (float)height);

    if (outppm) fprintf(outppm, "P6\n%d %d 255\n", width, height);
    if (outpgm) fprintf(outpgm, "P5\n%d %d 255\n", width, height);

    float  mind = +3.0f*maxdist;
    float  maxd = -3.0f*maxdist;

    for (int y = 0; y < height; y++, rowstart = float3_add3(rowstart, dy)) {
        float3  pos = rowstart;
        for (int x = 0; x < width; x++, pos = float3_add3(pos, dx)) {
            const float4  c = voxel_ray( Float4(eye.x, eye.y, eye.z, 0.0f),
                                         Float4(pos.x, pos.y, pos.z, 0.0f), maxdist);
            const int  r8 = (c.x <= 0.0f) ? 0 : (c.x < 1.0f) ? (int)(0.5f + 255.0f * c.x) : 255;
            const int  g8 = (c.y <= 0.0f) ? 0 : (c.y < 1.0f) ? (int)(0.5f + 255.0f * c.y) : 255;
            const int  b8 = (c.z <= 0.0f) ? 0 : (c.z < 1.0f) ? (int)(0.5f + 255.0f * c.z) : 255;
            const int  d8 = (c.w <= 0.0f) ? 0 : (c.w < maxdist) ? (int)(0.5f + 255.0f * c.w / maxdist) : 255;
            if (c.w < maxdist) {
                if (mind > c.w) mind = c.w;
                if (maxd < c.w) maxd = c.w;
            }

            if (outppm) {
                fputc(r8, outppm);
                fputc(g8, outppm);
                fputc(b8, outppm);
            }
            if (outpgm) {
                fputc(d8, outpgm);
            }
        }
        fprintf(stderr, "\rRow %d of %d completed.", y + 1, height);
        fflush(stderr);
    }

    if (outppm) fflush(outppm);
    if (outpgm) fflush(outpgm);
    fprintf(stderr, "\rRendering complete. Distances varied between %.6f and %.6f.\n", mind, maxd);
    fflush(stderr);
}

int main(int argc, char *argv[])
{
    FILE *ppm, *pgm;

    if (argc != 3 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        const char *arg0 = (argc > 0 && argv && argv[0] && argv[0][0]) ? argv[0] : "(this)";
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s [ -h | --help ]\n", arg0);
        fprintf(stderr, "       %s OUT.ppm DEPTH.pgm\n", arg0);
        fprintf(stderr, "\n");
        return EXIT_FAILURE;
    }

    voxel_size.x = 64;
    voxel_size.y = 64;
    voxel_size.z = 64;

    voxel_xstride = 1;
    voxel_ystride = (size_t)voxel_size.x;
    voxel_zstride = voxel_ystride * (size_t)voxel_size.y;
    const size_t  size = voxel_zstride * (size_t)voxel_size.y;

    voxel_cell = (unsigned char *)malloc(size);
    if (!voxel_cell) {
        fprintf(stderr, "Not enough memory for a %d x %d x %d voxel map.\n", voxel_size.x, voxel_size.y, voxel_size.z);
        return EXIT_FAILURE;
    }
    memset(voxel_cell, 0, size);

    /* Make all cell values transparent, */
    for (int i = 0; i < 256; i++) {
        for (int k = 0; k < 8; k++) {
            voxel_rgba[i][k] = Float4(0.0f, 0.0f, 0.0f, 0.0f);
        }
    }

    /* Cell type 1 faces are blue, red, and green; edges and vertices their mix. */
    voxel_rgba[1][1] = Float4(0.0f, 0.0f, 1.0f, 1.0f);
    voxel_rgba[1][2] = Float4(0.0f, 1.0f, 0.0f, 1.0f);
    voxel_rgba[1][4] = Float4(1.0f, 0.0f, 0.0f, 1.0f);
    voxel_rgba[1][3] = Float4(0.0f, 0.8f, 0.8f, 1.0f);
    voxel_rgba[1][5] = Float4(0.8f, 0.0f, 0.8f, 1.0f);
    voxel_rgba[1][6] = Float4(0.8f, 0.8f, 0.0f, 1.0f);
    voxel_rgba[1][7] = Float4(0.6f, 0.6f, 0.6f, 1.0f);

#ifdef SPHERE
    /* Create a shell at the center, minradius 10, maxradius 12 */
    {
        const int  cx = 32;
        const int  cy = 32;
        const int  cz = 32;
        const int  rrmin = 10*10;
        const int  rrmax = 12*12;

        for (int z = 0; z < voxel_size.z; z++) {
            const int  zz = (z-cz)*(z-cz);
            for (int y = 0; y < voxel_size.y; y++) {
                const int  zzyy = zz + (y-cy)*(y-cy);
                for (int x = 0; x < voxel_size.x; x++) {
                    const int  dd = zzyy + (x-cx)*(x-cx);
                    if (dd >= rrmin && dd < rrmax) {
                        voxel_cell[(size_t)x * voxel_xstride + (size_t)y * voxel_ystride + (size_t)z * voxel_zstride] = 1;
                    }
                }
            }
        }
    }
#else
    voxel_cell[voxel_xstride + voxel_ystride] = 1;
    voxel_cell[voxel_xstride] = 1;
    voxel_cell[voxel_ystride] = 1;
    voxel_cell[voxel_zstride] = 1;
#endif

    fprintf(stderr, "Constructed a %d x %d x %d voxel map.\n", voxel_size.x, voxel_size.y, voxel_size.z);

    ppm = fopen(argv[1], "wb");
    if (!ppm) {
        fprintf(stderr, "%s: %s.\n", argv[1], strerror(errno));
        return EXIT_FAILURE;
    }
    pgm = fopen(argv[2], "wb");
    if (!pgm) {
        fprintf(stderr, "%s: %s.\n", argv[2], strerror(errno));
        fclose(ppm);
        remove(argv[1]);
        return EXIT_FAILURE;
    }

#ifdef SPHERE
    renderPPM(ppm, pgm, 512, 384, Float3(0.0f, 0.0f, 0.0f), Float3(4.0f, 5.0f, 6.0f), Float3(4.0f, -2.0f, -1.0f), 128.0f);
#else
    renderPPM(ppm, pgm, 512, 512, Float3(-8.0f, -8.0f, -8.0f), Float3(1.0f, 4.0f, 4.0f), Float3(-1.0f, 0.0f, 1.0f), 9.0f);
#endif

    if (fclose(ppm)) {
        fprintf(stderr, "%s: Error closing file.\n", argv[1]);
        fclose(pgm);
        remove(argv[1]);
        remove(argv[2]);
        return EXIT_FAILURE;
    }
    if (fclose(pgm)) {
        fprintf(stderr, "%s: Error closing file.\n", argv[2]);
        remove(argv[1]);
        remove(argv[2]);
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Saved PPM image as '%s', and depth graymap as PGM image '%s'.\n", argv[1], argv[2]);
    return EXIT_SUCCESS;
}

There were issues in xnext, ynext, znext, xdelta, ydelta, and zdelta calculations, causing each voxel cell to be displayed as intersecting planes, instead of cube walls. It is important to note that when correct, they all have nonnegative w components, and that xnext.x, xdelta.x, ynext.y, ydelta.y, znext.z, and zdelta.z are all integers.

The code that makes the voxel space periodic is this:

        /* Position within the wraparound voxel space. */
        float4  temp = float4_floor4(pos);
        int4    posi = float4_int4(temp);

        /* Ensure posi is within the positive voxel space. */
        posi.x = posi.x % voxel_size.x; if (posi.x < 0) posi.x += voxel_size.x;
        posi.y = posi.y % voxel_size.y; if (posi.y < 0) posi.y += voxel_size.y;
        posi.z = posi.z % voxel_size.z; if (posi.z < 0) posi.z += voxel_size.z;

Floating-point vector modulo can also be used, but be careful of the semantics: usually, the remainder has the same sign as the argument, but here we want the positive remainder.

If you don't mind mirroring the voxel space, you can use something like posi = floor(modf(abs(pos), size)) instead in HLSL. Then, positive voxel space is periodic, but is mirrored with respect to each axial plane (x=0, y=0, and z=0).

To move each voxel cell wall towards the eye, I added the section

        /* Adjust cell coordinates so that each cell always defines an outer wall. */
        if ((intersection & 1) && (dir.x < 0.0f)) --posi.x;
        if ((intersection & 2) && (dir.y < 0.0f)) --posi.y;
        if ((intersection & 4) && (dir.z < 0.0f)) --posi.z;

You can replace this with a single vector addition or subtraction, if you create an array of 8 vectors (one for each possible intersection type).

Without this adjustment, looking towards negative x would show the y and z faces of the cell before the x face. This is because normally, the cells have only three faces, and to create a cube, you need to set four voxel cells: (0,0,0) has all three faces, (+1,0,0) has the opposite x face, (0,+1,0) has the opposite y face, and (0,0,+1) has the opposite z face.

What the above adjustment does, is shift the voxel cell walls so that they always face the eye. So, each voxel effectively has six faces.

This is what the results look like, if compiled with SPHERE defined (-DSPHERE):

Answer 3

For optimizing voxel tracing, I can think of two different approaches:

1. Multi-level mapping to skip larger empty volumes

2. Ganged initial tracing, so that a single tracer handles the common parts of many rays

---

Multi-level mapping means that you have an additional voxel map for blocks of k×k×k voxels. At minimum, each block is just a flag whether that block is empty or not. Point (x, y, z) belongs to block (floor(x/k), floor(y/k), floor(z/k)), so you wish to pick a k whose inverse you can express exactly with a floating-point number; for example, a power of two (k = 2, 4, 8, 16, 32, 64, ...).

Whenever the inner tracing function notices that the voxel block changes, it checks if the voxel block is empty. If it is empty, the inner tracing function can skip the voxel block entirely (using an inner loop that advances normally, but never examines the voxel data, just loops until it exits the current voxel block).

I don't think one can expect massive speedups with this, but it can be useful if the voxel block is actually either empty, or a reference to a look-up voxel map (which are all k×k×k voxels in size). This means that large voxel maps with lots of empty regions can fit in less memory.

Another option, making the tracing non-voxel-like, is to use large voxel cells (but small tracing depths, say on the order of 50 or so maximum). When one finds the nearest next position, the current and that next position forms a line through the current voxel; between two adjacent or opposite faces. If the contents of that voxel are represented using geometric primitives et cetera, one can check if the ray intersects anything within the voxel or not. Essentially, the voxel map is then just a way to group the possibly intersecting objects into subsets.

---

Ganged initial tracing means you have two different types of tracers.

The initial tracers take four rays, each with their own starting point, direction; and the screen coordinates for the set. For example, let S00, S01, S10, and S11 be their starting points (float3 or float4), D00, D01, D10, and D11 be their directions (float3 or float4), x0 be the screen x coordinate for S00 and S10, x1 be the screen x coordinate for S01 and S11, y0 be the screen y coordinate for S00 and S01, and y1 be the screen y coordinate for S10 and S11.

As long as the four rays stay in the same voxel, they are advanced by the same tracer. Whenever they diverge, the set is split according to their screen x and y coordinates. Whenever the set refers to only four connected pixels (i.e. x1 = x0 + 1 and y1 = y0 + 1), they are split into individual rays, and traced as described in my other answer as an individual ray.

The "trick" is that one cannot just use the average of the corresponding ray to calculate the added rays when splitting: they must be weighted per the coordinates.

For example, if we have x0 = 4 and x1 = 7, we can choose the split point x = 5 or x = 6. Whichever we pick, the new ray starting location is weighed by (x1 - x) / (x1 - x0) for the one corresponding to x0, and by (x - x0) / (x1 - x0) for the one corresponding to x1.

(Note, it may look like the weights are swapped, but they're not. Linear interpolation from x0 to x1 using 0 <= w <= 1 is (1-w)*x0 + w*x1, you see. Here, 1 - (x - x0) / (x1 - x0) = (x1 - x0 - x + x0) / (x1 - x0) = (x1 - x) / (x1 - x0).)

If the ganged (initial) tracer finds that the voxel cell – remember, it only progresses the rays if they stay within the same voxel cell – it has to "paint" a rectangular region of the display plane, from (x0,y0) to (x1,y1).

Now, in my other answer, the distance to the first voxel face the ray intersects with is calculated directly; that point saved in the xnext, ynext, and znext vectors. The ganged tracing (for a set of four rays) can be continued as long as the rays stay in the same or neighboring voxels. So, one option would be to gang rays that have the same $y$ coordinate, as long as their first intersections are in the same voxel. Assuming the rendered image is much wider than it is tall (which is typical for humans, we have a field of vision that is wider than it is tall), this could be a suitable compromise between complexity and speedup obtained.

(Then, one type of tracer takes two rays. As soon as the dot product of the ray current position integer coordinates exceeds 1, they are split. Otherwise, the rays are advanced as a pair. Splitting uses the split point x coordinate as the weight as described above. The new split ray may not be at a voxel face, though.)

I think I shall test the horizontally ganged tracer, and compare its speed to the individual one, and report here later if it was worth it.