XFEEDToolBox

using System.Buffers.Binary;
using System.IO.Compression;
using System.Runtime.InteropServices;
using System.Security.Cryptography;
using System.Text.Json;
using XFEEDToolBox.Core.Models;

namespace XFEEDToolBox.Core.Rendering;

/// <summary>The original integer GDM inputs, not a normalized overview image.</summary>
public sealed class GalaxyOriginalDensityMap
{
    private readonly uint[] massDensity;
    private readonly ushort[] dustSurface,milkyWaySurface,radialHeight;
    private readonly uint[] coefficients;
    public const int Width=2048,Height=2048;
    public const string PayloadSha256="55aa2bc69d9ae86a4262145ab5346a6212a0f77f9929debe228cae3c5b1298c5";
    public Coordinates Origin {get;}=new(49985,40985,24105);
    public IReadOnlyList<uint> Coefficients=>Array.AsReadOnly(coefficients);

    internal GalaxyOriginalDensityMap(JsonElement metadata,byte[] payload)
    {
        if(payload.Length!=41943040||metadata.GetProperty("payloadSha256").GetString()!=PayloadSha256||
           Convert.ToHexStringLower(SHA256.HashData(payload))!=PayloadSha256)
            throw new InvalidDataException("原银河密度数据校验失败。");
        coefficients=metadata.GetProperty("coefficients").EnumerateArray().Select(v=>v.GetUInt32()).ToArray();
        if(coefficients.Length!=14||coefficients[11]!=100000||coefficients[12]!=81910||coefficients[13]!=100000)
            throw new InvalidDataException("原银河密度定点系数无效。");
        massDensity=MemoryMarshal.Cast<byte,uint>(payload.AsSpan(0,16777216)).ToArray();
        dustSurface=MemoryMarshal.Cast<byte,ushort>(payload.AsSpan(16777216,8388608)).ToArray();
        milkyWaySurface=MemoryMarshal.Cast<byte,ushort>(payload.AsSpan(25165824,8388608)).ToArray();
        radialHeight=MemoryMarshal.Cast<byte,ushort>(payload.AsSpan(33554432,8388608)).ToArray();
    }

    /// <summary>Exact mass-zero address used by the original MilkyWayLayerJob.</summary>
    public static ulong CellAddress(float absoluteX,float absoluteY,float absoluteZ)
    {
        if(!float.IsFinite(absoluteX+absoluteY+absoluteZ))throw new ArgumentOutOfRangeException(nameof(absoluteX));
        var x=(ulong)(int)(Math.Clamp(absoluteX,0,99999.9f)*.1f);
        var y=(ulong)(int)(Math.Clamp(absoluteY,0,81909.9f)*.1f);
        var z=(ulong)(int)(Math.Clamp(absoluteZ,0,99999.9f)*.1f);
        return (((x<<13)|y)<<14|z)<<3;
    }

    private static void Cell(ulong address,out uint x,out uint y,out uint z)
    {
        if((address&7)!=0)throw new ArgumentOutOfRangeException(nameof(address),"此背景采样路径仅使用原始质量层零。");
        x=(uint)(address>>30)&0x3fff;y=(uint)(address>>17)&0x1fff;z=(uint)(address>>3)&0x3fff;
    }

    // 4CA5370: every intermediate multiplication/truncation is integer.
    public float StellarDensity(ulong address)
    {
        Cell(address,out var x,out _,out var z);
        var qx=(uint)((ulong)x*512*coefficients[1]>>31);
        var qz=(uint)((ulong)z*512*coefficients[4]>>31);
        var planar=Sample(massDensity,qx,qz);
        var height=HeightIntegral(address);
        var value=(((ulong)planar*coefficients[0]>>15)*height)>>15;
        return (float)value*(1f/65536);
    }

    // 4CA5940 -> 4CA5F30: the fixed-point coefficient is deliberately rounded
    // before multiplying the world-space cell centre, as in the original.
    public float MilkyWaySurface(ulong address)
    {
        Cell(address,out var x,out _,out var z);
        var sx=((ulong)Width<<30)/((ulong)coefficients[11]<<15);
        var sz=((ulong)Height<<30)/((ulong)coefficients[13]<<15);
        var qx=(uint)(((ulong)(5+x*10)*512*sx)>>15);
        var qz=(uint)(((ulong)(5+z*10)*512*sz)>>15);
        return Sample(milkyWaySurface,qx,qz)*(1f/32768);
    }

    // 4CA5240: average the four corners in the original dust integral image.
    public float DustDensity(ulong address)
    {
        Cell(address,out var x,out _,out var z);
        var x0=(uint)((ulong)x*4096*coefficients[5]>>34);
        var z0=(uint)((ulong)z*4096*coefficients[6]>>34);
        var x1=(uint)((ulong)(x+1)*coefficients[5]>>22);
        var z1=(uint)((ulong)(z+1)*coefficients[6]>>22);
        var sum=unchecked(Sample(dustSurface,x0,z0)+Sample(dustSurface,x0,z1)+Sample(dustSurface,x1,z0)+Sample(dustSurface,x1,z1));
        return (sum>>2)*(1f/32768);
    }

    // 4CA5A20: integrate the same radial/vertical kernel at all eight cell
    // corners; the map is a pre-integrated source, not a guessed scale height.
    public uint HeightIntegral(ulong address)
    {
        Cell(address,out var x,out var y,out var z);
        var cx=coefficients[11]/20;var cy=coefficients[12]/20;var cz=coefficients[13]/20;
        var ax=Distance(x,cx);var bx=Distance(x+1,cx);var az=Distance(z,cz);var bz=Distance(z+1,cz);
        var ay=Distance(y,cy);var by=Distance(y+1,cy);if(ay>by)(ay,by)=(by,ay);
        var qy0=(uint)((ulong)ay*16777216*coefficients[2]>>39);
        var qy1=(uint)((ulong)by*16777216*coefficients[2]>>39);
        var a=Integral(ax,az,qy0)-Integral(ax,az,qy1);
        var c=Integral(bx,az,qy0)-Integral(bx,az,qy1);
        var e=Integral(ax,bz,qy0)-Integral(ax,bz,qy1);
        var g=Integral(bx,bz,qy0)-Integral(bx,bz,qy1);
        return unchecked(a+c+e+g)>>3;
    }

    private uint Integral(uint x,uint z,uint y)
    {
        var radial=(ulong)(Math.Sqrt((double)x*x+(double)z*z)*32768);
        var qr=(uint)((radial*coefficients[3])>>37);
        return Sample(radialHeight,y,qr);
    }
    private static uint Distance(uint a,uint b)=>a>=b?a-b:b-a;

    private static uint Sample<T>(T[] pixels,uint qx,uint qz) where T:unmanaged
    {
        var x=qx>>9;var z=qz>>9;if(x+1>=Width||z+1>=Height)return 0;
        var wx=(qx&511)*64;var wz=(qz&511)*64;var i=(int)(z*Width+x);
        var top=((ulong)Value(pixels,i)*(32768-wx)>>15)+((ulong)Value(pixels,i+1)*wx>>15);
        var bottom=((ulong)Value(pixels,i+Width)*(32768-wx)>>15)+((ulong)Value(pixels,i+Width+1)*wx>>15);
        return (uint)(((ulong)(uint)top*(32768-wz)>>15)+((ulong)(uint)bottom*wz>>15));
    }
    private static uint Value<T>(T[] pixels,int index) where T:unmanaged=>
        typeof(T)==typeof(uint)?((uint[])(object)pixels)[index]:((ushort[])(object)pixels)[index];

    public uint SourceValue(int field,int x,int z)
    {
        if((uint)x>=Width||(uint)z>=Height)throw new ArgumentOutOfRangeException(nameof(x));
        var i=z*Width+x;
        return field switch {0x38=>massDensity[i],0x78=>dustSurface[i],0x138=>milkyWaySurface[i],0x140=>radialHeight[i],_=>throw new ArgumentOutOfRangeException(nameof(field))};
    }
}

public static class GalaxyOriginalDensity
{
    private static readonly Lazy<GalaxyOriginalDensityMap> reference=new(Read);
    public static GalaxyOriginalDensityMap Reference=>reference.Value;
    private static GalaxyOriginalDensityMap Read()
    {
        var assembly=typeof(GalaxyOriginalDensity).Assembly;
        using var jsonStream=assembly.GetManifestResourceStream("XFEEDToolBox.Core.Rendering.original-galaxy-density.json")!;
        using var json=JsonDocument.Parse(jsonStream);
        using var stream=assembly.GetManifestResourceStream("XFEEDToolBox.Core.Rendering.original-galaxy-density.bin.gz")!;
        using var gzip=new GZipStream(stream,CompressionMode.Decompress);var payload=new byte[41943040];gzip.ReadExactly(payload);
        if(gzip.ReadByte()!=-1)throw new InvalidDataException("原银河密度数据超出范围。");
        return new(json.RootElement,payload);
    }
}