Implement Gaussian anti-aliasing

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),
+	};
+}