add targeting probability density function

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

@@ -35,6 +35,7 @@ public class Examples {
         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);
@@ -212,6 +213,38 @@ public class Examples {
         );
     }
 
+    public static @NotNull Example getCornellBoxSphere(int height) {
+        if (height <= 0) height = 600;
+
+        var red = new LambertianMaterial(new Color(.65, .05, .05));
+        var white = new LambertianMaterial(new Color(.73, .73, .73));
+        var green = new LambertianMaterial(new Color(.12, .45, .15));
+        var light = new DiffuseLight(new Color(7.0, 7.0, 7.0));
+        var aluminum = new MetallicMaterial(new Color(0.8, 0.85, 0.88));
+        var glass = new DielectricMaterial(1.5);
+
+        return new Example(
+                new Scene(
+                        new Box(
+                                new AABB(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 AABB(new Vec3(0, 0, 0), new Vec3(165, 330, 165)),
+                                white, white, white, white, white, aluminum
+                        ).rotateY(Math.toRadians(15)).translate(new Vec3(265, 0, 295)),
+                        new Sphere(new Vec3(190, 90, 190), 90, glass)
+                ),
+                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;
 

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

@@ -17,17 +17,29 @@ public record Vec3(double x, double y, double z) {
         assert Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(z) : "x, y and z must be finite";
     }
 
+    /**
+     * {@return a uniformly random vector with components in the range [-1, 1)}
+     */
     public static @NotNull Vec3 random(@NotNull RandomGenerator random) {
-        return random(random, false);
-    }
-
-    public static @NotNull Vec3 random(@NotNull RandomGenerator random, boolean unit) {
-        var vec = new Vec3(
+        return new Vec3(
                 Math.fma(2, random.nextDouble(), -1),
                 Math.fma(2, random.nextDouble(), -1),
                 Math.fma(2, random.nextDouble(), -1)
         );
-        return unit ? vec.unit() : vec;
+    }
+
+    /**
+     * {@return a uniformly random unit vector}
+     */
+    public static @NotNull Vec3 random(@NotNull RandomGenerator random, boolean unit) {
+        if (!unit) return random(random);
+        Vec3 vec;
+        double squared;
+        do {
+            vec = random(random);
+            squared = vec.squared();
+        } while (squared > 1);
+        return vec.div(Math.sqrt(squared));
     }
 
     public static @NotNull Vec3 reflect(@NotNull Vec3 vec, @NotNull Vec3 normal) {

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

@@ -18,7 +18,13 @@ public interface Camera {
 
     /**
      * 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, @NotNull RandomGenerator random);
+    @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

@@ -79,13 +79,20 @@ public final class SimpleCamera implements Camera {
 
     /**
      * {@inheritDoc}
+     * @param x {@inheritDoc}
+     * @param y {@inheritDoc}
+     * @param i {@inheritDoc}
+     * @param j {@inheritDoc}
+     * @param n {@inheritDoc}
+     * @param random {@inheritDoc}
      */
-    public @NotNull Ray cast(int x, int y, @NotNull RandomGenerator random) {
+    @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);
 
         var origin = getRayOrigin(random);
-        var target = getRayTarget(x, y, random);
+        var target = getRayTarget(x, y, i, j, n, random);
         return new Ray(origin, target.minus(origin));
     }
 
@@ -98,8 +105,8 @@ public final class SimpleCamera implements Camera {
         if (blurRadius <= 0) return origin;
 
         while (true) {
-            var du = 2 * random.nextDouble() - 1;
-            var dv = 2 * random.nextDouble() - 1;
+            var du = Math.fma(2, random.nextDouble(), -1);
+            var dv = Math.fma(2, random.nextDouble(), -1);
             if (du * du + dv * dv >= 1) continue;
 
             var ru = blurRadius * du;
@@ -116,9 +123,10 @@ public final class SimpleCamera implements Camera {
     /**
      * {@return the target vector for a ray through the given pixel} The position is randomized within the pixel.
      */
-    private @NotNull Vec3 getRayTarget(int x, int y, @NotNull RandomGenerator random) {
-        double dx = x + random.nextDouble() - 0.5;
-        double dy = y + random.nextDouble() - 0.5;
+    private @NotNull Vec3 getRayTarget(int x, int y, int i, int j, int n, @NotNull RandomGenerator random) {
+        var factor = 1d / n;
+        var dx = x + Math.fma(factor, i + random.nextDouble(), -0.5);
+        var dy = y + Math.fma(factor, j + random.nextDouble(), -0.5);
         return pixel00.plus(pixelU.times(dx)).plus(pixelV.times(dy));
     }
 

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

@@ -32,7 +32,7 @@ public record DielectricMaterial(double refractionIndex, @NotNull Texture textur
                 .orElseGet(() -> Vec3.reflect(ray.direction(), hit.normal()));
 
         var attenuation = texture.get(hit);
-        return Optional.of(new ScatterResult(new Ray(hit.position(), newDirection), attenuation));
+        return Optional.of(new SpecularScatterResult(attenuation, new Ray(hit.position(), newDirection)));
     }
 
     private double reflectance(double cos) {

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

@@ -18,6 +18,7 @@ public final class DirectionalMaterial implements Material {
     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(
@@ -39,7 +40,7 @@ public final class DirectionalMaterial implements Material {
             if (back != null) return back.scatter(ray, hit, random);
         }
         // let the ray pass through without obstruction
-        return Optional.of(new ScatterResult(new Ray(ray.at(hit.t()), ray.direction()), Color.WHITE));
+        return Optional.of(new SpecularScatterResult(Color.WHITE, new Ray(ray.at(hit.t()), ray.direction())));
     }
 
     @Override
@@ -49,7 +50,7 @@ public final class DirectionalMaterial implements Material {
         } else {
             if (back != null) return back.emitted(hit);
         }
-        return Color.BLACK;
+        return Material.super.emitted(hit);
     }
 
     private record DirectionalTexture(@Nullable Texture front, @Nullable Texture back) implements Texture {

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

@@ -1,7 +1,7 @@
 package eu.jonahbauer.raytracing.render.material;
 
 import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.render.renderer.pdf.SphereProbabilityDensityFunction;
 import eu.jonahbauer.raytracing.render.texture.Color;
 import eu.jonahbauer.raytracing.render.texture.Texture;
 import eu.jonahbauer.raytracing.scene.HitResult;
@@ -13,7 +13,7 @@ import java.util.random.RandomGenerator;
 public record IsotropicMaterial(@NotNull Color albedo) implements Material {
     @Override
     public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
-        return Optional.of(new ScatterResult(new Ray(hit.position(), Vec3.random(random, true)), albedo()));
+        return Optional.of(new PdfScatterResult(albedo(), new SphereProbabilityDensityFunction()));
     }
 
     @Override

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

@@ -1,7 +1,7 @@
 package eu.jonahbauer.raytracing.render.material;
 
 import eu.jonahbauer.raytracing.math.Ray;
-import eu.jonahbauer.raytracing.math.Vec3;
+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;
@@ -17,11 +17,7 @@ public record LambertianMaterial(@NotNull Texture texture) implements Material {
 
     @Override
     public @NotNull Optional<ScatterResult> scatter(@NotNull Ray ray, @NotNull HitResult hit, @NotNull RandomGenerator random) {
-        var newDirection = hit.normal().plus(Vec3.random(random, true));
-        if (newDirection.isNearZero()) newDirection = hit.normal();
-
         var attenuation = texture.get(hit);
-        var scattered = new Ray(hit.position(), newDirection);
-        return Optional.of(new ScatterResult(scattered, attenuation));
+        return Optional.of(new PdfScatterResult(attenuation, new CosineProbabilityDensityFunction(hit.normal())));
     }
 }

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

@@ -1,6 +1,7 @@
 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.texture.Color;
 import eu.jonahbauer.raytracing.render.texture.Texture;
 import eu.jonahbauer.raytracing.scene.HitResult;
@@ -12,18 +13,59 @@ 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 Color#BLACK black}, i.e. no emission
+     */
     default @NotNull Color emitted(@NotNull HitResult hit) {
         return Color.BLACK;
     }
 
-    record ScatterResult(@NotNull Ray ray, @NotNull Color attenuation) {
-        public ScatterResult {
-            Objects.requireNonNull(ray, "ray");
+    /**
+     * 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 Color 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 Color attenuation, @NotNull ProbabilityDensityFunction pdf) implements ScatterResult {
+        public PdfScatterResult {
+            Objects.requireNonNull(attenuation, "attenuation");
+            Objects.requireNonNull(pdf, "pdf");
         }
     }
 }

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

@@ -28,6 +28,6 @@ public record MetallicMaterial(@NotNull Texture texture, double fuzz) implements
             newDirection = newDirection.unit().plus(Vec3.random(random, true).times(fuzz));
         }
         var attenuation = texture.get(hit);
-        return Optional.of(new ScatterResult(new Ray(hit.position(), newDirection), attenuation));
+        return Optional.of(new SpecularScatterResult(attenuation, new Ray(hit.position(), newDirection)));
     }
 }

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

@@ -2,6 +2,9 @@ package eu.jonahbauer.raytracing.render.renderer;
 
 import eu.jonahbauer.raytracing.math.Range;
 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.texture.Color;
 import eu.jonahbauer.raytracing.render.camera.Camera;
 import eu.jonahbauer.raytracing.render.canvas.Canvas;
@@ -14,7 +17,7 @@ import java.util.random.RandomGenerator;
 import java.util.stream.IntStream;
 
 public final class SimpleRenderer implements Renderer {
-    private final int samplesPerPixel;
+    private final int sqrtSamplesPerPixel;
     private final int maxDepth;
     private final double gamma;
 
@@ -30,7 +33,7 @@ public final class SimpleRenderer implements Renderer {
     }
 
     private SimpleRenderer(@NotNull Builder builder) {
-        this.samplesPerPixel = builder.samplesPerPixel;
+        this.sqrtSamplesPerPixel = (int) Math.sqrt(builder.samplesPerPixel);
         this.maxDepth = builder.maxDepth;
         this.gamma = builder.gamma;
 
@@ -38,6 +41,9 @@ public final class SimpleRenderer implements Renderer {
         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()) {
@@ -45,42 +51,66 @@ public final class SimpleRenderer implements Renderer {
         }
 
         if (iterative) {
-            var random = new Random();
+            renderIterative(camera, scene, canvas);
+        } else {
+            renderNonIterative(camera, scene, canvas);
+        }
+    }
 
-            // render one sample after the other
-            for (int i = 1 ; i <= samplesPerPixel; i++) {
-                var sample = i;
+    /**
+     * 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();
+
+        // 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, random);
+                        var ray = camera.cast(x, y, sif, sjf, sqrtSamplesPerPixel, random);
                         var c = getColor(scene, ray, random);
                         canvas.set(x, y, Color.average(canvas.get(x, y), c, sample));
                     }
                 });
             }
-
-            // apply gamma correction
-            getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
-                for (int x = 0; x < camera.getWidth(); x++) {
-                    canvas.set(x, y, Color.gamma(canvas.get(x, y), gamma));
-                }
-            });
-        } else {
-            var splittable = new SplittableRandom();
-            // render one pixel after the other
-            getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
-                var random = splittable.split();
-                for (int x = 0; x < camera.getWidth(); x++) {
-                    var color = Color.BLACK;
-                    for (int i = 1; i <= samplesPerPixel; i++) {
-                        var ray = camera.cast(x, y, random);
-                        var c = getColor(scene, ray, random);
-                        color = Color.average(color, c, i);
-                    }
-                    canvas.set(x, y, Color.gamma(color, gamma));
-                }
-            });
         }
+
+        // apply gamma correction
+        getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
+            for (int x = 0; x < camera.getWidth(); x++) {
+                canvas.set(x, y, Color.gamma(canvas.get(x, y), gamma));
+            }
+        });
+    }
+
+    /**
+     * 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();
+        // render one pixel after the other
+        getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
+            var random = splittable.split();
+            for (int x = 0; x < camera.getWidth(); x++) {
+                var color = Color.BLACK;
+                int i = 0;
+                for (int sj = 0; sj < sqrtSamplesPerPixel; sj++) {
+                    for (int si = 0; si < sqrtSamplesPerPixel; si++) {
+                        var ray = camera.cast(x, y, si, sj, sqrtSamplesPerPixel, random);
+                        var c = getColor(scene, ray, random);
+                        color = Color.average(color, c, ++i);
+                    }
+                }
+                canvas.set(x, y, Color.gamma(color, gamma));
+            }
+        });
     }
 
     /**
@@ -104,12 +134,29 @@ public final class SimpleRenderer implements Renderer {
             var hit = optional.get();
             var material = hit.material();
             var emitted = material.emitted(hit);
-            var scatter = material.scatter(ray, hit, random);
+            var result = material.scatter(ray, hit, random);
             color = Color.add(color, Color.multiply(attenuation, emitted));
 
-            if (scatter.isEmpty()) break;
-            attenuation = Color.multiply(attenuation, scatter.get().attenuation());
-            ray = scatter.get().ray();
+            if (result.isEmpty()) break;
+
+            switch (result.get()) {
+                case Material.SpecularScatterResult(var a, var scattered) -> {
+                    attenuation = Color.multiply(attenuation, a);
+                    ray = scattered;
+                }
+                case Material.PdfScatterResult(var a, var pdf) -> {
+                    if (scene.getLights() == null) {
+                        attenuation = Color.multiply(attenuation, a);
+                        ray = new Ray(hit.position(), pdf.generate(random));
+                    } else {
+                        var mixed = new MixtureProbabilityDensityFunction(new TargetingProbabilityDensityFunction(hit.position(), scene.getLights()), pdf);
+                        var direction = mixed.generate(random);
+                        var factor = pdf.value(direction) / mixed.value(direction);
+                        attenuation = Color.multiply(attenuation, Color.multiply(a, factor));
+                        ray = new Ray(hit.position(), direction);
+                    }
+                }
+            }
         }
 
         return color;

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.times(direction.unit());
+        return Math.max(0, cos / Math.PI);
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        var out = normal().plus(Vec3.random(random, true));
+        return out.isNearZero() ? normal() : out;
+    }
+}

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

@@ -0,0 +1,37 @@
+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 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) {
+        return weight * a.value(direction) + (1 - weight) * b.value(direction);
+    }
+
+    @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,11 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+public interface ProbabilityDensityFunction {
+    double value(@NotNull Vec3 direction);
+    @NotNull Vec3 generate(@NotNull RandomGenerator random);
+}

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

@@ -0,0 +1,19 @@
+package eu.jonahbauer.raytracing.render.renderer.pdf;
+
+import eu.jonahbauer.raytracing.math.Vec3;
+import org.jetbrains.annotations.NotNull;
+
+import java.util.random.RandomGenerator;
+
+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, true);
+    }
+}

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

@@ -0,0 +1,28 @@
+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.Objects;
+import java.util.random.RandomGenerator;
+
+/**
+ * A probability density function targeting a target.
+ */
+public record TargetingProbabilityDensityFunction(@NotNull Vec3 origin, @NotNull Target target) implements ProbabilityDensityFunction {
+    public TargetingProbabilityDensityFunction {
+        Objects.requireNonNull(origin, "origin");
+        Objects.requireNonNull(target, "target");
+    }
+
+    @Override
+    public double value(@NotNull Vec3 direction) {
+        return target.getProbabilityDensity(origin, direction);
+    }
+
+    @Override
+    public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
+        return target.getTargetingDirection(origin, random);
+    }
+}

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

@@ -26,6 +26,10 @@ public record Color(double r, double g, double b) implements Texture {
         return new Color(a.r() * b.r(), a.g() * b.g(), a.b() * b.b());
     }
 
+    public static @NotNull Color multiply(@NotNull Color a, double b) {
+        return new Color(a.r() * b, a.g() * b, a.b() * b);
+    }
+
     public static @NotNull Color add(@NotNull Color a, @NotNull Color b) {
         return new Color(a.r() + b.r(), a.g() + b.g(), a.b() + b.b());
     }
@@ -44,10 +48,11 @@ public record Color(double r, double g, double b) implements Texture {
     }
 
     public static @NotNull Color average(@NotNull Color current, @NotNull Color next, int index) {
+        var factor = 1d / index;
         return new Color(
-                current.r() + (next.r() - current.r()) / index,
-                current.g() + (next.g() - current.g()) / index,
-                current.b() + (next.b() - current.b()) / index
+                Math.fma(factor, next.r() - current.r(), current.r()),
+                Math.fma(factor, next.g() - current.g(), current.g()),
+                Math.fma(factor, next.b() - current.b(), current.b())
         );
     }
 

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

@@ -19,6 +19,9 @@ public interface Hittable {
      */
     @NotNull Optional<HitResult> hit(@NotNull Ray ray, @NotNull Range range);
 
+    /**
+     * {@return the axis-aligned bounding box of this hittable}
+     */
     @NotNull AABB getBoundingBox();
 
     default @NotNull Hittable translate(@NotNull Vec3 offset) {

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

@@ -6,6 +6,7 @@ import eu.jonahbauer.raytracing.render.texture.Color;
 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.List;
 import java.util.Objects;
@@ -14,6 +15,8 @@ public final class Scene extends HittableCollection {
     private final @NotNull HittableCollection objects;
     private final @NotNull SkyBox background;
 
+    private final @Nullable Target light;
+
     public Scene(@NotNull List<? extends @NotNull Hittable> objects) {
         this(Color.BLACK, objects);
     }
@@ -25,6 +28,8 @@ public final class Scene extends HittableCollection {
     public Scene(@NotNull SkyBox background, @NotNull List<? extends @NotNull Hittable> objects) {
         this.objects = new HittableBinaryTree(objects);
         this.background = Objects.requireNonNull(background);
+
+        this.light = (Target) objects.get(1);
     }
 
     public Scene(@NotNull Hittable @NotNull... objects) {
@@ -49,6 +54,10 @@ public final class Scene extends HittableCollection {
         return objects.getBoundingBox();
     }
 
+    public @Nullable Target getLights() {
+        return light;
+    }
+
     public @NotNull Color getBackgroundColor(@NotNull Ray ray) {
         return background.getColor(ray);
     }

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/Parallelogram.java

@@ -1,11 +1,17 @@
 package eu.jonahbauer.raytracing.scene.hittable2d;
 
 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.Target;
+import eu.jonahbauer.raytracing.scene.util.PdfUtil;
 import org.jetbrains.annotations.NotNull;
 
-public final class Parallelogram extends Hittable2D {
+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) {
@@ -22,4 +28,24 @@ public final class Parallelogram extends Hittable2D {
     public @NotNull AABB getBoundingBox() {
         return bbox;
     }
+
+    @Override
+    public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        var result = hit(new Ray(origin, direction), new Range(0.001, Double.POSITIVE_INFINITY));
+        if (result.isEmpty()) return 0;
+
+        var a = this.origin;
+        var b = this.origin.plus(u);
+        var c = this.origin.plus(v);
+        var d = b.plus(v);
+        var angle = PdfUtil.getSolidAngle(origin, a, b, d) + PdfUtil.getSolidAngle(origin, 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 this.origin.plus(u.times(alpha)).plus(v.times(beta)).minus(origin);
+    }
 }

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

@@ -7,13 +7,16 @@ 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 {
+public final class Box implements Hittable, Target {
     private final @NotNull AABB box;
     private final @Nullable Material @NotNull[] materials;
 
@@ -115,6 +118,61 @@ public final class Box implements Hittable {
         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), new Range(0.001, Double.POSITIVE_INFINITY)).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, true);
+
+        // 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();
+    }
+
     private enum Side {
         NEG_X(Vec3.UNIT_X.neg()),
         NEG_Y(Vec3.UNIT_Y.neg()),
@@ -155,5 +213,78 @@ public final class Box implements Hittable {
                 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);
+            var b = get(box, 0, 1);
+            var c = get(box, 1, 1);
+            var d = get(box, 1, 0);
+
+            return PdfUtil.getSolidAngle(pos, a, b, d) + PdfUtil.getSolidAngle(pos, c, b, d);
+        }
     }
 }

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

@@ -7,12 +7,14 @@ 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.Objects;
 import java.util.Optional;
+import java.util.random.RandomGenerator;
 
-public final class Sphere implements Hittable {
+public final class Sphere implements Hittable, Target {
     private final @NotNull Vec3 center;
     private final double radius;
     private final @NotNull Material material;
@@ -74,4 +76,25 @@ public final class Sphere implements Hittable {
     public @NotNull AABB getBoundingBox() {
         return bbox;
     }
+
+    @Override
+    public double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        if (hit(new Ray(origin, direction), new Range(0.001, Double.POSITIVE_INFINITY)).isEmpty()) 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);
+
+        Vec3 target;
+        do {
+            target = Vec3.random(random, true);
+        } while (target.times(direction) >= 0);
+
+        return target.times(radius).plus(center).minus(origin);
+    }
 }

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

@@ -5,8 +5,11 @@ 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 final class RotateY extends Transform {
     private final double cos;
     private final double sin;
@@ -49,41 +52,44 @@ public final class RotateY extends Transform {
         var origin = ray.origin();
         var direction = ray.direction();
 
-        var newOrigin = new Vec3(
-                cos * origin.x() - sin * origin.z(),
-                origin.y(),
-                sin * origin.x() + cos * origin.z()
-        );
-        var newDirection = new Vec3(
-                cos * direction.x() - sin * direction.z(),
-                direction.y(),
-                sin * direction.x() + cos * direction.z()
-        );
-
+        var newOrigin = transform(origin);
+        var newDirection = transform(direction);
         return new Ray(newOrigin, newDirection);
     }
 
     @Override
     protected @NotNull HitResult transform(@NotNull HitResult result) {
         var position = result.position();
-        var newPosition = new Vec3(
-                cos * position.x() + sin * position.z(),
-                position.y(),
-                - sin * position.x() + cos * position.z()
-        );
+        var newPosition = untransform(position);
 
         var normal = result.normal();
-        var newNormal = new Vec3(
-                cos * normal.x() + sin * normal.z(),
-                normal.y(),
-                -sin * normal.x() + cos * normal.z()
-        );
+        var newNormal = untransform(normal);
 
         return result.withPositionAndNormal(newPosition, newNormal);
     }
 
+    private @NotNull Vec3 transform(@NotNull Vec3 vec) {
+        return new Vec3(cos * vec.x() - sin * vec.z(), vec.y(), sin * vec.x() + cos * vec.z());
+    }
+
+    private @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 double getProbabilityDensity(@NotNull Vec3 origin, @NotNull Vec3 direction) {
+        if (!(object instanceof Target target)) throw new UnsupportedOperationException();
+        return target.getProbabilityDensity(transform(origin), transform(direction));
+    }
+
+    @Override
+    public @NotNull Vec3 getTargetingDirection(@NotNull Vec3 origin, @NotNull RandomGenerator random) {
+        if (!(object instanceof Target target)) throw new UnsupportedOperationException();
+        return untransform(target.getTargetingDirection(transform(origin), random));
+    }
 }

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

@@ -2,14 +2,17 @@ package eu.jonahbauer.raytracing.scene.transform;
 
 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 java.util.Objects;
 import java.util.Optional;
+import java.util.random.RandomGenerator;
 
-public abstract class Transform implements Hittable {
+public abstract class Transform implements Hittable, Target {
     protected final @NotNull Hittable object;
 
     protected Transform(@NotNull Hittable object) {

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

@@ -5,8 +5,11 @@ 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 final class Translate extends Transform {
     private final @NotNull Vec3 offset;
     private final @NotNull AABB bbox;
@@ -36,4 +39,16 @@ public final class Translate extends Transform {
     public @NotNull AABB getBoundingBox() {
         return bbox;
     }
+
+    @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/PdfUtil.java

@@ -0,0 +1,18 @@
+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();
+    }
+
+    public static double getSolidAngle(@NotNull Vec3 o, @NotNull Vec3 a, @NotNull Vec3 b, @NotNull Vec3 c) {
+        var i = a.minus(o).unit();
+        var j = b.minus(o).unit();
+        var k = c.minus(o).unit();
+
+        return 2 * Math.atan(Math.abs(i.times(j.cross(k))) / (1 + i.times(j) + j.times(k) + k.times(i)));
+    }
+}