acdream/src/AcDream.App/Rendering/FrustumCuller.cs

using System.Numerics;

namespace AcDream.App.Rendering;

/// <summary>
/// Six normalized view-frustum planes extracted from a View×Projection matrix.
/// Each plane is represented as (normal.X, normal.Y, normal.Z, distance) where
/// dot(normal, point) + distance >= 0 means the point is on the visible side.
/// </summary>
public readonly struct FrustumPlanes
{
    public readonly Vector4 Left;
    public readonly Vector4 Right;
    public readonly Vector4 Bottom;
    public readonly Vector4 Top;
    public readonly Vector4 Near;
    public readonly Vector4 Far;

    private FrustumPlanes(Vector4 left, Vector4 right, Vector4 bottom, Vector4 top, Vector4 near, Vector4 far)
    {
        Left   = left;
        Right  = right;
        Bottom = bottom;
        Top    = top;
        Near   = near;
        Far    = far;
    }

    /// <summary>
    /// Extracts the six frustum planes from a combined View×Projection matrix
    /// using the Gribb-Hartmann method. System.Numerics.Matrix4x4 is row-major,
    /// so rows are accessed directly via M{row}{col} fields.
    /// </summary>
    public static FrustumPlanes FromViewProjection(Matrix4x4 vp)
    {
        // System.Numerics.Matrix4x4 uses ROW-VECTOR convention: clip = worldPos * VP
        // So clip.x = dot(worldPos, col1), clip.w = dot(worldPos, col4), etc.
        // Gribb-Hartmann left plane: clip.x + clip.w >= 0 => dot(worldPos, col1 + col4) >= 0
        // Therefore we must operate on COLUMNS of the matrix (not rows).
        // Column vectors (M{row}{col}, so col1 has elements from col-index=1):
        var col1 = new Vector4(vp.M11, vp.M21, vp.M31, vp.M41);
        var col2 = new Vector4(vp.M12, vp.M22, vp.M32, vp.M42);
        var col3 = new Vector4(vp.M13, vp.M23, vp.M33, vp.M43);
        var col4 = new Vector4(vp.M14, vp.M24, vp.M34, vp.M44);

        // Gribb-Hartmann extraction for row-vector (System.Numerics) convention:
        var left   = Normalize(col4 + col1);
        var right  = Normalize(col4 - col1);
        var bottom = Normalize(col4 + col2);
        var top    = Normalize(col4 - col2);
        // NEAR is the one plane whose formula depends on the clip-space z range.
        // The familiar `col4 + col3` is the OpenGL form, for NDC z in [-1,1],
        // where the near plane is the locus of clip.z = -clip.w. Every acdream
        // projection comes from Matrix4x4.CreatePerspectiveFieldOfView (and
        // CreateOrthographic), whose NDC z range is [0,1] — there the near plane
        // is clip.z = 0, which is `col3` by itself. Using the GL form here put
        // the effective threshold at -n·f/(2f-n), about half the true near
        // distance, which only ever kept geometry the real frustum would drop.
        // See docs/ISSUES.md #248; same class as the PortalProjection near-test
        // bug at PortalProjection.cs:12-19. FAR is `col4 - col3` under BOTH
        // conventions and is deliberately left alone.
        var near   = Normalize(col3);
        var far    = Normalize(col4 - col3);

        return new FrustumPlanes(left, right, bottom, top, near, far);
    }

    private static Vector4 Normalize(Vector4 plane)
    {
        float length = MathF.Sqrt(plane.X * plane.X + plane.Y * plane.Y + plane.Z * plane.Z);
        return plane / length;
    }
}

/// <summary>
/// Conservative AABB-vs-frustum culling. Zero allocations; suitable for per-frame use.
/// </summary>
public static class FrustumCuller
{
    /// <summary>
    /// Returns true if the axis-aligned bounding box defined by
    /// [min, max] is potentially visible against the given frustum
    /// planes. Conservative: returns true for partial intersections.
    /// Zero allocations; suitable for per-frame use.
    /// </summary>
    public static bool IsAabbVisible(FrustumPlanes planes, Vector3 min, Vector3 max)
    {
        // For each plane, test the AABB's "most-positive vertex" —
        // the corner most in the direction of the plane normal. If
        // that corner is behind the plane, the entire box is outside.
        return TestPlane(planes.Left,   min, max)
            && TestPlane(planes.Right,  min, max)
            && TestPlane(planes.Bottom, min, max)
            && TestPlane(planes.Top,    min, max)
            && TestPlane(planes.Near,   min, max)
            && TestPlane(planes.Far,    min, max);
    }

    private static bool TestPlane(Vector4 plane, Vector3 min, Vector3 max)
    {
        // Pick the corner of the AABB most in the direction of the
        // plane normal (the "positive vertex").
        float px = plane.X >= 0 ? max.X : min.X;
        float py = plane.Y >= 0 ? max.Y : min.Y;
        float pz = plane.Z >= 0 ? max.Z : min.Z;

        // If the positive vertex is behind the plane, the box is
        // fully outside this half-space.
        return plane.X * px + plane.Y * py + plane.Z * pz + plane.W >= 0;
    }
}