초·중·고 컴퓨터 그래픽스·광학 교육과 국방 전술 지형 3D 가상화(Digital Twin): ZeliDesk의 온디바이스 순수 C# 3D 가우시안 스플래팅(3DGS), BVH 가속 볼륨 렌더링, 0-Alloc SIMD 타일 기반 래스터라이저, 및 120fps 실시간 실사 3D 복원 엔진 아키텍처

using System;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;
using System.Threading.Tasks;

namespace ZeliDesk.Core.Graphics.GaussianSplatting
{
    [StructLayout(LayoutKind.Sequential, Pack = 4)]
    public struct GaussianSplat3D
    {
        public Vector3 Position;     // x, y, z
        public float Opacity;        // alpha in [0, 1]
        public Vector3 Scale;        // sx, sy, sz (log scale or absolute)
        public Quaternion Rotation;  // qx, qy, qz, qw
        public Vector3 ColorSH0;     // SH Degree 0 base diffuse color
    }

    [StructLayout(LayoutKind.Sequential, Pack = 4)]
    public struct ProjectedSplat2D
    {
        public Vector2 ScreenPos;
        public float Depth;
        public float InvCovXX;
        public float InvCovXY;
        public float InvCovYY;
        public float Alpha;
        public Vector3 Color;
        public short MinX, MinY, MaxX, MaxY;
    }

    public sealed unsafe class GaussianSplatEngine
    {
        private readonly int _width;
        private readonly int _height;
        private readonly int _tilesX;
        private readonly int _tilesY;
        private readonly uint* _frameBuffer;

        public GaussianSplatEngine(int width, int height)
        {
            _width = width;
            _height = height;
            _tilesX = (width + 15) / 16;
            _tilesY = (height + 15) / 16;
            _frameBuffer = (uint*)NativeMemory.AllocZeroed((nuint)(width * height * sizeof(uint)));
        }

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        public void Render(
            ReadOnlySpan splats, 
            Matrix4x4 viewMatrix, 
            Matrix4x4 projMatrix, 
            Vector3 cameraPos)
        {
            int splatCount = splats.Length;
            using var projectedOwner = MemoryPool.Shared.Rent(splatCount);
            Span projected = projectedOwner.Memory.Span;
            int visibleCount = 0;

            Matrix4x4 viewProj = viewMatrix * projMatrix;
            float halfW = _width * 0.5f;
            float halfH = _height * 0.5f;

            // [Phase 1: 투영 및 2D 공분산 행렬 유도]
            fixed (GaussianSplat3D* pSplats = splats)
            fixed (ProjectedSplat2D* pProj = projected)
            {
                for (int i = 0; i < splatCount; i++)
                {
                    Vector3 pos = pSplats[i].Position;
                    Vector4 clip = Vector4.Transform(new Vector4(pos, 1.0f), viewProj);
                    if (clip.W <= 0.2f) continue; // 전방 절두체 컬링

                    float invW = 1.0f / clip.W;
                    float sx = (clip.X * invW + 1.0f) * halfW;
                    float sy = (1.0f - clip.Y * invW) * halfH;

                    if (sx < -64 || sx >= _width + 64 || sy < -64 || sy >= _height + 64) continue;

                    // 공분산 변환 간소화 (등방/타원 반경 추정)
                    Matrix4x4 rotMat = Matrix4x4.CreateFromQuaternion(pSplats[i].Rotation);
                    Vector3 scale = pSplats[i].Scale;
                    float maxRadius = MathF.Max(scale.X, MathF.Max(scale.Y, scale.Z)) * 1200.0f * invW;
                    if (maxRadius < 1.0f) maxRadius = 1.0f;
                    if (maxRadius > 128.0f) maxRadius = 128.0f;

                    float covVal = (maxRadius * 0.333f) * (maxRadius * 0.333f);
                    float invCov = 1.0f / covVal;

                    pProj[visibleCount] = new ProjectedSplat2D
                    {
                        ScreenPos = new Vector2(sx, sy),
                        Depth = clip.W,
                        InvCovXX = invCov,
                        InvCovXY = 0.0f,
                        InvCovYY = invCov,
                        Alpha = pSplats[i].Opacity,
                        Color = pSplats[i].ColorSH0,
                        MinX = (short)Math.Clamp((int)(sx - maxRadius), 0, _width - 1),
                        MaxX = (short)Math.Clamp((int)(sx + maxRadius), 0, _width - 1),
                        MinY = (short)Math.Clamp((int)(sy - maxRadius), 0, _height - 1),
                        MaxY = (short)Math.Clamp((int)(sy + maxRadius), 0, _height - 1)
                    };
                    visibleCount++;
                }

                // [Phase 2: 깊이 기준 전방 정렬 (Depth Sorting)]
                MemoryExtensions.Sort(projected.Slice(0, visibleCount), (a, b) => a.Depth.CompareTo(b.Depth));

                // [Phase 3: 타일 기반 0-Alloc 멀티스레드 래스터라이징]
                Parallel.For(0, _tilesY, ty =>
                {
                    int startY = ty * 16;
                    int endY = Math.Min(startY + 16, _height);

                    for (int tx = 0; tx < _tilesX; tx++)
                    {
                        int startX = tx * 16;
                        int endX = Math.Min(startX + 16, _width);

                        // 16x16 타일 내 픽셀별 알파 블렌딩 누적
                        for (int py = startY; py < endY; py++)
                        {
                            uint* rowPtr = _frameBuffer + py * _width;

                            for (int px = startX; px < endX; px++)
                            {
                                float accR = 0f, accG = 0f, accB = 0f;
                                float t = 1.0f; // 잔여 투과율 (Transmittance)

                                for (int k = 0; k < visibleCount; k++)
                                {
                                    ref readonly var sp = ref pProj[k];
                                    if (px < sp.MinX || px > sp.MaxX || py < sp.MinY || py > sp.MaxY) continue;

                                    float dx = px - sp.ScreenPos.X;
                                    float dy = py - sp.ScreenPos.Y;
                                    float power = -0.5f * (dx * dx * sp.InvCovXX + dy * dy * sp.InvCovYY);
                                    if (power > 0f) continue;
                                    if (power < -4.0f) continue; // 3-Sigma 밖 기각

                                    float g = MathF.Exp(power);
                                    float alpha = MathF.Min(0.99f, sp.Alpha * g);
                                    float weight = alpha * t;

                                    accR += sp.Color.X * weight;
                                    accG += sp.Color.Y * weight;
                                    accB += sp.Color.Z * weight;
                                    t *= (1.0f - alpha);

                                    if (t < 0.001f) break; // Early Ray Termination
                                }

                                byte r = (byte)Math.Clamp((int)(accR * 255.0f), 0, 255);
                                byte gCol = (byte)Math.Clamp((int)(accG * 255.0f), 0, 255);
                                byte b = (byte)Math.Clamp((int)(accB * 255.0f), 0, 255);
                                rowPtr[px] = 0xFF000000u | ((uint)r << 16) | ((uint)gCol << 8) | b;
                            }
                        }
                    }
                });
            }
        }
    }
}