improve performance by precomputing parts of the AABB intersection algorithm

A performance analysis showed, that the Double#isFinite checks in the Vec3 constructor add significant overhead without providing much value to the application.

README.md

@@ -0,0 +1,57 @@
+# raytracing
+
+Based on the series <a href="https://raytracing.github.io"><cite>Ray Tracing in One Weekend</cite></a>.
+
+## Scenes
+
+### simple
+![](./docs/simple.png)
+
+```
+java -jar raytracing.jar --samples 5000 --height 1080 SIMPLE
+```
+
+### spheres
+![](./docs/spheres.png)
+
+```
+java -jar raytracing.jar --samples 1000 --height 1080 SPHERES
+```
+### squares
+![](./docs/squares.png)
+
+```
+java -jar raytracing.jar --samples 500 --height 1200 SQUARES
+```
+
+### cornell box
+
+![](./docs/cornell.png)
+
+```
+java -jar raytracing.jar --samples 50000 --height 1200 CORNELL
+```
+
+### cornell box with smoke
+
+![](./docs/cornell_smoke.png)
+
+```
+java -jar raytracing.jar --samples 50000 --height 600 CORNELL_SMOKE
+```
+
+### diagramm
+
+![](./docs/diagramm.png)
+
+```
+java -jar raytracing.jar --samples 1000 --height 1080 DIAGRAMM
+```
+
+### a little bit of everything
+
+![](./docs/final.png)
+
+```
+java -jar raytracing.jar --samples 10000 --height 1200 FINAL
+```

build.gradle.kts

@@ -1,5 +1,16 @@
 plugins {
     id("java")
+    id("application")
+}
+
+java {
+    toolchain {
+        languageVersion = JavaLanguageVersion.of(22)
+    }
+}
+
+application {
+    mainClass = "eu.jonahbauer.raytracing.Main"
 }
 
 dependencies {

src/main/java/eu/jonahbauer/raytracing/Example.java

@@ -0,0 +1,8 @@
+package eu.jonahbauer.raytracing;
+
+import eu.jonahbauer.raytracing.render.camera.Camera;
+import eu.jonahbauer.raytracing.scene.Scene;
+import org.jetbrains.annotations.NotNull;
+
+public record Example(@NotNull Scene scene, @NotNull Camera camera) {
+}

src/main/java/eu/jonahbauer/raytracing/Examples.java

@@ -0,0 +1,415 @@
+package eu.jonahbauer.raytracing;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.color.ColorSpaces;
+import eu.jonahbauer.raytracing.render.spectrum.RGBAlbedoSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.RGBIlluminantSpectrum;
+import eu.jonahbauer.raytracing.render.texture.CheckerTexture;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.camera.SimpleCamera;
+import eu.jonahbauer.raytracing.render.material.*;
+import eu.jonahbauer.raytracing.render.texture.ImageTexture;
+import eu.jonahbauer.raytracing.render.texture.PerlinTexture;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import eu.jonahbauer.raytracing.scene.Scene;
+import eu.jonahbauer.raytracing.scene.SkyBox;
+import eu.jonahbauer.raytracing.scene.hittable2d.Parallelogram;
+import eu.jonahbauer.raytracing.scene.hittable3d.Box;
+import eu.jonahbauer.raytracing.scene.hittable3d.ConstantMedium;
+import eu.jonahbauer.raytracing.scene.hittable3d.Sphere;
+import eu.jonahbauer.raytracing.scene.util.HittableBinaryTree;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.*;
+import java.util.function.IntFunction;
+
+public class Examples {
+    private static final Map<String, IntFunction<Example>> REGISTRY = new HashMap<>();
+
+    private static void register(@NotNull String name, @NotNull IntFunction<Example> example) {
+        REGISTRY.put(name, example);
+    }
+
+    static {
+        register("SIMPLE", Examples::getSimpleScene);
+        register("SPHERES", Examples::getSpheres);
+        register("SQUARES", Examples::getSquares);
+        register("LIGHT", Examples::getLight);
+        register("CORNELL", Examples::getCornellBox);
+        register("CORNELL_SMOKE", Examples::getCornellBoxSmoke);
+        register("CORNELL_SPHERE", Examples::getCornellBoxSphere);
+        register("DIAGRAMM", Examples::getDiagramm);
+        register("EARTH", Examples::getEarth);
+        register("PERLIN", Examples::getPerlin);
+        register("FINAL", Examples::getFinal);
+    }
+
+    public static @NotNull IntFunction<Example> getByName(@NotNull String name) {
+        var out = REGISTRY.get(name);
+        if (out == null) throw new IllegalArgumentException("unknown example " + name + ", expected one of " + REGISTRY.keySet());
+        return out;
+    }
+
+    public static @NotNull Example getSimpleScene(int height) {
+        if (height <= 0) height = 675;
+        var cs = ColorSpaces.sRGB;
+        return new Example(
+                new Scene(getSkyBox(), List.of(
+                        new Sphere(new Vec3(0, -100.5, -1.0), 100.0, new LambertianMaterial(cs.albedo(0.8, 0.8, 0.0))),
+                        new Sphere(new Vec3(0, 0, -1.2), 0.5, new LambertianMaterial(cs.albedo(0.1, 0.2, 0.5))),
+                        new Sphere(new Vec3(-1.0, 0, -1.2), 0.5, new DielectricMaterial(1.5)),
+                        new Sphere(new Vec3(-1.0, 0, -1.2), 0.4, new DielectricMaterial(1 / 1.5)),
+                        new Sphere(new Vec3(1.0, 0, -1.2), 0.5, new MetallicMaterial(cs.albedo(0.8, 0.6, 0.2)))
+                )),
+                SimpleCamera.builder()
+                        .withImage(height * 16 / 9, height)
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getSpheres(int height) {
+        if (height <= 0) height = 675;
+
+        var cs = ColorSpaces.sRGB;
+
+        var rng = new Random(1);
+        var objects = new ArrayList<Hittable>();
+        objects.add(new Sphere(
+                new Vec3(0, -1000, 0), 1000,
+                new LambertianMaterial(new CheckerTexture(0.32,
+                        new RGBAlbedoSpectrum(cs, new ColorRGB(.2, .3, .1)),
+                        new RGBAlbedoSpectrum(cs, new ColorRGB(.9, .9, .9))
+                ))
+        ));
+
+        for (int a = -11; a < 11; a++) {
+            for (int b = -11; b < 11; b++) {
+                var center = new Vec3(a + 0.9 * rng.nextDouble(), 0.2, b + 0.9 * rng.nextDouble());
+                if (Vec3.distance(center, new Vec3(4, 0.2, 0)) <= 0.9) continue;
+
+                Material material;
+                var rnd = rng.nextDouble();
+                if (rnd < 0.8) {
+                    // diffuse
+                    var albedo = ColorRGB.random(rng).times(ColorRGB.random(rng));
+                    material = new LambertianMaterial(cs.albedo(albedo));
+                } else if (rnd < 0.95) {
+                    // metal
+                    var albedo = ColorRGB.random(rng, 0.5, 1.0);
+                    var fuzz = rng.nextDouble() * 0.5;
+                    material = new MetallicMaterial(cs.albedo(albedo), fuzz);
+                } else {
+                    // glass
+                    material = new DielectricMaterial(1.5);
+                }
+
+                objects.add(new Sphere(center, 0.2, material));
+            }
+        }
+
+        objects.add(new Sphere(new Vec3(0, 1, 0), 1.0, new DielectricMaterial(1.5)));
+        objects.add(new Sphere(new Vec3(-4, 1, 0), 1.0, new LambertianMaterial(cs.albedo(0.4, 0.2, 0.1))));
+        objects.add(new Sphere(new Vec3(4, 1, 0), 1.0, new MetallicMaterial(cs.albedo(0.7, 0.6, 0.5))));
+
+        var camera = SimpleCamera.builder()
+                .withImage(height * 16 / 9, height)
+                .withPosition(new Vec3(13, 2, 3))
+                .withTarget(new Vec3(0, 0, 0))
+                .withFieldOfView(Math.toRadians(20))
+                .withFocusDistance(10.0)
+                .withBlurAngle(Math.toRadians(0.6))
+                .build();
+
+        return new Example(new Scene(getSkyBox(), objects), camera);
+    }
+
+    public static @NotNull Example getSquares(int height) {
+        if (height <= 0) height = 600;
+        var cs = ColorSpaces.sRGB;
+        return new Example(
+                new Scene(getSkyBox(), List.of(
+                        new Parallelogram(new Vec3(-3, -2, 5), new Vec3(0, 0, -4), new Vec3(0, 4, 0), new LambertianMaterial(cs.albedo(1.0, 0.2, 0.2))),
+                        new Parallelogram(new Vec3(-2, -2, 0), new Vec3(4, 0, 0), new Vec3(0, 4, 0), new LambertianMaterial(cs.albedo(0.2, 1.0, 0.2))),
+                        new Parallelogram(new Vec3(3, -2, 1), new Vec3(0, 0, 4), new Vec3(0, 4, 0), new LambertianMaterial(cs.albedo(0.2, 0.2, 1.0))),
+                        new Parallelogram(new Vec3(-2, 3, 1), new Vec3(4, 0, 0), new Vec3(0, 0, 4), new LambertianMaterial(cs.albedo(1.0, 0.5, 0.0))),
+                        new Parallelogram(new Vec3(-2, -3, 5), new Vec3(4, 0, 0), new Vec3(0, 0, -4), new LambertianMaterial(cs.albedo(0.2, 0.8, 0.8)))
+                )),
+                SimpleCamera.builder()
+                        .withImage(height, height)
+                        .withFieldOfView(Math.toRadians(80))
+                        .withPosition(new Vec3(0, 0, 9))
+                        .withTarget(new Vec3(0, 0, 0))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getLight(int height) {
+        if (height <= 0) height = 225;
+        var cs = ColorSpaces.sRGB;
+        return new Example(
+                new Scene(List.of(
+                        new Sphere(new Vec3(0, -1000, 0), 1000, new LambertianMaterial(cs.albedo(0.2, 0.2, 0.2))),
+                        new Sphere(new Vec3(0, 2, 0), 2, new LambertianMaterial(cs.albedo(0.2, 0.2, 0.2))),
+                        new Parallelogram(new Vec3(3, 1, -2), new Vec3(2, 0, 0), new Vec3(0, 2, 0), new DiffuseLight(cs.illuminant(4.0))),
+                        new Sphere(new Vec3(0, 7, 0), 2, new DiffuseLight(cs.illuminant(4.0)))
+                )),
+                SimpleCamera.builder()
+                        .withImage(height * 16 / 9, height)
+                        .withFieldOfView(Math.toRadians(20))
+                        .withPosition(new Vec3(26, 3, 6))
+                        .withTarget(new Vec3(0, 2, 0))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getCornellBox(int height) {
+        if (height <= 0) height = 600;
+
+        var cs = ColorSpaces.sRGB;
+        var red = new LambertianMaterial(cs.albedo(.65, .05, .05));
+        var white = new LambertianMaterial(cs.albedo(.73, .73, .73));
+        var green = new LambertianMaterial(cs.albedo(.12, .45, .15));
+        var light = new DiffuseLight(cs.illuminant(15.0));
+
+        return new Example(
+                new Scene(List.of(
+                        new Box(new Vec3(0, 0, 0), new Vec3(555, 555, 555), white, white, red, green, white, null),
+                        new Parallelogram(new Vec3(343, 554, 332), new Vec3(-130, 0, 0), new Vec3(0, 0, -105), light),
+                        new Box(new Vec3(0, 0, 0), new Vec3(165, 330, 165), white)
+                                .rotateY(Math.toRadians(15))
+                                .translate(new Vec3(265, 0, 295)),
+                        new Box(new Vec3(0, 0, 0), new Vec3(165, 165, 165), white)
+                                .rotateY(Math.toRadians(-18))
+                                .translate(new Vec3(130, 0, 65))
+                )),
+                SimpleCamera.builder()
+                        .withImage(height, height)
+                        .withFieldOfView(Math.toRadians(40))
+                        .withPosition(new Vec3(278, 278, -800))
+                        .withTarget(new Vec3(278, 278, 0))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getCornellBoxSmoke(int height) {
+        if (height <= 0) height = 600;
+        var cs = ColorSpaces.sRGB;
+        var red = new LambertianMaterial(cs.albedo(.65, .05, .05));
+        var white = new LambertianMaterial(cs.albedo(.73, .73, .73));
+        var green = new LambertianMaterial(cs.albedo(.12, .45, .15));
+        var light = new DiffuseLight(cs.illuminant(15.0));
+
+        return new Example(
+                new Scene(List.of(
+                        new Box(new Vec3(0, 0, 0), new Vec3(555, 555, 555), white, white, red, green, white, null),
+                        new Parallelogram(new Vec3(113, 554, 127), new Vec3(330, 0, 0), new Vec3(0, 0, 305), light),
+                        new ConstantMedium(
+                                new Box(new Vec3(0, 0, 0), new Vec3(165, 330, 165), white)
+                                        .rotateY(Math.toRadians(15))
+                                        .translate(new Vec3(265, 0, 295)),
+                                0.01, new IsotropicMaterial(cs.albedo(ColorRGB.BLACK))
+                        ),
+                        new ConstantMedium(
+                                new Box(new Vec3(0, 0, 0), new Vec3(165, 165, 165), white)
+                                        .rotateY(Math.toRadians(-18))
+                                        .translate(new Vec3(130, 0, 65)),
+                                0.01, new IsotropicMaterial(cs.albedo(ColorRGB.WHITE))
+                        )
+                )),
+                SimpleCamera.builder()
+                        .withImage(height, height)
+                        .withFieldOfView(Math.toRadians(40))
+                        .withPosition(new Vec3(278, 278, -800))
+                        .withTarget(new Vec3(278, 278, 0))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getCornellBoxSphere(int height) {
+        if (height <= 0) height = 600;
+
+        var cs = ColorSpaces.sRGB;
+        var red = new LambertianMaterial(cs.albedo(.65, .05, .05));
+        var white = new LambertianMaterial(cs.albedo(.73, .73, .73));
+        var green = new LambertianMaterial(cs.albedo(.12, .45, .15));
+        var light = new DiffuseLight(cs.illuminant(7.0));
+        var glass = new DielectricMaterial(1.5);
+
+        var room = new Box(new Vec3(0, 0, 0), new Vec3(555, 555, 555), white, white, red, green, white, null);
+        var lamp = new Parallelogram(new Vec3(343, 554, 332), new Vec3(-130, 0, 0), new Vec3(0, 0, -105), light);
+        var box = new Box(new Vec3(0, 0, 0), new Vec3(165, 330, 165), white)
+                .rotateY(Math.toRadians(15))
+                .translate(new Vec3(265, 0, 295));
+        var sphere = new Sphere(new Vec3(190, 90, 190), 90, glass);
+
+        return new Example(
+                new Scene(List.of(room, box), List.of(lamp, sphere)),
+                SimpleCamera.builder()
+                        .withImage(height, height)
+                        .withFieldOfView(Math.toRadians(40))
+                        .withPosition(new Vec3(278, 278, -800))
+                        .withTarget(new Vec3(278, 278, 0))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getDiagramm(int height) {
+        if (height <= 0) height = 450;
+
+        var cs = ColorSpaces.sRGB;
+        record Partei(String name, ColorRGB color, double stimmen) { }
+        var data = List.of(
+                new Partei("CDU",       new ColorRGB(0x004B76), 18.9),
+                new Partei("SPD",       new ColorRGB(0xC0003C), 25.7),
+                new Partei("AfD",       new ColorRGB(0x80CDEC), 10.3),
+                new Partei("FDP",       new ColorRGB(0xF7BB3D), 11.5),
+                new Partei("DIE LINKE", new ColorRGB(0x5F316E), 4.9),
+                new Partei("GRÜNE",     new ColorRGB(0x00854A), 14.8),
+                new Partei("CSU",       new ColorRGB(0x0077B6), 5.2)
+        );
+        var white = new LambertianMaterial(cs.albedo(.99, .99, .99));
+
+        var count = data.size();
+        var size = 75d;
+        var spacing = 50d;
+
+        var x = (count + 1) * spacing + count * size + 1000;
+        var y = 500 + 1000;
+        var z = 2 * spacing + size + 1000;
+
+        var objects = new ArrayList<Hittable>();
+        objects.add(new Parallelogram(new Vec3(0, 0, 0), new Vec3(x, 0, 0), new Vec3(0, y, 0), white));
+        objects.add(new Parallelogram(new Vec3(0, 0, 0), new Vec3(0, 0, z), new Vec3(x, 0, 0), white));
+        objects.add(new Parallelogram(new Vec3(0, 0, 0), new Vec3(0, y, 0), new Vec3(0, 0, z), white));
+
+        for (int i = 0; i < data.size(); i++) {
+            var partei = data.get(i);
+            objects.add(new Box(
+                    new Vec3((i + 1) * spacing +       i * size, 0, spacing),
+                    new Vec3((i + 1) * spacing + (i + 1) * size, partei.stimmen() * 15, spacing + size),
+                    new DielectricMaterial(1.5, cs.albedo(partei.color()))
+            ));
+        }
+
+        return new Example(
+                new Scene(cs.illuminant().scale(1.25), objects),
+                SimpleCamera.builder()
+                        .withImage(height * 16 / 9, height)
+                        .withPosition(new Vec3(700, 250, 800))
+                        .withTarget(new Vec3(500, 200, 0))
+                        .withFieldOfView(Math.toRadians(40))
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getEarth(int height) {
+        if (height <= 0) height = 450;
+
+        return new Example(
+                new Scene(getSkyBox(), List.of(
+                        new Sphere(Vec3.ZERO, 2, new LambertianMaterial(new ImageTexture("earthmap.jpg", ColorSpaces.sRGB)))
+                )),
+                SimpleCamera.builder()
+                        .withImage(height * 16 / 9, height)
+                        .withFieldOfView(Math.toRadians(20))
+                        .withPosition(new Vec3(12, 0, 0))
+                        .withTarget(Vec3.ZERO)
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getPerlin(int height) {
+        if (height <= 0) height = 450;
+
+        var material = new LambertianMaterial(new PerlinTexture(4));
+
+        return new Example(
+                new Scene(getSkyBox(), List.of(
+                        new Sphere(new Vec3(0, -1000, 0), 1000, material),
+                        new Sphere(new Vec3(0, 2, 0), 2, material)
+                )),
+                SimpleCamera.builder()
+                        .withImage(height * 16 / 9, height)
+                        .withFieldOfView(Math.toRadians(20))
+                        .withPosition(new Vec3(13, 2, 3))
+                        .withTarget(Vec3.ZERO)
+                        .build()
+        );
+    }
+
+    public static @NotNull Example getFinal(int height) {
+        if (height <= 0) height = 400;
+
+        var cs = ColorSpaces.sRGB;
+        var objects = new ArrayList<Hittable>();
+        var random = new Random(1);
+
+        // boxes
+        var boxes = new ArrayList<Hittable>();
+        var ground = new LambertianMaterial(cs.albedo(0.48, 0.83, 0.53));
+        for (int i = 0; i < 20; i++) {
+            for (int j = 0; j < 20; j++) {
+                var w = 100.0;
+                var x0 = -1000.0 + i * w;
+                var z0 = -1000.0 + j * w;
+                var y0 = 0.0;
+                var x1 = x0 + w;
+                var y1 = random.nextInt(1, 101);
+                var z1 = z0 + w;
+                boxes.add(new Box(new Vec3(x0, y0, z0), new Vec3(x1, y1, z1), ground));
+            }
+        }
+        objects.add(new HittableBinaryTree(boxes));
+
+        // light
+        objects.add(new Parallelogram(
+                new Vec3(123, 554, 147), new Vec3(300, 0, 0), new Vec3(0, 0, 265),
+                new DiffuseLight(cs.illuminant(7.0))
+        ));
+
+        // spheres with different materials
+        objects.add(new Sphere(new Vec3(400, 400, 200), 50, new LambertianMaterial(cs.albedo(0.7, 0.3, 0.1))));
+        objects.add(new Sphere(new Vec3(260, 150, 45), 50, new DielectricMaterial(1.5)));
+        objects.add(new Sphere(new Vec3(0, 150, 145), 50, new MetallicMaterial(cs.albedo(0.8, 0.8, 0.9), 1.0)));
+
+        // glass sphere filled with gas
+        var boundary = new Sphere(new Vec3(360, 150, 145), 70, new DielectricMaterial(1.5));
+        objects.add(boundary);
+        objects.add(new ConstantMedium(boundary, 0.2, new IsotropicMaterial(cs.albedo(0.2, 0.4, 0.9))));
+
+        // put the world in a glass sphere
+        objects.add(new ConstantMedium(
+                new Sphere(new Vec3(0, 0, 0), 5000, new DielectricMaterial(1.5)),
+                0.0001, new IsotropicMaterial(cs.albedo(1.0, 1.0, 1.0))
+        ));
+
+        // textures spheres
+        objects.add(new Sphere(new Vec3(400, 200, 400), 100, new LambertianMaterial(new ImageTexture("earthmap.jpg", cs))));
+        objects.add(new Sphere(new Vec3(220, 280, 300), 80, new LambertianMaterial(new PerlinTexture(0.2))));
+
+        // box from spheres
+        var white = new LambertianMaterial(cs.albedo(.73, .73, .73));
+        var spheres = new ArrayList<Hittable>();
+        for (int j = 0; j < 1000; j++) {
+            spheres.add(new Sphere(new Vec3(random.nextDouble(165), random.nextDouble(165), random.nextDouble(165)), 10, white));
+        }
+        objects.add(new HittableBinaryTree(spheres).rotateY(Math.toRadians(15)).translate(new Vec3(-100, 270, 395)));
+
+        return new Example(
+                new Scene(objects),
+                SimpleCamera.builder()
+                        .withImage(height, height)
+                        .withFieldOfView(Math.toRadians(40))
+                        .withPosition(new Vec3(478, 278, -600))
+                        .withTarget(new Vec3(278, 278, 0))
+                        .build()
+        );
+    }
+
+    private static @NotNull SkyBox getSkyBox() {
+        return SkyBox.gradient(
+                new RGBIlluminantSpectrum(ColorSpaces.sRGB, new ColorRGB(0.5, 0.7, 1.0)),
+                ColorSpaces.sRGB.illuminant()
+        );
+    }
+}

src/main/java/eu/jonahbauer/raytracing/Main.java

@@ -1,38 +1,135 @@
 package eu.jonahbauer.raytracing;
 
-import eu.jonahbauer.raytracing.material.DielectricMaterial;
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
-import eu.jonahbauer.raytracing.material.LambertianMaterial;
+import eu.jonahbauer.raytracing.render.canvas.LiveCanvas;
-import eu.jonahbauer.raytracing.material.MetallicMaterial;
+import eu.jonahbauer.raytracing.render.canvas.XYZCanvas;
-import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.image.PNGImageWriter;
-import eu.jonahbauer.raytracing.render.Camera;
+import eu.jonahbauer.raytracing.render.renderer.SimpleRenderer;
-import eu.jonahbauer.raytracing.render.Color;
+import eu.jonahbauer.raytracing.render.color.ColorSpaces;
-import eu.jonahbauer.raytracing.render.ImageFormat;
+import org.jetbrains.annotations.NotNull;
-import eu.jonahbauer.raytracing.scene.Scene;
-import eu.jonahbauer.raytracing.scene.Sphere;
 
 import java.io.IOException;
+import java.nio.file.InvalidPathException;
 import java.nio.file.Path;
+import java.util.function.IntFunction;
 
 public class Main {
-    public static void main(String[] args) throws IOException {
+    public static final boolean DEBUG = false;
-        var scene = new Scene(
-                new Sphere(0, -100.5, - 1, 100, new LambertianMaterial(new Color(0.8, 0.8, 0.0))),
-                new Sphere(0, 0, - 1.2, 0.5, new LambertianMaterial(new Color(0.1, 0.2, 0.5))),
-                new Sphere(-1, 0, - 1, 0.5, new DielectricMaterial(1.5)),
-                new Sphere(-1, 0, - 1, 0.4, new DielectricMaterial(1 / 1.5)),
-                new Sphere(1, 0, - 1, 0.5, new MetallicMaterial(new Color(0.8, 0.6, 0.2), 1.0))
-        );
 
-        var camera = Camera.builder()
+    public static void main(String[] args) throws IOException {
-                .withImage(800, 450)
+        var config = Config.parse(args);
-                .withPosition(new Vec3(-2, 2, 1))
+        var example = config.example;
-                .withTarget(new Vec3(0, 0, -1))
+        var scene = example.scene();
-                .withFieldOfView(Math.toRadians(20))
+        var camera = example.camera();
-                .withFocusDistance(3.4)
+
-                .withBlurAngle(Math.toRadians(10))
+        var renderer = SimpleRenderer.builder()
+                .withSamplesPerPixel(config.samples)
+                .withSpectralSamples(config.spectralSamples)
+                .withMaxDepth(config.depth)
+                .withIterative(config.iterative)
+                .withParallel(config.parallel)
                 .build();
 
-        var image = camera.render(scene);
+        Canvas canvas;
-        ImageFormat.PNG.write(image, Path.of("scene-" + System.currentTimeMillis() + ".png"));
+        if (config.preview) {
+            var image = new LiveCanvas(new XYZCanvas(camera.getWidth(), camera.getHeight()), ColorSpaces.sRGB);
+            image.preview();
+            canvas = image;
+        } else {
+            canvas = new XYZCanvas(camera.getWidth(), camera.getHeight());
+        }
+
+        long time = System.nanoTime();
+        renderer.render(camera, scene, canvas);
+        System.out.printf("rendering finished after %dms", (System.nanoTime() - time) / 1_000_000);
+
+        PNGImageWriter.sRGB.write(canvas, config.path);
+    }
+
+    private record Config(@NotNull Example example, @NotNull Path path, boolean preview, boolean iterative, boolean parallel, int samples, int spectralSamples, int depth) {
+        public static @NotNull Config parse(@NotNull String @NotNull[] args) {
+            IntFunction<Example> example = null;
+            Path path = null;
+            boolean preview = true;
+            boolean iterative = false;
+            boolean parallel = false;
+            int samples = 1000;
+            int spectralSamples = 4;
+            int depth = 50;
+            int height = -1;
+
+            for (int i = 0; i < args.length; i++) {
+                switch (args[i]) {
+                    case "--output" -> {
+                        if (i + 1 == args.length) throw fail("missing value for parameter --output");
+                        try {
+                            path = Path.of(args[++i]);
+                        } catch (InvalidPathException ex) {
+                            throw fail("value " + args[i] + " is not a valid path");
+                        }
+                    }
+                    case "--preview" -> preview = true;
+                    case "--no-preview" -> preview = false;
+                    case "--iterative" -> iterative = true;
+                    case "--no-iterative" -> iterative = false;
+                    case "--parallel" -> parallel = true;
+                    case "--no-parallel" -> parallel = false;
+                    case "--samples" -> {
+                        if (i + 1 == args.length) throw fail("missing value for parameter --samples");
+                        try {
+                            samples = Integer.parseInt(args[++i]);
+                            if (samples <= 0) throw fail("samples must be positive");
+                        } catch (NumberFormatException ex) {
+                            throw fail("value " + args[i] + " is not a valid integer");
+                        }
+                    }
+                    case "--spectral-samples" -> {
+                        if (i + 1 == args.length) throw fail("missing value for parameter --spectral-samples");
+                        try {
+                            spectralSamples = Integer.parseInt(args[++i]);
+                            if (spectralSamples <= 0) throw fail("spectral samples must be positive");
+                        } catch (NumberFormatException ex) {
+                            throw fail("value " + args[i] + " is not a valid integer");
+                        }
+                    }
+                    case "--depth" -> {
+                        if (i + 1 == args.length) throw fail("missing value for parameter --depth");
+                        try {
+                            depth = Integer.parseInt(args[++i]);
+                            if (depth <= 0) throw fail("depth must be positive");
+                        } catch (NumberFormatException ex) {
+                            throw fail("value " + args[i] + " is not a valid integer");
+                        }
+                    }
+                    case "--height" -> {
+                        if (i + 1 == args.length) throw fail("missing value for parameter --height");
+                        try {
+                            height = Integer.parseInt(args[++i]);
+                            if (height <= 0) throw fail("height must be positive");
+                        } catch (NumberFormatException ex) {
+                            throw fail("value " + args[i] + " is not a valid integer");
+                        }
+                    }
+                    case String str when !str.startsWith("-") -> {
+                        try {
+                            example = Examples.getByName(str);
+                        } catch (IllegalArgumentException ex) {
+                            throw fail(ex.getMessage());
+                        }
+                    }
+                    default -> throw fail("unknown option " + args[i]);
+                }
+            }
+
+            if (example == null) example = Examples::getCornellBoxSmoke;
+            if (path == null) path = Path.of("scene-" + System.currentTimeMillis() + ".png");
+            return new Config(example.apply(height), path, preview, iterative, parallel, samples, spectralSamples, depth);
+        }
+
+        private static @NotNull RuntimeException fail(@NotNull String message) {
+            System.err.println(message);
+            System.exit(1);
+            return new RuntimeException();
+        }
     }
 }

src/main/java/eu/jonahbauer/raytracing/material/DielectricMaterial.java

@@ -1,32 +0,0 @@
-package eu.jonahbauer.raytracing.material;
-
-import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
-import eu.jonahbauer.raytracing.render.Color;
-import eu.jonahbauer.raytracing.scene.HitResult;
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Optional;
-
-public record DielectricMaterial(double refractionIndex) implements Material {
-    @Override
-    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit) {
-        var ri = hit.frontFace() ? (1 / refractionIndex) : refractionIndex;
-
-        var cosTheta = Math.min(- ray.direction().unit().times(hit.normal()), 1.0);
-        var reflectance = reflectance(cosTheta);
-        var reflect = reflectance > Math.random();
-
-        var newDirection = (reflect ? Optional.<Vec3>empty() : Vec3.refract(ray.direction(), hit.normal(), ri))
-                .orElseGet(() -> Vec3.reflect(ray.direction(), hit.normal()));
-
-        return Optional.of(new ScatterResult(new Ray(hit.position(), newDirection), Color.WHITE));
-    }
-
-    private double reflectance(double cos) {
-        // use schlick's approximation for reflectance
-        var r0 = (1 - refractionIndex) / (1 + refractionIndex);
-        r0 = r0 * r0;
-        return r0 + (1 - r0) * (1 - cos) * (1 - cos) * (1 - cos) * (1 - cos) * (1 - cos);
-    }
-}

src/main/java/eu/jonahbauer/raytracing/material/LambertianMaterial.java

@@ -1,25 +0,0 @@
-package eu.jonahbauer.raytracing.material;
-
-import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
-import eu.jonahbauer.raytracing.render.Color;
-import eu.jonahbauer.raytracing.scene.HitResult;
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Objects;
-import java.util.Optional;
-
-public record LambertianMaterial(@NotNull Color albedo) implements Material {
-    public LambertianMaterial {
-        Objects.requireNonNull(albedo, "albedo");
-    }
-
-    @Override
-    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit) {
-        var newDirection = hit.normal().plus(Vec3.random(true));
-        if (newDirection.isNearZero()) newDirection = hit.normal();
-
-        var scattered = new Ray(hit.position(), newDirection);
-        return Optional.of(new ScatterResult(scattered, albedo));
-    }
-}

src/main/java/eu/jonahbauer/raytracing/material/Material.java

@@ -1,21 +0,0 @@
-package eu.jonahbauer.raytracing.material;
-
-import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.render.Color;
-import eu.jonahbauer.raytracing.scene.HitResult;
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Objects;
-import java.util.Optional;
-
-public interface Material {
-
-    @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit);
-
-    record ScatterResult(@NotNull Ray ray, @NotNull Color attenuation) {
-        public ScatterResult {
-            Objects.requireNonNull(ray, "ray");
-            Objects.requireNonNull(attenuation, "attenuation");
-        }
-    }
-}

src/main/java/eu/jonahbauer/raytracing/material/MetallicMaterial.java

@@ -1,31 +0,0 @@
-package eu.jonahbauer.raytracing.material;
-
-import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
-import eu.jonahbauer.raytracing.render.Color;
-import eu.jonahbauer.raytracing.scene.HitResult;
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Objects;
-import java.util.Optional;
-
-public record MetallicMaterial(@NotNull Color albedo, double fuzz) implements Material {
-
-    public MetallicMaterial(@NotNull Color albedo) {
-        this(albedo, 0);
-    }
-
-    public MetallicMaterial {
-        Objects.requireNonNull(albedo, "albedo");
-        if (fuzz < 0 || !Double.isFinite(fuzz)) throw new IllegalArgumentException("fuzz must be non-negative");
-    }
-
-    @Override
-    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit) {
-        var newDirection = Vec3.reflect(ray.direction(), hit.normal());
-        if (fuzz > 0) {
-            newDirection = newDirection.unit().plus(Vec3.random(true).times(fuzz));
-        }
-        return Optional.of(new ScatterResult(new Ray(hit.position(), newDirection), albedo));
-    }
-}

src/main/java/eu/jonahbauer/raytracing/math/AABB.java

@@ -0,0 +1,129 @@
+package eu.jonahbauer.raytracing.math;
+
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Comparator;
+import java.util.List;
+import java.util.Optional;
+
+/**
+ * An axis-aligned bounding box.
+ */
+public record AABB(@NotNull Vec3 min, @NotNull Vec3 max) {
+    public static final AABB EMPTY = new AABB(Vec3.ZERO, Vec3.ZERO);
+    public static final Comparator<AABB> X_AXIS = Comparator.comparing(AABB::min, Comparator.comparingDouble(Vec3::x));
+    public static final Comparator<AABB> Y_AXIS = Comparator.comparing(AABB::min, Comparator.comparingDouble(Vec3::y));
+    public static final Comparator<AABB> Z_AXIS = Comparator.comparing(AABB::min, Comparator.comparingDouble(Vec3::z));
+
+    public AABB {
+        var a = min;
+        var b = max;
+        min = Vec3.min(a, b);
+        max = Vec3.max(a, b);
+    }
+
+    public AABB(@NotNull Range x, @NotNull Range y, @NotNull Range z) {
+        this(new Vec3(x.min(), y.min(), z.min()), new Vec3(x.max(), y.max(), z.max()));
+    }
+
+    public static @NotNull Optional<AABB> getBoundingBox(@NotNull List<? extends @NotNull Hittable> objects) {
+        var bbox = (AABB) null;
+        for (var object : objects) {
+            bbox = bbox == null ? object.getBoundingBox() : bbox.expand(object.getBoundingBox());
+        }
+        return Optional.ofNullable(bbox);
+    }
+
+    /**
+     * {@return the range of x values}
+     */
+    public @NotNull Range x() {
+        return new Range(min.x(), max.x());
+    }
+
+    /**
+     * {@return the range of y values}
+     */
+    public @NotNull Range y() {
+        return new Range(min.y(), max.y());
+    }
+
+    /**
+     * {@return the range of z values}
+     */
+    public @NotNull Range z() {
+        return new Range(min.z(), max.z());
+    }
+
+    /**
+     * {@return the center of this bounding box}
+     */
+    public @NotNull Vec3 center() {
+        return Vec3.average(min, max, 2);
+    }
+
+    /**
+     * Expands this bounding box to include the other bounding box.
+     * @param box a bounding box
+     * @return the expanded bounding box
+     */
+    public @NotNull AABB expand(@NotNull AABB box) {
+        return new AABB(Vec3.min(this.min, box.min), Vec3.max(this.max, box.max));
+    }
+
+    /**
+     * Tests whether the {@code ray} intersects this bounding box withing the {@code range}
+     * @param ray a ray
+     * @param range a range of valid {@code t}s
+     * @return {@code true} iff the ray intersects this bounding box, {@code false} otherwise
+     */
+    public boolean hit(@NotNull Ray ray, @NotNull Range range) {
+        var invDirection = ray.getInvDirection();
+        var negInvOrigin = ray.getNegInvOrigin();
+
+        var tminX = Math.fma(min.x(), invDirection.x(), negInvOrigin.x());
+        var tminY = Math.fma(min.y(), invDirection.y(), negInvOrigin.y());
+        var tminZ = Math.fma(min.z(), invDirection.z(), negInvOrigin.z());
+
+        var tmaxX = Math.fma(max.x(), invDirection.x(), negInvOrigin.x());
+        var tmaxY = Math.fma(max.y(), invDirection.y(), negInvOrigin.y());
+        var tmaxZ = Math.fma(max.z(), invDirection.z(), negInvOrigin.z());
+
+        var tlmax = max(
+                Math.min(tminX, tmaxX),
+                Math.min(tminY, tmaxY),
+                Math.min(tminZ, tmaxZ)
+        );
+        var tumin = min(
+                Math.max(tminX, tmaxX),
+                Math.max(tminY, tmaxY),
+                Math.max(tminZ, tmaxZ)
+        );
+
+        return tlmax < tumin && tumin >= range.min() && tlmax <= range.max();
+    }
+
+    private static double max(double a, double b, double c) {
+        return Math.max(a, Math.max(b, c));
+    }
+
+    private static double min(double a, double b, double c) {
+        return Math.min(a, Math.min(b, c));
+    }
+
+    /**
+     * Computes the {@code t} values of the intersections of a ray with the axis-aligned planes through a point.
+     * @param corner the point
+     * @param origin the origin point of the ray
+     * @param invDirection the {@linkplain Vec3#inv() inverted} direction of the ray
+     * @return a three-element array of the {@code t} values of the intersection with the yz-, xz- and xy-plane through {@code corner}
+     */
+    public static double @NotNull[] intersect(@NotNull Vec3 corner, @NotNull Vec3 origin, @NotNull Vec3 invDirection) {
+        return new double[] {
+                (corner.x() - origin.x()) * invDirection.x(),
+                (corner.y() - origin.y()) * invDirection.y(),
+                (corner.z() - origin.z()) * invDirection.z(),
+        };
+    }
+}

src/main/java/eu/jonahbauer/raytracing/math/IVec.java

@@ -0,0 +1,33 @@
+package eu.jonahbauer.raytracing.math;
+
+import org.jetbrains.annotations.Contract;
+import org.jetbrains.annotations.NotNull;
+
+/**
+ * A vector-like object that implements the standard mathematical operations
+ * @param <T> the type
+ */
+public interface IVec<T extends IVec<T>> {
+    @Contract(pure = true)
+    double get(int i);
+
+    @Contract(pure = true)
+    @NotNull T plus(@NotNull T other);
+
+    @Contract(pure = true)
+    @NotNull T minus(@NotNull T other);
+
+    @Contract(pure = true)
+    @NotNull T times(@NotNull T other);
+
+    @Contract(pure = true)
+    @NotNull T times(double d);
+
+    @Contract(pure = true)
+    default @NotNull T div(double d) {
+        return times(1 / d);
+    }
+
+    @Contract(pure = true)
+    double @NotNull[] toArray();
+}

src/main/java/eu/jonahbauer/raytracing/math/IVec3.java

@@ -0,0 +1,47 @@
+package eu.jonahbauer.raytracing.math;
+
+import org.jetbrains.annotations.Contract;
+import org.jetbrains.annotations.NotNull;
+
+/**
+ * A vector-like object with three components.
+ * @param <T> the type
+ */
+public interface IVec3<T extends Record & IVec3<T>> extends IVec<T> {
+    @Contract(pure = true)
+    default double component1() {
+        return toVec3().x();
+    }
+
+    @Contract(pure = true)
+    default double component2() {
+        return toVec3().y();
+    }
+
+    @Contract(pure = true)
+    default double component3() {
+        return toVec3().z();
+    }
+
+    @Override
+    @Contract(pure = true)
+    default double get(int i) {
+        return switch (i) {
+            case 0 -> component1();
+            case 1 -> component2();
+            case 2 -> component3();
+            default -> throw new IndexOutOfBoundsException(i);
+        };
+    }
+
+    @Contract(pure = true)
+    default @NotNull Vec3 toVec3() {
+        return new Vec3(component1(), component2(), component3());
+    }
+
+    @Override
+    @Contract(pure = true)
+    default double @NotNull [] toArray() {
+        return new double[] {component1(), component2(), component3()};
+    }
+}

src/main/java/eu/jonahbauer/raytracing/math/Matrix3.java

@@ -0,0 +1,255 @@
+package eu.jonahbauer.raytracing.math;
+
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+
+public record Matrix3(
+        double a11, double a12, double a13,
+        double a21, double a22, double a23,
+        double a31, double a32, double a33
+) {
+    public static @NotNull Matrix3 fromRows(@NotNull Vec3 @NotNull[] rows) {
+        if (rows.length != 3) throw new IllegalArgumentException();
+        return fromRows(rows[0], rows[1], rows[2]);
+    }
+
+    public static @NotNull Matrix3 fromRows(@NotNull Vec3 row0, @NotNull Vec3 row1, @NotNull Vec3 row2) {
+        return new Matrix3(
+                row0.x(), row0.y(), row0.z(),
+                row1.x(), row1.y(), row1.z(),
+                row2.x(), row2.y(), row2.z()
+        );
+    }
+
+    public static @NotNull Matrix3 fromColumns(@NotNull Vec3 @NotNull[] cols) {
+        if (cols.length != 3) throw new IllegalArgumentException();
+        return fromColumns(cols[0], cols[1], cols[2]);
+    }
+
+    public static @NotNull Matrix3 fromColumns(@NotNull Vec3 col0, @NotNull Vec3 col1, @NotNull Vec3 col2) {
+        return new Matrix3(
+                col0.x(), col1.x(), col2.x(),
+                col0.y(), col1.y(), col2.y(),
+                col0.z(), col1.z(), col2.z()
+        );
+    }
+
+    public static @NotNull Matrix3 fromArray(double @NotNull[] @NotNull[] array) {
+        return new Matrix3(
+                array[0][0], array[0][1], array[0][2],
+                array[1][0], array[1][1], array[1][2],
+                array[2][0], array[2][1], array[2][2]
+        );
+    }
+
+    public Matrix3() {
+        this(1, 1, 1);
+    }
+
+    public Matrix3(double a11, double a22, double a33) {
+        this(a11, 0, 0, 0, a22, 0, 0, 0, a33);
+    }
+
+    public @NotNull Matrix3 times(@NotNull Matrix3 other) {
+        return new Matrix3(
+                a11 * other.a11 + a12 * other.a21 + a13 * other.a31,
+                a11 * other.a12 + a12 * other.a22 + a13 * other.a32,
+                a11 * other.a13 + a12 * other.a23 + a13 * other.a33,
+                a21 * other.a11 + a22 * other.a21 + a23 * other.a31,
+                a21 * other.a12 + a22 * other.a22 + a23 * other.a32,
+                a21 * other.a13 + a22 * other.a23 + a23 * other.a33,
+                a31 * other.a11 + a32 * other.a21 + a33 * other.a31,
+                a31 * other.a12 + a32 * other.a22 + a33 * other.a32,
+                a31 * other.a13 + a32 * other.a23 + a33 * other.a33
+        );
+    }
+
+    public @NotNull Matrix3 times(double other) {
+        return new Matrix3(
+                a11 * other, a12 * other, a13 * other,
+                a21 * other, a22 * other, a23 * other,
+                a31 * other, a32 * other, a33 * other
+        );
+    }
+
+    public @NotNull Vec3 times(@NotNull Vec3 other) {
+        return new Vec3(
+                a11 * other.x() + a12 * other.y() + a13 * other.z(),
+                a21 * other.x() + a22 * other.y() + a23 * other.z(),
+                a31 * other.x() + a32 * other.y() + a33 * other.z()
+        );
+    }
+
+    public @NotNull Matrix3 plus(@NotNull Matrix3 other) {
+        return new Matrix3(
+                a11 + other.a11, a12 + other.a12, a13 + other.a13,
+                a21 + other.a21, a22 + other.a22, a23 + other.a23,
+                a31 + other.a31, a32 + other.a32, a33 + other.a33
+        );
+    }
+
+    public double det() {
+        return a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32
+                - a13 * a22 * a31 - a23 * a32 * a11 - a33 * a12 * a21;
+    }
+
+    public @NotNull Matrix3 invert() {
+        var det = det();
+        if (det == 0) throw new IllegalStateException();
+        var t = 1 / det;
+        return new Matrix3(
+                t * (Math.fma( a22, a33, -a23 * a32)),
+                t * (Math.fma(-a12, a33,  a13 * a32)),
+                t * (Math.fma( a12, a23, -a13 * a22)),
+                t * (Math.fma(-a21, a33,  a23 * a31)),
+                t * (Math.fma( a11, a33, -a13 * a31)),
+                t * (Math.fma(-a11, a23,  a13 * a21)),
+                t * (Math.fma( a21, a32, -a22 * a31)),
+                t * (Math.fma(-a11, a32,  a12 * a31)),
+                t * (Math.fma( a11, a22, -a12 * a21))
+        );
+    }
+
+    public @NotNull Vec3 column(int i) {
+        return switch (i) {
+            case 0 -> new Vec3(a11, a21, a31);
+            case 1 -> new Vec3(a12, a22, a32);
+            case 2 -> new Vec3(a13, a23, a33);
+            default -> throw new IndexOutOfBoundsException(i);
+        };
+    }
+
+    public @NotNull Vec3 @NotNull[] columns() {
+        return new Vec3[] {
+                new Vec3(a11, a21, a31),
+                new Vec3(a12, a22, a32),
+                new Vec3(a13, a23, a33)
+        };
+    }
+
+    public @NotNull Vec3 row(int i) {
+        return switch (i) {
+            case 0 -> new Vec3(a11, a12, a13);
+            case 1 -> new Vec3(a21, a22, a23);
+            case 2 -> new Vec3(a31, a32, a33);
+            default -> throw new IndexOutOfBoundsException(i);
+        };
+    }
+
+    public @NotNull Vec3 @NotNull[] rows() {
+        return new Vec3[] {
+                new Vec3(a11, a12, a13),
+                new Vec3(a21, a22, a23),
+                new Vec3(a31, a32, a33)
+        };
+    }
+
+    public double @NotNull[] @NotNull[] toArray() {
+        return new double[][] {
+                {a11, a12, a13},
+                {a21, a22, a23},
+                {a31, a32, a33}
+        };
+    }
+
+    public double get(int i, int j) {
+        Objects.checkIndex(i, 3);
+        Objects.checkIndex(j, 3);
+        var idx = 3 * i + j;
+        return switch (idx) {
+            case 0 -> a11;
+            case 1 -> a12;
+            case 2 -> a13;
+            case 3 -> a21;
+            case 4 -> a22;
+            case 5 -> a23;
+            case 6 -> a31;
+            case 7 -> a32;
+            case 8 -> a33;
+            default -> throw new AssertionError();
+        };
+    }
+
+    /**
+     * Performs lower-upper decomposition with partial pivoting (LUP decomposition) on {@code this} matrix.
+     * @param tolerance a small tolerance number to detect failure when the matrix is near degenerate
+     * @see <a href="https://en.wikipedia.org/w/index.php?title=LU_decomposition&oldid=1213102558#C_code_example">LU decomposition — Wikipedia, The Free Encyclopedia</a>
+     */
+    public @NotNull LUPDecomposition decompose(double tolerance) {
+        // unit permutation matrix
+        var perm = new int[] {0, 1, 2, 3};
+        var A = toArray();
+        var N = 3;
+
+        for (int i = 0; i < N; i++) {
+            double maxA = 0.0;
+            int imax = i;
+
+            for (int k = i; k < N; k++) {
+                double absA = Math.abs(A[k][i]);
+                if (absA > maxA) {
+                    maxA = absA;
+                    imax = k;
+                }
+            }
+
+            if (maxA < tolerance) throw new IllegalArgumentException("matrix is degenerate");
+
+            if (imax != i) {
+                // pivoting P
+                int j = perm[i];
+                perm[i] = perm[imax];
+                perm[imax] = j;
+
+                // pivoting rows of A
+                var ptr = A[i];
+                A[i] = A[imax];
+                A[imax] = ptr;
+
+                // counting pivots starting from N (for determinant)
+                perm[3]++;
+            }
+
+            for (int j = i + 1; j < N; j++) {
+                A[j][i] /= A[i][i];
+
+                for (int k = i + 1; k < N; k++) {
+                    A[j][k] -= A[j][i] * A[i][k];
+                }
+            }
+        }
+
+        return new LUPDecomposition(fromArray(A), perm);
+    }
+
+    public record LUPDecomposition(@NotNull Matrix3 matrix, int @NotNull[] permutation) {
+
+        /**
+         * Solves the equation {@code Ax = b} where {@code A} is the matrix that {@code this} decomposition was derived
+         * from.
+         * @param b the right hand side vector
+         * @return the solution vector
+         */
+        public @NotNull Vec3 solve(@NotNull Vec3 b) {
+            var N = 3;
+            var x = new double[N];
+            for (int i = 0; i < N; i++) {
+                x[i] = b.get(permutation[i]);
+
+                for (int k = 0; k < i; k++) {
+                    x[i] -= matrix.get(i, k) * x[k];
+                }
+            }
+
+            for (int i = N - 1; i >= 0; i--) {
+                for (int k = i + 1; k < N; k++) {
+                    x[i] -= matrix.get(i, k) * x[k];
+                }
+
+                x[i] /= matrix.get(i, i);
+            }
+            return new Vec3(x[0], x[1], x[2]);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/math/Range.java

@@ -16,4 +16,8 @@ public record Range(double min, double max) {
     public boolean surrounds(double value) {
         return min < value && value < max;
     }
+
+    public double size() {
+        return max - min;
+    }
 }

src/main/java/eu/jonahbauer/raytracing/math/Ray.java

@@ -1,21 +1,97 @@
 package eu.jonahbauer.raytracing.math;
 
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.scene.HitResult;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.Objects;
 
-public record Ray(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+public final class Ray {
-    public Ray {
+    private final @NotNull Vec3 origin;
-        Objects.requireNonNull(origin, "origin");
+    private final @NotNull Vec3 direction;
-        Objects.requireNonNull(direction, "direction");
+    private final @NotNull SampledWavelengths lambda;
+
+    private final @NotNull Vec3 inv;
+    private final @NotNull Vec3 negInvOrigin;
+
+    public Ray(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        this(origin, direction, SampledWavelengths.EMPTY);
+    }
+
+    public Ray(@NotNull Vec3 origin, @NotNull Vec3 direction, @NotNull SampledWavelengths lambda) {
+        this.origin = Objects.requireNonNull(origin, "origin");
+        this.direction = Objects.requireNonNull(direction, "direction");
+        this.lambda = Objects.requireNonNull(lambda, "lambda");
+
+        this.inv = direction.inv();
+        this.negInvOrigin = inv.neg().times(origin);
+    }
+
+    private Ray(@NotNull Vec3 origin, @NotNull Vec3 direction, @NotNull SampledWavelengths lambda, @NotNull Vec3 inv, @NotNull Vec3 negInvOrigin) {
+        this.origin = origin;
+        this.direction = direction;
+        this.lambda = lambda;
+        this.inv = inv;
+        this.negInvOrigin = negInvOrigin;
     }
 
     public @NotNull Vec3 at(double t) {
-        if (t < 0) throw new IllegalArgumentException("t must not be negative");
+        return Vec3.fma(t, direction, origin);
-        return new Vec3(
-                origin().x() + t * direction.x(),
-                origin().y() + t * direction.y(),
-                origin().z() + t * direction.z()
-        );
     }
+
+    public @NotNull Ray with(@NotNull HitResult hit, @NotNull Vec3 direction) {
+        return new Ray(hit.position(), direction, lambda);
+    }
+
+    public @NotNull Ray with(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        return new Ray(origin, direction, lambda);
+    }
+
+    public @NotNull Ray with(@NotNull SampledWavelengths lambda) {
+        return new Ray(origin, direction, lambda, inv, negInvOrigin);
+    }
+
+    public @NotNull Vec3 origin() {
+        return origin;
+    }
+
+    public @NotNull Vec3 direction() {
+        return direction;
+    }
+
+    public @NotNull SampledWavelengths lambda() {
+        return lambda;
+    }
+
+    public @NotNull Vec3 getInvDirection() {
+        return inv;
+    }
+
+    public @NotNull Vec3 getNegInvOrigin() {
+        return negInvOrigin;
+    }
+
+    @Override
+    public boolean equals(Object obj) {
+        if (obj == this) return true;
+        if (obj == null || obj.getClass() != this.getClass()) return false;
+        var that = (Ray) obj;
+        return Objects.equals(this.origin, that.origin) &&
+                Objects.equals(this.direction, that.direction) &&
+                Objects.equals(this.lambda, that.lambda);
+    }
+
+    @Override
+    public int hashCode() {
+        return Objects.hash(origin, direction, lambda);
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return "Ray[" +
+                "origin=" + origin + ", " +
+                "direction=" + direction + ", " +
+                "lambda=" + lambda + ']';
+    }
+
 }

src/main/java/eu/jonahbauer/raytracing/math/Vec3.java

@@ -3,39 +3,77 @@ package eu.jonahbauer.raytracing.math;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.Optional;
+import java.util.random.RandomGenerator;
 
-public record Vec3(double x, double y, double z) {
+import static eu.jonahbauer.raytracing.Main.DEBUG;
+
+public record Vec3(double x, double y, double z) implements IVec3<Vec3> {
     public static final Vec3 ZERO = new Vec3(0, 0, 0);
     public static final Vec3 UNIT_X = new Vec3(1, 0, 0);
     public static final Vec3 UNIT_Y = new Vec3(0, 1, 0);
     public static final Vec3 UNIT_Z = new Vec3(0, 0, 1);
 
     public Vec3 {
-        if (!Double.isFinite(x) || !Double.isFinite(y) || !Double.isFinite(z)) {
+        if (DEBUG) {
-            throw new IllegalArgumentException("x, y and z must be finite");
+            if (!Double.isFinite(x) || !Double.isFinite(y) || !Double.isFinite(z)) {
+                throw new IllegalArgumentException("x, y and z must be finite");
+            }
         }
     }
 
-    public static @NotNull Vec3 random() {
+    /**
-        return random(false);
+     * {@return a uniformly random unit vector}
+     */
+    public static @NotNull Vec3 random(@NotNull RandomGenerator random) {
+        double x, y, z;
+        double squared;
+        do {
+            x = Math.fma(2, random.nextDouble(), -1);
+            y = Math.fma(2, random.nextDouble(), -1);
+            z = Math.fma(2, random.nextDouble(), -1);
+            squared = x * x + y * y + z * z;
+        } while (squared > 1);
+        var factor = 1 / Math.sqrt(squared);
+        return new Vec3(x * factor, y * factor, z * factor);
     }
 
-    public static @NotNull Vec3 random(boolean unit) {
+    /**
-        var random = new Vec3(
+     * {@return a uniformly random unit vector on the opposite hemisphere of the given <code>direction</code>}
-                2 * Math.random() - 1,
+     */
-                2 * Math.random() - 1,
+    public static @NotNull Vec3 randomOppositeHemisphere(@NotNull RandomGenerator random, @NotNull Vec3 direction) {
-                2 * Math.random() - 1
+        double x, y, z;
-        );
+        double squared;
-        return unit ? random.unit() : random;
+        do {
+            x = Math.fma(2, random.nextDouble(), -1);
+            y = Math.fma(2, random.nextDouble(), -1);
+            z = Math.fma(2, random.nextDouble(), -1);
+            squared = x * x + y * y + z * z;
+        } while (squared > 1 || direction.x() * x + direction.y() * y + direction.z() * z >= 0);
+        var factor = 1 / Math.sqrt(squared);
+        return new Vec3(x * factor, y * factor, z * factor);
     }
 
+    /**
+     * Reflects a vector on the given {@code normal} vector.
+     * @param vec a vector
+     * @param normal the surface normal (must be a unit vector)
+     * @return the reflected vector
+     */
     public static @NotNull Vec3 reflect(@NotNull Vec3 vec, @NotNull Vec3 normal) {
-        return vec.minus(normal.times(2 * normal.times(vec)));
+        var factor = - 2 * normal.dot(vec);
+        return Vec3.fma(factor, normal, vec);
     }
 
+    /**
+     * Refracts a vector on the given {@code normal} vector.
+     * @param vec a vector
+     * @param normal the surface normal (must be a unit vector)
+     * @param ri the refractive index
+     * @return the refracted vector
+     */
     public static @NotNull Optional<Vec3> refract(@NotNull Vec3 vec, @NotNull Vec3 normal, double ri) {
         vec = vec.unit();
-        var cosTheta = Math.min(- vec.times(normal), 1.0);
+        var cosTheta = Math.min(- vec.dot(normal), 1.0);
         var sinTheta = Math.sqrt(1 - cosTheta * cosTheta);
         if (ri * sinTheta > 1) return Optional.empty();
 
@@ -44,59 +82,236 @@ public record Vec3(double x, double y, double z) {
         return Optional.of(rOutPerp.plus(rOutParallel));
     }
 
+    /**
+     * Rotates a vector around an {@code axis}.
+     * @param vec a vector
+     * @param axis the rotation axis
+     * @param angle the angle in radians
+     * @return the rotated vector
+     */
     public static @NotNull Vec3 rotate(@NotNull Vec3 vec, @NotNull Vec3 axis, double angle) {
         Vec3 vxp = axis.cross(vec);
         Vec3 vxvxp = axis.cross(vxp);
         return vec.plus(vxp.times(Math.sin(angle))).plus(vxvxp.times(1 - Math.cos(angle)));
     }
 
-    public @NotNull Vec3 plus(@NotNull Vec3 b) {
+    /**
-        return new Vec3(this.x + b.x, this.y + b.y, this.z + b.z);
+     * {@return the euclidean distance between two vectors}
+     * @param a a vector
+     * @param b another vector
+     */
+    public static double distance(@NotNull Vec3 a, @NotNull Vec3 b) {
+        return a.minus(b).length();
     }
-    
+
-    public @NotNull Vec3 minus(@NotNull Vec3 b) {
+    /**
-        return new Vec3(this.x - b.x, this.y - b.y, this.z - b.z);
+     * Computes a running average of vectors.
+     * @param current the current running average
+     * @param next the next vector
+     * @param index the one-based index of the next vector
+     * @return the new running average
+     */
+    public static @NotNull Vec3 average(@NotNull Vec3 current, @NotNull Vec3 next, int index) {
+        var factor = 1d / index;
+        return new Vec3(
+                Math.fma(factor, next.x() - current.x(), current.x()),
+                Math.fma(factor, next.y() - current.y(), current.y()),
+                Math.fma(factor, next.z() - current.z(), current.z())
+        );
     }
-    
+
-    public double times(@NotNull Vec3 b) {
+    /**
-        return this.x * b.x + this.y * b.y + this.z * b.z;
+     * {@return a component-wise maximum vector}
+     * @param a a vector
+     * @param b another vector
+     */
+    public static @NotNull Vec3 max(@NotNull Vec3 a, @NotNull Vec3 b) {
+        return new Vec3(
+                Math.max(a.x(), b.x()),
+                Math.max(a.y(), b.y()),
+                Math.max(a.z(), b.z())
+        );
     }
-    
+
-    public @NotNull Vec3 times(double b) {
+    /**
-        return new Vec3(this.x * b, this.y * b, this.z * b);
+     * {@return a component-wise minimum vector}
+     * @param a a vector
+     * @param b another vector
+     */
+    public static @NotNull Vec3 min(@NotNull Vec3 a, @NotNull Vec3 b) {
+        return new Vec3(
+                Math.min(a.x(), b.x()),
+                Math.min(a.y(), b.y()),
+                Math.min(a.z(), b.z())
+        );
+    }
+
+    /**
+     * {@return <code>a * b + c</code>}
+     * @param a scalar
+     * @param b a vector
+     * @param c another vector
+     */
+    public static @NotNull Vec3 fma(double a, @NotNull Vec3 b, @NotNull Vec3 c) {
+        return new Vec3(
+                Math.fma(a, b.x(), c.x()),
+                Math.fma(a, b.y(), c.y()),
+                Math.fma(a, b.z(), c.z())
+        );
+    }
+
+    public static @NotNull Vec3 fma(@NotNull Vec3 a, @NotNull Vec3 b, @NotNull Vec3 c) {
+        return new Vec3(
+                Math.fma(a.x(), b.x(), c.x()),
+                Math.fma(a.y(), b.y(), c.y()),
+                Math.fma(a.z(), b.z(), c.z())
+        );
+    }
+
+    public static double tripleProduct(@NotNull Vec3 a, @NotNull Vec3 b, @NotNull Vec3 c) {
+        return a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y - c.x * b.y * a.z - c.y * b.z * a.x - c.z * b.x * a.y;
+    }
+
+    /*
+     * Math
+     */
+
+    public @NotNull Vec3 plus(double x, double y, double z) {
+        return new Vec3(this.x + x, this.y + y, this.z + z);
+    }
+
+    public @NotNull Vec3 minus(double x, double y, double z) {
+        return new Vec3(this.x - x, this.y - y, this.z - z);
+    }
+
+    /**
+     * Adds a vector to this vector
+     * @param other a vector
+     * @return the sum of this and the other vector
+     */
+    @Override
+    public @NotNull Vec3 plus(@NotNull Vec3 other) {
+        return new Vec3(this.x + other.x, this.y + other.y, this.z + other.z);
+    }
+
+    /**
+     * Subtracts a vector from this vector
+     * @param other a vector
+     * @return the difference of this and the other vector
+     */
+    @Override
+    public @NotNull Vec3 minus(@NotNull Vec3 other) {
+        return new Vec3(this.x - other.x, this.y - other.y, this.z - other.z);
+    }
+
+    /**
+     * Multiplies this vector with a scalar
+     * @param t a scalar
+     * @return the product of this vector and the scalar
+     */
+    @Override
+    public @NotNull Vec3 times(double t) {
+        return new Vec3(this.x * t, this.y * t, this.z * t);
+    }
+
+    /**
+     * Multiplies this vector with another vector component-wise.
+     * @param other a vector
+     * @return the component-wise product of this vector and the other vector
+     */
+    @Override
+    public @NotNull Vec3 times(@NotNull Vec3 other) {
+        return new Vec3(this.x * other.x, this.y * other.y, this.z * other.z);
+    }
+
+    /**
+     * Divides this vector by a scalar
+     * @param t a scalar
+     * @return this vector divided by the scalar
+     */
+    @Override
+    public @NotNull Vec3 div(double t) {
+        return times(1 / t);
+    }
+
+    /**
+     * Computes the scalar product of this and another vector
+     * @param other a vector
+     * @return the scalar product
+     */
+    public double dot(@NotNull Vec3 other) {
+        return this.x * other.x + this.y * other.y + this.z * other.z;
     }
 
     public @NotNull Vec3 neg() {
         return new Vec3(-x, -y, -z);
     }
 
-    public @NotNull Vec3 cross(@NotNull Vec3 b) {
+    public @NotNull Vec3 inv() {
+        return new Vec3(1 / x, 1 / y, 1 / z);
+    }
+
+    public @NotNull Vec3 cross(@NotNull Vec3 other) {
         return new Vec3(
-                this.y() * b.z() - b.y() * this.z(),
+                Math.fma(this.y, other.z, - other.y * this.z),
-                this.z() * b.x() - b.z() * this.x(),
+                Math.fma(this.z, other.x, - other.z * this.x),
-                this.x() * b.y() - b.x() * this.y()
+                Math.fma(this.x, other.y, - other.x * this.y)
         );
     }
 
-    public @NotNull Vec3 div(double b) {
+    /**
-        return new Vec3(this.x / b, this.y / b, this.z / b);
+     * {@return the squared length of this vector}
-    }
+     */
-
     public double squared() {
         return this.x * this.x + this.y * this.y + this.z * this.z;
     }
-    
+
+    /**
+     * {@return the length of this vector}
+     */
     public double length() {
         return Math.sqrt(squared());
     }
 
+    /**
+     * {@return whether this vector is near zero}
+     */
     public boolean isNearZero() {
         var s = 1e-8;
         return Math.abs(x) < s && Math.abs(y) < s && Math.abs(z) < s;
     }
-    
+
+    /**
+     * {@return a unit vector with the same direction as this vector}
+     */
     public @NotNull Vec3 unit() {
-        return div(length());
+        var squared = squared();
+        if (squared == 1) return this;
+        return div(Math.sqrt(squared));
+    }
+
+    /*
+     * Accessors
+     */
+
+    @Override
+    public double component1() {
+        return x;
+    }
+
+    @Override
+    public double component2() {
+        return y;
+    }
+
+    @Override
+    public double component3() {
+        return z;
+    }
+
+    @Override
+    public @NotNull Vec3 toVec3() {
+        return this;
     }
 
     public @NotNull Vec3 withX(double x) {
@@ -110,4 +325,9 @@ public record Vec3(double x, double y, double z) {
     public @NotNull Vec3 withZ(double z) {
         return new Vec3(x, y, z);
     }
+
+    @Override
+    public @NotNull String toString() {
+        return "(" + x + "," + y + "," + z + ")";
+    }
 }

src/main/java/eu/jonahbauer/raytracing/render/Color.java

@@ -1,48 +0,0 @@
-package eu.jonahbauer.raytracing.render;
-
-import org.jetbrains.annotations.NotNull;
-
-public record Color(double r, double g, double b) {
-    public static final @NotNull Color BLACK = new Color(0.0, 0.0, 0.0);
-    public static final @NotNull Color WHITE = new Color(1.0, 1.0, 1.0);
-    public static final @NotNull Color SKY = new Color(0.5, 0.7, 1.0);
-    public static final @NotNull Color RED = new Color(1.0, 0.0, 0.0);
-    public static final @NotNull Color GREEN = new Color(0.0, 1.0, 0.0);
-    public static final @NotNull Color BLUE = new Color(0.0, 0.0, 1.0);
-
-    public static @NotNull Color lerp(@NotNull Color a, @NotNull Color b, double t) {
-        if (t < 0) return a;
-        if (t > 1) return b;
-        return new Color(
-                (1 - t) * a.r + t * b.r,
-                (1 - t) * a.g + t * b.g,
-                (1 - t) * a.b + t * b.b
-        );
-    }
-
-    public static @NotNull Color multiply(@NotNull Color a, @NotNull Color b) {
-        return new Color(a.r() * b.r(), a.g() * b.g(), a.b() * b.b());
-    }
-
-    public Color {
-        if (r < 0 || r > 1 || g < 0 || g > 1 || b < 0 || b > 1) {
-            throw new IllegalArgumentException("r, g and b must be in the range 0 to 1");
-        }
-    }
-    
-    public Color(int red, int green, int blue) {
-        this(red / 255f, green / 255f, blue / 255f);
-    }
-
-    public int red() {
-        return (int) (255.99 * r);
-    }
-
-    public int green() {
-        return (int) (255.99 * g);
-    }
-
-    public int blue() {
-        return (int) (255.99 * b);
-    }
-}

src/main/java/eu/jonahbauer/raytracing/render/Image.java

@@ -1,56 +0,0 @@
-package eu.jonahbauer.raytracing.render;
-
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Arrays;
-import java.util.Objects;
-import java.util.stream.Stream;
-
-public final class Image {
-    private final int width;
-    private final int height;
-
-    private final Color[][] data;
-
-    public Image(int width, int height) {
-        this.width = width;
-        this.height = height;
-
-        if (width <= 0) throw new IllegalArgumentException("width must be positive");
-        if (height <= 0) throw new IllegalArgumentException("height must be positive");
-
-        this.data = new Color[height][width];
-    }
-
-    public int width() {
-        return width;
-    }
-
-    public int height() {
-        return height;
-    }
-
-    public @NotNull Color get(int x, int y) {
-        Objects.checkIndex(x, width);
-        Objects.checkIndex(y, height);
-        return Objects.requireNonNullElse(this.data[y][x], Color.BLACK);
-    }
-
-    public @NotNull Stream<Color> pixels() {
-        return Arrays.stream(data).flatMap(Arrays::stream).map(c -> Objects.requireNonNullElse(c, Color.BLACK));
-    }
-
-    public void set(int x, int y, @NotNull Color color) {
-        Objects.checkIndex(x, width);
-        Objects.checkIndex(y, height);
-        this.data[y][x] = Objects.requireNonNull(color);
-    }
-
-    public void set(int x, int y, int red, int green, int blue) {
-        set(x, y, new Color(red, green, blue));
-    }
-
-    public void set(int x, int y, double r, double g, double b) {
-        set(x, y, new Color(r, g, b));
-    }
-}

src/main/java/eu/jonahbauer/raytracing/render/ImageFormat.java

@@ -1,130 +0,0 @@
-package eu.jonahbauer.raytracing.render;
-
-import org.jetbrains.annotations.NotNull;
-
-import java.io.*;
-import java.nio.charset.StandardCharsets;
-import java.nio.file.Files;
-import java.nio.file.Path;
-import java.util.zip.CRC32;
-import java.util.zip.CheckedOutputStream;
-import java.util.zip.DeflaterOutputStream;
-
-public enum ImageFormat {
-    PPM {
-        @Override
-        public void write(@NotNull Image image, @NotNull OutputStream out) throws IOException {
-            try (var writer = new BufferedWriter(new OutputStreamWriter(out, StandardCharsets.US_ASCII))) {
-                writer.write("P3\n");
-                writer.write(String.valueOf(image.width()));
-                writer.write(" ");
-                writer.write(String.valueOf(image.height()));
-                writer.write("\n255\n");
-
-                var it = image.pixels().iterator();
-                while (it.hasNext()) {
-                    var color = it.next();
-                    writer.write(String.valueOf(color.red()));
-                    writer.write(" ");
-                    writer.write(String.valueOf(color.green()));
-                    writer.write(" ");
-                    writer.write(String.valueOf(color.blue()));
-                    writer.write("\n");
-                }
-            }
-        }
-    },
-    PNG {
-        private static final byte[] MAGIC = new byte[] { (byte) 0x89, 0x50, 0x4E, 0x47, 0x0D, 0x0A, 0x1A, 0x0A };
-        private static final int IHDR_LENGTH = 13;
-        private static final int IHDR_TYPE = 0x49484452;
-        private static final int IDAT_TYPE = 0x49444154;
-        private static final int IEND_TYPE = 0x49454E44;
-
-        @Override
-        public void write(@NotNull Image image, @NotNull OutputStream out) throws IOException {
-            try (var data = new NoCloseDataOutputStream(out); var _ = data.closeable()) {
-                data.write(MAGIC);
-
-                writeIHDR(image, data);
-                writeIDAT(image, data);
-                writeIEND(image, data);
-            }
-        }
-
-        private void writeIHDR(@NotNull Image image, @NotNull DataOutputStream data) throws IOException {
-            data.writeInt(IHDR_LENGTH);
-            try (
-                    var crc = new CheckedOutputStream(data, new CRC32());
-                    var ihdr = new DataOutputStream(crc)
-            ) {
-                ihdr.writeInt(IHDR_TYPE);
-                ihdr.writeInt(image.width()); // image width
-                ihdr.writeInt(image.height()); // image height
-                ihdr.writeByte(8); // bit depth
-                ihdr.writeByte(2); // color type
-                ihdr.writeByte(0); // compression method
-                ihdr.writeByte(0); // filter method
-                ihdr.writeByte(0); // interlace method
-                ihdr.flush();
-                data.writeInt((int) crc.getChecksum().getValue());
-            }
-        }
-
-        private void writeIDAT(@NotNull Image image, @NotNull DataOutputStream data) throws IOException {
-            try (
-                    var baos = new ByteArrayOutputStream();
-                    var crc = new CheckedOutputStream(baos, new CRC32());
-                    var idat = new DataOutputStream(crc)
-            ) {
-                idat.writeInt(IDAT_TYPE);
-
-                try (var deflate = new DataOutputStream(new DeflaterOutputStream(idat))) {
-                    var pixels = image.pixels().iterator();
-                    for (int i = 0; pixels.hasNext(); i = (i + 1) % image.width()) {
-                        if (i == 0) deflate.writeByte(0); // filter type
-                        var pixel = pixels.next();
-                        deflate.writeByte(pixel.red());
-                        deflate.writeByte(pixel.green());
-                        deflate.writeByte(pixel.blue());
-                    }
-                }
-
-                var bytes = baos.toByteArray();
-                data.writeInt(bytes.length);
-                data.write(bytes);
-                data.writeInt((int) crc.getChecksum().getValue());
-            }
-        }
-
-        private void writeIEND(@NotNull Image image, @NotNull DataOutputStream data) throws IOException {
-            data.writeInt(0);
-            data.writeInt(IEND_TYPE);
-            data.writeInt(0);
-        }
-
-        private static class NoCloseDataOutputStream extends DataOutputStream {
-            public NoCloseDataOutputStream(OutputStream out) {
-                super(out);
-            }
-
-            @Override
-            public void close() {
-                // do nothing
-            }
-
-            public Closeable closeable() {
-                return super::close;
-            }
-        }
-    },
-    ;
-
-    public void write(@NotNull Image image, @NotNull Path path) throws IOException {
-        try (var out = Files.newOutputStream(path)) {
-            write(image, out);
-        }
-    }
-
-    public abstract void write(@NotNull Image image, @NotNull OutputStream out) throws IOException;
-}

src/main/java/eu/jonahbauer/raytracing/render/camera/Camera.java

@@ -0,0 +1,30 @@
+package eu.jonahbauer.raytracing.render.camera;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+public interface Camera {
+    /**
+     * {@return the width of this camera in pixels}
+     */
+    int getWidth();
+
+    /**
+     * {@return the height of this camera in pixels}
+     */
+    int getHeight();
+
+    /**
+     * Casts a ray through the given pixel.
+     * @param x the pixel x coordinate
+     * @param y the pixel y coordinate
+     * @param i the subpixel x coordinate
+     * @param j the subpixel y coordinate
+     * @param n the subpixel count (per side)
+     * @param random a random number generator
+     * @return a new ray
+     */
+    @NotNull Ray cast(int x, int y, int i, int j, int n, @NotNull RandomGenerator random);
+}

src/main/java/eu/jonahbauer/raytracing/render/camera/SimpleCamera.java

@@ -1,15 +1,14 @@
-package eu.jonahbauer.raytracing.render;
+package eu.jonahbauer.raytracing.render.camera;
 
-import eu.jonahbauer.raytracing.math.Range;
 import eu.jonahbauer.raytracing.math.Ray;
 import eu.jonahbauer.raytracing.math.Vec3;
-import eu.jonahbauer.raytracing.scene.Scene;
 import org.jetbrains.annotations.NotNull;
 import org.jetbrains.annotations.Nullable;
 
 import java.util.Objects;
+import java.util.random.RandomGenerator;
 
-public final class Camera {
+public final class SimpleCamera implements Camera {
     // image size
     private final int width;
     private final int height;
@@ -18,9 +17,6 @@ public final class Camera {
     private final @NotNull Vec3 origin;
 
     // rendering
-    private final int samplesPerPixel;
-    private final int maxDepth;
-    private final double gamma;
     private final double blurRadius;
 
     // internal properties
@@ -35,11 +31,11 @@ public final class Camera {
         return new Builder();
     }
 
-    public Camera() {
+    public static @NotNull Camera withDefaults() {
-        this(new Builder());
+        return new Builder().build();
     }
 
-    private Camera(@NotNull Builder builder) {
+    private SimpleCamera(@NotNull Builder builder) {
         this.width = builder.imageWidth;
         this.height = builder.imageHeight;
 
@@ -49,9 +45,6 @@ public final class Camera {
         this.origin = builder.position;
         var direction = (builder.direction == null ? builder.target.minus(builder.position).unit() : builder.direction);
 
-        this.samplesPerPixel = builder.samplePerPixel;
-        this.maxDepth = builder.maxDepth;
-        this.gamma = builder.gamma;
         this.blurRadius = Math.tan(0.5 * builder.blurAngle) * builder.focusDistance;
 
         // project direction the horizontal plane
@@ -67,87 +60,74 @@ public final class Camera {
 
         this.pixel00 = origin.plus(direction.times(builder.focusDistance))
                 .minus(u.times(0.5 * viewportWidth)).minus(v.times(0.5 * viewportHeight))
-                .plus(pixelU.div(2)).plus(pixelV.div(2));
+                .plus(pixelU.times(.5)).plus(pixelV.times(.5));
     }
 
-    public @NotNull Image render(@NotNull Scene scene) {
+    /**
-        var image = new Image(width, height);
+     * {@inheritDoc}
-
+     */
-        for (int y = 0; y < height; y++) {
+    public int getWidth() {
-            for (int x = 0; x < width; x++) {
+        return width;
-                var r = 0d;
-                var g = 0d;
-                var b = 0d;
-
-                for (int i = 0; i < samplesPerPixel; i++) {
-                    var ray = getRay(x, y);
-                    var color = getColor(scene, ray);
-                    r += color.r();
-                    g += color.g();
-                    b += color.b();
-                }
-
-                image.set(x, y, new Color(
-                        Math.pow(r / samplesPerPixel, 1 / gamma),
-                        Math.pow(g / samplesPerPixel, 1 / gamma),
-                        Math.pow(b / samplesPerPixel, 1 / gamma)
-                ));
-            }
-        }
-
-        return image;
     }
 
-    private @NotNull Ray getRay(int x, int y) {
+    /**
-        var origin = this.origin;
+     * {@inheritDoc}
-        if (blurRadius > 0) {
+     */
-            double bu, bv;
+    public int getHeight() {
-            do {
+        return height;
-                bu = 2 * Math.random() - 1;
-                bv = 2 * Math.random() - 1;
-            } while (bu * bu + bv * bv >= 1);
-
-            origin = origin.plus(u.times(blurRadius * bu)).plus(v.times(blurRadius * bv));
-        }
-
-        return new Ray(origin, getPixel(x, y).minus(origin));
     }
 
-    private @NotNull Vec3 getPixel(int x, int y) {
+    /**
+     * {@inheritDoc}
+     * @param x {@inheritDoc}
+     * @param y {@inheritDoc}
+     * @param i {@inheritDoc}
+     * @param j {@inheritDoc}
+     * @param n {@inheritDoc}
+     * @param random {@inheritDoc}
+     */
+    @Override
+    public @NotNull Ray cast(int x, int y, int i, int j, int n, @NotNull RandomGenerator random) {
         Objects.checkIndex(x, width);
         Objects.checkIndex(y, height);
 
-        double dx = x + Math.random() - 0.5;
+        var origin = getRayOrigin(random);
-        double dy = y + Math.random() - 0.5;
+        var target = getRayTarget(x, y, i, j, n, random);
-        return pixel00.plus(pixelU.times(dx)).plus(pixelV.times(dy));
+        return new Ray(origin, target.minus(origin));
     }
 
-    private @NotNull Color getColor(@NotNull Scene scene, @NotNull Ray ray) {
+    /**
-        return getColor(scene, ray, maxDepth);
+     * {@return the origin for a ray cast by this camera} The ray origin is randomized within a disk of
-    }
+     * radius {@link #blurRadius} centered on the camera position and perpendicular to the direction to simulate depth
+     * of field.
+     */
+    private @NotNull Vec3 getRayOrigin(@NotNull RandomGenerator random) {
+        if (blurRadius <= 0) return origin;
 
-    private @NotNull Color getColor(@NotNull Scene scene, @NotNull Ray ray, int depth) {
+        while (true) {
-        if (depth <= 0) return Color.BLACK;
+            var du = Math.fma(2, random.nextDouble(), -1);
+            var dv = Math.fma(2, random.nextDouble(), -1);
+            if (du * du + dv * dv >= 1) continue;
 
-        var optional = scene.hit(ray, new Range(0.001, Double.POSITIVE_INFINITY));
+            var ru = blurRadius * du;
-        if (optional.isPresent()) {
+            var rv = blurRadius * dv;
-            var hit = optional.get();
+
-            var material = hit.material();
+            return new Vec3(
-            return material.scatter(ray, hit)
+                    origin.x() + ru * u.x() + rv * v.x(),
-                    .map(scatter -> Color.multiply(scatter.attenuation(), getColor(scene, scatter.ray(), depth - 1)))
+                    origin.y() + ru * u.y() + rv * v.y(),
-                    .orElse(Color.BLACK);
+                    origin.z() + ru * u.z() + rv * v.z()
-        } else {
+            );
-            return getSkyboxColor(ray);
         }
     }
 
-    private static @NotNull Color getSkyboxColor(@NotNull Ray ray) {
+    /**
-        // altitude from -pi/2 to pi/2
+     * {@return the target vector for a ray through the given pixel} The position is randomized within the pixel.
-        var alt = Math.copySign(
+     */
-                Math.acos(ray.direction().withY(0).unit().times(ray.direction().unit())),
+    private @NotNull Vec3 getRayTarget(int x, int y, int i, int j, int n, @NotNull RandomGenerator random) {
-                ray.direction().y()
+        var factor = 1d / n;
-        );
+        var dx = x + Math.fma(factor, i + random.nextDouble(), -0.5);
-        return Color.lerp(Color.WHITE, Color.SKY, alt / Math.PI + 0.5);
+        var dy = y + Math.fma(factor, j + random.nextDouble(), -0.5);
+        return pixel00.plus(pixelU.times(dx)).plus(pixelV.times(dy));
     }
 
     public static class Builder {
@@ -163,9 +143,6 @@ public final class Camera {
         private double focusDistance = 10;
         private double blurAngle = 0.0;
 
-        private int samplePerPixel = 100;
-        private int maxDepth = 10;
-        private double gamma = 2.0;
 
         private Builder() {}
 
@@ -217,26 +194,8 @@ public final class Camera {
             return this;
         }
 
-        public @NotNull Builder withSamplesPerPixel(int samples) {
+        public @NotNull SimpleCamera build() {
-            if (samples <= 0) throw new IllegalArgumentException("samples must be positive");
+            return new SimpleCamera(this);
-            this.samplePerPixel = samples;
-            return this;
-        }
-
-        public @NotNull Builder withMaxDepth(int depth) {
-            if (depth <= 0) throw new IllegalArgumentException("depth must be positive");
-            this.maxDepth = depth;
-            return this;
-        }
-
-        public @NotNull Builder withGamma(double gamma) {
-            if (gamma <= 0 || !Double.isFinite(gamma)) throw new IllegalArgumentException("gamma must be positive");
-            this.gamma = gamma;
-            return this;
-        }
-
-        public @NotNull Camera build() {
-            return new Camera(this);
         }
 
     }

src/main/java/eu/jonahbauer/raytracing/render/canvas/Canvas.java

@@ -0,0 +1,37 @@
+package eu.jonahbauer.raytracing.render.canvas;
+
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import org.jetbrains.annotations.NotNull;
+
+public interface Canvas {
+    /**
+     * {@return the width of this canvas}
+     */
+    int getWidth();
+
+    /**
+     * {@return the height of this canvas}
+     */
+    int getHeight();
+
+    /**
+     * Adds a sample to this canvas
+     * @param x the pixel x coordinate
+     * @param y the pixel y coordinate
+     * @param n the index of the sample
+     * @param spectrum the sampled spectrum
+     * @param lambda the sampled wavelengths
+     */
+    void add(int x, int y, int n, @NotNull SampledSpectrum spectrum, @NotNull SampledWavelengths lambda);
+
+    /**
+     * {@return the color at a given pixel}
+     * @param x the pixel x coordinate
+     * @param y the pixel y coordinate
+     * @param cs the color space of the output
+     */
+    @NotNull ColorRGB getRGB(int x, int y, @NotNull ColorSpace cs);
+}

src/main/java/eu/jonahbauer/raytracing/render/canvas/LiveCanvas.java

@@ -0,0 +1,87 @@
+package eu.jonahbauer.raytracing.render.canvas;
+
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import org.jetbrains.annotations.NotNull;
+
+import javax.swing.*;
+import java.awt.*;
+import java.awt.event.WindowAdapter;
+import java.awt.event.WindowEvent;
+import java.awt.image.BufferedImage;
+
+public final class LiveCanvas implements Canvas {
+    private final @NotNull Canvas delegate;
+    private final @NotNull BufferedImage image;
+    private final @NotNull ColorSpace cs;
+
+    public LiveCanvas(@NotNull Canvas delegate, @NotNull ColorSpace cs) {
+        this.delegate = delegate;
+        this.image = new BufferedImage(delegate.getWidth(), delegate.getHeight(), BufferedImage.TYPE_INT_RGB);
+        this.cs = cs;
+    }
+
+    @Override
+    public int getWidth() {
+        return delegate.getWidth();
+    }
+
+    @Override
+    public int getHeight() {
+        return delegate.getHeight();
+    }
+
+    @Override
+    public void add(int x, int y, int n, @NotNull SampledSpectrum spectrum, @NotNull SampledWavelengths lambda) {
+        delegate.add(x, y, n, spectrum, lambda);
+        var color = cs.encode(delegate.getRGB(x, y, cs));
+        var rgb = color.red() << 16 | color.green() << 8 | color.blue();
+        image.setRGB(x, y, rgb);
+    }
+
+    @Override
+    public @NotNull ColorRGB getRGB(int x, int y, @NotNull ColorSpace cs) {
+        return delegate.getRGB(x, y, cs);
+    }
+
+    public @NotNull Thread preview() {
+        var frame = new JFrame();
+        frame.setDefaultCloseOperation(JFrame.DISPOSE_ON_CLOSE);
+        frame.setContentPane(new JPanel() {
+            {
+                setPreferredSize(new Dimension(image.getWidth(), image.getHeight()));
+            }
+
+            @Override
+            protected void paintComponent(Graphics g) {
+                g.drawImage(image, 0, 0, null);
+            }
+        });
+        frame.setResizable(false);
+        frame.pack();
+        frame.setVisible(true);
+
+        var update = Thread.ofVirtual().start(() -> {
+            while (!Thread.interrupted()) {
+                try {
+                    Thread.sleep(1000);
+                } catch (InterruptedException ex) {
+                    break;
+                }
+                frame.repaint();
+            }
+            frame.dispose();
+        });
+
+        frame.addWindowListener(new WindowAdapter() {
+            @Override
+            public void windowClosing(WindowEvent e) {
+                update.interrupt();
+            }
+        });
+
+        return update;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/canvas/RGBCanvas.java

@@ -0,0 +1,76 @@
+package eu.jonahbauer.raytracing.render.canvas;
+
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import org.jetbrains.annotations.NotNull;
+
+import java.awt.image.BufferedImage;
+import java.util.Objects;
+
+public final class RGBCanvas implements Canvas {
+    private final int width;
+    private final int height;
+
+    private final @NotNull ColorSpace cs;
+    private final @NotNull ColorRGB[][] data;
+
+    public RGBCanvas(int width, int height, @NotNull ColorSpace cs) {
+        this.width = width;
+        this.height = height;
+        this.cs = Objects.requireNonNull(cs);
+
+        if (width <= 0) throw new IllegalArgumentException("width must be positive");
+        if (height <= 0) throw new IllegalArgumentException("height must be positive");
+
+        this.data = new ColorRGB[height][width];
+    }
+
+    public RGBCanvas(@NotNull BufferedImage image, @NotNull ColorSpace cs) {
+        this(image.getWidth(), image.getHeight(), cs);
+
+        for (int y = 0; y < height; y++) {
+            for (int x = 0; x < width; x++) {
+                this.data[y][x] = new ColorRGB(image.getRGB(x, y));
+            }
+        }
+    }
+
+    @Override
+    public int getWidth() {
+        return width;
+    }
+
+    @Override
+    public int getHeight() {
+        return height;
+    }
+
+    @Override
+    public void add(int x, int y, int n, @NotNull SampledSpectrum spectrum, @NotNull SampledWavelengths lambda) {
+        assert x < width;
+        assert y < height;
+
+        var rgb = spectrum.toRGB(lambda, cs);
+        data[y][x] = ColorRGB.average(data[y][x], rgb, n);
+    }
+
+    @Override
+    public @NotNull ColorRGB getRGB(int x, int y, @NotNull ColorSpace cs) {
+        if (cs == this.cs) return get(x, y);
+        return cs.toRGB(this.cs.toXYZ(get(x, y)));
+    }
+
+    public @NotNull ColorRGB get(int x, int y) {
+        assert x < width;
+        assert y < height;
+        return Objects.requireNonNullElse(data[y][x], ColorRGB.BLACK);
+    }
+
+    public void set(int x, int y, @NotNull ColorRGB color) {
+        assert x < width;
+        assert y < height;
+        data[y][x] = Objects.requireNonNull(color);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/canvas/XYZCanvas.java

@@ -0,0 +1,74 @@
+package eu.jonahbauer.raytracing.render.canvas;
+
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.ColorXYZ;
+import org.jetbrains.annotations.NotNull;
+
+import java.awt.image.BufferedImage;
+import java.util.Objects;
+
+public final class XYZCanvas implements Canvas {
+    private final int width;
+    private final int height;
+
+    private final @NotNull ColorXYZ[][] data;
+
+    public XYZCanvas(int width, int height) {
+        this.width = width;
+        this.height = height;
+
+        if (width <= 0) throw new IllegalArgumentException("width must be positive");
+        if (height <= 0) throw new IllegalArgumentException("height must be positive");
+
+        this.data = new ColorXYZ[height][width];
+    }
+
+    public XYZCanvas(@NotNull BufferedImage image, @NotNull ColorSpace cs) {
+        this(image.getWidth(), image.getHeight());
+
+        for (int y = 0; y < height; y++) {
+            for (int x = 0; x < width; x++) {
+                data[y][x] = cs.toXYZ(new ColorRGB(image.getRGB(x, y)));
+            }
+        }
+    }
+
+    @Override
+    public int getWidth() {
+        return width;
+    }
+
+    @Override
+    public int getHeight() {
+        return height;
+    }
+
+    @Override
+    public void add(int x, int y, int n, @NotNull SampledSpectrum spectrum, @NotNull SampledWavelengths lambda) {
+        assert x < width;
+        assert y < height;
+
+        var xyz = spectrum.toXYZ(lambda);
+        data[y][x] = ColorXYZ.average(get(x, y), xyz, n);
+    }
+
+    public @NotNull ColorXYZ get(int x, int y) {
+        assert x < width;
+        assert y < height;
+        return Objects.requireNonNullElse(data[y][x], ColorXYZ.BLACK);
+    }
+
+    @Override
+    public @NotNull ColorRGB getRGB(int x, int y, @NotNull ColorSpace cs) {
+        return cs.toRGB(get(x, y));
+    }
+
+    public void set(int x, int y, @NotNull ColorXYZ color) {
+        assert x < width;
+        assert y < height;
+        data[y][x] = Objects.requireNonNull(color);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/Chromaticity.java

@@ -0,0 +1,9 @@
+package eu.jonahbauer.raytracing.render.color;
+
+/**
+ * A pair of chromaticity coordinates in the xyY color space
+ * @param x the x coordinate
+ * @param y the y coordinate
+ */
+public record Chromaticity(double x, double y) {
+}

src/main/java/eu/jonahbauer/raytracing/render/color/ColorRGB.java

@@ -0,0 +1,115 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.math.IVec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Random;
+
+public record ColorRGB(double r, double g, double b) implements IVec3<ColorRGB> {
+    public static final @NotNull ColorRGB BLACK = new ColorRGB(0.0, 0.0, 0.0);
+    public static final @NotNull ColorRGB WHITE = new ColorRGB(1.0, 1.0, 1.0);
+
+    public static @NotNull ColorRGB random(@NotNull Random random) {
+        return new ColorRGB(random.nextDouble(), random.nextDouble(), random.nextDouble());
+    }
+
+    public static @NotNull ColorRGB random(@NotNull Random random, double min, double max) {
+        var span = max - min;
+        return new ColorRGB(
+                Math.fma(random.nextDouble(), span, min),
+                Math.fma(random.nextDouble(), span, min),
+                Math.fma(random.nextDouble(), span, min)
+        );
+    }
+
+    public ColorRGB(int rgb) {
+        this(((rgb >> 16) & 0xFF) / 255d, ((rgb >> 8) & 0xFF) / 255d, (rgb & 0xFF) / 255d);
+    }
+
+    public ColorRGB {
+        assert Double.isFinite(r) : "r must be finite";
+        assert Double.isFinite(g) : "g must be finite";
+        assert Double.isFinite(b) : "b must be finite";
+    }
+
+    /*
+     * Math
+     */
+
+    public static @NotNull ColorRGB average(@NotNull ColorRGB current, @NotNull ColorRGB next, int index) {
+        return lerp(current, next, 1d / index);
+    }
+
+    public static @NotNull ColorRGB lerp(@NotNull ColorRGB a, @NotNull ColorRGB b, double t) {
+        if (t < 0) return a;
+        if (t > 1) return b;
+        return new ColorRGB(
+                Math.fma(t, b.r - a.r, a.r),
+                Math.fma(t, b.g - a.g, a.g),
+                Math.fma(t, b.b - a.b, a.b)
+        );
+    }
+
+    public static @NotNull ColorRGB fma(@NotNull ColorRGB a, @NotNull ColorRGB b, @NotNull ColorRGB c) {
+        return new ColorRGB(
+                Math.fma(a.r, b.r, c.r),
+                Math.fma(a.g, b.g, c.g),
+                Math.fma(a.b, b.b, c.b)
+        );
+    }
+
+    @Override
+    public @NotNull ColorRGB plus(@NotNull ColorRGB other) {
+        return new ColorRGB(r + other.r, g + other.g, b + other.b);
+    }
+
+    @Override
+    public @NotNull ColorRGB minus(@NotNull ColorRGB other) {
+        return new ColorRGB(r - other.r, g - other.g, b - other.b);
+    }
+
+    @Override
+    public @NotNull ColorRGB times(double d) {
+        return new ColorRGB(r * d, g * d, b * d);
+    }
+
+    @Override
+    public @NotNull ColorRGB times(@NotNull ColorRGB other) {
+        return new ColorRGB(r * other.r, g * other.g, b * other.b);
+    }
+
+    /*
+     * Accessors
+     */
+
+    public int red() {
+        return toInt(r);
+    }
+
+    public int green() {
+        return toInt(g);
+    }
+
+    public int blue() {
+        return toInt(b);
+    }
+
+    private static int toInt(double value) {
+        return Math.clamp((int) (255.99 * value), 0, 255);
+    }
+
+    @Override
+    public double component1() {
+        return r;
+    }
+
+    @Override
+    public double component2() {
+        return g;
+    }
+
+    @Override
+    public double component3() {
+        return b;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/ColorSpace.java

@@ -0,0 +1,164 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.math.Matrix3;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.*;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+
+/**
+ * An RGB color space.
+ */
+public final class ColorSpace {
+    private final @NotNull Chromaticity r;
+    private final @NotNull Chromaticity g;
+    private final @NotNull Chromaticity b;
+    private final @NotNull Chromaticity w;
+    private final @NotNull DenselySampledSpectrum illuminant;
+
+    private final @NotNull ColorXYZ R;
+    private final @NotNull ColorXYZ G;
+    private final @NotNull ColorXYZ B;
+    private final @NotNull ColorXYZ W;
+
+    private final @NotNull Matrix3 XYZfromRGB;
+    private final @NotNull Matrix3 RGBfromXYZ;
+    private final @NotNull SpectrumTable RGBtoSpectrumTable;
+    private final @NotNull TransferFunction transferFunction;
+
+    public ColorSpace(
+            @NotNull Chromaticity r, @NotNull Chromaticity g, @NotNull Chromaticity b,
+            @NotNull Spectrum illuminant, @NotNull SpectrumTable table, @NotNull TransferFunction transferFunction
+    ) {
+        this.r = Objects.requireNonNull(r, "r");
+        this.g = Objects.requireNonNull(g, "g");
+        this.b = Objects.requireNonNull(b, "b");
+        this.illuminant = new DenselySampledSpectrum(illuminant);
+        this.RGBtoSpectrumTable = table; // no null-check
+        this.transferFunction = transferFunction; // no null-check
+
+        this.W = illuminant.toXYZ();
+        this.w = W.xy();
+
+        this.R = new ColorXYZ(r);
+        this.G = new ColorXYZ(g);
+        this.B = new ColorXYZ(b);
+        var rgb = new Matrix3(
+                R.x(), G.x(), B.x(),
+                R.y(), G.y(), B.y(),
+                R.z(), G.z(), B.z()
+        );
+        var C = rgb.invert().times(W.toVec3());
+
+        this.XYZfromRGB = rgb.times(new Matrix3(C.x(), C.y(), C.z()));
+        this.RGBfromXYZ = XYZfromRGB.invert();
+    }
+
+    /*
+     * Conversions
+     */
+
+    public @NotNull ColorRGB toRGB(@NotNull ColorXYZ xyz) {
+        var out = RGBfromXYZ.times(xyz.toVec3());
+        return new ColorRGB(out.x(), out.y(), out.z());
+    }
+
+    public @NotNull ColorXYZ toXYZ(@NotNull ColorRGB rgb) {
+        var out = XYZfromRGB.times(rgb.toVec3());
+        return new ColorXYZ(out);
+    }
+
+    public @NotNull Vec3 toCIELab(@NotNull ColorRGB rgb) {
+        return toCIELab(toXYZ(rgb));
+    }
+
+    public @NotNull Vec3 toCIELab(@NotNull ColorXYZ xyz) {
+        return new Vec3(
+                116 * cieLabCbrt(xyz.y() / W.y()) - 16,
+                500 * (cieLabCbrt(xyz.x() / W.x()) - cieLabCbrt(xyz.y() / W.y())),
+                200 * (cieLabCbrt(xyz.y() / W.y()) - cieLabCbrt(xyz.z() / W.z()))
+        );
+    }
+
+    private static double cieLabCbrt(double x) {
+        var delta = 6.0 / 29.0;
+        if (x > delta * delta * delta) {
+            return Math.cbrt(x);
+        } else {
+            return x / (delta * delta * 3.0) + (4.0 / 29.0);
+        }
+    }
+
+    public @NotNull ColorRGB encode(@NotNull ColorRGB rgb) {
+        return transferFunction.encode(rgb);
+    }
+
+    public @NotNull ColorRGB decode(@NotNull ColorRGB rgb) {
+        return transferFunction.decode(rgb);
+    }
+
+    /*
+     * Spectrum
+     */
+
+    public @NotNull SigmoidPolynomial toPolynomial(@NotNull ColorRGB rgb) {
+        return RGBtoSpectrumTable.get(new ColorRGB(
+                Math.max(0, rgb.r()),
+                Math.max(0, rgb.g()),
+                Math.max(0, rgb.b())
+        ));
+    }
+
+    public @NotNull Spectrum toSpectrum(@NotNull ColorRGB rgb, @NotNull Spectrum.Type type) {
+        return switch (type) {
+            case ALBEDO -> new RGBAlbedoSpectrum(this, rgb);
+            case ILLUMINANT -> new RGBIlluminantSpectrum(this, rgb);
+            case UNBOUNDED -> new RGBUnboundedSpectrum(this, rgb);
+        };
+    }
+
+    public @NotNull Spectrum albedo(double r, double g, double b) {
+        return albedo(new ColorRGB(r, g, b));
+    }
+
+    public @NotNull Spectrum albedo(@NotNull ColorRGB rgb) {
+        return toSpectrum(rgb, Spectrum.Type.ALBEDO);
+    }
+
+    public @NotNull Spectrum illuminant(double intensity) {
+        return illuminant.scale(intensity);
+    }
+
+    public @NotNull Spectrum illuminant(double r, double g, double b) {
+        return illuminant(new ColorRGB(r, g, b));
+    }
+
+    public @NotNull Spectrum illuminant(@NotNull ColorRGB rgb) {
+        return toSpectrum(rgb, Spectrum.Type.ILLUMINANT);
+    }
+
+    /*
+     * Accessors
+     */
+
+    public @NotNull Chromaticity r() {
+        return r;
+    }
+
+    public @NotNull Chromaticity g() {
+        return g;
+    }
+
+    public @NotNull Chromaticity b() {
+        return b;
+    }
+
+    public @NotNull Chromaticity w() {
+        return w;
+    }
+
+    public @NotNull DenselySampledSpectrum illuminant() {
+        return illuminant;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/ColorSpaces.java

@@ -0,0 +1,44 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.IOException;
+import java.io.UncheckedIOException;
+import java.util.Objects;
+
+public final class ColorSpaces {
+    // Rec. ITU-R BT.709.3
+    public static final @NotNull ColorSpace sRGB = new ColorSpace(
+            new Chromaticity(0.6400, 0.3300),
+            new Chromaticity(0.3000, 0.6000),
+            new Chromaticity(0.1500, 0.0600),
+            Spectra.D65, read("sRGB_spectrum.bin"), TransferFunctions.sRGB
+    );
+    // P3-D65 (display)
+    public static final @NotNull ColorSpace DCI_P3 = new ColorSpace(
+            new Chromaticity(0.680, 0.320),
+            new Chromaticity(0.265, 0.690),
+            new Chromaticity(0.150, 0.060),
+            Spectra.D65, read("DCI_P3_spectrum.bin"), TransferFunctions.sRGB
+    );
+    // ITU-R Rec BT.2020
+    public static final @NotNull ColorSpace Rec2020 = new ColorSpace(
+            new Chromaticity(0.708, 0.292),
+            new Chromaticity(0.170, 0.797),
+            new Chromaticity(0.131, 0.046),
+            Spectra.D65, read("Rec2020_spectrum.bin"), null
+    );
+
+    private static @NotNull SpectrumTable read(@NotNull String name) {
+        try (var in = ColorSpaces.class.getResourceAsStream("/eu/jonahbauer/raytracing/colorspace/" + name)) {
+            return SpectrumTable.read(Objects.requireNonNull(in));
+        } catch (IOException e) {
+            throw new UncheckedIOException(e);
+        }
+    }
+
+    private ColorSpaces() {
+        throw new UnsupportedOperationException();
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/ColorXYZ.java

@@ -0,0 +1,111 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.math.IVec3;
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+/**
+ * A CIE XYZ color
+ */
+public record ColorXYZ(double x, double y, double z) implements IVec3<ColorXYZ> {
+    public static final double CIE_Y_INTEGRAL = 106.85689500000002;
+
+    public static final @NotNull ColorXYZ BLACK = new ColorXYZ(0, 0, 0);
+    public static final @NotNull ColorXYZ WHITE = new ColorXYZ(0, 1, 0);
+
+    public ColorXYZ(@NotNull Chromaticity chromaticity) {
+        this(chromaticity, 1);
+    }
+
+    public ColorXYZ(@NotNull Chromaticity chromaticity, double Y) {
+        this(
+                chromaticity.y() == 0 ? 0 : Y * chromaticity.x() / chromaticity.y(),
+                chromaticity.y() == 0 ? 0 : Y,
+                chromaticity.y() == 0 ? 0 : Y * (1 - chromaticity.x() - chromaticity.y()) / chromaticity.y()
+        );
+    }
+
+    public ColorXYZ(@NotNull Vec3 vec) {
+        this(vec.x(), vec.y(), vec.z());
+    }
+
+    public ColorXYZ {
+        assert Double.isFinite(x) : "x must be finite";
+        assert Double.isFinite(y) : "y must be finite";
+        assert Double.isFinite(z) : "z must be finite";
+    }
+
+    /*
+     * Math
+     */
+
+    public static @NotNull ColorXYZ average(@NotNull ColorXYZ a, @NotNull ColorXYZ b, int index) {
+        return lerp(a, b, 1d / index);
+    }
+
+    public static @NotNull ColorXYZ lerp(@NotNull ColorXYZ a, @NotNull ColorXYZ b, double t) {
+        if (t < 0) return a;
+        if (t > 1) return b;
+        return new ColorXYZ(
+                Math.fma(t, b.x - a.x, a.x),
+                Math.fma(t, b.y - a.y, a.y),
+                Math.fma(t, b.z - a.z, a.z)
+        );
+    }
+
+    public static @NotNull ColorXYZ fma(@NotNull ColorXYZ a, @NotNull ColorXYZ b, @NotNull ColorXYZ c) {
+        return new ColorXYZ(
+                Math.fma(a.x, b.x, c.x),
+                Math.fma(a.y, b.y, c.y),
+                Math.fma(a.z, b.z, c.z)
+        );
+    }
+
+    @Override
+    public @NotNull ColorXYZ plus(@NotNull ColorXYZ other) {
+        return new ColorXYZ(x + other.x, y + other.y, z + other.z);
+    }
+
+    @Override
+    public @NotNull ColorXYZ minus(@NotNull ColorXYZ other) {
+        return new ColorXYZ(x - other.x, y - other.y, z - other.z);
+    }
+
+    @Override
+    public @NotNull ColorXYZ times(@NotNull ColorXYZ other) {
+        return new ColorXYZ(x * other.x, y * other.y, z * other.z);
+    }
+
+    @Override
+    public @NotNull ColorXYZ times(double d) {
+        return new ColorXYZ(x * d, y * d, z * d);
+    }
+
+    /*
+     * Accessors
+     */
+
+    public double average() {
+        return (x + y + z) / 3;
+    }
+
+    public @NotNull Chromaticity xy() {
+        var factor = 1 / (x + y + z);
+        return new Chromaticity(factor * x, factor * y);
+    }
+
+    @Override
+    public double component1() {
+        return x;
+    }
+
+    @Override
+    public double component2() {
+        return y;
+    }
+
+    @Override
+    public double component3() {
+        return z;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/SigmoidPolynomial.java

@@ -0,0 +1,34 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+
+/**
+ * A function of the form {@code s(p(x))} where {@code p} is a polynomial of second degree and {@code s} is the sigmoid
+ * function <code>s(x) = 0.5 + x / (2 * sqrt(1 + x^2))</code>.
+ * <p>
+ * A function of this form is used to generate a {@link Spectrum} from an RGB value.
+ *
+ * @param c0 the coefficient of the quadratic term
+ * @param c1 the coefficient of the linear term
+ * @param c2 the coefficient of the constant term
+ */
+public record SigmoidPolynomial(double c0, double c1, double c2) {
+
+    public double get(double x) {
+        var p = Math.fma(Math.fma(c0, x, c1), x, c2);
+        if (!Double.isFinite(p)) return p > 0 ? 1 : 0;
+        return Math.fma(.5 * p, 1 / Math.sqrt(Math.fma(p, p, 1)), .5);
+    }
+
+    public double max() {
+        // evaluate at the edges
+        var result = Math.max(get(Spectrum.LAMBDA_MIN), get(Spectrum.LAMBDA_MAX));
+        var lambda = -c1 / (2 * c0);
+        if (lambda >= 360 && lambda <= 830) {
+            // evaluate at the vertex
+            return Math.max(result, get(lambda));
+        } else {
+            return result;
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/SpectrumTable.java

@@ -0,0 +1,142 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import org.jetbrains.annotations.NotNull;
+
+import java.io.*;
+import java.util.Arrays;
+
+/**
+ * A table of sigmoid polynomials used to convert between RGB values and spectra.
+ * <p>
+ * The rgb values are renormalized to xyz coordinates with {@code z} being the largest of the rgb components, and
+ * {@code x} and {@code y} being the other two rgb components divided by {@code z}. By this construction, {@code x},
+ * {@code y} and {@code z} are all in the range [0, 1] which allows for better use of the samples in a fixed grid.
+ * The {@code z} coordinate is additionally remapped using {@link #zNodes} to improve sampling at both ends of the scale.
+ * <p>
+ * The coefficients of the sigmoid functions are stored in a flattened five-dimensional array with indices as described
+ * in {@link #coefficients}.
+ */
+public final class SpectrumTable {
+    private final int resolution;
+
+    /**
+     * the remapped {@code z} values
+     */
+    private final double[] zNodes;
+
+    /**
+     * the stored coefficients as a flattened five-dimensional array with the following indices
+     * <ol>
+     *     <li>the component index of the biggest rgb component</li>
+     *     <li>the {@code z} coordinate</li>
+     *     <li>the {@code y} coordinate</li>
+     *     <li>the {@code x} coordinate</li>
+     *     <li>the coefficient index</li>
+     * </ol>
+     */
+    private final double[] coefficients;
+
+    public static void write(@NotNull SpectrumTable table, @NotNull OutputStream out) throws IOException {
+        var dos = new DataOutputStream(out);
+        dos.writeInt(table.resolution);
+        for (double z : table.zNodes) {
+            dos.writeDouble(z);
+        }
+        for (double c : table.coefficients) {
+            dos.writeDouble(c);
+        }
+        dos.flush();
+    }
+
+    public static @NotNull SpectrumTable read(@NotNull InputStream in) throws IOException {
+        var dis = new DataInputStream(in);
+        var resolution = dis.readInt();
+        var nodes = new double[resolution];
+        for (int i = 0; i < resolution; i++) {
+            nodes[i] = dis.readDouble();
+        }
+        var table = new double[3 * resolution * resolution * resolution * 3];
+        for (int i = 0; i < table.length; i++) {
+            table[i] = dis.readDouble();
+        }
+        return new SpectrumTable(resolution, nodes, table);
+    }
+
+    SpectrumTable(int resolution, double @NotNull[] zNodes, double[] coefficients) {
+        this.resolution = resolution;
+        this.zNodes = zNodes;
+        this.coefficients = coefficients;
+
+        // check input array lengths
+        if (zNodes.length != resolution) {
+            throw new IllegalArgumentException("length of zNodes must be equal to the RESOLUTION");
+        }
+        if (coefficients.length != 3 * resolution * resolution * resolution * 3) {
+            throw new IllegalArgumentException("coefficients length must be 3 * RESOLUTION * RESOLUTION * RESOLUTION * 3");
+        }
+        // check ascending zNodes
+        for (int i = 1; i < resolution; i++) {
+            if (zNodes[i - 1] >= zNodes[i]) {
+                throw new IllegalArgumentException("zNodes must be in increasing order");
+            }
+        }
+        if (zNodes[0] != 0.0 || zNodes[zNodes.length - 1] != 1.0) {
+            throw new IllegalArgumentException("zNodes must start with 0.0 and end with 1.0");
+        }
+    }
+
+    public @NotNull SigmoidPolynomial get(@NotNull ColorRGB color) {
+        // handle uniform rgb values
+        if (color.r() == color.g() && color.g() == color.b()) {
+            return new SigmoidPolynomial(0, 0, (color.r() - .5) / Math.sqrt(color.r() * (1 - color.r())));
+        }
+
+        // find maximum component and compute remapped component values
+        var max = color.r() > color.g()
+                ? (color.r() > color.b() ? 0 : 2)
+                : (color.g() > color.b() ? 1 : 2);
+
+        var z = color.get(max);
+        var x = color.get((max + 1) % 3) * (resolution - 1) / z;
+        var y = color.get((max + 2) % 3) * (resolution - 1) / z;
+
+        // compute integer indices and offsets for coefficient interpolation
+        int xi = Math.min((int) x, resolution - 2);
+        int yi = Math.min((int) y, resolution - 2);
+        int zi = Arrays.binarySearch(zNodes, z);
+        if (zi < 0) {
+            zi = -zi - 2;
+        } else if (zi > 0) {
+            zi = zi - 1;
+        }
+
+        var dx = x - xi;
+        var dy = y -yi;
+        var dz = (z - zNodes[zi]) / (zNodes[zi + 1] - zNodes[zi]);
+
+        // trilinearly interpolate sigmoid polynomial coefficients
+        var c = new double[3];
+        for (int i = 0; i < 3; i++) {
+            c[i] = lerp(dz,
+                        lerp(dy,
+                             lerp(dx, get(max, zi + 0, yi + 0, xi + 0, i), get(max, zi + 0, yi + 0, xi + 1, i)),
+                             lerp(dx, get(max, zi + 0, yi + 1, xi + 0, i), get(max, zi + 0, yi + 1, xi + 1, i))
+                        ),
+                        lerp(dy,
+                             lerp(dx, get(max, zi + 1, yi + 0, xi + 0, i), get(max, zi + 1, yi + 0, xi + 1, i)),
+                             lerp(dx, get(max, zi + 1, yi + 1, xi + 0, i), get(max, zi + 1, yi + 1, xi + 1, i))
+                        )
+            );
+        }
+
+        return new SigmoidPolynomial(c[0], c[1], c[2]);
+    }
+
+    private double get(int l, int z, int y, int x, int i) {
+        return coefficients[(((l * resolution + z) * resolution + y) * resolution + x) * 3 + i];
+    }
+
+    private static double lerp(double t, double a, double b) {
+        return Math.fma(t, b, Math.fma(-t, a, a));
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/SpectrumTableGenerator.java

@@ -0,0 +1,204 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import eu.jonahbauer.raytracing.math.Matrix3;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.util.Arrays;
+import java.util.Objects;
+import java.util.stream.IntStream;
+
+/**
+ * Generates a lookup table for RGB to spectrum conversion.
+ * <p>
+ * The spectrum for each RGB value is a {@link SigmoidPolynomial} with coefficients such that the round trip error
+ * from converting the RGB value to a spectrum and back is minimized.
+ * <p>
+ * <img src="doc-files/rgb2spectrum.png">
+ */
+public final class SpectrumTableGenerator {
+    private static final double EPSILON = 1e-4;
+    private static final int ITERATIONS = 15;
+
+    private final int resolution = 64;
+    private final @NotNull ColorSpace cs;
+
+    public static void main(String[] args) throws IOException {
+        var generator = new SpectrumTableGenerator(ColorSpaces.DCI_P3);
+        var table = generator.generate();
+
+        try (var out = Files.newOutputStream(Path.of("DCI_P3_spectrum.bin"))) {
+            SpectrumTable.write(table, out);
+        }
+    }
+
+    public SpectrumTableGenerator(@NotNull ColorSpace cs) {
+        this.cs = Objects.requireNonNull(cs);
+    }
+
+    public @NotNull SpectrumTable generate() {
+        var scale = new double[resolution];
+        for (int i = 0; i < scale.length; i++) {
+            var t = (double) i / (resolution - 1);
+            scale[i] = smoothstep(smoothstep(t));
+        }
+
+        var table = new double[3 * resolution * resolution * resolution * 3];
+
+        for (int l0 = 0; l0 < 3; l0++) {
+            var l = l0;
+            IntStream.range(0, resolution).parallel().forEach(i -> {
+                System.out.println("l = " + l + ", i = " + i);
+                var x = (double) i / (resolution - 1);
+                for (int j = 0; j < resolution; j++) {
+                    var y = (double) j / (resolution - 1);
+
+                    var start = resolution / 5;
+
+                    var c = new double[3];
+                    for (int k = start; k < resolution; k++) {
+                        var z = scale[k];
+                        var idx = ((((l * resolution + k) * resolution) + j) * resolution + i) * 3;
+                        var color = getColor(l, x, y, z);
+                        generate(color, c, table, idx);
+                    }
+
+                    Arrays.fill(c, 0);
+                    for (int k = start; k >= 0; --k) {
+                        var z = scale[k];
+                        var idx = ((((l * resolution + k) * resolution) + j) * resolution + i) * 3;
+                        var color = getColor(l, x, y, z);
+                        generate(color, c, table, idx);
+                    }
+                }
+            });
+        }
+        return new SpectrumTable(resolution, scale, table);
+    }
+
+    private void generate(@NotNull ColorRGB rgb, double @NotNull[] c, double @NotNull[] out, int offset) {
+        gaussNewton(rgb, c, ITERATIONS);
+        double c0 = 360.0, c1 = 1.0 / (830.0 - 360.0);
+        double A = c[0], B = c[1], C = c[2];
+        out[offset] = A * c1 * c1;
+        out[offset + 1] = B * c1 - 2 * A * c0 * c1 * c1;
+        out[offset + 2] = C - B * c0 * c1 + A * c0 * c0 * c1 * c1;
+    }
+
+    /**
+     * Use Gauss-Newton algorithm to calculate coefficients {@code c} of a {@link SigmoidPolynomial} such that the round
+     * trip error from converting the {@code rgb} value to a spectrum and back is minimized.
+     * @param rgb the input color
+     * @param c the coefficients, used as initial values and output
+     * @param it the number of iterations
+     */
+    private void gaussNewton(@NotNull ColorRGB rgb, double @NotNull[] c, int it) {
+        var bestQuality = Double.POSITIVE_INFINITY;
+        var bestCoefficients = new double[3];
+
+        for (int i = 0; i < it; ++i) {
+            var polynomial = new SigmoidPolynomial(c[0], c[1], c[2]);
+            var residual = getResidual(rgb, polynomial);
+            var jacobian = getJacobian(rgb, polynomial);
+
+            var delta = jacobian.decompose(1e-15).solve(residual);
+            for (int j = 0; j < 3; ++j) {
+                c[j] -= delta.get(j);
+            }
+
+            // catch runaway
+            double max = Math.max(Math.max(c[0], c[1]), c[2]);
+            if (max > 200) {
+                for (int j = 0; j < 3; ++j) {
+                    c[j] *= 200 / max;
+                }
+            }
+
+            var quality = residual.squared();
+            if (quality <= 1e-6) {
+                return;
+            } else if (quality < bestQuality) {
+                bestQuality = quality;
+                System.arraycopy(c, 0, bestCoefficients, 0, 3);
+            }
+        }
+
+        System.arraycopy(bestCoefficients, 0, c, 0, 3);
+    }
+
+    /**
+     * Calculates the Jacobian matrix of the {@code polynomial}.
+     */
+    private @NotNull Matrix3 getJacobian(@NotNull ColorRGB rgb, @NotNull SigmoidPolynomial polynomial) {
+        var jac = new double[3][3];
+
+        // central finite difference coefficients for first derivative with sixth-order accuracy
+        var factors = new double[] { -1d/60, 3d/20, -3d/4, 0, 3d/4, -3d/20,	1d/60 };
+
+        for (int i = 0; i < 3; i++) {
+            var derivative = Vec3.ZERO;
+            for (int d = - factors.length / 2, j = 0; j < factors.length; d++, j++) {
+                if (factors[j] == 0) continue;
+                var tmp = switch (i) {
+                    case 0 -> new SigmoidPolynomial(polynomial.c0() + d * EPSILON, polynomial.c1(), polynomial.c2());
+                    case 1 -> new SigmoidPolynomial(polynomial.c0(), polynomial.c1() + d * EPSILON, polynomial.c2());
+                    case 2 -> new SigmoidPolynomial(polynomial.c0(), polynomial.c1(), polynomial.c2() + d * EPSILON);
+                    default -> throw new AssertionError();
+                };
+                var r = getResidual(rgb, tmp);
+                derivative = Vec3.fma(factors[j], r, derivative);
+            }
+
+            for (int j = 0; j < 3; j++) {
+                jac[j][i] = derivative.get(j) / EPSILON;
+            }
+        }
+
+        return new Matrix3(
+                jac[0][0], jac[0][1], jac[0][2],
+                jac[1][0], jac[1][1], jac[1][2],
+                jac[2][0], jac[2][1], jac[2][2]
+        );
+    }
+
+    /**
+     * Calculates the difference between the RGB color and the result of converting the RGB color to a spectrum using
+     * the given coefficients, illuminating it with the color space's standard illuminant, and converting it back to an
+     * RBG color. The output is a vector in CIE Lab color space.
+     */
+    private @NotNull Vec3 getResidual(@NotNull ColorRGB rgb, @NotNull SigmoidPolynomial polynomial) {
+        var out = new SigmoidPolynomialSpectrum(polynomial, cs).toXYZ();
+        return cs.toCIELab(rgb).minus(cs.toCIELab(out));
+    }
+
+
+    private static double smoothstep(double x) {
+        return x * x * (3.0 - 2.0 * x);
+    }
+
+    private static @NotNull ColorRGB getColor(int l, double x, double y, double z) {
+        var rgb = new double[3];
+        rgb[l] = z;
+        rgb[(l + 1) % 3] = x * z;
+        rgb[(l + 2) % 3] = y * z;
+        return new ColorRGB(rgb[0], rgb[1], rgb[2]);
+    }
+
+    private record SigmoidPolynomialSpectrum(@NotNull SigmoidPolynomial polynomial, @NotNull ColorSpace cs) implements Spectrum {
+
+        @Override
+        public double max() {
+            return polynomial.max();
+        }
+
+        @Override
+        public double get(double lambda) {
+            var l = (lambda - Spectrum.LAMBDA_MIN) / (Spectrum.LAMBDA_MAX - Spectrum.LAMBDA_MIN);
+            return polynomial.get(l) * cs.illuminant().get(lambda);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/color/TransferFunction.java

@@ -0,0 +1,8 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import org.jetbrains.annotations.NotNull;
+
+public interface TransferFunction {
+    @NotNull ColorRGB decode(@NotNull ColorRGB rgb);
+    @NotNull ColorRGB encode(@NotNull ColorRGB rgb);
+}

src/main/java/eu/jonahbauer/raytracing/render/color/TransferFunctions.java

@@ -0,0 +1,50 @@
+package eu.jonahbauer.raytracing.render.color;
+
+import org.jetbrains.annotations.NotNull;
+
+public final class TransferFunctions {
+    public static final @NotNull TransferFunction sRGB = new ComponentTransferFunction() {
+        @Override
+        protected double encode(double value) {
+            if (value <= 0.0031308) return 12.92 * value;
+            return 1.055 * Math.pow(value, 1. / 2.4) - 0.055;
+        }
+
+        @Override
+        protected double decode(double value) {
+            if (value <= 0.04045) return value / 12.92;
+            return Math.pow((value + 0.055) / 1.055, 2.4d);
+        }
+    };
+
+    public static final @NotNull TransferFunction LINEAR = new TransferFunction() {
+        @Override
+        public @NotNull ColorRGB encode(@NotNull ColorRGB rgb) {
+            return rgb;
+        }
+
+        @Override
+        public @NotNull ColorRGB decode(@NotNull ColorRGB rgb) {
+            return rgb;
+        }
+    };
+
+    private TransferFunctions() {
+        throw new UnsupportedOperationException();
+    }
+
+    private abstract static class ComponentTransferFunction implements TransferFunction {
+        @Override
+        public final @NotNull ColorRGB decode(@NotNull ColorRGB rgb) {
+            return new ColorRGB(decode(rgb.r()), decode(rgb.g()), decode(rgb.b()));
+        }
+
+        @Override
+        public final @NotNull ColorRGB encode(@NotNull ColorRGB rgb) {
+            return new ColorRGB(encode(rgb.r()), encode(rgb.g()), encode(rgb.b()));
+        }
+
+        protected abstract double encode(double value);
+        protected abstract double decode(double value);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/image/ImageWriter.java

@@ -0,0 +1,19 @@
+package eu.jonahbauer.raytracing.render.image;
+
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.IOException;
+import java.io.OutputStream;
+import java.nio.file.Files;
+import java.nio.file.Path;
+
+public interface ImageWriter {
+    default void write(@NotNull Canvas image, @NotNull Path path) throws IOException {
+        try (var out = Files.newOutputStream(path)) {
+            write(image, out);
+        }
+    }
+
+    void write(@NotNull Canvas canvas, @NotNull OutputStream out) throws IOException;
+}

src/main/java/eu/jonahbauer/raytracing/render/image/PNGImageWriter.java

@@ -0,0 +1,107 @@
+package eu.jonahbauer.raytracing.render.image;
+
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.ColorSpaces;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.*;
+import java.util.Objects;
+import java.util.zip.CRC32;
+import java.util.zip.CheckedOutputStream;
+import java.util.zip.DeflaterOutputStream;
+
+public class PNGImageWriter implements ImageWriter {
+    public static final @NotNull PNGImageWriter sRGB = new PNGImageWriter(ColorSpaces.sRGB);
+
+    private static final byte[] MAGIC = new byte[] { (byte) 0x89, 0x50, 0x4E, 0x47, 0x0D, 0x0A, 0x1A, 0x0A };
+    private static final int IHDR_LENGTH = 13;
+    private static final int IHDR_TYPE = 0x49484452;
+    private static final int IDAT_TYPE = 0x49444154;
+    private static final int IEND_TYPE = 0x49454E44;
+    private static final int IEND_CRC = 0xAE426082;
+
+    private final @NotNull ColorSpace cs;
+
+    public PNGImageWriter(@NotNull ColorSpace cs) {
+        this.cs = Objects.requireNonNull(cs, "cs");
+    }
+
+    @Override
+    public void write(@NotNull Canvas image, @NotNull OutputStream out) throws IOException {
+        try (var data = new NoCloseDataOutputStream(out); var _ = data.closeable()) {
+            data.write(MAGIC);
+
+            writeIHDR(image, data);
+            writeIDAT(image, data);
+            writeIEND(image, data);
+        }
+    }
+
+    private void writeIHDR(@NotNull Canvas image, @NotNull DataOutputStream data) throws IOException {
+        data.writeInt(IHDR_LENGTH);
+        try (
+                var crc = new CheckedOutputStream(data, new CRC32());
+                var ihdr = new DataOutputStream(crc)
+        ) {
+            ihdr.writeInt(IHDR_TYPE);
+            ihdr.writeInt(image.getWidth()); // image width
+            ihdr.writeInt(image.getHeight()); // image height
+            ihdr.writeByte(8); // bit depth
+            ihdr.writeByte(2); // color type
+            ihdr.writeByte(0); // compression method
+            ihdr.writeByte(0); // filter method
+            ihdr.writeByte(0); // interlace method
+            ihdr.flush();
+            data.writeInt((int) crc.getChecksum().getValue());
+        }
+    }
+
+    private void writeIDAT(@NotNull Canvas image, @NotNull DataOutputStream data) throws IOException {
+        try (
+                var baos = new ByteArrayOutputStream();
+                var crc = new CheckedOutputStream(baos, new CRC32());
+                var idat = new DataOutputStream(crc)
+        ) {
+            idat.writeInt(IDAT_TYPE);
+
+            try (var deflate = new DataOutputStream(new DeflaterOutputStream(idat))) {
+                for (int y = 0; y < image.getHeight(); y++) {
+                    deflate.writeByte(0); // filter type
+                    for (int x = 0; x < image.getWidth(); x++) {
+                        var pixel = cs.encode(image.getRGB(x, y, cs));
+                        deflate.writeByte(pixel.red());
+                        deflate.writeByte(pixel.green());
+                        deflate.writeByte(pixel.blue());
+                    }
+                }
+            }
+
+            var bytes = baos.toByteArray();
+            data.writeInt(bytes.length - 4); // don't include type in length
+            data.write(bytes);
+            data.writeInt((int) crc.getChecksum().getValue());
+        }
+    }
+
+    private void writeIEND(@NotNull Canvas image, @NotNull DataOutputStream data) throws IOException {
+        data.writeInt(0);
+        data.writeInt(IEND_TYPE);
+        data.writeInt(IEND_CRC);
+    }
+
+    private static class NoCloseDataOutputStream extends DataOutputStream {
+        public NoCloseDataOutputStream(OutputStream out) {
+            super(out);
+        }
+
+        @Override
+        public void close() {
+            // do nothing
+        }
+
+        public Closeable closeable() {
+            return super::close;
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/image/PPMImageWriter.java

@@ -0,0 +1,46 @@
+package eu.jonahbauer.raytracing.render.image;
+
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.ColorSpaces;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.BufferedWriter;
+import java.io.IOException;
+import java.io.OutputStream;
+import java.io.OutputStreamWriter;
+import java.nio.charset.StandardCharsets;
+import java.util.Objects;
+
+public class PPMImageWriter implements ImageWriter {
+    public static final PPMImageWriter sRGB = new PPMImageWriter(ColorSpaces.sRGB);
+
+    private final @NotNull ColorSpace cs;
+
+    public PPMImageWriter(@NotNull ColorSpace cs) {
+        this.cs = Objects.requireNonNull(cs, "cs");
+    }
+
+    @Override
+    public void write(@NotNull Canvas image, @NotNull OutputStream out) throws IOException {
+        try (var writer = new BufferedWriter(new OutputStreamWriter(out, StandardCharsets.US_ASCII))) {
+            writer.write("P3\n");
+            writer.write(String.valueOf(image.getWidth()));
+            writer.write(" ");
+            writer.write(String.valueOf(image.getHeight()));
+            writer.write("\n255\n");
+
+            for (int y = 0; y < image.getHeight(); y++) {
+                for (int x = 0; x < image.getWidth(); x++) {
+                    var color = cs.encode(image.getRGB(x, y, cs));
+                    writer.write(String.valueOf(color.red()));
+                    writer.write(" ");
+                    writer.write(String.valueOf(color.green()));
+                    writer.write(" ");
+                    writer.write(String.valueOf(color.blue()));
+                    writer.write("\n");
+                }
+            }
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/DielectricMaterial.java

@@ -0,0 +1,85 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record DielectricMaterial(@NotNull RefractiveIndex ri, @NotNull Texture texture) implements Material {
+    public DielectricMaterial {
+        Objects.requireNonNull(ri, "ri");
+        Objects.requireNonNull(texture, "texture");
+    }
+
+    public DielectricMaterial(@NotNull RefractiveIndex ri) {
+        this(ri, Spectra.WHITE);
+    }
+
+    public DielectricMaterial(double ri) {
+        this(new ConstantRefractiveIndex(ri), Spectra.WHITE);
+    }
+
+    public DielectricMaterial(double ri, @NotNull Texture texture) {
+        this(new ConstantRefractiveIndex(ri), texture);
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        var ri = switch (this.ri) {
+            case ConstantRefractiveIndex(var x) -> x;
+            case RefractiveIndex x -> x.get(ray.lambda().collapse());
+        };
+        if (hit.isFrontFace()) ri = 1 / ri;
+
+        var cosTheta = Math.min(- ray.direction().unit().dot(hit.normal()), 1.0);
+        var reflectance = reflectance(cosTheta, ri);
+        var reflect = reflectance > random.nextDouble();
+
+        var newDirection = (reflect ? Optional.<Vec3>empty() : Vec3.refract(ray.direction(), hit.normal(), ri))
+                .orElseGet(() -> Vec3.reflect(ray.direction(), hit.normal()));
+
+        var attenuation = texture.get(hit);
+        return Optional.of(new SpecularScatterResult(attenuation, ray.with(hit, newDirection)));
+    }
+
+    private double reflectance(double cos, double ri) {
+        // use schlick's approximation for reflectance
+        var r0 = (1 - ri) / (1 + ri);
+        r0 = r0 * r0;
+        return r0 + (1 - r0) * (1 - cos) * (1 - cos) * (1 - cos) * (1 - cos) * (1 - cos);
+    }
+
+    @FunctionalInterface
+    public interface RefractiveIndex {
+        double get(double lambda);
+    }
+
+    public record ConstantRefractiveIndex(double ri) implements RefractiveIndex {
+
+        @Override
+        public double get(double lambda) {
+            return ri;
+        }
+    }
+
+    public record SellmeierRefractiveIndex(
+            double B1, double B2, double B3,
+            double C1, double C2, double C3
+    ) implements RefractiveIndex {
+
+        @Override
+        public double get(double lambda) {
+            var l2 = lambda * lambda * 1E-6; // square and convert to µm
+            var x = 1 + B1 * l2 / (l2 - C1)
+                      + B2 * l2 / (l2 - C2)
+                      + B3 * l2 / (l2 - C3);
+            return Math.sqrt(x);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/DiffuseLight.java

@@ -0,0 +1,27 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record DiffuseLight(@NotNull Texture texture) implements Material {
+    public DiffuseLight {
+        Objects.requireNonNull(texture, "texture");
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        return Optional.empty();
+    }
+
+    @Override
+    public @NotNull Spectrum emitted(@NotNull HitResult hit) {
+        return texture.get(hit);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/DirectionalMaterial.java

@@ -0,0 +1,69 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
+
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public final class DirectionalMaterial implements Material {
+    private final @Nullable Material front;
+    private final @Nullable Material back;
+
+    private final @NotNull Texture texture;
+
+    public DirectionalMaterial(@Nullable Material front, @Nullable Material back) {
+        if (front == null && back == null) throw new IllegalArgumentException("front and back must not both be null");
+        this.front = front;
+        this.back = back;
+        this.texture = new DirectionalTexture(
+                front != null ? front.texture() : null,
+                back != null ? back.texture() : null
+        );
+    }
+
+    @Override
+    public @NotNull Texture texture() {
+        return texture;
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        if (hit.isFrontFace()) {
+            if (front != null) return front.scatter(ray, hit, random);
+        } else {
+            if (back != null) return back.scatter(ray, hit, random);
+        }
+        // let the ray pass through without obstruction
+        return Optional.of(new SpecularScatterResult(Spectra.WHITE, ray.with(hit, ray.direction())));
+    }
+
+    @Override
+    public @NotNull Spectrum emitted(@NotNull HitResult hit) {
+        if (hit.isFrontFace()) {
+            if (front != null) return front.emitted(hit);
+        } else {
+            if (back != null) return back.emitted(hit);
+        }
+        return Material.super.emitted(hit);
+    }
+
+    private record DirectionalTexture(@Nullable Texture front, @Nullable Texture back) implements Texture {
+
+        @Override
+        public @NotNull Spectrum get(double u, double v, @NotNull Vec3 p) {
+            throw new UnsupportedOperationException();
+        }
+
+        @Override
+        public boolean isUVRequired() {
+            return front() != null && front().isUVRequired() || back() != null && back().isUVRequired();
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/IsotropicMaterial.java

@@ -0,0 +1,28 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.renderer.pdf.SphereProbabilityDensityFunction;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record IsotropicMaterial(@NotNull Spectrum albedo) implements Material {
+    public IsotropicMaterial {
+        Objects.requireNonNull(albedo, "albedo");
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        return Optional.of(new PdfScatterResult(albedo(), new SphereProbabilityDensityFunction()));
+    }
+
+    @Override
+    public @NotNull Texture texture() {
+        return albedo();
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/LambertianMaterial.java

@@ -0,0 +1,23 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.renderer.pdf.CosineProbabilityDensityFunction;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record LambertianMaterial(@NotNull Texture texture) implements Material {
+    public LambertianMaterial {
+        Objects.requireNonNull(texture, "texture");
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        var attenuation = texture.get(hit);
+        return Optional.of(new PdfScatterResult(attenuation, new CosineProbabilityDensityFunction(hit.normal())));
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/Material.java

@@ -0,0 +1,73 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.renderer.pdf.ProbabilityDensityFunction;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public interface Material {
+
+    /**
+     * {@return the texture associated with this material}
+     */
+    @NotNull Texture texture();
+
+    /**
+     * Scatters a light ray after it hit a surface.
+     * @param ray the incoming light ray
+     * @param hit information about the light ray hitting some object
+     * @param random a random number generator
+     * @return a {@code ScatterResult} if the ray is scattered or an {@linkplain Optional#empty() empty optional} if the
+     * ray is absorbed.
+     */
+    @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random);
+
+    /**
+     * {@return the color emitted for a given hit}
+     * @implSpec the default implementation returns {@linkplain ColorRGB#BLACK black}, i.e. no emission
+     */
+    default @NotNull Spectrum emitted(@NotNull HitResult hit) {
+        return Spectra.BLACK;
+    }
+
+    /**
+     * The result of a {@linkplain Material#scatter(Ray, HitResult, RandomGenerator) scattering operation}.
+     */
+    sealed interface ScatterResult {}
+
+    /**
+     * The result of a specular {@linkplain #scatter(Ray, HitResult, RandomGenerator) scattering operation}. A
+     * specular is a scattering operation with a very small number of possible scattered rays (like a
+     * perfect reflection which only has one possible scattered ray).
+     * @param attenuation the attenuation of the scattered light ray
+     * @param ray the scattered light ray
+     */
+    record SpecularScatterResult(@NotNull Spectrum attenuation, @NotNull Ray ray) implements ScatterResult {
+        public SpecularScatterResult {
+            Objects.requireNonNull(attenuation, "attenuation");
+            Objects.requireNonNull(ray, "ray");
+        }
+    }
+
+    /**
+     * The result of a probability density function based
+     * {@linkplain #scatter(Ray, HitResult, RandomGenerator) scattering operation}. A probability density function
+     * based scattering operation uses a probability density function to determine the scatter direction.
+     * @param attenuation the attenuation of the scattered light ray
+     * @param pdf the probability density function
+     */
+    record PdfScatterResult(@NotNull Spectrum attenuation, @NotNull ProbabilityDensityFunction pdf) implements ScatterResult {
+        public PdfScatterResult {
+            Objects.requireNonNull(attenuation, "attenuation");
+            Objects.requireNonNull(pdf, "pdf");
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/Materials.java

@@ -0,0 +1,25 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.render.color.ColorSpaces;
+import eu.jonahbauer.raytracing.render.material.DielectricMaterial.SellmeierRefractiveIndex;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.render.texture.CheckerTexture;
+import org.jetbrains.annotations.NotNull;
+
+public final class Materials {
+    public static final @NotNull Material GLASS = new DielectricMaterial(new SellmeierRefractiveIndex(
+            1.0361212, 0.231792344, 1.01046945,
+            6.00069867E-3, 2.00179144E-2, 1.03560653E2
+    ));
+
+    public static final @NotNull Material MIRROR = new MetallicMaterial(ColorSpaces.sRGB.albedo(0.7, 0.7, 0.7));
+
+    public static final @NotNull Material DEBUG = new DirectionalMaterial(
+            new LambertianMaterial(new CheckerTexture(50.0, ColorSpaces.sRGB.albedo(1.0, 0.0, 1.0), Spectra.BLACK)),
+            new LambertianMaterial(new CheckerTexture(50.0, ColorSpaces.sRGB.albedo(1.0, 1.0, 1.0), Spectra.BLACK))
+    );
+
+    private Materials() {
+        throw new UnsupportedOperationException();
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/material/MetallicMaterial.java

@@ -0,0 +1,33 @@
+package eu.jonahbauer.raytracing.render.material;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record MetallicMaterial(@NotNull Texture texture, double fuzz) implements Material {
+
+    public MetallicMaterial(@NotNull Texture texture) {
+        this(texture, 0);
+    }
+
+    public MetallicMaterial {
+        Objects.requireNonNull(texture, "texture");
+        if (fuzz < 0 || !Double.isFinite(fuzz)) throw new IllegalArgumentException("fuzz must be non-negative");
+    }
+
+    @Override
+    public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
+        var newDirection = Vec3.reflect(ray.direction(), hit.normal());
+        if (fuzz > 0) {
+            newDirection = Vec3.fma(fuzz, Vec3.random(random), newDirection.unit());
+        }
+        var attenuation = texture.get(hit);
+        return Optional.of(new SpecularScatterResult(attenuation, ray.with(hit, newDirection)));
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/Renderer.java

@@ -0,0 +1,20 @@
+package eu.jonahbauer.raytracing.render.renderer;
+
+import eu.jonahbauer.raytracing.render.camera.Camera;
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
+import eu.jonahbauer.raytracing.render.canvas.XYZCanvas;
+import eu.jonahbauer.raytracing.scene.Scene;
+import org.jetbrains.annotations.NotNull;
+
+public interface Renderer {
+    default @NotNull Canvas render(@NotNull Camera camera, @NotNull Scene scene) {
+        var image = new XYZCanvas(camera.getWidth(), camera.getHeight());
+        render(camera, scene, image);
+        return image;
+    }
+
+    /**
+     * Renders the {@code scene} as seen by the {@code camera} to the {@code canvas}.
+     */
+    void render(@NotNull Camera camera, @NotNull Scene scene, @NotNull Canvas canvas);
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/SimpleRenderer.java

@@ -0,0 +1,260 @@
+package eu.jonahbauer.raytracing.render.renderer;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.render.renderer.pdf.TargetingProbabilityDensityFunction;
+import eu.jonahbauer.raytracing.render.renderer.pdf.MixtureProbabilityDensityFunction;
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.camera.Camera;
+import eu.jonahbauer.raytracing.render.canvas.Canvas;
+import eu.jonahbauer.raytracing.render.spectrum.SampledWavelengths;
+import eu.jonahbauer.raytracing.scene.Scene;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Random;
+import java.util.SplittableRandom;
+import java.util.random.RandomGenerator;
+import java.util.stream.IntStream;
+
+import static eu.jonahbauer.raytracing.Main.DEBUG;
+
+public final class SimpleRenderer implements Renderer {
+    private final int sqrtSamplesPerPixel;
+    private final int maxDepth;
+
+    private final int spectralSamples;
+    private final SampledSpectrum black;
+    private final SampledSpectrum white;
+
+    private final boolean parallel;
+    private final boolean iterative;
+
+    public static @NotNull Builder builder() {
+        return new Builder();
+    }
+
+    public static @NotNull Renderer withDefaults() {
+        return new Builder().build();
+    }
+
+    private SimpleRenderer(@NotNull Builder builder) {
+        this.sqrtSamplesPerPixel = (int) Math.sqrt(builder.samplesPerPixel);
+        this.maxDepth = builder.maxDepth;
+
+        this.spectralSamples = builder.spectralSamples;
+        this.black = new SampledSpectrum(spectralSamples, 0);
+        this.white = new SampledSpectrum(spectralSamples, 1);
+
+        this.parallel = builder.parallel;
+        this.iterative = builder.iterative;
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public void render(@NotNull Camera camera, @NotNull Scene scene, @NotNull Canvas canvas) {
+        if (canvas.getWidth() != camera.getWidth() || canvas.getHeight() != camera.getHeight()) {
+            throw new IllegalArgumentException("sizes of camera and canvas are different");
+        }
+
+        if (iterative) {
+            renderIterative(camera, scene, canvas);
+        } else {
+            renderNonIterative(camera, scene, canvas);
+        }
+    }
+
+    /**
+     * Renders the {@code scene} as seen by the {@code camera} to the {@code canvas}, taking one sample per pixel at
+     * a time and updating the canvas after each sample.
+     */
+    private void renderIterative(@NotNull Camera camera, @NotNull Scene scene, @NotNull Canvas canvas) {
+        var random = new Random(0);
+
+        // render one sample after the other
+        int i = 0;
+        for (int sj = 0; sj < sqrtSamplesPerPixel; sj++) {
+            for (int si = 0; si < sqrtSamplesPerPixel; si++) {
+                var sample = ++i;
+                var sif = si;
+                var sjf = sj;
+                getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
+                    for (int x = 0; x < camera.getWidth(); x++) {
+                        var ray = camera.cast(x, y, sif, sjf, sqrtSamplesPerPixel, random);
+                        var c = getColor(scene, ray, random);
+                        canvas.add(x, y, sample, c, ray.lambda());
+                    }
+                });
+            }
+        }
+    }
+
+    /**
+     * Renders the {@code scene} as seen by the {@code camera} to the {@code canvas}, taking some amount of samples
+     * per pixel and updating the canvas after each pixel.
+     */
+    private void renderNonIterative(@NotNull Camera camera, @NotNull Scene scene, @NotNull Canvas canvas) {
+        var splittable = new SplittableRandom(0);
+        // render one pixel after the other
+        getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
+            var random = splittable.split();
+            for (int x = 0; x < camera.getWidth(); x++) {
+                int i = 0;
+                for (int sj = 0; sj < sqrtSamplesPerPixel; sj++) {
+                    for (int si = 0; si < sqrtSamplesPerPixel; si++) {
+                        var lambda = SampledWavelengths.uniform(random.nextDouble(), spectralSamples);
+                        var ray = camera.cast(x, y, si, sj, sqrtSamplesPerPixel, random).with(lambda);
+                        if (DEBUG) {
+                            System.out.println("Casting ray " + ray + " through pixel (" + x + "," + y + ") at subpixel (" + si + "," + sj + ")...");
+                        }
+                        var c = getColor(scene, ray, random);
+                        canvas.add(x, y, ++i, c, ray.lambda());
+                    }
+                }
+            }
+        });
+    }
+
+    /**
+     * {@return the color of the given ray in the given scene}
+     */
+    private @NotNull SampledSpectrum getColor(@NotNull Scene scene, @NotNull Ray ray, @NotNull RandomGenerator random) {
+        return getColor0(scene, ray, maxDepth, random);
+    }
+
+    private @NotNull SampledSpectrum getColor0(@NotNull Scene scene, @NotNull Ray ray, int depth, @NotNull RandomGenerator random) {
+        var color = black;
+        var attenuation = white;
+
+        while (depth-- > 0) {
+            var optional = scene.hit(ray, random);
+            if (optional.isEmpty()) {
+                var background = scene.getBackgroundColor(ray);
+                color = SampledSpectrum.fma(attenuation, background, color);
+                if (DEBUG) {
+                    System.out.println("  Hit background: " + background);
+                }
+                break;
+            }
+
+            var hit = optional.get();
+            if (DEBUG) {
+                System.out.println("  Hit " + hit.target() + " at t=" + hit.t() + " (" + hit.position() + ")");
+            }
+            var material = hit.material();
+            var emitted = material.emitted(hit).sample(ray.lambda());
+            if (DEBUG && !black.equals(emitted)) {
+                System.out.println("  Emitted: " + emitted);
+            }
+
+            var result = material.scatter(ray, hit, random);
+            color = SampledSpectrum.fma(attenuation, emitted, color);
+
+            if (result.isEmpty()) {
+                if (DEBUG) {
+                    System.out.println("  Absorbed");
+                }
+                break;
+            }
+
+            switch (result.get()) {
+                case Material.SpecularScatterResult(var a, var scattered) -> {
+                    attenuation = attenuation.times(a.sample(ray.lambda()));
+                    ray = scattered;
+
+                    if (DEBUG) {
+                        System.out.println("  Specular scattering with albedo " + a);
+                    }
+                }
+                case Material.PdfScatterResult(var a, var pdf) -> {
+                    if (scene.getTargets() == null) {
+                        attenuation = attenuation.times(a.sample(ray.lambda()));
+                        ray = new Ray(hit.position(), pdf.generate(random), ray.lambda());
+
+                        if (DEBUG) {
+                            System.out.println("  Pdf scattering with albedo " + a);
+                        }
+                    } else {
+                        var mixed = new MixtureProbabilityDensityFunction(new TargetingProbabilityDensityFunction(hit.position(), scene.getTargets()), pdf, 0.5);
+                        var direction = mixed.generate(random).unit();
+
+                        var idealPdf = pdf.value(direction);
+                        var actualPdf = mixed.value(direction);
+                        if (actualPdf == 0) break; // when actualPdf is 0, the ray should have never been generated by mixed.generate
+
+                        var factor = idealPdf / actualPdf;
+
+                        attenuation = attenuation.times(a.sample(ray.lambda()).times(factor));
+                        ray = new Ray(hit.position(), direction, ray.lambda());
+
+                        if (DEBUG) {
+                            System.out.println("  Pdf scattering with albedo " + a + " and factor " + factor);
+                        }
+                    }
+                }
+            }
+
+            if (DEBUG) {
+                System.out.println("  Combined color is " + color);
+                System.out.println("  Combined attenuation is " + attenuation);
+                System.out.println("  New ray is " + ray);
+            }
+        }
+
+        if (DEBUG) {
+            System.out.println("  Final color is " + color);
+        }
+
+        return color;
+    }
+
+    /**
+     * {@return a stream of the pixels in a canvas with the given size} The pixels {@code x} and {@code y} coordinate
+     * are encoded in the longs lower and upper 32 bits respectively.
+     */
+    private static @NotNull IntStream getScanlineStream(int height, boolean parallel) {
+        var stream = IntStream.range(0, height).map(i -> height - i - 1);
+        return parallel ? stream.parallel() : stream;
+    }
+
+    public static class Builder {
+        private int samplesPerPixel = 100;
+        private int maxDepth = 10;
+        private int spectralSamples = 4;
+        private boolean parallel = true;
+        private boolean iterative = false;
+
+        public @NotNull Builder withSamplesPerPixel(int samples) {
+            if (samples <= 0) throw new IllegalArgumentException("samples must be positive");
+            this.samplesPerPixel = samples;
+            return this;
+        }
+
+        public @NotNull Builder withMaxDepth(int depth) {
+            if (depth <= 0) throw new IllegalArgumentException("depth must be positive");
+            this.maxDepth = depth;
+            return this;
+        }
+
+        public @NotNull Builder withSpectralSamples(int samples) {
+            if (samples <= 0) throw new IllegalArgumentException("samples must be positive");
+            this.spectralSamples = samples;
+            return this;
+        }
+
+        public @NotNull Builder withParallel(boolean parallel) {
+            this.parallel = parallel;
+            return this;
+        }
+
+        public @NotNull Builder withIterative(boolean iterative) {
+            this.iterative = iterative;
+            return this;
+        }
+
+        public @NotNull SimpleRenderer build() {
+            return new SimpleRenderer(this);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/pdf/CosineProbabilityDensityFunction.java

@@ -0,0 +1,27 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.random.RandomGenerator;
+
+public record CosineProbabilityDensityFunction(@NotNull Vec3 normal) implements ProbabilityDensityFunction {
+
+    public CosineProbabilityDensityFunction {
+        Objects.requireNonNull(normal, "normal");
+        normal = normal.unit();
+    }
+
+    @Override
+    public double value(@NotNull Vec3 direction) {
+        var cos = normal.dot(direction);
+        return Math.max(0, cos / Math.PI);
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        var out = normal().plus(Vec3.random(random));
+        return out.isNearZero() ? normal() : out;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/pdf/MixtureProbabilityDensityFunction.java

@@ -0,0 +1,46 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.random.RandomGenerator;
+
+/**
+ * Mixes between two probability density functions (pdf) using a weight. When the weight is closer to zero, the
+ * influence of the second pdf is stronger. When the weight is closer to one, the influence of the first pdf is stronger.
+ * @param a the first probability density function
+ * @param b the second probability density function
+ * @param weight a weight in the range [0, 1]
+ */
+public record MixtureProbabilityDensityFunction(
+        @NotNull ProbabilityDensityFunction a,
+        @NotNull ProbabilityDensityFunction b,
+        double weight
+) implements ProbabilityDensityFunction {
+    public MixtureProbabilityDensityFunction(@NotNull ProbabilityDensityFunction a, @NotNull ProbabilityDensityFunction b) {
+        this(a, b, 0.5);
+    }
+
+    public MixtureProbabilityDensityFunction {
+        Objects.requireNonNull(a);
+        Objects.requireNonNull(b);
+        weight = Math.clamp(weight, 0, 1);
+    }
+
+    @Override
+    public double value(@NotNull Vec3 direction) {
+        var v = a.value(direction);
+        var w = b.value(direction);
+        return Math.fma(weight, v, Math.fma(-weight, w, w));
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        if (random.nextDouble() < weight) {
+            return a.generate(random);
+        } else {
+            return b.generate(random);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/pdf/ProbabilityDensityFunction.java

@@ -0,0 +1,25 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+/**
+ * A probability density function used for sampling random directions when scattering a ray.
+ */
+public interface ProbabilityDensityFunction {
+
+    /**
+     * {@return the value of this probability density function at the given point}
+     * @param direction the direction
+     */
+    double value(@NotNull Vec3 direction);
+
+    /**
+     * Generates a random direction based on this probability density function.
+     * @param random a random number generator
+     * @return the random direction
+     */
+    @NotNull Vec3 generate(@NotNull RandomGenerator random);
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/pdf/SphereProbabilityDensityFunction.java

@@ -0,0 +1,22 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+/**
+ * A probability density function sampling the sphere uniformly.
+ */
+public record SphereProbabilityDensityFunction() implements ProbabilityDensityFunction {
+
+    @Override
+    public double value(@NotNull Vec3 direction) {
+        return 1 / (4 * Math.PI);
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        return Vec3.random(random);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/renderer/pdf/TargetingProbabilityDensityFunction.java

@@ -0,0 +1,41 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.Target;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Objects;
+import java.util.random.RandomGenerator;
+
+/**
+ * A probability density function targeting a target.
+ * @see Target
+ */
+public final class TargetingProbabilityDensityFunction implements ProbabilityDensityFunction {
+    private final @NotNull Vec3 origin;
+    private final @NotNull List<@NotNull Target> targets;
+
+    public TargetingProbabilityDensityFunction(@NotNull Vec3 origin, @NotNull List<@NotNull Target> targets) {
+        this.origin = Objects.requireNonNull(origin, "origin");
+        this.targets = new ArrayList<>(targets);
+    }
+
+    @Override
+    public double value(@NotNull Vec3 direction) {
+        var weight = 1d / targets.size();
+        var sum = 0.0;
+
+        for (var target : targets) {
+            sum = Math.fma(weight, target.getProbabilityDensity(origin, direction), sum);
+        }
+
+        return sum;
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        return targets.get(random.nextInt(targets.size())).getTargetingDirection(origin, random);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/BlackbodySpectrum.java

@@ -0,0 +1,45 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+public final class BlackbodySpectrum implements Spectrum {
+    /**
+     * the speed of light in m/s
+     */
+    private static final double c = 299792458d;
+
+    /**
+     * the planck constant in m^2*kg/s
+     */
+    private static final double h = 6.62607015E-34;
+
+    /**
+     * the boltzmann constant in m^2*kg/s^2/K
+     */
+    private static final double k = 1.380649E-23;
+
+    /**
+     * wien's displacement constant in m*K
+     */
+    private static final double b = 2.897771995E-3;
+
+    private final double T;
+    private final double factor;
+
+    public BlackbodySpectrum(double T) {
+        if (T < 0) throw new IllegalArgumentException("T must be non-negative");
+        this.T = T;
+        this.factor = 1 / get(b / T);
+    }
+
+    @Override
+    public double max() {
+        return 1;
+    }
+
+    @Override
+    public double get(double lambda) {
+        lambda *= 1E-9;
+        var l2 = lambda * lambda;
+        var x = h * c / (lambda * k * T);
+        return 2 * h * c * c / (l2 * l2 * lambda) / (Math.exp(x) - 1) * factor;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/ConstantSpectrum.java

@@ -0,0 +1,17 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+/**
+ * A constant spectrum.
+ * @param c the constant value
+ */
+public record ConstantSpectrum(double c) implements Spectrum {
+    @Override
+    public double max() {
+        return c;
+    }
+
+    @Override
+    public double get(double lambda) {
+        return c;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/DenselySampledSpectrum.java

@@ -0,0 +1,51 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Arrays;
+
+/**
+ * A spectrum sampled in one nanometer intervals.
+ */
+public final class DenselySampledSpectrum implements Spectrum {
+    private final double[] samples;
+    private final int min;
+
+    private final double max;
+
+    public DenselySampledSpectrum(@NotNull Spectrum spectrum) {
+        this(spectrum, LAMBDA_MIN, LAMBDA_MAX);
+    }
+
+    public DenselySampledSpectrum(@NotNull Spectrum spectrum, int min, int max) {
+        if (max - min + 1 <= 0) throw new IllegalArgumentException("samples must not be empty");
+        this.samples = new double[max - min + 1];
+        var maxValue = 0d;
+        for (int lambda = min, i = 0; lambda <= max; lambda++, i++) {
+            var sample = spectrum.get(lambda);
+            if (sample > maxValue) maxValue = sample;
+            this.samples[i] = sample;
+        }
+        this.min = min;
+        this.max = maxValue;
+    }
+
+    public DenselySampledSpectrum(double @NotNull[] samples, int lambdaMin) {
+        if (samples.length == 0) throw new IllegalArgumentException("samples must not be empty");
+        this.samples = Arrays.copyOf(samples, samples.length);
+        this.min = lambdaMin;
+        this.max = Arrays.stream(this.samples).max().orElseThrow();
+    }
+
+    @Override
+    public double max() {
+        return max;
+    }
+
+    @Override
+    public double get(double lambda) {
+        int offset = (int) Math.round(lambda) - min;
+        if (offset < 0 || offset >= samples.length) return 0;
+        return samples[offset];
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/PiecewiseLinearSpectrum.java

@@ -0,0 +1,48 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import java.util.Arrays;
+
+public final class PiecewiseLinearSpectrum implements Spectrum {
+    private final double[] lambdas;
+    private final double[] values;
+    private final double max;
+
+    public PiecewiseLinearSpectrum(double[] lambdas, double[] values) {
+        if (lambdas.length != values.length) {
+            throw new IllegalArgumentException("lambdas and values must have the same length");
+        }
+
+        this.lambdas = Arrays.copyOf(lambdas, lambdas.length);
+        this.values = Arrays.copyOf(values, values.length);
+
+        var max = 0d;
+        for (int i = 1; i < this.lambdas.length; i++) {
+            if (this.lambdas[i] <= this.lambdas[i - 1]) {
+                throw new IllegalArgumentException("lambdas must be in increasing order");
+            }
+            if (this.values[i] < 0) {
+                throw new IllegalArgumentException("values must be non-negative");
+            } else if (this.values[i] > max) {
+                max = this.values[i];
+            }
+        }
+        this.max = max;
+    }
+
+    @Override
+    public double max() {
+        return max;
+    }
+
+    @Override
+    public double get(double lambda) {
+        if (lambdas.length == 0 || lambda < lambdas[0] || lambda > lambdas[lambdas.length - 1]) return 0;
+        if (lambda == lambdas[lambdas.length - 1]) return values[values.length - 1];
+
+        var i = Arrays.binarySearch(lambdas, lambda);
+        if (i < 0) i = -i - 1;
+
+        var t = (lambda - lambdas[i]) / (lambdas[i + 1] - lambdas[i]);
+        return (1 - t) * values[i] + t * values[i + 1];
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/RGBAlbedoSpectrum.java

@@ -0,0 +1,27 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.SigmoidPolynomial;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import org.jetbrains.annotations.NotNull;
+
+public final class RGBAlbedoSpectrum implements Spectrum {
+    private final @NotNull SigmoidPolynomial polynomial;
+
+    public RGBAlbedoSpectrum(@NotNull ColorSpace cs, @NotNull ColorRGB rgb) {
+        if (rgb.r() < 0 || rgb.r() > 1 || rgb.g() < 0 || rgb.g() > 1 || rgb.b() < 0 || rgb.b() > 1) {
+            throw new IllegalArgumentException();
+        }
+        this.polynomial = cs.toPolynomial(rgb);
+    }
+
+    @Override
+    public double max() {
+        return polynomial.max();
+    }
+
+    @Override
+    public double get(double lambda) {
+        return polynomial.get(lambda);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/RGBIlluminantSpectrum.java

@@ -0,0 +1,36 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.SigmoidPolynomial;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import org.jetbrains.annotations.NotNull;
+
+/**
+ * A spectrum based on an RGB color used as an illuminant. The spectrum is adjusted to account for the color space's
+ * standard illuminant.
+ */
+public final class RGBIlluminantSpectrum implements Spectrum {
+    private final double scale;
+    private final @NotNull SigmoidPolynomial polynomial;
+    private final @NotNull Spectrum illuminant;
+
+    public RGBIlluminantSpectrum(@NotNull ColorSpace cs, @NotNull ColorRGB rgb) {
+        if (rgb.r() < 0 || rgb.g() < 0 || rgb.b() < 0) {
+            throw new IllegalArgumentException();
+        }
+        var max = Math.max(rgb.r(), Math.max(rgb.g(), rgb.b()));
+        this.scale = 2 * max;
+        this.polynomial = cs.toPolynomial(scale != 0 ? rgb.div(scale) : ColorRGB.BLACK);
+        this.illuminant = cs.illuminant();
+    }
+
+    @Override
+    public double max() {
+        return scale * polynomial.max() * illuminant.max();
+    }
+
+    @Override
+    public double get(double lambda) {
+        return scale * polynomial.get(lambda) * illuminant.get(lambda);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/RGBUnboundedSpectrum.java

@@ -0,0 +1,30 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.SigmoidPolynomial;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import org.jetbrains.annotations.NotNull;
+
+public final class RGBUnboundedSpectrum implements Spectrum {
+    private final double scale;
+    private final @NotNull SigmoidPolynomial polynomial;
+
+    public RGBUnboundedSpectrum(@NotNull ColorSpace cs, @NotNull ColorRGB rgb) {
+        if (rgb.r() < 0 || rgb.g() < 0 || rgb.b() < 0) {
+            throw new IllegalArgumentException();
+        }
+        var max = Math.max(rgb.r(), Math.max(rgb.g(), rgb.b()));
+        this.scale = 2 * max;
+        this.polynomial = cs.toPolynomial(scale != 0 ? rgb.div(scale) : ColorRGB.BLACK);
+    }
+
+    @Override
+    public double max() {
+        return scale * polynomial.max();
+    }
+
+    @Override
+    public double get(double lambda) {
+        return scale * polynomial.get(lambda);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/SampledSpectrum.java

@@ -0,0 +1,139 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.math.IVec;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.ColorXYZ;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Arrays;
+
+// TODO use Vector API to parallelize operations
+public final class SampledSpectrum implements IVec<SampledSpectrum> {
+    private final double @NotNull[] values;
+
+    public SampledSpectrum(@NotNull SampledWavelengths lambdas, @NotNull Spectrum spectrum) {
+        var values = new double[lambdas.size()];
+        for (int i = 0; i < values.length; i++) {
+            values[i] = spectrum.get(lambdas.get(i));
+        }
+        this.values = values;
+    }
+
+    public SampledSpectrum(int count, double value) {
+        var values = new double[count];
+        Arrays.fill(values, value);
+        this.values = values;
+    }
+
+    private SampledSpectrum(double @NotNull[] values) {
+        this.values = values;
+    }
+
+    /*
+     * Math
+     */
+
+    public static @NotNull SampledSpectrum fma(@NotNull SampledSpectrum a, @NotNull SampledSpectrum b, @NotNull SampledSpectrum c) {
+        var out = new double[a.values.length];
+        for (int i = 0; i < a.values.length; i++) {
+            out[i] = Math.fma(a.values[i], b.values[i], c.values[i]);
+        }
+        return new SampledSpectrum(out);
+    }
+
+    public static @NotNull SampledSpectrum lerp(@NotNull SampledSpectrum a, @NotNull SampledSpectrum b, double t) {
+        if (t < 0) return a;
+        if (t > 1) return b;
+        var out = new double[a.values.length];
+        for (int i = 0; i < a.values.length; i++) {
+            out[i] = Math.fma(t, b.values[i] - a.values[i], a.values[i]);
+        }
+        return new SampledSpectrum(out);
+    }
+
+    @Override
+    public @NotNull SampledSpectrum plus(@NotNull SampledSpectrum other) {
+        var out = new double[other.values.length];
+        for (int i = 0; i < other.values.length; i++) {
+            out[i] = values[i] + other.values[i];
+        }
+        return new SampledSpectrum(out);
+    }
+
+    @Override
+    public @NotNull SampledSpectrum minus(@NotNull SampledSpectrum other) {
+        var out = new double[other.values.length];
+        for (int i = 0; i < other.values.length; i++) {
+            out[i] = values[i] - other.values[i];
+        }
+        return new SampledSpectrum(out);
+    }
+
+    @Override
+    public @NotNull SampledSpectrum times(@NotNull SampledSpectrum other) {
+        var out = new double[other.values.length];
+        for (int i = 0; i < other.values.length; i++) {
+            out[i] = values[i] * other.values[i];
+        }
+        return new SampledSpectrum(out);
+    }
+
+    @Override
+    public @NotNull SampledSpectrum times(double d) {
+        var out = new double[values.length];
+        for (int i = 0; i < values.length; i++) {
+            out[i] = values[i] * d;
+        }
+        return new SampledSpectrum(out);
+    }
+
+    @Override
+    public double @NotNull [] toArray() {
+        return Arrays.copyOf(values, values.length);
+    }
+
+    /*
+     * Accessors
+     */
+
+    @Override
+    public double get(int index) {
+        return values[index];
+    }
+
+    public int size() {
+        return values.length;
+    }
+
+    /*
+     * Object
+     */
+
+    @Override
+    public boolean equals(Object obj) {
+        return obj instanceof SampledSpectrum o && Arrays.equals(values, o.values);
+    }
+
+    @Override
+    public int hashCode() {
+        return Arrays.hashCode(values);
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return "SampledSpectrum[values=" + Arrays.toString(values) + "]";
+    }
+
+    /*
+     * Conversions
+     */
+
+    public @NotNull ColorXYZ toXYZ(@NotNull SampledWavelengths lambdas) {
+        return lambdas.toXYZ(this);
+    }
+
+    public @NotNull ColorRGB toRGB(@NotNull SampledWavelengths lambdas, @NotNull ColorSpace cs) {
+        return cs.toRGB(toXYZ(lambdas));
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/SampledWavelengths.java

@@ -0,0 +1,120 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.render.color.ColorXYZ;
+import org.jetbrains.annotations.Contract;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Arrays;
+
+/**
+ * A set of sampled wavelength that can be tracked together.
+ */
+public final class SampledWavelengths {
+    public static final SampledWavelengths EMPTY = new SampledWavelengths(new double[0], new double[0]);
+
+    private final double @NotNull[] lambdas;
+    private final double @NotNull[] pdf;
+
+    public static @NotNull SampledWavelengths uniform(double rng, int count) {
+        return uniform(rng, count, Spectrum.LAMBDA_MIN, Spectrum.LAMBDA_MAX);
+    }
+
+    public static @NotNull SampledWavelengths uniform(double rng, int count, double min, double max) {
+        var lambdas = new double[count];
+
+        // choose first sample at random
+        lambdas[0] = (1 - rng) * min + rng * max;
+
+        // choose next samples in equal intervals, wrapping if necessary
+        var delta = (max - min) / count;
+        for (int i = 1; i < count; i++) {
+            lambdas[i] = lambdas[i - 1] + delta;
+            if (lambdas[i] > max) {
+                lambdas[i] = min + (lambdas[i] - max);
+            }
+        }
+
+        var pdf = new double[count];
+        Arrays.fill(pdf, 1 / (max - min));
+        return new SampledWavelengths(lambdas, pdf);
+    }
+
+    private SampledWavelengths(double @NotNull[] lambdas, double @NotNull[] pdf) {
+        this.lambdas = lambdas;
+        this.pdf = pdf;
+    }
+
+    @Contract(pure = true)
+    public double get(int index) {
+        return lambdas[index];
+    }
+
+    @Contract(pure = true)
+    public int size() {
+        return lambdas.length;
+    }
+
+    /**
+     * Terminates the secondary wavelengths. This method should be called after a wavelength-dependent operation like
+     * refraction that makes it incorrect to track multiple wavelengths together.
+     */
+    @Contract(mutates = "this")
+    public double collapse() {
+        if (pdf.length >= 2 || pdf[1] != 0) {
+            Arrays.fill(pdf, 1, pdf.length, 0d);
+            pdf[0] /= pdf.length;
+        }
+        return lambdas[0];
+    }
+
+    /*
+     * Object
+     */
+
+    @Override
+    public boolean equals(Object obj) {
+        return obj instanceof SampledWavelengths o && Arrays.equals(lambdas, o.lambdas) && Arrays.equals(pdf, o.pdf);
+    }
+
+    @Override
+    public int hashCode() {
+        return 31 * Arrays.hashCode(lambdas) + Arrays.hashCode(pdf);
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return "SampledWavelengths[lambdas=" + Arrays.toString(lambdas) + ", pdf=" + Arrays.toString(pdf) + "]";
+    }
+
+    /*
+     * Conversions
+     */
+
+    @NotNull
+    ColorXYZ toXYZ(@NotNull SampledSpectrum spectrum) {
+        var x = Spectra.X.sample(this);
+        var y = Spectra.Y.sample(this);
+        var z = Spectra.Z.sample(this);
+
+        return new ColorXYZ(
+                toXYZ0(spectrum, x) / ColorXYZ.CIE_Y_INTEGRAL,
+                toXYZ0(spectrum, y) / ColorXYZ.CIE_Y_INTEGRAL,
+                toXYZ0(spectrum, z) / ColorXYZ.CIE_Y_INTEGRAL
+        );
+    }
+
+    private double toXYZ0(@NotNull SampledSpectrum spectrum, @NotNull SampledSpectrum cie) {
+        var avg = 0d;
+        for (int i = 0; i < spectrum.size(); i++) {
+            var pdf = this.pdf[i];
+            double value;
+            if (pdf == 0) {
+                value = 0;
+            } else {
+                value = spectrum.get(i) * cie.get(i) / pdf;
+            }
+            avg = Math.fma(1d / (i + 1), value - avg, avg);
+        }
+        return avg;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/ScaledSpectrum.java

@@ -0,0 +1,15 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import org.jetbrains.annotations.NotNull;
+
+public record ScaledSpectrum(@NotNull Spectrum spectrum, double scale) implements Spectrum {
+    @Override
+    public double max() {
+        return spectrum.max() * scale;
+    }
+
+    @Override
+    public double get(double lambda) {
+        return spectrum.get(lambda) * scale;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/Spectra.java

@@ -0,0 +1,410 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.render.color.ColorXYZ;
+import org.jetbrains.annotations.NotNull;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.InputStreamReader;
+import java.io.UncheckedIOException;
+import java.nio.charset.StandardCharsets;
+import java.util.ArrayList;
+
+public final class Spectra {
+    private static final String PATH_PREFIX = "/eu/jonahbauer/raytracing/spectrum/";
+
+    /**
+     * the CIE XYZ color matching curve for X
+     */
+    public static final Spectrum X = new DenselySampledSpectrum(new PiecewiseLinearSpectrum(CIE_XYZ.CIE_lambda, CIE_XYZ.CIE_X));
+    /**
+     * the CIE XYZ color matching curve for Y
+     */
+    public static final Spectrum Y = new DenselySampledSpectrum(new PiecewiseLinearSpectrum(CIE_XYZ.CIE_lambda, CIE_XYZ.CIE_Y));;
+    /**
+     * the CIE XYZ color matching curve for Z
+     */
+    public static final Spectrum Z = new DenselySampledSpectrum(new PiecewiseLinearSpectrum(CIE_XYZ.CIE_lambda, CIE_XYZ.CIE_Z));;
+    /**
+     * the CIE standard illuminant D50
+     * @see <a href="https://doi.org/10.25039/CIE.DS.hjfjmt59">CIE 2022, CIE standard illuminant D65, International Commission on Illumination (CIE), Vienna, Austria, DOI: 10.25039/CIE.DS.hjfjmt59</a>
+     */
+    public static final Spectrum D50 = read("CIE_std_illum_D50.csv", true);
+    /**
+     * the CIE standard illuminant D65
+     * @see <a href="https://doi.org/10.25039/CIE.DS.etgmuqt5">CIE 2022, Relative spectral power distributions of CIE standard illuminants A, D65 and D50 (wavelengths in standard air) (data table), International Commission on Illumination (CIE), Vienna, Austria, DOI:10.25039/CIE.DS.etgmuqt5</a>
+     */
+    public static final Spectrum D65 = read("CIE_std_illum_D65.csv", true);
+
+    public static final Spectrum BLACK = new ConstantSpectrum(0);
+    public static final Spectrum WHITE = new ConstantSpectrum(1);
+
+    private static @NotNull Spectrum read(@NotNull String path, boolean normalize) {
+        var lambda = new ArrayList<Double>();
+        var values = new ArrayList<Double>();
+
+        try (
+                var is = Spectra.class.getResourceAsStream(PATH_PREFIX + path);
+                var in = new BufferedReader(new InputStreamReader(is, StandardCharsets.US_ASCII))
+        ) {
+            String line;
+            while ((line = in.readLine()) != null) {
+                var parts = line.split(",");
+                lambda.add(Double.parseDouble(parts[0]));
+                values.add(Double.parseDouble(parts[1]));
+            }
+        } catch (IOException e) {
+            throw new UncheckedIOException(e);
+        }
+
+        var pls = new PiecewiseLinearSpectrum(
+                lambda.stream().mapToDouble(Double::doubleValue).toArray(),
+                values.stream().mapToDouble(Double::doubleValue).toArray()
+        );
+        if (normalize) {
+            return pls.scale(ColorXYZ.CIE_Y_INTEGRAL / Util.innerProduct(pls, Spectra.Y));
+        } else {
+            return pls;
+        }
+    }
+
+    private Spectra() {
+        throw new UnsupportedOperationException();
+    }
+
+    /**
+     * @see <a href="https://doi.org/10.25039/CIE.DS.xvudnb9b">CIE 2018, CIE 1931 colour-matching functions , 2 degree observer (data table),
+     * International Commission on Illumination (CIE), Vienna, Austria,
+     * DOI:10.25039/CIE.DS.xvudnb9b</a>
+     */
+    private static final class CIE_XYZ {
+        private static final double[] CIE_lambda = {
+                360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376,
+                377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393,
+                394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410,
+                411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427,
+                428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444,
+                445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461,
+                462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478,
+                479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495,
+                496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512,
+                513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529,
+                530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546,
+                547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563,
+                564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580,
+                581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597,
+                598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614,
+                615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631,
+                632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648,
+                649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665,
+                666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682,
+                683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699,
+                700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716,
+                717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733,
+                734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750,
+                751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767,
+                768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784,
+                785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801,
+                802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818,
+                819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830
+        };
+
+
+        private static final double[] CIE_X = {
+                // CIE X function values
+                0.0001299000,   0.0001458470,   0.0001638021,   0.0001840037,   0.0002066902,
+                0.0002321000,   0.0002607280,   0.0002930750,   0.0003293880,   0.0003699140,
+                0.0004149000,   0.0004641587,   0.0005189860,   0.0005818540,   0.0006552347,
+                0.0007416000,   0.0008450296,   0.0009645268,   0.001094949,    0.001231154,
+                0.001368000,    0.001502050,    0.001642328,    0.001802382,    0.001995757,
+                0.002236000,    0.002535385,    0.002892603,    0.003300829,    0.003753236,
+                0.004243000,    0.004762389,    0.005330048,    0.005978712,    0.006741117,
+                0.007650000,    0.008751373,    0.01002888,     0.01142170,     0.01286901,
+                0.01431000,     0.01570443,     0.01714744,     0.01878122,     0.02074801,
+                0.02319000,     0.02620736,     0.02978248,     0.03388092,     0.03846824,
+                0.04351000,     0.04899560,     0.05502260,     0.06171880,     0.06921200,
+                0.07763000,     0.08695811,     0.09717672,     0.1084063,      0.1207672,
+                0.1343800,      0.1493582,      0.1653957,      0.1819831,      0.1986110,
+                0.2147700,      0.2301868,      0.2448797,      0.2587773,      0.2718079,
+                0.2839000,      0.2949438,      0.3048965,      0.3137873,      0.3216454,
+                0.3285000,      0.3343513,      0.3392101,      0.3431213,      0.3461296,
+                0.3482800,      0.3495999,      0.3501474,      0.3500130,      0.3492870,
+                0.3480600,      0.3463733,      0.3442624,      0.3418088,      0.3390941,
+                0.3362000,      0.3331977,      0.3300411,      0.3266357,      0.3228868,
+                0.3187000,      0.3140251,      0.3088840,      0.3032904,      0.2972579,
+                0.2908000,      0.2839701,      0.2767214,      0.2689178,      0.2604227,
+                0.2511000,      0.2408475,      0.2298512,      0.2184072,      0.2068115,
+                0.1953600,      0.1842136,      0.1733273,      0.1626881,      0.1522833,
+                0.1421000,      0.1321786,      0.1225696,      0.1132752,      0.1042979,
+                0.09564000,     0.08729955,     0.07930804,     0.07171776,     0.06458099,
+                0.05795001,     0.05186211,     0.04628152,     0.04115088,     0.03641283,
+                0.03201000,     0.02791720,     0.02414440,     0.02068700,     0.01754040,
+                0.01470000,     0.01216179,     0.009919960,    0.007967240,    0.006296346,
+                0.004900000,    0.003777173,    0.002945320,    0.002424880,    0.002236293,
+                0.002400000,    0.002925520,    0.003836560,    0.005174840,    0.006982080,
+                0.009300000,    0.01214949,     0.01553588,     0.01947752,     0.02399277,
+                0.02910000,     0.03481485,     0.04112016,     0.04798504,     0.05537861,
+                0.06327000,     0.07163501,     0.08046224,     0.08973996,     0.09945645,
+                0.1096000,      0.1201674,      0.1311145,      0.1423679,      0.1538542,
+                0.1655000,      0.1772571,      0.1891400,      0.2011694,      0.2133658,
+                0.2257499,      0.2383209,      0.2510668,      0.2639922,      0.2771017,
+                0.2904000,      0.3038912,      0.3175726,      0.3314384,      0.3454828,
+                0.3597000,      0.3740839,      0.3886396,      0.4033784,      0.4183115,
+                0.4334499,      0.4487953,      0.4643360,      0.4800640,      0.4959713,
+                0.5120501,      0.5282959,      0.5446916,      0.5612094,      0.5778215,
+                0.5945000,      0.6112209,      0.6279758,      0.6447602,      0.6615697,
+                0.6784000,      0.6952392,      0.7120586,      0.7288284,      0.7455188,
+                0.7621000,      0.7785432,      0.7948256,      0.8109264,      0.8268248,
+                0.8425000,      0.8579325,      0.8730816,      0.8878944,      0.9023181,
+                0.9163000,      0.9297995,      0.9427984,      0.9552776,      0.9672179,
+                0.9786000,      0.9893856,      0.9995488,      1.0090892,      1.0180064,
+                1.0263000,      1.0339827,      1.0409860,      1.0471880,      1.0524667,
+                1.0567000,      1.0597944,      1.0617992,      1.0628068,      1.0629096,
+                1.0622000,      1.0607352,      1.0584436,      1.0552244,      1.0509768,
+                1.0456000,      1.0390369,      1.0313608,      1.0226662,      1.0130477,
+                1.0026000,      0.9913675,      0.9793314,      0.9664916,      0.9528479,
+                0.9384000,      0.9231940,      0.9072440,      0.8905020,      0.8729200,
+                0.8544499,      0.8350840,      0.8149460,      0.7941860,      0.7729540,
+                0.7514000,      0.7295836,      0.7075888,      0.6856022,      0.6638104,
+                0.6424000,      0.6215149,      0.6011138,      0.5811052,      0.5613977,
+                0.5419000,      0.5225995,      0.5035464,      0.4847436,      0.4661939,
+                0.4479000,      0.4298613,      0.4120980,      0.3946440,      0.3775333,
+                0.3608000,      0.3444563,      0.3285168,      0.3130192,      0.2980011,
+                0.2835000,      0.2695448,      0.2561184,      0.2431896,      0.2307272,
+                0.2187000,      0.2070971,      0.1959232,      0.1851708,      0.1748323,
+                0.1649000,      0.1553667,      0.1462300,      0.1374900,      0.1291467,
+                0.1212000,      0.1136397,      0.1064650,      0.09969044,     0.09333061,
+                0.08740000,     0.08190096,     0.07680428,     0.07207712,     0.06768664,
+                0.06360000,     0.05980685,     0.05628216,     0.05297104,     0.04981861,
+                0.04677000,     0.04378405,     0.04087536,     0.03807264,     0.03540461,
+                0.03290000,     0.03056419,     0.02838056,     0.02634484,     0.02445275,
+                0.02270000,     0.02108429,     0.01959988,     0.01823732,     0.01698717,
+                0.01584000,     0.01479064,     0.01383132,     0.01294868,     0.01212920,
+                0.01135916,     0.01062935,     0.009938846,    0.009288422,    0.008678854,
+                0.008110916,    0.007582388,    0.007088746,    0.006627313,    0.006195408,
+                0.005790346,    0.005409826,    0.005052583,    0.004717512,    0.004403507,
+                0.004109457,    0.003833913,    0.003575748,    0.003334342,    0.003109075,
+                0.002899327,    0.002704348,    0.002523020,    0.002354168,    0.002196616,
+                0.002049190,    0.001910960,    0.001781438,    0.001660110,    0.001546459,
+                0.001439971,    0.001340042,    0.001246275,    0.001158471,    0.001076430,
+                0.0009999493,   0.0009287358,   0.0008624332,   0.0008007503,   0.0007433960,
+                0.0006900786,   0.0006405156,   0.0005945021,   0.0005518646,   0.0005124290,
+                0.0004760213,   0.0004424536,   0.0004115117,   0.0003829814,   0.0003566491,
+                0.0003323011,   0.0003097586,   0.0002888871,   0.0002695394,   0.0002515682,
+                0.0002348261,   0.0002191710,   0.0002045258,   0.0001908405,   0.0001780654,
+                0.0001661505,   0.0001550236,   0.0001446219,   0.0001349098,   0.0001258520,
+                0.0001174130,   0.0001095515,   0.0001022245,   0.00009539445,  0.00008902390,
+                0.00008307527,  0.00007751269,  0.00007231304,  0.00006745778,  0.00006292844,
+                0.00005870652,  0.00005477028,  0.00005109918,  0.00004767654,  0.00004448567,
+                0.00004150994,  0.00003873324,  0.00003614203,  0.00003372352,  0.00003146487,
+                0.00002935326,  0.00002737573,  0.00002552433,  0.00002379376,  0.00002217870,
+                0.00002067383,  0.00001927226,  0.00001796640,  0.00001674991,  0.00001561648,
+                0.00001455977,  0.00001357387,  0.00001265436,  0.00001179723,  0.00001099844,
+                0.00001025398,  0.000009559646, 0.000008912044, 0.000008308358, 0.000007745769,
+                0.000007221456, 0.000006732475, 0.000006276423, 0.000005851304, 0.000005455118,
+                0.000005085868, 0.000004741466, 0.000004420236, 0.000004120783, 0.000003841716,
+                0.000003581652, 0.000003339127, 0.000003112949, 0.000002902121, 0.000002705645,
+                0.000002522525, 0.000002351726, 0.000002192415, 0.000002043902, 0.000001905497,
+                0.000001776509, 0.000001656215, 0.000001544022, 0.000001439440, 0.000001341977,
+                0.000001251141
+        };
+
+        private static final double[] CIE_Y = {
+                // CIE Y function values
+                0.000003917000,  0.000004393581,  0.000004929604,  0.000005532136,  0.000006208245,
+                0.000006965000,  0.000007813219,  0.000008767336,  0.000009839844,  0.00001104323,
+                0.00001239000,   0.00001388641,   0.00001555728,   0.00001744296,   0.00001958375,
+                0.00002202000,   0.00002483965,   0.00002804126,   0.00003153104,   0.00003521521,
+                0.00003900000,   0.00004282640,   0.00004691460,   0.00005158960,   0.00005717640,
+                0.00006400000,   0.00007234421,   0.00008221224,   0.00009350816,   0.0001061361,
+                0.0001200000,    0.0001349840,    0.0001514920,    0.0001702080,    0.0001918160,
+                0.0002170000,    0.0002469067,    0.0002812400,    0.0003185200,    0.0003572667,
+                0.0003960000,    0.0004337147,    0.0004730240,    0.0005178760,    0.0005722187,
+                0.0006400000,    0.0007245600,    0.0008255000,    0.0009411600,    0.001069880,
+                0.001210000,     0.001362091,     0.001530752,     0.001720368,     0.001935323,
+                0.002180000,     0.002454800,     0.002764000,     0.003117800,     0.003526400,
+                0.004000000,     0.004546240,     0.005159320,     0.005829280,     0.006546160,
+                0.007300000,     0.008086507,     0.008908720,     0.009767680,     0.01066443,
+                0.01160000,      0.01257317,      0.01358272,      0.01462968,      0.01571509,
+                0.01684000,      0.01800736,      0.01921448,      0.02045392,      0.02171824,
+                0.02300000,      0.02429461,      0.02561024,      0.02695857,      0.02835125,
+                0.02980000,      0.03131083,      0.03288368,      0.03452112,      0.03622571,
+                0.03800000,      0.03984667,      0.04176800,      0.04376600,      0.04584267,
+                0.04800000,      0.05024368,      0.05257304,      0.05498056,      0.05745872,
+                0.06000000,      0.06260197,      0.06527752,      0.06804208,      0.07091109,
+                0.07390000,      0.07701600,      0.08026640,      0.08366680,      0.08723280,
+                0.09098000,      0.09491755,      0.09904584,      0.1033674,       0.1078846,
+                0.1126000,       0.1175320,       0.1226744,       0.1279928,       0.1334528,
+                0.1390200,       0.1446764,       0.1504693,       0.1564619,       0.1627177,
+                0.1693000,       0.1762431,       0.1835581,       0.1912735,       0.1994180,
+                0.2080200,       0.2171199,       0.2267345,       0.2368571,       0.2474812,
+                0.2586000,       0.2701849,       0.2822939,       0.2950505,       0.3085780,
+                0.3230000,       0.3384021,       0.3546858,       0.3716986,       0.3892875,
+                0.4073000,       0.4256299,       0.4443096,       0.4633944,       0.4829395,
+                0.5030000,       0.5235693,       0.5445120,       0.5656900,       0.5869653,
+                0.6082000,       0.6293456,       0.6503068,       0.6708752,       0.6908424,
+                0.7100000,       0.7281852,       0.7454636,       0.7619694,       0.7778368,
+                0.7932000,       0.8081104,       0.8224962,       0.8363068,       0.8494916,
+                0.8620000,       0.8738108,       0.8849624,       0.8954936,       0.9054432,
+                0.9148501,       0.9237348,       0.9320924,       0.9399226,       0.9472252,
+                0.9540000,       0.9602561,       0.9660074,       0.9712606,       0.9760225,
+                0.9803000,       0.9840924,       0.9874812,       0.9903128,       0.9928116,
+                0.9949501,       0.9967108,       0.9980983,       0.9991120,       0.9997482,
+                1.0000000,       0.9998567,       0.9993046,       0.9983255,       0.9968987,
+                0.9950000,       0.9926005,       0.9897426,       0.9864444,       0.9827241,
+                0.9786000,       0.9740837,       0.9691712,       0.9638568,       0.9581349,
+                0.9520000,       0.9454504,       0.9384992,       0.9311628,       0.9234576,
+                0.9154000,       0.9070064,       0.8982772,       0.8892048,       0.8797816,
+                0.8700000,       0.8598613,       0.8493920,       0.8386220,       0.8275813,
+                0.8163000,       0.8047947,       0.7930820,       0.7811920,       0.7691547,
+                0.7570000,       0.7447541,       0.7324224,       0.7200036,       0.7074965,
+                0.6949000,       0.6822192,       0.6694716,       0.6566744,       0.6438448,
+                0.6310000,       0.6181555,       0.6053144,       0.5924756,       0.5796379,
+                0.5668000,       0.5539611,       0.5411372,       0.5283528,       0.5156323,
+                0.5030000,       0.4904688,       0.4780304,       0.4656776,       0.4534032,
+                0.4412000,       0.4290800,       0.4170360,       0.4050320,       0.3930320,
+                0.3810000,       0.3689184,       0.3568272,       0.3447768,       0.3328176,
+                0.3210000,       0.3093381,       0.2978504,       0.2865936,       0.2756245,
+                0.2650000,       0.2547632,       0.2448896,       0.2353344,       0.2260528,
+                0.2170000,       0.2081616,       0.1995488,       0.1911552,       0.1829744,
+                0.1750000,       0.1672235,       0.1596464,       0.1522776,       0.1451259,
+                0.1382000,       0.1315003,       0.1250248,       0.1187792,       0.1127691,
+                0.1070000,       0.1014762,       0.09618864,      0.09112296,      0.08626485,
+                0.08160000,      0.07712064,      0.07282552,      0.06871008,      0.06476976,
+                0.06100000,      0.05739621,      0.05395504,      0.05067376,      0.04754965,
+                0.04458000,      0.04175872,      0.03908496,      0.03656384,      0.03420048,
+                0.03200000,      0.02996261,      0.02807664,      0.02632936,      0.02470805,
+                0.02320000,      0.02180077,      0.02050112,      0.01928108,      0.01812069,
+                0.01700000,      0.01590379,      0.01483718,      0.01381068,      0.01283478,
+                0.01192000,      0.01106831,      0.01027339,      0.009533311,     0.008846157,
+                0.008210000,     0.007623781,     0.007085424,     0.006591476,     0.006138485,
+                0.005723000,     0.005343059,     0.004995796,     0.004676404,     0.004380075,
+                0.004102000,     0.003838453,     0.003589099,     0.003354219,     0.003134093,
+                0.002929000,     0.002738139,     0.002559876,     0.002393244,     0.002237275,
+                0.002091000,     0.001953587,     0.001824580,     0.001703580,     0.001590187,
+                0.001484000,     0.001384496,     0.001291268,     0.001204092,     0.001122744,
+                0.001047000,     0.0009765896,    0.0009111088,    0.0008501332,    0.0007932384,
+                0.0007400000,    0.0006900827,    0.0006433100,    0.0005994960,    0.0005584547,
+                0.0005200000,    0.0004839136,    0.0004500528,    0.0004183452,    0.0003887184,
+                0.0003611000,    0.0003353835,    0.0003114404,    0.0002891656,    0.0002684539,
+                0.0002492000,    0.0002313019,    0.0002146856,    0.0001992884,    0.0001850475,
+                0.0001719000,    0.0001597781,    0.0001486044,    0.0001383016,    0.0001287925,
+                0.0001200000,    0.0001118595,    0.0001043224,    0.00009733560,   0.00009084587,
+                0.00008480000,   0.00007914667,   0.00007385800,   0.00006891600,   0.00006430267,
+                0.00006000000,   0.00005598187,   0.00005222560,   0.00004871840,   0.00004544747,
+                0.00004240000,   0.00003956104,   0.00003691512,   0.00003444868,   0.00003214816,
+                0.00003000000,   0.00002799125,   0.00002611356,   0.00002436024,   0.00002272461,
+                0.00002120000,   0.00001977855,   0.00001845285,   0.00001721687,   0.00001606459,
+                0.00001499000,   0.00001398728,   0.00001305155,   0.00001217818,   0.00001136254,
+                0.00001060000,   0.000009885877,  0.000009217304,  0.000008592362,  0.000008009133,
+                0.000007465700,  0.000006959567,  0.000006487995,  0.000006048699,  0.000005639396,
+                0.000005257800,  0.000004901771,  0.000004569720,  0.000004260194,  0.000003971739,
+                0.000003702900,  0.000003452163,  0.000003218302,  0.000003000300,  0.000002797139,
+                0.000002607800,  0.000002431220,  0.000002266531,  0.000002113013,  0.000001969943,
+                0.000001836600,  0.000001712230,  0.000001596228,  0.000001488090,  0.000001387314,
+                0.000001293400,  0.000001205820,  0.000001124143,  0.000001048009,  0.0000009770578,
+                0.0000009109300, 0.0000008492513, 0.0000007917212, 0.0000007380904, 0.0000006881098,
+                0.0000006415300, 0.0000005980895, 0.0000005575746, 0.0000005198080, 0.0000004846123,
+                0.0000004518100
+        };
+
+        private static final double[] CIE_Z = {
+                // CIE Z function values
+                0.0006061000,   0.0006808792,   0.0007651456,   0.0008600124,   0.0009665928,
+                0.001086000,    0.001220586,    0.001372729,    0.001543579,    0.001734286,
+                0.001946000,    0.002177777,    0.002435809,    0.002731953,    0.003078064,
+                0.003486000,    0.003975227,    0.004540880,    0.005158320,    0.005802907,
+                0.006450001,    0.007083216,    0.007745488,    0.008501152,    0.009414544,
+                0.01054999,     0.01196580,     0.01365587,     0.01558805,     0.01773015,
+                0.02005001,     0.02251136,     0.02520288,     0.02827972,     0.03189704,
+                0.03621000,     0.04143771,     0.04750372,     0.05411988,     0.06099803,
+                0.06785001,     0.07448632,     0.08136156,     0.08915364,     0.09854048,
+                0.1102000,      0.1246133,      0.1417017,      0.1613035,      0.1832568,
+                0.2074000,      0.2336921,      0.2626114,      0.2947746,      0.3307985,
+                0.3713000,      0.4162091,      0.4654642,      0.5196948,      0.5795303,
+                0.6456000,      0.7184838,      0.7967133,      0.8778459,      0.9594390,
+                1.0390501,      1.1153673,      1.1884971,      1.2581233,      1.3239296,
+                1.3856000,      1.4426352,      1.4948035,      1.5421903,      1.5848807,
+                1.6229600,      1.6564048,      1.6852959,      1.7098745,      1.7303821,
+                1.7470600,      1.7600446,      1.7696233,      1.7762637,      1.7804334,
+                1.7826000,      1.7829682,      1.7816998,      1.7791982,      1.7758671,
+                1.7721100,      1.7682589,      1.7640390,      1.7589438,      1.7524663,
+                1.7441000,      1.7335595,      1.7208581,      1.7059369,      1.6887372,
+                1.6692000,      1.6475287,      1.6234127,      1.5960223,      1.5645280,
+                1.5281000,      1.4861114,      1.4395215,      1.3898799,      1.3387362,
+                1.2876400,      1.2374223,      1.1878243,      1.1387611,      1.0901480,
+                1.0419000,      0.9941976,      0.9473473,      0.9014531,      0.8566193,
+                0.8129501,      0.7705173,      0.7294448,      0.6899136,      0.6521049,
+                0.6162000,      0.5823286,      0.5504162,      0.5203376,      0.4919673,
+                0.4651800,      0.4399246,      0.4161836,      0.3938822,      0.3729459,
+                0.3533000,      0.3348578,      0.3175521,      0.3013375,      0.2861686,
+                0.2720000,      0.2588171,      0.2464838,      0.2347718,      0.2234533,
+                0.2123000,      0.2011692,      0.1901196,      0.1792254,      0.1685608,
+                0.1582000,      0.1481383,      0.1383758,      0.1289942,      0.1200751,
+                0.1117000,      0.1039048,      0.09666748,     0.08998272,     0.08384531,
+                0.07824999,     0.07320899,     0.06867816,     0.06456784,     0.06078835,
+                0.05725001,     0.05390435,     0.05074664,     0.04775276,     0.04489859,
+                0.04216000,     0.03950728,     0.03693564,     0.03445836,     0.03208872,
+                0.02984000,     0.02771181,     0.02569444,     0.02378716,     0.02198925,
+                0.02030000,     0.01871805,     0.01724036,     0.01586364,     0.01458461,
+                0.01340000,     0.01230723,     0.01130188,     0.01037792,     0.009529306,
+                0.008749999,    0.008035200,    0.007381600,    0.006785400,    0.006242800,
+                0.005749999,    0.005303600,    0.004899800,    0.004534200,    0.004202400,
+                0.003900000,    0.003623200,    0.003370600,    0.003141400,    0.002934800,
+                0.002749999,    0.002585200,    0.002438600,    0.002309400,    0.002196800,
+                0.002100000,    0.002017733,    0.001948200,    0.001889800,    0.001840933,
+                0.001800000,    0.001766267,    0.001737800,    0.001711200,    0.001683067,
+                0.001650001,    0.001610133,    0.001564400,    0.001513600,    0.001458533,
+                0.001400000,    0.001336667,    0.001270000,    0.001205000,    0.001146667,
+                0.001100000,    0.001068800,    0.001049400,    0.001035600,    0.001021200,
+                0.001000000,    0.0009686400,   0.0009299200,   0.0008868800,   0.0008425600,
+                0.0008000000,   0.0007609600,   0.0007236800,   0.0006859200,   0.0006454400,
+                0.0006000000,   0.0005478667,   0.0004916000,   0.0004354000,   0.0003834667,
+                0.0003400000,   0.0003072533,   0.0002831600,   0.0002654400,   0.0002518133,
+                0.0002400000,   0.0002295467,   0.0002206400,   0.0002119600,   0.0002021867,
+                0.0001900000,   0.0001742133,   0.0001556400,   0.0001359600,   0.0001168533,
+                0.0001000000,   0.00008613333,  0.00007460000,  0.00006500000,  0.00005693333,
+                0.00004999999,  0.00004416000,  0.00003948000,  0.00003572000,  0.00003264000,
+                0.00003000000,  0.00002765333,  0.00002556000,  0.00002364000,  0.00002181333,
+                0.00002000000,  0.00001813333,  0.00001620000,  0.00001420000,  0.00001213333,
+                0.00001000000,  0.000007733333, 0.000005400000, 0.000003200000, 0.000001333333,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000,
+                0.000000000000
+        };
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/Spectrum.java

@@ -0,0 +1,79 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.color.ColorXYZ;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.texture.Texture;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.SkyBox;
+import org.jetbrains.annotations.NotNull;
+
+public interface Spectrum extends Texture, SkyBox {
+    int LAMBDA_MIN = 360;
+    int LAMBDA_MAX = 830;
+
+    /**
+     * {@return the maximum value of <code>this</code> spectrum over the range of wavelengths}
+     */
+    double max();
+
+    /**
+     * {@return the value of <code>this</code> spectrum at a given wavelength}
+     * @param lambda the wavelength in nanometers
+     */
+    double get(double lambda);
+
+    default @NotNull Spectrum scale(double scale) {
+        return new ScaledSpectrum(this, scale);
+    }
+
+    default @NotNull SampledSpectrum sample(@NotNull SampledWavelengths lambdas) {
+        return new SampledSpectrum(lambdas, this);
+    }
+
+    default @NotNull ColorXYZ toXYZ() {
+        return new ColorXYZ(
+                Util.innerProduct(Spectra.X, this) / ColorXYZ.CIE_Y_INTEGRAL,
+                Util.innerProduct(Spectra.Y, this) / ColorXYZ.CIE_Y_INTEGRAL,
+                Util.innerProduct(Spectra.Z, this) / ColorXYZ.CIE_Y_INTEGRAL
+        );
+    }
+
+    default @NotNull ColorRGB toRGB(@NotNull ColorSpace cs) {
+        return cs.toRGB(toXYZ());
+    }
+
+    /*
+     * Texture
+     */
+
+    @Override
+    default @NotNull Spectrum get(@NotNull HitResult hit) {
+        return this;
+    }
+
+    @Override
+    default @NotNull Spectrum get(double u, double v, @NotNull Vec3 p) {
+        return this;
+    }
+
+    @Override
+    default boolean isUVRequired() {
+        return false;
+    }
+
+    /*
+     * SkyBox
+     */
+
+    @Override
+    default @NotNull SampledSpectrum getColor(@NotNull Ray ray) {
+        return this.sample(ray.lambda());
+    }
+
+    enum Type {
+        ALBEDO, ILLUMINANT, UNBOUNDED
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/spectrum/Util.java

@@ -0,0 +1,17 @@
+package eu.jonahbauer.raytracing.render.spectrum;
+
+import org.jetbrains.annotations.NotNull;
+
+final class Util {
+    private Util() {
+        throw new UnsupportedOperationException();
+    }
+
+    public static double innerProduct(@NotNull Spectrum f, @NotNull Spectrum g) {
+        var integral = 0.0;
+        for (var lambda = Spectrum.LAMBDA_MIN; lambda <= Spectrum.LAMBDA_MAX; lambda++) {
+            integral += f.get(lambda) * g.get(lambda);
+        }
+        return integral;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/texture/CheckerTexture.java

@@ -0,0 +1,22 @@
+package eu.jonahbauer.raytracing.render.texture;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import org.jetbrains.annotations.NotNull;
+
+public record CheckerTexture(double scale, @NotNull Texture even, @NotNull Texture odd) implements Texture {
+
+    @Override
+    public @NotNull Spectrum get(double u, double v, @NotNull Vec3 p) {
+        var x = (int) Math.floor(p.x() / scale);
+        var y = (int) Math.floor(p.y() / scale);
+        var z = (int) Math.floor(p.z() / scale);
+        var even = (x + y + z) % 2 == 0;
+        return even ? even().get(u, v, p) : odd().get(u, v, p);
+    }
+
+    @Override
+    public boolean isUVRequired() {
+        return even.isUVRequired() || odd.isUVRequired();
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/texture/ImageTexture.java

@@ -0,0 +1,55 @@
+package eu.jonahbauer.raytracing.render.texture;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.color.ColorRGB;
+import eu.jonahbauer.raytracing.render.color.ColorSpace;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import org.jetbrains.annotations.NotNull;
+
+import javax.imageio.ImageIO;
+import java.awt.image.BufferedImage;
+import java.io.IOException;
+import java.io.UncheckedIOException;
+import java.util.Objects;
+
+public final class ImageTexture implements Texture {
+    private static final String PATH_PREFIX = "/eu/jonahbauer/raytracing/textures/";
+
+    private final int width;
+    private final int height;
+    private final @NotNull Spectrum[][] spectra;
+
+    public ImageTexture(@NotNull BufferedImage image, @NotNull ColorSpace cs, @NotNull Spectrum.Type type) {
+        this.width = image.getWidth();
+        this.height = image.getHeight();
+        this.spectra = new Spectrum[height][width];
+
+        for (int y = 0; y < height; y++) {
+            for (int x = 0; x < width; x++) {
+                var rgb = cs.decode(new ColorRGB(image.getRGB(x, y)));
+                spectra[y][x] = cs.toSpectrum(rgb, type);
+            }
+        }
+    }
+
+    public ImageTexture(@NotNull String path, @NotNull ColorSpace cs) {
+        this(read(path), cs, Spectrum.Type.ALBEDO);
+    }
+
+    private static @NotNull BufferedImage read(@NotNull String path) {
+        try (var in = Objects.requireNonNull(ImageTexture.class.getResourceAsStream(PATH_PREFIX + path))) {
+            return ImageIO.read(in);
+        } catch (IOException ex) {
+            throw new UncheckedIOException(ex);
+        }
+    }
+
+    @Override
+    public @NotNull Spectrum get(double u, double v, @NotNull Vec3 p) {
+        u = Math.clamp(u, 0, 1);
+        v = 1 - Math.clamp(v, 0, 1);
+        int x = (int) (u * (width - 1));
+        int y = (int) (v * (height - 1));
+        return spectra[y][x];
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/texture/PerlinTexture.java

@@ -0,0 +1,154 @@
+package eu.jonahbauer.raytracing.render.texture;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Random;
+import java.util.function.DoubleFunction;
+import java.util.random.RandomGenerator;
+
+public final class PerlinTexture implements Texture {
+    private static final int POINT_COUNT = 256;
+    private static final @NotNull Random RANDOM = new Random();
+    private static final @NotNull DoubleFunction<Spectrum> GREYSCALE = Spectra.WHITE::scale;
+
+    private final double scale;
+    private final int turbulence;
+    private final @NotNull DoubleFunction<Spectrum> color;
+
+    private final int mask;
+    private final Vec3[] randvec;
+    private final int[] permX;
+    private final int[] permY;
+    private final int[] permZ;
+
+    public PerlinTexture() {
+        this(1.0);
+    }
+
+    public PerlinTexture(double scale) {
+        this(scale, 7);
+    }
+
+    public PerlinTexture(double scale, int turbulence) {
+        this(scale, turbulence, GREYSCALE);
+    }
+
+    public PerlinTexture(double scale, int turbulence, @NotNull DoubleFunction<Spectrum> color) {
+        this(scale, turbulence, color, POINT_COUNT, RANDOM);
+    }
+
+    public PerlinTexture(
+            double scale, int turbulence, @NotNull DoubleFunction<Spectrum> color,
+            int count, @NotNull RandomGenerator random
+    ) {
+        if ((count & (count - 1)) != 0) throw new IllegalArgumentException("count must be a power of two");
+        if (turbulence <= 0) throw new IllegalArgumentException("turbulence must be positive");
+
+        this.scale = scale;
+        this.turbulence = turbulence;
+        this.color = Objects.requireNonNull(color, "color");
+
+        this.mask = count - 1;
+        this.randvec = new Vec3[count];
+        for (int i = 0; i < count; i++) {
+            this.randvec[i] = Vec3.random(random);
+        }
+        this.permX = generatePerm(count, random);
+        this.permY = generatePerm(count, random);
+        this.permZ = generatePerm(count, random);
+    }
+
+    private static int @NotNull[] generatePerm(int count, @NotNull RandomGenerator random) {
+        int[] p = new int[count];
+        for (int i = 0; i < count; i++) {
+            p[i] = i;
+        }
+        permutate(p, random);
+        return p;
+    }
+
+    private static void permutate(int @NotNull[] p, @NotNull RandomGenerator random) {
+        for (int i = p.length - 1; i > 0; i--) {
+            int target = random.nextInt(i);
+            int tmp = p[i];
+            p[i] = p[target];
+            p[target] = tmp;
+        }
+    }
+
+    public double getNoise(@NotNull Vec3 p) {
+        var x = p.x() * scale;
+        var y = p.y() * scale;
+        var z = p.z() * scale;
+
+        var u = x - Math.floor(x);
+        var v = y - Math.floor(y);
+        var w = z - Math.floor(z);
+
+        int i = (int) Math.floor(x);
+        int j = (int) Math.floor(y);
+        int k = (int) Math.floor(z);
+
+        var c = new Vec3[8];
+        for (int di = 0; di < 2; di++) {
+            for (int dj = 0; dj < 2; dj++) {
+                for (int dk = 0; dk < 2; dk++) {
+                    c[di << 2 | dj << 1 | dk] = randvec[permX[(i + di) & mask] ^ permY[(j + dj) & mask] ^ permZ[(k + dk) & mask]];
+                }
+            }
+        }
+
+        return interpolate(c, u, v, w);
+    }
+
+    public double getNoise(@NotNull Vec3 p, int depth) {
+        var accum = 0.0;
+        var temp = p;
+        var weight = 1.0;
+
+        for (int i = 0; i < depth; i++) {
+            accum = Math.fma(weight, getNoise(temp), accum);
+            weight *= 0.5;
+            temp = temp.times(2);
+        }
+
+        return accum;
+    }
+
+    @Override
+    public @NotNull Spectrum get(double u, double v, @NotNull Vec3 p) {
+        var noise = getNoise(p, turbulence);
+        var t = Math.fma(0.5, Math.sin(Math.PI * noise), 0.5);
+        return color.apply(t);
+    }
+
+    private static double interpolate(Vec3[] c, double u, double v, double w) {
+        var uu = u * u * Math.fma(-2, u, 3);
+        var vv = v * v * Math.fma(-2, v, 3);
+        var ww = w * w * Math.fma(-2, w, 3);
+
+        var accum = 0.0;
+        for (int i = 0; i < 2; i++) {
+            for (int j = 0; j < 2; j++) {
+                for (int k = 0; k < 2; k++) {
+                    var vec = c[i << 2 | j << 1 | k];
+                    var dot = (u - i) * vec.x() + (v - j) * vec.y() + (w - k) * vec.z();
+                    accum += Math.fma(i, uu, (1 - i) * (1 - uu))
+                             * Math.fma(j, vv, (1 - j) * (1 - vv))
+                             * Math.fma(k, ww, (1 - k) * (1 - ww))
+                             * dot;
+                }
+            }
+        }
+        return accum;
+    }
+
+    @Override
+    public boolean isUVRequired() {
+        return false;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/render/texture/Texture.java

@@ -0,0 +1,33 @@
+package eu.jonahbauer.raytracing.render.texture;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import org.jetbrains.annotations.NotNull;
+
+public interface Texture {
+    /**
+     * {@return the color of <code>this</code> texture for a hit}
+     */
+    default @NotNull Spectrum get(@NotNull HitResult hit) {
+        return get(hit.u(), hit.v(), hit.position());
+    }
+
+    /**
+     * {@return the color of <code>this</code> texture at the specified position}
+     * @param u the texture u coordinate
+     * @param v the texture v coordinate
+     * @param p the position
+     */
+    @NotNull Spectrum get(double u, double v, @NotNull Vec3 p);
+
+    /**
+     * Returns whether {@link #get(double, double, Vec3)} uses the {@code u} and/or {@code v} parameters.
+     * When a texture indicates that the {@code u} and {@code v} coordinates are not required, the calculation may be
+     * skipped and {@link Double#NaN} will be passed.
+     * @return whether {@link #get(double, double, Vec3)} uses the {@code u} and/or {@code v} parameters
+     */
+    default boolean isUVRequired() {
+        return true;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/HitResult.java

@@ -1,22 +1,32 @@
 package eu.jonahbauer.raytracing.scene;
 
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
 import eu.jonahbauer.raytracing.math.Vec3;
-import eu.jonahbauer.raytracing.material.Material;
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.render.texture.Texture;
 import org.jetbrains.annotations.NotNull;
 
-import java.util.Objects;
+import java.util.random.RandomGenerator;
 
+/**
+ * The result of a {@linkplain Hittable#hit(Ray, Range, RandomGenerator) hit}.
+ * @param t the {@code t} value at which the hit occurs
+ * @param position the position of the hit
+ * @param normal the surface normal at the hit position
+ * @param target the hit target (for debug purposes only)
+ * @param material the material of the surface
+ * @param u the texture u coordinate (or {@code Double.NaN} if the {@linkplain Material#texture() material's texture} does {@linkplain Texture#isUVRequired() not depend} on the uv-coordinates)
+ * @param v the texture v coordinate (or {@code Double.NaN} if the {@linkplain Material#texture() material's texture} does {@linkplain Texture#isUVRequired() not depend} on the uv-coordinates)
+ * @param isFrontFace whether the front or the back of the surface was it
+ */
 public record HitResult(
-        double t,
+        double t, @NotNull Vec3 position, @NotNull Vec3 normal, @NotNull Hittable target,
-        @NotNull Vec3 position,
+        @NotNull Material material, double u, double v, boolean isFrontFace
-        @NotNull Vec3 normal,
-        @NotNull Material material,
-        boolean frontFace
 ) implements Comparable<HitResult> {
-    public HitResult {
+
-        if (t < 0 || !Double.isFinite(t)) throw new IllegalArgumentException("t must be non-negative");
+    public @NotNull HitResult withPositionAndNormal(@NotNull Vec3 position, @NotNull Vec3 normal) {
-        Objects.requireNonNull(position, "position");
+        return new HitResult(t, position, normal, target, material, u, v, isFrontFace);
-        normal = normal.unit();
     }
 
     @Override

src/main/java/eu/jonahbauer/raytracing/scene/Hittable.java

@@ -1,17 +1,49 @@
 package eu.jonahbauer.raytracing.scene;
 
+import eu.jonahbauer.raytracing.math.AABB;
 import eu.jonahbauer.raytracing.math.Range;
 import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.transform.RotateY;
+import eu.jonahbauer.raytracing.scene.transform.Translate;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.Optional;
+import java.util.random.RandomGenerator;
 
 public interface Hittable {
+    @NotNull Range FORWARD = new Range(0.001, Double.POSITIVE_INFINITY);
 
     /**
-     * {@return the value <code>t</code> such that <code>ray.at(t)</code> is the intersection of this shaped closest to
+     * @see #hit(Ray, Range, RandomGenerator)
-     * the ray origin, or <code>Double.NaN</code> if the ray does not intersect this shape}
-     * @param ray a ray
      */
-    @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range);
+    default @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull RandomGenerator random) {
+        return hit(ray, FORWARD, random);
+    }
+
+    /**
+     * Tests whether the {@code ray} intersects {@code this} hittable.
+     * <p>
+     * The second parameter {@code range} allows the implementation to skip unnecessary calculations if it can
+     * determine that a hit (if any) will fall outside the valid range of {@code t}s. The returned hit may still be
+     * outside the valid range and has to be checked by the caller.
+     * @param ray a ray
+     * @param range the range of valid {@code t}s
+     * @return the result of the hit test, containing (among others) the value {@code t} such that {@code ray.at(t)} is
+     *         a point on {@code this} hittable
+     */
+    @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @NotNull RandomGenerator random);
+
+    /**
+     * {@return the axis-aligned bounding box of this hittable}
+     */
+    @NotNull AABB getBoundingBox();
+
+    default @NotNull Hittable translate(@NotNull Vec3 offset) {
+        return Translate.create(this, offset);
+    }
+
+    default @NotNull Hittable rotateY(double angle) {
+        return RotateY.create(this, angle);
+    }
 }

src/main/java/eu/jonahbauer/raytracing/scene/Scene.java

@@ -1,31 +1,62 @@
 package eu.jonahbauer.raytracing.scene;
 
-import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.AABB;
 import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.Spectra;
+import eu.jonahbauer.raytracing.scene.util.HittableBinaryTree;
+import eu.jonahbauer.raytracing.scene.util.HittableCollection;
 import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
 
+import java.util.ArrayList;
 import java.util.List;
-import java.util.Optional;
+import java.util.Objects;
+import java.util.random.RandomGenerator;
 
-public record Scene(@NotNull List<@NotNull Hittable> objects) implements Hittable {
+public final class Scene extends HittableCollection {
+    private final @NotNull HittableCollection objects;
+    private final @NotNull SkyBox background;
 
-    public Scene {
+    private final @Nullable List<@NotNull Target> targets;
-        objects = List.copyOf(objects);
+
+    public Scene(@NotNull List<? extends @NotNull Hittable> objects) {
+        this(objects, null);
     }
 
-    public Scene(@NotNull Hittable @NotNull ... objects) {
+    public Scene(@NotNull List<? extends @NotNull Hittable> objects, @Nullable List<? extends @NotNull Target> targets) {
-        this(List.of(objects));
+        this(Spectra.BLACK, objects, targets);
     }
 
-    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range) {
+    public Scene(@NotNull SkyBox background, @NotNull List<? extends @NotNull Hittable> objects) {
-        var result = (HitResult) null;
+        this(background, objects, null);
-        for (var object : objects) {
+    }
-            var r = object.hit(ray, range);
+
-            if (r.isPresent() && range.surrounds(r.get().t())) {
+    public Scene(@NotNull SkyBox background, @NotNull List<? extends @NotNull Hittable> objects, @Nullable List<? extends @NotNull Target> targets) {
-                result = r.get();
+        var list = new ArrayList<Hittable>(objects.size() + (targets != null ? targets.size() : 0));
-                range = new Range(range.min(), result.t());
+        list.addAll(objects);
-            }
+        if (targets != null) list.addAll(targets);
-        }
+        this.objects = new HittableBinaryTree(list);
-        return Optional.ofNullable(result);
+        this.background = Objects.requireNonNull(background);
+
+        this.targets = targets != null ? List.copyOf(targets) : null;
+    }
+
+    @Override
+    public void hit(@NotNull Ray ray, @NotNull State state, @NotNull RandomGenerator random) {
+        objects.hit(ray, state, random);
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return objects.getBoundingBox();
+    }
+
+    public @Nullable List<@NotNull Target> getTargets() {
+        return targets;
+    }
+
+    public @NotNull SampledSpectrum getBackgroundColor(@NotNull Ray ray) {
+        return background.getColor(ray);
     }
 }

src/main/java/eu/jonahbauer/raytracing/scene/SkyBox.java

@@ -0,0 +1,27 @@
+package eu.jonahbauer.raytracing.scene;
+
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.render.spectrum.SampledSpectrum;
+import eu.jonahbauer.raytracing.render.spectrum.Spectrum;
+import org.jetbrains.annotations.NotNull;
+
+@FunctionalInterface
+public interface SkyBox {
+    @NotNull SampledSpectrum getColor(@NotNull Ray ray);
+
+    static @NotNull SkyBox gradient(@NotNull Spectrum top, @NotNull Spectrum bottom) {
+        return ray -> {
+            // altitude from -pi/2 to pi/2
+            var alt = Math.copySign(
+                    Math.acos(ray.direction().withY(0).unit().dot(ray.direction().unit())),
+                    ray.direction().y()
+            );
+
+            return SampledSpectrum.lerp(
+                    top.sample(ray.lambda()),
+                    bottom.sample(ray.lambda()),
+                    alt / Math.PI + 0.5
+            );
+        };
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/Sphere.java

@@ -1,58 +0,0 @@
-package eu.jonahbauer.raytracing.scene;
-
-import eu.jonahbauer.raytracing.material.Material;
-import eu.jonahbauer.raytracing.math.Range;
-import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
-import org.jetbrains.annotations.NotNull;
-
-import java.util.Objects;
-import java.util.Optional;
-
-public record Sphere(@NotNull Vec3 center, double radius, @NotNull Material material) implements Hittable {
-
-    public Sphere {
-        Objects.requireNonNull(center, "center");
-        Objects.requireNonNull(material, "material");
-        if (radius <= 0 || !Double.isFinite(radius)) throw new IllegalArgumentException("radius must be positive");
-    }
-
-    public Sphere(double x, double y, double z, double r, @NotNull Material material) {
-        this(new Vec3(x, y, z), r, material);
-    }
-
-    @Override
-    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range) {
-        var oc = ray.origin().minus(center());
-
-        var a = ray.direction().squared();
-        var h = ray.direction().times(oc);
-        var c = oc.squared() - radius * radius;
-
-        var discriminant = h * h - a * c;
-        if (discriminant < 0) return Optional.empty();
-
-        var sd = Math.sqrt(discriminant);
-
-        double t = (- h - sd) / a;
-        if (!range.surrounds(t)) t = (- h + sd) / a;
-        if (!range.surrounds(t)) return Optional.empty();
-
-        var position = ray.at(t);
-        var normal = position.minus(center);
-        var frontFace = normal.times(ray.direction()) < 0;
-        return Optional.of(new HitResult(t, position, frontFace ? normal : normal.times(-1), material, frontFace));
-    }
-
-    public @NotNull Sphere withCenter(@NotNull Vec3 center) {
-        return new Sphere(center, radius, material);
-    }
-
-    public @NotNull Sphere withCenter(double x, double y, double z) {
-        return withCenter(new Vec3(x, y, z));
-    }
-
-    public @NotNull Sphere withRadius(double radius) {
-        return new Sphere(center, radius, material);
-    }
-}

src/main/java/eu/jonahbauer/raytracing/scene/Target.java

@@ -0,0 +1,39 @@
+package eu.jonahbauer.raytracing.scene;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.renderer.pdf.TargetingProbabilityDensityFunction;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+/**
+ * An interface for objects that can be targeted. A target can construct randomly distributed directions in which
+ * it will be hit from a given origin.
+ * @see TargetingProbabilityDensityFunction
+ */
+public interface Target extends Hittable {
+    /**
+     * Returns the probability density for a direction as sampled by {@link #getTargetingDirection(Vec3, RandomGenerator)}.
+     * @param origin the origin
+     * @param direction the direction
+     * @return the probability density for a direction as sampled by {@link #getTargetingDirection(Vec3, RandomGenerator)}
+     */
+    double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction);
+
+    /**
+     * {@return a vector targeting this hittable from the <code>origin</code>} The vector is chosen randomly.
+     * @param origin the origin
+     * @param random a random number generator
+     */
+    @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random);
+
+    @Override
+    default @NotNull Target translate(@NotNull Vec3 offset) {
+        return (Target) Hittable.super.translate(offset);
+    }
+
+    @Override
+    default @NotNull Target rotateY(double angle) {
+        return (Target) Hittable.super.rotateY(angle);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable2d/Ellipse.java

@@ -0,0 +1,25 @@
+package eu.jonahbauer.raytracing.scene.hittable2d;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.Material;
+import org.jetbrains.annotations.NotNull;
+
+public final class Ellipse extends Hittable2D {
+    private final @NotNull AABB bbox;
+
+    public Ellipse(@NotNull Vec3 origin, @NotNull Vec3 u, @NotNull Vec3 v, @NotNull Material material) {
+        super(origin, u, v, material);
+        this.bbox = new AABB(origin.minus(u).minus(v), origin.plus(u).plus(v));
+    }
+
+    @Override
+    protected boolean isInterior(double alpha, double beta) {
+        return alpha * alpha + beta * beta < 1;
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable2d/Hittable2D.java

@@ -0,0 +1,93 @@
+package eu.jonahbauer.raytracing.scene.hittable2d;
+
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public abstract class Hittable2D implements Hittable {
+    protected final @NotNull Vec3 origin;
+    protected final @NotNull Vec3 u;
+    protected final @NotNull Vec3 v;
+    private final @NotNull Material material;
+
+    // internal
+    private final @NotNull Vec3 normal;
+    private final double d;
+    private final @NotNull Vec3 w;
+
+    protected Hittable2D(@NotNull Vec3 origin, @NotNull Vec3 u, @NotNull Vec3 v, @NotNull Material material) {
+        this.origin = Objects.requireNonNull(origin);
+        this.u = Objects.requireNonNull(u);
+        this.v = Objects.requireNonNull(v);
+        this.material = Objects.requireNonNull(material);
+
+        var n = u.cross(v);
+        if (n.squared() < 1e-8) throw new IllegalArgumentException();
+        this.normal = n.unit();
+        this.d = origin.dot(normal);
+        this.w = n.div(n.squared());
+    }
+
+    @Override
+    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @Nullable RandomGenerator random) {
+        var denominator = ray.direction().dot(normal);
+        if (Math.abs(denominator) < 1e-8) return Optional.empty(); // parallel
+
+        var t = (d - ray.origin().dot(normal)) / denominator;
+        if (!range.surrounds(t)) return Optional.empty();
+
+        var position = ray.at(t);
+        var p = position.minus(origin);
+
+        var alpha = w.dot(p.cross(v));
+        var beta = w.dot(u.cross(p));
+        if (!isInterior(alpha, beta)) return Optional.empty();
+
+        var frontFace = denominator < 0;
+        return Optional.of(new HitResult(
+                t, position, frontFace ? normal : normal.neg(), this,
+                material, alpha, beta, frontFace
+        ));
+    }
+
+    protected double hit0(@NotNull Ray ray, @NotNull Range range) {
+        var denominator = ray.direction().dot(normal);
+        if (Math.abs(denominator) < 1e-8) return Double.NaN; // parallel
+
+        var t = (d - ray.origin().dot(normal)) / denominator;
+        if (!range.surrounds(t)) return Double.NaN;
+
+        var position = ray.at(t);
+        var p = position.minus(origin);
+
+        var alpha = Vec3.tripleProduct(w, p, v);
+        var beta = Vec3.tripleProduct(w, u, p);
+        if (!isInterior(alpha, beta)) return Double.NaN;
+
+        return t;
+    }
+
+    protected @NotNull Vec3 get(double alpha, double beta) {
+        return new Vec3(
+                Math.fma(beta, v.x(), Math.fma(alpha, u.x(), origin.x())),
+                Math.fma(beta, v.y(), Math.fma(alpha, u.y(), origin.y())),
+                Math.fma(beta, v.z(), Math.fma(alpha, u.z(), origin.z()))
+        );
+    }
+
+    protected abstract boolean isInterior(double alpha, double beta);
+
+    @Override
+    public @NotNull String toString() {
+        return this.getClass().getSimpleName() + "(origin=" + origin + ", u=" + u + ", v=" + v + ")";
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable2d/Parallelogram.java

@@ -0,0 +1,50 @@
+package eu.jonahbauer.raytracing.scene.hittable2d;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.scene.Target;
+import eu.jonahbauer.raytracing.scene.util.PdfUtil;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+public final class Parallelogram extends Hittable2D implements Target {
+    private final @NotNull AABB bbox;
+
+    public Parallelogram(@NotNull Vec3 origin, @NotNull Vec3 u, @NotNull Vec3 v, @NotNull Material material) {
+        super(origin, u, v, material);
+        this.bbox = new AABB(origin, origin.plus(u).plus(v));
+    }
+
+    @Override
+    protected boolean isInterior(double alpha, double beta) {
+        return 0 <= alpha && alpha < 1 && 0 <= beta && beta < 1;
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+
+    @Override
+    public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        if (Double.isNaN(hit0(new Ray(origin, direction), FORWARD))) return 0;
+
+        var o = this.origin.minus(origin);
+        var a = o.unit();
+        var b = o.plus(u).unit();
+        var c = o.plus(v).unit();
+        var d = o.plus(u).plus(v).unit();
+        var angle = PdfUtil.getSolidAngle(a, b, d) + PdfUtil.getSolidAngle(c, b, d);
+        return 1 / angle;
+    }
+
+    @Override
+    public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+        var alpha = random.nextDouble();
+        var beta = random.nextDouble();
+        return get(alpha, beta).minus(origin);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable2d/Triangle.java

@@ -0,0 +1,25 @@
+package eu.jonahbauer.raytracing.scene.hittable2d;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.Material;
+import org.jetbrains.annotations.NotNull;
+
+public final class Triangle extends Hittable2D {
+    private final @NotNull AABB bbox;
+
+    public Triangle(@NotNull Vec3 origin, @NotNull Vec3 u, @NotNull Vec3 v, @NotNull Material material) {
+        super(origin, u, v, material);
+        this.bbox = new AABB(origin, origin.plus(u).plus(v));
+    }
+
+    @Override
+    protected boolean isInterior(double alpha, double beta) {
+        return 0 <= alpha && 0 <= beta && alpha + beta <= 1;
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable3d/Box.java

@@ -0,0 +1,304 @@
+package eu.jonahbauer.raytracing.scene.hittable3d;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import eu.jonahbauer.raytracing.scene.Target;
+import eu.jonahbauer.raytracing.scene.util.PdfUtil;
+import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public final class Box implements Hittable, Target {
+    private final @NotNull AABB box;
+    private final @Nullable Material @NotNull[] materials;
+
+    public Box(@NotNull Vec3 a, @NotNull Vec3 b, @NotNull Material material) {
+        this(new AABB(a, b), material);
+    }
+
+    public Box(@NotNull AABB box, @NotNull Material material) {
+        this(box, Objects.requireNonNull(material, "material"), material, material, material, material, material);
+    }
+
+    public Box(
+            @NotNull Vec3 a, @NotNull Vec3 b,
+            @Nullable Material top, @Nullable Material bottom,
+            @Nullable Material left, @Nullable Material right,
+            @Nullable Material front, @Nullable Material back
+    ) {
+        this(new AABB(a, b), top, bottom, left, right, front, back);
+    }
+
+    public Box(
+            @NotNull AABB box,
+            @Nullable Material top, @Nullable Material bottom,
+            @Nullable Material left, @Nullable Material right,
+            @Nullable Material front, @Nullable Material back
+    ) {
+        this.box = Objects.requireNonNull(box, "box");
+        this.materials = new Material[] { left, bottom, back, right, top, front };
+    }
+
+    @Override
+    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @Nullable RandomGenerator random) {
+        // based on AABB#hit with additional detection of the side hit
+        var origin = ray.origin();
+        var direction = ray.direction();
+        var invDirection = direction.inv();
+
+        var tmin = AABB.intersect(box.min(), origin, invDirection);
+        var tmax = AABB.intersect(box.max(), origin, invDirection);
+
+        double tlmax = Double.NEGATIVE_INFINITY;
+        double tumin = Double.POSITIVE_INFINITY;
+
+        Side entry = null;
+        Side exit = null;
+
+        for (int i = 0; i < 3; i++) {
+            if (direction.get(i) >= 0) {
+                if (tmin[i] > tlmax) {
+                    tlmax = tmin[i];
+                    entry = Side.NEGATIVE[i];
+                }
+                if (tmax[i] < tumin) {
+                    tumin = tmax[i];
+                    exit = Side.POSITIVE[i];
+                }
+            } else {
+                if (tmax[i] > tlmax) {
+                    tlmax = tmax[i];
+                    entry = Side.POSITIVE[i];
+                }
+                if (tmin[i] < tumin) {
+                    tumin = tmin[i];
+                    exit = Side.NEGATIVE[i];
+                }
+            }
+        }
+
+        if (tlmax < tumin && tumin >= range.min() && tlmax <= range.max()) {
+            assert entry != null && exit != null;
+            return hit0(tlmax, tumin, entry, exit, ray, range);
+        } else {
+            return Optional.empty();
+        }
+    }
+
+    private @NotNull Optional<HitResult> hit0(double tmin, double tmax, @NotNull Side entry, @NotNull Side exit, @NotNull Ray ray, @NotNull Range range) {
+        double t;
+        Side side;
+        boolean frontFace;
+        Material material;
+        Vec3 normal;
+        if (range.surrounds(tmin) && materials[entry.ordinal()] != null) {
+            t = tmin;
+            side = entry;
+            frontFace = true;
+            material = materials[entry.ordinal()];
+            normal = side.normal;
+        } else if (range.surrounds(tmax) && materials[exit.ordinal()] != null) {
+            t = tmax;
+            side = exit;
+            frontFace = false;
+            material = materials[exit.ordinal()];
+            normal = side.normal.neg();
+        } else {
+            return Optional.empty();
+        }
+
+        var position = ray.at(t);
+        var uv = material.texture().isUVRequired();
+        var u = uv ? side.getTextureU(box, position) : Double.NaN;
+        var v = uv ? side.getTextureV(box, position) : Double.NaN;
+        return Optional.of(new HitResult(t, position, normal, this, material, u, v, frontFace));
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return box;
+    }
+
+    @Override
+    public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        if (contains(origin)) return 1 / (4 * Math.PI);
+        if (hit(new Ray(origin, direction), null).isEmpty()) return 0;
+
+        var solidAngle = 0d;
+        for (var s : Side.values()) {
+            if (!s.isExterior(box, origin)) continue;
+            solidAngle += s.getSolidAngle(box, origin);
+        }
+
+        return 1 / solidAngle;
+    }
+
+    @Override
+    public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+        if (contains(origin)) return Vec3.random(random);
+
+        // determine sides facing the origin and their solid angles
+        int visible = 0;
+
+        // at most three faces are visible
+        Side[] sides = new Side[3];
+        double[] solidAngle = new double[3];
+
+        double accumSolidAngle = 0;
+        for (var s : Side.values()) {
+            if (!s.isExterior(box, origin)) continue;
+
+            var sa = s.getSolidAngle(box, origin);
+            accumSolidAngle += sa;
+            sides[visible] = s;
+            solidAngle[visible] = accumSolidAngle;
+            visible++;
+        }
+
+        // choose a random side facing the origin based on their relative solid angles
+        var r = random.nextDouble() * solidAngle[visible - 1];
+        for (int j = 0; j < visible; j++) {
+            if (r < solidAngle[j]) {
+                // choose a random point on that side
+                var target = sides[j].random(box, random);
+                return target.minus(origin);
+            }
+        }
+
+        throw new AssertionError();
+    }
+
+    private boolean contains(@NotNull Vec3 point) {
+        return box.min().x() < point.x() && point.x() < box.max().x()
+               && box.min().y() < point.y() && point.y() < box.max().y()
+               && box.min().z() < point.z() && point.z() < box.max().z();
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return "Box(min=" + box.min() + ", max=" + box.max() + ")";
+    }
+
+    private enum Side {
+        NEG_X(Vec3.UNIT_X.neg()),
+        NEG_Y(Vec3.UNIT_Y.neg()),
+        NEG_Z(Vec3.UNIT_Z.neg()),
+        POS_X(Vec3.UNIT_X),
+        POS_Y(Vec3.UNIT_Y),
+        POS_Z(Vec3.UNIT_Z),
+        ;
+
+        private static final Side[] NEGATIVE = new Side[] {Side.NEG_X, Side.NEG_Y, Side.NEG_Z};
+        private static final Side[] POSITIVE = new Side[] {Side.POS_X, Side.POS_Y, Side.POS_Z};
+
+        private final @NotNull Vec3 normal;
+
+        Side(@NotNull Vec3 normal) {
+            this.normal = Objects.requireNonNull(normal, "normal");
+        }
+
+        /**
+         * {@return the texture u coordinate for a position on this side of the box}
+         */
+        public double getTextureU(@NotNull AABB box, @NotNull Vec3 pos) {
+            return switch (this) {
+                case NEG_X -> (pos.z() - box.min().z()) / (box.max().z() - box.min().z());
+                case POS_X -> (box.max().z() - pos.z()) / (box.max().z() - box.min().z());
+                case NEG_Y, POS_Y, POS_Z -> (pos.x() - box.min().x()) / (box.max().x() - box.min().x());
+                case NEG_Z -> (box.max().x() - pos.x()) / (box.max().x() - box.min().x());
+            };
+        }
+
+        /**
+         * {@return the texture v coordinate for a position on this side of the box}
+         */
+        public double getTextureV(@NotNull AABB box, @NotNull Vec3 pos) {
+            return switch (this) {
+                case NEG_X, POS_X, NEG_Z, POS_Z -> (pos.y() - box.min().y()) / (box.max().y() - box.min().y());
+                case NEG_Y -> (pos.z() - box.min().z()) / (box.max().z() - box.min().z());
+                case POS_Y -> (box.max().z() - pos.z()) / (box.max().z() - box.min().z());
+            };
+        }
+
+        /**
+         * {@return whether the given position is outside of the box only considering <code>this</code> side}
+         */
+        public boolean isExterior(@NotNull AABB box, @NotNull Vec3 pos) {
+            return switch (this) {
+                case NEG_X -> pos.x() < box.min().x();
+                case NEG_Y -> pos.y() < box.min().y();
+                case NEG_Z -> pos.z() < box.min().z();
+                case POS_X -> pos.x() > box.max().x();
+                case POS_Y -> pos.y() > box.max().y();
+                case POS_Z -> pos.z() > box.max().z();
+            };
+        }
+
+        /**
+         * {@return the point on <code>this</code> side of the <code>box</code> with the given <code>u</code>-<code>v</code>-coordinates}
+         */
+        public @NotNull Vec3 get(@NotNull AABB box, double u, double v) {
+            return switch (this) {
+                case NEG_X -> new Vec3(
+                        box.min().x(),
+                        Math.fma(v, box.max().y() - box.min().y(), box.min().y()),
+                        Math.fma(u, box.max().z() - box.min().z(), box.min().z())
+                );
+                case NEG_Y -> new Vec3(
+                        Math.fma(u, box.max().x() - box.min().x(), box.min().x()),
+                        box.min().y(),
+                        Math.fma(v, box.max().z() - box.min().z(), box.min().z())
+                );
+                case NEG_Z -> new Vec3(
+                        Math.fma(u, box.min().x() - box.max().x(), box.max().x()),
+                        Math.fma(v, box.max().y() - box.min().y(), box.min().y()),
+                        box.min().z()
+                );
+                case POS_X -> new Vec3(
+                        box.max().x(),
+                        Math.fma(v, box.max().y() - box.min().y(), box.min().y()),
+                        Math.fma(u, box.min().z() - box.max().z(), box.max().z())
+                );
+                case POS_Y -> new Vec3(
+                        Math.fma(u, box.max().x() - box.min().x(), box.min().x()),
+                        box.max().y(),
+                        Math.fma(v, box.min().z() - box.max().z(), box.max().z())
+                );
+                case POS_Z -> new Vec3(
+                        Math.fma(u, box.max().x() - box.min().x(), box.min().x()),
+                        Math.fma(v, box.max().y() - box.min().y(), box.min().y()),
+                        box.max().z()
+                );
+            };
+        }
+
+        /**
+         * {@return a random point on <code>this</code> side of the <code>box</code>}
+         */
+        public @NotNull Vec3 random(@NotNull AABB box, @NotNull RandomGenerator random) {
+            var u = random.nextDouble();
+            var v = random.nextDouble();
+            return get(box, u, v);
+        }
+
+        /**
+         * {@return the solid angle covered by <code>this</code> side of the <code>box</code> when viewed from <code>pos</code>}
+         */
+        public double getSolidAngle(@NotNull AABB box, @NotNull Vec3 pos) {
+            var a = get(box, 0, 0).minus(pos).unit();
+            var b = get(box, 0, 1).minus(pos).unit();
+            var c = get(box, 1, 1).minus(pos).unit();
+            var d = get(box, 1, 0).minus(pos).unit();
+
+            return PdfUtil.getSolidAngle(a, b, d) + PdfUtil.getSolidAngle(c, b, d);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable3d/ConstantMedium.java

@@ -0,0 +1,43 @@
+package eu.jonahbauer.raytracing.scene.hittable3d;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.material.IsotropicMaterial;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public record ConstantMedium(@NotNull Hittable boundary, double density, @NotNull IsotropicMaterial material) implements Hittable {
+
+    @Override
+    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @NotNull RandomGenerator random) {
+        var hit1 = boundary.hit(ray, Range.UNIVERSE, random);
+        if (hit1.isEmpty()) return Optional.empty();
+
+        var hit2 = boundary.hit(ray, new Range(hit1.get().t() + 0.0001, Double.POSITIVE_INFINITY), random);
+        if (hit2.isEmpty()) return Optional.empty();
+
+        var tmin = Math.max(range.min(), hit1.get().t());
+        var tmax = Math.min(range.max(), hit2.get().t());
+        if (tmin >= tmax) return Optional.empty();
+        if (tmin < 0) tmin = 0;
+
+        var length = ray.direction().length();
+        var distance = length * (tmax - tmin);
+        var hitDistance = - Math.log(random.nextDouble()) / density;
+        if (hitDistance > distance) return Optional.empty();
+
+        var t = tmin + hitDistance / length;
+        return Optional.of(new HitResult(t, ray.at(t), Vec3.UNIT_X, this, material, 0, 0, true)); // arbitrary normal, u, v and isFrontFace
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return boundary().getBoundingBox();
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/hittable3d/Sphere.java

@@ -0,0 +1,118 @@
+package eu.jonahbauer.raytracing.scene.hittable3d;
+
+import eu.jonahbauer.raytracing.render.material.Material;
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import eu.jonahbauer.raytracing.scene.Target;
+import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public final class Sphere implements Hittable, Target {
+    private final @NotNull Vec3 center;
+    private final double radius;
+    private final @NotNull Material material;
+
+    private final @NotNull AABB bbox;
+
+    private final @NotNull Vec3 normalizedCenter;
+    private final double invRadius;
+
+    public Sphere(@NotNull Vec3 center, double radius, @NotNull Material material) {
+        this.center = Objects.requireNonNull(center, "center");
+        this.material = Objects.requireNonNull(material, "material");
+        if (radius <= 0 || !Double.isFinite(radius)) throw new IllegalArgumentException("radius must be positive");
+        this.radius = radius;
+
+        this.invRadius = 1 / radius;
+        this.normalizedCenter = this.center.times(-this.invRadius);
+
+        this.bbox = new AABB(
+                center.minus(radius, radius, radius),
+                center.plus(radius, radius, radius)
+        );
+    }
+
+    @Override
+    public @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @Nullable RandomGenerator random) {
+        var t = hit0(ray, range);
+        if (Double.isNaN(t)) return Optional.empty();
+
+        var position = ray.at(t);
+        var normal = Vec3.fma(invRadius, position, normalizedCenter);
+        var frontFace = normal.dot(ray.direction()) < 0;
+
+        double u;
+        double v;
+        if (material.texture().isUVRequired()) {
+            var theta = Math.acos(-normal.y());
+            var phi = Math.atan2(-normal.z(), normal.x()) + Math.PI;
+            u = phi / (2 * Math.PI);
+            v = theta / Math.PI;
+        } else {
+            u = Double.NaN;
+            v = Double.NaN;
+        }
+
+        return Optional.of(new HitResult(
+                t, position, frontFace ? normal : normal.neg(), this,
+                material, u, v, frontFace
+        ));
+    }
+
+    private double hit0(@NotNull Ray ray, @NotNull Range range) {
+        var oc = ray.origin().minus(center);
+
+        var a = ray.direction().squared();
+        var h = ray.direction().dot(oc);
+        var c = oc.squared() - radius * radius;
+
+        var discriminant = h * h - a * c;
+        if (discriminant < 0) return Double.NaN;
+
+        var sd = Math.sqrt(discriminant);
+
+        double t = (- h - sd) / a;
+        if (!range.surrounds(t)) t = (- h + sd) / a;
+        if (!range.surrounds(t)) return Double.NaN;
+
+        return t;
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+
+    @Override
+    public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        if (Double.isNaN(hit0(new Ray(origin, direction), FORWARD))) return 0;
+
+        var cos = Math.sqrt(1 - radius * radius / (center.minus(origin).squared()));
+        var solidAngle = 2 * Math.PI * (1 - cos);
+        return 1 / solidAngle;
+    }
+
+    @Override
+    public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+        var direction = center.minus(origin);
+        var out = Vec3.randomOppositeHemisphere(random, direction);
+        return new Vec3(
+                Math.fma(radius, out.x(), center.x() - origin.x()),
+                Math.fma(radius, out.y(), center.y() - origin.y()),
+                Math.fma(radius, out.z(), center.z() - origin.z())
+        );
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return "Sphere(center=" + center + ", radius=" + radius + ")";
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/transform/RotateY.java

@@ -0,0 +1,113 @@
+package eu.jonahbauer.raytracing.scene.transform;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import eu.jonahbauer.raytracing.scene.Target;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+public sealed class RotateY extends Transform {
+    private final double cos;
+    private final double sin;
+
+    private final @NotNull AABB bbox;
+
+    public static @NotNull RotateY create(@NotNull Hittable object, double angle) {
+        if (object instanceof Target) {
+            return new RotateYTarget(object, angle);
+        } else {
+            return new RotateY(object, angle);
+        }
+    }
+
+    private RotateY(@NotNull Hittable object, double angle) {
+        super(object);
+        this.cos = Math.cos(angle);
+        this.sin = Math.sin(angle);
+
+        var bbox = object.getBoundingBox();
+
+        var min = new Vec3(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE);
+        var max = new Vec3(- Double.MAX_VALUE, - Double.MAX_VALUE, - Double.MAX_VALUE);
+
+        for (int i = 0; i < 2; i++) {
+            for (int j = 0; j < 2; j++) {
+                for (int k = 0; k < 2; k++) {
+                    var x = i * bbox.max().x() + (1 - i) * bbox.min().x();
+                    var y = j * bbox.max().y() + (1 - j) * bbox.min().y();
+                    var z = k * bbox.max().z() + (1 - k) * bbox.min().z();
+
+                    var newx =  cos * x + sin * z;
+                    var newz = -sin * x + cos * z;
+
+                    var temp = new Vec3(newx, y, newz);
+
+                    min = Vec3.min(min, temp);
+                    max = Vec3.max(max, temp);
+                }
+            }
+        }
+
+        this.bbox = new AABB(min, max);
+    }
+
+    @Override
+    protected final @NotNull Ray transform(@NotNull Ray ray) {
+        var origin = ray.origin();
+        var direction = ray.direction();
+
+        var newOrigin = transform(origin);
+        var newDirection = transform(direction);
+        return new Ray(newOrigin, newDirection, ray.lambda());
+    }
+
+    @Override
+    protected final @NotNull HitResult transform(@NotNull HitResult result) {
+        var position = result.position();
+        var newPosition = untransform(position);
+
+        var normal = result.normal();
+        var newNormal = untransform(normal);
+
+        return result.withPositionAndNormal(newPosition, newNormal);
+    }
+
+    protected final @NotNull Vec3 transform(@NotNull Vec3 vec) {
+        return new Vec3(cos * vec.x() - sin * vec.z(), vec.y(), sin * vec.x() + cos * vec.z());
+    }
+
+    protected final @NotNull Vec3 untransform(@NotNull Vec3 vec) {
+        return new Vec3(cos * vec.x() + sin * vec.z(), vec.y(), - sin * vec.x() + cos * vec.z());
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return object + " rotated by " + Math.toDegrees(Math.atan2(sin, cos)) + "° around the y axis";
+    }
+
+    private static final class RotateYTarget extends RotateY implements Target {
+
+        private RotateYTarget(@NotNull Hittable object, double angle) {
+            super(object, angle);
+        }
+
+        @Override
+        public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+            return ((Target) object).getProbabilityDensity(transform(origin), transform(direction));
+        }
+
+        @Override
+        public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+            return untransform(((Target) object).getTargetingDirection(transform(origin), random));
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/transform/Transform.java

@@ -0,0 +1,28 @@
+package eu.jonahbauer.raytracing.scene.transform;
+
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public abstract class Transform implements Hittable {
+    protected final @NotNull Hittable object;
+
+    protected Transform(@NotNull Hittable object) {
+        this.object = Objects.requireNonNull(object);
+    }
+
+    protected abstract @NotNull Ray transform(@NotNull Ray ray);
+
+    protected abstract @NotNull HitResult transform(@NotNull HitResult result);
+
+    @Override
+    public final @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @NotNull RandomGenerator random) {
+        return object.hit(transform(ray), range, random).map(this::transform);
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/transform/Translate.java

@@ -0,0 +1,73 @@
+package eu.jonahbauer.raytracing.scene.transform;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import eu.jonahbauer.raytracing.scene.Target;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+public sealed class Translate extends Transform {
+    protected final @NotNull Vec3 offset;
+    private final @NotNull AABB bbox;
+
+    public static @NotNull Translate create(@NotNull Hittable object, @NotNull Vec3 offset) {
+        if (object instanceof Target) {
+            return new TranslateTarget(object, offset);
+        } else {
+            return new Translate(object, offset);
+        }
+    }
+
+    private Translate(@NotNull Hittable object, @NotNull Vec3 offset) {
+        super(object);
+        this.offset = offset;
+
+        var bbox = object.getBoundingBox();
+        this.bbox = new AABB(
+                bbox.min().plus(offset),
+                bbox.max().plus(offset)
+        );
+    }
+
+    @Override
+    protected final @NotNull Ray transform(@NotNull Ray ray) {
+        return new Ray(ray.origin().minus(offset), ray.direction(), ray.lambda());
+    }
+
+    @Override
+    protected final @NotNull HitResult transform(@NotNull HitResult result) {
+        return result.withPositionAndNormal(result.position().plus(offset), result.normal());
+    }
+
+    @Override
+    public final @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+
+    @Override
+    public @NotNull String toString() {
+        return object + " translated by " + offset;
+    }
+
+    private static final class TranslateTarget extends Translate implements Target {
+        private TranslateTarget(@NotNull Hittable object, @NotNull Vec3 offset) {
+            super(object, offset);
+        }
+
+        @Override
+        public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+            if (!(object instanceof Target target)) throw new UnsupportedOperationException();
+            return target.getProbabilityDensity(origin.minus(offset), direction);
+        }
+
+        @Override
+        public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+            if (!(object instanceof Target target)) throw new UnsupportedOperationException();
+            return target.getTargetingDirection(origin.minus(offset), random);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/util/HittableBinaryTree.java

@@ -0,0 +1,68 @@
+package eu.jonahbauer.raytracing.scene.util;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
+
+import java.util.Comparator;
+import java.util.List;
+import java.util.random.RandomGenerator;
+
+public final class HittableBinaryTree extends HittableCollection {
+    private final @Nullable Hittable left;
+    private final @Nullable Hittable right;
+    private final @NotNull AABB bbox;
+
+    public HittableBinaryTree(@NotNull List<? extends @NotNull Hittable> objects) {
+        bbox = AABB.getBoundingBox(objects).orElse(AABB.EMPTY);
+        if (objects.isEmpty()) {
+            left = null;
+            right = null;
+        } else if (objects.size() == 1) {
+            left = objects.getFirst();
+            right = null;
+        } else if (objects.size() == 2) {
+            left = objects.getFirst();
+            right = objects.getLast();
+        } else {
+            var x = bbox.x().size();
+            var y = bbox.y().size();
+            var z = bbox.z().size();
+            Comparator<AABB> comparator;
+            if (x > y && x > z) {
+                comparator = AABB.X_AXIS;
+            } else if (y > z) {
+                comparator = AABB.Y_AXIS;
+            } else {
+                comparator = AABB.Z_AXIS;
+            }
+            var sorted = objects.stream().sorted(Comparator.comparing(Hittable::getBoundingBox, comparator)).toList();
+            var size = sorted.size();
+
+            left = new HittableBinaryTree(sorted.subList(0, size / 2));
+            right = new HittableBinaryTree(sorted.subList(size / 2, size));
+        }
+    }
+
+    @Override
+    public void hit(@NotNull Ray ray, @NotNull State state, @NotNull RandomGenerator random) {
+        if (!bbox.hit(ray, state.getRange())) return;
+        if (left instanceof HittableCollection coll) {
+            coll.hit(ray, state, random);
+        } else if (left != null) {
+            hit(state, ray, left, random);
+        }
+        if (right instanceof HittableCollection coll) {
+            coll.hit(ray, state, random);
+        } else if (right != null) {
+            hit(state, ray, right, random);
+        }
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/util/HittableCollection.java

@@ -0,0 +1,53 @@
+package eu.jonahbauer.raytracing.scene.util;
+
+import eu.jonahbauer.raytracing.math.Range;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.scene.HitResult;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.Objects;
+import java.util.Optional;
+import java.util.random.RandomGenerator;
+
+public abstract class HittableCollection implements Hittable {
+
+    @Override
+    public final @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range, @NotNull RandomGenerator random) {
+        var state = new State(range);
+        hit(ray, state, random);
+        return state.getResult();
+    }
+
+    public abstract void hit(@NotNull Ray ray, @NotNull State state, @NotNull RandomGenerator random);
+
+    protected static boolean hit(@NotNull State state, @NotNull Ray ray, @NotNull Hittable object, @NotNull RandomGenerator random) {
+        var r = object.hit(ray, state.range, random);
+        if (r.isPresent()) {
+            if (state.range.surrounds(r.get().t())){
+                state.result = r.get();
+                state.range = new Range(state.range.min(), state.result.t());
+            }
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+    public static class State {
+        private @NotNull Range range;
+        private HitResult result;
+
+        private State(@NotNull Range range) {
+            this.range = Objects.requireNonNull(range);
+        }
+
+        public @NotNull Range getRange() {
+            return range;
+        }
+
+        private @NotNull Optional<HitResult> getResult() {
+            return Optional.ofNullable(result);
+        }
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/util/HittableList.java

@@ -0,0 +1,34 @@
+package eu.jonahbauer.raytracing.scene.util;
+
+import eu.jonahbauer.raytracing.math.AABB;
+import eu.jonahbauer.raytracing.math.Ray;
+import eu.jonahbauer.raytracing.scene.Hittable;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.random.RandomGenerator;
+
+public final class HittableList extends HittableCollection {
+    private final @NotNull List<Hittable> objects;
+    private final @NotNull AABB bbox;
+
+    public HittableList(@NotNull List<? extends @NotNull Hittable> objects) {
+        this.objects = new ArrayList<>(objects);
+        this.bbox = AABB.getBoundingBox(this.objects).orElse(AABB.EMPTY);
+    }
+
+    public HittableList(@NotNull Hittable @NotNull... objects) {
+        this(List.of(objects));
+    }
+
+    @Override
+    public void hit(@NotNull Ray ray, @NotNull State state, @NotNull RandomGenerator random) {
+        objects.forEach(object -> hit(state, ray, object, random));
+    }
+
+    @Override
+    public @NotNull AABB getBoundingBox() {
+        return bbox;
+    }
+}

src/main/java/eu/jonahbauer/raytracing/scene/util/PdfUtil.java

@@ -0,0 +1,19 @@
+package eu.jonahbauer.raytracing.scene.util;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+public final class PdfUtil {
+    private PdfUtil() {
+        throw new UnsupportedOperationException();
+    }
+
+    /**
+     * {@return the solid angle of the triangle abc as seen from the origin} The vectors {@code a}, {@code b} and {@code c}
+     * must be unit vectors.
+     */
+    public static double getSolidAngle(@NotNull Vec3 a, @NotNull Vec3 b, @NotNull Vec3 c) {
+        var angle = 2 * Math.atan(Math.abs(Vec3.tripleProduct(a, b, c)) / (1 + a.dot(b) + b.dot(c) + c.dot(a)));
+        return angle < 0 ? 2 * Math.PI + angle : angle;
+    }
+}