diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..85ee419 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,113 @@ +name: CI + +on: + push: + branches: [ main ] + pull_request: + branches: [ main ] + +jobs: + test: + name: Test + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, windows-latest, macos-latest] + zig-version: [master] + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Setup Zig + uses: mlugg/setup-zig@v2 + with: + version: ${{ matrix.zig-version }} + + - name: Run tests + run: zig build test + + benchmark: + name: Benchmark + runs-on: ubuntu-latest + if: github.event_name == 'pull_request' + permissions: + pull-requests: write + contents: read + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Setup Zig + uses: mlugg/setup-zig@v2 + with: + version: master + + - name: Run benchmark + run: | + echo "Running benchmarks..." + zig build bench -Doptimize=ReleaseFast > benchmark_results.txt 2>&1 || true + echo "Benchmark completed" + continue-on-error: true + + - name: Read benchmark results + id: benchmark + run: | + if [ -f benchmark_results.txt ]; then + echo "results<> $GITHUB_OUTPUT + cat benchmark_results.txt >> $GITHUB_OUTPUT + echo "EOF" >> $GITHUB_OUTPUT + else + echo "results=No benchmark results found" >> $GITHUB_OUTPUT + fi + + - name: Comment PR with benchmark results + uses: actions/github-script@v7 + with: + script: | + const results = `${{ steps.benchmark.outputs.results }}`; + + const comment = `## ๐Ÿš€ Benchmark Results + +
+ Click to view benchmark results + + \`\`\` + ${results} + \`\`\` + +
+ + > Benchmarks run on: \`ubuntu-latest\` with \`ReleaseFast\` optimization + > + > To run benchmarks locally: \`zig build bench -Doptimize=ReleaseFast\` + `; + + // Check if we already commented on this PR + const { data: comments } = await github.rest.issues.listComments({ + owner: context.repo.owner, + repo: context.repo.repo, + issue_number: context.issue.number, + }); + + const botComment = comments.find(comment => + comment.user.type === 'Bot' && + comment.body.includes('๐Ÿš€ Benchmark Results') + ); + + if (botComment) { + // Update existing comment + await github.rest.issues.updateComment({ + owner: context.repo.owner, + repo: context.repo.repo, + comment_id: botComment.id, + body: comment + }); + } else { + // Create new comment + await github.rest.issues.createComment({ + owner: context.repo.owner, + repo: context.repo.repo, + issue_number: context.issue.number, + body: comment + }); + } diff --git a/.gitignore b/.gitignore index cbcffab..28e7464 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .zig-cache/ zig-out/ +zig-pkg/ .vscode/ ignore_* \ No newline at end of file diff --git a/README.md b/README.md index 555d840..4a1603f 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,101 @@ -This library is a simple linear algebra library. +# Zig Math Library (ZML) -This library was written for my own use, and will be extended and optimized as I need it. +[![CI](https://github.com/flying-swallow/zig-linear-algebra/actions/workflows/ci.yml/badge.svg)](https://github.com/flying-swallow/zig-linear-algebra/actions/workflows/ci.yml) +A high-performance linear algebra library for Zig, providing vector, matrix, quaternion, and geometry operations. +## Features + +- **Vectors** (`zml.vec`) โ€” norm, dot, cross, normalize, reflect, distance, angle, swizzle, fused `sin_cos`, and more. +- **Matrices** (`zml.Mat`) โ€” generic column-major `Mat(T, cols, rows)` with multiply, transforms, etc. +- **Quaternions** (`zml.quat`) โ€” rotation, slerp/nlerp, axis-angle and Euler conversions. +- **Geometry** (`zml.geom`) โ€” AABB, sphere, plane, capsule, OBB, ray, frustum, overlap/containment tests. +- Extras: `zml.scalar` (clamp/lerp/smoothstep), `zml.color`, `zml.packing`, `zml.random`. +- Built on Zig's native `@Vector` SIMD types โ€” **but functions also accept plain arrays** (see below). + +## Installation + +```zig +zig fetch --save "git+https://github.com/flying-swallow/zig-linear-algebra.git" +``` + +Then in your `build.zig`: + +```zig +const zml = b.dependency("zml", .{ + .target = target, + .optimize = optimize, +}); + +exe.root_module.addImport("zml", zml.module("zml")); +``` + +## Usage + +```zig +const zml = @import("zml"); + +const a = zml.Vec3f32{ 1, 2, 3 }; // @Vector(3, f32) +const b = zml.Vec3f32{ 4, 5, 6 }; + +const d = zml.vec.dot(a, b); // f32 = 32 +const c = zml.vec.cross(a, b); // @Vector(3, f32) +const n = zml.vec.normalize(a); // @Vector(3, f32) +const r = zml.vec.sin_cos(a); // .{ .sin_out, .cos_out } +``` + +## `@Vector` and array inputs + +Every length-generic vector function accepts **both** a native `@Vector(N, T)` and a plain +`[N]T` array, and the return **preserves the input's container kind** (array in โ†’ array out, +vector in โ†’ vector out): + +```zig +const arr = [3]f32{ 1, 2, 3 }; +const d = zml.vec.dot(arr, arr); // f32 +const nz = zml.vec.normalize(arr); // [3]f32 +``` + +Arrays are meant for arbitrarily long data. Rather than coercing a large array into one wide +`@Vector(N, T)` โ€” which makes LLVM emit a single enormous SIMD instruction and bloats the binary +โ€” array inputs are processed in **chunks sized to the CPU's native SIMD width** +(`std.simd.suggestVectorLengthForCpu`) via a compact loop. `@Vector` inputs keep the single-op +path (you chose that width explicitly). + +For example, `norm` over a `[245]f32` compiles to **324 bytes** via the chunked path, versus +**2,529 bytes** for the equivalent `@Vector(245, f32)` op (`-OReleaseSmall`) โ€” and it is faster +too (see below). + +## Benchmarks + +Run the benchmark suite (uses [zBench](https://github.com/hendriknielaender/zBench)): + +```sh +zig build bench -Doptimize=ReleaseFast +``` + +Representative results (`ReleaseFast`; absolute numbers are machine-dependent โ€” the point is the +array-vs-`@Vector` ratio at large `N`): + +| Operation (N = 245) | array (chunked) | `@Vector(N)` (single op) | speedup | +| -------------------------- | --------------: | -----------------------: | ------: | +| `norm` | 37 ns | 235 ns | 6.3ร— | +| `dot` | 40 ns | 313 ns | 7.8ร— | +| `normalize` | 57 ns | 261 ns | 4.6ร— | + +| Sin/Cos over 256 elements | time/run | +| -------------------------- | ---------: | +| scalar `std.math` loop | 695 ns | +| `sin_cos` (`@Vector`) | 106 ns | +| `sin_cos` (array, chunked) | 139 ns | + +For the length-generic reductions/maps, the chunked array path is both **smaller and faster** +than one wide `@Vector` op at large `N`: the wide op forces LLVM into a slow, bloated instruction +sequence, while the native-width loop stays compact and vectorized. `sin_cos` is pure +element-wise, so the `@Vector` form already vectorizes cleanly and the two are comparable. + +## Testing + +```sh +zig build test +``` diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 78e5dc5..0000000 --- a/TODO.md +++ /dev/null @@ -1,11 +0,0 @@ -## TODO: (in this order~) - -- [x] rework how precision of output types is handled for vectors -- [x] rework how vectors handle integer types -- [x] add doc comments for Vector operations -- [x] add more operations to Vectors -- [x] rework how precision of output types is handled for vectors again -- [x] Add more operations to Matrices -- [ ] add testing to matrices operations -- [ ] add install guide -- [ ] add usage guide diff --git a/bench/bench.zig b/bench/bench.zig new file mode 100644 index 0000000..3ab25e2 --- /dev/null +++ b/bench/bench.zig @@ -0,0 +1,217 @@ +const std = @import("std"); +const zbench = @import("zbench"); +const zml = @import("zml"); + +fn benchmark_multiply(comptime size: usize) type { + return struct { + const Matrix512 = zml.Mat(f32, size, size); + a: Matrix512, + b: Matrix512, + + fn init() @This() { + var input: Matrix512 = undefined; + for (&input.items, 0..) |*row, i| { + for (row, 0..) |*col, j| { + col.* = @floatFromInt(i * size + j); + } + } + return .{ + .a = input, + .b = input, + }; + } + + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, Matrix512.mul, .{ self.a, self.b })); + } + }; +} + +fn bench_sin_cos_fused(comptime size: usize) type { + return struct { + angles: @Vector(size, f32), + + fn init() @This() { + var val: [size]f32 = undefined; + for (0..size) |k| { + val[k] = @as(f32, @floatFromInt(k)) * 0.01; + } + return .{ .angles = val }; + } + + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.sin_cos, .{self.angles})); + } + }; +} + +fn benchmark_sin_cos_system(comptime size: usize) type { + return struct { + angles: [size]f32, + + fn init() @This() { + var val: [size]f32 = undefined; + for (0..size) |k| { + val[k] = @as(f32, @floatFromInt(k)) * 0.01; + } + return .{ .angles = val }; + } + + pub fn run(self: *@This(), _: std.mem.Allocator) void { + var sin_val: [size]f32 = undefined; + var cos_val: [size]f32 = undefined; + for(self.angles, 0..) |angle, i| { + sin_val[i] = std.math.sin(angle); + cos_val[i] = std.math.cos(angle); + } + std.mem.doNotOptimizeAway(sin_val); + std.mem.doNotOptimizeAway(cos_val); + } + }; +} + +// Fill a length-`size` array with distinct positive values (kept away from zero +// so reductions/normalize stay well-conditioned). +fn ramp(comptime size: usize) [size]f32 { + var val: [size]f32 = undefined; + for (0..size) |k| val[k] = @as(f32, @floatFromInt(k % 251)) * 0.03 + 1.0; + return val; +} + +// `sin_cos` on a plain array: processed in native-width SIMD chunks (simd.zig). +fn bench_sin_cos_array(comptime size: usize) type { + return struct { + angles: [size]f32, + fn init() @This() { + return .{ .angles = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.sin_cos, .{self.angles})); + } + }; +} + +// The remaining factories each come in an `_array` (chunked) and a `_vector` +// (single wide @Vector op) flavor over the SAME data, so a run shows the +// chunked path's cost against the wide-vector path it replaces for arrays. + +fn bench_norm_array(comptime size: usize) type { + return struct { + data: [size]f32, + fn init() @This() { + return .{ .data = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.norm_adv, .{ self.data, 32 })); + } + }; +} + +fn bench_norm_vector(comptime size: usize) type { + return struct { + data: @Vector(size, f32), + fn init() @This() { + return .{ .data = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.norm_adv, .{ self.data, 32 })); + } + }; +} + +fn bench_dot_array(comptime size: usize) type { + return struct { + a: [size]f32, + b: [size]f32, + fn init() @This() { + return .{ .a = ramp(size), .b = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.dot, .{ self.a, self.b })); + } + }; +} + +fn bench_dot_vector(comptime size: usize) type { + return struct { + a: @Vector(size, f32), + b: @Vector(size, f32), + fn init() @This() { + return .{ .a = ramp(size), .b = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.dot, .{ self.a, self.b })); + } + }; +} + +fn bench_normalize_array(comptime size: usize) type { + return struct { + data: [size]f32, + fn init() @This() { + return .{ .data = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.normalize, .{self.data})); + } + }; +} + +fn bench_normalize_vector(comptime size: usize) type { + return struct { + data: @Vector(size, f32), + fn init() @This() { + return .{ .data = ramp(size) }; + } + pub fn run(self: *@This(), _: std.mem.Allocator) void { + std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.normalize, .{self.data})); + } + }; +} + +pub fn main() !void { + const io = std.Io.Threaded.global_single_threaded.io(); + const stdout: std.Io.File = .stdout(); + + var bench = zbench.Benchmark.init(std.heap.page_allocator, .{}); + defer bench.deinit(); + + try bench.addParam("Multiple 4x4 matrix multiplication", &benchmark_multiply(4).init(), .{ + .iterations = 256, + }); + try bench.addParam("Multiple 512x512 matrix multiplication", &benchmark_multiply(512).init(), .{ + .iterations = 256, + }); + try bench.addParam("Sin/Cos", &benchmark_sin_cos_system(256).init(), .{ + .iterations = 256, + }); + try bench.addParam("Sin/Cos vectorized", &bench_sin_cos_fused(256).init(), .{ + .iterations = 256, + }); + try bench.addParam("Sin/Cos array (chunked)", &bench_sin_cos_array(256).init(), .{ + .iterations = 256, + }); + + // Chunked array path vs single wide-@Vector op, at the motivating length 245. + try bench.addParam("norm [245] array (chunked)", &bench_norm_array(245).init(), .{ + .iterations = 4096, + }); + try bench.addParam("norm @Vector(245)", &bench_norm_vector(245).init(), .{ + .iterations = 4096, + }); + try bench.addParam("dot [245] array (chunked)", &bench_dot_array(245).init(), .{ + .iterations = 4096, + }); + try bench.addParam("dot @Vector(245)", &bench_dot_vector(245).init(), .{ + .iterations = 4096, + }); + try bench.addParam("normalize [245] array (chunked)", &bench_normalize_array(245).init(), .{ + .iterations = 4096, + }); + try bench.addParam("normalize @Vector(245)", &bench_normalize_vector(245).init(), .{ + .iterations = 4096, + }); + + try bench.run(io, stdout); +} + diff --git a/build.zig b/build.zig index 7b2fac6..6d70eb0 100644 --- a/build.zig +++ b/build.zig @@ -4,7 +4,7 @@ pub fn build(b: *std.Build) void { const target = b.standardTargetOptions(.{}); const optimize = b.standardOptimizeOption(.{}); - const root_module = b.addModule("zla", .{ + const root_module = b.addModule("zml", .{ .root_source_file = b.path("src/root.zig"), .target = target, .optimize = optimize, @@ -17,10 +17,28 @@ pub fn build(b: *std.Build) void { const test_step = b.step("test", "Run unit tests"); test_step.dependOn(&run_tests.step); } - - { - const tests_check = b.addTest(.{ .name = "check", .root_module = root_module }); - const check = b.step("check", "Check if exe and tests compile"); - check.dependOn(&tests_check.step); + + if (b.isRoot()) { + if (b.lazyDependency("zbench", .{ + .target = target, + .optimize = optimize, + })) |zbench_dep| { + const bench = b.addExecutable(.{ + .name = "bench", + .root_module = b.createModule(.{ + .root_source_file = b.path("bench/bench.zig"), + .target = target, + .optimize = optimize, + .imports = &.{ + .{ .name = "zml", .module = root_module }, + .{ .name = "zbench", .module = zbench_dep.module("zbench") }, + }, + }), + }); + const bench_step = b.step("bench", "run benchmark"); + const bench_cmd = b.addRunArtifact(bench); + bench_step.dependOn(&bench_cmd.step); + } } + } diff --git a/build.zig.zon b/build.zig.zon index deb2757..b0c361c 100644 --- a/build.zig.zon +++ b/build.zig.zon @@ -1,8 +1,14 @@ .{ - .name = .zla, - .version = "0.0.2", - .fingerprint = 0x55da41f1f0b5a9ce, // Changing this has security and trust implications. - .dependencies = .{}, + .name = .zml, + .version = "0.1.0", + .fingerprint = 0x32700c0d44ccd0f5, // Changing this has security and trust implications. + .dependencies = .{ + .zbench = .{ + .url = "git+https://github.com/hendriknielaender/zBench.git#f12a2262bcd9d300207890555c577da5e29de453", + .hash = "zbench-0.14.1-YTdc781OAQCR6wxqkgAwm9bG4TnFGzm8c9aE-ZPiTAU1", + .lazy = true, + }, + }, .paths = .{ "build.zig", "build.zig.zon", diff --git a/src/benchmark.zig b/src/benchmark.zig deleted file mode 100644 index 1e79dfa..0000000 --- a/src/benchmark.zig +++ /dev/null @@ -1,154 +0,0 @@ -const std = @import("std"); - -const module = @import("module.zig"); -const Matrix = module.Mat; - -test "benchmark matrix multiplication 256x256 * 256x256" { - const size = 256; - const Matrix512 = Matrix(f32, size, size); - const a = matrixMaker(size, size, 1); - const b = matrixMaker(size, size, 1); - - var timer = try std.time.Timer.start(); - - const n = 10; - for (0..n) |_| { - std.mem.doNotOptimizeAway(@call(.never_inline, Matrix512.mul, .{ &a, &b })); - } - - const elapsed = timer.read(); - std.debug.print("--------------------------------\n", .{}); - std.debug.print("matrix multiplication benchmark:\n", .{}); - std.debug.print("operation: 256x256 * 256x256 matrix multiplication\n", .{}); - std.debug.print("time per operation: {d} ns\n", .{ - @as(f64, @floatFromInt(elapsed)) / n, - }); - std.debug.print("time per operation: {d} ms\n", .{ - @as(f64, @floatFromInt(elapsed)) / n / 1000000, - }); -} - -test "benchmark matrix multiplication 8x8 * 8x8" { - const size = 8; - const Matrix8 = Matrix(f32, size, size); - const a = matrixMaker(size, size, 1); - const b = matrixMaker(size, size, 1); - - var timer = try std.time.Timer.start(); - - const n = 10; - for (0..n) |_| { - std.mem.doNotOptimizeAway(@call(.never_inline, Matrix8.mul, .{ &a, &b })); - } - - const elapsed = timer.read(); - std.debug.print("--------------------------------\n", .{}); - std.debug.print("matrix multiplication benchmark:\n", .{}); - std.debug.print("operation: 8x8 * 8x8 matrix multiplication\n", .{}); - std.debug.print("time per operation: {d} ns\n", .{ - @as(f64, @floatFromInt(elapsed)) / n, - }); - std.debug.print("time per operation: {d} ms\n", .{ - @as(f64, @floatFromInt(elapsed)) / n / 1000000, - }); -} - -test "benchmark matrix multiplication 4x4 * 4x4" { - const Mat4 = Matrix(f32, 4, 4); - - const n = 1e4; - const nn = 1e4; - - const as = try std.testing.allocator.alloc(Mat4, n); - defer std.testing.allocator.free(as); - for (as, 0..) |*item, i| item.* = matrixMaker(4, 4, i); - const bs = try std.testing.allocator.alloc(Mat4, n); - defer std.testing.allocator.free(bs); - for (bs, 0..) |*item, i| item.* = matrixMaker(4, 4, i); - - var timer = try std.time.Timer.start(); - - for (0..nn) |_| { - for (as, bs) |*a, b| { - a.* = @call(.always_inline, Mat4.mul, .{ a, b }); - } - } - - const elapsed = timer.read(); - std.debug.print("--------------------------------\n", .{}); - std.debug.print("matrix multiplication benchmark:\n", .{}); - std.debug.print("operation: 4x4 * 4x4 matrix multiplication\n", .{}); - std.debug.print("time per operation: {d} ns\n", .{ - @as(f64, @floatFromInt(elapsed)) / (n * nn), - }); -} - -test "benchmark matrix multiplication 3x3 * 3x3" { - const Mat3 = Matrix(f32, 3, 3); - - const n = 1e4; - const nn = 1e4; - - const as = try std.testing.allocator.alloc(Mat3, n); - defer std.testing.allocator.free(as); - for (as, 0..) |*item, i| item.* = matrixMaker(3, 3, i); - const bs = try std.testing.allocator.alloc(Mat3, n); - defer std.testing.allocator.free(bs); - for (bs, 0..) |*item, i| item.* = matrixMaker(3, 3, i); - - var timer = try std.time.Timer.start(); - - for (0..nn) |_| { - for (as, bs) |*a, b| { - a.* = @call(.always_inline, Mat3.mul, .{ a, b }); - } - } - - const elapsed = timer.read(); - std.debug.print("--------------------------------\n", .{}); - std.debug.print("matrix multiplication benchmark:\n", .{}); - std.debug.print("operation: 3x3 * 3x3 matrix multiplication\n", .{}); - std.debug.print("time per operation: {d} ns\n", .{ - @as(f64, @floatFromInt(elapsed)) / (n * nn), - }); -} - -test "benchmark matrix multiplication 2x2 * 2x2" { - const Mat2 = Matrix(f32, 2, 2); - - const n = 1e4; - const nn = 1e4; - - const as = try std.testing.allocator.alloc(Mat2, n); - defer std.testing.allocator.free(as); - for (as, 0..) |*item, i| item.* = matrixMaker(2, 2, i); - const bs = try std.testing.allocator.alloc(Mat2, n); - defer std.testing.allocator.free(bs); - for (bs, 0..) |*item, i| item.* = matrixMaker(2, 2, i); - - var timer = try std.time.Timer.start(); - - for (0..nn) |_| { - for (as, bs) |*a, b| { - a.* = @call(.always_inline, Mat2.mul, .{ a, &b }); - } - } - - const elapsed = timer.read(); - std.debug.print("--------------------------------\n", .{}); - std.debug.print("matrix multiplication benchmark:\n", .{}); - std.debug.print("operation: 2x2 * 2x2 matrix multiplication\n", .{}); - std.debug.print("time per operation: {d} ns\n", .{ - @as(f64, @floatFromInt(elapsed)) / (n * nn), - }); -} - -fn matrixMaker(comptime rows: usize, comptime cols: usize, added: usize) Matrix(f32, rows, cols) { - var result = Matrix(f32, rows, cols).fromCM(undefined); - for (&result.items, 0..) |*row, i| { - for (row, 0..) |*col, j| { - col.* = @floatFromInt(i * cols + j + added); - } - } - return result; -} diff --git a/src/color.zig b/src/color.zig new file mode 100644 index 0000000..0e49a2d --- /dev/null +++ b/src/color.zig @@ -0,0 +1,134 @@ +const std = @import("std"); + +const Vec3 = @Vector(3, f32); + +fn map3(v: Vec3, comptime f: fn (f32) f32) Vec3 { + return .{ f(v[0]), f(v[1]), f(v[2]) }; +} + +// --------------------------------------------------------------------------- +// sRGB transfer functions (IEC 61966-2-1) +// --------------------------------------------------------------------------- + +pub fn srgb_to_linear(c: f32) f32 { + return if (c <= 0.04045) c / 12.92 else std.math.pow(f32, (c + 0.055) / 1.055, 2.4); +} + +pub fn linear_to_srgb(c: f32) f32 { + return if (c <= 0.0031308) c * 12.92 else 1.055 * std.math.pow(f32, c, 1.0 / 2.4) - 0.055; +} + +pub fn srgb_to_linear3(c: Vec3) Vec3 { + return map3(c, srgb_to_linear); +} + +pub fn linear_to_srgb3(c: Vec3) Vec3 { + return map3(c, linear_to_srgb); +} + +// --------------------------------------------------------------------------- +// Luminance (Rec. 709 linear-light weights) +// --------------------------------------------------------------------------- + +pub fn luminance(rgb: Vec3) f32 { + return @reduce(.Add, rgb * Vec3{ 0.2126, 0.7152, 0.0722 }); +} + +// --------------------------------------------------------------------------- +// RGB <-> YCoCg (reversible luma/chroma transform) +// --------------------------------------------------------------------------- + +pub fn rgb_to_ycocg(rgb: Vec3) Vec3 { + return .{ + 0.25 * rgb[0] + 0.5 * rgb[1] + 0.25 * rgb[2], + 0.5 * rgb[0] - 0.5 * rgb[2], + -0.25 * rgb[0] + 0.5 * rgb[1] - 0.25 * rgb[2], + }; +} + +pub fn ycocg_to_rgb(ycocg: Vec3) Vec3 { + const y = ycocg[0]; + const co = ycocg[1]; + const cg = ycocg[2]; + return .{ y + co - cg, y + cg, y - co - cg }; +} + +// --------------------------------------------------------------------------- +// Tone mapping +// --------------------------------------------------------------------------- + +pub fn reinhard(x: Vec3) Vec3 { + return x / (@as(Vec3, @splat(1)) + x); +} + +pub fn reinhard_inverse(x: Vec3) Vec3 { + return x / (@as(Vec3, @splat(1)) - x); +} + +/// Narkowicz's fitted ACES filmic tone mapping curve. +pub fn aces_film(x: Vec3) Vec3 { + const a: Vec3 = @splat(2.51); + const b: Vec3 = @splat(0.03); + const c: Vec3 = @splat(2.43); + const d: Vec3 = @splat(0.59); + const e: Vec3 = @splat(0.14); + const num = x * (a * x + b); + const den = x * (c * x + d) + e; + return @min(@max(num / den, @as(Vec3, @splat(0))), @as(Vec3, @splat(1))); +} + +// --------------------------------------------------------------------------- +// PQ / SMPTE ST 2084 (Rec. 2100). Linear input normalized so 1.0 == 10000 nits. +// --------------------------------------------------------------------------- + +const pq_m1 = 0.1593017578125; +const pq_m2 = 78.84375; +const pq_c1 = 0.8359375; +const pq_c2 = 18.8515625; +const pq_c3 = 18.6875; + +pub fn linear_to_pq(l: f32) f32 { + const lm1 = std.math.pow(f32, @max(l, 0), pq_m1); + return std.math.pow(f32, (pq_c1 + pq_c2 * lm1) / (1.0 + pq_c3 * lm1), pq_m2); +} + +pub fn pq_to_linear(n: f32) f32 { + const nm2 = std.math.pow(f32, @max(n, 0), 1.0 / pq_m2); + return std.math.pow(f32, @max(nm2 - pq_c1, 0) / (pq_c2 - pq_c3 * nm2), 1.0 / pq_m1); +} + +test "srgb round trip" { + for ([_]f32{ 0, 0.001, 0.04, 0.5, 1 }) |c| { + try std.testing.expectApproxEqAbs(c, linear_to_srgb(srgb_to_linear(c)), 1e-5); + } + const v = Vec3{ 0.1, 0.5, 0.9 }; + const back = linear_to_srgb3(srgb_to_linear3(v)); + try std.testing.expectApproxEqAbs(v[0], back[0], 1e-5); + try std.testing.expectApproxEqAbs(v[2], back[2], 1e-5); +} + +test "luminance" { + try std.testing.expectApproxEqAbs(@as(f32, 1), luminance(.{ 1, 1, 1 }), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.7152), luminance(.{ 0, 1, 0 }), 1e-6); +} + +test "ycocg round trip" { + const rgb = Vec3{ 0.2, 0.6, 0.9 }; + const back = ycocg_to_rgb(rgb_to_ycocg(rgb)); + try std.testing.expectApproxEqAbs(rgb[0], back[0], 1e-6); + try std.testing.expectApproxEqAbs(rgb[1], back[1], 1e-6); + try std.testing.expectApproxEqAbs(rgb[2], back[2], 1e-6); +} + +test "reinhard round trip" { + const x = Vec3{ 0.1, 1.0, 4.0 }; + const back = reinhard_inverse(reinhard(x)); + try std.testing.expectApproxEqAbs(x[0], back[0], 1e-5); + try std.testing.expectApproxEqAbs(x[2], back[2], 1e-4); +} + +test "pq round trip" { + for ([_]f32{ 0.0, 0.01, 0.1, 0.5, 1.0 }) |l| { + try std.testing.expectApproxEqAbs(l, pq_to_linear(linear_to_pq(l)), 1e-4); + } +} diff --git a/src/geometry.zig b/src/geometry.zig new file mode 100644 index 0000000..1f5d295 --- /dev/null +++ b/src/geometry.zig @@ -0,0 +1,246 @@ +const std = @import("std"); +const vec = @import("vector.zig"); + +pub const Primative = enum { + AABB, + Plane, + Sphere, + OrientedBox, + Capsule, + Frustum, +}; + +/// Result of testing a primitive against a bounding volume such as a frustum. +pub const IntersectionState = enum { + outside, + inside, + partial, +}; + +pub const AABB = @import("geometry/aabb.zig").AABB; +pub const AABBf32 = AABB(f32); +pub const AABBf64 = AABB(f64); + +pub const Plane = @import("geometry/plane.zig").Plane; +pub const Planef32 = Plane(f32); +pub const Planef64 = Plane(f64); + +pub const Sphere = @import("geometry/sphere.zig").Sphere; +pub const Capsule = @import("geometry/capsule.zig").Capsule; +pub const OrientedBoundedBox = @import("geometry/obb.zig").OrientedBoundedBox; +pub const Frustum = @import("geometry/frustum.zig").Frustum; +pub const Frustumf32 = Frustum(f32); +pub const Frustumf64 = Frustum(f64); + + + +pub const overlap = @import("geometry/overlap.zig"); +pub const ray = @import("geometry/ray.zig"); +pub const contains = @import("geometry/contains.zig"); + +test { + _ = @import("geometry/aabb.zig"); + _ = @import("geometry/plane.zig"); + _ = @import("geometry/sphere.zig"); + _ = @import("geometry/capsule.zig"); + _ = @import("geometry/obb.zig"); + _ = @import("geometry/frustum.zig"); + _ = @import("geometry/overlap.zig"); + _ = @import("geometry/ray.zig"); + _ = @import("geometry/contains.zig"); +} + + +//pub fn get_vector_from_buffer(comptime T: type, vertex_index: usize, buffer: []align(4) const u8, byte_offset: usize, byte_stride: usize) T { +// const vector_len = switch(@typeInfo(T)) { +// .vector => |v| v.len, +// .array => |array| array.len, +// else => @compileError("Expected a vector type got: " ++ @typeName(T)), +// }; +// var arr: [vector_len]std.meta.Child(T) = undefined; +// const element_size = @sizeOf(std.meta.Child(T)); +// +// // Copy bytes and cast to the correct type +// const total_bytes = vector_len * element_size; +// const begin = byte_offset + vertex_index * byte_stride; +// const byte_slice = buffer[begin..begin + total_bytes]; +// @memcpy(std.mem.asBytes(&arr), byte_slice); +// return arr; +//} +// +//test "vector_from_buffer - basic 2D f32 vector extraction" { +// // Create a buffer with 2D f32 vectors: [1.0, 2.0], [3.0, 4.0], [5.0, 6.0] +// const data = [_]f32{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; +// const buffer = std.mem.sliceAsBytes(&data); +// +// // Extract first vector (index 0) +// const vec1 = get_vector_from_buffer(@Vector(2, f32), 0, buffer, 0, 2 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 1.0), vec1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), vec1[1]); +// +// // Extract second vector (index 1) +// const vec2 = get_vector_from_buffer(@Vector(2, f32), 1, buffer, 0, 2 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 3.0), vec2[0]); +// try std.testing.expectEqual(@as(f32, 4.0), vec2[1]); +// +// // Extract third vector (index 2) +// const vec3 = get_vector_from_buffer(@Vector(2, f32), 2, buffer, 0, 2 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 5.0), vec3[0]); +// try std.testing.expectEqual(@as(f32, 6.0), vec3[1]); +//} +// +//test "vector_from_buffer - 3D f32 vector extraction" { +// // Create a buffer with 3D f32 vectors: [1.0, 2.0, 3.0], [4.0, 5.0, 6.0] +// var data = [_]f32{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; +// const buffer = std.mem.sliceAsBytes(data[0..]); +// +// // Extract first 3D vector +// const vec1 = get_vector_from_buffer(@Vector(3, f32), 0, buffer, 0, 3 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 1.0), vec1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), vec1[1]); +// try std.testing.expectEqual(@as(f32, 3.0), vec1[2]); +// +// // Extract second 3D vector +// const vec2 = get_vector_from_buffer(@Vector(3, f32), 1, buffer, 0, 3 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 4.0), vec2[0]); +// try std.testing.expectEqual(@as(f32, 5.0), vec2[1]); +// try std.testing.expectEqual(@as(f32, 6.0), vec2[2]); +//} +// +//test "vector_from_buffer - 4D f32 vector extraction" { +// // Create a buffer with 4D f32 vectors: [1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0] +// const buffer = std.mem.asBytes(&[_]f32{ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 }); +// +// // Extract first 4D vector +// const vec1 = get_vector_from_buffer(@Vector(4, f32), 0, buffer, 0, 4 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 1.0), vec1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), vec1[1]); +// try std.testing.expectEqual(@as(f32, 3.0), vec1[2]); +// try std.testing.expectEqual(@as(f32, 4.0), vec1[3]); +// +// // Extract second 4D vector +// const vec2 = get_vector_from_buffer(@Vector(4, f32), 1, buffer, 0, 4 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 5.0), vec2[0]); +// try std.testing.expectEqual(@as(f32, 6.0), vec2[1]); +// try std.testing.expectEqual(@as(f32, 7.0), vec2[2]); +// try std.testing.expectEqual(@as(f32, 8.0), vec2[3]); +//} +// +//test "vector_from_buffer - integer vectors" { +// // Create a buffer with i32 vectors: [1, 2], [3, 4], [5, 6] +// const buffer = std.mem.asBytes(&[_]i32{ 1, 2, 3, 4, 5, 6 }); +// +// // Extract vectors +// const vec1 = get_vector_from_buffer(@Vector(2, i32), 0, buffer, 0, 2 * @sizeOf(i32)); +// try std.testing.expectEqual(@as(i32, 1), vec1[0]); +// try std.testing.expectEqual(@as(i32, 2), vec1[1]); +// +// const vec2 = get_vector_from_buffer(@Vector(2, i32), 1, buffer, 0, 2 * @sizeOf(i32)); +// try std.testing.expectEqual(@as(i32, 3), vec2[0]); +// try std.testing.expectEqual(@as(i32, 4), vec2[1]); +//} +// +//test "vector_from_buffer - with byte offset" { +// // Create a buffer with some padding at the beginning +// const data = [_]f32{ 999.0, 888.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; // First two are padding +// const buffer = std.mem.asBytes(&data); +// +// // Start reading from offset (skip the padding) +// const offset = 2 * @sizeOf(f32); +// +// const vec1 = get_vector_from_buffer(@Vector(2, f32), 0, buffer, offset, 2 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 1.0), vec1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), vec1[1]); +// +// const vec2 = get_vector_from_buffer(@Vector(2, f32), 1, buffer, offset, 2 * @sizeOf(f32)); +// try std.testing.expectEqual(@as(f32, 3.0), vec2[0]); +// try std.testing.expectEqual(@as(f32, 4.0), vec2[1]); +//} +// +//test "vector_from_buffer - interleaved data with larger stride" { +// // Create interleaved position and normal data: [pos_x, pos_y, norm_x, norm_y, pos_x, pos_y, norm_x, norm_y, ...] +// const data = [_]f32{ +// 1.0, 2.0, 0.1, 0.2, // vertex 0: pos(1,2), normal(0.1,0.2) +// 3.0, 4.0, 0.3, 0.4, // vertex 1: pos(3,4), normal(0.3,0.4) +// 5.0, 6.0, 0.5, 0.6 // vertex 2: pos(5,6), normal(0.5,0.6) +// }; +// const buffer = std.mem.asBytes(&data); +// +// // Extract positions (stride = 4 * f32, offset = 0) +// const stride = 4 * @sizeOf(f32); +// +// const pos1 = get_vector_from_buffer(@Vector(2, f32), 0, buffer, 0, stride); +// try std.testing.expectEqual(@as(f32, 1.0), pos1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), pos1[1]); +// +// const pos2 = get_vector_from_buffer(@Vector(2, f32), 1, buffer, 0, stride); +// try std.testing.expectEqual(@as(f32, 3.0), pos2[0]); +// try std.testing.expectEqual(@as(f32, 4.0), pos2[1]); +// +// // Extract normals (stride = 4 * f32, offset = 2 * f32) +// const normal_offset = 2 * @sizeOf(f32); +// +// const norm1 = get_vector_from_buffer(@Vector(2, f32), 0, buffer, normal_offset, stride); +// try std.testing.expectEqual(@as(f32, 0.1), norm1[0]); +// try std.testing.expectEqual(@as(f32, 0.2), norm1[1]); +// +// const norm2 = get_vector_from_buffer(@Vector(2, f32), 1, buffer, normal_offset, stride); +// try std.testing.expectEqual(@as(f32, 0.3), norm2[0]); +// try std.testing.expectEqual(@as(f32, 0.4), norm2[1]); +//} +// +//test "vector_from_buffer - f64 vectors" { +// // Test with double precision floats +// const buffer = std.mem.asBytes(&[_]f64{ 1.5, 2.5, 3.5, 4.5 }); +// +// const vec1 = get_vector_from_buffer(@Vector(2, f64), 0, buffer, 0, 2 * @sizeOf(f64)); +// try std.testing.expectEqual(@as(f64, 1.5), vec1[0]); +// try std.testing.expectEqual(@as(f64, 2.5), vec1[1]); +// +// const vec2 = get_vector_from_buffer(@Vector(2, f64), 1, buffer, 0, 2 * @sizeOf(f64)); +// try std.testing.expectEqual(@as(f64, 3.5), vec2[0]); +// try std.testing.expectEqual(@as(f64, 4.5), vec2[1]); +//} +// +//test "vector_from_buffer - mixed data types with padding" { +// // Test extracting from a buffer where vectors are not tightly packed +// const VtxData = extern struct { +// pos: @Vector(3, f32), +// padding1: u32, // 4 bytes padding +// normal: @Vector(3, f32), +// padding2: u32, // 4 bytes padding +// }; +// +// const vertices = [_]VtxData{ +// .{ .pos = .{ 1.0, 2.0, 3.0 }, .padding1 = 0, .normal = .{ 0.1, 0.2, 0.3 }, .padding2 = 0 }, +// .{ .pos = .{ 4.0, 5.0, 6.0 }, .padding1 = 0, .normal = .{ 0.4, 0.5, 0.6 }, .padding2 = 0 }, +// }; +// +// const buffer = std.mem.asBytes(&vertices); +// const vertex_stride = @sizeOf(VtxData); +// +// // Extract positions +// const pos_offset = @offsetOf(VtxData, "pos"); +// const pos1 = get_vector_from_buffer(@Vector(3, f32), 0, buffer, pos_offset, vertex_stride); +// try std.testing.expectEqual(@as(f32, 1.0), pos1[0]); +// try std.testing.expectEqual(@as(f32, 2.0), pos1[1]); +// try std.testing.expectEqual(@as(f32, 3.0), pos1[2]); +// +// const pos2 = get_vector_from_buffer(@Vector(3, f32), 1, buffer, pos_offset, vertex_stride); +// try std.testing.expectEqual(@as(f32, 4.0), pos2[0]); +// try std.testing.expectEqual(@as(f32, 5.0), pos2[1]); +// try std.testing.expectEqual(@as(f32, 6.0), pos2[2]); +// +// // Extract normals +// const normal_offset = @offsetOf(VtxData, "normal"); +// const norm1 = get_vector_from_buffer(@Vector(3, f32), 0, buffer, normal_offset, vertex_stride); +// try std.testing.expectEqual(@as(f32, 0.1), norm1[0]); +// try std.testing.expectEqual(@as(f32, 0.2), norm1[1]); +// try std.testing.expectEqual(@as(f32, 0.3), norm1[2]); +// +// const norm2 = get_vector_from_buffer(@Vector(3, f32), 1, buffer, normal_offset, vertex_stride); +// try std.testing.expectEqual(@as(f32, 0.4), norm2[0]); +// try std.testing.expectEqual(@as(f32, 0.5), norm2[1]); +// try std.testing.expectEqual(@as(f32, 0.6), norm2[2]); +//} +// diff --git a/src/geometry/aabb.zig b/src/geometry/aabb.zig new file mode 100644 index 0000000..4da14a2 --- /dev/null +++ b/src/geometry/aabb.zig @@ -0,0 +1,497 @@ +const std = @import("std"); +const zml = @import("../root.zig"); +const geometry = zml.geom; + +pub fn AABB(comptime T: type) type { + return struct { + pub const child: type = T; + pub const primative_type: geometry.Primative = .AABB; + const Self = @This(); + + min: @Vector(3, T), + max: @Vector(3, T), + + pub const empty = Self { + .min = .{0,0,0}, + .max = .{0,0,0}, + }; + + pub const InvDirection = struct { + inv_direction: @Vector(3, T), // 1 / ray direction + is_parallel: @Vector(3, bool), // true if ray direction is close to zero + + pub fn from_direction(direction: @Vector(3, T)) @This() { + return .{ + .is_parallel = @abs(direction) < @as(@Vector(3, T), @splat(1.0e-20)), + .inv_direction = @as(@Vector(3, T), @splat(1)) / direction, + }; + } + }; + + /// Get the closest point on or in this box to point + pub fn get_closest_point(self: Self, point: @Vector(3, T)) @Vector(3, T) { + return @min(@max(point, self.min), self.max); + } + + pub fn get_sqr_distance_to(self: Self, point: @Vector(3, T)) T { + return zml.vec.norm_sqr(self.get_closest_point(point) - point); + } + + pub fn from_two_points(p1: @Vector(3, T), p2: @Vector(3, T)) @This() { + std.debug.assert(@reduce(.And, p1 <= p2)); + return @This(){ + .min = @min(p1, p2), + .max = @max(p1, p2), + }; + } + + pub fn translate(self: Self, translation: @Vector(3, T)) Self { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return .{ + .min = self.min + translation, + .max = self.max + translation, + }; + } + + pub fn transform(self: Self, mat: zml.Mat(T, 4, 4)) Self { + const pos = mat.position(); + var new_min = pos; + var new_max = pos; + + inline for (0..3) |c| { + const m_col = mat.column(c); + const col: @Vector(3, T) = .{ m_col[0], m_col[1], m_col[2] }; + const a = col * @as(@Vector(3, T), @splat(self.min[c])); + const b = col * @as(@Vector(3, T), @splat(self.max[c])); + new_min += @min(a, b); + new_max += @max(a, b); + } + return .{ .min = new_min, .max = new_max }; + } + + pub fn get_center(self: @This()) @Vector(3, T) { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return (self.min + self.max) * @as(@Vector(3, T), @splat(0.5)); + } + + pub fn get_size(self: @This()) @Vector(3, T) { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return self.max - self.min; + } + + // Calculate the support vector for this convex shape + pub fn get_support(self: @This(), direction: @Vector(3, T)) @Vector(3, T) { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return @select(T, direction < @as(@Vector(3, T), @splat(0)), self.max, self.min); + } + + pub fn encapsulate_aabb(self: Self, a: Self) Self { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return .{ .min = @min(self.min, a.min), .max = @max(self.max, a.max) }; + } + + pub fn intersect(self: Self, inRHS: Self) Self { + std.debug.assert(@reduce(.And, self.min <= self.max)); + return Self{ + .min = @max(self.min, inRHS.min), + .max = @min(self.max, inRHS.max), + }; + } + + pub fn expand_by(self: *Self, in: @Vector(3, T)) void { + std.debug.assert(@reduce(.And, self.min <= self.max)); + self.min -= in; + self.max += in; + } + + pub fn surface_area(self: Self) T { + std.debug.assert(@reduce(.And, self.min <= self.max)); + const extent = self.max - self.min; + return 2 * zml.vec.dot(zml.vec.swizzle(extent, "yzz"), zml.vec.swizzle(extent, "xxy")); + } + + pub fn volume(self: Self) T { + std.debug.assert(@reduce(.And, self.min <= self.max)); + const extent = self.max - self.min; + return @reduce(.Mul, extent); + } + }; +} + +test "aabb_transform" { + const AABBf32 = AABB(f32); + const aabb = AABBf32.from_two_points(.{ 0.0, 0.0, 0.0 }, .{ 1.0, 1.0, 1.0 }); + var translation: zml.Mat(f32, 4, 4) = .identity; + translation = translation.translate(.{ 1.0, 2.0, 3.0 }); + const transformed = aabb.transform(translation); + try std.testing.expect(zml.vec.is_close_default(transformed.min, .{ 1.0, 2.0, 3.0 })); + try std.testing.expect(zml.vec.is_close_default(transformed.max, .{ 2.0, 3.0, 4.0 })); +} + +test "surface_area" { + const AABBf32 = AABB(f32); + const aabb = AABBf32.from_two_points(.{ 0.0, 0.0, 0.0 }, .{ 1.0, 2.0, 3.0 }); + try std.testing.expectApproxEqRel(aabb.surface_area(), 22.0, 1.0e-6); +} + +test "volumn" { + const AABBf32 = AABB(f32); + const aabb = AABBf32.from_two_points(.{ 0.0, 0.0, 0.0 }, .{ 1.0, 2.0, 3.0 }); + try std.testing.expectApproxEqRel(aabb.volume(), 6.0, 1.0e-6); +} + +test "InvDirection" { + const AABBf32 = AABB(f32); + const dir = @Vector(3, f32){ 1.0, 0.0, -1.0 }; + const invDir = AABBf32.InvDirection.from_direction(dir); + try std.testing.expectEqual(invDir.is_parallel, @Vector(3, bool){ false, true, false }); + try std.testing.expectApproxEqRel(invDir.inv_direction[0], 1.0, 1.0e-6); + try std.testing.expectApproxEqRel(invDir.inv_direction[1], 0.0, 1.0e-6); + try std.testing.expectApproxEqRel(invDir.inv_direction[2], -1.0, 1.0e-6); +} + +test "aabb_get_support" { + const aabb: AABB(f32) = .from_two_points(.{ 0.0, 0.0, 0.0 }, .{ 1.0, 1.0, 1.0 }); + + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ 0.5774, 0.5774, 0.5774 }), .{ 0.0, 0.0, 0.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ 0.5774, 0.5774, -0.5774 }), .{ 0.0, 0.0, 1.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ 0.5774, -0.5774, 0.5774 }), .{ 0.0, 1.0, 0.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ 0.5774, -0.5774, -0.5774 }), .{ 0.0, 1.0, 1.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ -0.5774, 0.5774, 0.5774 }), .{ 1.0, 0.0, 0.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ -0.5774, 0.5774, -0.5774 }), .{ 1.0, 0.0, 1.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ -0.5774, -0.5774, 0.5774 }), .{ 1.0, 1.0, 0.0 })); + try std.testing.expect(zml.vec.is_close_default(aabb.get_support(.{ -0.5774, -0.5774, -0.5774 }), .{ 1.0, 1.0, 1.0 })); +} + +//const Plane = struct { +// normal: @Vector(3, f32), +// distance: f32, +// +// pub fn getNormal(self: @This()) @Vector(3, f32) { +// return self.normal; +// } +// +// pub fn signedDistance(self: @This(), point: @Vector(3, f32)) f32 { +// return @reduce(.Add, self.normal * point) + self.distance; +// } +//}; +// +//const AABox = @This(); +//min: @Vector(3, f32), +//max: @Vector(3, f32), +// +///// Create box from 2 points +//pub fn from_two_points(p1: @Vector(3, f32), p2: @Vector(3, f32)) AABox { +// return AABox{ +// .min = @min(p1, p2), +// .max = @min(p1, p2), +// }; +//} +// +///// Create box from indexed triangle +//pub fn from_triangles(vertices: []f32, indexes: []u32) AABox { +// var box = from_two_points(vertices[indexes[0]], vertices[indexes[1]]); +// box.encapsulate(vertices[indexes[2]]); +// return box; +//} +// +///// Get bounding box of size FLT_MAX +//pub fn biggest() AABox { +// // Max half extent of AABox is 0.5 * FLT_MAX so that getSize() remains finite +// const half_max = 0.5 * math.floatMax(f32); +// return AABox{ +// .min = vec3Replicate(-half_max), +// .max = vec3Replicate(half_max), +// }; +//} +// +///// Comparison operators +//pub fn eql(self: AABox, other: AABox) bool { +// return @reduce(.And, self.min == other.min) and @reduce(.And, self.max == other.max); +//} +// +//pub fn neql(self: AABox, other: AABox) bool { +// return !self.eql(other); +//} +// +///// Reset the bounding box to an empty bounding box +//pub fn setEmpty(self: *AABox) void { +// self.min = vec3Replicate(math.floatMax(f32)); +// self.max = vec3Replicate(-math.floatMax(f32)); +//} +// +///// Check if the bounding box is valid (max >= min) +//pub fn isValid(self: AABox) bool { +// return self.min[0] <= self.max[0] and +// self.min[1] <= self.max[1] and +// self.min[2] <= self.max[2]; +//} +// +///// Encapsulate point in bounding box +//pub fn encapsulate(self: *AABox, pos: Vec3) void { +// self.min = vec3Min(self.min, pos); +// self.max = vec3Max(self.max, pos); +//} +// +///// Encapsulate bounding box in bounding box +//pub fn encapsulateBox(self: *AABox, other: AABox) void { +// self.min = vec3Min(self.min, other.min); +// self.max = vec3Max(self.max, other.max); +//} +// +///// Encapsulate triangle in bounding box +//pub fn encapsulateTriangle(self: *AABox, triangle: Triangle) void { +// self.encapsulate(Vec3{ triangle.v[0][0], triangle.v[0][1], triangle.v[0][2] }); +// self.encapsulate(Vec3{ triangle.v[1][0], triangle.v[1][1], triangle.v[1][2] }); +// self.encapsulate(Vec3{ triangle.v[2][0], triangle.v[2][1], triangle.v[2][2] }); +//} +// +///// Encapsulate indexed triangle in bounding box +//pub fn encapsulateIndexedTriangle(self: *AABox, vertices: VertexList, triangle: IndexedTriangle) void { +// for (triangle.idx) |idx| { +// self.encapsulate(vertices[idx]); +// } +//} +// +///// Intersect this bounding box with other, returns the intersection +//pub fn intersect(self: AABox, other: AABox) AABox { +// return AABox{ +// .min = vec3Max(self.min, other.min), +// .max = vec3Min(self.max, other.max), +// }; +//} +// +///// Make sure that each edge of the bounding box has a minimal length +//pub fn ensureMinimalEdgeLength(self: *AABox, min_edge_length: f32) void { +// const min_length = vec3Replicate(min_edge_length); +// const size = self.max - self.min; +// const mask = vec3Less(size, min_length); +// self.max = vec3Select(mask, self.min + min_length, self.max); +//} +// +///// Widen the box on both sides by vector +//pub fn expandBy(self: *AABox, vector: Vec3) void { +// self.min -= vector; +// self.max += vector; +//} +// +///// Get center of bounding box +//pub fn getCenter(self: AABox) Vec3 { +// return (self.min + self.max) * vec3Replicate(0.5); +//} +// +///// Get extent of bounding box (half of the size) +//pub fn getExtent(self: AABox) Vec3 { +// return (self.max - self.min) * vec3Replicate(0.5); +//} +// +///// Get size of bounding box +//pub fn getSize(self: AABox) Vec3 { +// return self.max - self.min; +//} +// +///// Get surface area of bounding box +//pub fn getSurfaceArea(self: AABox) f32 { +// const extent = self.max - self.min; +// return 2.0 * (extent[0] * extent[1] + extent[0] * extent[2] + extent[1] * extent[2]); +//} +// +///// Get volume of bounding box +//pub fn getVolume(self: AABox) f32 { +// const extent = self.max - self.min; +// return extent[0] * extent[1] * extent[2]; +//} +// +///// Check if this box contains another box +//pub fn contains(self: AABox, other: AABox) bool { +// return testAllXYZTrue(vec3LessOrEqual(self.min, other.min)) and +// testAllXYZTrue(vec3GreaterOrEqual(self.max, other.max)); +//} +// +///// Check if this box contains a point +//pub fn containsPoint(self: AABox, point: Vec3) bool { +// return testAllXYZTrue(vec3LessOrEqual(self.min, point)) and +// testAllXYZTrue(vec3GreaterOrEqual(self.max, point)); +//} +// +///// Check if this box contains a double precision point +//pub fn containsPointD(self: AABox, point: DVec3) bool { +// const point_f32 = Vec3{ @floatCast(point[0]), @floatCast(point[1]), @floatCast(point[2]) }; +// return self.containsPoint(point_f32); +//} +// +///// Check if this box overlaps with another box +//pub fn overlaps(self: AABox, other: AABox) bool { +// return !testAnyXYZTrue(vec3Greater(self.min, other.max)) and +// !testAnyXYZTrue(vec3Less(self.max, other.min)); +//} +// +///// Check if this box overlaps with a plane +//pub fn overlapsPlane(self: AABox, plane: Plane) bool { +// const normal = plane.getNormal(); +// const dist_normal = plane.signedDistance(self.getSupport(normal)); +// const dist_min_normal = plane.signedDistance(self.getSupport(-normal)); +// return dist_normal * dist_min_normal <= 0.0; // If both support points are on the same side of the plane we don't overlap +//} +// +///// Translate bounding box +//pub fn translate(self: *AABox, translation: Vec3) void { +// self.min += translation; +// self.max += translation; +//} +// +///// Translate bounding box with double precision +//pub fn translateD(self: *AABox, translation: DVec3) void { +// const min_d = DVec3{ self.min[0], self.min[1], self.min[2] } + translation; +// const max_d = DVec3{ self.max[0], self.max[1], self.max[2] } + translation; +// +// // Round down for min, round up for max to ensure conservative bounds +// self.min = Vec3{ @floatCast(min_d[0]), @floatCast(min_d[1]), @floatCast(min_d[2]) }; +// self.max = Vec3{ @floatCast(max_d[0]), @floatCast(max_d[1]), @floatCast(max_d[2]) }; +//} +// +///// Transform bounding box by 4x4 matrix +//pub fn transformed(self: AABox, matrix: Mat44) AABox { +// // Start with the translation of the matrix +// var new_min = Vec3{ matrix[3][0], matrix[3][1], matrix[3][2] }; +// var new_max = new_min; +// +// // Now find the extreme points by considering the product of the min and max with each column of matrix +// var c: u32 = 0; +// while (c < 3) : (c += 1) { +// const col = Vec3{ matrix[0][c], matrix[1][c], matrix[2][c] }; +// +// const a = col * vec3Replicate(self.min[c]); +// const b = col * vec3Replicate(self.max[c]); +// +// new_min += vec3Min(a, b); +// new_max += vec3Max(a, b); +// } +// +// return AABox{ .min = new_min, .max = new_max }; +//} +// +///// Transform bounding box by double precision 4x4 matrix +//pub fn transformedD(self: AABox, matrix: DMat44) AABox { +// // Extract rotation part as f32 matrix +// var rotation: Mat44 = undefined; +// for (0..3) |i| { +// for (0..3) |j| { +// rotation[i][j] = @floatCast(matrix[i][j]); +// } +// } +// rotation[3] = .{ 0, 0, 0, 1 }; +// +// var result = self.transformed(rotation); +// const translation = DVec3{ matrix[3][0], matrix[3][1], matrix[3][2] }; +// result.translateD(translation); +// return result; +//} +// +///// Scale this bounding box, can handle non-uniform and negative scaling +//pub fn scaled(self: AABox, scale: Vec3) AABox { +// return fromTwoPoints(self.min * scale, self.max * scale); +//} +// +///// Calculate the support vector for this convex shape +//pub fn getSupport(self: AABox, direction: Vec3) Vec3 { +// return vec3Select(vec3Less(direction, vec3Zero()), self.min, self.max); +//} +// +///// Get the vertices of the face that faces direction the most +//pub fn getSupportingFace(self: AABox, direction: Vec3, vertices: *[4]Vec3) void { +// const axis = getHighestComponentIndex(direction); +// +// if (direction[axis] < 0.0) { +// switch (axis) { +// 0 => { +// vertices[0] = Vec3{ self.max[0], self.min[1], self.min[2] }; +// vertices[1] = Vec3{ self.max[0], self.max[1], self.min[2] }; +// vertices[2] = Vec3{ self.max[0], self.max[1], self.max[2] }; +// vertices[3] = Vec3{ self.max[0], self.min[1], self.max[2] }; +// }, +// 1 => { +// vertices[0] = Vec3{ self.min[0], self.max[1], self.min[2] }; +// vertices[1] = Vec3{ self.min[0], self.max[1], self.max[2] }; +// vertices[2] = Vec3{ self.max[0], self.max[1], self.max[2] }; +// vertices[3] = Vec3{ self.max[0], self.max[1], self.min[2] }; +// }, +// 2 => { +// vertices[0] = Vec3{ self.min[0], self.min[1], self.max[2] }; +// vertices[1] = Vec3{ self.max[0], self.min[1], self.max[2] }; +// vertices[2] = Vec3{ self.max[0], self.max[1], self.max[2] }; +// vertices[3] = Vec3{ self.min[0], self.max[1], self.max[2] }; +// }, +// else => unreachable, +// } +// } else { +// switch (axis) { +// 0 => { +// vertices[0] = Vec3{ self.min[0], self.min[1], self.min[2] }; +// vertices[1] = Vec3{ self.min[0], self.min[1], self.max[2] }; +// vertices[2] = Vec3{ self.min[0], self.max[1], self.max[2] }; +// vertices[3] = Vec3{ self.min[0], self.max[1], self.min[2] }; +// }, +// 1 => { +// vertices[0] = Vec3{ self.min[0], self.min[1], self.min[2] }; +// vertices[1] = Vec3{ self.max[0], self.min[1], self.min[2] }; +// vertices[2] = Vec3{ self.max[0], self.min[1], self.max[2] }; +// vertices[3] = Vec3{ self.min[0], self.min[1], self.max[2] }; +// }, +// 2 => { +// vertices[0] = Vec3{ self.min[0], self.min[1], self.min[2] }; +// vertices[1] = Vec3{ self.min[0], self.max[1], self.min[2] }; +// vertices[2] = Vec3{ self.max[0], self.max[1], self.min[2] }; +// vertices[3] = Vec3{ self.max[0], self.min[1], self.min[2] }; +// }, +// else => unreachable, +// } +// } +//} +// +///// Get the closest point on or in this box to point +//pub fn getClosestPoint(self: AABox, point: Vec3) Vec3 { +// return vec3Min(vec3Max(point, self.min), self.max); +//} +// +///// Get the squared distance between point and this box (will be 0 if point is inside the box) +//pub fn getSqDistanceTo(self: AABox, point: Vec3) f32 { +// const closest = self.getClosestPoint(point); +// return vec3LengthSq(closest - point); +//} +// +//// Tests +//test "AABox basic functionality" { +// const testing = std.testing; +// +// // Test creation from two points +// const p1 = Vec3{ 0, 0, 0 }; +// const p2 = Vec3{ 1, 1, 1 }; +// const box = fromTwoPoints(p1, p2); +// +// try testing.expectEqual(Vec3{ 0, 0, 0 }, box.min); +// try testing.expectEqual(Vec3{ 1, 1, 1 }, box.max); +// +// // Test center calculation +// const center = box.getCenter(); +// try testing.expectEqual(Vec3{ 0.5, 0.5, 0.5 }, center); +// +// // Test size calculation +// const size = box.getSize(); +// try testing.expectEqual(Vec3{ 1, 1, 1 }, size); +// +// // Test contains point +// try testing.expect(box.containsPoint(Vec3{ 0.5, 0.5, 0.5 })); +// try testing.expect(!box.containsPoint(Vec3{ 2, 2, 2 })); +//} +// +//test "AABox encapsulation" { +// const testing = std.testing; +// +// var box = fromTwoPoints(Vec3{ 0, 0, 0 }, Vec3{ 1, 1, 1 }); +// box.encapsulate(Vec3{ 2, 0.5, 0.5 }); +// +// try testing.expectEqual(Vec3{ 0, 0, 0 }, box.min); +// try testing.expectEqual(Vec3{ 2, 1, 1 }, box.max); +//} +// diff --git a/src/geometry/buffer_view.zig b/src/geometry/buffer_view.zig new file mode 100644 index 0000000..28251dc --- /dev/null +++ b/src/geometry/buffer_view.zig @@ -0,0 +1,9 @@ +const std = @import("std"); +const meta = @import("../meta.zig"); + +pub fn BufferView(comptime T: type) type { + return struct { + pub const inner_type = [meta.lengthOf(T)]meta.Child(T); + pub const num_elements = meta.lengthOf(T); + }; +} diff --git a/src/geometry/capsule.zig b/src/geometry/capsule.zig new file mode 100644 index 0000000..84088be --- /dev/null +++ b/src/geometry/capsule.zig @@ -0,0 +1,26 @@ +const geometry = @import("../geometry.zig"); +const zml = @import("../root.zig"); + +pub fn Capsule(comptime T: type) type { + return struct { + pub const child: type = T; + pub const primative_type = geometry.Primative.Capsule; + const Self = @This(); + + hemisphere_centers: [2]@Vector(3, T), // the two hemisphere centers + radius: T, // radius of the capsule + + pub fn center(self: Self) @Vector(3, T) { + return (self.hemisphere_centers[0] + self.hemisphere_centers[1]) * @as(@Vector(3, T), @splat(0.5)); + } + + pub fn get_cylinder_height(self: Self) T { + return zml.vec.distance(self.hemisphere_centers[0], self.hemisphere_centers[1]); + } + + pub fn get_total_height(self: Self) T { + return self.get_cylinder_height() + self.radius * @as(T, 2); + } + }; +} + diff --git a/src/geometry/contains.zig b/src/geometry/contains.zig new file mode 100644 index 0000000..d7dd555 --- /dev/null +++ b/src/geometry/contains.zig @@ -0,0 +1,39 @@ +const std = @import("std"); +const zml = @import("../root.zig"); + +pub fn aabb_contains_point(a: anytype, pt: @Vector(3, @TypeOf(a).child)) bool { + comptime { + std.debug.assert(@TypeOf(a).primative_type == .AABB); + } + return @reduce(.And, (pt >= a.min) & (pt <= a.max)); +} + +//pub fn capsule_contains_point(a: anytype, pt: @Vector(3, @TypeOf(a).child)) bool { +// comptime { +// std.debug.assert(@TypeOf(a).primative_type == .Capsule); +// } +// const ab = a.hemisphere_centers[1] - a.hemisphere_centers[0]; +// const t = @max(@min(zml.vec.dot(pt - a.hemisphere_centers[0], ab) / zml.vec.dot(ab, ab), @as(@TypeOf(a).child, 1)), @as(@TypeOf(a).child, 0)); +// const closest_point = a.hemisphere_centers[0] + ab * @as(@Vector(3, @TypeOf(a).child), t); +// return zml.vec.distance_sqr(pt, closest_point) <= a.radius * a.radius; +//} + +test aabb_contains_point { + const aabb: zml.geom.AABB(f32) = .from_two_points(.{ -1, -1, -1 }, .{ 1, 1, 1 }); + try std.testing.expect(aabb_contains_point(aabb, .{ 0, 0, 0 })); + try std.testing.expect(!aabb_contains_point(aabb, .{ 2, 0, 0 })); + try std.testing.expect(!aabb_contains_point(aabb, .{ 0, -2, 0 })); + try std.testing.expect(!aabb_contains_point(aabb, .{ 0, 0, 2 })); +} + +//test capsule_contains_point { +// const capsule: zml.geom.Capsule(f32) = .{ +// .hemisphere_centers = .{ .{ 0, 0, -1 }, .{ 0, 0, 1 } }, +// .radius = 1, +// }; +// try std.testing.expect(capsule_contains_point(capsule, .{ 0, 0, 0 })); +// try std.testing.expect(capsule_contains_point(capsule, .{ 1, 0, 0 })); +// try std.testing.expect(!capsule_contains_point(capsule, .{ 2, 0, 0 })); +// try std.testing.expect(!capsule_contains_point(capsule, .{ 0, -2, 0 })); +// try std.testing.expect(!capsule_contains_point(capsule, .{ 0, 0, 3 })); +//} diff --git a/src/geometry/cylinder.zig b/src/geometry/cylinder.zig new file mode 100644 index 0000000..e69de29 diff --git a/src/geometry/frustum.zig b/src/geometry/frustum.zig new file mode 100644 index 0000000..653b072 --- /dev/null +++ b/src/geometry/frustum.zig @@ -0,0 +1,120 @@ +const std = @import("std"); +const zml = @import("../root.zig"); +const Mat = @import("../matrix.zig").Mat; +const vec = @import("../vector.zig"); +const geometry = @import("../geometry.zig"); +const Plane = @import("plane.zig").Plane; + +/// A view frustum represented by its six bounding planes with inward-pointing +/// normals: left, right, bottom, top, near, far. Built from a view-projection +/// matrix using the Gribb-Hartmann method for a right-handed, [0, 1]-depth +/// projection (the convention produced by `Mat.perspective`). +pub fn Frustum(comptime T: type) type { + return struct { + pub const child: type = T; + pub const primative_type: geometry.Primative = .Frustum; + const Self = @This(); + + planes: [6]Plane(T), + + fn make_plane(v: @Vector(4, T)) Plane(T) { + const n = @Vector(3, T){ v[0], v[1], v[2] }; + const inv_len = 1.0 / vec.norm(n); + return .{ + .normal = n * @as(@Vector(3, T), @splat(inv_len)), + .c = v[3] * inv_len, + }; + } + + /// Extract the six frustum planes from a (model-)view-projection matrix. + pub fn from_view_projection(mvp: Mat(T, 4, 4)) Self { + const r0 = mvp.row(0); + const r1 = mvp.row(1); + const r2 = mvp.row(2); + const r3 = mvp.row(3); + return .{ + .planes = .{ + make_plane(r3 + r0), // left + make_plane(r3 - r0), // right + make_plane(r3 + r1), // bottom + make_plane(r3 - r1), // top + make_plane(r2), // near (zero-to-one depth) + make_plane(r3 - r2), // far + }, + }; + } + + pub fn intersect_sphere_state(self: Self, sphere: anytype) geometry.IntersectionState { + comptime std.debug.assert(@TypeOf(sphere).primative_type == .Sphere); + var inside = true; + for (self.planes) |p| { + const d = p.signed_distance(sphere.center); + if (d < -sphere.radius) return .outside; + if (d < sphere.radius) inside = false; + } + return if (inside) .inside else .partial; + } + + pub fn intersect_sphere(self: Self, sphere: anytype) bool { + return self.intersect_sphere_state(sphere) != .outside; + } + + pub fn intersect_aabb_state(self: Self, aabb: anytype) geometry.IntersectionState { + comptime std.debug.assert(@TypeOf(aabb).primative_type == .AABB); + var inside = true; + for (self.planes) |p| { + // AABB.get_support(d) returns the corner minimizing dot(d, corner), + // so get_support(-normal) is the corner furthest along +normal. + if (p.signed_distance(aabb.get_support(-p.normal)) < 0) return .outside; + if (p.signed_distance(aabb.get_support(p.normal)) < 0) inside = false; + } + return if (inside) .inside else .partial; + } + + pub fn intersect_aabb(self: Self, aabb: anytype) bool { + return self.intersect_aabb_state(aabb) != .outside; + } + + pub fn intersect_capsule_state(self: Self, capsule: anytype) geometry.IntersectionState { + comptime std.debug.assert(@TypeOf(capsule).primative_type == .Capsule); + var inside = true; + for (self.planes) |p| { + const d0 = p.signed_distance(capsule.hemisphere_centers[0]); + const d1 = p.signed_distance(capsule.hemisphere_centers[1]); + if (@max(d0, d1) < -capsule.radius) return .outside; + if (@min(d0, d1) < capsule.radius) inside = false; + } + return if (inside) .inside else .partial; + } + + pub fn intersect_capsule(self: Self, capsule: anytype) bool { + return self.intersect_capsule_state(capsule) != .outside; + } + }; +} + +test "frustum culling" { + const M = Mat(f32, 4, 4).orthographic(-1, 1, -1, 1, 1, 10); + const f = Frustum(f32).from_view_projection(M); + const Sphere = zml.geom.Sphere; + const AABB = zml.geom.AABB; + const Capsule = zml.geom.Capsule; + const State = geometry.IntersectionState; + + // Spheres + try std.testing.expectEqual(State.inside, f.intersect_sphere_state(Sphere(f32).from_center_radius(.{ 0, 0, -5 }, 0.1))); + try std.testing.expectEqual(State.outside, f.intersect_sphere_state(Sphere(f32).from_center_radius(.{ 5, 0, -5 }, 0.1))); + try std.testing.expectEqual(State.outside, f.intersect_sphere_state(Sphere(f32).from_center_radius(.{ 0, 0, -0.5 }, 0.1))); + try std.testing.expectEqual(State.partial, f.intersect_sphere_state(Sphere(f32).from_center_radius(.{ 1, 0, -5 }, 0.5))); + try std.testing.expect(f.intersect_sphere(Sphere(f32).from_center_radius(.{ 0, 0, -5 }, 0.1))); + try std.testing.expect(!f.intersect_sphere(Sphere(f32).from_center_radius(.{ 5, 0, -5 }, 0.1))); + + // AABBs + try std.testing.expectEqual(State.inside, f.intersect_aabb_state(AABB(f32).from_two_points(.{ -0.5, -0.5, -6 }, .{ 0.5, 0.5, -4 }))); + try std.testing.expectEqual(State.outside, f.intersect_aabb_state(AABB(f32).from_two_points(.{ 2, 2, -6 }, .{ 3, 3, -4 }))); + try std.testing.expectEqual(State.partial, f.intersect_aabb_state(AABB(f32).from_two_points(.{ 0.5, 0.5, -6 }, .{ 1.5, 1.5, -4 }))); + + // Capsules + try std.testing.expect(f.intersect_capsule(Capsule(f32){ .hemisphere_centers = .{ .{ 0, 0, -4 }, .{ 0, 0, -6 } }, .radius = 0.2 })); + try std.testing.expect(!f.intersect_capsule(Capsule(f32){ .hemisphere_centers = .{ .{ 5, 0, -4 }, .{ 5, 0, -6 } }, .radius = 0.2 })); +} diff --git a/src/geometry/obb.zig b/src/geometry/obb.zig new file mode 100644 index 0000000..5296486 --- /dev/null +++ b/src/geometry/obb.zig @@ -0,0 +1,12 @@ +const Mat = @import("../matrix.zig").Mat; +const geometry = @import("../geometry.zig"); + +pub fn OrientedBoundedBox(comptime T: type) type { + return struct { + pub const primative_type = geometry.Primative.OrientedBox; + pub const child: type = T; + orientation: Mat(4,4, T), + half_extent: @Vector(3, T), + }; +} + diff --git a/src/geometry/overlap.zig b/src/geometry/overlap.zig new file mode 100644 index 0000000..c196fcc --- /dev/null +++ b/src/geometry/overlap.zig @@ -0,0 +1,155 @@ +const std = @import("std"); +const zml = @import("../root.zig"); +const Mat = @import("../matrix.zig").Mat; +const geometry = @import("../geometry.zig"); +const Sphere = @import("sphere.zig").Sphere; +const vector = @import("../vector.zig"); + +// aabb +pub fn overlap_sphere_sphere(a: anytype, b: anytype) bool { + comptime { + std.debug.assert(@TypeOf(a).child == @TypeOf(b).child); + std.debug.assert(@TypeOf(a).primative_type == .Sphere); + std.debug.assert(@TypeOf(b).primative_type == .Sphere); + } + return vector.norm_sqr(a.center - b.center) <= (a.radius + b.radius) * (a.radius + b.radius); +} + +pub fn overlap_aabb_aabb(a: anytype, b: anytype) bool { + comptime { + std.debug.assert(@TypeOf(a).child == @TypeOf(b).child); + std.debug.assert(@TypeOf(a).primative_type == .AABB); + std.debug.assert(@TypeOf(b).primative_type == .AABB); + } + return !@reduce(.Or, (a.min > b.max) | (a.max < b.min)); +} + +pub fn overlap_aabb_plane(a: anytype, b: anytype) bool { + comptime { + std.debug.assert(@TypeOf(a).child == @TypeOf(b).child); + std.debug.assert(@TypeOf(a).primative_type == .AABB); + std.debug.assert(@TypeOf(b).primative_type == .Plane); + } + const dist_normal = b.signed_distance(a.get_support(b.normal)); + const dist_min_normal = b.signed_distance(a.get_support(-b.normal)); + return dist_normal * dist_min_normal <= 0; +} + +pub fn overlap_aabb_aabb_4(a: anytype, minX: @Vector(4, @TypeOf(a).inner_type), maxX: @Vector(4, @TypeOf(a).inner_type), minY: @Vector(4, @TypeOf(a).inner_type), maxY: @Vector(4, @TypeOf(a).inner_type), minZ: @Vector(4, @TypeOf(a).inner_type), maxZ: @Vector(4, @TypeOf(a).inner_type)) @Vector(4, bool) { + comptime { + std.debug.assert(@TypeOf(a).primative_type == .AABB); + } + const box1_minx = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.min[0])); + const box1_miny = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.min[1])); + const box1_minz = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.min[2])); + const box1_maxx = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.max[0])); + const box1_maxy = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.max[1])); + const box1_maxz = @as(@Vector(4, @TypeOf(a).inner_type), @splat(a.max[2])); + + const nooverlap_x = (box1_minx > maxX) | (box1_maxx < minX); + const nooverlap_y = (box1_miny > maxY) | (box1_maxy < minY); + const nooverlap_z = (box1_minz > maxZ) | (box1_maxz < minZ); + return !(nooverlap_x | nooverlap_y | nooverlap_z); +} + +/// Axis-aligned box vs triangle overlap using the separating-axis theorem +/// (Akenine-Mรถller). The triangle is given by its three vertices. +pub fn overlap_aabb_triangle( + aabb: anytype, + t0: @Vector(3, @TypeOf(aabb).child), + t1: @Vector(3, @TypeOf(aabb).child), + t2: @Vector(3, @TypeOf(aabb).child), +) bool { + comptime std.debug.assert(@TypeOf(aabb).primative_type == .AABB); + const T = @TypeOf(aabb).child; + const V = @Vector(3, T); + + const c = aabb.get_center(); + const h = aabb.get_size() * @as(V, @splat(0.5)); + const v0 = t0 - c; + const v1 = t1 - c; + const v2 = t2 - c; + const edges = [3]V{ v1 - v0, v2 - v1, v0 - v2 }; + const verts = [3]V{ v0, v1, v2 }; + + // 9 axes: cross products of the triangle edges with the box axes. + inline for (0..3) |ei| { + inline for (0..3) |ai| { + var unit: V = .{ 0, 0, 0 }; + unit[ai] = 1; + const axis = vector.cross(edges[ei], unit); + var pmin: T = std.math.inf(T); + var pmax: T = -std.math.inf(T); + inline for (0..3) |k| { + const p = vector.dot(axis, verts[k]); + pmin = @min(pmin, p); + pmax = @max(pmax, p); + } + const r = vector.dot(h, @abs(axis)); + if (pmin > r or pmax < -r) return false; + } + } + + // 3 box face normals: the triangle's AABB must overlap the box. + inline for (0..3) |i| { + const tmin = @min(v0[i], @min(v1[i], v2[i])); + const tmax = @max(v0[i], @max(v1[i], v2[i])); + if (tmin > h[i] or tmax < -h[i]) return false; + } + + // 1 axis: the triangle plane normal. + const n = vector.cross(edges[0], edges[1]); + if (@abs(vector.dot(n, v0)) > vector.dot(h, @abs(n))) return false; + + return true; +} + +// generic overlap function +pub inline fn overlap(a: anytype, b: anytype) bool { + const a_primative: geometry.Primative = @TypeOf(a).primative_type; + const b_primative: geometry.Primative = @TypeOf(b).primative_type; + if (a_primative == .Sphere and b_primative == .Sphere) return overlap_sphere_sphere(a, b); + if (a_primative == .AABB and b_primative == .AABB) return overlap_aabb_aabb(a, b); + @compileError("Unsupported primative overlap: " ++ @typeName(a_primative) ++ " " ++ @typeName(b_primative)); +} + +test overlap_aabb_triangle { + const box = geometry.AABB(f32).from_two_points(.{ -1, -1, -1 }, .{ 1, 1, 1 }); + // Triangle slicing through the box. + try std.testing.expect(overlap_aabb_triangle(box, .{ -2, 0, 0 }, .{ 2, 0, 0 }, .{ 0, 2, 0 })); + // Triangle with a vertex inside the box. + try std.testing.expect(overlap_aabb_triangle(box, .{ 0, 0, 0 }, .{ 5, 0, 0 }, .{ 0, 5, 0 })); + // Triangle far outside the box. + try std.testing.expect(!overlap_aabb_triangle(box, .{ 5, 5, 5 }, .{ 6, 5, 5 }, .{ 5, 6, 5 })); + // Triangle separated by a face plane (all z > 1). + try std.testing.expect(!overlap_aabb_triangle(box, .{ 0, 0, 2 }, .{ 1, 0, 3 }, .{ 0, 1, 2.5 })); +} + +test overlap_sphere_sphere { + const s1: Sphere(f32) = .from_center_radius(.{ 0, 0, 0 }, 1); + const s2: Sphere(f32) = .from_center_radius(.{ 0, 0, 1.5 }, 1); + const s3: Sphere(f32) = .from_center_radius(.{ 0, 0, 3 }, 1); + + try std.testing.expect(overlap_sphere_sphere(s1, s2)); + try std.testing.expect(!overlap_sphere_sphere(s1, s3)); + // Symmetric + try std.testing.expect(overlap(s1, s2)); + try std.testing.expect(overlap(s2, s1)); +} + +test overlap_aabb_aabb { + const aabb1: geometry.AABB(f32) = .from_two_points(.{ 0, 0, 0 }, .{ 2, 2, 2 }); + const aabb2: geometry.AABB(f32) = .from_two_points(.{ 1, 1, 1 }, .{ 3, 3, 3 }); // Overlapping + const aabb3: geometry.AABB(f32) = .from_two_points(.{ 5, 5, 5 }, .{ 7, 7, 7 }); // Non-overlapping + const aabb4: geometry.AABB(f32) = .from_two_points(.{ 2, 0, 0 }, .{ 4, 2, 2 }); // Edge touching + + try std.testing.expect(overlap_aabb_aabb(aabb1, aabb2)); // Overlapping boxes should return true + try std.testing.expect(!overlap_aabb_aabb(aabb1, aabb3)); // Non-overlapping boxes should return false + try std.testing.expect(overlap_aabb_aabb(aabb1, aabb4)); // Edge touching boxes should return true + + + // Symmetric - test the generic overlap function + try std.testing.expect(overlap(aabb1, aabb2)); + try std.testing.expect(!overlap(aabb1, aabb3)); + try std.testing.expect(overlap(aabb2, aabb1)); +} diff --git a/src/geometry/plane.zig b/src/geometry/plane.zig new file mode 100644 index 0000000..0721854 --- /dev/null +++ b/src/geometry/plane.zig @@ -0,0 +1,126 @@ +const std = @import("std"); +const vec = @import("../vector.zig"); +const Mat = @import("../matrix.zig").Mat; + +const zml = @import("../root.zig"); +const geometry = zml.geom; + +pub fn Plane(comptime T: type) type { + return struct { + pub const child: type = T; + pub const primative_type: geometry.Primative = .Plane; + + normal: @Vector(3, T), // normal vector + c: T, // constant + + const Self = @This(); + + pub fn normal_and_constant(self: Self) @Vector(4, T) { + return .{ self.normal[0], self.normal[1], self.normal[2], self.c }; + } + + pub fn from_point_and_normal(point: @Vector(3, T), normal: @Vector(3, T)) Self { + return .{ .normal = normal, .c = -vec.dot(normal, point) }; + } + + pub fn from_points_ccw(a: @Vector(3, T), b: @Vector(3, T), c: @Vector(3, T)) Self { + return from_point_and_normal(a, vec.normalize(vec.cross((b - a), (c - a)))); + } + + pub fn offset(self: Self, distance: T) Self { + return .{ + .normal = self.normal, + .c = self.c - distance, + }; + } + + pub fn transform(self: Self, m: Mat(T, 4, 4)) Self { + const transformed_normal = m.extract(3, 3).mul(vec.to_mat(self.normal)); + return .{ .normal = transformed_normal.column(0), .c = self.c - vec.dot(m.position(), transformed_normal.column(0)) }; + } + + pub fn scaled(self: Self, scale: @Vector(3, T)) Self { + const scaled_normal = self.normal / scale; + const scaled_normal_length = vec.norm(scaled_normal); + return .{ + .normal = self.normal / @as(@Vector(3, T), @splat(scaled_normal_length)), + .c = self.c / scaled_normal_length, + }; + } + + // distance point to plane + pub fn signed_distance(self: Self, pt: @Vector(3, T)) T { + return vec.dot(self.normal, pt) + self.c; + } + + pub fn project_point_plane(self: Self, pt: @Vector(3, T)) @Vector(3, T) { + return pt - self.normal * @as(@Vector(3, T), @splat(self.signed_distance(pt))); + } + + pub fn intersect_plane(p1: Self, p2: Self, p3: Self) ?@Vector(3, T) { + // We solve the equation: + // |ax, ay, az, aw| | x | | 0 | + // |bx, by, bz, bw| * | y | = | 0 | + // |cx, cy, cz, cw| | z | | 0 | + // |0, 0, 0, 1| | 1 | | 1 | + // Where normal of plane 1 = (ax, ay, az), plane constant of 1 = aw, normal of plane 2 = (bx, by, bz) etc. + // This involves inverting the matrix and multiplying it with [0, 0, 0, 1] + + // Fetch the normals and plane constants for the three planes + const a = p1.normal_and_constant(); + const b = p2.normal_and_constant(); + const c = p3.normal_and_constant(); + + const denom = vec.dot(p1.normal, vec.cross(p2.normal, p3.normal)); + if (denom == 0) return null; + + // The numerator is: + // [aw*(bz*cy-by*cz)+ay*(bw*cz-bz*cw)+az*(by*cw-bw*cy)] + // [aw*(bx*cz-bz*cx)+ax*(bz*cw-bw*cz)+az*(bw*cx-bx*cw)] + // [aw*(by*cx-bx*cy)+ax*(bw*cy-by*cw)+ay*(bx*cw-bw*cx)] + const numerator = + @as(@Vector(3, T), @splat(p1.c)) * (vec.swizzle(b, "zxy") * vec.swizzle(c, "yzx") - vec.swizzle(b, "yzx") * vec.swizzle(c, "zxy")) + vec.swizzle(a, "yxx") * (vec.swizzle(b, "wzw") * vec.swizzle(c, "zwy") - vec.swizzle(b, "zwy") * vec.swizzle(c, "wzw")) + vec.swizzle(a, "zzy") * (vec.swizzle(b, "ywx") * vec.swizzle(c, "wxw") - vec.swizzle(b, "wxw") * vec.swizzle(c, "ywx")); + return numerator / @as(@Vector(3, T), @splat(denom)); + } + }; +} + +test "transformed" { + const transform = Mat(f32, 4, 4) + .rotate(.identity, 0.1 * std.math.pi, vec.normalize(@Vector(3, f32){ 1, 2, 3 })) + .translate(@Vector(3, f32){ 5, -7, 9 }); + + const point = @Vector(3, f32){ 11.0, 13.0, 15.0 }; + const normal = vec.normalize(@Vector(3, f32){ -3.0, 5.0, -7.0 }); + + const p1 = Plane(f32).from_point_and_normal(point, normal).transform(transform); + const p2 = Plane(f32).from_point_and_normal(vec.extract(transform.mul( + vec.to_mat(@Vector(4, f32){ point[0], point[1], point[2], 1.0 })).column(0), 3), + transform.extract(3, 3).mul(vec.to_mat(normal)).column(0)); + + try std.testing.expectApproxEqAbs(p1.normal[0], p2.normal[0], 0.000001); + try std.testing.expectApproxEqAbs(p1.normal[1], p2.normal[1], 0.000001); + try std.testing.expectApproxEqAbs(p1.normal[2], p2.normal[2], 0.000001); +} + +test "signed_distance" { + const plane = Plane(f32).from_point_and_normal(.{ 0, 2, 0 }, .{ 0, 1, 0 }); + try std.testing.expectApproxEqRel(plane.signed_distance(.{ 5, 7, 0 }), 5.0, 0.0001); + try std.testing.expectApproxEqRel(plane.signed_distance(.{ 5, -3, 0 }), -5.0, 0.0001); +} + +test "intersect_plane" { + const p1 = Plane(f32).from_point_and_normal(.{ 0, 2, 0 }, .{ 0, 1, 0 }); + const p2 = Plane(f32).from_point_and_normal(.{ 3, 0, 0 }, .{ 1, 0, 0 }); + const p3 = Plane(f32).from_point_and_normal(.{ 0, 0, 4 }, .{ 0, 0, 1 }); + { + const point = Plane(f32).intersect_plane(p1, p2, p3); + try std.testing.expect(point != null); + try std.testing.expectEqual(point.?, .{ 3, 2, 4 }); + } + { + const p4 = Plane(f32).from_point_and_normal(.{ 0, 3, 0 }, .{ 0, 1, 0 }); + const point = Plane(f32).intersect_plane(p1, p2, p4); + try std.testing.expect(point == null); + } +} diff --git a/src/geometry/ray.zig b/src/geometry/ray.zig new file mode 100644 index 0000000..f5f88cc --- /dev/null +++ b/src/geometry/ray.zig @@ -0,0 +1,264 @@ +const std = @import("std"); +const zml = @import("../root.zig"); + +pub fn InvDirection(comptime T: type) type { + return struct { + inv_direction: @Vector(3, T), // 1 / ray direction + is_parallel: @Vector(3, bool), // true if ray direction is close to zero + + pub fn from_direction(direction: @Vector(3, T)) @This() { + return .{ + .is_parallel = @abs(direction) < @as(@Vector(3, T), @splat(1.0e-20)), + .inv_direction = @as(@Vector(3, T), @splat(1)) / direction, + }; + } + }; +} + +//pub fn ray_cylinder(cylinder: anytype, origin: @Vector(3, @TypeOf(cylinder).child), direction: @Vector(3, @TypeOf(cylinder).child)) std.meta.Float(@bitSizeOf(@TypeOf(cylinder).child)) { +// comptime { +// std.debug.assert(@TypeOf(cylinder).primative_type == .Cylinder); // ensure cylinder type +// } +// +// const orgin_xz = @select(@TypeOf(cylinder).child, .{ true, false, true }, origin, @as(@Vector(3, @TypeOf(cylinder).child), @splat(0))); +// const origin_xz_len_sq = zml.vec.norm_sqr(orgin_xz); +// const r_sq = cylinder.radius * cylinder.radius; +// if (origin_xz_len_sq > r_sq) { +// // Ray starts outside the infinite cylinder +// // Solve: |RayOrigin_xz + fraction * RayDirection_xz|^2 = r^2 to find fraction +// const direction_xz = @select(@TypeOf(cylinder).child, .{ true, false, true }, direction, @as(@Vector(3, @TypeOf(cylinder).child), @splat(0))); +// const a = zml.vec.norm_sqr(direction_xz); +// const b = 2 * zml.vec.dot(orgin_xz, direction_xz); +// const c = origin_xz_len_sq - r_sq; +// const root_terms = zml.find_roots(@TypeOf(cylinder).child, a, b, c); +// if (root_terms.num_roots == 0) { +// return std.math.floatMax(@TypeOf(cylinder).child); +// } +// // Get fraction corresponding to the ray entering the circle +// const fraction = @min(root_terms.roots[0], root_terms.roots[1]); +// if (fraction >= 0) { +// return fraction; +// } +// } else { +// return 0.0; +// } +// return std.math.floatMax(@TypeOf(cylinder).child); +//} + +// Ray - Axis Aligned Bounding Box intersection +// Note: Can return negative t values if the ray origin is inside the AABB +// return std.math.floatMax(T) if no hit +pub fn ray_aabb(aabb: anytype, origin: @Vector(3, @TypeOf(aabb).child), invDir: InvDirection(@TypeOf(aabb).child)) std.meta.Float(@bitSizeOf(@TypeOf(aabb).child)) { + comptime { + std.debug.assert(@TypeOf(aabb).primative_type == .AABB); // ensure aabb type + } + const flt_min = @as(@Vector(3, @TypeOf(aabb).child), @splat(-std.math.floatMax(@TypeOf(aabb).child))); + const flt_max = @as(@Vector(3, @TypeOf(aabb).child), @splat(std.math.floatMax(@TypeOf(aabb).child))); + + // Test against all three axes simultaneously. + const t1 = (aabb.min - origin) * invDir.inv_direction; + const t2 = (aabb.max - origin) * invDir.inv_direction; + + // Compute the max of min(t1,t2) and the min of max(t1,t2) ensuring we don't + // use the results from any directions parallel to the slab. + var t_min = @select(@TypeOf(aabb).child, invDir.is_parallel, flt_min, @min(t1, t2)); + var t_max = @select(@TypeOf(aabb).child, invDir.is_parallel, flt_max, @max(t1, t2)); + + // t_min.xyz = maximum(t_min.x, t_min.y, t_min.z); + t_min = @max(t_min, @shuffle(@TypeOf(aabb).child, t_min, t_min, [3]u8{ 1, 2, 0 })); + t_min = @max(t_min, @shuffle(@TypeOf(aabb).child, t_min, t_min, [3]u8{ 2, 0, 1 })); + + // t_max.xyz = minimum(t_max.x, t_max.y, t_max.z); + t_max = @min(t_max, @shuffle(@TypeOf(aabb).child, t_max, t_max, [3]u8{ 1, 2, 0 })); + t_max = @min(t_max, @shuffle(@TypeOf(aabb).child, t_max, t_max, [3]u8{ 2, 0, 1 })); + + // if (t_min > t_max) return FLT_MAX; + var no_intersection = t_min > t_max; + + // if (t_max < 0.0f) return FLT_MAX; + no_intersection = no_intersection | (t_max < @as(@TypeOf(t_max), @splat(0))); + + // if (inInvDirection.mIsParallel && !(Min <= inOrigin && inOrigin <= Max)) return FLT_MAX; else return t_min; + const no_parallel_overlap = (origin < aabb.min) | (origin > aabb.max); + no_intersection = no_intersection | (invDir.is_parallel & no_parallel_overlap); + no_intersection = no_intersection | @as(@TypeOf(no_intersection), @splat(no_intersection[1])); + no_intersection = no_intersection | @as(@TypeOf(no_intersection), @splat(no_intersection[2])); + + return @select(@TypeOf(aabb).child, no_intersection, flt_max, t_min)[0]; +} + +// Ray - Axis Aligned Bounding Box intersection +// Note: Can return negative t values if the ray origin is inside the AABB +// return -std.math.floatMax(T) and std.math.floatMax(T) if no hit +pub fn ray_aabb_with_enter_exit(aabb: anytype, origin: @Vector(3, @TypeOf(aabb).child), invDir: InvDirection(@TypeOf(aabb).child)) struct { + min: std.meta.Float(@bitSizeOf(@TypeOf(aabb).child)), + max: std.meta.Float(@bitSizeOf(@TypeOf(aabb).child)) +} { + comptime { + std.debug.assert(@TypeOf(aabb).primative_type == .AABB); // ensure aabb type + } + const flt_min = @as(@Vector(3, @TypeOf(aabb).child), @splat(-std.math.floatMax(@TypeOf(aabb).child))); + const flt_max = @as(@Vector(3, @TypeOf(aabb).child), @splat(std.math.floatMax(@TypeOf(aabb).child))); + + // Test against all three axes simultaneously. + const t1 = (aabb.min - origin) * invDir.inv_direction; + const t2 = (aabb.max - origin) * invDir.inv_direction; + + // Compute the max of min(t1,t2) and the min of max(t1,t2) ensuring we don't + // use the results from any directions parallel to the slab. + var t_min = @select(@TypeOf(aabb).child, invDir.is_parallel, flt_min, @min(t1, t2)); + var t_max = @select(@TypeOf(aabb).child, invDir.is_parallel, flt_max, @max(t1, t2)); + + // t_min.xyz = maximum(t_min.x, t_min.y, t_min.z); + t_min = @max(t_min, @shuffle(@TypeOf(aabb).child, t_min, t_min, [3]u8{ 1, 2, 0 })); + t_min = @max(t_min, @shuffle(@TypeOf(aabb).child, t_min, t_min, [3]u8{ 2, 0, 1 })); + + // t_max.xyz = minimum(t_max.x, t_max.y, t_max.z); + t_max = @min(t_max, @shuffle(@TypeOf(aabb).child, t_max, t_max, [3]u8{ 1, 2, 0 })); + t_max = @min(t_max, @shuffle(@TypeOf(aabb).child, t_max, t_max, [3]u8{ 2, 0, 1 })); + + // if (t_min > t_max) return FLT_MAX; + var no_intersection = t_min > t_max; + + // if (t_max < 0.0f) return FLT_MAX; + no_intersection = no_intersection | (t_max < @as(@TypeOf(t_max), @splat(0))); + + // if (inInvDirection.mIsParallel && !(Min <= inOrigin && inOrigin <= Max)) return FLT_MAX; else return t_min; + const no_parallel_overlap = (origin < aabb.min) | (origin > aabb.max); + no_intersection = no_intersection | (invDir.is_parallel & no_parallel_overlap); + no_intersection = no_intersection | @as(@TypeOf(no_intersection), @splat(no_intersection[1])); + no_intersection = no_intersection | @as(@TypeOf(no_intersection), @splat(no_intersection[2])); + + return .{ + .min = @select(@TypeOf(aabb).child, no_intersection, flt_max, t_min)[0], + .max = @select(@TypeOf(aabb).child, no_intersection, flt_max, t_max)[0], + }; + +} + +/// Intersect ray with triangle, returns closest point or FLT_MAX if no hit (branch less version) +/// Adapted from: http://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm +pub fn ray_triangle(comptime T: type, origin: @Vector(3, T), direction: @Vector(3, T), v0: @Vector(3, T), v1: @Vector(3, T), v2: @Vector(3, T)) T { + const epsilon: @Vector(3, T) = @as(@Vector(3, T), @splat(1.0e-12)); + + const zero: @Vector(3, T) = @as(@Vector(3, T), @splat(0)); + const one: @Vector(3, T) = @as(@Vector(3, T), @splat(1)); + + // Find vectors for two edges sharing v0 + const e1 = v1 - v0; + const e2 = v2 - v0; + + // Begin calculating determinant - also used to calculate u parameter + const p = zml.vec.cross(direction, e2); + + // If determinant is near zero, ray lies in plane of triangle + var det = @as(@Vector(3, T), @splat(zml.vec.dot(e1, p))); + + // Check if determinant is near zero + const det_near_zero = @abs(det) < epsilon; + + // when the determinant is near zero, return no intersection + det = @select(T, det_near_zero, one, det); + + // Calculate distance from v0 to ray origin + const s = origin - v0; + + // Calculate u parameter and test bounds + const u = @as(@Vector(3, T), @splat(zml.vec.dot(s, p))) / det; + + // Prepare to test v parameter + const q = zml.vec.cross(s, e1); + + // Calculate v parameter and test bounds + const v = @as(@Vector(3, T), @splat(zml.vec.dot(direction, q))) / det; + + // get intersection point + const t = @as(@Vector(3, T), @splat(zml.vec.dot(e2, q))) / det; + + const no_intersection = + (det_near_zero | (u < zero)) | ((v < zero) | ((u + v) > one)) | (t < zero); + + return @select(T, no_intersection, @as(@Vector(3, T), @splat(std.math.floatMax(T))), t)[0]; +} + +test ray_aabb { + const aabb: zml.geom.AABB(f32) = .from_two_points(.{ -1, -1, -1 }, .{ 1, 1, 1 }); + inline for (0..3) |axis| { + { + // Ray starting in the center of the box, pointing high + const origin = @Vector(3, f32){ 0, 0, 0 }; + var direction = @Vector(3, f32){ 0, 0, 0 }; + direction[axis] = 1; + const fraction = ray_aabb(aabb, origin, .from_direction(direction)); + try std.testing.expectApproxEqRel(-1.0, fraction, 1.0e-6); + } + } +} + +test ray_triangle { + const v0 = @Vector(3, f32){ 0, 0, 0 }; + const v1 = @Vector(3, f32){ 1, 0, 0 }; + const v2 = @Vector(3, f32){ 0, 1, 0 }; + + { + // Ray starting above the triangle, pointing down + const origin = @Vector(3, f32){ 0.25, 0.25, 1 }; + const direction = @Vector(3, f32){ 0, 0, -1 }; + const fraction = ray_triangle(f32, origin, direction, v0, v1, v2); + try std.testing.expectApproxEqRel(1.0, fraction, 1.0e-6); + } + { + // Ray starting below the triangle, pointing up + const origin = @Vector(3, f32){ 0.25, 0.25, -1 }; + const direction = @Vector(3, f32){ 0, 0, 1 }; + const fraction = ray_triangle(f32, origin, direction, v0, v1, v2); + try std.testing.expectApproxEqRel(1.0, fraction, 1.0e-6); + } + { + // Ray starting to the side of the triangle pointing away + const origin = @Vector(3, f32){ -1, -1, 0 }; + const direction = @Vector(3, f32){ -1, -1, 0 }; + const fraction = ray_triangle(f32, origin, direction, v0, v1, v2); + try std.testing.expectEqual(std.math.floatMax(f32), fraction); + } + { + // Ray starting to the side of the triangle pointing towards + const origin = @Vector(3, f32){ -1, -1, 0 }; + const direction = @Vector(3, f32){ 1, 1, 0 }; + const fraction = ray_triangle(f32, origin, direction, v0, v1, v2); + try std.testing.expectEqual(std.math.floatMax(f32), fraction); + } +} + +//test ray_cylinder { +// const cylinder: zml.geom.Cylinder(f32) = .from_two_points_radius(.{ 0, 0, 0 }, .{ 0, 0, 2 }, 1.0); +// { +// // Ray starting outside the cylinder, pointing towards +// const origin = @Vector(3, f32){ 2, 0, 1 }; +// const direction = @Vector(3, f32){ -1, 0, 0 }; +// const fraction = ray_cylinder(cylinder, origin, direction); +// try std.testing.expectApproxEqRel(1.0, fraction, 1.0e-6); +// } +// { +// // Ray starting outside the cylinder, pointing away +// const origin = @Vector(3, f32){ 2, 0, 1 }; +// const direction = @Vector(3, f32){ 1, 0, 0 }; +// const fraction = ray_cylinder(cylinder, origin, direction); +// try std.testing.expectEqual(std.math.floatMax(f32), fraction); +// } +// { +// // Ray starting inside the cylinder +// const origin = @Vector(3, f32){ 0.5, 0, 1 }; +// const direction = @Vector(3, f32){ 1, 0, 0 }; +// const fraction = ray_cylinder(cylinder, origin, direction); +// try std.testing.expectApproxEqRel(0.0, fraction, 1.0e-6); +// } +//} + +//pub fn ray_aabb(origin: anytype) { +// +//} + +//pub fn ray_intersect_aabb() { +// +//} diff --git a/src/geometry/segment.zig b/src/geometry/segment.zig new file mode 100644 index 0000000..e69de29 diff --git a/src/geometry/sphere.zig b/src/geometry/sphere.zig new file mode 100644 index 0000000..fce39d3 --- /dev/null +++ b/src/geometry/sphere.zig @@ -0,0 +1,25 @@ +const zml = @import("../root.zig"); +const geometry = zml.geom; + +pub fn Sphere(comptime T: type) type { + return struct { + pub const child: type = T; + pub const primative_type: geometry.Primative = .Sphere; + + pub const Self = @This(); + pub const empty: Self = .{ + .center = @Vector(3, T){0, 0, 0} , + .radius = 0, + }; + + center: @Vector(3, T), + radius: T, + + pub fn from_center_radius(center: @Vector(3, T), radius: T) Sphere(T) { + return .{ + .center = center, + .radius = radius, + }; + } + }; +} diff --git a/src/geometry/test.zig b/src/geometry/test.zig new file mode 100644 index 0000000..e69de29 diff --git a/src/matrix.zig b/src/matrix.zig index fbb7c1a..3e87f8c 100644 --- a/src/matrix.zig +++ b/src/matrix.zig @@ -1,33 +1,32 @@ const std = @import("std"); const vec = @import("./root.zig").vec; +pub const Mat4f32 = Mat(f32, 4, 4); +pub const Mat4f64 = Mat(f64, 4, 4); + /// column major generic matrix type pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type { return extern struct { const Self = @This(); - - items: [cols][rows]T, - - // comptime rows: comptime_int = rows, TODO: open an issue about comptime fields not being allowed on extern structs - // comptime cols: comptime_int = cols, - // comptime Type: type = T, - pub const rows: comptime_int = rows_; pub const cols: comptime_int = cols_; pub const Type: type = T; + pub const is_square: bool = rows == cols; - pub inline fn fromCM(values: [cols][rows]T) Self { + items: [cols][rows]T, + + pub inline fn from_column_major_array(values: [cols][rows]T) Self { return .{ .items = values }; } /// performs a transpose operation, usefull for more human readable mat literals - pub inline fn fromRM(values: [rows][cols]T) Self { - return Mat(T, rows, cols).fromCM(values).transpose(); + pub inline fn from_row_major_array(values: [rows][cols]T) Self { + return Mat(T, rows, cols).from_column_major_array(values).transpose(); } // generated code seams fast but tbd pub fn transpose(self: Self) Mat(T, rows, cols) { - var result: Mat(T, rows, cols) = .fromCM(undefined); + var result: Mat(T, rows, cols) = .from_column_major_array(undefined); for (0..cols) |c| { for (0..rows) |r| { result.items[r][c] = self.items[c][r]; @@ -37,21 +36,22 @@ pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type } /// Scalar multiplication - pub fn scalarMul(self: Self, scalar: T) Self { + pub fn scalar_mul(self: Self, scalar: T) Self { const items: [rows * cols]T = @bitCast(self.items); var result_items: [rows * cols]T = undefined; for (&result_items, items) |*result_item, item| { result_item.* = item * scalar; } - return .fromCM(@bitCast(result_items)); + return .from_column_major_array(@bitCast(result_items)); } - // `other` can be a scalar or a matrix with the same number of rows as the numbers of columns of `self` + // `other` must be a matrix with the same number of rows as the numbers of columns of `self` pub fn mul(self: Self, other: anytype) Mat(T, @TypeOf(other).cols, Self.rows) { - if (Self.cols != @TypeOf(other).rows) @compileError("number of columns of self must be equal to number of rows of other"); - if (Self.Type != @TypeOf(other).Type) @compileError("type of self must be equal to Type of other"); + comptime { + std.debug.assert(Self.cols == @TypeOf(other).rows); + std.debug.assert(Self.Type == @TypeOf(other).Type); + } const Wt = @Vector(rows, T); - var result: Mat(T, @TypeOf(other).cols, Self.rows) = undefined; for (0..@TypeOf(result).cols) |i| { result.items[i] = @as(Wt, self.items[0]) * @as(Wt, @splat(other.items[i][0])); @@ -62,10 +62,6 @@ pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type return result; } - pub inline fn selfMul(self: *Self, other: anytype) void { - self.* = self.mul(other); - } - pub fn add(self: Self, other: Self) Self { var result = Self{ .items = undefined }; for (0..cols) |c| { @@ -76,8 +72,52 @@ pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type return result; } - pub inline fn selfAdd(self: *Self, other: Self) void { - self.* = self.add(other); + pub inline fn modify_row(self: Self, index: usize, v: anytype) Self { + comptime { + std.debug.assert(@typeInfo(@TypeOf(v)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(v)).vector.len <= cols); + } + var result = self; + for (0..@typeInfo(@TypeOf(v)).vector.len) |i| { + result.items[i][index] = v[i]; + } + return result; + } + + pub inline fn modify_column(self: Self, index: usize, v: anytype) Self { + comptime { + std.debug.assert(@typeInfo(@TypeOf(v)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(v)).vector.len <= rows); + } + + var result = self; + for (0..@typeInfo(@TypeOf(v)).vector.len) |i| { + result.items[index][i] = v[i]; + } + return result; + } + + pub inline fn column(self: Self, index: usize) @Vector(rows, T) { + return self.items[index]; + } + + pub inline fn row(self: Self, index: usize) @Vector(cols, T) { + var result: @Vector(cols, T) = undefined; + inline for (0..cols) |c| { + result[c] = self.items[c][index]; + } + return result; + } + + pub fn extract(self: Self, comptime sub_col: usize, comptime sub_row: usize) Mat(T, sub_col, sub_row) { + if (sub_col > cols or sub_row > rows) @compileError("sub matrix dimensions must be less than or equal to matrix dimensions"); + var result: Mat(T, sub_col, sub_row) = .from_column_major_array(undefined); + for (0..sub_col) |c| { + for (0..sub_row) |r| { + result.items[c][r] = self.items[c][r]; + } + } + return result; } pub fn sub(self: Self, other: Self) Self { @@ -90,155 +130,327 @@ pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type return result; } - pub inline fn selfSub(self: *Self, other: Self) void { - self.* = self.sub(other); - } - /// create a perspective projection matrix pub fn perspective(fovy: T, aspect: T, near: T, far: T) Self { - if (rows != 4 or cols != 4) @compileError("Perspective matrix must be 4x4"); + if (rows == 4 and cols == 4) { + const tanHalfFovy = std.math.tan(fovy / 2); - const tanHalfFovy = std.math.tan(fovy / 2); + var result: Self = .zero; + result.items[0][0] = 1.0 / (aspect * tanHalfFovy); + result.items[1][1] = 1.0 / tanHalfFovy; + result.items[2][2] = far / (near - far); + result.items[2][3] = -1.0; + result.items[3][2] = -(far * near) / (far - near); - var result: Self = .zero; - result.items[0][0] = 1.0 / (aspect * tanHalfFovy); - result.items[1][1] = 1.0 / tanHalfFovy; - result.items[2][2] = far / (near - far); - result.items[2][3] = -1.0; - result.items[3][2] = -(far * near) / (far - near); - - return result; + return result; + } + unreachable; } /// create a look-at view matrix pub fn lookAt(eye: @Vector(3, T), center: @Vector(3, T), up: @Vector(3, T)) Self { - if (rows != 4 or cols != 4) @compileError("Look-at matrix must be 4x4"); + if (rows == 4 and cols == 4) { + const f = vec.normalize(center - eye); + const s = vec.normalize(vec.cross(f, up)); + const u = vec.cross(s, f); - const f = vec.normalize(center - eye); - const s = vec.normalize(vec.cross(f, up)); - const u = vec.cross(s, f); - - var result: Self = .identity; - result.items[0][0] = s[0]; - result.items[1][0] = s[1]; - result.items[2][0] = s[2]; + var result: Self = .identity; + result.items[0][0] = s[0]; + result.items[1][0] = s[1]; + result.items[2][0] = s[2]; - result.items[0][1] = u[0]; - result.items[1][1] = u[1]; - result.items[2][1] = u[2]; + result.items[0][1] = u[0]; + result.items[1][1] = u[1]; + result.items[2][1] = u[2]; - result.items[0][2] = -f[0]; - result.items[1][2] = -f[1]; - result.items[2][2] = -f[2]; + result.items[0][2] = -f[0]; + result.items[1][2] = -f[1]; + result.items[2][2] = -f[2]; - result.items[3][0] = -vec.dot(s, eye); - result.items[3][1] = -vec.dot(u, eye); - result.items[3][2] = vec.dot(f, eye); + result.items[3][0] = -vec.dot(s, eye); + result.items[3][1] = -vec.dot(u, eye); + result.items[3][2] = vec.dot(f, eye); - return result; + return result; + } + unreachable; } pub fn translate(self: Self, vector: @Vector(rows - 1, T)) Self { - if (rows != cols) @compileError("Transform matrix must be square"); - var result = self; - result.items[cols - 1][0 .. rows - 1].* = self.items[cols - 1][0 .. rows - 1].* + vector; - return result; - } - - pub inline fn selfTranslate(self: *Self, vector: @Vector(rows - 1, T)) void { - self.* = self.translate(vector); + if (rows == cols) { + var result = self; + result.items[cols - 1][0 .. rows - 1].* = self.items[cols - 1][0 .. rows - 1].* + vector; + return result; + } + unreachable; } pub inline fn position(self: Self) @Vector(rows - 1, T) { - if (rows != cols) @compileError("Transform matrix must be square"); - return self.items[cols - 1][0 .. rows - 1].*; + if (rows == cols) { + return self.items[cols - 1][0 .. rows - 1].*; + } + unreachable; } /// Scaling transform matrix pub fn scale(self: Self, factors: @Vector(rows - 1, T)) Self { - if (rows != cols) @compileError("Transform matrix must be square"); + if (rows == cols) { + var result = self; + inline for (0..rows - 1) |i| { + result.items[i][0 .. rows - 1].* = self.items[i][0 .. rows - 1].* * @as(@Vector(rows - 1, T), @splat(factors[i])); + } + return result; + } + unreachable; + } - var result = self; - for (0..rows - 1) |i| { - result.items[i][0 .. cols - 1].* = self.items[i][0 .. cols - 1].* * factors; + pub fn rotate(self: Self, angle: T, axis: @Vector(3, T)) Self { + if (rows == 4 and cols == 4) { + const a = vec.normalize(axis); + const c = std.math.cos(angle); + const s = std.math.sin(angle); + const t = 1.0 - c; + + const rot: Self = .{ + .items = .{ + .{ t * a[0] * a[0] + c, t * a[0] * a[1] + s * a[2], t * a[0] * a[2] - s * a[1], 0 }, + .{ t * a[0] * a[1] - s * a[2], t * a[1] * a[1] + c, t * a[1] * a[2] + s * a[0], 0 }, + .{ t * a[0] * a[2] + s * a[1], t * a[1] * a[2] - s * a[0], t * a[2] * a[2] + c, 0 }, + .{ 0, 0, 0, 1 }, + }, + }; + + return self.mul(rot); } + unreachable; + } + /// The (cols-1)x(rows-1) sub matrix formed by deleting the given column and row. + pub fn minor(self: Self, comptime del_row: usize, comptime del_col: usize) Mat(T, cols - 1, rows - 1) { + var result: Mat(T, cols - 1, rows - 1) = undefined; + var cc: usize = 0; + inline for (0..cols) |c| { + if (c == del_col) continue; + var rr: usize = 0; + inline for (0..rows) |r| { + if (r == del_row) continue; + result.items[cc][rr] = self.items[c][r]; + rr += 1; + } + cc += 1; + } return result; } - pub inline fn selfScale(self: *Self, vector: @Vector(rows - 1, T)) void { - self.* = self.scale(vector); + /// The determinant of a square matrix (cofactor expansion along the first row). + pub fn determinant(self: Self) T { + comptime std.debug.assert(is_square); + if (rows == 1) return self.items[0][0]; + if (rows == 2) return self.items[0][0] * self.items[1][1] - self.items[1][0] * self.items[0][1]; + var det: T = 0; + inline for (0..cols) |c| { + const sign: T = if (c % 2 == 0) 1 else -1; + det += sign * self.items[c][0] * self.minor(0, c).determinant(); + } + return det; } - pub fn rotate(self: Self, angle: T, axis: @Vector(rows, T)) Self { - if (rows != cols) @compileError("Transform matrix must be square"); - if (rows != 4 or cols != 4) @compileError("unsuported dimensions, only suports 4x4"); + /// The inverse of a square matrix via the adjugate (cofactor) method. + /// Asserts the matrix is square; the determinant must be non-zero. + pub fn inverse(self: Self) Self { + comptime std.debug.assert(is_square and rows >= 2); + const inv_det = 1.0 / self.determinant(); + var result: Self = undefined; + // A^-1[row=i, col=j] = cofactor(j, i) / det (adjugate transposes the cofactor matrix) + inline for (0..rows) |i| { + inline for (0..cols) |j| { + const sign: T = if ((i + j) % 2 == 0) 1 else -1; + result.items[j][i] = sign * self.minor(j, i).determinant() * inv_det; + } + } + return result; + } - const a = vec.normalize(axis); - const c = std.math.cos(angle); - const s = std.math.sin(angle); - const t = 1.0 - c; + /// Fast inverse for a rigid transform (rotation + translation, no scale/shear): + /// transpose the 3x3 rotation and negate the rotated translation. 4x4 only. + pub fn inverse_ortho(self: Self) Self { + comptime std.debug.assert(rows == 4 and cols == 4); + var result: Self = .identity; + inline for (0..3) |c| { + inline for (0..3) |r| { + result.items[c][r] = self.items[r][c]; + } + } + const t = self.position(); + inline for (0..3) |r| { + var s: T = 0; + inline for (0..3) |k| { + s += self.items[r][k] * t[k]; + } + result.items[3][r] = -s; + } + return result; + } - const rot: Self = .{ - .items = .{ - .{ t * a[0] * a[0] + c, t * a[0] * a[1] + s * a[2], t * a[0] * a[2] - s * a[1], 0 }, - .{ t * a[0] * a[1] - s * a[2], t * a[1] * a[1] + c, t * a[1] * a[2] + s * a[0], 0 }, - .{ t * a[0] * a[2] + s * a[1], t * a[1] * a[2] - s * a[0], t * a[2] * a[2] + c, 0 }, - .{ 0, 0, 0, 1 }, - }, - }; + /// Build a rotation matrix from a quaternion (x, y, z, w). Works for 3x3 and 4x4. + pub fn from_quat(q: @Vector(4, T)) Self { + comptime std.debug.assert(is_square and (rows == 3 or rows == 4)); + const x = q[0]; + const y = q[1]; + const z = q[2]; + const w = q[3]; + var result: Self = .identity; + result.items[0][0] = 1 - 2 * (y * y + z * z); + result.items[0][1] = 2 * (x * y + w * z); + result.items[0][2] = 2 * (x * z - w * y); + result.items[1][0] = 2 * (x * y - w * z); + result.items[1][1] = 1 - 2 * (x * x + z * z); + result.items[1][2] = 2 * (y * z + w * x); + result.items[2][0] = 2 * (x * z + w * y); + result.items[2][1] = 2 * (y * z - w * x); + result.items[2][2] = 1 - 2 * (x * x + y * y); + return result; + } - return self.mul(rot); + /// Right-handed orthographic projection mapping depth to [0, 1] (matches `perspective`). 4x4 only. + pub fn orthographic(left: T, right: T, bottom: T, top: T, near: T, far: T) Self { + if (rows == 4 and cols == 4) { + var result: Self = .identity; + result.items[0][0] = 2.0 / (right - left); + result.items[1][1] = 2.0 / (top - bottom); + result.items[2][2] = -1.0 / (far - near); + result.items[3][0] = -(right + left) / (right - left); + result.items[3][1] = -(top + bottom) / (top - bottom); + result.items[3][2] = -near / (far - near); + return result; + } + unreachable; } - pub inline fn selfRotate(self: *Self, angle: T, axis: @Vector(3, T)) void { - self.* = self.rotate(angle, axis); + /// Right-handed perspective projection with an infinite far plane, depth range [0, 1]. 4x4 only. + pub fn perspective_infinite(fovy: T, aspect: T, near: T) Self { + if (rows == 4 and cols == 4) { + const tan_half_fovy = std.math.tan(fovy / 2); + var result: Self = .zero; + result.items[0][0] = 1.0 / (aspect * tan_half_fovy); + result.items[1][1] = 1.0 / tan_half_fovy; + result.items[2][2] = -1.0; + result.items[2][3] = -1.0; + result.items[3][2] = -near; + return result; + } + unreachable; } - pub const identity = blk: { - if (rows != cols) @compileError("Identity matrix must be square"); - var result: Self = .zero; - for (0..cols) |i| { - result.items[i][i] = 1; + /// Right-handed reversed-Z perspective (near -> 1, far -> 0), depth range [0, 1]. + /// Reversed-Z maximizes floating-point depth precision. 4x4 only. + pub fn perspective_reverse_z(fovy: T, aspect: T, near: T, far: T) Self { + if (rows == 4 and cols == 4) { + const tan_half_fovy = std.math.tan(fovy / 2); + var result: Self = .zero; + result.items[0][0] = 1.0 / (aspect * tan_half_fovy); + result.items[1][1] = 1.0 / tan_half_fovy; + result.items[2][2] = near / (far - near); + result.items[2][3] = -1.0; + result.items[3][2] = (far * near) / (far - near); + return result; + } + unreachable; + } + + /// Per-axis scale factors, taken from the lengths of the first three basis columns. + pub fn get_scale(self: Self) @Vector(3, T) { + comptime std.debug.assert(rows >= 3 and cols >= 3); + var s: @Vector(3, T) = undefined; + inline for (0..3) |i| { + s[i] = vec.norm(@Vector(3, T){ self.items[i][0], self.items[i][1], self.items[i][2] }); + } + return s; + } + + pub const Decomposed = struct { + translation: @Vector(3, T), + rotation: @Vector(4, T), // quaternion (x, y, z, w) + scale: @Vector(3, T), + }; + + /// Decompose an affine transform into translation, rotation (quaternion) and scale. + /// Assumes no shear and non-negative scale. 4x4 only. + pub fn decompose(self: Self) Decomposed { + comptime std.debug.assert(rows == 4 and cols == 4); + const s = self.get_scale(); + var rot = self; + inline for (0..3) |i| { + const inv = 1.0 / s[i]; + inline for (0..3) |r| { + rot.items[i][r] = self.items[i][r] * inv; + } } - break :blk result; + return .{ + .translation = self.position(), + .rotation = @import("quat.zig").from_matrix(rot), + .scale = s, + }; + } + + pub const Projection = struct { + near: T, + far: T, + fovy: T, + aspect: T, }; - pub const zero: Self = .fromCM(@splat(@splat(0))); + /// Recover the parameters of a right-handed [0, 1]-depth perspective projection + /// (as produced by `perspective`). 4x4 only. + pub fn decompose_projection(self: Self) Projection { + comptime std.debug.assert(rows == 4 and cols == 4); + const tan_half = 1.0 / self.items[1][1]; + const a = self.items[2][2]; + const b = self.items[3][2]; + const near = b / a; + return .{ + .near = near, + .far = a * near / (1.0 + a), + .fovy = 2.0 * std.math.atan(tan_half), + .aspect = self.items[1][1] / self.items[0][0], + }; + } + + pub const identity = blk: { + if (rows == cols) { + var result: Self = .zero; + for (0..cols) |i| { + result.items[i][i] = 1; + } + break :blk result; + } + unreachable; + }; - // TODO: doc - pub fn format( - self: Self, - comptime fmt: []const u8, - options: std.fmt.FormatOptions, - writer: anytype, - ) !void { - const separator_indx = comptime std.mem.lastIndexOfScalar(u8, fmt, '|'); - const scalar_fmt = if (separator_indx) |i| fmt[i + 1 ..] else fmt; - const grid_fmt = if (separator_indx) |i| fmt[0..i] else "g"; + pub const zero: Self = .from_column_major_array(@splat(@splat(0))); - if (comptime std.mem.eql(u8, grid_fmt, "l")) { - try writer.writeAll("["); + pub fn format(self: @This(), writer: *std.Io.Writer) std.Io.Writer.Error!void { + var max_widths: [cols]usize = @splat(0); + for (0..cols) |c| { for (0..rows) |r| { - try writer.writeAll("["); - for (0..cols) |c| { - try std.fmt.formatType(self.items[c][r], scalar_fmt, options, writer, 0); - if (c < cols - 1) try writer.writeAll(", "); - } - if (r < rows - 1) try writer.writeAll(", "); - try writer.writeAll("]"); + const len = std.fmt.count("{d}", .{self.items[c][r]}); + max_widths[c] = @max(max_widths[c], len); } - try writer.writeAll("]"); - } else if (comptime std.mem.eql(u8, grid_fmt, "g")) { - for (0..rows) |r| { - try writer.writeAll("["); - for (0..cols) |c| { - try std.fmt.formatType(self.items[c][r], scalar_fmt, options, writer, 0); - if (c < cols - 1) try writer.writeAll(", "); + } + + for (0..rows) |r| { + try writer.writeAll("["); + for (0..cols) |c| { + const len = std.fmt.count("{d}", .{self.items[c][r]}); + for (0..max_widths[c] - len) |_| { + try writer.writeByte(' '); } - try writer.writeAll("]\n"); + try writer.print("{d}", .{self.items[c][r]}); + if (c < cols - 1) try writer.writeAll(", "); } - } else @compileError("invalid matrix fmt specificer, specificer must be `l` or `g` but is `" ++ grid_fmt ++ "`"); + try writer.writeByte(']'); + if (r != rows - 1) try writer.writeByte('\n'); + } } pub fn eql(self: Self, other: Self) bool { @@ -254,30 +466,63 @@ pub fn Mat(comptime T: type, comptime cols_: usize, comptime rows_: usize) type }; } +test "format" { + const c: Mat(f32, 3, 3) = .from_row_major_array(.{ + .{ 9, 12, 15 }, + .{ 19, 26, 33 }, + .{ 29, 40, 51 }, + }); + var buff: [128]u8 = undefined; + const result = try std.fmt.bufPrint(&buff, "{f}", .{c}); + try std.testing.expectEqualStrings( + \\[ 9, 12, 15] + \\[19, 26, 33] + \\[29, 40, 51] + , result); +} + test "translate" { const c = Mat(f32, 4, 4).identity.translate(.{ 1, 2, 3 }); - const excpected_c: Mat(f32, 4, 4) = .fromRM(.{ + const expected: Mat(f32, 4, 4) = .from_row_major_array(.{ .{ 1, 0, 0, 1 }, .{ 0, 1, 0, 2 }, .{ 0, 0, 1, 3 }, .{ 0, 0, 0, 1 }, }); - try std.testing.expectEqual(excpected_c, c); + try std.testing.expectEqual(expected, c); +} + +test "modify_column" { + var m = Mat(f32, 4, 4).zero; + m = m.modify_column(0, @Vector(4, f32){ 1, 2, 3, 4 }); + try std.testing.expectEqual(m.column(0), .{ 1, 2, 3, 4 }); + m = m.modify_column(0, @Vector(3, f32){ 5, 6, 7 }); + try std.testing.expectEqual(m.column(0), .{ 5, 6, 7, 4 }); + + m = m.modify_column(0, @Vector(4, f32){ 8, 9, 10, 0 }); + m = m.modify_column(1, @Vector(4, f32){ 11, 12, 13, 0 }); + m = m.modify_column(2, @Vector(4, f32){ 14, 15, 16, 0 }); + m = m.modify_column(3, @Vector(4, f32){ 17, 18, 19, 0 }); + + try std.testing.expectEqual(m.column(0), .{ 8, 9, 10, 0 }); + try std.testing.expectEqual(m.column(1), .{ 11, 12, 13, 0 }); + try std.testing.expectEqual(m.column(2), .{ 14, 15, 16, 0 }); + try std.testing.expectEqual(m.position(), .{ 17, 18, 19 }); } test "mul" { { - const a = Mat(f32, 2, 2).fromCM(.{ + const a = Mat(f32, 2, 2).from_column_major_array(.{ .{ 1, 2 }, .{ 3, 4 }, }); - const b = Mat(f32, 2, 2).fromCM(.{ + const b = Mat(f32, 2, 2).from_column_major_array(.{ .{ 5, 6 }, .{ 7, 8 }, }); const c = a.mul(b); - const excpected_c = Mat(f32, 2, 2).fromCM(.{ + const excpected_c = Mat(f32, 2, 2).from_column_major_array(.{ .{ 23, 34 }, .{ 31, 46 }, }); @@ -285,18 +530,18 @@ test "mul" { } { - const a = Mat(f32, 2, 3).fromCM(.{ + const a = Mat(f32, 2, 3).from_column_major_array(.{ .{ 1, 2, 3 }, .{ 4, 5, 6 }, }); - const b = Mat(f32, 3, 2).fromCM(.{ + const b = Mat(f32, 3, 2).from_column_major_array(.{ .{ 1, 2 }, .{ 3, 4 }, .{ 5, 6 }, }); const c = a.mul(b); - const excpected_c = Mat(f32, 3, 3).fromCM(.{ + const excpected_c = Mat(f32, 3, 3).from_column_major_array(.{ .{ 9, 12, 15 }, .{ 19, 26, 33 }, .{ 29, 40, 51 }, @@ -305,13 +550,13 @@ test "mul" { } { - const a = Mat(f32, 4, 4).fromCM(.{ + const a = Mat(f32, 4, 4).from_column_major_array(.{ .{ 1, 2, 3, 4 }, .{ 5, 6, 7, 8 }, .{ 9, 10, 11, 12 }, .{ 13, 14, 15, 16 }, }); - const b = Mat(f32, 4, 4).fromCM(.{ + const b = Mat(f32, 4, 4).from_column_major_array(.{ .{ 17, 18, 19, 20 }, .{ 21, 22, 23, 24 }, .{ 25, 26, 27, 28 }, @@ -319,7 +564,7 @@ test "mul" { }); const c = a.mul(b); - const excpected_c = Mat(f32, 4, 4).fromCM(.{ + const excpected_c = Mat(f32, 4, 4).from_column_major_array(.{ .{ 538, 612, 686, 760 }, .{ 650, 740, 830, 920 }, .{ 762, 868, 974, 1080 }, @@ -328,3 +573,186 @@ test "mul" { try std.testing.expectEqual(excpected_c, c); } } + +test "scale" { + { + const mat: Mat(f32, 4, 4) = .from_row_major_array(.{ + .{ 1, 0, 0, 5 }, + .{ 0, 1, 0, 6 }, + .{ 0, 0, 1, 7 }, + .{ 0, 0, 0, 1 }, + }); + const scaled = mat.scale(.{ 2, 3, 4 }); + + const expected: Mat(f32, 4, 4) = .from_row_major_array(.{ + .{ 2, 0, 0, 5 }, + .{ 0, 3, 0, 6 }, + .{ 0, 0, 4, 7 }, + .{ 0, 0, 0, 1 }, + }); + try std.testing.expectEqual(expected, scaled); + } + + { + const mat: Mat(f32, 3, 3) = .from_row_major_array(.{ + .{ 0.707, -0.707, 0 }, + .{ 0.707, 0.707, 0 }, + .{ 0, 0, 1 }, + }); + + const scaled = mat.scale(.{ 2, 3 }); + + const expected: Mat(f32, 3, 3) = .from_row_major_array(.{ + .{ 1.414, -2.121, 0 }, + .{ 1.414, 2.121, 0 }, + .{ 0, 0, 1 }, + }); + + for (0..3) |c| { + for (0..3) |r| { + try std.testing.expectApproxEqAbs(expected.items[c][r], scaled.items[c][r], 0.001); + } + } + } +} + +fn mul_vec4(m: Mat(f32, 4, 4), v: @Vector(4, f32)) @Vector(4, f32) { + var out: @Vector(4, f32) = .{ 0, 0, 0, 0 }; + inline for (0..4) |c| { + out += @as(@Vector(4, f32), m.items[c]) * @as(@Vector(4, f32), @splat(v[c])); + } + return out; +} + +fn expect_mat_close(comptime C: usize, comptime R: usize, expected: Mat(f32, C, R), actual: Mat(f32, C, R), tol: f32) !void { + inline for (0..C) |c| { + inline for (0..R) |r| { + try std.testing.expectApproxEqAbs(expected.items[c][r], actual.items[c][r], tol); + } + } +} + +test "determinant" { + try std.testing.expectEqual(@as(f32, 1), Mat(f32, 4, 4).identity.determinant()); + + const m2: Mat(f32, 2, 2) = .from_row_major_array(.{ + .{ 1, 2 }, + .{ 3, 4 }, + }); + try std.testing.expectApproxEqAbs(@as(f32, -2), m2.determinant(), 1e-6); + + const m3: Mat(f32, 3, 3) = .from_row_major_array(.{ + .{ 2, 0, 1 }, + .{ 1, 3, 2 }, + .{ 1, 0, 1 }, + }); + try std.testing.expectApproxEqAbs(@as(f32, 3), m3.determinant(), 1e-6); +} + +test "inverse" { + const m: Mat(f32, 4, 4) = .from_row_major_array(.{ + .{ 4, 1, 0, 0 }, + .{ 1, 3, 1, 0 }, + .{ 0, 1, 2, 1 }, + .{ 0, 0, 1, 2 }, + }); + try expect_mat_close(4, 4, .identity, m.mul(m.inverse()), 1e-4); + + const m3: Mat(f32, 3, 3) = .from_row_major_array(.{ + .{ 2, 0, 1 }, + .{ 1, 3, 2 }, + .{ 1, 0, 1 }, + }); + try expect_mat_close(3, 3, .identity, m3.mul(m3.inverse()), 1e-4); +} + +test "inverse_ortho" { + const quat = @import("quat.zig"); + const q = quat.from_rotation(@Vector(3, f32){ 0, 1, 0 }, 0.7); + var m = Mat(f32, 4, 4).from_quat(q); + m = m.translate(.{ 3, -2, 5 }); + + try expect_mat_close(4, 4, .identity, m.mul(m.inverse_ortho()), 1e-5); + // For a rigid transform inverse_ortho must agree with the general inverse. + try expect_mat_close(4, 4, m.inverse(), m.inverse_ortho(), 1e-4); +} + +test "from_quat" { + const quat = @import("quat.zig"); + const axis = @Vector(3, f32){ 0, 0, 1 }; + const angle: f32 = std.math.pi / 3.0; + const q = quat.from_rotation(axis, angle); + // from_quat must match the equivalent axis-angle rotation matrix. + try expect_mat_close(4, 4, Mat(f32, 4, 4).identity.rotate(angle, axis), Mat(f32, 4, 4).from_quat(q), 1e-5); +} + +test "orthographic" { + const m = Mat(f32, 4, 4).orthographic(-2, 2, -1, 1, 0.5, 10); + // near corner (left, bottom, -near) -> NDC (-1, -1, 0) + const near_corner = mul_vec4(m, .{ -2, -1, -0.5, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, -1), near_corner[0], 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, -1), near_corner[1], 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 0), near_corner[2], 1e-5); + // far plane -> z_ndc 1 + const far_pt = mul_vec4(m, .{ 0, 0, -10, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, 1), far_pt[2], 1e-5); +} + +test "perspective_reverse_z" { + const near: f32 = 0.1; + const far: f32 = 100.0; + const m = Mat(f32, 4, 4).perspective_reverse_z(std.math.pi / 2.0, 1.0, near, far); + const at_near = mul_vec4(m, .{ 0, 0, -near, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, 1), at_near[2] / at_near[3], 1e-4); + const at_far = mul_vec4(m, .{ 0, 0, -far, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, 0), at_far[2] / at_far[3], 1e-4); +} + +test "perspective_infinite" { + const near: f32 = 0.1; + const m = Mat(f32, 4, 4).perspective_infinite(std.math.pi / 2.0, 1.0, near); + const at_near = mul_vec4(m, .{ 0, 0, -near, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, 0), at_near[2] / at_near[3], 1e-4); + const far_away = mul_vec4(m, .{ 0, 0, -1.0e9, 1 }); + try std.testing.expectApproxEqAbs(@as(f32, 1), far_away[2] / far_away[3], 1e-3); +} + +test "get_scale" { + const s = Mat(f32, 4, 4).identity.scale(.{ 2, 3, 4 }).get_scale(); + try std.testing.expectApproxEqAbs(@as(f32, 2), s[0], 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 3), s[1], 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 4), s[2], 1e-5); +} + +test "decompose" { + const quat = @import("quat.zig"); + const q = quat.from_rotation(vec.normalize(@Vector(3, f32){ 1, 1, 0 }), 0.5); + var m = Mat(f32, 4, 4).from_quat(q); + m = m.scale(.{ 2, 3, 4 }); + m = m.translate(.{ 4, -1, 7 }); + + const d = m.decompose(); + try std.testing.expectApproxEqAbs(@as(f32, 2), d.scale[0], 1e-4); + try std.testing.expectApproxEqAbs(@as(f32, 3), d.scale[1], 1e-4); + try std.testing.expectApproxEqAbs(@as(f32, 4), d.scale[2], 1e-4); + try std.testing.expectApproxEqAbs(@as(f32, 4), d.translation[0], 1e-4); + try std.testing.expectApproxEqAbs(@as(f32, 7), d.translation[2], 1e-4); + + // Recomposing must reproduce the original matrix. + var m2 = Mat(f32, 4, 4).from_quat(d.rotation); + m2 = m2.scale(d.scale); + m2 = m2.translate(d.translation); + try expect_mat_close(4, 4, m, m2, 1e-3); +} + +test "decompose_projection" { + const fovy: f32 = std.math.pi / 3.0; + const aspect: f32 = 1.6; + const near: f32 = 0.2; + const far: f32 = 75.0; + const d = Mat(f32, 4, 4).perspective(fovy, aspect, near, far).decompose_projection(); + try std.testing.expectApproxEqAbs(near, d.near, 1e-3); + try std.testing.expectApproxEqAbs(far, d.far, 5e-2); + try std.testing.expectApproxEqAbs(fovy, d.fovy, 1e-4); + try std.testing.expectApproxEqAbs(aspect, d.aspect, 1e-4); +} diff --git a/src/meta.zig b/src/meta.zig new file mode 100644 index 0000000..b20f916 --- /dev/null +++ b/src/meta.zig @@ -0,0 +1,65 @@ +const std = @import("std"); + +/// Element/scalar type of a vector or array type. +pub fn Child(comptime T: type) type { + return switch (@typeInfo(T)) { + .vector, .array => std.meta.Child(T), + else => @compileError("Expected vector or array type, got: " ++ @typeName(T)), + }; +} + +/// Number of elements in a vector or array type. +pub fn lengthOf(comptime T: type) comptime_int { + return switch (@typeInfo(T)) { + .vector => |v| v.len, + .array => |a| a.len, + else => @compileError("Expected vector or array type, got: " ++ @typeName(T)), + }; +} + +/// The `@Vector` type equivalent to a vector-or-array type. Used to coerce inputs +/// into a SIMD vector for the actual math. +pub fn AsVector(comptime T: type) type { + return @Vector(lengthOf(T), Child(T)); +} + +/// A container of the SAME kind (array vs vector) as `T`, but with `len` elements. +/// Used so a function's return shape matches its input shape. +pub fn Reshape(comptime T: type, comptime len: usize) type { + return switch (@typeInfo(T)) { + .vector => @Vector(len, Child(T)), + .array => [len]Child(T), + else => @compileError("Expected vector or array type, got: " ++ @typeName(T)), + }; +} + +test Child { + try std.testing.expect(Child(@Vector(3, f32)) == f32); + try std.testing.expect(Child([3]f32) == f32); + try std.testing.expect(Child([4]i32) == i32); +} + +test lengthOf { + try std.testing.expectEqual(3, lengthOf(@Vector(3, f32))); + try std.testing.expectEqual(4, lengthOf([4]f64)); +} + +test AsVector { + try std.testing.expect(AsVector(@Vector(3, f32)) == @Vector(3, f32)); + try std.testing.expect(AsVector([3]f32) == @Vector(3, f32)); + // An array coerces to its equivalent @Vector. + const arr = [3]f32{ 1, 2, 3 }; + const v: AsVector(@TypeOf(arr)) = arr; + try std.testing.expectEqual(@as(f32, 6), @reduce(.Add, v)); +} + +test Reshape { + // Same kind is preserved, length can change. + try std.testing.expect(Reshape(@Vector(4, f32), 2) == @Vector(2, f32)); + try std.testing.expect(Reshape([4]f32, 2) == [2]f32); + try std.testing.expect(Reshape([3]i32, 3) == [3]i32); + // A @Vector coerces back to its equivalent array. + const v = @Vector(3, f32){ 1, 2, 3 }; + const arr: Reshape([3]f32, 3) = v; + try std.testing.expectEqual(@as(f32, 2), arr[1]); +} diff --git a/src/packing.zig b/src/packing.zig new file mode 100644 index 0000000..b877649 --- /dev/null +++ b/src/packing.zig @@ -0,0 +1,145 @@ +const std = @import("std"); + +// --------------------------------------------------------------------------- +// Bit reinterpretation +// --------------------------------------------------------------------------- + +pub fn asuint(x: f32) u32 { + return @bitCast(x); +} + +pub fn asfloat(x: u32) f32 { + return @bitCast(x); +} + +// --------------------------------------------------------------------------- +// UNORM / SNORM (fixed-point normalized integers of arbitrary bit width) +// --------------------------------------------------------------------------- + +/// Pack a value in [0, 1] into a `bits`-wide unsigned normalized integer. +pub fn pack_unorm(comptime bits: u16, x: f32) u32 { + const maxv: f32 = @floatFromInt((@as(u64, 1) << bits) - 1); + return @intFromFloat(@round(std.math.clamp(x, 0, 1) * maxv)); +} + +/// Unpack a `bits`-wide unsigned normalized integer back into [0, 1]. +pub fn unpack_unorm(comptime bits: u16, v: u32) f32 { + const maxv: f32 = @floatFromInt((@as(u64, 1) << bits) - 1); + return @as(f32, @floatFromInt(v)) / maxv; +} + +/// Pack a value in [-1, 1] into a `bits`-wide signed normalized integer. +pub fn pack_snorm(comptime bits: u16, x: f32) i32 { + const maxv: f32 = @floatFromInt((@as(u64, 1) << (bits - 1)) - 1); + return @intFromFloat(@round(std.math.clamp(x, -1, 1) * maxv)); +} + +/// Unpack a `bits`-wide signed normalized integer back into [-1, 1]. +pub fn unpack_snorm(comptime bits: u16, v: i32) f32 { + const maxv: f32 = @floatFromInt((@as(u64, 1) << (bits - 1)) - 1); + return @max(@as(f32, @floatFromInt(v)) / maxv, -1.0); +} + +// --------------------------------------------------------------------------- +// Half precision (IEEE binary16). Zig has a native f16, so these are bit casts. +// --------------------------------------------------------------------------- + +pub fn pack_f16(x: f32) u16 { + return @bitCast(@as(f16, @floatCast(x))); +} + +pub fn unpack_f16(v: u16) f32 { + return @floatCast(@as(f16, @bitCast(v))); +} + +// --------------------------------------------------------------------------- +// Packed float formats derived from the binary16 layout (5-bit exponent, bias 15) +// --------------------------------------------------------------------------- + +/// R11G11B10F: 11-bit (5e6m) R and G, 10-bit (5e5m) B, no sign. Negatives clamp to 0. +/// Mantissa bits are truncated, so this is a lossy pack. +pub fn pack_ufloat_11_11_10(rgb: @Vector(3, f32)) u32 { + const r: u32 = f32_to_ufloat(rgb[0], 6); + const g: u32 = f32_to_ufloat(rgb[1], 6); + const b: u32 = f32_to_ufloat(rgb[2], 5); + return r | (g << 11) | (b << 22); +} + +pub fn unpack_ufloat_11_11_10(v: u32) @Vector(3, f32) { + return .{ + ufloat_to_f32(v & 0x7FF, 6), + ufloat_to_f32((v >> 11) & 0x7FF, 6), + ufloat_to_f32((v >> 22) & 0x3FF, 5), + }; +} + +fn f32_to_ufloat(x: f32, comptime mantissa_bits: u4) u32 { + if (!(x > 0)) return 0; // <= 0 or NaN + const h: u16 = @bitCast(@as(f16, @floatCast(x))); + const exp: u32 = (h >> 10) & 0x1F; + const mant: u32 = (h >> (10 - mantissa_bits)) & ((@as(u32, 1) << mantissa_bits) - 1); + return (exp << mantissa_bits) | mant; +} + +fn ufloat_to_f32(v: u32, comptime mantissa_bits: u4) f32 { + const exp: u16 = @intCast((v >> mantissa_bits) & 0x1F); + const mant: u16 = @intCast(v & ((@as(u32, 1) << mantissa_bits) - 1)); + const h: u16 = (exp << 10) | (mant << (10 - mantissa_bits)); + return @floatCast(@as(f16, @bitCast(h))); +} + +// --------------------------------------------------------------------------- +// FP8 (OCP E5M2): 1 sign, 5 exponent, 2 mantissa โ€” the top 8 bits of binary16. +// --------------------------------------------------------------------------- + +pub fn pack_f8_e5m2(x: f32) u8 { + const h: u16 = @bitCast(@as(f16, @floatCast(x))); + return @intCast(h >> 8); +} + +pub fn unpack_f8_e5m2(v: u8) f32 { + return @floatCast(@as(f16, @bitCast(@as(u16, v) << 8))); +} + +test "unorm/snorm round trip" { + for ([_]f32{ 0, 0.25, 0.5, 0.75, 1 }) |x| { + try std.testing.expectApproxEqAbs(x, unpack_unorm(8, pack_unorm(8, x)), 1.0 / 255.0); + try std.testing.expectApproxEqAbs(x, unpack_unorm(16, pack_unorm(16, x)), 1.0 / 65535.0); + } + for ([_]f32{ -1, -0.5, 0, 0.5, 1 }) |x| { + try std.testing.expectApproxEqAbs(x, unpack_snorm(8, pack_snorm(8, x)), 1.0 / 127.0); + try std.testing.expectApproxEqAbs(x, unpack_snorm(16, pack_snorm(16, x)), 1.0 / 32767.0); + } + // Out-of-range values saturate. + try std.testing.expectEqual(@as(f32, 1), unpack_unorm(8, pack_unorm(8, 5))); + try std.testing.expectEqual(@as(f32, 0), unpack_unorm(8, pack_unorm(8, -5))); +} + +test "f16 round trip" { + for ([_]f32{ 0, 1, -2, 0.5, 100.25 }) |x| { + try std.testing.expectEqual(x, unpack_f16(pack_f16(x))); + } +} + +test "asfloat/asuint" { + try std.testing.expectEqual(@as(f32, 1.5), asfloat(asuint(1.5))); + try std.testing.expectEqual(@as(u32, 0x3F800000), asuint(1.0)); +} + +test "ufloat_11_11_10 round trip" { + // Values with exact low-order mantissa survive the pack. + const rgb = @Vector(3, f32){ 1.0, 2.0, 0.5 }; + const back = unpack_ufloat_11_11_10(pack_ufloat_11_11_10(rgb)); + try std.testing.expectApproxEqAbs(rgb[0], back[0], 1e-6); + try std.testing.expectApproxEqAbs(rgb[1], back[1], 1e-6); + try std.testing.expectApproxEqAbs(rgb[2], back[2], 1e-6); + // Negatives clamp to zero. + const neg = unpack_ufloat_11_11_10(pack_ufloat_11_11_10(.{ -1, 0, 0 })); + try std.testing.expectEqual(@as(f32, 0), neg[0]); +} + +test "f8_e5m2 round trip" { + for ([_]f32{ 0, 0.5, 1, 1.5, 2, -2 }) |x| { + try std.testing.expectEqual(x, unpack_f8_e5m2(pack_f8_e5m2(x))); + } +} diff --git a/src/quat.zig b/src/quat.zig new file mode 100644 index 0000000..3338662 --- /dev/null +++ b/src/quat.zig @@ -0,0 +1,389 @@ +const std = @import("std"); +const vector = @import("vector.zig"); +const meta = @import("meta.zig"); +const zml = @import("root.zig"); + +pub const Quat4f32 = @Vector(4, f32); +pub const Quat4f64 = @Vector(4, f64); + +fn map_to_vector(a: anytype) meta.AsVector(@TypeOf(a)) { + return a; +} + +fn Quat(comptime T: type) type { + return @Vector(4, T); +} + +pub fn identity(comptime T: type) @Vector(4, T) { + return .{ 0, 0, 0, 1 }; +} + +pub fn zero(comptime T: type) @Vector(4, T) { + return .{ 0, 0, 0, 0 }; +} + +pub fn x_axis(comptime T: type) @Vector(3, T) { + return .{ 1, 0, 0 }; +} + +pub fn y_axis(comptime T: type) @Vector(3, T) { + return .{ 0, 1, 0 }; +} + +pub fn z_axis(comptime T: type) @Vector(3, T) { + return .{ 0, 0, 1 }; +} + +// Create quaternion that rotates a vector from the direction of inFrom to the direction of inTo along the shortest path +// @see https://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm +pub fn from_to(from: anytype, to: @TypeOf(from)) @Vector(4, std.meta.Child(@TypeOf(from))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(from)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(from)).vector.len == 3); + } + + //Uses (inFrom = v1, inTo = v2): + + //angle = arcos(v1 . v2 / |v1||v2|) + //axis = normalize(v1 x v2) + + //Quaternion is then: + + //s = sin(angle / 2) + //x = axis.x * s + //y = axis.y * s + //z = axis.z * s + //w = cos(angle / 2) + + //Using identities: + + //sin(2 * a) = 2 * sin(a) * cos(a) + //cos(2 * a) = cos(a)^2 - sin(a)^2 + //sin(a)^2 + cos(a)^2 = 1 + + //This reduces to: + + //x = (v1 x v2).x + //y = (v1 x v2).y + //z = (v1 x v2).z + //w = |v1||v2| + v1 . v2 + + //which then needs to be normalized because the whole equation was multiplied by 2 cos(angle / 2) + const len_v1_v2 = std.math.sqrt(vector.norm_sqr(from) * vector.norm_sqr(to)); + const w = len_v1_v2 + vector.dot(from, to); + if (w == 0.0) { + if (len_v1_v2 == 0.0) { + return identity(std.meta.Child(@TypeOf(from))); + } else { + const norm_perp = vector.norm_perpendicular(from); + return .{ norm_perp[0], norm_perp[1], norm_perp[2], 0 }; + } + } + const v = vector.cross(from, to); + return vector.normalize(@Vector(4, std.meta.Child(@TypeOf(from))){ v[0], v[1], v[2], w }); +} + +pub fn get_twist(inAxis: anytype) @Vector(4, std.meta.Child(@TypeOf(inAxis))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(inAxis)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(inAxis)).vector.len == 4); + } + + const dir = vector.dot(vector.extract(inAxis, 3), inAxis) * inAxis; + const twist: @Vector(4, std.meta.Child(@TypeOf(inAxis))) = .{ dir[0], dir[1], dir[2], inAxis[3] }; + const twist_len = vector.norm_sqr(twist); + if (twist_len == 0.0) { + return twist / @as(@Vector(4, std.meta.Child(@TypeOf(inAxis))), @splat(std.math.sqrt(twist_len))); + } + return identity(std.meta.Child(@TypeOf(inAxis))); +} + +//TODO: optimize with simd +pub fn mul(a: anytype, b: @TypeOf(a)) @TypeOf(a) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(a)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(a)).vector.len == 4); + } + const inner_a = map_to_vector(a); + const inner_b = map_to_vector(b); + + const lx = inner_a[0]; + const ly = inner_a[1]; + const lz = inner_a[2]; + const lw = inner_a[3]; + + const rx = inner_b[0]; + const ry = inner_b[1]; + const rz = inner_b[2]; + const rw = inner_b[3]; + + return @TypeOf(a){ lw * rx + lx * rw + ly * rz - lz * ry, lw * ry - lx * rz + ly * rw + lz * rx, lw * rz + lx * ry - ly * rx + lz * rw, lw * rw - lx * rx - ly * ry - lz * rz }; +} + +pub fn norm(q: anytype) std.meta.Float(@bitSizeOf(std.meta.Child(@TypeOf(q)))) { + return vector.norm(q); +} + +pub fn to_axis_angle(q: anytype) struct { + axis: @Vector(3, std.meta.Child(@TypeOf(q))), + angle: std.meta.Child(@TypeOf(q)), +} { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + const C = std.meta.Child(@TypeOf(q)); + const qw_clamped = std.math.clamp(q[3], -1.0, 1.0); + const angle = 2.0 * std.math.acos(qw_clamped); + const s = std.math.sqrt(1.0 - qw_clamped * qw_clamped); + if (s < 0.001) { // axis direction is arbitrary when the angle is close to zero + return .{ .axis = .{ 1, 0, 0 }, .angle = angle }; + } + return .{ .axis = @Vector(3, C){ q[0] / s, q[1] / s, q[2] / s }, .angle = angle }; +} + +/// Rotate a 3D vector by a (unit) quaternion. Uses the branchless +/// t = 2ยท(u ร— v); v' = v + wยทt + u ร— t form (u = q.xyz). +pub fn rotate_vector(q: anytype, v: @Vector(3, std.meta.Child(@TypeOf(q)))) @Vector(3, std.meta.Child(@TypeOf(q))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + const C = std.meta.Child(@TypeOf(q)); + const u = @Vector(3, C){ q[0], q[1], q[2] }; + const t = @as(@Vector(3, C), @splat(2)) * vector.cross(u, v); + return v + @as(@Vector(3, C), @splat(q[3])) * t + vector.cross(u, t); +} + +/// Rotate a 3D vector by the inverse of a (unit) quaternion. +pub fn rotate_vector_inv(q: anytype, v: @Vector(3, std.meta.Child(@TypeOf(q)))) @Vector(3, std.meta.Child(@TypeOf(q))) { + return rotate_vector(conjugate(q), v); +} + +/// Convert a quaternion (x, y, z, w) to a 4x4 rotation matrix. +pub fn to_matrix(q: anytype) zml.Mat(std.meta.Child(@TypeOf(q)), 4, 4) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + return zml.Mat(std.meta.Child(@TypeOf(q)), 4, 4).from_quat(q); +} + +/// Extract a quaternion from the rotation part of a matrix (Shepperd's method). +pub fn from_matrix(m: anytype) @Vector(4, @TypeOf(m).Type) { + // element(row, col) is stored column-major as items[col][row] + const m00 = m.items[0][0]; + const m01 = m.items[1][0]; + const m02 = m.items[2][0]; + const m10 = m.items[0][1]; + const m11 = m.items[1][1]; + const m12 = m.items[2][1]; + const m20 = m.items[0][2]; + const m21 = m.items[1][2]; + const m22 = m.items[2][2]; + const trace = m00 + m11 + m22; + if (trace > 0) { + const s = @sqrt(trace + 1.0) * 2.0; // s = 4 * qw + return .{ (m21 - m12) / s, (m02 - m20) / s, (m10 - m01) / s, 0.25 * s }; + } else if (m00 > m11 and m00 > m22) { + const s = @sqrt(1.0 + m00 - m11 - m22) * 2.0; // s = 4 * qx + return .{ 0.25 * s, (m01 + m10) / s, (m02 + m20) / s, (m21 - m12) / s }; + } else if (m11 > m22) { + const s = @sqrt(1.0 + m11 - m00 - m22) * 2.0; // s = 4 * qy + return .{ (m01 + m10) / s, 0.25 * s, (m12 + m21) / s, (m02 - m20) / s }; + } else { + const s = @sqrt(1.0 + m22 - m00 - m11) * 2.0; // s = 4 * qz + return .{ (m02 + m20) / s, (m12 + m21) / s, 0.25 * s, (m10 - m01) / s }; + } +} + +pub fn conjugate(q: anytype) @TypeOf(q) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + return q * @Vector(4, std.meta.Child(@TypeOf(q))){ -1, -1, -1, 1 }; +} + +pub fn inverse(q: anytype) @TypeOf(q) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + return conjugate(q) / @as(@TypeOf(q), @splat(vector.norm(q))); +} + +pub fn slerp(a: anytype, b: anytype, factor: std.meta.Child(@TypeOf(a))) Quat(std.meta.Child(@TypeOf(a))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(a)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(a)).vector.len == 4); + std.debug.assert(@typeInfo(@TypeOf(b)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(b)).vector.len == 4); + std.debug.assert(@TypeOf(a) == @TypeOf(b)); + } + const inner_a = map_to_vector(a); + const inner_b = map_to_vector(b); + const delta: std.meta.Child(@TypeOf(a)) = 0.0001; + + var sign_scale1: std.meta.Child(@TypeOf(a)) = 1.0; + var cos_omega = vector.dot(inner_a, inner_b); + + if (cos_omega < 0.0) { + cos_omega = -cos_omega; + sign_scale1 = -1.0; + } + + // Calculate coefficients + var scale0: std.meta.Child(@TypeOf(a)) = undefined; + var scale1: std.meta.Child(@TypeOf(a)) = undefined; + if (1.0 - cos_omega > delta) { + // Standard case (slerp) + const omega = std.math.acos(cos_omega); + const sin_omega = std.math.sin(omega); + scale0 = std.math.sin((1.0 - factor) * omega) / sin_omega; + scale1 = sign_scale1 * std.math.sin(factor * omega) / sin_omega; + } else { + // Quaternions are very close so we can do a linear interpolation + scale0 = 1.0 - factor; + scale1 = sign_scale1 * factor; + } + + return vector.normalize(@as(Quat(std.meta.Child(@TypeOf(a))), @splat(scale0)) * inner_a + + @as(Quat(std.meta.Child(@TypeOf(a))), @splat(scale1)) * inner_b); +} + +pub fn lerp(a: anytype, b: anytype, factor: std.meta.Child(@TypeOf(a))) Quat(std.meta.Child(@TypeOf(a))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(a)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(a)).vector.len == 4); + std.debug.assert(@typeInfo(@TypeOf(b)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(b)).vector.len == 4); + std.debug.assert(@TypeOf(a) == @TypeOf(b)); + } + + return @as(Quat(std.meta.Child(@TypeOf(a))), @splat(1.0 - factor)) * map_to_vector(a) + + @as(Quat(std.meta.Child(@TypeOf(a))), @splat(factor)) * map_to_vector(b); +} + +pub fn from_eular_angles(inAngles: anytype) @Vector(4, std.meta.Child(@TypeOf(inAngles))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(inAngles)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(inAngles)).vector.len == 3); + } + + const half = @as(@TypeOf(inAngles), @splat(0.5)) * inAngles; + const res = vector.sin_cos(half); + + const cx = res.cos_out[0]; + const sx = res.sin_out[0]; + const cy = res.cos_out[1]; + const sy = res.sin_out[1]; + const cz = res.cos_out[2]; + const sz = res.sin_out[2]; + + return .{ cz * sx * cy - sz * cx * sy, cz * cx * sy + sz * sx * cy, sz * cx * cy - cz * sx * sy, cz * cx * cy + sz * sx * sy }; +} + +pub fn from_rotation(axis: anytype, angle: std.meta.Child(@TypeOf(axis))) @Vector(4, std.meta.Child(@TypeOf(axis))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(axis)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(axis)).vector.len == 3); + } + const in_axis = map_to_vector(axis); + std.debug.assert(vector.is_normalized_default(in_axis)); + return .{ in_axis[0] * std.math.sin(angle * 0.5), in_axis[1] * std.math.sin(angle * 0.5), in_axis[2] * std.math.sin(angle * 0.5), std.math.cos(angle * 0.5) }; +} + +pub fn to_eular_angles(q: anytype) @Vector(3, std.meta.Child(@TypeOf(q))) { + comptime { + std.debug.assert(@typeInfo(@TypeOf(q)) == .vector); + std.debug.assert(@typeInfo(@TypeOf(q)).vector.len == 4); + } + + const ysqr = q[1] * q[1]; + + // roll (x-axis rotation) + const t0 = 2.0 * (q[3] * q[0] + q[1] * q[2]); + const t1 = 1.0 - 2.0 * (q[0] * q[0] + ysqr); + const roll = std.math.atan2(t0, t1); + + // pitch (y-axis rotation) + var t2 = 2.0 * (q[3] * q[1] - q[2] * q[0]); + t2 = if (t2 > 1.0) 1.0 else if (t2 < -1.0) -1.0 else t2; + const pitch = std.math.asin(t2); + + // yaw (z-axis rotation) + const t3 = 2.0 * (q[3] * q[2] + q[0] * q[1]); + const t4 = 1.0 - 2.0 * (ysqr + q[2] * q[2]); + const yaw = std.math.atan2(t3, t4); + + return .{ roll, pitch, yaw }; +} + +test from_to { + try std.testing.expect(vector.is_close_default(from_to(@Vector(3, f32){ 10, 0, 0 }, @Vector(3, f32){ 20, 0, 0 }), identity(f32))); +} + +test mul { + try std.testing.expect(vector.is_close_default(mul(Quat4f32{ 0, 1, 0, 0 }, Quat4f32{ 1, 0, 0, 0 }), Quat4f32{ 0, 0, -1, 0 })); + try std.testing.expect(vector.is_close_default(mul(Quat4f32{ 1, 0, 0, 0 }, Quat4f32{ 0, 1, 0, 0 }), Quat4f32{ 0, 0, 1, 0 })); + try std.testing.expect(vector.is_close_default(mul(Quat4f32{ 2, 3, 4, 1 }, Quat4f32{ 6, 7, 8, 5 }), Quat4f32{ 12, 30, 24, -60 })); +} + +test slerp { + const v1 = identity(f32); + const v2: Quat4f32 = from_rotation(zml.Vec3f32{ 1, 0, 0 }, 0.99 * std.math.pi); + try std.testing.expect(vector.is_close_default(slerp(v1, v2, 0.25), from_rotation(x_axis(f32), 0.25 * 0.99 * std.math.pi))); + + const v3 = vector.normalize(Quat4f32{ 1, 2, 3, 4 }); + try std.testing.expect(vector.is_close_default(slerp(v3, -v3, 0.5), v3)); +} + +test lerp { + const v1: Quat4f32 = .{ 1, 2, 3, 4 }; + const v2: Quat4f32 = .{ 5, 6, 7, 8 }; + try std.testing.expect(vector.is_close_default(lerp(v1, v2, 0.25), Quat4f32{ 2, 3, 4, 5 })); +} + +//test to_eular_angles { +// var qx: Quat4f32 = from_eular_angles(from_rotation(x_axis(f32), std.math.degreesToRadians(-10))); +// var qy: Quat4f32 = from_eular_angles(from_rotation(y_axis(f32), std.math.degreesToRadians(-20))); +// var qz: Quat4f32 = from_eular_angles(from_rotation(z_axis(f32), std.math.degreesToRadians(-30))); +//} + +test rotate_vector { + const q = from_rotation(z_axis(f32), std.math.pi / 2.0); + // 90ยฐ about +z takes +x to +y (right-handed, active rotation). + try std.testing.expect(vector.is_close(rotate_vector(q, @Vector(3, f32){ 1, 0, 0 }), .{ 0, 1, 0 }, 1e-6)); + // inverse rotation undoes it + const v = @Vector(3, f32){ 0.3, 0.5, -0.2 }; + try std.testing.expect(vector.is_close(rotate_vector_inv(q, rotate_vector(q, v)), v, 1e-6)); +} + +test to_matrix { + const q = vector.normalize(from_eular_angles(@Vector(3, f32){ 0.3, -0.6, 1.1 })); + const m = to_matrix(q); + const v = @Vector(3, f32){ 0.2, -0.7, 0.5 }; + // Applying the matrix (as a column vector) must match rotate_vector. + var mv: @Vector(3, f32) = .{ 0, 0, 0 }; + inline for (0..3) |c| { + mv += @Vector(3, f32){ m.items[c][0], m.items[c][1], m.items[c][2] } * @as(@Vector(3, f32), @splat(v[c])); + } + try std.testing.expect(vector.is_close(mv, rotate_vector(q, v), 1e-5)); +} + +test from_matrix { + const q = vector.normalize(from_eular_angles(@Vector(3, f32){ 0.3, -0.6, 1.1 })); + var back = from_matrix(to_matrix(q)); + // q and -q are the same rotation; align signs before comparing. + if (vector.dot(q, back) < 0) back = -back; + try std.testing.expect(vector.is_close(back, q, 1e-5)); +} + +test to_axis_angle { + const axis = vector.normalize(@Vector(3, f32){ 1, 2, 3 }); + const angle: f32 = 1.2; + const aa = to_axis_angle(from_rotation(axis, angle)); + try std.testing.expectApproxEqAbs(angle, aa.angle, 1e-5); + try std.testing.expect(vector.is_close(aa.axis, axis, 1e-5)); +} diff --git a/src/random.zig b/src/random.zig new file mode 100644 index 0000000..24633f9 --- /dev/null +++ b/src/random.zig @@ -0,0 +1,154 @@ +const std = @import("std"); + +const INV_2_32: f32 = 2.3283064365386963e-10; // 1 / 2^32 + +// --------------------------------------------------------------------------- +// Integer hashing / stateless PRNG steps +// --------------------------------------------------------------------------- + +/// Wellons' "lowbias32" integer hash โ€” a strong bit mixer for u32 keys. +pub fn hash_u32(x0: u32) u32 { + var x = x0; + x ^= x >> 16; + x *%= 0x7feb352d; + x ^= x >> 15; + x *%= 0x846ca68b; + x ^= x >> 16; + return x; +} + +/// Fold a new value into a running hash (boost-style combiner). +pub fn hash_combine(seed: u32, v: u32) u32 { + return seed ^ (hash_u32(v) +% 0x9e3779b9 +% (seed << 6) +% (seed >> 2)); +} + +/// Numerical Recipes linear congruential step. +pub fn lcg_next(state: u32) u32 { + return state *% 1664525 +% 1013904223; +} + +/// Marsaglia xorshift32 step (state must be non-zero). +pub fn xorshift32(state: u32) u32 { + var x = state; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return x; +} + +/// Map a u32 to a float in [0, 1) using the top 24 bits. +pub fn uint_to_float01(u: u32) f32 { + return @as(f32, @floatFromInt(u >> 8)) * (1.0 / 16777216.0); +} + +// --------------------------------------------------------------------------- +// Low-discrepancy sequences +// --------------------------------------------------------------------------- + +/// The van der Corput / Halton radical inverse of `index` in the given base. +pub fn halton(index: u32, base: u32) f32 { + var f: f32 = 1; + var r: f32 = 0; + var i = index; + const bf: f32 = @floatFromInt(base); + while (i > 0) { + f /= bf; + r += f * @as(f32, @floatFromInt(i % base)); + i /= base; + } + return r; +} + +fn radical_inverse_base2(bits0: u32) f32 { + var bits = bits0; + bits = (bits << 16) | (bits >> 16); + bits = ((bits & 0x55555555) << 1) | ((bits & 0xAAAAAAAA) >> 1); + bits = ((bits & 0x33333333) << 2) | ((bits & 0xCCCCCCCC) >> 2); + bits = ((bits & 0x0F0F0F0F) << 4) | ((bits & 0xF0F0F0F0) >> 4); + bits = ((bits & 0x00FF00FF) << 8) | ((bits & 0xFF00FF00) >> 8); + return @as(f32, @floatFromInt(bits)) * INV_2_32; +} + +/// The i-th of N points of the Hammersley set on the unit square. +pub fn hammersley_2d(i: u32, n: u32) @Vector(2, f32) { + return .{ @as(f32, @floatFromInt(i)) / @as(f32, @floatFromInt(n)), radical_inverse_base2(i) }; +} + +/// Weyl additive recurrence in 1D (golden-ratio increment). +pub fn weyl_1d(n: u32) f32 { + return @as(f32, @floatFromInt(n *% 2654435769)) * INV_2_32; +} + +/// The R2 low-discrepancy sequence (Roberts 2018), a 2D generalization of the +/// golden-ratio sequence. +pub fn weyl_2d(n: u32) @Vector(2, f32) { + return .{ + @as(f32, @floatFromInt(n *% 3242174889)) * INV_2_32, + @as(f32, @floatFromInt(n *% 2447445414)) * INV_2_32, + }; +} + +// --------------------------------------------------------------------------- +// Morton / Z-order curve (16 bits per axis) +// --------------------------------------------------------------------------- + +fn part1by1(x0: u32) u32 { + var x = x0 & 0x0000FFFF; + x = (x | (x << 8)) & 0x00FF00FF; + x = (x | (x << 4)) & 0x0F0F0F0F; + x = (x | (x << 2)) & 0x33333333; + x = (x | (x << 1)) & 0x55555555; + return x; +} + +fn compact1by1(x0: u32) u32 { + var x = x0 & 0x55555555; + x = (x | (x >> 1)) & 0x33333333; + x = (x | (x >> 2)) & 0x0F0F0F0F; + x = (x | (x >> 4)) & 0x00FF00FF; + x = (x | (x >> 8)) & 0x0000FFFF; + return x; +} + +pub fn morton2d_encode(x: u16, y: u16) u32 { + return part1by1(x) | (part1by1(y) << 1); +} + +pub fn morton2d_decode(code: u32) struct { x: u16, y: u16 } { + return .{ .x = @intCast(compact1by1(code)), .y = @intCast(compact1by1(code >> 1)) }; +} + +test "hash determinism" { + try std.testing.expectEqual(hash_u32(42), hash_u32(42)); + try std.testing.expect(hash_u32(1) != hash_u32(2)); + try std.testing.expect(hash_combine(0, 5) != hash_combine(0, 6)); +} + +test "uint_to_float01 range" { + try std.testing.expectEqual(@as(f32, 0), uint_to_float01(0)); + const f = uint_to_float01(0xFFFFFFFF); + try std.testing.expect(f >= 0 and f < 1); +} + +test "halton base 2" { + try std.testing.expectApproxEqAbs(@as(f32, 0.5), halton(1, 2), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.25), halton(2, 2), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.75), halton(3, 2), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 1.0 / 3.0), halton(1, 3), 1e-6); +} + +test "hammersley/weyl range" { + try std.testing.expectEqual(@Vector(2, f32){ 0, 0 }, hammersley_2d(0, 16)); + const w = weyl_2d(7); + try std.testing.expect(w[0] >= 0 and w[0] < 1 and w[1] >= 0 and w[1] < 1); +} + +test "morton round trip" { + const code = morton2d_encode(12345, 54321); + const d = morton2d_decode(code); + try std.testing.expectEqual(@as(u16, 12345), d.x); + try std.testing.expectEqual(@as(u16, 54321), d.y); + // Interleaving: (1,0) -> bit0, (0,1) -> bit1. + try std.testing.expectEqual(@as(u32, 1), morton2d_encode(1, 0)); + try std.testing.expectEqual(@as(u32, 2), morton2d_encode(0, 1)); +} diff --git a/src/root.zig b/src/root.zig index 2721e25..baab608 100644 --- a/src/root.zig +++ b/src/root.zig @@ -1,16 +1,96 @@ +const matrix = @import("matrix.zig"); +const std = @import("std"); + pub const vec = @import("vector.zig"); -pub const Mat = @import("matrix.zig").Mat; +pub const quat = @import("quat.zig"); +pub const Mat = matrix.Mat; + +pub const Vec4f32 = vec.Vec4f32; +pub const Vec3f32 = vec.Vec3f32; +pub const Vec2f32 = vec.Vec2f32; + +pub const Vec4f64 = vec.Vec4f64; +pub const Vec3f64 = vec.Vec3f64; +pub const Vec2f64 = vec.Vec2f64; + +pub const Quat4f32 = quat.Quat4f32; +pub const Quat4f64 = quat.Quat4f64; + +pub const Mat4f32 = matrix.Mat4f32; +pub const Mat4f64 = matrix.Mat4f64; + +pub const Mat3f32 = Mat(f32, 3, 3); +pub const Mat3f64 = Mat(f64, 3, 3); + +pub const geom = @import("geometry.zig"); +pub const meta = @import("meta.zig"); +pub const simd = @import("simd.zig"); +pub const scalar = @import("scalar.zig"); +pub const packing = @import("packing.zig"); +pub const color = @import("color.zig"); +pub const random = @import("random.zig"); test { - @import("std").testing.refAllDeclsRecursive(@This()); + _ = @import("vector.zig"); + _ = @import("matrix.zig"); + _ = @import("quat.zig"); + _ = @import("geometry.zig"); + _ = @import("meta.zig"); + _ = @import("simd.zig"); + _ = @import("scalar.zig"); + _ = @import("packing.zig"); + _ = @import("color.zig"); + _ = @import("random.zig"); } -const std = @import("std"); - -pub fn toRadians(degrees: anytype) @TypeOf(degrees) { +pub fn to_radians(degrees: anytype) @TypeOf(degrees) { return degrees * (std.math.pi / 180); } -pub fn toDegrees(radians: anytype) @TypeOf(radians) { +pub fn to_degrees(radians: anytype) @TypeOf(radians) { return radians * (1 / (std.math.pi / 180)); } + +pub fn find_roots(comptime T: type, a: T, b: T, c: T) struct { + num_roots: u8, + roots: [2]T, +} { + // Check if this is a linear equation + if (a == 0) { + // Check if this is a constant equation + if (b == 0) + return .{ + .num_roots = 0, + .roots = .{ 0, 0 }, + }; + + // Linear equation with 1 solution + const r1 = -c / b; + return .{ + .num_roots = 1, + .roots = .{ r1, r1 }, + }; + } + + // See Numerical Recipes in C, Chapter 5.6 Quadratic and Cubic Equations + const det: T = (b * b) - 4 * a * c; + if (det < 0) + return .{ + .num_roots = 0, + .roots = .{ 0, 0 }, + }; + + const q: T = (b + std.math.sign(b) * std.math.sqrt(det)) / -2; + const r1 = q / a; + if (q == 0) { + return .{ + .num_roots = 1, + .roots = .{ r1, r1 }, + }; + } + const r2 = c / q; + return .{ + .num_roots = 2, + .roots = .{ r1, r2 }, + }; +} diff --git a/src/scalar.zig b/src/scalar.zig new file mode 100644 index 0000000..ff79b09 --- /dev/null +++ b/src/scalar.zig @@ -0,0 +1,81 @@ +const std = @import("std"); + +/// Scalar element type of a scalar-or-vector type (`f32` -> `f32`, `@Vector(n, f32)` -> `f32`). +fn Elem(comptime V: type) type { + return switch (@typeInfo(V)) { + .vector => |info| info.child, + else => V, + }; +} + +/// Broadcast a scalar to `V` (splat for vectors, identity for scalars). +fn splatLike(comptime V: type, value: Elem(V)) V { + return switch (@typeInfo(V)) { + .vector => @splat(value), + else => value, + }; +} + +/// Clamp `x` to the inclusive range [lo, hi]. Works on scalars and float vectors. +pub fn clamp(x: anytype, lo: Elem(@TypeOf(x)), hi: Elem(@TypeOf(x))) @TypeOf(x) { + const V = @TypeOf(x); + return @min(@max(x, splatLike(V, lo)), splatLike(V, hi)); +} + +/// Clamp `x` to [0, 1]. +pub fn saturate(x: anytype) @TypeOf(x) { + return clamp(x, 0, 1); +} + +/// Linear interpolation: a + (b - a) * t. +pub fn lerp(a: anytype, b: @TypeOf(a), t: Elem(@TypeOf(a))) @TypeOf(a) { + const V = @TypeOf(a); + return a + (b - a) * splatLike(V, t); +} + +/// Multiply-add: a * b + c. +pub fn madd(a: anytype, b: @TypeOf(a), c: @TypeOf(a)) @TypeOf(a) { + return a * b + c; +} + +/// 0 if x < edge, else 1 (component-wise). `edge` and `x` share the same type. +pub fn step(edge: anytype, x: @TypeOf(edge)) @TypeOf(edge) { + const V = @TypeOf(edge); + return switch (@typeInfo(V)) { + .vector => @select(Elem(V), x >= edge, splatLike(V, 1), splatLike(V, 0)), + else => @as(V, if (x >= edge) 1 else 0), + }; +} + +/// Saturated linear ramp from edge0 to edge1. +pub fn linearstep(edge0: anytype, edge1: @TypeOf(edge0), x: @TypeOf(edge0)) @TypeOf(edge0) { + return saturate((x - edge0) / (edge1 - edge0)); +} + +/// Hermite smoothstep: 0 below edge0, 1 above edge1, smooth in between. +pub fn smoothstep(edge0: anytype, edge1: @TypeOf(edge0), x: @TypeOf(edge0)) @TypeOf(edge0) { + const V = @TypeOf(edge0); + const t = saturate((x - edge0) / (edge1 - edge0)); + return t * t * (splatLike(V, 3) - splatLike(V, 2) * t); +} + +test "scalar helpers" { + try std.testing.expectApproxEqAbs(@as(f32, 0.5), lerp(@as(f32, 0), @as(f32, 1), 0.5), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 7), madd(@as(f32, 2), @as(f32, 3), @as(f32, 1)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0), saturate(@as(f32, -3)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 1), saturate(@as(f32, 3)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.75), clamp(@as(f32, 5), 0.25, 0.75), 1e-6); + + try std.testing.expectApproxEqAbs(@as(f32, 0), smoothstep(@as(f32, 0), @as(f32, 1), @as(f32, 0)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 1), smoothstep(@as(f32, 0), @as(f32, 1), @as(f32, 1)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.5), smoothstep(@as(f32, 0), @as(f32, 1), @as(f32, 0.5)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0.25), linearstep(@as(f32, 0), @as(f32, 4), @as(f32, 1)), 1e-6); + + try std.testing.expectApproxEqAbs(@as(f32, 1), step(@as(f32, 0.5), @as(f32, 0.7)), 1e-6); + try std.testing.expectApproxEqAbs(@as(f32, 0), step(@as(f32, 0.5), @as(f32, 0.2)), 1e-6); + + // Vector variants. + try std.testing.expectEqual(@Vector(2, f32){ 1, 2 }, lerp(@Vector(2, f32){ 0, 0 }, @Vector(2, f32){ 2, 4 }, 0.5)); + try std.testing.expectEqual(@Vector(2, f32){ 0, 1 }, saturate(@Vector(2, f32){ -1, 2 })); + try std.testing.expectEqual(@Vector(2, f32){ 0, 1 }, step(@Vector(2, f32){ 0.5, 0.5 }, @Vector(2, f32){ 0.2, 0.7 })); +} diff --git a/src/simd.zig b/src/simd.zig new file mode 100644 index 0000000..9762d12 --- /dev/null +++ b/src/simd.zig @@ -0,0 +1,260 @@ +//! Chunked SIMD kernels for ARRAY inputs. +//! +//! `vector.zig` accepts both `@Vector(N,T)` and `[N]T`. A `@Vector` input is +//! processed as a single native op (the caller chose that width). An array, +//! however, may be arbitrarily long (`[245]f32`); coercing it to one wide +//! `@Vector(245,f32)` makes LLVM emit one enormous SIMD op and bloats the +//! binary. These kernels instead walk an array in chunks sized by the CPU's +//! native SIMD width (`laneCount`), using a runtime loop over full chunks plus +//! one comptime-sized remainder chunk, so no wide vector is ever materialized. + +const std = @import("std"); +const builtin = @import("builtin"); +const meta = @import("meta.zig"); + +/// Native SIMD width (lane count) for element type `T` on the target CPU, or 1 +/// when the target has no vector unit for `T`. This is the chunk size used by +/// every kernel below. +pub fn laneCount(comptime T: type) comptime_int { + return std.simd.suggestVectorLengthForCpu(T, builtin.cpu) orelse 1; +} + +/// Sum of `a[i] * b[i]`, accumulated in compute type `C`. +/// +/// `C` must be either `Child(a)` itself (identity coercion) or a WIDER float +/// (element-wise float widen) โ€” never an int->float change, which is illegal on +/// a vector. The int path keeps `C == T` and the caller applies `@floatFromInt` +/// once to the scalar result, matching `vector.norm_sqr_adv`. +pub fn sumProducts(comptime C: type, a: anytype, b: anytype) C { + const T = meta.Child(@TypeOf(a)); + const N = comptime meta.lengthOf(@TypeOf(a)); + const W = comptime laneCount(T); + const full = comptime (N / W) * W; // largest multiple of W <= N + const rem = comptime N % W; // leftover lanes + + // Bind to const locals so `x[i..]` (runtime i) is addressable. + const ba = a; + const bb = b; + + var acc: @Vector(W, C) = @splat(0); + var i: usize = 0; + while (i < full) : (i += W) { // runtime loop -> bounded code size + const ca: @Vector(W, T) = ba[i..][0..W].*; + const cb: @Vector(W, T) = bb[i..][0..W].*; + const va: @Vector(W, C) = ca; // identity or float-widen + const vb: @Vector(W, C) = cb; + acc += va * vb; + } + var total: C = @reduce(.Add, acc); + if (rem != 0) { // comptime guard: @Vector(0, C) is illegal + const ca: @Vector(rem, T) = ba[full..][0..rem].*; + const cb: @Vector(rem, T) = bb[full..][0..rem].*; + const va: @Vector(rem, C) = ca; + const vb: @Vector(rem, C) = cb; + total += @reduce(.Add, va * vb); + } + return total; +} + +/// Sum of `a[i] * a[i]`, accumulated in compute type `C`. +pub fn sumSquares(comptime C: type, a: anytype) C { + return sumProducts(C, a, a); +} + +/// Apply width-generic `op(chunk, ctx)` to every native-width chunk of `a`, +/// returning a fresh array of the same element type and length. `op` must be +/// generic over the chunk width: `fn (chunk: anytype, ctx) @TypeOf(chunk)`. +pub fn mapUnary(a: anytype, ctx: anytype, comptime op: anytype) [meta.lengthOf(@TypeOf(a))]meta.Child(@TypeOf(a)) { + const T = meta.Child(@TypeOf(a)); + const N = comptime meta.lengthOf(@TypeOf(a)); + const W = comptime laneCount(T); + const full = comptime (N / W) * W; + const rem = comptime N % W; + + const ba = a; + var out: [N]T = undefined; + var i: usize = 0; + while (i < full) : (i += W) { + const c: @Vector(W, T) = ba[i..][0..W].*; + out[i..][0..W].* = op(c, ctx); // @Vector -> [W]T store coercion + } + if (rem != 0) { + const c: @Vector(rem, T) = ba[full..][0..rem].*; + out[full..][0..rem].* = op(c, ctx); + } + return out; +} + +/// Binary version of `mapUnary`: `op(chunkA, chunkB, ctx) @TypeOf(chunkA)`. +pub fn mapBinary(a: anytype, b: anytype, ctx: anytype, comptime op: anytype) [meta.lengthOf(@TypeOf(a))]meta.Child(@TypeOf(a)) { + const T = meta.Child(@TypeOf(a)); + const N = comptime meta.lengthOf(@TypeOf(a)); + const W = comptime laneCount(T); + const full = comptime (N / W) * W; + const rem = comptime N % W; + + const ba = a; + const bb = b; + var out: [N]T = undefined; + var i: usize = 0; + while (i < full) : (i += W) { + const ca: @Vector(W, T) = ba[i..][0..W].*; + const cb: @Vector(W, T) = bb[i..][0..W].*; + out[i..][0..W].* = op(ca, cb, ctx); + } + if (rem != 0) { + const ca: @Vector(rem, T) = ba[full..][0..rem].*; + const cb: @Vector(rem, T) = bb[full..][0..rem].*; + out[full..][0..rem].* = op(ca, cb, ctx); + } + return out; +} + +/// Element-wise `a - b` returned as an ARRAY, so a following reduction chunks it +/// instead of materializing one wide `@Vector`. Used by distance / is_close. +pub fn sub(a: anytype, b: anytype) [meta.lengthOf(@TypeOf(a))]meta.Child(@TypeOf(a)) { + return mapBinary(a, b, {}, struct { + fn op(x: anytype, y: anytype, _: void) @TypeOf(x) { + return x - y; + } + }.op); +} + +/// Width-generic sin/cos core (cephes-derived). Operates on ONE `@Vector` of any +/// width and float element type. Called by both the `@Vector` public path and +/// the per-chunk array path so the math lives in exactly one place. +/// +/// Implementation based on sinf.c from the cephes library, combining sinf and +/// cosf, changing octants to quadrants and vectorizing it. +/// Original implementation by Stephen L. Moshier (See: http://www.moshier.net/) +pub fn sinCosVec(in: anytype) struct { sin: @TypeOf(in), cos: @TypeOf(in) } { + const FVec = @TypeOf(in); + const info = switch (@typeInfo(FVec)) { + .vector => |v| v, + else => @compileError("sinCosVec expects a @Vector, got: " ++ @typeName(FVec)), + }; + const w = info.len; + const F = info.child; + const bits = switch (@typeInfo(F)) { + .float => |f| f.bits, + else => @compileError("sinCosVec only supports floating point vectors"), + }; + const UVec = @Vector(w, switch (bits) { + 32 => u32, + 64 => u64, + else => @compileError("Unsupported float size"), + }); + const last_bit = switch (bits) { + 32 => 0x80000000, + 64 => 0x8000000000000000, + else => unreachable, + }; + const num_bits = bits; + + // Make argument positive and remember sign for sin only since cos is symmetric around x (highest bit of a float is the sign bit) + var sin_sign = @as(UVec, @bitCast(in)) & @as(UVec, @splat(last_bit)); + var x: FVec = @bitCast(@as(UVec, @bitCast(in)) ^ sin_sign); + + // x / (PI / 2) rounded to nearest int gives us the quadrant closest to x + const quadrant: UVec = @intFromFloat(@as(FVec, @splat(0.6366197723675814)) * x + @as(FVec, @splat(0.5))); + const float_quadrant: FVec = @floatFromInt(quadrant); + + // Make x relative to the closest quadrant using a two step Cody-Waite argument reduction. + // After this we have x in the range [-PI / 4, PI / 4]. + x = ((x - float_quadrant * @as(FVec, @splat(1.5703125))) - float_quadrant * @as(FVec, @splat(0.0004837512969970703125))) - float_quadrant * @as(FVec, @splat(7.549789948768648e-8)); + + // Calculate x2 = x^2 + const x2 = x * x; + + // Taylor expansion: + // Cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! + ... = (((x2/8!- 1/6!) * x2 + 1/4!) * x2 - 1/2!) * x2 + 1 + const taylor_cos = ((@as(FVec, @splat(2.443315711809948e-5)) * x2 - @as(FVec, @splat(1.388731625493765e-3))) * x2 + @as(FVec, @splat(4.166664568298827e-2))) * x2 * x2 - @as(FVec, @splat(0.5)) * x2 + @as(FVec, @splat(1.0)); + // Sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ... = ((-x2/7! + 1/5!) * x2 - 1/3!) * x2 * x + x + const taylor_sin = ((@as(FVec, @splat(-1.9515295891e-4)) * x2 + @as(FVec, @splat(8.3321608736e-3))) * x2 - @as(FVec, @splat(1.6666654611e-1))) * x2 * x + x; + + // The lowest 2 bits of quadrant indicate the quadrant that we are in. Extract + // them to determine which Taylor expansion to use and the signs. + const bit1: UVec = quadrant << @as(UVec, @splat(num_bits - 1)); // bit 0 + const bit2: UVec = (quadrant << @as(UVec, @splat(num_bits - 2))) & @as(UVec, @splat(last_bit)); // bit 1 + + // Select which one of the results is sin and which one is cos based on bit1 + const s = @select(F, bit1 > @as(UVec, @splat(0)), taylor_cos, taylor_sin); + const c = @select(F, bit1 > @as(UVec, @splat(0)), taylor_sin, taylor_cos); + + sin_sign = sin_sign ^ bit2; + const cos_sign = bit1 ^ bit2; + + return .{ + .sin = @as(FVec, @bitCast(@as(UVec, @bitCast(s)) ^ sin_sign)), + .cos = @as(FVec, @bitCast(@as(UVec, @bitCast(c)) ^ cos_sign)), + }; +} + +/// Chunked array driver for sin/cos: returns two `[N]T` arrays. +pub fn sinCos(a: anytype) struct { + sin: [meta.lengthOf(@TypeOf(a))]meta.Child(@TypeOf(a)), + cos: [meta.lengthOf(@TypeOf(a))]meta.Child(@TypeOf(a)), +} { + const T = meta.Child(@TypeOf(a)); + const N = comptime meta.lengthOf(@TypeOf(a)); + const W = comptime laneCount(T); + const full = comptime (N / W) * W; + const rem = comptime N % W; + + const ba = a; + var s: [N]T = undefined; + var c: [N]T = undefined; + var i: usize = 0; + while (i < full) : (i += W) { + const chunk: @Vector(W, T) = ba[i..][0..W].*; + const r = sinCosVec(chunk); + s[i..][0..W].* = r.sin; + c[i..][0..W].* = r.cos; + } + if (rem != 0) { + const chunk: @Vector(rem, T) = ba[full..][0..rem].*; + const r = sinCosVec(chunk); + s[full..][0..rem].* = r.sin; + c[full..][0..rem].* = r.cos; + } + return .{ .sin = s, .cos = c }; +} + +test laneCount { + // Always at least 1 lane so chunking is well-defined on any target. + try std.testing.expect(laneCount(f32) >= 1); + try std.testing.expect(laneCount(i8) >= 1); +} + +test sumSquares { + // Chunked reduction over a remainder-forcing length. + const a = [_]f32{ 1, 2, 3, 4, 5, 6, 7 }; + try std.testing.expectApproxEqAbs(@as(f32, 140), sumSquares(f32, a), 1e-3); + // Integer accumulation is exact and associative. + const ai = [_]i32{ 1, 2, 3, 4, 5 }; + try std.testing.expectEqual(@as(i32, 55), sumSquares(i32, ai)); +} + +test sumProducts { + const a = [_]f32{ 1, 2, 3, 4, 5 }; + const b = [_]f32{ 5, 4, 3, 2, 1 }; + try std.testing.expectApproxEqAbs(@as(f32, 35), sumProducts(f32, a, b), 1e-3); +} + +test sub { + const a = [_]f32{ 4, 6, 8 }; + const b = [_]f32{ 1, 2, 3 }; + const d = sub(a, b); + try std.testing.expect(@TypeOf(d) == [3]f32); + try std.testing.expectEqual([3]f32{ 3, 4, 5 }, d); +} + +test mapUnary { + const a = [_]f32{ 1, 2, 3, 4, 5 }; + const out = mapUnary(a, @as(f32, 2), struct { + fn op(c: anytype, s: f32) @TypeOf(c) { + return c * @as(@TypeOf(c), @splat(s)); + } + }.op); + try std.testing.expectEqual([5]f32{ 2, 4, 6, 8, 10 }, out); +} diff --git a/src/testing.zig b/src/testing.zig new file mode 100644 index 0000000..0f05536 --- /dev/null +++ b/src/testing.zig @@ -0,0 +1,10 @@ +const vector = @import("vector.zig"); +const std = @import("std"); + +pub fn expect_is_close(vec: anytype, other: @TypeOf(vec), epsilon: std.meta.Child(@TypeOf(vec))) !void{ + if(vector.distance(vec, other) > @as(std.meta.Child(@TypeOf(vec)), epsilon)) { + std.debug.print("Vectors not close enough: {} vs {}\n", .{vec, other}); + return error.TestExpectedApproxIsClose; + } +} + diff --git a/src/vector.zig b/src/vector.zig index 098d352..9ea3f52 100644 --- a/src/vector.zig +++ b/src/vector.zig @@ -1,10 +1,40 @@ +// Jolt Physics Library (https://github.com/jrouwe/JoltPhysics) +// SPDX-FileCopyrightText: 2021 Jorrit Rouwe +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2025 Michael Pollind const std = @import("std"); - const Float = std.meta.Float; -pub fn info(T: type) std.builtin.Type.Vector { - if (@typeInfo(T) != .vector) @compileError("Expected a @Vector type got: " ++ @typeName(T)); - return @typeInfo(T).vector; +const Mat = @import("matrix.zig").Mat; +const testing = @import("testing.zig"); +const meta = @import("meta.zig"); +const simd = @import("simd.zig"); + +pub const Vec4f32 = @Vector(4, f32); +pub const Vec3f32 = @Vector(3, f32); +pub const Vec2f32 = @Vector(2, f32); + +pub const Vec4f64 = @Vector(4, f64); +pub const Vec3f64 = @Vector(3, f64); +pub const Vec2f64 = @Vector(2, f64); + +pub fn to_mat(vec: anytype) Mat(meta.Child(@TypeOf(vec)), 1, meta.lengthOf(@TypeOf(vec))) { + var result: Mat(meta.Child(@TypeOf(vec)), 1, meta.lengthOf(@TypeOf(vec))) = .zero; + inline for (0..meta.lengthOf(@TypeOf(vec))) |i| { + result.items[0][i] = vec[i]; + } + return result; +} + +pub fn extract(vec: anytype, comptime len: usize) meta.Reshape(@TypeOf(vec), len) { + comptime { + std.debug.assert(len <= meta.lengthOf(@TypeOf(vec))); + } + var result: @Vector(len, meta.Child(@TypeOf(vec))) = undefined; + inline for (0..len) |i| { + result[i] = vec[i]; + } + return result; } /// Returns a new vector with the components swizzled. @@ -22,8 +52,9 @@ pub fn info(T: type) std.builtin.Type.Vector { /// const v2 = v.sw("wxx"); /// ``` /// here v2 is equal to `Vec(3, f32).init(.{ 4, 1, 1 });` -pub fn sw(vec: anytype, comptime components: []const u8) @Vector(components.len, info(@TypeOf(vec)).child) { - const T = info(@TypeOf(vec)).child; +pub fn swizzle(vec: anytype, comptime components: []const u8) meta.Reshape(@TypeOf(vec), components.len) { + const T = meta.Child(@TypeOf(vec)); + const v: meta.AsVector(@TypeOf(vec)) = vec; comptime var mask: [components.len]u8 = undefined; comptime var i: usize = 0; inline for (components) |c| { @@ -39,37 +70,81 @@ pub fn sw(vec: anytype, comptime components: []const u8) @Vector(components.len, return @shuffle( T, - vec, + v, @as(@Vector(1, T), undefined), mask, ); } +pub fn norm_sqr_adv(vec: anytype, comptime precision: u8) Float(precision) { + const T = meta.Child(@TypeOf(vec)); + // Array inputs are processed in native-width SIMD chunks (see simd.zig) so a + // large array never materializes one huge @Vector. + if (@typeInfo(@TypeOf(vec)) == .array) { + if (@typeInfo(T) == .int) return @floatFromInt(simd.sumSquares(T, vec)); + const C = if (precision > @bitSizeOf(T)) Float(precision) else T; + return @floatCast(simd.sumSquares(C, vec)); + } + const v: meta.AsVector(@TypeOf(vec)) = vec; + if (@typeInfo(T) == .int) { + return @as( + Float(precision), + @floatFromInt(@reduce( + .Add, + v * v, + )), + ); + } else { + const items: blk: { + if (precision > @bitSizeOf(T)) { + break :blk @Vector(meta.lengthOf(@TypeOf(vec)), Float(precision)); + } else { + break :blk meta.AsVector(@TypeOf(vec)); + } + } = v; + + return @floatCast(@reduce( + .Add, + items * items, + )); + } +} + +pub inline fn norm_sqr(vec: anytype) Float(@bitSizeOf(meta.Child(@TypeOf(vec)))) { + return norm_sqr_adv(vec, @bitSizeOf(meta.Child(@TypeOf(vec)))); +} + /// Returns the norm of the vector as a Float. /// /// the precsion parameter is the number of bits of the output. /// the precision of the calculations will match the precision of the output type. -pub fn normAdv(vec: anytype, comptime precision: u8) Float(precision) { - const T = info(@TypeOf(vec)).child; - +pub fn norm_adv(vec: anytype, comptime precision: u8) Float(precision) { + const T = meta.Child(@TypeOf(vec)); + if (@typeInfo(@TypeOf(vec)) == .array) { + if (@typeInfo(T) == .int) + return std.math.sqrt(@as(Float(precision), @floatFromInt(simd.sumSquares(T, vec)))); + const C = if (precision > @bitSizeOf(T)) Float(precision) else T; + return @floatCast(std.math.sqrt(simd.sumSquares(C, vec))); + } + const v: meta.AsVector(@TypeOf(vec)) = vec; if (@typeInfo(T) == .int) { return std.math.sqrt( @as( Float(precision), @floatFromInt(@reduce( .Add, - vec * vec, + v * v, )), ), ); } else { const items: blk: { if (precision > @bitSizeOf(T)) { - break :blk @Vector(info(@TypeOf(vec)).len, Float(precision)); + break :blk @Vector(meta.lengthOf(@TypeOf(vec)), Float(precision)); } else { - break :blk @TypeOf(vec); + break :blk meta.AsVector(@TypeOf(vec)); } - } = vec; + } = v; return @floatCast(std.math.sqrt( @reduce( @@ -83,36 +158,73 @@ pub fn normAdv(vec: anytype, comptime precision: u8) Float(precision) { /// Returns the norm of the vector as a Float with the default precision. /// the precision of the output is the number of bits of the output. /// see `normAdv` for more information. -pub inline fn norm(vec: anytype) Float(@bitSizeOf(info(@TypeOf(vec)).child)) { - return normAdv(vec, @bitSizeOf(info(@TypeOf(vec)).child)); +pub inline fn norm(vec: anytype) Float(@bitSizeOf(meta.Child(@TypeOf(vec)))) { + return norm_adv(vec, @bitSizeOf(meta.Child(@TypeOf(vec)))); } /// Returns a new vector with the same direction as the original vector, but with a norm closest to 1. pub fn normalize(vec: anytype) @TypeOf(vec) { - const self_norm = switch (@typeInfo(info(@TypeOf(vec)).child)) { - .float => norm(vec), - .int => @as(info(@TypeOf(vec)).child, @intFromFloat(norm(vec))), + const T = meta.Child(@TypeOf(vec)); + if (@typeInfo(@TypeOf(vec)) == .array) { + const self_norm = switch (@typeInfo(T)) { + .float => norm(vec), + .int => @as(T, @intFromFloat(norm(vec))), + else => unreachable, + }; + return simd.mapUnary(vec, self_norm, struct { + fn op(c: anytype, d: T) @TypeOf(c) { + return c / @as(@TypeOf(c), @splat(d)); + } + }.op); + } + const v: meta.AsVector(@TypeOf(vec)) = vec; + const self_norm = switch (@typeInfo(T)) { + .float => norm(v), + .int => @as(T, @intFromFloat(norm(v))), else => unreachable, }; - return vec / @as( - @TypeOf(vec), + return v / @as( + meta.AsVector(@TypeOf(vec)), @splat(self_norm), ); } /// dot product of two vectors. -pub fn dot(vec: anytype, other: @TypeOf(vec)) info(@TypeOf(vec)).child { - return @reduce(.Add, vec * other); +pub fn dot(vec: anytype, other: @TypeOf(vec)) meta.Child(@TypeOf(vec)) { + if (@typeInfo(@TypeOf(vec)) == .array) + return simd.sumProducts(meta.Child(@TypeOf(vec)), vec, other); + const a: meta.AsVector(@TypeOf(vec)) = vec; + const b: meta.AsVector(@TypeOf(vec)) = other; + return @reduce(.Add, a * b); +} + +pub fn norm_perpendicular(a: anytype) @TypeOf(a) { + comptime { + std.debug.assert(meta.lengthOf(@TypeOf(a)) == 3); + } + const V = meta.AsVector(@TypeOf(a)); + const v: V = a; + if (@abs(v[0]) > @abs(v[1])) { + const len = norm(swizzle(v, "xz")); + return V{ v[2], 0.0, -v[0] } / @as(V, @splat(len)); + } else { + const len = norm(swizzle(v, "yz")); + return V{ 0.0, v[2], -v[1] } / @as(V, @splat(len)); + } } /// Returns the cross product of two vectors. -pub fn cross(self: anytype, other: @TypeOf(self)) @TypeOf(self) { - if (info(@TypeOf(self)).len != 3) @compileError("vector must have three elements for cross() to be defined"); - const T = info(@TypeOf(self)).child; - const self1 = @shuffle(T, self, self, [3]u8{ 1, 2, 0 }); - const self2 = @shuffle(T, self, self, [3]u8{ 2, 0, 1 }); - const other1 = @shuffle(T, other, other, [3]u8{ 1, 2, 0 }); - const other2 = @shuffle(T, other, other, [3]u8{ 2, 0, 1 }); +pub fn cross(a: anytype, b: @TypeOf(a)) @TypeOf(a) { + comptime { + std.debug.assert(meta.lengthOf(@TypeOf(a)) == 3); + } + const T = meta.Child(@TypeOf(a)); + const av: meta.AsVector(@TypeOf(a)) = a; + const bv: meta.AsVector(@TypeOf(a)) = b; + const self1 = @shuffle(T, av, av, [3]u8{ 1, 2, 0 }); + const self2 = @shuffle(T, av, av, [3]u8{ 2, 0, 1 }); + const other1 = @shuffle(T, bv, bv, [3]u8{ 1, 2, 0 }); + const other2 = @shuffle(T, bv, bv, [3]u8{ 2, 0, 1 }); return self1 * other2 - self2 * other1; } @@ -121,42 +233,168 @@ pub fn cross(self: anytype, other: @TypeOf(self)) @TypeOf(self) { /// /// the precsion parameter is the number of bits of the Vector type T. /// the precision of the calculations will match the precision of the output type. -pub fn distanceAdv(vec: anytype, other: @TypeOf(vec), comptime precision: u8) Float(precision) { - return normAdv(vec - other, precision); +pub fn distance_adv(a: anytype, b: @TypeOf(a), comptime precision: u8) Float(precision) { + if (@typeInfo(@TypeOf(a)) == .array) + return norm_adv(simd.sub(a, b), precision); // sub -> [N]T, norm_adv chunks it + const av: meta.AsVector(@TypeOf(a)) = a; + const bv: meta.AsVector(@TypeOf(a)) = b; + return norm_adv(av - bv, precision); } /// Returns the distance between two vectors. /// /// the precision of the output is the number of bits of T. /// see `distanceAdv` for more information. -pub inline fn distance(vec: anytype, other: @TypeOf(vec)) Float(@bitSizeOf(info(@TypeOf(vec)).child)) { - return distanceAdv(vec, other, @bitSizeOf(info(@TypeOf(vec)).child)); +pub inline fn distance(vec: anytype, other: @TypeOf(vec)) Float(@bitSizeOf(meta.Child(@TypeOf(vec)))) { + return distance_adv(vec, other, @bitSizeOf(meta.Child(@TypeOf(vec)))); } /// Returns the angle between two vectors. /// /// the precsion parameter is the number of bits of the Vector type T. /// the precision of the calculations will match the precision of the output type. -pub fn angleAdv(vec: anytype, other: @TypeOf(vec), comptime precision: u8) Float(precision) { - return std.math.acos(dot(vec, other) / (normAdv(vec, precision) * normAdv(other, precision))); +pub fn angle_adv(vec: anytype, other: @TypeOf(vec), comptime precision: u8) Float(precision) { + return std.math.acos(dot(vec, other) / (norm_adv(vec, precision) * norm_adv(other, precision))); } /// Returns the angle between two vectors. -pub inline fn angle(vec: anytype, other: @TypeOf(vec)) info(@TypeOf(vec)).child { - return angleAdv(vec, other, @bitSizeOf(info(@TypeOf(vec)).child)); +pub inline fn angle(a: anytype, b: @TypeOf(a)) meta.Child(@TypeOf(a)) { + return angle_adv(a, b, @bitSizeOf(meta.Child(@TypeOf(a)))); } /// Returns a new vector that is the reflection of the original vector on the given normal. pub fn reflect(vec: anytype, normal: @TypeOf(vec)) @TypeOf(vec) { - const dot_product = dot(vec, normal); - return vec - (normal * - @as(@TypeOf(normal), @splat(2)) * - @as(@TypeOf(normal), @splat(dot_product))); + const T = meta.Child(@TypeOf(vec)); + if (@typeInfo(@TypeOf(vec)) == .array) { + const dp = dot(vec, normal); // chunked dot + return simd.mapBinary(vec, normal, dp, struct { + fn op(cv: anytype, cn: anytype, d: T) @TypeOf(cv) { + const two: @TypeOf(cv) = @splat(2); + const dd: @TypeOf(cv) = @splat(d); + return cv - cn * two * dd; + } + }.op); + } + const V = meta.AsVector(@TypeOf(vec)); + const v: V = vec; + const n: V = normal; + const dot_product = dot(v, n); + return v - (n * + @as(V, @splat(2)) * + @as(V, @splat(dot_product))); +} + +/// Simultaneously computes the sine and cosine of each component. +/// The cephes-derived core lives in `simd.sinCosVec`; array inputs are processed +/// in native-width SIMD chunks via `simd.sinCos`, `@Vector` inputs in one op. +pub fn sin_cos(input: anytype) struct { sin_out: @TypeOf(input), cos_out: @TypeOf(input) } { + if (@typeInfo(@TypeOf(input)) == .array) { + const r = simd.sinCos(input); + return .{ .sin_out = r.sin, .cos_out = r.cos }; + } + const r = simd.sinCosVec(input); // input is a @Vector -> width == its len + return .{ .sin_out = r.sin, .cos_out = r.cos }; } /// Returns a new vector with a direction closest to the original vector, but with a magnitude scaled by the given value. -pub inline fn scale(vec: anytype, value: info(@TypeOf(vec)).child) @TypeOf(vec) { - return vec * @as(@TypeOf(vec), @splat(value)); +pub inline fn scale(a: anytype, value: meta.Child(@TypeOf(a))) @TypeOf(a) { + if (@typeInfo(@TypeOf(a)) == .array) + return simd.mapUnary(a, value, struct { + fn op(c: anytype, s: meta.Child(@TypeOf(a))) @TypeOf(c) { + return c * @as(@TypeOf(c), @splat(s)); + } + }.op); + const V = meta.AsVector(@TypeOf(a)); + const v: V = a; + return v * @as(V, @splat(value)); +} + +pub fn is_close_default(a: anytype, b: @TypeOf(a)) bool { + return is_close(a, b, 1.0e-12); +} + +pub fn is_close(a: anytype, b: @TypeOf(a), max_distance_sqr: Float(@bitSizeOf(meta.Child(@TypeOf(a))))) bool { + if (@typeInfo(@TypeOf(a)) == .array) + return norm_sqr(simd.sub(a, b)) <= max_distance_sqr; + const av: meta.AsVector(@TypeOf(a)) = a; + const bv: meta.AsVector(@TypeOf(a)) = b; + return norm_sqr(av - bv) <= max_distance_sqr; +} + +pub fn is_normalized_default(a: anytype) bool { + return is_normalized(a, 1.0e-6); +} + +pub fn is_normalized(a: anytype, tolerance: Float(@bitSizeOf(meta.Child(@TypeOf(a))))) bool { + return @abs(norm_sqr(a) - 1.0) <= tolerance; +} + +/// Encode a unit vector into a 2-component octahedral representation in [-1, 1]^2. +/// Useful for compact normal storage. See `decode_oct` for the inverse. +pub fn encode_oct(n: anytype) @Vector(2, meta.Child(@TypeOf(n))) { + const T = meta.Child(@TypeOf(n)); + const l1 = @abs(n[0]) + @abs(n[1]) + @abs(n[2]); + var p = @Vector(2, T){ n[0] / l1, n[1] / l1 }; + if (n[2] <= 0) { + // Capture components first: `p = .{...}` builds in place, so reading + // p[0]/p[1] on the RHS would otherwise see partially-updated values. + const px = p[0]; + const py = p[1]; + const sx: T = if (px >= 0) 1 else -1; + const sy: T = if (py >= 0) 1 else -1; + p = .{ (1 - @abs(py)) * sx, (1 - @abs(px)) * sy }; + } + return p; +} + +/// Decode a 2-component octahedral representation back into a unit vector. +pub fn decode_oct(e: anytype) @Vector(3, meta.Child(@TypeOf(e))) { + const T = meta.Child(@TypeOf(e)); + var n = @Vector(3, T){ e[0], e[1], 1 - @abs(e[0]) - @abs(e[1]) }; + const t = @max(-n[2], 0); + n[0] += if (n[0] >= 0) -t else t; + n[1] += if (n[1] >= 0) -t else t; + return normalize(n); +} + +/// Build a right-handed orthonormal basis (tangent, bitangent) around a unit +/// normal using Duff et al.'s branchless method. `n` is the third basis vector. +pub fn build_basis(n: anytype) struct { tangent: @TypeOf(n), bitangent: @TypeOf(n) } { + const T = meta.Child(@TypeOf(n)); + const s: T = std.math.copysign(@as(T, 1), n[2]); + const a = -1.0 / (s + n[2]); + const b = n[0] * n[1] * a; + return .{ + .tangent = .{ 1.0 + s * n[0] * n[0] * a, s * b, -s * n[0] }, + .bitangent = .{ b, s + n[1] * n[1] * a, -n[1] }, + }; +} + +test encode_oct { + const dirs = [_]@Vector(3, f32){ + .{ 0, 0, 1 }, + .{ 0, 0, -1 }, + .{ 1, 0, 0 }, + .{ 0, -1, 0 }, + .{ 1, 2, 3 }, + .{ -3, 2, -1 }, + .{ 0.5, -0.5, -0.7 }, + .{ -2, -2, -2 }, + }; + for (dirs) |d| { + const dn = normalize(d); + try std.testing.expect(is_close(decode_oct(encode_oct(dn)), dn, 1e-6)); + } +} + +test build_basis { + const n = normalize(@Vector(3, f32){ 0.3, -0.5, 0.8 }); + const basis = build_basis(n); + try std.testing.expectApproxEqAbs(@as(f32, 0), dot(basis.tangent, n), 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 0), dot(basis.bitangent, n), 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 0), dot(basis.tangent, basis.bitangent), 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 1), norm(basis.tangent), 1e-5); + try std.testing.expectApproxEqAbs(@as(f32, 1), norm(basis.bitangent), 1e-5); } test scale { @@ -165,32 +403,35 @@ test scale { try std.testing.expectEqual(@Vector(2, f32){ 0.5, 1 }, scale(v, 0.5)); } -test sw { +test is_close { + try std.testing.expect(is_close(@Vector(4, f32){ 1, 2, 3, 4 }, @Vector(4, f32){ 1.001, 2.001, 3.001, 4.001 }, 1.0e-4)); + try std.testing.expect(!is_close(@Vector(4, f32){ 1, 2, 3, 4 }, @Vector(4, f32){ 1.001, 2.001, 3.001, 4.001 }, 1.0e-6)); +} + +test swizzle { const v = @Vector(2, f32){ 1, 2 }; - try std.testing.expectEqual(@as(f32, 2), sw(v, "yx")[0]); - try std.testing.expectEqual(@as(f32, 1), sw(v, "yx")[1]); + try std.testing.expectEqual(@as(f32, 2), swizzle(v, "yx")[0]); + try std.testing.expectEqual(@as(f32, 1), swizzle(v, "yx")[1]); const v2 = @Vector(3, f32){ 1, 2, 3 }; const v2_expected = @Vector(3, f32){ 2, 3, 1 }; - try std.testing.expectEqual(v2_expected, sw(v2, "yzx")); + try std.testing.expectEqual(v2_expected, swizzle(v2, "yzx")); } -test angleAdv { +test angle_adv { const v1 = @Vector(2, f32){ 1, 0 }; const v2 = @Vector(2, f32){ 0, 1 }; const expected: f64 = @as(f64, std.math.pi) / @as(f64, 2); - try std.testing.expectApproxEqAbs(expected, angleAdv(v1, v2, 32), 0.0000001); - try std.testing.expectApproxEqAbs(expected, angleAdv(v1, v2, 64), 0.00000000001); + try std.testing.expectApproxEqAbs(expected, angle_adv(v1, v2, 32), 0.0000001); + try std.testing.expectApproxEqAbs(expected, angle_adv(v1, v2, 64), 0.00000000001); } test normalize { // Test with floating point vector - const v_f32 = @Vector(3, f32){ 3, 4, 0 }; - const normalized_f32 = normalize(v_f32); + const normalized_f32 = normalize(@Vector(3, f32){ 3, 4, 0 }); try std.testing.expectApproxEqAbs(@as(f32, 0.6), normalized_f32[0], 0.0001); try std.testing.expectApproxEqAbs(@as(f32, 0.8), normalized_f32[1], 0.0001); try std.testing.expectApproxEqAbs(@as(f32, 0), normalized_f32[2], 0.0001); try std.testing.expectApproxEqAbs(@as(f32, 1), norm(normalized_f32), 0.0001); - // Test with integer vector const v_i32 = @Vector(2, i32){ 3, 4 }; const normalized_i32 = normalize(v_i32); @@ -255,26 +496,26 @@ test reflect { try std.testing.expectApproxEqAbs(expected2[2], result[2], 0.0001); } -test normAdv { +test norm_adv { const v = @Vector(2, f64){ 1, 0 }; - try std.testing.expectEqual(@as(f32, 1), normAdv(v, 32)); + try std.testing.expectEqual(@as(f32, 1), norm_adv(v, 32)); const v2 = @Vector(2, f16){ 1, 2 }; - const result = normAdv(v2, 64); + const result = norm_adv(v2, 64); try std.testing.expect(@TypeOf(result) == f64); const expected = 2.2360679774997896964091736687; - try std.testing.expectEqual(@as(f16, expected), normAdv(v2, 16)); - try std.testing.expectEqual(@as(f32, expected), normAdv(v2, 32)); - try std.testing.expectEqual(@as(f64, expected), normAdv(v2, 64)); + try std.testing.expectEqual(@as(f16, expected), norm_adv(v2, 16)); + try std.testing.expectEqual(@as(f32, expected), norm_adv(v2, 32)); + try std.testing.expectEqual(@as(f64, expected), norm_adv(v2, 64)); const v3 = @Vector(2, i8){ 1, 2 }; - try std.testing.expect(@TypeOf(normAdv(v3, 64)) == f64); - try std.testing.expectEqual(@as(f16, expected), normAdv(v3, 16)); - try std.testing.expectEqual(@as(f32, expected), normAdv(v3, 32)); - try std.testing.expectEqual(@as(f64, expected), normAdv(v3, 64)); + try std.testing.expect(@TypeOf(norm_adv(v3, 64)) == f64); + try std.testing.expectEqual(@as(f16, expected), norm_adv(v3, 16)); + try std.testing.expectEqual(@as(f32, expected), norm_adv(v3, 32)); + try std.testing.expectEqual(@as(f64, expected), norm_adv(v3, 64)); const v4 = @Vector(2, u8){ 1, 2 }; - try std.testing.expect(@TypeOf(normAdv(v4, 64)) == f64); - try std.testing.expectEqual(@as(f16, expected), normAdv(v4, 16)); - try std.testing.expectEqual(@as(f32, expected), normAdv(v4, 32)); - try std.testing.expectEqual(@as(f64, expected), normAdv(v4, 64)); + try std.testing.expect(@TypeOf(norm_adv(v4, 64)) == f64); + try std.testing.expectEqual(@as(f16, expected), norm_adv(v4, 16)); + try std.testing.expectEqual(@as(f32, expected), norm_adv(v4, 32)); + try std.testing.expectEqual(@as(f64, expected), norm_adv(v4, 64)); } test norm { @@ -287,3 +528,189 @@ test norm { const v3 = @Vector(4, i32){ 1, 2, 2, 0 }; try std.testing.expectApproxEqAbs(@as(f32, 3), norm(v3), 0.0001); } + +test sin_cos { + { + const input = @Vector(2, f32){ 0, std.math.pi * 0.5 }; + const res = sin_cos(input); + try testing.expect_is_close(res.cos_out, @Vector(2, f32){ 1, 0 }, 0.0001); + try testing.expect_is_close(res.sin_out, @Vector(2, f32){ 0, 1.0 }, 0.0001); + } + { + const input = @Vector(4, f32){ 0, std.math.pi * 0.5, std.math.pi, -0.5 * std.math.pi }; + const res = sin_cos(input); + try testing.expect_is_close(res.cos_out, @Vector(4, f32){ 1, 0, -1, 0 }, 0.0001); + try testing.expect_is_close(res.sin_out, @Vector(4, f32){ 0, 1.0, 0, -1.0 }, 0.0001); + } + var ms: f64 = 0.0; + var mc: f64 = 0.0; + + var x: f32 = -100.0 * std.math.pi; + while (x < 100.0 * std.math.pi) : (x += 1.0e-3) { + const xv = @as(Vec4f32, @splat(x)) + Vec4f32{ 0.0e-4, 2.5e-4, 5.0e-4, 7.5e-4 }; + const res2 = sin_cos(xv); + inline for (0..3) |i| { + const s1 = std.math.sin(xv[i]); + const s2 = res2.sin_out[i]; + ms = @max(ms, @abs(s1 - s2)); + + const c1 = std.math.cos(xv[i]); + const c2 = res2.cos_out[i]; + mc = @max(mc, @abs(c1 - c2)); + } + } + try std.testing.expect(ms < 1.0e-7); + try std.testing.expect(mc < 1.0e-7); +} + +test "array inputs are accepted and preserve shape" { + // Same-length vector-returning: array in -> array out. + { + const result = normalize([3]f32{ 3, 4, 0 }); + try std.testing.expect(@TypeOf(result) == [3]f32); + try std.testing.expectApproxEqAbs(@as(f32, 0.6), result[0], 0.0001); + try std.testing.expectApproxEqAbs(@as(f32, 0.8), result[1], 0.0001); + try std.testing.expectApproxEqAbs(@as(f32, 0), result[2], 0.0001); + } + // @Vector in -> @Vector out (unchanged behavior). + { + const result = normalize(@Vector(3, f32){ 3, 4, 0 }); + try std.testing.expect(@TypeOf(result) == @Vector(3, f32)); + } + // Scalar-returning: container-independent. + { + const d = dot([3]f32{ 1, 2, 3 }, [3]f32{ 4, 5, 6 }); + try std.testing.expect(@TypeOf(d) == f32); + try std.testing.expectEqual(@as(f32, 32), d); + try std.testing.expectApproxEqAbs(@as(f32, 5), norm([2]f32{ 3, 4 }), 0.0001); + } + // cross: array in -> array out. + { + const result = cross([3]f32{ 1, 0, 0 }, [3]f32{ 0, 1, 0 }); + try std.testing.expect(@TypeOf(result) == [3]f32); + try std.testing.expectApproxEqAbs(@as(f32, 1), result[2], 0.0001); + } + // reflect: array in -> array out. + { + const result = reflect([2]f32{ 1, -1 }, [2]f32{ 0, 1 }); + try std.testing.expect(@TypeOf(result) == [2]f32); + try std.testing.expectApproxEqAbs(@as(f32, 1), result[1], 0.0001); + } + // scale: array in -> array out. + { + const result = scale([2]f32{ 1, 2 }, 2); + try std.testing.expect(@TypeOf(result) == [2]f32); + try std.testing.expectApproxEqAbs(@as(f32, 4), result[1], 0.0001); + } + // Length-changing: extract/swizzle preserve the container kind. + { + const e = extract([4]f32{ 1, 2, 3, 4 }, 2); + try std.testing.expect(@TypeOf(e) == [2]f32); + try std.testing.expectEqual(@as(f32, 1), e[0]); + try std.testing.expectEqual(@as(f32, 2), e[1]); + + const s = swizzle([2]f32{ 1, 2 }, "yx"); + try std.testing.expect(@TypeOf(s) == [2]f32); + try std.testing.expectEqual(@as(f32, 2), s[0]); + try std.testing.expectEqual(@as(f32, 1), s[1]); + } + // distance / angle over arrays. + { + try std.testing.expectApproxEqAbs(@as(f32, 5), distance([2]f32{ 1, 1 }, [2]f32{ 4, 5 }), 0.0001); + try std.testing.expectApproxEqAbs(@as(f32, std.math.pi / 2.0), angle([2]f32{ 1, 0 }, [2]f32{ 0, 1 }), 0.0001); + } + // is_close over arrays. + { + try std.testing.expect(is_close([2]f32{ 1, 2 }, [2]f32{ 1.0001, 2.0001 }, 1.0e-4)); + } + // sin_cos: array in -> array out. + { + const res = sin_cos([2]f32{ 0, std.math.pi * 0.5 }); + try std.testing.expect(@TypeOf(res.sin_out) == [2]f32); + try testing.expect_is_close(@as(@Vector(2, f32), res.cos_out), @Vector(2, f32){ 1, 0 }, 0.0001); + try testing.expect_is_close(@as(@Vector(2, f32), res.sin_out), @Vector(2, f32){ 0, 1.0 }, 0.0001); + } +} + +test "chunked array path matches single-@Vector path" { + @setEvalBranchQuota(10000); + // For each length (incl. remainder-forcing and the motivating 245), the + // chunked array kernels must agree with the single-op @Vector path on the + // same data. Pure element-wise ops (scale, sin_cos) must match exactly; + // reductions and reduction-dependent maps agree to within float reordering. + // Inner fill/compare loops are RUNTIME `for` (not `inline`) so the unrolled + // outer length loop doesn't blow the comptime branch quota. + const lengths = [_]comptime_int{ 3, 8, 16, 17, 64, 65, 77, 245 }; + inline for (lengths) |N| { + // Strictly-positive, moderate f32 data so reduction magnitudes stay nonzero. + var a: [N]f32 = undefined; + var b: [N]f32 = undefined; + for (0..N) |i| { + a[i] = @as(f32, @floatFromInt((i % 13) + 1)) * 0.5; // 0.5 .. 6.5 + b[i] = @as(f32, @floatFromInt((i % 7) + 1)) * 0.25; // 0.25 .. 1.75 + } + const va: @Vector(N, f32) = a; + const vb: @Vector(N, f32) = b; + + // Reductions -> scalar: relative tolerance absorbs the sum reordering. + try std.testing.expectApproxEqRel(norm_sqr(va), norm_sqr(a), 1e-4); + try std.testing.expectApproxEqRel(norm(va), norm(a), 1e-4); + try std.testing.expectApproxEqRel(dot(va, vb), dot(a, b), 1e-4); + try std.testing.expectApproxEqRel(distance(va, vb), distance(a, b), 1e-4); + + // scale: pure per-lane multiply -> bit-identical. + try std.testing.expectEqual(@as([N]f32, scale(va, 2.0)), scale(a, 2.0)); + + // sin_cos: pure per-lane transcendental -> bit-identical. + const sc_v = sin_cos(va); + const sc_a = sin_cos(a); + try std.testing.expectEqual(@as([N]f32, sc_v.sin_out), sc_a.sin_out); + try std.testing.expectEqual(@as([N]f32, sc_v.cos_out), sc_a.cos_out); + + // normalize: components are O(1/sqrt(N)); only the norm scalar reorders. + const nz_v: [N]f32 = normalize(va); + const nz_a = normalize(a); + for (0..N) |i| try std.testing.expectApproxEqAbs(nz_v[i], nz_a[i], 1e-4); + + // reflect about a UNIT normal keeps magnitudes O(input) so an abs tol works. + const unit_v: [N]f32 = normalize(vb); + const unit_a = normalize(b); + const rf_v: [N]f32 = reflect(va, @as(@Vector(N, f32), unit_v)); + const rf_a = reflect(a, unit_a); + for (0..N) |i| try std.testing.expectApproxEqAbs(rf_v[i], rf_a[i], 1e-3); + + // is_close over arrays agrees with the vector path. + try std.testing.expectEqual(is_close(va, vb, 1000.0), is_close(a, b, 1000.0)); + } + + // Integer reductions/maps are exact (integer add is associative mod 2^bits). + inline for (.{ 3, 17, 64, 245 }) |N| { + var ai: [N]i32 = undefined; + var bi: [N]i32 = undefined; + for (0..N) |i| { + ai[i] = @intCast((i % 5) + 1); + bi[i] = @intCast((i % 3) + 1); + } + const vai: @Vector(N, i32) = ai; + const vbi: @Vector(N, i32) = bi; + try std.testing.expectEqual(norm_sqr(vai), norm_sqr(ai)); + try std.testing.expectEqual(dot(vai, vbi), dot(ai, bi)); + try std.testing.expectEqual(@as([N]i32, normalize(vai)), normalize(ai)); + } + + // f16 with precision-64 widening: both paths accumulate in f64, so they match. + inline for (.{ 17, 64 }) |N| { + var af: [N]f16 = undefined; + for (0..N) |i| af[i] = @floatFromInt((i % 4) + 1); + const vaf: @Vector(N, f16) = af; + try std.testing.expectApproxEqRel(norm_adv(vaf, 64), norm_adv(af, 64), 1e-6); + } + + // f64 reductions. + inline for (.{ 64, 245 }) |N| { + var ad: [N]f64 = undefined; + for (0..N) |i| ad[i] = @as(f64, @floatFromInt((i % 11) + 1)) * 0.3; + const vad: @Vector(N, f64) = ad; + try std.testing.expectApproxEqRel(norm(vad), norm(ad), 1e-12); + } +}