Fancy benchmarks

benches/Cargo.toml

@@ -0,0 +1,10 @@
+[package]
+name = "benches"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+bvh = { path = ".." }
+clap = { version = "4.5.27", features = ["derive"] }
+nalgebra = "0.33.2"
+rand = "0.9.0"

benches/plot.py

@@ -0,0 +1,75 @@
+import matplotlib.pyplot as plt
+from matplotlib.colors import LinearSegmentedColormap
+import numpy as np
+import math
+from subprocess import Popen, PIPE
+
+# number of benchmarks per power of 10
+resolution = 4
+# benchmark 1..=10^max rays/triangles
+max = 4
+
+all_benchmarks = [
+    "traverse",
+    "traverse_iterator",
+    "nearest_traverse_iterator",
+    "nearest_child_traverse_iterator"
+]
+
+cmap = LinearSegmentedColormap.from_list("speedup", ["red", "white", "green"], N=256)
+
+def bench(rays, triangles, samples, benchmark):
+    proc = Popen([
+        "cargo",
+        "run",
+        "--release",
+        "--",
+        "--triangles",
+        str(triangles),
+        "--rays",
+        str(rays),
+        "--samples",
+        str(samples),
+        "--benchmark",
+        benchmark
+    ], stdout=PIPE)
+    stdout, _ = proc.communicate()
+    return float(stdout)
+
+def plot(colormaps):
+    mesh = []
+    for x in range(0, max * resolution + 1):
+        mesh.append(round(math.pow(10, x * (1.0 / resolution))))
+
+    n = len(all_benchmarks)
+    fig, axs = plt.subplots(1, n, figsize=(n * 2 + 2, 3),
+                            layout='constrained', squeeze=False)
+
+    for (ax, benchmark) in zip(axs.flat, all_benchmarks):
+        data = np.empty((max * resolution + 1, max * resolution + 1))
+
+        with np.nditer(data, flags=['multi_index'], op_flags=['writeonly']) as it:
+            for x in it:
+                rays = mesh[it.multi_index[0]]
+                triangles = mesh[it.multi_index[1]]
+                samples = 51
+                if triangles < 50:
+                    # sample more due to counteract observed noise
+                    samples = 501
+                elif triangles >= 100 and rays >= 100:
+                    # these run longer so less intrinsic noise; save time
+                    # by sampling less
+                    samples = 3
+                x[...] = bench(rays, triangles, samples, benchmark)
+
+        ax.set_title(benchmark + "\nBVH Speedup Coefficient")
+        ax.set_xlabel("Triangles")
+        ax.set_ylabel("Rays")
+        ax.set_xscale('log')
+        ax.set_yscale('log')
+        psm = ax.pcolormesh(mesh, mesh, data, cmap=cmap, rasterized=True, vmin=0, vmax=2, shading="gouraud")
+        fig.colorbar(psm, ax=ax)
+
+    plt.show()
+
+plot([None])

benches/src/main.rs

@@ -0,0 +1,220 @@
+use std::{hint::black_box, time::Instant};
+
+use bvh::{
+    aabb::{Aabb, Bounded},
+    bounding_hierarchy::BHShape,
+    bvh::Bvh,
+    ray::Ray,
+};
+use clap::{Parser, ValueEnum};
+use nalgebra::{Point3, Vector3};
+use rand::{rng, Rng};
+
+#[derive(Parser)]
+#[command(version, about, long_about = None)]
+struct Cli {
+    #[arg(long)]
+    rays: usize,
+    #[arg(long)]
+    triangles: usize,
+    #[arg(long)]
+    samples: usize,
+    #[arg(long)]
+    benchmark: Benchmark,
+}
+
+#[derive(ValueEnum, Debug, Clone)]
+#[clap(rename_all = "snake_case")]
+enum Benchmark {
+    Traverse,
+    TraverseIterator,
+    NearestTraverseIterator,
+    NearestChildTraverseIterator,
+}
+
+fn main() {
+    let cli = Cli::parse();
+    let mut rng = rng();
+
+    let mut samples = Vec::new();
+    let mut rays = Vec::new();
+    let mut triangles = Vec::new();
+
+    for i in 0..cli.samples {
+        rays.clear();
+        triangles.clear();
+
+        for _ in 0..cli.rays {
+            rays.push(Ray::<f32, 3>::new(
+                Point3::new(
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                ),
+                Vector3::new(
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                ),
+            ));
+        }
+
+        for _ in 0..cli.triangles {
+            let center = Point3::new(
+                rng.random_range(-1000.0..=1000.0),
+                rng.random_range(-1000.0..=1000.0),
+                rng.random_range(-1000.0..=1000.0),
+            );
+            let a = center
+                + Vector3::new(
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                );
+            let b = center
+                + Vector3::new(
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                );
+            let c = center
+                + Vector3::new(
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                    rng.random_range(-1.0..=1.0),
+                );
+            triangles.push(Triangle {
+                a,
+                b,
+                c,
+                aabb: Aabb::empty().grow(&a).grow(&b).grow(&c),
+                bh_node_index: usize::MAX,
+            });
+        }
+
+        let mut brute_force_duration = f64::NAN;
+        let mut bvh_duration = f64::NAN;
+
+        let mut measure_brute_force = |triangles: &[Triangle]| {
+            let start_brute_force = Instant::now();
+            for ray in &rays {
+                match cli.benchmark {
+                    Benchmark::Traverse | Benchmark::TraverseIterator => {
+                        black_box(
+                            black_box(&triangles)
+                                .iter()
+                                .filter(|triangle| triangle.intersect(&ray))
+                                .count(),
+                        );
+                    }
+                    Benchmark::NearestTraverseIterator
+                    | Benchmark::NearestChildTraverseIterator => {
+                        let mut results = black_box(
+                            black_box(&triangles)
+                                .iter()
+                                .filter(|triangle| triangle.intersect(&ray))
+                                .collect::<Vec<_>>(),
+                        );
+                        results.sort_by_cached_key(|triangle| {
+                            (ray.intersection_slice_for_aabb(&triangle.aabb).unwrap().0 * 1000.0)
+                                as usize
+                        });
+                        black_box(results);
+                    }
+                }
+            }
+            brute_force_duration = start_brute_force.elapsed().as_secs_f64();
+        };
+
+        let mut measure_bvh = |triangles: &mut Vec<Triangle>| {
+            let start_bvh = Instant::now();
+            let bvh = Bvh::build(black_box(triangles));
+            for ray in &rays {
+                match cli.benchmark {
+                    Benchmark::Traverse => {
+                        black_box(
+                            bvh.traverse(black_box(ray), black_box(&triangles))
+                                .into_iter()
+                                .filter(|triangle| triangle.intersect(&ray))
+                                .count(),
+                        );
+                    }
+                    Benchmark::TraverseIterator => {
+                        black_box(
+                            bvh.traverse_iterator(black_box(ray), black_box(&triangles))
+                                .filter(|triangle| triangle.intersect(&ray))
+                                .count(),
+                        );
+                    }
+                    Benchmark::NearestTraverseIterator => {
+                        black_box(
+                            bvh.nearest_traverse_iterator(black_box(ray), black_box(&triangles))
+                                .filter(|triangle| triangle.intersect(&ray))
+                                .count(),
+                        );
+                    }
+                    Benchmark::NearestChildTraverseIterator => {
+                        black_box(
+                            bvh.nearest_child_traverse_iterator(
+                                black_box(ray),
+                                black_box(&triangles),
+                            )
+                            .filter(|triangle| triangle.intersect(&ray))
+                            .count(),
+                        );
+                    }
+                }
+            }
+            bvh_duration = start_bvh.elapsed().as_secs_f64();
+        };
+
+        // Flip order to minimize bias due to caching.
+        if i % 2 == 0 {
+            measure_bvh(&mut triangles);
+            measure_brute_force(&triangles);
+        } else {
+            measure_brute_force(&triangles);
+            measure_bvh(&mut triangles);
+        }
+
+        let bvh_speedup = brute_force_duration / bvh_duration;
+        samples.push(bvh_speedup);
+    }
+
+    samples.sort_by(|a, b| a.partial_cmp(b).unwrap());
+
+    // Median.
+    println!("{}", samples[cli.samples / 2]);
+}
+
+struct Triangle {
+    a: Point3<f32>,
+    b: Point3<f32>,
+    c: Point3<f32>,
+    aabb: Aabb<f32, 3>,
+    bh_node_index: usize,
+}
+
+impl Triangle {
+    fn intersect(&self, ray: &Ray<f32, 3>) -> bool {
+        ray.intersects_triangle(&self.a, &self.b, &self.c)
+            .distance
+            .is_finite()
+    }
+}
+
+impl Bounded<f32, 3> for Triangle {
+    fn aabb(&self) -> Aabb<f32, 3> {
+        self.aabb
+    }
+}
+
+impl BHShape<f32, 3> for Triangle {
+    fn bh_node_index(&self) -> usize {
+        self.bh_node_index
+    }
+
+    fn set_bh_node_index(&mut self, val: usize) {
+        self.bh_node_index = val;
+    }
+}

justfile

@@ -33,6 +33,10 @@ bench:
 bench_simd:
     cargo bench --features bench,simd
 
+# WIP/temporary
+bench_fancy:
+    cd benches; python3 plot.py
+
 # fuzz the library
 fuzz:
     cargo fuzz run fuzz