#pragma once

// Standalone-only port of the research field's round, winding tunnel profile.
// The cubic centerline, variable radius and elliptical vertical section follow
// GeneratedTunnelEdgeFromColumnProfile in VoxelRedefinedDensityField.cpp.
namespace FStandaloneResearchTunnels
{
    inline constexpr float WorldCenterZVoxels=-80.0f;

    inline uint32 Mix(uint32 X)
    {
        X ^= X >> 16;
        X *= 0x7FEB352Du;
        X ^= X >> 15;
        X *= 0x846CA68Bu;
        X ^= X >> 16;
        return X;
    }

    inline uint32 EdgeHash(int32 X, int32 Y, int32 Axis, uint32 Seed)
    {
        uint32 H = Seed ^ 0x9E3779B9u;
        H = Mix(H ^ static_cast<uint32>(X) * 0x85EBCA6Bu);
        H = Mix(H ^ static_cast<uint32>(Y) * 0xC2B2AE35u);
        return Mix(H ^ static_cast<uint32>(Axis) * 0x27D4EB2Fu);
    }

    inline float Unit(uint32 X)
    {
        return static_cast<float>(Mix(X) >> 8) * (1.0f / 16777216.0f);
    }

    inline uint32 NoiseHash(int32 X,int32 Y,int32 Z,uint32 Seed)
    {
        uint32 H=Mix(Seed^0xB5297A4Du);
        H=Mix(H^static_cast<uint32>(X)*0x68E31DA4u);
        H=Mix(H^static_cast<uint32>(Y)*0x1B56C4E9u);
        return Mix(H^static_cast<uint32>(Z)*0x9E3779B9u);
    }

    inline float NoiseValue(int32 X,int32 Y,int32 Z,uint32 Seed)
    {
        return static_cast<float>(NoiseHash(X,Y,Z,Seed)>>8)*(2.0f/16777215.0f)-1.0f;
    }

    inline float Fade(float T)
    {
        return T*T*(3.0f-2.0f*T);
    }

    inline float ValueNoise(float X,float Y,float Z,int32 Wavelength,uint32 Seed)
    {
        const float Inv=1.0f/static_cast<float>(Wavelength);
        const float FX=X*Inv,FY=Y*Inv,FZ=Z*Inv;
        const int32 X0=FMath::FloorToInt(FX),Y0=FMath::FloorToInt(FY),Z0=FMath::FloorToInt(FZ);
        const float TX=Fade(FX-static_cast<float>(X0));
        const float TY=Fade(FY-static_cast<float>(Y0));
        const float TZ=Fade(FZ-static_cast<float>(Z0));
        const float Z00=FMath::Lerp(NoiseValue(X0,Y0,Z0,Seed),NoiseValue(X0+1,Y0,Z0,Seed),TX);
        const float Z10=FMath::Lerp(NoiseValue(X0,Y0+1,Z0,Seed),NoiseValue(X0+1,Y0+1,Z0,Seed),TX);
        const float Z01=FMath::Lerp(NoiseValue(X0,Y0,Z0+1,Seed),NoiseValue(X0+1,Y0,Z0+1,Seed),TX);
        const float Z11=FMath::Lerp(NoiseValue(X0,Y0+1,Z0+1,Seed),NoiseValue(X0+1,Y0+1,Z0+1,Seed),TX);
        return FMath::Lerp(FMath::Lerp(Z00,Z10,TY),FMath::Lerp(Z01,Z11,TY),TZ);
    }

    inline float WallRelief(float X,float Y,float Z,uint32 Seed)
    {
        return 0.70f*ValueNoise(X,Y,Z,8,Seed^0xA511E9B3u)
            +1.05f*ValueNoise(X,Y,Z,16,Seed^0x63D83595u)
            +1.45f*ValueNoise(X,Y,Z,32,Seed^0xB5297A4Du);
    }

    inline float SmoothMax(float A,float B,float Width)
    {
        const float H=FMath::Max(Width-FMath::Abs(A-B),0.0f)/Width;
        return FMath::Max(A,B)+(H*H)*(Width*0.25f);
    }

    inline float Apply(float BaseDensity, float X, float Y, float Z, uint32 Seed,
                       bool bEnabled, bool bEnhancedWalls)
    {
        if (!bEnabled) return BaseDensity;
        constexpr float CenterZ = WorldCenterZVoxels;
        if (FMath::Abs(Z-CenterZ) > 36.0f) return BaseDensity;

        constexpr float Pitch = 144.0f;
        constexpr float RadiusBase = 28.0f;
        constexpr float HalfHeight = 26.0f;
        const int32 BaseCellX = FMath::FloorToInt(X/Pitch);
        const int32 BaseCellY = FMath::FloorToInt(Y/Pitch);
        float Best = -FLT_MAX;
        bool bHasEdge = false;
        const auto CellCenter = [&](int32 CX, int32 CY, float& OutX, float& OutY)
        {
            const uint32 H = EdgeHash(CX,CY,3,Seed);
            OutX = (static_cast<float>(CX)+0.5f)*Pitch
                + (Unit(H^0xA511E9B3u)-0.5f)*Pitch*0.18f;
            OutY = (static_cast<float>(CY)+0.5f)*Pitch
                + (Unit(H^0x63D83595u)-0.5f)*Pitch*0.18f;
        };
        const auto EvaluateEdge = [&](int32 CX, int32 CY, int32 Axis)
        {
            const uint32 H = EdgeHash(CX,CY,Axis,Seed);
            if (Unit(H^0x68E31DA4u) > 0.56f) return;
            bHasEdge = true;
            const int32 NX = CX + (Axis == 0 ? 1 : 0);
            const int32 NY = CY + (Axis == 1 ? 1 : 0);
            float AX,AY,BX,BY;
            CellCenter(CX,CY,AX,AY);
            CellCenter(NX,NY,BX,BY);

            const float ABX = BX-AX, ABY = BY-AY;
            const float LengthSquared = ABX*ABX+ABY*ABY;
            float T = LengthSquared > 1.0e-8f
                ? FMath::Clamp(((X-AX)*ABX+(Y-AY)*ABY)/LengthSquared,0.0f,1.0f)
                : 0.0f;
            const float Length = FMath::Sqrt(LengthSquared);
            const float NXrm = Length > 1.0e-4f ? -ABY/Length : 0.0f;
            const float NYrm = Length > 1.0e-4f ? ABX/Length : 0.0f;
            const float BendA = (Unit(H^0xB5297A4Du)*2.0f-1.0f)*Pitch*0.16f;
            const float BendB = (Unit(H^0x1B56C4E9u)*2.0f-1.0f)*Pitch*0.16f;
            const float C1X=AX+ABX/3.0f+NXrm*BendA, C1Y=AY+ABY/3.0f+NYrm*BendA;
            const float C2X=AX+ABX*(2.0f/3.0f)+NXrm*BendB;
            const float C2Y=AY+ABY*(2.0f/3.0f)+NYrm*BendB;

            for (int32 Iteration=0;Iteration<4;++Iteration)
            {
                const float U=T, O=1.0f-U;
                const float PX=O*O*O*AX+3.0f*O*O*U*C1X+3.0f*O*U*U*C2X+U*U*U*BX;
                const float PY=O*O*O*AY+3.0f*O*O*U*C1Y+3.0f*O*U*U*C2Y+U*U*U*BY;
                const float DX=3.0f*O*O*(C1X-AX)+6.0f*O*U*(C2X-C1X)+3.0f*U*U*(BX-C2X);
                const float DY=3.0f*O*O*(C1Y-AY)+6.0f*O*U*(C2Y-C1Y)+3.0f*U*U*(BY-C2Y);
                const float DDX=6.0f*O*(C2X-2.0f*C1X+AX)+6.0f*U*(BX-2.0f*C2X+C1X);
                const float DDY=6.0f*O*(C2Y-2.0f*C1Y+AY)+6.0f*U*(BY-2.0f*C2Y+C1Y);
                const float Den=DX*DX+DY*DY+(PX-X)*DDX+(PY-Y)*DDY;
                if (FMath::Abs(Den)>1.0e-4f)
                    T=FMath::Clamp(U-((PX-X)*DX+(PY-Y)*DY)/Den,0.0f,1.0f);
            }

            const float O=1.0f-T;
            const float PX=O*O*O*AX+3.0f*O*O*T*C1X+3.0f*O*T*T*C2X+T*T*T*BX;
            const float PY=O*O*O*AY+3.0f*O*O*T*C1Y+3.0f*O*T*T*C2Y+T*T*T*BY;
            const float DX=X-PX,DY=Y-PY;
            const float Distance=FMath::Sqrt(DX*DX+DY*DY);
            const float Phase=Unit(H^0xA511E9B3u)*6.2831853071795864769f;
            const float Radius=RadiusBase*(0.82f+0.36f*Unit(H^0x63D83595u))
                *(1.0f+0.22f*VoxelMath::DetSin(6.2831853071795864769f*T+Phase));
            const float VerticalWave=4.0f*VoxelMath::DetSin(3.14159265358979323846f*T)
                *VoxelMath::DetSin(6.2831853071795864769f*T+Phase);
            const float VerticalOffset=Z-(CenterZ+VerticalWave);
            const float HorizontalDistance=FMath::Max(0.0f,Distance);
            const float Normalized=FMath::Sqrt(
                FMath::Square(HorizontalDistance/FMath::Max(Radius,1.0f))
                +FMath::Square(VerticalOffset/FMath::Max(HalfHeight,1.0f)));
            const float EdgeDensity=(1.0f-Normalized)*FMath::Min(Radius,HalfHeight);
            Best=bEnhancedWalls&&Best>-FLT_MAX
                ?SmoothMax(Best,EdgeDensity,4.0f):FMath::Max(Best,EdgeDensity);
        };

        for (int32 CY=BaseCellY-1;CY<=BaseCellY+1;++CY)
        for (int32 CX=BaseCellX-1;CX<=BaseCellX+1;++CX)
        {
            EvaluateEdge(CX,CY,0);
            EvaluateEdge(CX,CY,1);
        }
        if(!bHasEdge)return BaseDensity;
        if(bEnhancedWalls)Best+=WallRelief(X,Y,Z,Seed);
        // The wider base/round-field blend removes tiny marching-cubes wedges where a round
        // branch meets an existing room or passage wall. Relief remains active on the tunnel
        // wall; this only rounds the structural union boundary.
        return bEnhancedWalls?SmoothMax(BaseDensity,Best,8.0f):FMath::Max(BaseDensity,Best);
    }
}