feat: add more math functions

Signed-off-by: Michael Pollind <mpollind@gmail.com>

Author: pollend

README.md

@@ -1,4 +1,4 @@
-# Zig Linear Algebra (ZLA)
+# Zig Math Library (ZML)
 
 [![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)
 
@@ -13,11 +13,11 @@ zig fetch --save "git+https://github.com/flying-swallow/zig-linear-algebra.git"
 Then in your `build.zig`:
 
 ```zig
-const zla = b.dependency("zla", .{
+const zml = b.dependency("zml", .{
     .target = target,
     .optimize = optimize,
 });
 
-exe.root_module.addImport("zla", zla.module("zla"));
+exe.root_module.addImport("zml", zml.module("zml"));
 ```
 

bench/bench.zig

@@ -1,10 +1,10 @@
 const std = @import("std");
 const zbench = @import("zbench");
-const zla = @import("zla");
+const zml = @import("zml");
 
 fn benchmark_multiply(comptime size: usize) type {
     return struct {
-        const Matrix512 = zla.Mat(f32, size, size);
+        const Matrix512 = zml.Mat(f32, size, size);
         a: Matrix512,
         b: Matrix512,
 
@@ -40,7 +40,7 @@ fn bench_sin_cos_fused(comptime size: usize) type {
         }
 
         pub fn run(self: *@This(), _: std.mem.Allocator) void {
-           std.mem.doNotOptimizeAway(@call(.never_inline, zla.vec.sin_cos, .{self.angles}));
+           std.mem.doNotOptimizeAway(@call(.never_inline, zml.vec.sin_cos, .{self.angles}));
         }
     };
 }

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,
@@ -30,7 +30,7 @@ pub fn build(b: *std.Build) void {
                     .target = target,
                     .optimize = optimize,
                     .imports = &.{
-                        .{ .name = "zla", .module = root_module },
+                        .{ .name = "zml", .module = root_module },
                         .{ .name = "zbench", .module = zbench_dep.module("zbench") },
                     },
                 }),

build.zig.zon

@@ -1,7 +1,7 @@
 .{
-    .name = .zla,
-    .version = "0.0.2",
-    .fingerprint = 0x55da41f1f0b5a9ce, // Changing this has security and trust implications.
+    .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",

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);
+    }
+}

src/geometry.zig

@@ -7,6 +7,14 @@ pub const Primative = enum {
     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;
@@ -20,6 +28,9 @@ 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);
 
 
 
@@ -33,6 +44,7 @@ test {
     _ = @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");

src/geometry/aabb.zig

@@ -1,6 +1,6 @@
 const std = @import("std");
-const zla = @import("../root.zig");
-const geometry = zla.geom;
+const zml = @import("../root.zig");
+const geometry = zml.geom;
 
 pub fn AABB(comptime T: type) type {
     return struct {
@@ -34,7 +34,7 @@ pub fn AABB(comptime T: type) type {
         }
 
         pub fn get_sqr_distance_to(self: Self, point: @Vector(3, T)) T {
-            return zla.vec.norm_sqr(self.get_closest_point(point) - point);
+            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() {
@@ -53,7 +53,7 @@ pub fn AABB(comptime T: type) type {
             };
         }
 
-        pub fn transform(self: Self, mat: zla.Mat(T, 4, 4)) Self {
+        pub fn transform(self: Self, mat: zml.Mat(T, 4, 4)) Self {
             const pos = mat.position();
             var new_min = pos;
             var new_max = pos;
@@ -107,7 +107,7 @@ pub fn AABB(comptime T: type) type {
         pub fn surface_area(self: Self) T {
             std.debug.assert(@reduce(.And, self.min <= self.max));
             const extent = self.max - self.min;
-            return 2 * zla.vec.dot(zla.vec.swizzle(extent, "yzz"), zla.vec.swizzle(extent, "xxy"));
+            return 2 * zml.vec.dot(zml.vec.swizzle(extent, "yzz"), zml.vec.swizzle(extent, "xxy"));
         }
 
         pub fn volume(self: Self) T {
@@ -121,11 +121,11 @@ pub fn AABB(comptime T: type) type {
 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: zla.Mat(f32, 4, 4) = .identity;
+    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(zla.vec.is_close_default(transformed.min, .{ 1.0, 2.0, 3.0 }));
-    try std.testing.expect(zla.vec.is_close_default(transformed.max, .{ 2.0, 3.0, 4.0 }));
+    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" {
@@ -153,14 +153,14 @@ test "InvDirection" {
 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(zla.vec.is_close_default(aabb.get_support(.{ 0.5774, 0.5774, 0.5774 }), .{ 0.0, 0.0, 0.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ 0.5774, 0.5774, -0.5774 }), .{ 0.0, 0.0, 1.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ 0.5774, -0.5774, 0.5774 }), .{ 0.0, 1.0, 0.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ 0.5774, -0.5774, -0.5774 }), .{ 0.0, 1.0, 1.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ -0.5774, 0.5774, 0.5774 }), .{ 1.0, 0.0, 0.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ -0.5774, 0.5774, -0.5774 }), .{ 1.0, 0.0, 1.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ -0.5774, -0.5774, 0.5774 }), .{ 1.0, 1.0, 0.0 }));
-    try std.testing.expect(zla.vec.is_close_default(aabb.get_support(.{ -0.5774, -0.5774, -0.5774 }), .{ 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 {

src/geometry/buffer_view.zig

@@ -1,16 +1,9 @@
 const std = @import("std");
+const meta = @import("../meta.zig");
 
 pub fn BufferView(comptime T: type) type {
     return struct {
-        pub const inner_type = switch (@typeInfo(T)) { 
-            .vector => |v| [v.len]std.meta.Child(T),
-            .array => |a| [a.len]std.meta.Child(T),
-            else => @compileError("Expected a vector or array type, got: " ++ @typeName(T)),
-        };
-        pub const num_elements = switch (@typeInfo(T)) {
-            .vector => |v| v.len,
-            .array => |a| a.len,
-            else => @compileError("Expected a vector or array type, got: " ++ @typeName(T)),
-        };
+        pub const inner_type = [meta.lengthOf(T)]meta.Child(T);
+        pub const num_elements = meta.lengthOf(T);
     };
 }

src/geometry/capsule.zig

@@ -1,5 +1,5 @@
 const geometry = @import("../geometry.zig");
-const zla = @import("../root.zig");
+const zml = @import("../root.zig");
 
 pub fn Capsule(comptime T: type) type {
     return struct {
@@ -15,7 +15,7 @@ pub fn Capsule(comptime T: type) type {
         }
 
         pub fn get_cylinder_height(self: Self) T {
-            return zla.vec.distance(self.hemisphere_centers[0], self.hemisphere_centers[1]);
+            return zml.vec.distance(self.hemisphere_centers[0], self.hemisphere_centers[1]);
         }
 
         pub fn get_total_height(self: Self) T {

src/geometry/contains.zig

@@ -1,5 +1,5 @@
 const std = @import("std");
-const zla = @import("../root.zig");
+const zml = @import("../root.zig");
 
 pub fn aabb_contains_point(a: anytype, pt: @Vector(3, @TypeOf(a).child)) bool {
     comptime {
@@ -13,21 +13,21 @@ pub fn aabb_contains_point(a: anytype, pt: @Vector(3, @TypeOf(a).child)) bool {
 //        std.debug.assert(@TypeOf(a).primative_type == .Capsule);
 //    }
 //    const ab = a.hemisphere_centers[1] - a.hemisphere_centers[0];
-//    const t = @max(@min(zla.vec.dot(pt - a.hemisphere_centers[0], ab) / zla.vec.dot(ab, ab), @as(@TypeOf(a).child, 1)), @as(@TypeOf(a).child, 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 zla.vec.distance_sqr(pt, closest_point) <= a.radius * a.radius;
+//    return zml.vec.distance_sqr(pt, closest_point) <= a.radius * a.radius;
 //} 
 
 test aabb_contains_point {
-    const aabb: zla.geom.AABB(f32) = .from_two_points(.{ -1, -1, -1 }, .{ 1, 1, 1 });
+    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: zla.geom.Capsule(f32) = .{
+//    const capsule: zml.geom.Capsule(f32) = .{
 //        .hemisphere_centers = .{ .{ 0, 0, -1 }, .{ 0, 0, 1 } },
 //        .radius = 1,
 //    };