Editor/Services/ArchFaceOverlap.cs
using System;
using System.Collections.Generic;
using System.Linq;
using HalfEdgeMesh;
using Sandbox;

namespace Sunless.Architecture;

// One built face, with the plane it lies in already keyed - the key is how a face finds the only faces it could
// ever be covered by, and it is what an incremental build uses to decide which parts still have to be compared.
sealed class ArchOverlapFace {
	public PolygonMesh Mesh { get; init; }
	public FaceHandle Handle { get; init; }
	public Vector3[] Corners { get; init; }
	public Vector3 Normal { get; init; }
	public Vector3 Centre { get; init; }
	public float Area { get; init; }
	public int Part { get; init; }
	public (long, long, long, long) Plane { get; init; }
}

public sealed class ArchFaceOverlap {
	const float PlaneAngle = 0.02f;
	const float PlaneOffset = 0.06f;

	public int RemoveCoveredFaces( IReadOnlyList<ArchBuiltPart> parts ) {
		return Resolve( Collect( parts, null ) );
	}

	// Faces in EMISSION order, and a subset is a filter rather than a re-ordering: the pass breaks a mutual cover
	// by area but breaks everything else by which face it reached first, so a reordered subset could kill the other
	// half of a pair and a warm build would stop matching a cold one.
	internal List<ArchOverlapFace> Collect( IReadOnlyList<ArchBuiltPart> parts, IReadOnlySet<int> only ) {
		var faces = new List<ArchOverlapFace>();

		for ( var partIndex = 0; partIndex < parts.Count; partIndex++ ) {
			if ( only is not null && !only.Contains( partIndex ) ) {
				continue;
			}

			faces.AddRange( Faces( parts[partIndex], partIndex ) );
		}

		return faces;
	}

	internal static IEnumerable<ArchOverlapFace> Faces( ArchBuiltPart part, int partIndex ) {
		var mesh = part.Canvas.Finish();

		foreach ( var handle in mesh.FaceHandles.ToList() ) {
			var corners = mesh.GetFaceVertexPositions( handle, part.Canvas.Projection ).ToArray();

			if ( corners.Length < 3 ) {
				continue;
			}

			mesh.ComputeFaceNormal( handle, out var local );

			var normal = part.Canvas.Projection.NormalToWorld( local ).Normal;
			var centre = corners.Aggregate( Vector3.Zero, ( total, point ) => total + point ) / corners.Length;

			yield return new ArchOverlapFace {
				Mesh = mesh,
				Handle = handle,
				Corners = corners,
				Normal = normal,
				Centre = centre,
				Area = Area( corners ),
				Part = partIndex,
				Plane = PlaneKey( normal, centre )
			};
		}
	}

	internal int Resolve( List<ArchOverlapFace> faces ) {
		var groups = faces.GroupBy( face => face.Plane );
		var dead = new Dictionary<PolygonMesh, HashSet<FaceHandle>>();

		foreach ( var group in groups ) {
			var candidates = group.ToList();
			var extents = candidates.Select( Bounds ).ToArray();

			foreach ( var (first, second) in Pairs( extents ) ) {
				var a = candidates[first];
				var b = candidates[second];

				if ( a.Part == b.Part && ReferenceEquals( a.Mesh, b.Mesh ) && (IsDead( dead, a ) || IsDead( dead, b )) ) {
					continue;
				}

				var remove = Redundant( a, b );

				if ( remove is null || IsDead( dead, remove ) ) {
					continue;
				}

				(dead.TryGetValue( remove.Mesh, out var handles ) ? handles : dead[remove.Mesh] = new HashSet<FaceHandle>()).Add( remove.Handle );
			}
		}

		foreach ( var entry in dead ) {
			entry.Key.RemoveFaces( entry.Value.ToList() );
		}

		return dead.Values.Sum( handles => handles.Count );
	}

	// Pairs whose extents meet, in nested-walk order — order is load-bearing for tie-breaking.
	internal static List<(int First, int Second)> Pairs( BBox[] extents ) {
		var pairs = new List<(int, int)>();

		if ( extents.Length < 2 ) {
			return pairs;
		}

		var axis = Widest( extents );
		var order = Enumerable.Range( 0, extents.Length ).OrderBy( index => Along( extents[index].Mins, axis ) ).ToArray();

		for ( var lead = 0; lead < order.Length; lead++ ) {
			var first = order[lead];
			var reach = Along( extents[first].Maxs, axis ) + Slack;

			for ( var trail = lead + 1; trail < order.Length; trail++ ) {
				var second = order[trail];

				// Ordered by where each extent BEGINS, so once one begins past where this one ends, so does every
				// one behind it.
				if ( Along( extents[second].Mins, axis ) > reach ) {
					break;
				}

				if ( Touches( extents[first], extents[second] ) ) {
					pairs.Add( first < second ? (first, second) : (second, first) );
				}
			}
		}

		pairs.Sort();

		return pairs;
	}

	static int Widest( BBox[] extents ) {
		var nearest = extents[0].Mins;
		var furthest = extents[0].Mins;

		foreach ( var extent in extents ) {
			nearest = Vector3.Min( nearest, extent.Mins );
			furthest = Vector3.Max( furthest, extent.Mins );
		}

		var spread = furthest - nearest;

		if ( spread.x >= spread.y && spread.x >= spread.z ) {
			return 0;
		}

		return spread.y >= spread.z ? 1 : 2;
	}

	static float Along( Vector3 point, int axis ) => axis switch { 0 => point.x, 1 => point.y, _ => point.z };

	internal static (long, long, long, long) PlaneKey( Vector3 normal, Vector3 centre ) {
		var canonical = Canonical( normal );

		return (
			(long)MathF.Round( canonical.x / PlaneAngle ),
			(long)MathF.Round( canonical.y / PlaneAngle ),
			(long)MathF.Round( canonical.z / PlaneAngle ),
			(long)MathF.Round( Vector3.Dot( canonical, centre ) / PlaneOffset ));
	}

	// Only a face its partner covers WHOLE - a partial dip leaves the rest out in the open.
	static ArchOverlapFace Redundant( ArchOverlapFace a, ArchOverlapFace b ) {
		var normal = Canonical( a.Normal );
		var flatA = Project( a.Corners, normal, a.Centre );
		var flatB = Project( b.Corners, normal, a.Centre );

		var aCovered = Covers( flatB, flatA );
		var bCovered = Covers( flatA, flatB );

		if ( aCovered && bCovered ) {
			return a.Area <= b.Area ? a : b;
		}

		return aCovered ? a : bCovered ? b : null;
	}

	internal static BBox Bounds( ArchOverlapFace face ) => BBox.FromPoints( face.Corners );

	const float Slack = 0.5f;

	internal static bool Touches( BBox a, BBox b ) {
		return a.Mins.x <= b.Maxs.x + Slack && a.Maxs.x >= b.Mins.x - Slack
			&& a.Mins.y <= b.Maxs.y + Slack && a.Maxs.y >= b.Mins.y - Slack
			&& a.Mins.z <= b.Maxs.z + Slack && a.Maxs.z >= b.Mins.z - Slack;
	}

	static bool Covers( IReadOnlyList<Vector2> cover, IReadOnlyList<Vector2> face ) {
		return Samples( face ).All( point => Contains( cover, point ) );
	}

	static bool IsDead( Dictionary<PolygonMesh, HashSet<FaceHandle>> dead, ArchOverlapFace face ) {
		return dead.TryGetValue( face.Mesh, out var handles ) && handles.Contains( face.Handle );
	}

	internal static Vector3 Canonical( Vector3 normal ) {
		if ( normal.x < -0.001f || normal.x is > -0.001f and < 0.001f && normal.y < -0.001f || normal.x is > -0.001f and < 0.001f && normal.y is > -0.001f and < 0.001f && normal.z < 0f ) {
			return -normal;
		}

		return normal;
	}

	static Vector2[] Project( IReadOnlyList<Vector3> corners, Vector3 normal, Vector3 origin ) {
		var result = new Vector2[corners.Count];
		var axis = MathF.Abs( normal.x ) >= MathF.Abs( normal.y ) && MathF.Abs( normal.x ) >= MathF.Abs( normal.z ) ? 0 : MathF.Abs( normal.y ) >= MathF.Abs( normal.z ) ? 1 : 2;

		for ( var index = 0; index < corners.Count; index++ ) {
			var point = corners[index] - origin;
			result[index] = axis switch {
				0 => new Vector2( point.y, point.z ),
				1 => new Vector2( point.x, point.z ),
				_ => new Vector2( point.x, point.y )
			};
		}

		return result;
	}

	static IEnumerable<Vector2> Samples( IReadOnlyList<Vector2> corners ) {
		var centre = corners.Aggregate( Vector2.Zero, ( total, point ) => total + point ) / corners.Count;
		yield return centre;

		foreach ( var corner in corners ) {
			yield return Vector2.Lerp( corner, centre, 0.25f );
			yield return Vector2.Lerp( corner, centre, 0.6f );
		}
	}

	static bool Contains( IReadOnlyList<Vector2> loop, Vector2 point ) => ArchFootprint.Contains( loop, point );

	static float Area( IReadOnlyList<Vector3> corners ) {
		var normal = Vector3.Zero;

		for ( var index = 0; index < corners.Count; index++ ) {
			normal += Vector3.Cross( corners[index], corners[(index + 1) % corners.Count] );
		}

		return normal.Length * 0.5f;
	}
}