minor performance improvements

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);
@@ -50,14 +51,13 @@ public class Examples {
     public static @NotNull Example getSimpleScene(int height) {
         if (height <= 0) height = 675;
         return new Example(
-                new Scene(
-                        getSkyBox(),
+                new Scene(getSkyBox(), List.of(
                         new Sphere(new Vec3(0, -100.5, -1.0), 100.0, new LambertianMaterial(new Color(0.8, 0.8, 0.0))),
                         new Sphere(new Vec3(0, 0, -1.2), 0.5, new LambertianMaterial(new Color(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(new Color(0.8, 0.6, 0.2), 0.0))
-                ),
+                )),
                 SimpleCamera.builder()
                         .withImage(height * 16 / 9, height)
                         .build()
@@ -118,14 +118,13 @@ public class Examples {
     public static @NotNull Example getSquares(int height) {
         if (height <= 0) height = 600;
         return new Example(
-                new Scene(
-                        getSkyBox(),
+                new Scene(getSkyBox(), List.of(
                         new Parallelogram(new Vec3(-3, -2, 5), new Vec3(0, 0, -4), new Vec3(0, 4, 0), new LambertianMaterial(new Color(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(new Color(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(new Color(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(new Color(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(new Color(0.2, 0.8, 0.8)))
-                ),
+                )),
                 SimpleCamera.builder()
                         .withImage(height, height)
                         .withFieldOfView(Math.toRadians(80))
@@ -138,12 +137,12 @@ public class Examples {
     public static @NotNull Example getLight(int height) {
         if (height <= 0) height = 225;
         return new Example(
-                new Scene(
+                new Scene(List.of(
                         new Sphere(new Vec3(0, -1000, 0), 1000, new LambertianMaterial(new Color(0.2, 0.2, 0.2))),
                         new Sphere(new Vec3(0, 2, 0), 2, new LambertianMaterial(new Color(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(new Color(4.0, 4.0, 4.0))),
                         new Sphere(new Vec3(0, 7, 0), 2, new DiffuseLight(new Color(4.0, 4.0, 4.0)))
-                ),
+                )),
                 SimpleCamera.builder()
                         .withImage(height * 16 / 9, height)
                         .withFieldOfView(Math.toRadians(20))
@@ -162,15 +161,16 @@ public class Examples {
         var light = new DiffuseLight(new Color(15.0, 15.0, 15.0));
 
         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 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))
-                ),
+                        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))
@@ -188,21 +188,49 @@ public class Examples {
         var light = new DiffuseLight(new Color(7.0, 7.0, 7.0));
 
         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 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)),
+                                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(Color.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)),
+                                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(Color.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 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 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))
@@ -264,10 +292,9 @@ public class Examples {
         if (height <= 0) height = 450;
 
         return new Example(
-                new Scene(
-                        getSkyBox(),
+                new Scene(getSkyBox(), List.of(
                         new Sphere(Vec3.ZERO, 2, new LambertianMaterial(new ImageTexture("/earthmap.jpg")))
-                ),
+                )),
                 SimpleCamera.builder()
                         .withImage(height * 16 / 9, height)
                         .withFieldOfView(Math.toRadians(20))
@@ -283,11 +310,10 @@ public class Examples {
         var material = new LambertianMaterial(new PerlinTexture(4));
 
         return new Example(
-                new Scene(
-                        getSkyBox(),
+                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))

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

@@ -13,6 +13,8 @@ import java.nio.file.Path;
 import java.util.function.IntFunction;
 
 public class Main {
+    public static final boolean DEBUG = false;
+
     public static void main(String[] args) throws IOException {
         var config = Config.parse(args);
         var example = config.example;

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

@@ -35,26 +35,49 @@ public record AABB(@NotNull Vec3 min, @NotNull Vec3 max) {
         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 origin = ray.origin();
         var direction = ray.direction();
@@ -85,6 +108,13 @@ public record AABB(@NotNull Vec3 min, @NotNull Vec3 max) {
         return tlmax < tumin && tumin >= range.min() && tlmax <= range.max();
     }
 
+    /**
+     * 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(),

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

@@ -11,16 +11,6 @@ public record Ray(@NotNull Vec3 origin, @NotNull Vec3 direction) {
     }
 
     public @NotNull Vec3 at(double t) {
-        return new Vec3(
-                Math.fma(t, direction.x(), origin.x()),
-                Math.fma(t, direction.y(), origin.y()),
-                Math.fma(t, direction.z(), origin.z())
-        );
-    }
-
-    public int vmask() {
-        return (direction().x() < 0 ? 1 : 0)
-                | (direction().y() < 0 ? 2 : 0)
-                | (direction().z() < 0 ? 4 : 0);
+        return Vec3.fma(t, direction, origin);
     }
 }

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

@@ -5,6 +5,8 @@ import org.jetbrains.annotations.NotNull;
 import java.util.Optional;
 import java.util.random.RandomGenerator;
 
+import static eu.jonahbauer.raytracing.Main.DEBUG;
+
 public record Vec3(double x, double y, double z) {
     public static final Vec3 ZERO = new Vec3(0, 0, 0);
     public static final Vec3 MAX = new Vec3(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE);
@@ -14,43 +16,63 @@ public record Vec3(double x, double y, double z) {
     public static final Vec3 UNIT_Z = new Vec3(0, 0, 1);
 
     public Vec3 {
-        assert Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(z) : "x, y and z must be finite";
+        if (DEBUG) {
+            if (!Double.isFinite(x) || !Double.isFinite(y) || !Double.isFinite(z)) {
+                throw new IllegalArgumentException("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 new Vec3(
-                Math.fma(2, random.nextDouble(), -1),
-                Math.fma(2, random.nextDouble(), -1),
-                Math.fma(2, random.nextDouble(), -1)
-        );
     }
 
     /**
      * {@return a uniformly random unit vector}
      */
-    public static @NotNull Vec3 random(@NotNull RandomGenerator random, boolean unit) {
-        if (!unit) return random(random);
-        Vec3 vec;
+    public static @NotNull Vec3 random(@NotNull RandomGenerator random) {
+        double x, y, z;
         double squared;
         do {
-            vec = random(random);
-            squared = vec.squared();
+            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);
-        return vec.div(Math.sqrt(squared));
+        var factor = 1 / Math.sqrt(squared);
+        return new Vec3(x * factor, y * factor, z * factor);
     }
 
+    /**
+     * {@return a uniformly random unit vector on the opposite hemisphere of the given <code>direction</code>}
+     */
+    public static @NotNull Vec3 randomOppositeHemisphere(@NotNull RandomGenerator random, @NotNull Vec3 direction) {
+        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 || 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) {
         var factor = - 2 * normal.times(vec);
-        return new Vec3(
-                Math.fma(factor, normal.x(), vec.x()),
-                Math.fma(factor, normal.y(), vec.y()),
-                Math.fma(factor, normal.z(), vec.z())
-        );
+        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);
@@ -62,16 +84,35 @@ 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)));
     }
 
+    /**
+     * {@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();
     }
 
+    /**
+     * 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(
@@ -81,6 +122,11 @@ public record Vec3(double x, double y, double 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()),
@@ -89,6 +135,11 @@ public record Vec3(double x, double y, double z) {
         );
     }
 
+    /**
+     * {@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()),
@@ -97,6 +148,24 @@ public record Vec3(double x, double y, double 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 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;
+    }
+
     public @NotNull Vec3 plus(double x, double y, double z) {
         return new Vec3(this.x + x, this.y + y, this.z + z);
     }
@@ -105,59 +174,115 @@ public record Vec3(double x, double y, double z) {
         return new Vec3(this.x - x, this.y - y, this.z - z);
     }
 
-    public @NotNull Vec3 plus(@NotNull Vec3 b) {
-        return new Vec3(this.x + b.x, this.y + b.y, this.z + b.z);
+    /**
+     * Adds a vector to this vector
+     * @param other a vector
+     * @return the sum of this and the other vector
+     */
+    public @NotNull Vec3 plus(@NotNull Vec3 other) {
+        return new Vec3(this.x + other.x, this.y + other.y, this.z + other.z);
     }
 
-    public @NotNull Vec3 minus(@NotNull Vec3 b) {
-        return new Vec3(this.x - b.x, this.y - b.y, this.z - b.z);
+    /**
+     * Subtracts a vector from this vector
+     * @param other a vector
+     * @return the difference of this and the other vector
+     */
+    public @NotNull Vec3 minus(@NotNull Vec3 other) {
+        return new Vec3(this.x - other.x, this.y - other.y, this.z - other.z);
     }
 
-    public double times(@NotNull Vec3 b) {
-        return this.x * b.x + this.y * b.y + this.z * b.z;
+    /**
+     * Computes the scalar product of this and another vector
+     * @param other a vector
+     * @return the scalar product
+     */
+    public double times(@NotNull Vec3 other) {
+        return this.x * other.x + this.y * other.y + this.z * other.z;
     }
 
-    public @NotNull Vec3 times(double b) {
-        return new Vec3(this.x * b, this.y * b, this.z * b);
+    /**
+     * Multiplies this vector with a scalar
+     * @param t a scalar
+     * @return the product of this vector and the scalar
+     */
+    public @NotNull Vec3 times(double t) {
+        return new Vec3(this.x * t, this.y * t, this.z * t);
     }
 
+    /**
+     * Negates this vector.
+     * {@return the negated vector}
+     */
     public @NotNull Vec3 neg() {
         return new Vec3(-x, -y, -z);
     }
 
+    /**
+     * Inverts each component of this vector.
+     * @return the inverted vector.
+     */
     public @NotNull Vec3 inv() {
         return new Vec3(1 / x, 1 / y, 1 / z);
     }
 
-    public @NotNull Vec3 cross(@NotNull Vec3 b) {
+    /**
+     * Computes the cross-product of this and another vector
+     * @param other a vector
+     * @return the cross-product
+     */
+    public @NotNull Vec3 cross(@NotNull Vec3 other) {
         return new Vec3(
-                this.y() * b.z() - b.y() * this.z(),
-                this.z() * b.x() - b.z() * this.x(),
-                this.x() * b.y() - b.x() * this.y()
+                Math.fma(this.y, other.z, - other.y * this.z),
+                Math.fma(this.z, other.x, - other.z * this.x),
+                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);
+    /**
+     * Divides this vector by a scalar
+     * @param t a scalar
+     * @return this vector divided by the scalar
+     */
+    public @NotNull Vec3 div(double t) {
+        return new Vec3(this.x / t, this.y / t, this.z / t);
     }
 
+    /**
+     * {@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));
     }
 
+    /**
+     * {@return the n-th component of this vector}
+     * @param axis the component index
+     */
     public double get(int axis) {
         return switch (axis) {
             case 0 -> x;
@@ -178,4 +303,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/camera/SimpleCamera.java

@@ -60,7 +60,7 @@ public final class SimpleCamera implements 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));
     }
 
     /**

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

@@ -25,7 +25,7 @@ public record MetallicMaterial(@NotNull Texture texture, double fuzz) implements
     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 = newDirection.unit().plus(Vec3.random(random, true).times(fuzz));
+            newDirection = Vec3.fma(fuzz, Vec3.random(random), newDirection.unit());
         }
         var attenuation = texture.get(hit);
         return Optional.of(new SpecularScatterResult(attenuation, new Ray(hit.position(), newDirection)));

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

@@ -3,7 +3,7 @@ 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.HittableProbabilityDensityFunction;
+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;
@@ -16,6 +16,8 @@ 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;
@@ -62,7 +64,7 @@ public final class SimpleRenderer implements Renderer {
      * 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();
+        var random = new Random(0);
 
         // render one sample after the other
         int i = 0;
@@ -94,7 +96,7 @@ public final class SimpleRenderer implements Renderer {
      * 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();
+        var splittable = new SplittableRandom(0);
         // render one pixel after the other
         getScanlineStream(camera.getHeight(), parallel).forEach(y -> {
             var random = splittable.split();
@@ -104,6 +106,9 @@ public final class SimpleRenderer implements Renderer {
                 for (int sj = 0; sj < sqrtSamplesPerPixel; sj++) {
                     for (int si = 0; si < sqrtSamplesPerPixel; si++) {
                         var ray = camera.cast(x, y, si, sj, sqrtSamplesPerPixel, random);
+                        if (DEBUG) {
+                            System.out.println("Casting ray " + ray + " through pixel (" + x + "," + y + ") at subpixel (" + si + "," + sj + ")...");
+                        }
                         var c = getColor(scene, ray, random);
                         color = Color.average(color, c, ++i);
                     }
@@ -125,34 +130,83 @@ public final class SimpleRenderer implements Renderer {
         var attenuation = Color.WHITE;
 
         while (depth-- > 0) {
-            var optional = scene.hit(ray, new Range(0.001, Double.POSITIVE_INFINITY));
+            var optional = scene.hit(ray);
             if (optional.isEmpty()) {
-                color = Color.add(color, Color.multiply(attenuation, scene.getBackgroundColor(ray)));
+                var background = scene.getBackgroundColor(ray);
+                color = Color.add(color, Color.multiply(attenuation, background));
+                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);
+            if (DEBUG && !Color.BLACK.equals(emitted)) {
+                System.out.println("  Emitted: " + emitted);
+            }
+
             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;
+            if (result.isEmpty()) {
+                if (DEBUG) {
+                    System.out.println("  Absorbed");
+                }
+                break;
+            }
 
             switch (result.get()) {
                 case Material.SpecularScatterResult(var a, var scattered) -> {
                     attenuation = Color.multiply(attenuation, a);
                     ray = scattered;
+
+                    if (DEBUG) {
+                        System.out.println("  Specular scattering with albedo " + a);
+                    }
                 }
                 case Material.PdfScatterResult(var a, var pdf) -> {
+                    if (scene.getTargets() == null) {
                         attenuation = Color.multiply(attenuation, a);
                         ray = new Ray(hit.position(), pdf.generate(random));
+
+                        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 = Color.multiply(attenuation, Color.multiply(a, factor));
+                        ray = new Ray(hit.position(), direction);
+
+                        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;
     }

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

@@ -15,13 +15,13 @@ public record CosineProbabilityDensityFunction(@NotNull Vec3 normal) implements
 
     @Override
     public double value(@NotNull Vec3 direction) {
-        var cos = normal.times(direction.unit());
+        var cos = normal.times(direction);
         return Math.max(0, cos / Math.PI);
     }
 
     @Override
     public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
-        var out = normal().plus(Vec3.random(random, true));
+        var out = normal().plus(Vec3.random(random));
         return out.isNearZero() ? normal() : out;
     }
 }

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

@@ -6,6 +6,13 @@ 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,
@@ -23,7 +30,9 @@ public record MixtureProbabilityDensityFunction(
 
     @Override
     public double value(@NotNull Vec3 direction) {
-        return weight * a.value(direction) + (1 - weight) * b.value(direction);
+        var v = a.value(direction);
+        var w = b.value(direction);
+        return Math.fma(weight, v, Math.fma(-weight, w, w));
     }
 
     @Override

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

@@ -5,7 +5,21 @@ 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

@@ -5,6 +5,9 @@ import org.jetbrains.annotations.NotNull;
 
 import java.util.random.RandomGenerator;
 
+/**
+ * A probability density function sampling the sphere uniformly.
+ */
 public record SphereProbabilityDensityFunction() implements ProbabilityDensityFunction {
 
     @Override
@@ -14,6 +17,6 @@ public record SphereProbabilityDensityFunction() implements ProbabilityDensityFu
 
     @Override
     public @NotNull Vec3 generate(@NotNull RandomGenerator random) {
-        return Vec3.random(random, true);
+        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/texture/CheckerTexture.java

@@ -5,7 +5,6 @@ import org.jetbrains.annotations.NotNull;
 
 public record CheckerTexture(double scale, @NotNull Texture even, @NotNull Texture odd) implements Texture {
 
-
     @Override
     public @NotNull Color get(double u, double v, @NotNull Vec3 p) {
         var x = (int) Math.floor(p.x() / scale);

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

@@ -1,11 +1,15 @@
 package eu.jonahbauer.raytracing.render.texture;
 
+import eu.jonahbauer.raytracing.math.Ray;
 import eu.jonahbauer.raytracing.math.Vec3;
+import eu.jonahbauer.raytracing.scene.SkyBox;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.Random;
 
-public record Color(double r, double g, double b) implements Texture {
+import static eu.jonahbauer.raytracing.Main.DEBUG;
+
+public record Color(double r, double g, double b) implements Texture, SkyBox {
     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 RED = new Color(1.0, 0.0, 0.0);
@@ -26,6 +30,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());
     }
@@ -78,6 +86,17 @@ public record Color(double r, double g, double b) implements Texture {
         this(red / 255f, green / 255f, blue / 255f);
     }
 
+    public Color {
+        if (DEBUG) {
+            if (!Double.isFinite(r) || !Double.isFinite(g) || !Double.isFinite(b)) {
+                throw new IllegalArgumentException("r, g and b must be finite");
+            }
+            if (r < 0 || g < 0 || b < 0) {
+                throw new IllegalArgumentException("r, g and b must be non-negative");
+            }
+        }
+    }
+
     public int red() {
         return toInt(r);
     }
@@ -95,6 +114,11 @@ public record Color(double r, double g, double b) implements Texture {
         return this;
     }
 
+    @Override
+    public @NotNull Color getColor(@NotNull Ray ray) {
+        return this;
+    }
+
     @Override
     public boolean isUVRequired() {
         return false;

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

@@ -53,7 +53,7 @@ public final class PerlinTexture implements Texture {
         this.mask = count - 1;
         this.randvec = new Vec3[count];
         for (int i = 0; i < count; i++) {
-            this.randvec[i] = Vec3.random(random, true);
+            this.randvec[i] = Vec3.random(random);
         }
         this.permX = generatePerm(count, random);
         this.permY = generatePerm(count, random);

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

@@ -5,12 +5,27 @@ 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 Color 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 Color 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,13 +1,27 @@
 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.render.material.Material;
+import eu.jonahbauer.raytracing.render.texture.Texture;
 import org.jetbrains.annotations.NotNull;
 
 import java.util.Objects;
 
+/**
+ * The result of a {@linkplain Hittable#hit(Ray, Range) 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, @NotNull Vec3 position, @NotNull Vec3 normal,
+        double t, @NotNull Vec3 position, @NotNull Vec3 normal, @NotNull Hittable target,
         @NotNull Material material, double u, double v, boolean isFrontFace
 ) implements Comparable<HitResult> {
     public HitResult {
@@ -16,7 +30,7 @@ public record HitResult(
     }
 
     public @NotNull HitResult withPositionAndNormal(@NotNull Vec3 position, @NotNull Vec3 normal) {
-        return new HitResult(t, position, normal, material, u, v, isFrontFace);
+        return new HitResult(t, position, normal, target, material, u, v, isFrontFace);
     }
 
     @Override

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

@@ -11,21 +11,38 @@ import org.jetbrains.annotations.NotNull;
 import java.util.Optional;
 
 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
-     * the ray origin, or <code>Double.NaN</code> if the ray does not intersect this shape}
+     * @see #hit(Ray, Range)
+     */
+    default @NotNull Optional<HitResult> hit(@NotNull Ray ray) {
+        return hit(ray, FORWARD);
+    }
+
+    /**
+     * 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);
 
+    /**
+     * {@return the axis-aligned bounding box of this hittable}
+     */
     @NotNull AABB getBoundingBox();
 
     default @NotNull Hittable translate(@NotNull Vec3 offset) {
-        return new Translate(this, offset);
+        return Translate.create(this, offset);
     }
 
     default @NotNull Hittable rotateY(double angle) {
-        return new RotateY(this, angle);
+        return RotateY.create(this, angle);
     }
 }

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

@@ -6,7 +6,9 @@ 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.ArrayList;
 import java.util.List;
 import java.util.Objects;
 
@@ -14,29 +16,28 @@ public final class Scene extends HittableCollection {
     private final @NotNull HittableCollection objects;
     private final @NotNull SkyBox background;
 
+    private final @Nullable List<@NotNull Target> targets;
+
     public Scene(@NotNull List<? extends @NotNull Hittable> objects) {
-        this(Color.BLACK, objects);
+        this(objects, null);
     }
 
-    public Scene(@NotNull Color background, @NotNull List<? extends @NotNull Hittable> objects) {
-        this(SkyBox.solid(background), objects);
+    public Scene(@NotNull List<? extends @NotNull Hittable> objects, @Nullable List<? extends @NotNull Target> targets) {
+        this(Color.BLACK, objects, targets);
     }
 
     public Scene(@NotNull SkyBox background, @NotNull List<? extends @NotNull Hittable> objects) {
-        this.objects = new HittableBinaryTree(objects);
+        this(background, objects, null);
+    }
+
+    public Scene(@NotNull SkyBox background, @NotNull List<? extends @NotNull Hittable> objects, @Nullable List<? extends @NotNull Target> targets) {
+        var list = new ArrayList<Hittable>(objects.size() + (targets != null ? targets.size() : 0));
+        list.addAll(objects);
+        if (targets != null) list.addAll(targets);
+        this.objects = new HittableBinaryTree(list);
         this.background = Objects.requireNonNull(background);
-    }
 
-    public Scene(@NotNull Hittable @NotNull... objects) {
-        this(List.of(objects));
-    }
-
-    public Scene(@NotNull Color background, @NotNull Hittable @NotNull... objects) {
-        this(background, List.of(objects));
-    }
-
-    public Scene(@NotNull SkyBox background, @NotNull Hittable @NotNull... objects) {
-        this(background, List.of(objects));
+        this.targets = targets != null ? List.copyOf(targets) : null;
     }
 
     @Override
@@ -49,6 +50,10 @@ public final class Scene extends HittableCollection {
         return objects.getBoundingBox();
     }
 
+    public @Nullable List<@NotNull Target> getTargets() {
+        return targets;
+    }
+
     public @NotNull Color getBackgroundColor(@NotNull Ray ray) {
         return background.getColor(ray);
     }

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

@@ -18,8 +18,4 @@ public interface SkyBox {
             return Color.lerp(bottom, top, alt / Math.PI + 0.5);
         };
     }
-
-    static @NotNull SkyBox solid(@NotNull Color color) {
-        return _ -> color;
-    }
 }

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

@@ -52,10 +52,40 @@ public abstract class Hittable2D implements Hittable {
 
         var frontFace = denominator < 0;
         return Optional.of(new HitResult(
-                t, position, frontFace ? normal : normal.neg(),
+                t, position, frontFace ? normal : normal.neg(), this,
                 material, alpha, beta, frontFace
         ));
     }
 
+    protected double hit0(@NotNull Ray ray, @NotNull Range range) {
+        var denominator = ray.direction().times(normal);
+        if (Math.abs(denominator) < 1e-8) return Double.NaN; // parallel
+
+        var t = (d - ray.origin().times(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

@@ -1,11 +1,16 @@
 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;
 
-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 +27,24 @@ public final class Parallelogram extends Hittable2D {
     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/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;
 
@@ -25,6 +28,15 @@ public final class Box implements Hittable {
         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,
@@ -107,7 +119,7 @@ public final class Box implements Hittable {
         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, material, u, v, frontFace));
+        return Optional.of(new HitResult(t, position, normal, this, material, u, v, frontFace));
     }
 
     @Override
@@ -115,6 +127,66 @@ 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)).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()),
@@ -155,5 +227,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).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

@@ -32,7 +32,7 @@ public record ConstantMedium(@NotNull Hittable boundary, double density, @NotNul
         if (hitDistance > distance) return Optional.empty();
 
         var t = tmin + hitDistance / length;
-        return Optional.of(new HitResult(t, ray.at(t), Vec3.UNIT_X, material, 0, 0, true)); // arbitrary normal, u, v and isFrontFace
+        return Optional.of(new HitResult(t, ray.at(t), Vec3.UNIT_X, this, material, 0, 0, true)); // arbitrary normal, u, v and isFrontFace
     }
 
     @Override

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

@@ -7,24 +7,32 @@ 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;
 
     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)
@@ -33,23 +41,11 @@ public final class Sphere implements Hittable {
 
     @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 t = hit0(ray, range);
+        if (Double.isNaN(t)) return Optional.empty();
 
         var position = ray.at(t);
-        var normal = position.minus(center).div(radius);
+        var normal = Vec3.fma(invRadius, position, normalizedCenter);
         var frontFace = normal.times(ray.direction()) < 0;
 
         double u;
@@ -65,13 +61,57 @@ public final class Sphere implements Hittable {
         }
 
         return Optional.of(new HitResult(
-                t, position, frontFace ? normal : normal.neg(),
+                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().times(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

@@ -5,15 +5,26 @@ 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;
 
-public final class RotateY extends Transform {
+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 RotateY(@NotNull Hittable object, double angle) {
+    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);
@@ -45,45 +56,58 @@ public final class RotateY extends Transform {
     }
 
     @Override
-    protected @NotNull Ray transform(@NotNull Ray ray) {
+    protected final @NotNull Ray transform(@NotNull Ray ray) {
         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) {
+    protected final @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);
     }
 
+    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

@@ -2,12 +2,15 @@ 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 {
     protected final @NotNull Hittable object;

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

@@ -5,13 +5,24 @@ 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;
 
-public final class Translate extends Transform {
-    private final @NotNull Vec3 offset;
+import java.util.random.RandomGenerator;
+
+public sealed class Translate extends Transform {
+    protected final @NotNull Vec3 offset;
     private final @NotNull AABB bbox;
 
-    public Translate(@NotNull Hittable object, @NotNull Vec3 offset) {
+    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;
 
@@ -23,17 +34,40 @@ public final class Translate extends Transform {
     }
 
     @Override
-    protected @NotNull Ray transform(@NotNull Ray ray) {
+    protected final @NotNull Ray transform(@NotNull Ray ray) {
         return new Ray(ray.origin().minus(offset), ray.direction());
     }
 
     @Override
-    protected @NotNull HitResult transform(@NotNull HitResult result) {
+    protected final @NotNull HitResult transform(@NotNull HitResult result) {
         return result.withPositionAndNormal(result.position().plus(offset), result.normal());
     }
 
     @Override
-    public @NotNull AABB getBoundingBox() {
+    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/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.times(b) + b.times(c) + c.times(a)));
+        return angle < 0 ? 2 * Math.PI + angle : angle;
+    }
+}