Here is a full C99/C11 implementation of a simple voxel renderer:
// SPDX-License-Identifier: CC0-1.0
//
// Compile using e.g.
// gcc -Wall -Wextra -O2 this.c -lm -o this
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <errno.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;
/* 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;
}
/*
* Voxel map
*/
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)
{
const float eps = 0.5f / 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 > eps) {
xnext = float4_add4(pos, float4_mul1(dir, 1.0f - posf.x));
xdelta = float4_div1(dir, dir.x);
} else
if (dir.x < -eps) {
xnext = float4_sub4(pos, float4_mul1(dir, 0.0f - posf.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 > eps) {
ynext = float4_add4(pos, float4_mul1(dir, 1.0f - posf.y));
ydelta = float4_div1(dir, dir.y);
} else
if (dir.y < -eps) {
ynext = float4_sub4(pos, float4_mul1(dir, 0.0f - posf.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 > eps) {
znext = float4_add4(pos, float4_mul1(dir, 1.0f - posf.z));
zdelta = float4_div1(dir, dir.z);
} else
if (dir.z < -eps) {
znext = float4_sub4(pos, float4_mul1(dir, 0.0f - posf.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);
}
if (xnext.w < 0.0f || xdelta.w < 0.0f) fprintf(stderr, "Warning: xnext.w = %.6f, xdelta.w = %.6f\n", xnext.w, xdelta.w);
if (ynext.w < 0.0f || ydelta.w < 0.0f) fprintf(stderr, "Warning: ynext.w = %.6f, ydelta.w = %.6f\n", ynext.w, ydelta.w);
if (znext.w < 0.0f || zdelta.w < 0.0f) fprintf(stderr, "Warning: znext.w = %.6f, zdelta.w = %.6f\n", znext.w, zdelta.w);
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)
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.
*/
/* 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 (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);
/* 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;
}
}
}
}
}
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;
}
renderPPM(ppm, pgm, 320, 160, Float3(68.0f, 32.0f, 68.0f), Float3(4.0f,4.0f,4.0f), Float3(4.0f, -4.0f, 0.0f), 256.0f);
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;
}
Here is a full C99/C11 implementation of a simple voxel renderer.
Most of the code is supportFirst, let's define vectors.h, for HLSL-like float3/, float4/, and int4 types, so that the code is easier to port. (I don't have anything that can compile and run HLSL, myself.)support:
Compile the program, then run it giving two file names as command-line parameters. The first one will be a NetPBM PPM image of the scene, and the second one a NetPBM PGM depth map:
// 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 */
The voxel space here is periodic, so although thereHere is only one sphere in the voxel mapexample renderer, the rendering has multiple copies of it.render.c:
Note how the nearest sphere – in the lower right corner of the rendered image – has green faces on it, while they obviously shouldn't exist on that side of the sphere? They're better visible if you increase the image size, too. What is going on?
// 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;
}
SinceThere were issues in xnext, ynext, znext, xdelta, ydelta, and zdelta calculations, causing each voxel cell has only three faces (and not sixto be displayed as intersecting planes, like ainstead of cube has), we see the faces we expectwalls. It is important to see onlynote that when viewing in a negative direction (where unit ray direction vector hascorrect, they all nonpositivehave nonnegative w components). When we view in other directions, we do not see the surface face firstand that xnext.x, but the faces between nonempty voxel cells!xdelta.x, ynext.y, ydelta.y, znext.z, and zdelta.z are all integers.
In 2D,The code that makes the equivalent 2×2 block has extra protruding linesvoxel 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;
In 3DFloating-point vector modulo can also be used, those extra faces are like vanes or plates standing up on the outsidebut be careful of the surface!semantics: usually, the remainder has the same sign as the argument, but here we want the positive remainder.
Simply putIf 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 we wantyou 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 drawcreate a nice voxel boxcube, weyou need to populateset four voxel cells: $$\begin{array}{ccl} \text{Cell} & \text{Value} & \text{Description} \\ \hline (0,0,0) & 7 & \text{x, y, and z faces} \\ (1,0,0) & 1 & \text{x face} \\ (0,1,0) & 2 & \text{y face} \\ (1,1,1) & 4 & \text{z face} \\ \end{array}$$ (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.
Alternatively, you can adjustWhat the cell index calculation, by preparing an array of eight vectors (constant for each rayabove adjustment does, depends only on the direction; each component is either 0 or -1), that are added to the current position just before the cell index calculation. (Essentially, each intersection type has their own adjustment.) That way we can shift the x faces, the y faces, and the z faces separately from the others for the voxel grid, depending on the ray direction;cell walls so that they always face the eye. So, each cell thenvoxel effectively has SIXsix faces.
This is what the results look like, if compiled with SPHERE defined (-DSPHERE):
