Implement Gaussian anti-aliasing

include/rng.h

@@ -10,8 +10,7 @@ typedef struct {
 } rng_t;
 
 rng_t rng_init(unsigned seed);
-
-double rng_double(rng_t *rng);
 vec3_t rng_vec3(rng_t *rng);
+vec3_t rng_gaussian_xy(rng_t *rng, double stddev);
 
 #endif

src/camera.c

@@ -4,12 +4,14 @@
 #include "rng.h"
 
 #include <float.h>
-#include <stdio.h>
 #include <threads.h>
 
 #define MAX_ITER 10
 #define MIN_T 1e-6
+
 #define NSAMPLES 100
+#define SAMPLE_WEIGHT (1.0 / (double)NSAMPLES)
+#define SAMPLE_STDDEV 0.333
 
 #define NTHREADS 40
 
@@ -76,18 +78,23 @@ static int render_thread(void *arg)
 		    = vec3_add(camera->pix_origin, vec3_scale(camera->y_step, y));
 		for (unsigned x = 0; x < w; ++x) {
 			const vec3_t pix = vec3_add(row, vec3_scale(camera->x_step, x));
-			const ray_t ray = {
-				.orig = camera->pos,
-				.dir = vec3_unit(pix),
-			};
 
 			vec3_t colour = black;
 			for (unsigned i = 0; i < NSAMPLES; ++i) {
+				const vec3_t jitter
+				    = rng_gaussian_xy(&slice->rng, SAMPLE_STDDEV);
+				const vec3_t offset = vec3_add(
+				    vec3_scale(camera->x_step, jitter.x),
+				    vec3_scale(camera->y_step, jitter.y));
 
-				const double weight = 1.0 / NSAMPLES;
+				const ray_t ray = {
+					.orig = camera->pos,
+					.dir = vec3_unit(vec3_add(pix, offset)),
+				};
 				const vec3_t sample = trace(
 				    ray, slice->scene, slice->scene_count, &slice->rng);
-				colour = vec3_add(colour, vec3_scale(sample, weight));
+
+				colour = vec3_add(colour, vec3_scale(sample, SAMPLE_WEIGHT));
 			}
 
 			setpix(colour, slice->pixels + (w * y + x));

src/rng.c

@@ -1,8 +1,14 @@
 #include "rng.h"
 
 #include <limits.h>
+#include <math.h>
+#include <stdio.h>
 #include <sys/time.h>
 
+#ifndef M_PI
+#define M_PI 3.14159265258979323846264
+#endif
+
 static uint32_t next(rng_t *rng)
 {
 	uint32_t x = rng->state;
@@ -12,6 +18,16 @@ static uint32_t next(rng_t *rng)
 	return rng->state = x;
 }
 
+static double canonical(rng_t *rng)
+{
+	return (double)next(rng) / (double)UINT32_MAX;
+}
+
+static double disc(rng_t *rng)
+{
+	return 2.0 * canonical(rng) - 1.0;
+}
+
 rng_t rng_init(unsigned seed)
 {
 	struct timeval tv;
@@ -19,13 +35,24 @@ rng_t rng_init(unsigned seed)
 	return (rng_t) { .state = tv.tv_usec + seed };
 }
 
-double rng_double(rng_t *rng)
-{
-	return 2.0 * (double)next(rng) / (double)UINT32_MAX - 1.0;
-}
-
 vec3_t rng_vec3(rng_t *rng)
 {
-	const vec3_t v = { rng_double(rng), rng_double(rng), rng_double(rng) };
+	const vec3_t v = { disc(rng), disc(rng), disc(rng) };
 	return vec3_unit(v);
 }
+
+vec3_t rng_gaussian_xy(rng_t *rng, double stddev)
+{
+	const double r1 = canonical(rng);
+	double r2;
+	do
+		r2 = canonical(rng);
+	while (r2 == 0);
+
+	const double theta = 2.0 * M_PI * r1;
+	const double mag = stddev * sqrt(-2.0 * log(r2));
+	return (vec3_t) {
+		.x = mag * cos(theta),
+		.y = mag * sin(theta),
+	};
+}