namespace AdvancedMappingTools;

/// <summary>
/// Cuts one mesh out of another. The boolean of the engine is clean for simple convex shapes but
/// breaks on a curved cutter, so this does it itself: both meshes become BSP trees, each is
/// clipped by the other, and the pieces are joined into one closed mesh again. Pieces of two
/// faces rarely end on the same points, so corners are added along the edges until every
/// edge has a partner. Faces the cut does not touch stay exactly as they were, with their
/// material and texture placement. It works in double precision, a curved tunnel makes many
/// nearly parallel planes.
/// </summary>
public static class MeshCutter
{
	const double PlaneEpsilon = 5e-4;
	const double WeldDistance = 0.01;
	const double OnEdgeDistance = 0.01;
	const double CrumbArea = 0.01;
	const int CrumbPasses = 4;
	const double SameDirection = 500.0;
	const double TurnEpsilon = 1e-7;
	const double StraightCorner = 0.02;
	const int JoinPasses = 4;
	const int BlisterPasses = 4;
	const int NeedlePasses = 8;
	const int CrackPasses = 4;
	const int CrackPoints = 96;
	const double CrackLength = 64.0;
	const double NeedleWidth = 0.01;
	const double NeedleEdge = 0.05;
	const double BlisterArea = 16.0;
	const int StackSize = 256 * 1024 * 1024;

	readonly record struct Vec( double X, double Y, double Z )
	{
		public static Vec operator +( Vec a, Vec b ) => new( a.X + b.X, a.Y + b.Y, a.Z + b.Z );
		public static Vec operator -( Vec a, Vec b ) => new( a.X - b.X, a.Y - b.Y, a.Z - b.Z );
		public static Vec operator *( Vec a, double s ) => new( a.X * s, a.Y * s, a.Z * s );

		public double Dot( Vec b ) => X * b.X + Y * b.Y + Z * b.Z;
		public Vec Cross( Vec b ) => new( Y * b.Z - Z * b.Y, Z * b.X - X * b.Z, X * b.Y - Y * b.X );
		public double Length => Math.Sqrt( Dot( this ) );
		public Vec Normal => this * (1.0 / Math.Max( Length, 1e-300 ));

		public static Vec From( Vector3 v ) => new( v.x, v.y, v.z );
		public Vector3 ToVector() => new( (float)X, (float)Y, (float)Z );
	}

	/// <summary>
	/// A record of how a triangle was divided. Tree planes run on forever and split far more than
	/// the cut needs: where both halves of a split survive, they become one again.
	/// <see cref="Face"/> is the face it came from in the mesh being cut, or -1 for the cutter.
	/// </summary>
	sealed class Fragment
	{
		public Vec[] Points;
		public int Face;
		public int Plate;
		public Vec Direction;
		public Fragment First;
		public Fragment Second;
		public bool Kept;

		/// <summary>
		/// True when all of this fragment survived. Halves that did are put back together here.
		/// </summary>
		public bool Settle()
		{
			if ( First is null ) return Kept;

			var first = First.Settle();
			var second = Second.Settle();

			if ( !first || !second ) return false;

			First = null;
			Second = null;
			Kept = true;

			return true;
		}

		public void Collect( List<Fragment> into )
		{
			if ( First is not null )
			{
				First.Collect( into );
				Second.Collect( into );
			}
			else if ( Kept )
			{
				into.Add( this );
			}
		}
	}

	/// <summary>
	/// A convex piece of a face on its way through the trees.
	/// </summary>
	sealed class Piece
	{
		public Vec[] Points;
		public Vec Normal;
		public double Offset;
		public Fragment Fragment;
		public bool Inverted;

		public Piece Flipped()
		{
			var points = (Vec[])Points.Clone();
			Array.Reverse( points );

			return new Piece { Points = points, Normal = Normal * -1.0, Offset = -Offset, Fragment = Fragment, Inverted = !Inverted };
		}

		public Piece Part( Vec[] points, Fragment fragment ) =>
			new() { Points = points, Normal = Normal, Offset = Offset, Fragment = fragment, Inverted = Inverted };

		/// <summary>
		/// A record for a part of this piece, in the direction of the original face.
		/// </summary>
		public Fragment Record( Vec[] points )
		{
			var stored = (Vec[])points.Clone();

			if ( Inverted ) Array.Reverse( stored );

			return new Fragment { Points = stored, Face = Fragment.Face, Plate = Fragment.Plate, Direction = Fragment.Direction };
		}
	}

	const int Coplanar = 0;
	const int Front = 1;
	const int Back = 2;
	const int Spanning = 3;

	static void Split( Vec normal, double offset, Piece piece, List<Piece> coplanarFront, List<Piece> coplanarBack, List<Piece> front, List<Piece> back )
	{
		var count = piece.Points.Length;
		var sides = new int[count];
		var kind = 0;

		for ( var i = 0; i < count; i++ )
		{
			var t = normal.Dot( piece.Points[i] ) - offset;
			var side = t < -PlaneEpsilon ? Back : t > PlaneEpsilon ? Front : Coplanar;

			kind |= side;
			sides[i] = side;
		}

		switch ( kind )
		{
			case Coplanar:
				(normal.Dot( piece.Normal ) > 0.0 ? coplanarFront : coplanarBack).Add( piece );
				return;

			case Front:
				front.Add( piece );
				return;

			case Back:
				back.Add( piece );
				return;
		}

		var f = new List<Vec>( count + 2 );
		var b = new List<Vec>( count + 2 );

		for ( var i = 0; i < count; i++ )
		{
			var j = (i + 1) % count;
			var si = sides[i];
			var sj = sides[j];
			var vi = piece.Points[i];
			var vj = piece.Points[j];

			if ( si != Back ) f.Add( vi );
			if ( si != Front ) b.Add( vi );

			if ( (si | sj) != Spanning ) continue;

			var t = (offset - normal.Dot( vi )) / normal.Dot( vj - vi );
			var v = vi + (vj - vi) * t;

			f.Add( v );
			b.Add( v );
		}

		// Only a real division is recorded, a sliver on one side leaves the piece as it was.
		if ( f.Count < 3 || b.Count < 3 )
		{
			(f.Count >= 3 ? front : back).Add( piece );
			return;
		}

		var frontPoints = f.ToArray();
		var backPoints = b.ToArray();
		var frontRecord = piece.Record( frontPoints );
		var backRecord = piece.Record( backPoints );

		piece.Fragment.First = frontRecord;
		piece.Fragment.Second = backRecord;

		front.Add( piece.Part( frontPoints, frontRecord ) );
		back.Add( piece.Part( backPoints, backRecord ) );
	}

	sealed class Node
	{
		Vec _normal;
		double _offset;
		bool _hasPlane;
		Node _front;
		Node _back;
		readonly List<Piece> _pieces = [];

		public void Build( List<Piece> pieces )
		{
			if ( pieces.Count == 0 ) return;

			if ( !_hasPlane )
			{
				_normal = pieces[0].Normal;
				_offset = pieces[0].Offset;
				_hasPlane = true;
			}

			var front = new List<Piece>();
			var back = new List<Piece>();

			foreach ( var piece in pieces )
				Split( _normal, _offset, piece, _pieces, _pieces, front, back );

			if ( front.Count > 0 ) (_front ??= new Node()).Build( front );
			if ( back.Count > 0 ) (_back ??= new Node()).Build( back );
		}

		public void Invert()
		{
			for ( var i = 0; i < _pieces.Count; i++ )
				_pieces[i] = _pieces[i].Flipped();

			_normal *= -1.0;
			_offset = -_offset;

			_front?.Invert();
			_back?.Invert();

			(_front, _back) = (_back, _front);
		}

		public List<Piece> Clip( List<Piece> pieces )
		{
			if ( !_hasPlane ) return [.. pieces];

			var front = new List<Piece>();
			var back = new List<Piece>();

			foreach ( var piece in pieces )
				Split( _normal, _offset, piece, front, back, front, back );

			if ( _front is not null ) front = _front.Clip( front );

			if ( _back is not null ) back = _back.Clip( back );
			else back = [];

			front.AddRange( back );

			return front;
		}

		public void ClipTo( Node other )
		{
			var clipped = other.Clip( _pieces );

			_pieces.Clear();
			_pieces.AddRange( clipped );

			_front?.ClipTo( other );
			_back?.ClipTo( other );
		}

		public void Collect( List<Piece> into )
		{
			into.AddRange( _pieces );

			_front?.Collect( into );
			_back?.Collect( into );
		}
	}

	static Piece MakePiece( Vec a, Vec b, Vec c, int face, int plate )
	{
		var normal = (b - a).Cross( c - a );

		if ( normal.Length < 1e-9 ) return null;

		normal = normal.Normal;

		return new Piece { Points = [a, b, c], Normal = normal, Offset = normal.Dot( a ), Fragment = new Fragment { Points = [a, b, c], Face = face, Plate = plate, Direction = normal } };
	}

	static List<Piece> Pieces( PolygonMesh mesh, Transform transform, bool keepFaces, out List<HalfEdgeMesh.FaceHandle> handles )
	{
		var pieces = new List<Piece>();

		handles = [.. mesh.FaceHandles];

		if ( handles.Count > 0 && mesh.CreateFace( handles[0], transform, Color.White ) is null )
			mesh.Rebuild();

		for ( var f = 0; f < handles.Count; f++ )
		{
			var vertices = mesh.CreateFace( handles[f], transform, Color.White );

			if ( vertices is null ) continue;

			for ( var i = 0; i + 2 < vertices.Length; i += 3 )
			{
				var piece = MakePiece( Vec.From( vertices[i].Position ), Vec.From( vertices[i + 1].Position ), Vec.From( vertices[i + 2].Position ), keepFaces ? f : -1, keepFaces ? f : -1 - f );

				if ( piece is not null ) pieces.Add( piece );
			}
		}

		return pieces;
	}

	/// <summary>
	/// A fixed shuffle. Faces of a swept shape lie in rings, and a tree built in that order is
	/// as deep as the mesh is large and splits far more than it has to.
	/// </summary>
	static List<Piece> Shuffled( List<Piece> pieces )
	{
		var result = new List<Piece>( pieces );
		var state = 2463534242u;

		for ( var i = result.Count - 1; i > 0; i-- )
		{
			state ^= state << 13;
			state ^= state >> 17;
			state ^= state << 5;

			var j = (int)(state % (uint)(i + 1));

			(result[i], result[j]) = (result[j], result[i]);
		}

		return result;
	}

	/// <summary>
	/// A closed mesh has a partner for every edge. Only then has it an inside, and only then
	/// does a cut through it get walls of its own. A sheet just gets a hole.
	/// </summary>
	static bool IsClosed( PolygonMesh mesh ) => !mesh.HalfEdgeHandles.Any( mesh.IsEdgeOpen );

	/// <summary>
	/// The mesh with the cutter taken out of it, or nothing if the two do not touch. The cutter
	/// is given in its own space, <paramref name="relative"/> takes it into that of the mesh.
	/// New faces inside the cut get <paramref name="material"/>.
	/// </summary>
	public static PolygonMesh Subtract( PolygonMesh mesh, PolygonMesh cutter, Transform relative, Material material )
	{
		// Bounds of the vertices themselves. A mesh made in the editor carries a transform of its own
		// for its texture mapping, which CalculateBounds would add; the cut works in vertex space.
		var meshBounds = BBox.FromPoints( mesh.VertexHandles.Select( v => mesh.GetVertexPosition( v ) ) );
		var cutterBounds = BBox.FromPoints( cutter.VertexHandles.Select( v => relative.PointToWorld( cutter.GetVertexPosition( v ) ) ) );

		if ( !meshBounds.Overlaps( cutterBounds ) ) return null;

		var solid = IsClosed( mesh );
		var body = Pieces( mesh, new Transform( Vector3.Zero, Rotation.Identity ), true, out var faces );
		var tool = Pieces( cutter, relative, false, out _ );
		var roots = body.Concat( tool ).Select( x => x.Fragment ).ToList();

		// Trees go as deep as they like, so this runs on a stack that is large enough for it.
		var worker = new System.Threading.Thread( () =>
		{
			foreach ( var piece in Cut( body, tool, solid ) )
				piece.Fragment.Kept = true;

			foreach ( var root in roots )
				root.Settle();
		}, StackSize );

		worker.Start();
		worker.Join();

		return Assemble( mesh, faces, roots, material );
	}

	/// <summary>
	/// A cut that passes a corner by a hair can come out with a crack. Moved by less than anyone
	/// can see the same cut misses that corner, so where the result is not clean it is tried
	/// again a little beside itself, and the cleanest result counts.
	/// </summary>
	public static PolygonMesh SubtractClean( PolygonMesh mesh, PolygonMesh cutter, Transform relative, Material material )
	{
		var best = default( PolygonMesh );
		var fewest = int.MaxValue;
		var closed = IsClosed( mesh );

		foreach ( var nudge in Nudges )
		{
			var moved = relative;

			moved.Position += nudge;

			var result = Subtract( mesh, cutter, moved, material );

			if ( result is null ) return best;

			var report = MeshDiagnostics.Analyze( result );
			var faults = report.Errors + report.Warnings;

			// A hole in an open sheet has an open edge all around it, as it should.
			if ( !closed ) faults -= report.Count( MeshIssueKind.OpenEdge ) + report.Count( MeshIssueKind.UnweldedSeam );

			if ( faults < fewest )
			{
				fewest = faults;
				best = result;
			}

			if ( faults <= 0 ) break;
		}

		return best;
	}

	static readonly Vector3[] Nudges =
	[
		Vector3.Zero,
		new( 0.031f, 0.017f, 0.023f ),
		new( -0.027f, 0.041f, -0.013f ),
		new( 0.019f, -0.036f, 0.044f ),
		new( -0.043f, -0.021f, -0.037f ),
		new( 0.053f, 0.047f, -0.029f )
	];

	static List<Piece> Cut( List<Piece> body, List<Piece> tool, bool solid )
	{
		// Subtracting is adding the cutter to the inverted mesh and inverting that.
		var a = new Node();
		var b = new Node();

		a.Build( Shuffled( body ) );
		b.Build( Shuffled( tool ) );

		a.Invert();
		a.ClipTo( b );

		var kept = new List<Piece>();

		if ( solid )
		{
			b.ClipTo( a );
			b.Invert();
			b.ClipTo( a );
			b.Invert();
			b.Collect( kept );
		}

		a.Collect( kept );

		return kept;
	}

	sealed class Ring
	{
		public int[] Points;
		public int Face;
		public int Plate;
		public Vec Direction;
	}

	/// <summary>
	/// Rebuild a mesh from what is left. A face that lost nothing comes back as it was, the rest
	/// is pieces. Every piece then gets a corner wherever another piece ends on one of its edges.
	/// </summary>
	static PolygonMesh Assemble( PolygonMesh source, List<HalfEdgeMesh.FaceHandle> faces, List<Fragment> roots, Material material )
	{
		var touched = new HashSet<int>();
		var pieces = new List<Fragment>();

		foreach ( var root in roots )
		{
			if ( root.Face >= 0 && !root.Kept ) touched.Add( root.Face );

			root.Collect( pieces );
		}

		// Nothing was taken away and nothing was added: the cutter passed by.
		if ( touched.Count == 0 && !pieces.Any( x => x.Face < 0 ) ) return null;

		var points = new List<Vec>();
		var grid = new Dictionary<(long, long, long), int>();

		int Point( Vec p )
		{
			var key = ((long)Math.Round( p.X / WeldDistance ), (long)Math.Round( p.Y / WeldDistance ), (long)Math.Round( p.Z / WeldDistance ));

			for ( var x = -1; x <= 1; x++ )
				for ( var y = -1; y <= 1; y++ )
					for ( var z = -1; z <= 1; z++ )
					{
						if ( !grid.TryGetValue( (key.Item1 + x, key.Item2 + y, key.Item3 + z), out var near ) ) continue;
						if ( (points[near] - p).Length <= WeldDistance ) return near;
					}

			grid[key] = points.Count;
			points.Add( p );

			return points.Count - 1;
		}

		var rings = new List<Ring>();

		void AddRing( IEnumerable<Vec> positions, int face, int plate, Vec direction )
		{
			var ring = new List<int>();

			foreach ( var position in positions )
			{
				var id = Point( position );

				if ( ring.Count == 0 || ring[^1] != id ) ring.Add( id );
			}

			if ( ring.Count > 1 && ring[^1] == ring[0] ) ring.RemoveAt( ring.Count - 1 );
			if ( ring.Count >= 3 ) rings.Add( new Ring { Points = [.. ring], Face = face, Plate = plate, Direction = direction } );
		}

		for ( var f = 0; f < faces.Count; f++ )
		{
			if ( touched.Contains( f ) ) continue;

			AddRing( source.GetFaceVertices( faces[f] ).Select( x => Vec.From( source.GetVertexPosition( x ) ) ), f, f, new Vec() );
		}

		foreach ( var piece in pieces )
		{
			if ( piece.Face >= 0 && !touched.Contains( piece.Face ) ) continue;

			// What is left of the cutter looks into the cut.
			AddRing( piece.Face >= 0 ? piece.Points : piece.Points.Reverse(), piece.Face, piece.Plate, piece.Face >= 0 ? piece.Direction : piece.Direction * -1.0 );
		}

		FillEdges( rings, points );

		foreach ( var ring in rings )
			ring.Points = WithoutSpikes( ring.Points );

		rings.RemoveAll( x => x.Points.Length < 3 );

		CollapseCrumbs( rings, points );

		rings = Merge( rings, points );
		rings = [.. rings.SelectMany( x => x.Face >= 0 && !touched.Contains( x.Face ) ? [x] : Split( x, points ) )];
		rings = WithoutBlisters( rings, points );
		rings = WithoutNeedles( rings, points );

		for ( var pass = 0; pass < CrackPasses; pass++ )
		{
			var before = rings.Count;

			rings = WithoutCracks( rings, points, IsClosed( source ) ? double.MaxValue : CrackLength );
			rings = WithoutNeedles( rings, points );

			if ( rings.Count == before ) break;
		}

		rings = WithoutPillows( rings );

		return Build( source, faces, rings, points, material );
	}

	/// <summary>
	/// Two faces on the same corners, back to back, have nothing between them and nothing to show.
	/// </summary>
	static List<Ring> WithoutPillows( List<Ring> rings )
	{
		var gone = new HashSet<Ring>();

		foreach ( var same in rings.GroupBy( x => string.Join( ' ', x.Points.OrderBy( p => p ) ) ).Where( x => x.Count() > 1 ) )
		{
			var group = same.ToList();

			for ( var i = 0; i < group.Count; i++ )
			{
				for ( var j = i + 1; j < group.Count; j++ )
				{
					var a = group[i].Points;
					var b = group[j].Points;

					if ( gone.Contains( group[i] ) || gone.Contains( group[j] ) ) continue;

					var at = Array.IndexOf( b, a[0] );

					if ( b[(at + b.Length - 1) % b.Length] != a[1] ) continue;

					gone.Add( group[i] );
					gone.Add( group[j] );
				}
			}
		}

		return gone.Count == 0 ? rings : [.. rings.Where( x => !gone.Contains( x ) )];
	}

	/// <summary>
	/// What is still open at this point is a crack a few points around, between edges that missed
	/// each other. It gets a face of its own. Most of those have no width and go again as needles,
	/// which is how a corner on one side of a crack finally reaches the other side. A body that
	/// was closed has no open edges at all, and the real edge of an open sheet is far longer.
	/// </summary>
	static List<Ring> WithoutCracks( List<Ring> rings, List<Vec> points, double longest )
	{
		var used = new HashSet<(int, int)>();

		foreach ( var ring in rings )
		{
			for ( var i = 0; i < ring.Points.Length; i++ )
				used.Add( (ring.Points[i], ring.Points[(i + 1) % ring.Points.Length]) );
		}

		var open = used.Where( x => !used.Contains( (x.Item2, x.Item1) ) ).ToList();

		if ( open.Count == 0 ) return rings;

		var arriving = new Dictionary<int, List<int>>();

		foreach ( var (from, to) in open )
		{
			if ( !arriving.TryGetValue( to, out var list ) ) arriving[to] = list = [];

			list.Add( from );
		}

		var done = new HashSet<(int, int)>();
		var added = new List<Ring>();

		foreach ( var (first, second) in open )
		{
			if ( !done.Add( (first, second) ) ) continue;

			// The new face runs against the open edges.
			var loop = new List<int> { second, first };
			var current = first;
			var around = (points[first] - points[second]).Length;
			var closed = false;

			while ( loop.Count < CrackPoints && around < longest )
			{
				if ( !arriving.TryGetValue( current, out var options ) ) break;

				var next = options.FindIndex( x => !done.Contains( (x, current) ) );

				if ( next < 0 ) break;

				var point = options[next];

				done.Add( (point, current) );
				around += (points[point] - points[current]).Length;

				if ( point == second )
				{
					closed = true;
					break;
				}

				loop.Add( point );
				current = point;
			}

			if ( closed && around < longest )
				Divide( loop, -1, points, added );
		}

		return [.. rings, .. added];
	}

	/// <summary>
	/// Triangles without any width are left where three points lie on a line. With a short edge such
	/// a needle gives up that edge. Without one its middle point lies on the long edge, so the
	/// face on the other side of that edge takes the point over, and the needle has no reason to exist.
	/// </summary>
	static List<Ring> WithoutNeedles( List<Ring> rings, List<Vec> points )
	{
		for ( var pass = 0; pass < NeedlePasses; pass++ )
		{
			var edges = new Dictionary<(int, int), Ring>();

			foreach ( var ring in rings )
			{
				for ( var i = 0; i < ring.Points.Length; i++ )
					edges[(ring.Points[i], ring.Points[(i + 1) % ring.Points.Length])] = ring;
			}

			var moved = new Dictionary<int, int>();
			var gone = new HashSet<Ring>();
			var grown = new HashSet<Ring>();

			foreach ( var ring in rings )
			{
				if ( ring.Points.Length != 3 || gone.Contains( ring ) || grown.Contains( ring ) ) continue;

				var sides = Enumerable.Range( 0, 3 )
					.Select( i => (A: ring.Points[i], B: ring.Points[(i + 1) % 3], Other: ring.Points[(i + 2) % 3]) )
					.OrderBy( x => (points[x.A] - points[x.B]).Length )
					.ToArray();

				var longest = sides[2];
				var length = (points[longest.A] - points[longest.B]).Length;

				var area = Area( ring.Points, points );

				if ( area >= CrumbArea && area * 2.0 / Math.Max( length, 1e-12 ) >= NeedleWidth ) continue;

				var shortest = sides[0];

				if ( (points[shortest.A] - points[shortest.B]).Length < NeedleEdge )
				{
					if ( !moved.ContainsKey( shortest.A ) && !moved.ContainsValue( shortest.A ) && !moved.ContainsKey( shortest.B ) )
						moved[shortest.A] = shortest.B;

					continue;
				}

				if ( !edges.TryGetValue( (longest.B, longest.A), out var neighbour ) || gone.Contains( neighbour ) || grown.Contains( neighbour ) ) continue;

				var at = Array.IndexOf( neighbour.Points, longest.B );

				if ( at < 0 || neighbour.Points[(at + 1) % neighbour.Points.Length] != longest.A ) continue;

				var taken = neighbour.Points.ToList();

				taken.Insert( at + 1, longest.Other );
				neighbour.Points = [.. taken];

				grown.Add( neighbour );
				gone.Add( ring );
			}

			if ( moved.Count == 0 && gone.Count == 0 ) break;

			rings = [.. rings.Where( x => !gone.Contains( x ) )];

			foreach ( var ring in rings )
				ring.Points = WithoutSpikes( [.. ring.Points.Select( x => moved.TryGetValue( x, out var to ) ? to : x )] );

			rings = [.. rings.Where( x => x.Points.Length >= 3 ).SelectMany( x => grown.Contains( x ) ? Split( x, points ) : [x] )];
			rings = WithoutBlisters( rings, points );
		}

		return rings;
	}

	/// <summary>
	/// Where a cut runs along a face at the width of a hair, the two sides can leave a blister: thin
	/// faces back to back, closed in on themselves and hanging on the real surface by edges
	/// that then have four faces instead of two. Such an edge has no place in a mesh.
	/// A real part of a body never hangs on like that, and these are far too small to be one.
	/// </summary>
	static List<Ring> WithoutBlisters( List<Ring> rings, List<Vec> points )
	{
		for ( var pass = 0; pass < BlisterPasses; pass++ )
		{
			var uses = new Dictionary<(int, int), List<int>>();

			for ( var r = 0; r < rings.Count; r++ )
			{
				var ring = rings[r].Points;

				for ( var i = 0; i < ring.Length; i++ )
				{
					var key = (ring[i], ring[(i + 1) % ring.Length]);

					if ( !uses.TryGetValue( key, out var list ) ) uses[key] = list = [];

					list.Add( r );
				}
			}

			var parents = Enumerable.Range( 0, rings.Count ).ToArray();

			int Find( int x )
			{
				while ( parents[x] != x ) x = parents[x] = parents[parents[x]];

				return x;
			}

			var pinched = new HashSet<int>();

			foreach ( var ((a, b), forward) in uses )
			{
				uses.TryGetValue( (b, a), out var backward );

				if ( forward.Count == 1 && backward?.Count == 1 )
				{
					parents[Find( forward[0] )] = Find( backward[0] );
				}
				else if ( forward.Count > 1 || (backward?.Count ?? 0) > 1 )
				{
					pinched.UnionWith( forward );

					if ( backward is not null ) pinched.UnionWith( backward );
				}
			}

			if ( pinched.Count == 0 ) return rings;

			var areas = new Dictionary<int, double>();

			for ( var r = 0; r < rings.Count; r++ )
				areas[Find( r )] = areas.GetValueOrDefault( Find( r ) ) + Area( rings[r].Points, points );

			var blisters = pinched.Select( Find ).Where( x => areas[x] < BlisterArea ).ToHashSet();

			if ( blisters.Count == 0 ) return rings;

			rings = [.. rings.Where( ( _, r ) => !blisters.Contains( Find( r ) ) )];
		}

		return rings;
	}

	/// <summary>
	/// Tree planes run on forever, so a face comes out of the trees in many pieces that were
	/// divided for no reason, and those are where slivers come from. All pieces of a face in one plane
	/// become one again: what they share goes, and what is left is the true outline of what
	/// survived the cut, with any holes it has.
	/// </summary>
	static List<Ring> Merge( List<Ring> rings, List<Vec> points )
	{
		var result = new List<Ring>( rings.Count );

		foreach ( var plate in rings.GroupBy( x => (x.Plate, (long)Math.Round( x.Direction.X * SameDirection ), (long)Math.Round( x.Direction.Y * SameDirection ), (long)Math.Round( x.Direction.Z * SameDirection )) ) )
		{
			var group = plate.ToList();

			if ( group.Count == 1 )
			{
				result.AddRange( group );
				continue;
			}

			result.AddRange( Outlines( group, points ) ?? group );
		}

		return result;
	}

	/// <summary>
	/// The outlines of a group of rings in one plane as new rings, or nothing if they do not make any.
	/// </summary>
	static List<Ring> Outlines( List<Ring> group, List<Vec> points )
	{
		var edges = new HashSet<(int From, int To)>();

		foreach ( var ring in group )
		{
			for ( var i = 0; i < ring.Points.Length; i++ )
			{
				var a = ring.Points[i];
				var b = ring.Points[(i + 1) % ring.Points.Length];

				if ( edges.Remove( (b, a) ) ) continue;
				if ( !edges.Add( (a, b) ) ) return null;
			}
		}

		if ( edges.Count < 3 ) return null;

		var first = group[0];
		var normal = first.Direction;
		var side = Math.Abs( normal.X ) < 0.7 ? new Vec( 1.0, 0.0, 0.0 ) : new Vec( 0.0, 1.0, 0.0 );
		var axisU = normal.Cross( side ).Normal;
		var axisV = normal.Cross( axisU ).Normal;

		(double X, double Y) Flat( int point ) => (points[point].Dot( axisU ), points[point].Dot( axisV ));

		var leaving = new Dictionary<int, List<int>>();

		foreach ( var (from, to) in edges )
		{
			if ( !leaving.TryGetValue( from, out var list ) ) leaving[from] = list = [];

			list.Add( to );
		}

		var loops = new List<List<int>>();

		while ( edges.Count > 0 )
		{
			var (start, next) = edges.First();
			var loop = new List<int> { start };
			var before = start;
			var current = next;

			edges.Remove( (start, next) );
			leaving[start].Remove( next );

			while ( current != start )
			{
				loop.Add( current );

				if ( !leaving.TryGetValue( current, out var options ) || options.Count == 0 ) return null;

				// Where two outlines touch in one point this keeps each of them to itself.
				var chosen = options[0];

				if ( options.Count > 1 )
				{
					var (px, py) = Flat( before );
					var (cx, cy) = Flat( current );
					var incoming = Math.Atan2( cy - py, cx - px );

					chosen = options.OrderBy( option =>
					{
						var (nx, ny) = Flat( option );
						var turn = Math.Atan2( ny - cy, nx - cx ) - incoming;

						while ( turn <= -Math.PI ) turn += Math.PI * 2.0;
						while ( turn > Math.PI ) turn -= Math.PI * 2.0;

						return turn;
					} ).First();
				}

				edges.Remove( (current, chosen) );
				options.Remove( chosen );

				before = current;
				current = chosen;
			}

			if ( loop.Count >= 3 ) loops.Add( loop );
		}

		double Signed( List<int> loop )
		{
			var sum = 0.0;

			for ( var i = 0; i < loop.Count; i++ )
			{
				var (ax, ay) = Flat( loop[i] );
				var (bx, by) = Flat( loop[(i + 1) % loop.Count] );

				sum += ax * by - bx * ay;
			}

			return sum * 0.5;
		}

		var outers = loops.Where( x => Signed( x ) > 0.0 ).ToList();
		var holes = loops.Where( x => Signed( x ) <= 0.0 ).ToList();

		if ( outers.Count == 0 ) return null;

		var result = new List<Ring>();
		var wanted = 0.0;
		var made = 0.0;

		double Corners( int a, int b, int c )
		{
			var (ax, ay) = Flat( a );
			var (bx, by) = Flat( b );
			var (cx, cy) = Flat( c );

			return ((bx - ax) * (cy - ay) - (by - ay) * (cx - ax)) * 0.5;
		}

		foreach ( var outer in outers )
		{
			var loop = outer;

			wanted += Signed( outer );

			foreach ( var hole in holes.Where( h => Encloses( outer, Flat( h[0] ), Flat ) ).OrderByDescending( h => Signed( h ) ).ToList() )
			{
				wanted += Signed( hole );
				loop = Bridged( loop, hole, Flat );

				if ( loop is null ) return null;

				holes.Remove( hole );
			}

			var plain = Plain( loop, Flat );
			var triangles = plain ? [[.. loop]] : EarClip( loop, Flat );

			if ( triangles is null ) return null;

			foreach ( var ring in triangles )
			{
				for ( var i = 1; i + 1 < ring.Length; i++ )
				{
					var area = Corners( ring[0], ring[i], ring[i + 1] );

					if ( area <= 0.0 ) return null;

					made += area;
				}

				result.Add( new Ring { Points = ring, Face = first.Face, Plate = first.Plate, Direction = first.Direction } );
			}
		}

		// What comes out has to cover what went in, or the outlines were not what they seemed.
		if ( holes.Count != 0 || Math.Abs( made - wanted ) > Math.Max( wanted * 1e-4, 1e-3 ) ) return null;

		return result;
	}

	static bool Encloses( List<int> loop, (double X, double Y) point, Func<int, (double X, double Y)> flat )
	{
		var inside = false;

		for ( var i = 0; i < loop.Count; i++ )
		{
			var (ax, ay) = flat( loop[i] );
			var (bx, by) = flat( loop[(i + 1) % loop.Count] );

			if ( ay > point.Y != by > point.Y && point.X < ax + (point.Y - ay) / (by - ay) * (bx - ax) )
				inside = !inside;
		}

		return inside;
	}

	/// <summary>
	/// A loop that is all outward corners, which a face can take as it is.
	/// </summary>
	static bool Plain( List<int> loop, Func<int, (double X, double Y)> flat )
	{
		for ( var i = 0; i < loop.Count; i++ )
		{
			var (ax, ay) = flat( loop[(i + loop.Count - 1) % loop.Count] );
			var (bx, by) = flat( loop[i] );
			var (cx, cy) = flat( loop[(i + 1) % loop.Count] );
			var first = Math.Sqrt( (bx - ax) * (bx - ax) + (by - ay) * (by - ay) );
			var second = Math.Sqrt( (cx - bx) * (cx - bx) + (cy - by) * (cy - by) );

			if ( ((bx - ax) * (cy - by) - (by - ay) * (cx - bx)) / Math.Max( first * second, 1e-12 ) < StraightCorner ) return false;
		}

		return true;
	}

	/// <summary>
	/// A hole joined to its outline by a cut from one to the other, so the two are one loop.
	/// </summary>
	static List<int> Bridged( List<int> outer, List<int> hole, Func<int, (double X, double Y)> flat )
	{
		bool Crosses( (double X, double Y) a, (double X, double Y) b, List<int> loop, int skipA, int skipB )
		{
			for ( var i = 0; i < loop.Count; i++ )
			{
				var u = loop[i];
				var v = loop[(i + 1) % loop.Count];

				if ( u == skipA || v == skipA || u == skipB || v == skipB ) continue;

				var (ux, uy) = flat( u );
				var (vx, vy) = flat( v );
				var d1 = (b.X - a.X) * (uy - a.Y) - (b.Y - a.Y) * (ux - a.X);
				var d2 = (b.X - a.X) * (vy - a.Y) - (b.Y - a.Y) * (vx - a.X);
				var d3 = (vx - ux) * (a.Y - uy) - (vy - uy) * (a.X - ux);
				var d4 = (vx - ux) * (b.Y - uy) - (vy - uy) * (b.X - ux);

				if ( d1 * d2 <= 0.0 && d3 * d4 <= 0.0 ) return true;
			}

			return false;
		}

		var best = (Outer: -1, Hole: -1);
		var shortest = double.MaxValue;

		for ( var o = 0; o < outer.Count; o++ )
		{
			var a = flat( outer[o] );

			for ( var h = 0; h < hole.Count; h++ )
			{
				var b = flat( hole[h] );
				var length = (a.X - b.X) * (a.X - b.X) + (a.Y - b.Y) * (a.Y - b.Y);

				if ( length >= shortest ) continue;
				if ( Crosses( a, b, outer, outer[o], hole[h] ) || Crosses( a, b, hole, outer[o], hole[h] ) ) continue;

				shortest = length;
				best = (o, h);
			}
		}

		if ( best.Outer < 0 ) return null;

		var joined = new List<int>( outer.Count + hole.Count + 2 );

		for ( var i = 0; i <= best.Outer; i++ ) joined.Add( outer[i] );
		for ( var i = 0; i <= hole.Count; i++ ) joined.Add( hole[(best.Hole + i) % hole.Count] );
		for ( var i = best.Outer; i < outer.Count; i++ ) joined.Add( outer[i] );

		return joined;
	}

	/// <summary>
	/// Cut a loop into triangles, the best shaped corners first. A corner of a triangle may only go
	/// when no other point of the loop lies in it or on its edges, so every corner on a straight
	/// edge stays a corner of something and no neighbour loses its partner.
	/// </summary>
	static List<int[]> EarClip( List<int> loop, Func<int, (double X, double Y)> flat )
	{
		var left = new List<int>( loop );
		var result = new List<int[]>( loop.Count );

		while ( left.Count > 3 )
		{
			var chosen = -1;
			var bestShape = 0.0;

			for ( var i = 0; i < left.Count; i++ )
			{
				var previous = left[(i + left.Count - 1) % left.Count];
				var corner = left[i];
				var following = left[(i + 1) % left.Count];

				if ( previous == following ) continue;

				var a = flat( previous );
				var b = flat( corner );
				var c = flat( following );
				var area = (b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X);
				var longest = Math.Max( Math.Max( (b.X - a.X) * (b.X - a.X) + (b.Y - a.Y) * (b.Y - a.Y), (c.X - b.X) * (c.X - b.X) + (c.Y - b.Y) * (c.Y - b.Y) ), (a.X - c.X) * (a.X - c.X) + (a.Y - c.Y) * (a.Y - c.Y) );
				var shape = area / Math.Max( longest, 1e-12 );

				if ( shape <= TurnEpsilon || shape <= bestShape ) continue;

				var blocked = false;
				var ab = Math.Sqrt( (b.X - a.X) * (b.X - a.X) + (b.Y - a.Y) * (b.Y - a.Y) );
				var bc = Math.Sqrt( (c.X - b.X) * (c.X - b.X) + (c.Y - b.Y) * (c.Y - b.Y) );
				var ca = Math.Sqrt( (a.X - c.X) * (a.X - c.X) + (a.Y - c.Y) * (a.Y - c.Y) );

				for ( var k = 0; k < left.Count && !blocked; k++ )
				{
					var other = left[k];

					if ( other == previous || other == corner || other == following ) continue;

					var p = flat( other );
					var d1 = ((b.X - a.X) * (p.Y - a.Y) - (b.Y - a.Y) * (p.X - a.X)) / Math.Max( ab, 1e-9 );
					var d2 = ((c.X - b.X) * (p.Y - b.Y) - (c.Y - b.Y) * (p.X - b.X)) / Math.Max( bc, 1e-9 );
					var d3 = ((a.X - c.X) * (p.Y - c.Y) - (a.Y - c.Y) * (p.X - c.X)) / Math.Max( ca, 1e-9 );

					blocked = d1 >= -OnEdgeDistance && d2 >= -OnEdgeDistance && d3 >= -OnEdgeDistance;
				}

				if ( blocked ) continue;

				bestShape = shape;
				chosen = i;
			}

			if ( chosen < 0 ) return null;

			result.Add( [left[(chosen + left.Count - 1) % left.Count], left[chosen], left[(chosen + 1) % left.Count]] );
			left.RemoveAt( chosen );
		}

		result.Add( [.. left] );

		return result;
	}

	/// <summary>
	/// A piece with a corner on a straight line has to be divided here: left to itself it
	/// can fall into a triangle without area.
	/// </summary>
	static IEnumerable<Ring> Split( Ring ring, List<Vec> points )
	{
		var loop = ring.Points;

		bool Straight( int i )
		{
			var a = points[loop[(i + loop.Length - 1) % loop.Length]];
			var b = points[loop[i]];
			var c = points[loop[(i + 1) % loop.Length]];

			return (b - a).Normal.Cross( (c - b).Normal ).Length < StraightCorner;
		}

		if ( !Enumerable.Range( 0, loop.Length ).Any( Straight ) ) return [ring];

		var result = new List<Ring>();

		Divide( [.. loop], ring.Face, points, result );

		return result;
	}

	/// <summary>
	/// Divide a flat piece along the diagonal that leaves both sides as far from thin as they can be.
	/// A diagonal between points of one straight edge leaves a side without area and so never wins.
	/// </summary>
	static void Divide( List<int> loop, int face, List<Vec> points, List<Ring> into )
	{
		if ( loop.Count <= 3 )
		{
			if ( loop.Count == 3 ) into.Add( new Ring { Points = [.. loop], Face = face, Plate = face } );

			return;
		}

		var count = loop.Count;
		var bestFirst = 0;
		var bestSecond = 2;
		var bestShape = double.MinValue;

		for ( var i = 0; i < count; i++ )
		{
			for ( var j = i + 2; j < count; j++ )
			{
				if ( i == 0 && j == count - 1 ) continue;

				var shape = Math.Min( Fatness( loop, i, j, points ), Fatness( loop, j, i + count, points ) );

				if ( shape > bestShape )
				{
					bestShape = shape;
					bestFirst = i;
					bestSecond = j;
				}
			}
		}

		var first = new List<int>();
		var second = new List<int>();

		for ( var k = bestFirst; k <= bestSecond; k++ ) first.Add( loop[k] );
		for ( var k = bestSecond; k <= bestFirst + count; k++ ) second.Add( loop[k % count] );

		Divide( first, face, points, into );
		Divide( second, face, points, into );
	}

	/// <summary>
	/// The area of a stretch of a loop closed by its last point, against the size of it.
	/// </summary>
	static double Fatness( List<int> loop, int from, int to, List<Vec> points )
	{
		var count = loop.Count;
		var origin = points[loop[from % count]];
		var sum = new Vec();
		var longest = 0.0;

		for ( var k = from + 1; k < to; k++ )
		{
			sum += (points[loop[k % count]] - origin).Cross( points[loop[(k + 1) % count]] - origin );
			longest = Math.Max( longest, (points[loop[k % count]] - origin).Length );
		}

		longest = Math.Max( longest, (points[loop[to % count]] - origin).Length );

		return sum.Length / Math.Max( longest * longest, 1e-12 );
	}

	/// <summary>
	/// Where a cut passes a corner very closely it leaves crumbs, faces far too small to be seen
	/// or used. Each one gives up its shortest edge, which takes the face away and leaves its
	/// neighbours connected.
	/// </summary>
	static void CollapseCrumbs( List<Ring> rings, List<Vec> points )
	{
		for ( var pass = 0; pass < CrumbPasses; pass++ )
		{
			var moved = new Dictionary<int, int>();

			foreach ( var ring in rings )
			{
				if ( Area( ring.Points, points ) >= CrumbArea ) continue;
				if ( ring.Points.Any( moved.ContainsKey ) ) continue;

				var (from, to) = Enumerable.Range( 0, ring.Points.Length )
					.Select( i => (A: ring.Points[i], B: ring.Points[(i + 1) % ring.Points.Length]) )
					.OrderBy( x => (points[x.A] - points[x.B]).Length )
					.First();

				if ( !moved.ContainsValue( from ) ) moved[from] = to;
			}

			if ( moved.Count == 0 ) return;

			foreach ( var ring in rings )
				ring.Points = WithoutSpikes( [.. ring.Points.Select( x => moved.TryGetValue( x, out var to ) ? to : x )] );

			rings.RemoveAll( x => x.Points.Length < 3 );
		}
	}

	static double Area( int[] ring, List<Vec> points )
	{
		var sum = new Vec();

		for ( var i = 1; i + 1 < ring.Length; i++ )
			sum += (points[ring[i]] - points[ring[0]]).Cross( points[ring[i + 1]] - points[ring[0]] );

		return sum.Length * 0.5;
	}

	/// <summary>
	/// A thin piece can lose all its width to welding and then runs out and straight back.
	/// </summary>
	static int[] WithoutSpikes( int[] ring )
	{
		var points = ring.ToList();
		var changed = true;

		while ( changed && points.Count >= 3 )
		{
			changed = false;

			for ( var i = 0; i < points.Count; i++ )
			{
				var before = points[(i + points.Count - 1) % points.Count];
				var after = points[(i + 1) % points.Count];

				if ( before == points[i] )
				{
					points.RemoveAt( i );
					changed = true;
					break;
				}

				if ( before != after ) continue;

				points.RemoveAt( i );
				points.RemoveAt( i % points.Count );
				changed = true;
				break;
			}
		}

		return [.. points];
	}

	/// <summary>
	/// Give every edge the points that lie on it. Two pieces next to each other can be split in
	/// different places, their shared edge then has corners on one side only.
	/// </summary>
	static void FillEdges( List<Ring> rings, List<Vec> points )
	{
		var order = Enumerable.Range( 0, points.Count ).OrderBy( i => points[i].X ).ToArray();
		var xs = order.Select( i => points[i].X ).ToArray();

		foreach ( var ring in rings )
		{
			var filled = new List<int>( ring.Points.Length );

			for ( var i = 0; i < ring.Points.Length; i++ )
			{
				var a = ring.Points[i];
				var b = ring.Points[(i + 1) % ring.Points.Length];

				filled.Add( a );

				var pa = points[a];
				var pb = points[b];
				var delta = pb - pa;
				var length = delta.Length;

				if ( length < 1e-9 ) continue;

				var direction = delta * (1.0 / length);
				var minX = Math.Min( pa.X, pb.X ) - OnEdgeDistance;
				var maxX = Math.Max( pa.X, pb.X ) + OnEdgeDistance;
				var first = Array.BinarySearch( xs, minX );

				if ( first < 0 ) first = ~first;

				var found = new List<(double T, int Point)>();

				for ( var k = first; k < order.Length && xs[k] <= maxX; k++ )
				{
					var c = order[k];

					if ( c == a || c == b ) continue;

					var offset = points[c] - pa;
					var t = offset.Dot( direction );

					if ( t <= OnEdgeDistance || t >= length - OnEdgeDistance ) continue;
					if ( (offset - direction * t).Length > OnEdgeDistance ) continue;

					found.Add( (t, c) );
				}

				found.Sort( ( x, y ) => x.T.CompareTo( y.T ) );

				foreach ( var (_, point) in found )
					filled.Add( point );
			}

			ring.Points = [.. filled];
		}
	}

	/// <summary>
	/// The rings in an order where each one shares an edge with one before it.
	/// </summary>
	static List<Ring> InOrder( List<Ring> rings )
	{
		var byEdge = new Dictionary<(int, int), List<int>>();

		for ( var r = 0; r < rings.Count; r++ )
		{
			var ring = rings[r].Points;

			for ( var i = 0; i < ring.Length; i++ )
			{
				var a = ring[i];
				var b = ring[(i + 1) % ring.Length];
				var key = a < b ? (a, b) : (b, a);

				if ( !byEdge.TryGetValue( key, out var list ) ) byEdge[key] = list = [];

				list.Add( r );
			}
		}

		var seen = new bool[rings.Count];
		var ordered = new List<Ring>( rings.Count );
		var queue = new Queue<int>();

		for ( var start = 0; start < rings.Count; start++ )
		{
			if ( seen[start] ) continue;

			seen[start] = true;
			queue.Enqueue( start );

			while ( queue.Count > 0 )
			{
				var current = queue.Dequeue();
				var ring = rings[current].Points;

				ordered.Add( rings[current] );

				for ( var i = 0; i < ring.Length; i++ )
				{
					var a = ring[i];
					var b = ring[(i + 1) % ring.Length];

					foreach ( var next in byEdge[a < b ? (a, b) : (b, a)] )
					{
						if ( seen[next] ) continue;

						seen[next] = true;
						queue.Enqueue( next );
					}
				}
			}
		}

		return ordered;
	}

	static PolygonMesh Build( PolygonMesh source, List<HalfEdgeMesh.FaceHandle> faces, List<Ring> rings, List<Vec> points, Material material )
	{
		var mesh = new PolygonMesh();

		mesh.SetTransform( source.Transform );

		var handles = new HalfEdgeMesh.VertexHandle?[points.Count];
		var newFaces = new List<HalfEdgeMesh.FaceHandle>();

		void Finish( HalfEdgeMesh.FaceHandle face, Ring ring )
		{
			if ( ring.Face >= 0 )
			{
				var from = faces[ring.Face];

				source.GetFaceTextureParameters( from, out var axisU, out var axisV, out var scale );

				mesh.SetFaceMaterial( face, source.GetFaceMaterial( from ) );
				mesh.SetFaceTextureParameters( face, axisU, axisV, scale );
				mesh.SetTextureOffset( face, source.GetTextureOffset( from ) );
			}
			else
			{
				mesh.SetFaceMaterial( face, material );
				newFaces.Add( face );
			}
		}

		// A face can only join at a corner that still has a free side, and what is free
		// changes with every face that joins. What did not fit gets another turn after the rest.
		var waiting = InOrder( rings );

		for ( var pass = 0; pass < JoinPasses && waiting.Count > 0; pass++ )
		{
			var later = new List<Ring>();

			foreach ( var ring in waiting )
			{
				var loop = new HalfEdgeMesh.VertexHandle[ring.Points.Length];

				for ( var i = 0; i < loop.Length; i++ )
					loop[i] = handles[ring.Points[i]] ??= mesh.AddVertex( points[ring.Points[i]].ToVector() );

				var face = mesh.AddFace( loop );

				if ( face.IsValid ) Finish( face, ring );
				else later.Add( ring );
			}

			if ( later.Count == waiting.Count ) break;

			waiting = later;
		}

		foreach ( var ring in waiting )
		{
			var face = mesh.AddFace( mesh.AddVertices( [.. ring.Points.Select( x => points[x].ToVector() )] ) );

			if ( face.IsValid ) Finish( face, ring );
		}

		mesh.SetSmoothingAngle( 40.0f );
		mesh.ComputeFaceTextureCoordinatesFromParameters();

		foreach ( var face in newFaces )
			mesh.TextureAlignToGrid( Transform.Zero, face );

		return mesh;
	}
}