Ocean.cs
using System;
using System.Collections.Generic;
using Sandbox.Rendering;

namespace OceanLib;

/// <summary>
/// Simulates the open ocean on the GPU, after Tessendorf's FFT method. One wind spectrum (Phillips, see
/// <see cref="OceanSpectrum"/>) is split by wave length across three cascades of different tile sizes - the swell, the
/// waves on it and the chop on those - so the sea has detail from hundreds of metres down to a hand's width, and the
/// tiles, being of unrelated sizes, never visibly repeat. Each cascade's spectrum is moved to the current time and
/// turned into a surface by an inverse FFT in compute shaders, giving two maps per cascade:
/// <see cref="Cascade.Displacement"/> (choppy x and y shift, height, shear) and <see cref="Cascade.Derivatives"/>
/// (slopes and stretch, mip-mapped for distant shading). The sea depends only on these settings and the time, so every
/// machine sees the same waves.
///
/// The sea state - wind, wave height, look - comes from an <see cref="OceanPreset"/>. The ocean can hold two at once
/// and blend between them (<see cref="BlendTo"/>, <see cref="SetBlend"/>): both are built once and kept, and the blend
/// is a number the GPU and the CPU replica each apply as they go, so a calm can build into a storm without anything
/// being rebuilt.
///
/// For physics, <see cref="HeightAt"/> and <see cref="WaterVelocityAt"/> read a CPU replica of the swell
/// (<see cref="OceanWaves"/>). The GPU work runs from a command list on the main camera at the end of each frame, for
/// the time of the next one, timed as "Ocean FFT" in the GPU profiler.
/// </summary>
[Title( "Ocean" ), Category( "Ocean" ), Icon( "waves" )]
public sealed class Ocean : Component
{
	/// <summary>Waves per side of a cascade's tile. Fixed by the FFT shader's shared memory.</summary>
	public const int Size = 256;

	public const int CascadeCount = 3;

	/// <summary>s&amp;box units (inches) per metre. Every length the ocean exposes is in metres unless it says otherwise.</summary>
	public const float UnitsPerMetre = 39.37f;

	/// <summary>The sea state to start from. Empty is a moderate sea (<see cref="OceanSea.Default"/>).</summary>
	[Property, Group( "Sea" )] public OceanPreset Preset { get; set; }

	/// <summary>Which way the wind blows, in degrees from +x.</summary>
	[Property, Group( "Sea" ), Range( 0, 360 )] public float WindAngle { get; set; } = 0.0f;

	/// <summary>How tightly the waves follow the wind: 2 is Tessendorf's, higher lines them up more.</summary>
	[Property, Group( "Sea" ), Range( 1, 8 )] public float Directionality { get; set; } = 2.0f;

	/// <summary>Waves much shorter than this, in metres, are damped out.</summary>
	[Property, Group( "Sea" ), Range( 0, 2 )] public float SmallWaveCutoff { get; set; } = 0.1f;

	[Property, Group( "Sea" )] public int Seed { get; set; } = 1;

	/// <summary>Metres covered by each cascade's tile, longest first. Unrelated sizes keep the tiles from lining up.</summary>
	[Property, Group( "Cascades" )] public Vector3 CascadeLengths { get; set; } = new( 500.0f, 97.0f, 19.3f );

	/// <summary>
	/// Where one cascade hands over to the next, in waves per tile of the finer one: each cascade keeps the waves too
	/// long for the next to hold this many of.
	/// </summary>
	[Property, Group( "Cascades" ), Range( 2, 16 )] public float BandSplit { get; set; } = 12.0f;

	/// <summary>Seconds after which the sea repeats exactly - long enough that nobody notices.</summary>
	[Property, Group( "Time" )] public float LoopPeriod { get; set; } = 200.0f;

	[Property, Group( "Time" ), Range( 0, 4 )] public float TimeScale { get; set; } = 1.0f;

	/// <summary>One tile size of the sea and the maps the FFT writes for it.</summary>
	public sealed class Cascade
	{
		/// <summary>Metres covered by the tile.</summary>
		public float TileLength { get; init; }

		/// <summary>Choppy x and y shift and height, in metres, and the shear dDx/dy in w.</summary>
		public Texture Displacement { get; init; }

		/// <summary>Slope of the height in x and y, and the stretch of the shift along x and y. Mip-mapped.</summary>
		public Texture Derivatives { get; init; }

		internal Texture SpectrumA { get; init; }
		internal Texture SpectrumB { get; init; }
	}

	/// <summary>One sea state, built: each cascade's initial spectrum, on the CPU and on the GPU.</summary>
	public sealed class SpectrumSet
	{
		public OceanSpectrum[] Spectra { get; init; }
		public OceanSea Sea { get; init; }
		internal Texture[] Textures { get; init; }

		internal static SpectrumSet Create( OceanSpectrum[] spectra, OceanSea sea )
		{
			var textures = new Texture[spectra.Length];

			for ( int i = 0; i < spectra.Length; i++ )
			{
				textures[i] = Texture.Create( Size, Size, ImageFormat.RGBA32323232F )
					.WithUAVBinding()
					.WithData<Vector4>( spectra[i].Texels.AsSpan() )
					.WithName( $"Ocean Initial Spectrum {i}" )
					.Finish();
			}

			return new SpectrumSet { Spectra = spectra, Sea = sea, Textures = textures };
		}

		/// <summary>The sea part way between two states, exactly as the blend shows it.</summary>
		internal static SpectrumSet Lerp( SpectrumSet a, SpectrumSet b, float t )
		{
			var spectra = new OceanSpectrum[a.Spectra.Length];

			for ( int i = 0; i < spectra.Length; i++ )
				spectra[i] = OceanSpectrum.Lerp( a.Spectra[i], b.Spectra[i], t );

			return Create( spectra, OceanSea.Lerp( a.Sea, b.Sea, t ) );
		}
	}

	/// <summary>The cascades, longest waves first.</summary>
	public Cascade[] Cascades { get; private set; } = Array.Empty<Cascade>();

	/// <summary>The two sea states being blended, and how far from the first to the second.</summary>
	public SpectrumSet From { get; private set; }
	public SpectrumSet To { get; private set; }
	public float Blend { get; private set; }

	/// <summary>The sea as it is now, between the two states being blended.</summary>
	public OceanSea Sea => From is null ? StartSea : OceanSea.Lerp( From.Sea, To.Sea, Blend );

	/// <summary>The time, wrapped to the loop period, the maps currently show.</summary>
	public float OceanTime { get; private set; }

	/// <summary>The sea replayed on the CPU, for physics: see <see cref="OceanWaves"/>.</summary>
	public OceanWaves Waves { get; private set; }

	/// <summary>The first active ocean in the scene: what floating things float on.</summary>
	public static Ocean Current { get; private set; }

	/// <summary>
	/// The sea's time at this very moment, wrapped to the loop period: unlike <see cref="OceanTime"/>, which steps
	/// once a frame, it moves on between the physics steps of a frame too.
	/// </summary>
	public float PhysicsTime => (float)((m_Time + (Time.Now - m_TimeTakenAt) * TimeScale) % LoopPeriod + LoopPeriod) % LoopPeriod;

	/// <summary>Which way the wind blows, as a unit vector across the sea.</summary>
	public Vector3 WindDirection => Rotation.FromYaw( WindAngle ).Forward;

	/// <summary>The height of the sea's surface over a point in the world, in world units. The mean level is z = 0.</summary>
	public float HeightAt( Vector3 world )
	{
		if ( Waves is null )
			return 0.0f;

		return Waves.HeightAt( ToMetres( world ), PhysicsTime, Sea.Choppiness ) * UnitsPerMetre;
	}

	/// <summary>
	/// Where the water that lies over a point at rest has been carried by the waves, in world units: sideways along
	/// the orbit and up or down.
	/// </summary>
	public Vector3 DisplacementAt( Vector3 world )
	{
		if ( Waves is null )
			return Vector3.Zero;

		return Waves.Displacement( ToMetres( world ), PhysicsTime, Sea.Choppiness ) * UnitsPerMetre;
	}

	/// <summary>
	/// How fast the water at the surface near a point is moving, in world units per second: the orbit each bit of
	/// water swings round as a wave passes - forward on the crest, back in the trough - which is what carries a
	/// floating thing along. From the swell's displacement a moment apart.
	/// </summary>
	public Vector3 WaterVelocityAt( Vector3 world )
	{
		if ( Waves is null )
			return Vector3.Zero;

		// Without the choppiness: that pinches the crests for the eye, and would carry the water round its orbit nearly
		// twice as fast as a real sea does.
		const float dt = 0.05f;
		var metres = ToMetres( world );
		float t = PhysicsTime;
		var now = Waves.Displacement( metres, t, 1.0f );
		var next = Waves.Displacement( metres, t + dt, 1.0f );
		return (next - now) / dt * UnitsPerMetre;
	}

	/// <summary>
	/// Moves the sea to another state over some seconds, from wherever it is now - even part way through another
	/// blend. Zero seconds switches at once.
	/// </summary>
	public void BlendTo( OceanPreset preset, float seconds )
	{
		if ( !m_Driven )
			m_From = StartSea;

		// Part way through a blend: what shows now becomes the start of the new one, so nothing jumps.
		if ( From is not null && Blend > 0.0f )
		{
			var snapshot = SpectrumSet.Lerp( From, To, Blend );
			DropSnapshot();
			m_Snapshot = snapshot;
			m_From = snapshot.Sea;
		}

		m_To = preset?.Sea ?? OceanSea.Default;
		m_Blend = 0.0f;
		m_Rate = seconds > 0.0f ? 1.0f / seconds : float.PositiveInfinity;
		m_Driven = true;
	}

	/// <summary>
	/// Holds the sea at a fixed point between two states: 0 is all of the first, 1 all of the second. For a slider or
	/// anything else that drives the blend itself every frame.
	/// </summary>
	public void SetBlend( OceanPreset from, OceanPreset to, float blend )
	{
		DropSnapshot();
		m_From = from?.Sea ?? OceanSea.Default;
		m_To = to?.Sea ?? OceanSea.Default;
		m_Blend = blend.Clamp( 0.0f, 1.0f );
		m_Rate = 0.0f;
		m_Driven = true;
	}

	/// <summary>The tile length of a cascade, in metres.</summary>
	public float LengthOf( int cascade ) => cascade switch { 0 => CascadeLengths.x, 1 => CascadeLengths.y, _ => CascadeLengths.z };

	/// <summary>Hands the maps and their tile sizes to anything drawn or dispatched with the ocean's shaders.</summary>
	public void Bind( RenderAttributes attributes )
	{
		for ( int i = 0; i < Cascades.Length; i++ )
		{
			attributes.Set( $"OceanDisplacement{i}", Cascades[i].Displacement );
			attributes.Set( $"OceanDerivatives{i}", Cascades[i].Derivatives );
		}

		attributes.Set( "OceanTileLengths", CascadeLengths * UnitsPerMetre );
		attributes.Set( "OceanUnitsPerMetre", UnitsPerMetre );
	}

	/// <summary>The same, for a command list.</summary>
	public void Bind( CommandList commands )
	{
		for ( int i = 0; i < Cascades.Length; i++ )
		{
			commands.Attributes.Set( $"OceanDisplacement{i}", Cascades[i].Displacement );
			commands.Attributes.Set( $"OceanDerivatives{i}", Cascades[i].Derivatives );
		}

		commands.Attributes.Set( "OceanTileLengths", CascadeLengths * UnitsPerMetre );
		commands.Attributes.Set( "OceanUnitsPerMetre", UnitsPerMetre );
	}

	private ComputeShader m_TimeShader;
	private ComputeShader m_FftShader;
	private ComputeShader m_AssembleShader;
	private CameraCommands m_Commands;
	private double m_Time;
	private double m_TimeTakenAt;

	// Every sea state built so far, and the CPU replica for every pair blended so far.
	private readonly Dictionary<OceanSea, SpectrumSet> m_Sets = new();
	private readonly Dictionary<(SpectrumSet, SpectrumSet), OceanWaves> m_Replicas = new();
	private int m_LayoutHash;

	// What BlendTo or SetBlend asked for; until either is called, the preset is the one sea state.
	private bool m_Driven;
	private OceanSea m_From, m_To;
	private float m_Blend, m_Rate;

	// The sea as it stood when a blend was interrupted: built from the two it was between, not from settings.
	private SpectrumSet m_Snapshot;

	private static readonly ProfilingSampler s_Profile = new( "Ocean FFT" );

	private OceanSea StartSea => Preset?.Sea ?? OceanSea.Default;

	private static Vector2 ToMetres( Vector3 world ) => new Vector2( world.x, world.y ) / UnitsPerMetre;

	private void DropSnapshot()
	{
		if ( m_Snapshot is null )
			return;

		foreach ( var key in m_Replicas.Keys.Where( k => k.Item1 == m_Snapshot ).ToList() )
			m_Replicas.Remove( key );

		m_Snapshot = null;
	}

	protected override void OnEnabled()
	{
		m_TimeShader = new ComputeShader( "ocean/shaders/ocean_timespectrum_cs.shader" );
		m_FftShader = new ComputeShader( "ocean/shaders/ocean_fft_cs.shader" );
		m_AssembleShader = new ComputeShader( "ocean/shaders/ocean_assemble_cs.shader" );

		m_Commands = new CameraCommands( "Ocean FFT", Stage.AfterUI, 0 );
		m_LayoutHash = 0;
		Refresh();
		Current ??= this;
	}

	protected override void OnDisabled()
	{
		if ( Current == this )
			Current = null;

		m_Commands?.Dispose();
		m_Commands = null;
	}

	private static Texture CreateTarget( string name, ImageFormat format, bool mips = false )
	{
		var builder = Texture.Create( Size, Size, format ).WithUAVBinding().WithName( name );
		return (mips ? builder.WithMips() : builder).Finish();
	}

	/// <summary>Brings the cascades, the two sea states and the CPU replica up to date with the settings.</summary>
	private void Refresh()
	{
		// The tiles and the spectrum's shape set every sea state: when they change, everything built so far is stale.
		int layout = HashCode.Combine( CascadeLengths, BandSplit, Seed, LoopPeriod, WindAngle, Directionality, SmallWaveCutoff );

		if ( layout != m_LayoutHash || Cascades.Length != CascadeCount )
		{
			m_LayoutHash = layout;
			m_Sets.Clear();
			m_Replicas.Clear();

			// An interrupted blend's snapshot was built from the old shape: start again from where it stood.
			if ( m_Snapshot is not null )
			{
				m_From = m_Snapshot.Sea;
				m_Snapshot = null;
			}

			var cascades = new Cascade[CascadeCount];
			for ( int i = 0; i < CascadeCount; i++ )
			{
				var old = i < Cascades.Length ? Cascades[i] : null;
				cascades[i] = new Cascade
				{
					TileLength = LengthOf( i ),
					SpectrumA = old?.SpectrumA ?? CreateTarget( $"Ocean Spectrum A {i}", ImageFormat.RGBA32323232F ),
					SpectrumB = old?.SpectrumB ?? CreateTarget( $"Ocean Spectrum B {i}", ImageFormat.RGBA32323232F ),
					// Half floats: a few millimetres of precision on waves of metres, half the bandwidth to read them.
					Displacement = old?.Displacement ?? CreateTarget( $"Ocean Displacement {i}", ImageFormat.RGBA16161616F ),
					Derivatives = old?.Derivatives ?? CreateTarget( $"Ocean Derivatives {i}", ImageFormat.RGBA16161616F, mips: true ),
				};
			}

			Cascades = cascades;
		}

		if ( !m_Driven )
			m_From = m_To = StartSea;

		From = m_Snapshot ?? SetFor( m_From );
		To = SetFor( m_To );
		Blend = From == To ? 0.0f : m_Blend;

		if ( !m_Replicas.TryGetValue( (From, To), out var replica ) )
			m_Replicas[(From, To)] = replica = new OceanWaves( From.Spectra[0], To.Spectra[0], LoopPeriod );

		Waves = replica;
		Waves.Blend = Blend;
	}

	/// <summary>The built sea state for some settings, building it the first time.</summary>
	private SpectrumSet SetFor( OceanSea sea )
	{
		if ( m_Sets.TryGetValue( sea, out var set ) )
			return set;

		var settings = new OceanSpectrum.Settings( sea.WindSpeed, WindAngle, Directionality, SmallWaveCutoff, Seed );
		var spectra = new OceanSpectrum[CascadeCount];
		double variance = 0.0;

		for ( int i = 0; i < CascadeCount; i++ )
		{
			// Each cascade takes the waves too long for the next one to hold BandSplit of; the last takes the rest.
			float minK = i == 0 ? 0.0f : 2.0f * MathF.PI / LengthOf( i ) * BandSplit;
			float maxK = i == CascadeCount - 1 ? float.MaxValue : 2.0f * MathF.PI / LengthOf( i + 1 ) * BandSplit;

			spectra[i] = new OceanSpectrum( Size, LengthOf( i ), minK, maxK, i, settings );
			variance += spectra[i].RawVariance;
		}

		float scale = OceanSpectrum.ScaleFor( sea.SignificantHeight, variance );

		for ( int i = 0; i < CascadeCount; i++ )
			spectra[i].Scale( scale );

		set = SpectrumSet.Create( spectra, sea );
		m_Sets[sea] = set;
		return set;
	}

	protected override void OnUpdate()
	{
		if ( m_Commands is null )
			return;

		AdvanceBlend();
		Refresh();
		m_Commands.Follow( Scene );

		// The maps are built at the end of this frame, for the moment the next one will show.
		m_Time += Time.Delta * TimeScale;
		m_TimeTakenAt = Time.Now;
		OceanTime = (float)((m_Time + Time.Delta * TimeScale) % LoopPeriod);
		Current ??= this;

		// The CPU replica's next grid gets started on a worker thread now, so it's ready before the physics needs it.
		Waves.Prefetch( PhysicsTime );

		Record();
	}

	private void AdvanceBlend()
	{
		if ( m_Rate <= 0.0f )
			return;

		m_Blend = float.IsInfinity( m_Rate ) ? 1.0f : m_Blend.Approach( 1.0f, Time.Delta * m_Rate );

		if ( m_Blend < 1.0f )
			return;

		// Arrived: the target is now the one sea state.
		m_From = m_To;
		DropSnapshot();
		m_Blend = 0.0f;
		m_Rate = 0.0f;
	}

	private void Record()
	{
		var cl = m_Commands.List;
		cl.Reset();

		using ( cl.ProfileScope( s_Profile ) )
		{
			cl.Attributes.Set( "OceanSize", Size );
			cl.Attributes.Set( "OceanTime", OceanTime );
			cl.Attributes.Set( "LoopPeriod", LoopPeriod );
			cl.Attributes.Set( "Choppiness", Sea.Choppiness );
			cl.Attributes.Set( "SpectrumBlend", Blend );

			for ( int i = 0; i < Cascades.Length; i++ )
			{
				var cascade = Cascades[i];

				// 1. The spectrum at this time, blended between the two sea states.
				cl.Attributes.Set( "InitialSpectrum", From.Textures[i] );
				cl.Attributes.Set( "InitialSpectrumB", To.Textures[i] );
				cl.Attributes.Set( "SpectrumA", cascade.SpectrumA );
				cl.Attributes.Set( "SpectrumB", cascade.SpectrumB );
				cl.Attributes.Set( "TileLength", cascade.TileLength );
				cl.DispatchCompute( m_TimeShader, Size, Size, 1 );
				cl.UavBarrier( cascade.SpectrumA );
				cl.UavBarrier( cascade.SpectrumB );

				// 2. Inverse FFT along x, then along y: one thread group per line.
				for ( int direction = 0; direction < 2; direction++ )
				{
					cl.Attributes.Set( "Direction", direction );
					cl.DispatchCompute( m_FftShader, Size, Size, 1 );
					cl.UavBarrier( cascade.SpectrumA );
					cl.UavBarrier( cascade.SpectrumB );
				}

				// 3. The maps the surface reads, the slopes mip-mapped so distant waves shade without shimmering.
				cl.Attributes.Set( "Displacement", cascade.Displacement );
				cl.Attributes.Set( "Derivatives", cascade.Derivatives );
				cl.DispatchCompute( m_AssembleShader, Size, Size, 1 );
				cl.UavBarrier( cascade.Displacement );
				cl.UavBarrier( cascade.Derivatives );
				cl.GenerateMipMaps( cascade.Derivatives, Graphics.DownsampleMethod.Box );
			}
		}
	}
}