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):
