using System;
using System.Collections.Generic;
using Sandbox;

namespace TerrainLines;

// A distance-based height profile, independent of point density and the transverse shoulder.
internal sealed class LineProfile
{
	private readonly double[] _heights;
	private readonly float _length;
	private readonly bool[] _roundJoins;

	internal LineProfile( TerrainStorage storage, IReadOnlyList<Vector3> points, float length, float evenRamp = 0, float slopeLimit = 90 )
	{
		_length = length;
		int intervals = (int)Math.Clamp( Math.Ceiling( length / (storage.TerrainSize / storage.Resolution) ), 8, 65536 );
		_heights = new double[intervals + 1];
		float travelled = 0;
		int segment = 1;
		float segmentLength = LineRaster.Distance2D( points[0], points[1] );
		for ( int i = 0; i <= intervals; i++ )
		{
			float distance = length * i / intervals;
			while ( segment < points.Count - 1 && travelled + segmentLength < distance )
			{
				travelled += segmentLength;
				segment++;
				segmentLength = LineRaster.Distance2D( points[segment - 1], points[segment] );
			}
			float t = segmentLength > 0 ? Math.Clamp( (distance - travelled) / segmentLength, 0, 1 ) : 0;
			var position = points[segment - 1] + (points[segment] - points[segment - 1]) * t;
			_heights[i] = LineRaster.SampleHeight( storage, position.x, position.y );
		}

		// Four implicit diffusion passes spread abrupt rises over roughly a tenth of the stroke.
		// Fixed endpoints preserve arrival heights at this stage; linear grades are unchanged.
		double lambda = Math.Pow( intervals * 0.1, 2 ) / 8;
		var upper = new double[intervals];
		var rhs = new double[intervals];
		for ( int i = 1; i < intervals; i++ )
			upper[i] = -lambda / (1 + 2 * lambda + lambda * upper[i - 1]);
		for ( int pass = 0; pass < 4; pass++ )
		{
			for ( int i = 1; i < intervals; i++ )
			{
				double value = _heights[i];
				if ( i == 1 ) value += lambda * _heights[0];
				if ( i == intervals - 1 ) value += lambda * _heights[intervals];
				rhs[i] = (value + lambda * rhs[i - 1]) / (1 + 2 * lambda + lambda * upper[i - 1]);
			}
			for ( int i = intervals - 1; i > 0; i-- )
				_heights[i] = rhs[i] - (i < intervals - 1 ? upper[i] * _heights[i + 1] : 0);
		}
		double start = _heights[0], end = _heights[^1];
		for ( int i = 0; i <= intervals; i++ )
		{
			double grade = start + (end - start) * i / intervals;
			_heights[i] += (grade - _heights[i]) * Math.Clamp( evenRamp, 0, 1 );
		}
		if ( slopeLimit < 90 )
		{
			var original = (double[])_heights.Clone();
			double rise = Math.Tan( Math.Clamp( slopeLimit, 0, 90 ) * Math.PI / 180 ) * length / intervals;
			// The smallest raised envelope with this grade bound. Both sweeps are needed for hills and valleys.
			for ( int i = 1; i <= intervals; i++ ) _heights[i] = Math.Max( _heights[i], _heights[i - 1] - rise );
			for ( int i = intervals - 1; i >= 0; i-- ) _heights[i] = Math.Max( _heights[i], _heights[i + 1] - rise );
			_roundJoins = new bool[intervals + 1];
			for ( int i = 1; i < intervals; i++ )
				_roundJoins[i] = _heights[i - 1] > original[i - 1] + 1e-8 || _heights[i] > original[i] + 1e-8 || _heights[i + 1] > original[i + 1] + 1e-8;
		}
	}

	internal float Sample( float distance )
	{
		float t = Math.Clamp( distance / _length, 0, 1 );
		float sample = t * (_heights.Length - 1);
		int knot = (int)MathF.Round( sample );
		if ( knot > 0 && knot < _heights.Length - 1 && (_roundJoins?[knot] ?? false) )
		{
			// Round only affected joins over one sample interval. The derivative is a convex blend
			// of the two compliant grades, so rounding cannot introduce a steeper section.
			double u = sample - knot + 0.5;
			double a = (_heights[knot - 1] + _heights[knot]) * 0.5;
			double b = _heights[knot];
			double c = (_heights[knot] + _heights[knot + 1]) * 0.5;
			return (float)((1 - u) * (1 - u) * a + 2 * (1 - u) * u * b + u * u * c);
		}
		int index = Math.Min( (int)sample, _heights.Length - 2 );
		return (float)(_heights[index] + (_heights[index + 1] - _heights[index]) * (sample - index));
	}
}