cpu_voxel_ray_tracing_in_370_lines_of_c89

/** Me: Can we stop and get Teardown? 
 * Mom: we have Teardown at home.
 * Teardown at home:
 * 
 * Note that there are many bugs (and obviously missing features (eg. reflections)). Compiled using TCC (this should comply with ANSI C89) **/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#define PI 3.14159
#define FOV 90.0f*(PI/180.0f)

#define MAX(a,b) (((a)>(b))?(a):(b))
#define MIN(a,b) (((a)<(b))?(a):(b))

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

float mag(float_vec vec)
{
	return sqrt(fabs(vec.x) + fabs(vec.y) + fabs(vec.z));
}

float_vec norm_vec(float_vec unnorm)
{
	float scale = 1.0f / mag(unnorm);

	float_vec normed = { scale * unnorm.x, scale * unnorm.y, scale * unnorm.z };

	return normed;
}

float_vec mul_vec_scl(float scl, float_vec sec)
{
	float_vec new = { scl * sec.x, scl * sec.y, scl * sec.z };
	return new;
}

float_vec add(float_vec fir, float_vec sec)
{
	float_vec new = { fir.x + sec.x, fir.y + sec.y, fir.z + sec.z };
	return new;
}

float dot(float_vec fir, float_vec sec)
{
	return (fir.x * sec.x) + (fir.y * sec.y) + (fir.z * sec.z);
}

float phong(float ia, float id, float is, float ka, float kd, float ks, float alpha, float_vec norm, float_vec raydir, float_vec light_dir)
{
	float il = 0;

	float_vec rm = add(mul_vec_scl(2.0f*dot(light_dir, norm), norm), mul_vec_scl(-1.0f, light_dir));
	il += ka*ia;

	il += (id*kd*dot(light_dir, norm));
	il += (is*ks*(float)pow((double)(MAX(0.0f, dot(rm, raydir))), (double)alpha));
	return il;
}

typedef struct
{
	float_vec pos;
	float_vec dir;
} ray;

typedef struct
{
	int mat;
	float_vec hit_pos;
	float_vec norm;
} rayhit;

/** Helper function for dda **/
void sign(float_vec val, int *out[3])
{
	int signedarr[3] = { ((val.x > 0) ? 1 : -1), ((val.y > 0) ? 1 : -1), ((val.z > 0) ? 1 : -1) };
	memcpy(out, signedarr, sizeof(int)*3);
}

/** Helper function for dda **/
void maxcalc(ray curr_ray, int sign[3], float *out[3])
{
	int curr_next[3] = { ((sign[0] > 0) ? (ceil(curr_ray.pos.x)) : (floor(curr_ray.pos.x))), ((sign[1] > 0) ? (ceil(curr_ray.pos.y)) : (floor(curr_ray.pos.y))), ((sign[2] > 0) ? (ceil(curr_ray.pos.z)) : (floor(curr_ray.pos.z))) };
	float maxs[3] = { (fabs(((float)curr_next[0]) - curr_ray.pos.x) / curr_ray.dir.x), (fabs(((float)curr_next[1]) - curr_ray.pos.y) / curr_ray.dir.y), (fabs(((float)curr_next[2]) - curr_ray.pos.z) / curr_ray.dir.z) };
	memcpy(out, maxs, sizeof(float)*3);
}

/** Helper function for dda **/
void deltacalc(float_vec val, int sign[3], float *out[3])
{
	float deltas[3] = { ((sign[0] > 0) ? (1 / val.x) : (1 / -val.x)), ((sign[1] > 0) ? (1 / val.y) : (1 / -val.y)), ((sign[2] > 0) ? (1 / val.z) : (1 / -val.z)) };
	memcpy(out, deltas, sizeof(float)*3);
}

/** See: "A Fast Voxel Traversal Algorithm for Ray Tracing" by Amanatides and Andrew Woo (http://www.cse.yorku.ca/~amana/research/grid.pdf) **/
rayhit dda(ray curr_ray, int scene[64][64][64])
{
	/** Initialization **/
	rayhit current_hit;

	int int_pos[3] = { round(curr_ray.pos.x), round(curr_ray.pos.y), round(curr_ray.pos.z) };

	int step[3];
	float tmax[3];
	float tdelta[3];

	int collided_mat = 0;
	float len = 1;
	float_vec norm = { 0.0f, 0.0f, 0.0f };
	
	sign(curr_ray.dir, &step);
	maxcalc(curr_ray, step, &tmax);
	deltacalc(curr_ray.dir, step, &tdelta);

	while((collided_mat == 0))
	{
		if (tmax[0] < tmax[1])
		{
			if (tmax[0] < tmax[2])
			{
				int_pos[0] += step[0];

				if ((int_pos[0] > 63) || (int_pos[0] < 0))
				{
					break;
				}
				tmax[0] += tdelta[0];
				len = (float)(tmax[0]);

				float_vec temp = { step[0], 0.0f, 0.0f };
				norm = temp;
			}
			else
			{
				int_pos[2] += step[2];

				if ((int_pos[2] > 63) || (int_pos[2] < 0))
				{
					break;
				}
				
				tmax[2] += tdelta[2];
				len = (float)(tmax[2]);

				float_vec temp = { 0.0f, 0.0f, step[2] };
				norm = temp;
			}
		}
		else
		{
			if (tmax[1] < tmax[2])
			{
				int_pos[1] += step[1];

				if ((int_pos[1] > 63) || (int_pos[1] < 0))
				{
					break;
				}
				
				tmax[1] += tdelta[1];
				len = (float)(tmax[1]);

				float_vec temp = { 0.0f, step[1], 0.0f };
				norm = temp;
			}
			else
			{
				int_pos[2] += step[2];

				if ((int_pos[2] > 63) || (int_pos[2] < 0))
				{
					break;
				}
				
				tmax[2] += tdelta[2];
				len = (float)(tmax[2]);

				float_vec temp = { 0.0f, 0.0f, step[2] };
				norm = temp;
			}
		}

		collided_mat = scene[int_pos[0]][int_pos[1]][int_pos[2]];
	};

	current_hit.mat = collided_mat;
	current_hit.hit_pos = add(curr_ray.pos, mul_vec_scl(len, curr_ray.dir));
	current_hit.norm = norm;

	return current_hit;
}

int main()
{
	int scene[64][64][64];
	float_vec cam = { 32.0f, 32.0f, 63.0f };

	ray init_rays[128][128];
	int colors[128][128][3];

	FILE *img;

	printf("Starting \n");

	/** Scene construction **/
	for (int x = 0; x < 64; x++)
	{
		for (int y = 0; y < 64; y++)
		{
			for (int z = 0; z < 64; z++)
			{
				if ((42 > x && x > 22) && (42 > y && y > 22) && (42 > z && z > 22))
				{
					scene[x][y][z] = 1;
				}
				else
				{
					scene[x][y][z] = 0;
				}

				if (y == 0)
				{
					scene[x][y][z] = 2;
				}
				
				if (z == 0)
				{
					scene[x][y][z] = 3;
				}

				if (x == 63)
				{
					scene[x][y][z] = 4;
				}
			}
		}
	}
	printf("Scene constructed \n");


	for (int x = 0; x < 128; x++)
	{
		for (int y = 0; y < 128; y++)
		{
			init_rays[x][y].dir.x = (((2.0f*((((float)x)+0.5f)/128.0f))-1.0f)*tan(FOV/2.0f));
			init_rays[x][y].dir.y = (1.0f-(2.0f*((((float)y)+0.5f)/128.0f)))*tan(FOV/2.0f);
			init_rays[x][y].dir.z = -1.0f;

			init_rays[x][y].dir = norm_vec(init_rays[x][y].dir);

			init_rays[x][y].pos = cam;
		}
	}

	for (int x = 0; x < 128; x++)
	{
		for (int y = 0; y < 128; y++)
		{
			rayhit hit = dda(init_rays[x][y], scene);
			
			float_vec light_dir = { -0.8f, 0.7f, 0.8f };
			float_vec elipsion = mul_vec_scl(-0.511f, hit.norm);
			float_vec offset_pos = add(elipsion, hit.hit_pos);
			ray shadow_ray = { offset_pos, light_dir };
			rayhit sha = dda(shadow_ray, scene);
			
			light_dir = norm_vec(light_dir);
			if (sha.mat != 0)
			{
				int il = (int)round(255.0f*0.2f*0.4f);

				if (hit.mat == 4)
				{
					colors[x][y][0] = il;
					colors[x][y][1] = 0;
					colors[x][y][2] = 0;
				}
				else if (hit.mat == 3)
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = il;
					colors[x][y][2] = 0;
				}
				else if (hit.mat == 2)
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = 0;
					colors[x][y][2] = il;
				}
				else if (hit.mat == 1)
				{
					colors[x][y][0] = il;
					colors[x][y][1] = il;
					colors[x][y][2] = il;
				}
				else
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = 0;
					colors[x][y][2] = 0;
				}
			}
			else
			{
				int il = (int)MIN(255.0f, round(255.0f*phong(0.2f, 0.8f, 0.9f, 0.4f, 0.1f, 0.9f, 0.7f, hit.norm, norm_vec(add(hit.hit_pos, mul_vec_scl(-1.0f, cam))), mul_vec_scl(-1.0f, light_dir))));

				if (hit.mat == 4)
				{
					colors[x][y][0] = il;
					colors[x][y][1] = 0;
					colors[x][y][2] = 0;
				}
				else if (hit.mat == 3)
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = il;
					colors[x][y][2] = 0;
				}
				else if (hit.mat == 2)
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = 0;
					colors[x][y][2] = il;
				}
				else if (hit.mat == 1)
				{
					colors[x][y][0] = il;
					colors[x][y][1] = il;
					colors[x][y][2] = il;
				}
				else
				{
					colors[x][y][0] = 0;
					colors[x][y][1] = 0;
					colors[x][y][2] = 0;
				}
			}
		}
	}

	img = fopen("./out.ppm", "w");

	fputs("P3\n", img);
	fputs("128 128\n", img);
	fputs("255\n", img);
	for (int x = 0; x < 128; x++)
	{
		for (int y = 0; y < 128; y++)
		{
			for (int c = 0; c < 3; c++)
			{
				fprintf(img, "%d ", colors[x][y][c]);
			}
		}
		fputs("\n", img);
	}
	fclose(img);

	printf("Rendering completed \n");
	return 0;
}