// Geometry functions
// Copyright (C) 2002  David Rosen
// Copyright (C) 2023, 2025  Nguyễn Gia Phong
//
// This file is part of Black Shades.
//
// Black Shades is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// Black Shades is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Black Shades.  If not, see <https://www.gnu.org/licenses/>.

const Child = std.meta.Child;
const degreesToRadians = std.math.degreesToRadians;
const floatEps = std.math.floatEps;
const std = @import("std");

pub fn sqr(x: anytype) @TypeOf(x) {
    return x * x;
}

fn dot(u: anytype, v: @TypeOf(u)) Child(@TypeOf(u)) {
    return @reduce(.Add, u * v);
}

pub const XYZ = extern struct { x: f32, y: f32, z: f32 };

export fn sqrlen(v: XYZ) f32 {
    const u: @Vector(3, f32) = @bitCast(v);
    return dot(u, u);
}

pub fn norm(v: anytype) Child(@TypeOf(v)) {
    return @sqrt(dot(v, v));
}

export fn len(v: XYZ) f32 {
    const u: @Vector(3, f32) = @bitCast(v);
    return norm(u);
}

pub export fn crossProduct(u: XYZ, v: XYZ) XYZ {
    return .{
        .x = u.y * v.z - u.z * v.y,
        .y = u.z * v.x - u.x * v.z,
        .z = u.x * v.y - u.y * v.x,
    };
}

pub fn splat(comptime n: comptime_int, x: anytype) @Vector(n, @TypeOf(x)) {
    return @splat(x);
}

pub export fn normalize(v: XYZ) XYZ {
    const u: @Vector(3, f32) = @bitCast(v);
    const d = norm(u);
    return if (d == 0) v else @bitCast(u / splat(3, d));
}

export fn reflect(v: XYZ, n: XYZ) XYZ {
    const u: @Vector(3, f32) = @bitCast(v);
    const m: @Vector(3, f32) = @bitCast(n);
    return @bitCast(u - m * splat(3, dot(u, m) * 2));
}

pub fn rotate2d(i: *f32, j: *f32, a: f32) void {
    if (a == 0) return;
    const x = i.*;
    const y = j.*;
    i.* = x * @cos(a) - y * @sin(a);
    j.* = x * @sin(a) + y * @cos(a);
}

export fn rotate(v: XYZ, deg_x: f32, deg_y: f32, deg_z: f32) XYZ {
    var u = v;
    // TODO: optimize
    rotate2d(&u.x, &u.y, degreesToRadians(deg_z));
    rotate2d(&u.z, &u.x, degreesToRadians(deg_y));
    rotate2d(&u.y, &u.z, degreesToRadians(deg_x));
    return u;
}

pub export fn segCrossSphere(a: XYZ, b: XYZ, i: XYZ, r: f32) bool {
    const p: @Vector(3, f32) = @bitCast(a);
    const q: @Vector(3, f32) = @bitCast(b);
    const c: @Vector(3, f32) = @bitCast(i);
    if (@reduce(.Or, @max(p, q) < c - splat(3, r))) return false;
    if (@reduce(.Or, @min(p, q) > c + splat(3, r))) return false;
    // https://en.wikipedia.org/wiki/Line–sphere_intersection
    const d = q - p; // line's direction
    const u = d / splat(3, norm(d)); // unit vector
    return sqr(dot(u, (p - c))) >= @reduce(.Add, sqr(p - c)) - sqr(r);
}

pub export fn segCrossTrigon(start: XYZ, end: XYZ,
                             p_a: *const XYZ, p_b: *const XYZ, p_c: *const XYZ,
                             normal: *const XYZ, intersection: *XYZ) bool {
    const p: @Vector(3, f32) = @bitCast(start);
    const q: @Vector(3, f32) = @bitCast(end);
    const n: @Vector(3, f32) = @bitCast(normal.*);
    const denom = dot(q - p, n);
    if (@abs(denom) < floatEps(f32))
        return false; // parallel segment and triangle

    const a: @Vector(3, f32) = @bitCast(p_a.*);
    const mu = (dot(a, n) - dot(p, n)) / denom;
    const i = p + (q - p) * splat(3, mu);
    if (mu < 0 or mu > 1)
        return false; // intersection not within segment

    // Check if intersection is in the triangle
    const n_abs = @abs(n);
    const n_max = @reduce(.Max, n_abs);
    const k: struct { usize, usize } = if (n_max == n_abs[0])
        .{ 1, 2 }
    else if (n_max == n_abs[1])
        .{ 0, 2 }
    else if (n_max == n_abs[2])
        .{ 0, 1 }
    else unreachable;

    const b: @Vector(3, f32) = @bitCast(p_b.*);
    const c: @Vector(3, f32) = @bitCast(p_c.*);
    const u = @Vector(3, f32){ i[k[0]], b[k[0]], c[k[0]] } - splat(3, a[k[0]]);
    const v = @Vector(3, f32){ i[k[1]], b[k[1]], c[k[1]] } - splat(3, a[k[1]]);
    intersection.* = @bitCast(i);

    if (@abs(u[1]) < floatEps(f32)) {
        const s = u[0] / u[2];
        if (s >= 0 and s <= 1) {
            const t = (v[0] - s * v[2]) / v[1];
            if (t >= 0 and s + t <= 1)
                return true;
        }
    } else {
        const s = (v[0] * u[1] - u[0] * v[1]) / (v[2] * u[1] - u[2] * v[1]);
        if (s >= 0 and s <= 1) {
            const t = (u[0] - s * u[2]) / u[1];
            if (t >= 0 and s + t <= 1)
                return true;
        }
    }
    return false;
}

fn transpose(comptime n: comptime_int, m: [n]@Vector(n, f32)) @TypeOf(m) {
    const flat: @Vector(sqr(n), f32) = @bitCast(m);
    return @bitCast(@shuffle(f32, flat, undefined, blk: {
        var v: @Vector(sqr(n), i32) = undefined;
        for (0..n) |i| {
            for (0..n) |j|
                v[i * n + j] = @intCast(j * n + i);
        }
        break :blk v;
    }));
}

export fn setFrustum(frustum: *[6]@Vector(4, f32),
                     p: *const [4]@Vector(4, f32),
                     mv: *const [4]@Vector(4, f32)) void {
    var mvp: [4]@Vector(4, f32) = undefined;
    for (&mvp, transpose(4, p.*)) |*u, p_col| {
        for (mv, 0..4) |mv_row, i|
            u[i] = dot(p_col, mv_row);
    } // matrix multiplication

    frustum.* = .{
        mvp[3] - mvp[0], mvp[3] + mvp[0], // right & left planes
        mvp[3] - mvp[1], mvp[3] + mvp[1], // bottom & top planes
        mvp[3] - mvp[2], mvp[3] + mvp[2], // far & near planes
    };
    for (frustum) |*plane| // normalize
        plane.* /= @splat(norm(plane.* * @Vector(4, f32){ 1, 1, 1, 0 }));
}

export fn cubeInFrustum(frustum: *const [6]@Vector(4, f32),
                        x: f32, y: f32, z: f32, size: f32) bool {
    const delta = [_]f32{ -size, size };
    loop: for (frustum) |*plane| {
        for (delta) |dx| for (delta) |dy| for (delta) |dz|
            if (dot(plane.*, @Vector(4, f32){ x + dx, y + dy, z + dz, 1 }) > 0)
                continue :loop;
        return false;
    }
    return true;
}

export fn sphereInFrustum(frustum: *const [6]@Vector(4, f32),
                          x: f32, y: f32, z: f32, r: f32) bool {
    for (frustum) |*plane|
        if (dot(plane.*, @Vector(4, f32){ x, y, z, 1 }) <= -r)
            return false;
    return true;
}