Add dielectric material

demo.c

@@ -21,6 +21,7 @@ static const obj_t scene[] = {
 	SPHERE(0.3, -0.7, -1.8, 0.3, LAMBERTIAN(0.9, 0.6, 0.2)),
 	SPHERE(-1.6, -0.8, -2.0, 0.2, REFLECTIVE(1.0, 0.9, 0.4)),
 	SPHERE(-6.0, 5.0, -5.0, 6.0, REFLECTIVE(0.9, 0.9, 0.9)),
+	SPHERE(-0.7, -0.75, -1.5, 0.25, DIELECTRIC(1.5)),
 	SPHERE(0.0, -1001.0, 0.0, 1000.0, LAMBERTIAN(0.3, 0.3, 0.3)),
 };
 

include/material.h

@@ -16,6 +16,12 @@
 		.params = { .reflective = { .tint = { r, g, b } } }, \
 	}
 
+#define DIELECTRIC(e) \
+	{ \
+		.scatter = scatter_dielectric, \
+		.params = { .dielectric = { .eta = e } }, \
+	}
+
 typedef struct {
 	vec3_t point, normal;
 	double t;
@@ -30,9 +36,14 @@ typedef struct {
 	vec3_t tint;
 } reflective_params_t;
 
+typedef struct {
+	double eta;
+} dielectric_params_t;
+
 typedef union {
 	lambertian_params_t lambertian;
 	reflective_params_t reflective;
+	dielectric_params_t dielectric;
 } material_params_t;
 
 typedef bool scatter_fn_t(
@@ -50,5 +61,8 @@ bool scatter_lambertian(
 bool scatter_reflective(
     material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
     vec3_t *atten_out);
+bool scatter_dielectric(
+    material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
+    vec3_t *atten_out);
 
 #endif

include/rng.h

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

include/vec3.h

@@ -11,6 +11,7 @@ vec3_t vec3_sub(vec3_t v, vec3_t u);
 double vec3_dot(vec3_t v, vec3_t u);
 vec3_t vec3_cross(vec3_t v, vec3_t u);
 vec3_t vec3_hadamard(vec3_t v, vec3_t u);
+double vec3_len_squared(vec3_t v);
 double vec3_len(vec3_t v);
 vec3_t vec3_unit(vec3_t v);
 

src/material.c

@@ -1,5 +1,24 @@
 #include "material.h"
 
+#include <math.h>
+#include <stdio.h>
+
+static vec3_t reflect(vec3_t v, vec3_t n)
+{
+	const vec3_t dn = vec3_scale(n, -2 * vec3_dot(n, v));
+	return vec3_unit(vec3_add(v, dn));
+}
+
+static vec3_t refract(vec3_t v, vec3_t n, double eta_ratio)
+{
+	const double cos_theta = fmin(-vec3_dot(v, n), 1.0);
+	const vec3_t perp
+	    = vec3_scale(vec3_add(v, vec3_scale(n, cos_theta)), eta_ratio);
+	const vec3_t para
+	    = vec3_scale(n, -sqrt(fabs(1.0 - vec3_len_squared(perp))));
+	return vec3_unit(vec3_add(perp, para));
+}
+
 bool scatter_lambertian(
     material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
     vec3_t *atten_out)
@@ -17,10 +36,30 @@ bool scatter_reflective(
 {
 	(void)rng;
 
-	const vec3_t dn
-	    = vec3_scale(hit.normal, -2 * vec3_dot(hit.normal, ray->dir));
 	ray->orig = hit.point;
-	ray->dir = vec3_unit(vec3_add(ray->dir, dn));
+	ray->dir = reflect(ray->dir, hit.normal);
 	*atten_out = params.reflective.tint;
 	return true;
 }
+
+bool scatter_dielectric(
+    material_params_t params, hit_t hit, rng_t *rng, ray_t *ray,
+    vec3_t *atten_out)
+{
+	const double eta_ratio
+	    = hit.front ? 1.0 / params.dielectric.eta : params.dielectric.eta;
+	const vec3_t n = hit.front ? hit.normal : vec3_scale(hit.normal, -1);
+	const double cos_theta = fmin(-vec3_dot(ray->dir, n), 1.0);
+	const double sin_theta = sqrt(1.0 - cos_theta * cos_theta);
+	const double r = pow((1 - eta_ratio) / (1 + eta_ratio), 2);
+	const double reflectance = r + (1 - r) * pow(1 - cos_theta, 5);
+
+	ray->orig = hit.point;
+	if (sin_theta * eta_ratio > 1.0 || rng_canon(rng) < reflectance)
+		ray->dir = reflect(ray->dir, n);
+	else
+		ray->dir = refract(ray->dir, n, eta_ratio);
+	*atten_out = (vec3_t) { 1.0, 1.0, 1.0 };
+
+	return true;
+}

src/rng.c

@@ -18,14 +18,9 @@ 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;
+	return 2.0 * rng_canon(rng) - 1.0;
 }
 
 rng_t rng_init(unsigned seed)
@@ -35,6 +30,11 @@ rng_t rng_init(unsigned seed)
 	return (rng_t) { .state = tv.tv_usec + seed };
 }
 
+double rng_canon(rng_t *rng)
+{
+	return (double)next(rng) / (double)UINT32_MAX;
+}
+
 vec3_t rng_vec3(rng_t *rng)
 {
 	const vec3_t v = { disc(rng), disc(rng), disc(rng) };
@@ -43,10 +43,10 @@ vec3_t rng_vec3(rng_t *rng)
 
 vec3_t rng_gaussian_xy(rng_t *rng, double stddev)
 {
-	const double r1 = canonical(rng);
+	const double r1 = rng_canon(rng);
 	double r2;
 	do
-		r2 = canonical(rng);
+		r2 = rng_canon(rng);
 	while (r2 == 0);
 
 	const double theta = 2.0 * M_PI * r1;

src/vec3.c

@@ -49,9 +49,14 @@ vec3_t vec3_hadamard(vec3_t v, vec3_t u)
 	};
 }
 
+double vec3_len_squared(vec3_t v)
+{
+	return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
 double vec3_len(vec3_t v)
 {
-	return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
+	return sqrt(vec3_len_squared(v));
 }
 
 vec3_t vec3_unit(vec3_t v)