VoxelForge/Standalone/Source/StandaloneResearchTunnels.h

#pragma once

#include <atomic>

// 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 std::atomic<bool> Enabled{false};
    inline std::atomic<float> CenterZVoxels{0.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 void Set(bool bEnabled, float CenterZ)
    {
        CenterZVoxels.store(CenterZ,std::memory_order_relaxed);
        Enabled.store(bEnabled,std::memory_order_release);
    }

    inline float Apply(float BaseDensity, float X, float Y, float Z, uint32 Seed)
    {
        if (!Enabled.load(std::memory_order_acquire)) return BaseDensity;
        const float CenterZ = CenterZVoxels.load(std::memory_order_relaxed);
        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*FMath::Sin(6.2831853071795864769f*T+Phase));
            const float VerticalWave=4.0f*FMath::Sin(3.14159265358979323846f*T)
                *FMath::Sin(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)));
            Best=FMath::Max(Best,(1.0f-Normalized)*FMath::Min(Radius,HalfHeight));
        };

        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);
        }
        return bHasEdge ? FMath::Max(BaseDensity,Best) : BaseDensity;
    }
}