namespace Sandbox;

/// <summary>
/// How the voxel shell is sealed before mesh extraction.
/// Splat scans produce thin shells with holes, the fill mode determines how those holes are patched to create watertight collision.
/// </summary>
public enum CollisionFillMode
{
	/// <summary>
	/// No fill, only the raw voxel shell is used. Small holes remain.
	/// </summary>
	None,

	/// <summary>
	/// For outdoor/terrain scenes. Walks each vertical (XY) column upward from the bottom and fills everything below the first solid voxel,
	/// creating a solid ground volume. Holes in the floor disappear because the entire volume beneath the surface becomes solid.
	/// </summary>
	FloorFill,

	/// <summary>
	/// For interior/room scenes. Dilates the solid grid to bridge wall gaps, then flood-fills inward from the bounding-box boundary.
	/// Everything reachable from outside becomes solid, leaving only the enclosed interior as empty space. Produces thick watertight walls with no holes.
	/// </summary>
	ExternalFill,
}

/// <summary>
/// Generates a voxelized collision mesh from Gaussian splat positions.
/// Attach alongside <see cref="GaussianSplatRenderer"/>, it will automatically
/// pick up the loaded splat data and rebuild collision when it changes.
/// </summary>
[Title( "Gaussian Splat Collider" )]
[Category( "Physics" )]
[Icon( "blur_on" )]
public sealed class GaussianSplatCollider : Collider
{
	/// <summary>
	/// Mesh shapes are always concave (non-convex triangle soup), they can only be static.
	/// </summary>
	public override bool IsConcave => true;

	float _voxelSize = 4.0f;

	/// <summary>
	/// Size of each voxel in world units (inches). Smaller = more accurate but more triangles.
	/// </summary>
	[Property, Group( "Voxel" ), Range( 1f, 32f )]
	public float VoxelSize
	{
		get => _voxelSize;
		set
		{
			if ( _voxelSize == value ) return;
			_voxelSize = value;
			RebuildMesh();
		}
	}

	/// <summary>
	/// Minimum cumulative opacity a voxel must contain to be considered solid.
	/// Each splat contributes its opacity (0..1) to the voxel it occupies.
	/// Higher values filter out stray low-opacity splats in empty space.
	/// </summary>
	[Property, Group( "Voxel" ), Range( 0.1f, 50f )]
	public float MinOccupancy
	{
		get => _minOccupancy;
		set
		{
			if ( _minOccupancy == value ) return;
			_minOccupancy = value;
			RebuildMesh();
		}
	}

	float _minOccupancy = 3.0f;

	/// <summary>
	/// Minimum per-splat opacity required for a splat to contribute to collision.
	/// Splats below this are completely ignored by the voxelizer.
	/// This is separate from the render MinOpacity, collision typically needs a much higher threshold to avoid invisible floaters.
	/// </summary>
	[Property, Group( "Voxel" ), Range( 0f, 1f )]
	public float MinOpacity
	{
		get => _minOpacity;
		set
		{
			if ( _minOpacity == value ) return;
			_minOpacity = value;
			RebuildMesh();
		}
	}

	float _minOpacity = 0.4f;

	/// <summary>
	/// When enabled, only the single largest connected cluster of solid voxels is kept, all smaller floating islands are removed.
	/// This is the most effective way to eliminate stray collision in empty space.
	/// </summary>
	[Property, Group( "Voxel" )]
	public bool KeepLargestClusterOnly
	{
		get => _keepLargestClusterOnly;
		set
		{
			if ( _keepLargestClusterOnly == value ) return;
			_keepLargestClusterOnly = value;
			RebuildMesh();
		}
	}

	bool _keepLargestClusterOnly = true;

	/// <summary>
	/// How the thin voxel shell is sealed before mesh extraction.
	/// FloorFill for outdoor/terrain, ExternalFill for rooms, None to skip.
	/// </summary>
	[Property, Group( "Fill" )]
	public CollisionFillMode FillMode
	{
		get => _fillMode;
		set
		{
			if ( _fillMode == value ) return;
			_fillMode = value;
			RebuildMesh();
		}
	}

	CollisionFillMode _fillMode = CollisionFillMode.FloorFill;

	/// <summary>
	/// Dilation distance in voxels used by ExternalFill to bridge small wall gaps
	/// before the outside flood-fill. Increase if walls have noisy holes the fill
	/// leaks through. Only used when FillMode is ExternalFill.
	/// </summary>
	[Property, Group( "Fill" ), Range( 0, 8 )]
	public int DilationRadius
	{
		get => _dilationRadius;
		set
		{
			if ( _dilationRadius == value ) return;
			_dilationRadius = value;
			RebuildMesh();
		}
	}

	int _dilationRadius = 2;

	/// <summary>
	/// For FloorFill: minimum number of solid neighbors (out of 8 in XY plane)
	/// an XY column must have to be eligible for floor-filling.
	/// Prevents filling columns in large open exterior areas. Only used when FillMode is FloorFill.
	/// </summary>
	[Property, Group( "Fill" ), Range( 0, 8 )]
	public int FloorFillNeighborThreshold
	{
		get => _floorFillNeighborThreshold;
		set
		{
			if ( _floorFillNeighborThreshold == value ) return;
			_floorFillNeighborThreshold = value;
			RebuildMesh();
		}
	}

	int _floorFillNeighborThreshold = 1;

	/// <summary>
	/// Number of Laplacian smoothing iterations to apply to the collision mesh.
	/// 0 = axis-aligned greedy mesh ("faces" mode, fastest, staircase edges).
	/// Higher values produce smoother slopes and contours at the cost of more triangles and slight collision inaccuracy. 8-16 is a good starting range.
	/// </summary>
	[Property, Group( "Mesh" ), Range( 0, 64 )]
	public int SmoothIterations
	{
		get => _smoothIterations;
		set
		{
			if ( _smoothIterations == value ) return;
			_smoothIterations = value;
			RebuildMesh();
		}
	}

	int _smoothIterations = 0;

	// Cached source data for rebuilding when settings change
	SceneGaussianSplatObject.SplatPosition[] _cachedPositions;
	float[] _cachedOpacities;

	// Cached mesh data
	Vector3[] _meshVertices;
	int[] _meshIndices;

	/// <summary>
	/// Triangle count of the generated collision mesh (read-only).
	/// </summary>
	[Property, Group( "Info" ), ReadOnly]
	public int TriangleCount => (_meshIndices?.Length ?? 0) / 3;

	protected override IEnumerable<PhysicsShape> CreatePhysicsShapes( PhysicsBody targetBody, Transform local )
	{
		if ( _meshVertices is null || _meshIndices is null || _meshVertices.Length == 0 )
			yield break;

		// Apply the GameObject's world scale to the mesh vertices so the
		// collision matches the rendered splat cloud at any transform scale.
		var scale = WorldScale;
		var needsScale = scale != Vector3.One;
		var vertices = _meshVertices;

		if ( needsScale )
		{
			vertices = new Vector3[_meshVertices.Length];
			for ( int i = 0; i < vertices.Length; i++ )
				vertices[i] = _meshVertices[i] * scale;
		}

		var shape = targetBody.AddMeshShape( vertices, _meshIndices );
		if ( shape.IsValid() )
			yield return shape;
	}

	/// <summary>
	/// Rebuild the collision mesh from the given splat positions and per-splat opacities.
	/// Called by <see cref="GaussianSplatRenderer"/> when splat data is loaded.
	/// </summary>
	public void BuildFromPositions( SceneGaussianSplatObject.SplatPosition[] positions, float[] opacities )
	{
		_cachedPositions = positions;
		_cachedOpacities = opacities;
		RebuildMesh();
	}

	/// <summary>
	/// Clear collision data.
	/// </summary>
	public void Clear()
	{
		_cachedPositions = null;
		_cachedOpacities = null;
		_meshVertices = null;
		_meshIndices = null;
		Rebuild();
	}

	void RebuildMesh()
	{
		if ( _cachedPositions is null || _cachedPositions.Length == 0 )
		{
			_meshVertices = null;
			_meshIndices = null;
			Rebuild();
			return;
		}

		BuildVoxelMesh( _cachedPositions, _cachedOpacities, _voxelSize, _minOccupancy, _minOpacity,
			_keepLargestClusterOnly, _fillMode, _dilationRadius, _floorFillNeighborThreshold,
			_smoothIterations, out _meshVertices, out _meshIndices );
		Rebuild();
	}

	/// <summary>
	/// Re-generate the mesh when voxel settings change (called via property setters → Rebuild)
	/// The base Collider.Rebuild calls RebuildImmediately which calls CreatePhysicsShapes,
	/// so we just need the mesh data to be current before that happens.
	/// </summary>
	protected override void RebuildImmediately()
	{
		// If we have cached source data and the mesh needs regenerating,
		// do it before the base class tries to create physics shapes.
		if ( _cachedPositions is not null && _meshVertices is null )
		{
			BuildVoxelMesh( _cachedPositions, _cachedOpacities, _voxelSize, _minOccupancy, _minOpacity,
				_keepLargestClusterOnly, _fillMode, _dilationRadius, _floorFillNeighborThreshold,
				_smoothIterations, out _meshVertices, out _meshIndices );
		}

		base.RebuildImmediately();
	}

	/// <summary>
	/// Full voxel collision pipeline following the splat-transform approach:
	/// 1. Accumulate opacity per voxel, threshold to solid set
	/// 2. Keep largest connected cluster (remove floaters)
	/// 3. Shell fill (FloorFill or ExternalFill) to seal holes
	/// 4. Greedy mesh extraction
	/// </summary>
	static void BuildVoxelMesh(
		ReadOnlySpan<SceneGaussianSplatObject.SplatPosition> positions,
		ReadOnlySpan<float> opacities,
		float voxelSize,
		float minOccupancy,
		float minOpacity,
		bool keepLargestClusterOnly,
		CollisionFillMode fillMode,
		int dilationRadius,
		int floorFillNeighborThreshold,
		int smoothIterations,
		out Vector3[] outVertices,
		out int[] outIndices )
	{
		float invSize = 1f / voxelSize;
		bool hasOpacities = opacities.Length == positions.Length;

		// Step 1: Accumulate opacity per voxel cell, skipping splats below the opacity floor
		var voxelDensity = new Dictionary<(int x, int y, int z), float>( positions.Length / 4 );

		for ( int i = 0; i < positions.Length; i++ )
		{
			float opacity = hasOpacities ? opacities[i] : 1f;

			if ( opacity < minOpacity )
				continue;

			var p = positions[i].Position;
			var key = (
				(int)MathF.Floor( p.x * invSize ),
				(int)MathF.Floor( p.y * invSize ),
				(int)MathF.Floor( p.z * invSize )
			);

			if ( voxelDensity.TryGetValue( key, out float existing ) )
				voxelDensity[key] = existing + opacity;
			else
				voxelDensity[key] = opacity;
		}

		// Step 2: Build solid voxel set from density threshold
		var solidVoxels = new HashSet<(int x, int y, int z)>( voxelDensity.Count );
		foreach ( var (key, density) in voxelDensity )
		{
			if ( density >= minOccupancy )
				solidVoxels.Add( key );
		}

		if ( solidVoxels.Count == 0 )
		{
			outVertices = Array.Empty<Vector3>();
			outIndices = Array.Empty<int>();
			return;
		}

		// Step 3: Keep only the largest connected cluster to remove floaters
		if ( keepLargestClusterOnly )
		{
			KeepLargestCluster( solidVoxels );
		}

		if ( solidVoxels.Count == 0 )
		{
			outVertices = Array.Empty<Vector3>();
			outIndices = Array.Empty<int>();
			return;
		}

		// Step 4: Shell fill, seal holes in the thin voxel surface
		switch ( fillMode )
		{
			case CollisionFillMode.FloorFill:
				ApplyFloorFill( solidVoxels, floorFillNeighborThreshold );
				break;
			case CollisionFillMode.ExternalFill:
				ApplyExternalFill( solidVoxels, dilationRadius );
				break;
		}

		// Step 5: Mesh extraction, greedy (faces) or smoothed
		if ( smoothIterations > 0 )
		{
			ExtractSmoothedMesh( solidVoxels, voxelSize, smoothIterations, out outVertices, out outIndices );
		}
		else
		{
			ExtractSurfaceMesh( solidVoxels, voxelSize, out outVertices, out outIndices );
		}
	}

	/// <summary>
	/// Floor fill for outdoor/terrain scenes.
	/// For each (X,Y) column that has at least one solid voxel and enough solid neighbors, fills everything below
	/// the highest solid voxel in that column. This turns a thin ground shell into a solid volume with no holes.
	/// </summary>
	static void ApplyFloorFill( HashSet<(int x, int y, int z)> solidVoxels, int neighborThreshold )
	{
		// Compute bounds
		int minX = int.MaxValue, minY = int.MaxValue, minZ = int.MaxValue;
		int maxX = int.MinValue, maxY = int.MinValue, maxZ = int.MinValue;

		foreach ( var (vx, vy, vz) in solidVoxels )
		{
			if ( vx < minX ) minX = vx; if ( vx > maxX ) maxX = vx;
			if ( vy < minY ) minY = vy; if ( vy > maxY ) maxY = vy;
			if ( vz < minZ ) minZ = vz; if ( vz > maxZ ) maxZ = vz;
		}

		// Build a set of XY columns that contain at least one solid voxel,
		// tracking the highest Z in each column
		var columnMaxZ = new Dictionary<(int x, int y), int>();
		foreach ( var (vx, vy, vz) in solidVoxels )
		{
			var col = (vx, vy);
			if ( columnMaxZ.TryGetValue( col, out int existing ) )
				columnMaxZ[col] = Math.Max( existing, vz );
			else
				columnMaxZ[col] = vz;
		}

		// 8-connected XY neighbor offsets for the proximity check
		ReadOnlySpan<(int dx, int dy)> xyDirs =
		[
			(1, 0), (-1, 0), (0, 1), (0, -1),
			(1, 1), (1, -1), (-1, 1), (-1, -1)
		];

		// For each column, fill from minZ up to the highest solid voxel
		// only if the column has enough neighboring columns with solid voxels
		foreach ( var (col, topZ) in columnMaxZ )
		{
			if ( neighborThreshold > 0 )
			{
				int solidNeighborCols = 0;
				foreach ( var (dx, dy) in xyDirs )
				{
					if ( columnMaxZ.ContainsKey( (col.x + dx, col.y + dy) ) )
						solidNeighborCols++;
				}

				if ( solidNeighborCols < neighborThreshold )
					continue;
			}

			// Fill the entire column below the top solid voxel
			for ( int z = minZ; z <= topZ; z++ )
			{
				solidVoxels.Add( (col.x, col.y, z) );
			}
		}
	}

	/// <summary>
	/// External fill for interior/room scenes. Dilates the solid set to bridge
	/// small wall gaps, then flood-fills from the bounding-box boundary inward.
	/// Everything reachable from outside becomes solid, leaving only the enclosed
	/// interior as empty space. Produces thick watertight walls.
	/// </summary>
	static void ApplyExternalFill( HashSet<(int x, int y, int z)> solidVoxels, int dilationRadius )
	{
		// Compute bounds (with padding for flood boundary)
		int minX = int.MaxValue, minY = int.MaxValue, minZ = int.MaxValue;
		int maxX = int.MinValue, maxY = int.MinValue, maxZ = int.MinValue;

		foreach ( var (vx, vy, vz) in solidVoxels )
		{
			if ( vx < minX ) minX = vx; if ( vx > maxX ) maxX = vx;
			if ( vy < minY ) minY = vy; if ( vy > maxY ) maxY = vy;
			if ( vz < minZ ) minZ = vz; if ( vz > maxZ ) maxZ = vz;
		}

		// Step A: Dilate — expand solid set by dilationRadius voxels to bridge gaps
		if ( dilationRadius > 0 )
		{
			var dilated = new HashSet<(int x, int y, int z)>( solidVoxels );

			// Iterative 6-connected dilation, one voxel per pass
			ReadOnlySpan<(int dx, int dy, int dz)> dirs =
			[
				(1, 0, 0), (-1, 0, 0),
				(0, 1, 0), (0, -1, 0),
				(0, 0, 1), (0, 0, -1)
			];

			for ( int pass = 0; pass < dilationRadius; pass++ )
			{
				var frontier = new List<(int x, int y, int z)>();
				foreach ( var cell in dilated )
				{
					foreach ( var (dx, dy, dz) in dirs )
					{
						var n = (cell.x + dx, cell.y + dy, cell.z + dz);
						if ( !dilated.Contains( n ) )
							frontier.Add( n );
					}
				}

				foreach ( var n in frontier )
					dilated.Add( n );
			}

			// Use the dilated set for the flood fill (not the original)
			// but we'll mark voxels in the original solidVoxels set
			var dilatedForFlood = dilated;

			// Step B: Flood-fill from the bounding-box boundary, treating dilated voxels as walls.
			// Everything NOT reachable = interior void. Mark exterior as solid.
			int pad = 2; // padding around bounds for the flood seed layer
			int bMinX = minX - pad, bMinY = minY - pad, bMinZ = minZ - pad;
			int bMaxX = maxX + pad, bMaxY = maxY + pad, bMaxZ = maxZ + pad;

			var exterior = new HashSet<(int x, int y, int z)>();
			var queue = new Queue<(int x, int y, int z)>();

			// Seed all boundary cells that aren't dilated-solid
			for ( int x = bMinX; x <= bMaxX; x++ )
			{
				for ( int y = bMinY; y <= bMaxY; y++ )
				{
					for ( int z = bMinZ; z <= bMaxZ; z++ )
					{
						bool isBoundary = x == bMinX || x == bMaxX
							|| y == bMinY || y == bMaxY
							|| z == bMinZ || z == bMaxZ;

						if ( isBoundary && !dilatedForFlood.Contains( (x, y, z) ) )
						{
							if ( exterior.Add( (x, y, z) ) )
								queue.Enqueue( (x, y, z) );
						}
					}
				}
			}

			// Flood-fill exterior
			while ( queue.Count > 0 )
			{
				var current = queue.Dequeue();

				foreach ( var (dx, dy, dz) in dirs )
				{
					var n = (current.x + dx, current.y + dy, current.z + dz);

					// Stay within padded bounds
					if ( n.Item1 < bMinX || n.Item1 > bMaxX
						|| n.Item2 < bMinY || n.Item2 > bMaxY
						|| n.Item3 < bMinZ || n.Item3 > bMaxZ )
						continue;

					// Can't flood through dilated walls
					if ( dilatedForFlood.Contains( n ) )
						continue;

					if ( exterior.Add( n ) )
						queue.Enqueue( n );
				}
			}

			// Step C: Everything inside the padded bounds that is NOT exterior
			// and NOT already solid → mark as solid. This fills the void outside
			// the room, making walls infinitely thick.
			for ( int x = minX; x <= maxX; x++ )
			{
				for ( int y = minY; y <= maxY; y++ )
				{
					for ( int z = minZ; z <= maxZ; z++ )
					{
						if ( !exterior.Contains( (x, y, z) ) )
							solidVoxels.Add( (x, y, z) );
					}
				}
			}
		}
	}

	/// <summary>
	/// 6-connected flood fill to find all connected components, then keep only
	/// the single largest one. This is the nuclear option for floating collision —
	/// it guarantees only the main structure has collision, regardless of how many
	/// stray splats exist elsewhere.
	/// </summary>
	static void KeepLargestCluster( HashSet<(int x, int y, int z)> solidVoxels )
	{
		var remaining = new HashSet<(int x, int y, int z)>( solidVoxels );
		var queue = new Queue<(int x, int y, int z)>();
		List<(int x, int y, int z)> largestCluster = null;
		var currentCluster = new List<(int x, int y, int z)>();

		ReadOnlySpan<(int dx, int dy, int dz)> dirs =
		[
			(1, 0, 0), (-1, 0, 0),
			(0, 1, 0), (0, -1, 0),
			(0, 0, 1), (0, 0, -1)
		];

		while ( remaining.Count > 0 )
		{
			var seed = remaining.First();
			currentCluster.Clear();
			queue.Clear();
			queue.Enqueue( seed );
			remaining.Remove( seed );

			while ( queue.Count > 0 )
			{
				var current = queue.Dequeue();
				currentCluster.Add( current );

				foreach ( var (dx, dy, dz) in dirs )
				{
					var neighbor = (current.x + dx, current.y + dy, current.z + dz);
					if ( remaining.Remove( neighbor ) )
					{
						queue.Enqueue( neighbor );
					}
				}
			}

			if ( largestCluster is null || currentCluster.Count > largestCluster.Count )
			{
				largestCluster = new List<(int x, int y, int z)>( currentCluster );
			}
		}

		solidVoxels.Clear();

		if ( largestCluster is not null )
		{
			foreach ( var v in largestCluster )
				solidVoxels.Add( v );
		}
	}

	/// <summary>
	/// Greedy meshing: for each face direction, slice through the voxel grid one layer
	/// at a time. Within each slice, merge coplanar exposed faces into maximal rectangles.
	/// This dramatically reduces triangle count for flat surfaces (floors, walls, ceilings).
	/// </summary>
	static void ExtractSurfaceMesh(
		HashSet<(int x, int y, int z)> solidVoxels,
		float voxelSize,
		out Vector3[] outVertices,
		out int[] outIndices )
	{
		// Compute axis-aligned bounding box of the solid voxels
		int minX = int.MaxValue, minY = int.MaxValue, minZ = int.MaxValue;
		int maxX = int.MinValue, maxY = int.MinValue, maxZ = int.MinValue;

		foreach ( var (vx, vy, vz) in solidVoxels )
		{
			if ( vx < minX ) minX = vx; if ( vx > maxX ) maxX = vx;
			if ( vy < minY ) minY = vy; if ( vy > maxY ) maxY = vy;
			if ( vz < minZ ) minZ = vz; if ( vz > maxZ ) maxZ = vz;
		}

		int sizeX = maxX - minX + 1;
		int sizeY = maxY - minY + 1;
		int sizeZ = maxZ - minZ + 1;

		var vertices = new List<Vector3>( solidVoxels.Count );
		var indices = new List<int>( solidVoxels.Count );

		// Process each of the 6 face directions.
		// For each direction, we iterate layer-by-layer along the normal axis
		// and greedily merge exposed faces within each 2D slice.
		//
		// Face directions: 0=+X, 1=-X, 2=+Y, 3=-Y, 4=+Z, 5=-Z
		// For each face we define:
		//   axis: the normal axis (0=X, 1=Y, 2=Z)
		//   sign: +1 or -1 (which side of the voxel)
		//   u,v: the two tangent axes for the 2D slice

		for ( int face = 0; face < 6; face++ )
		{
			int axis = face / 2;      // 0,0,1,1,2,2
			int sign = (face % 2 == 0) ? 1 : -1;

			// Tangent axes for the 2D slice
			int uAxis = (axis + 1) % 3;
			int vAxis = (axis + 2) % 3;

			// Bounds along each axis
			int axisMin = axis switch { 0 => minX, 1 => minY, _ => minZ };
			int axisMax = axis switch { 0 => maxX, 1 => maxY, _ => maxZ };
			int uMin = uAxis switch { 0 => minX, 1 => minY, _ => minZ };
			int uMax = uAxis switch { 0 => maxX, 1 => maxY, _ => maxZ };
			int vMin = vAxis switch { 0 => minX, 1 => minY, _ => minZ };
			int vMax = vAxis switch { 0 => maxX, 1 => maxY, _ => maxZ };
			int uSize = uMax - uMin + 1;
			int vSize = vMax - vMin + 1;

			// Mask of exposed faces in this slice
			var mask = new bool[uSize * vSize];

			for ( int layer = axisMin; layer <= axisMax; layer++ )
			{
				// Build the mask for this slice
				int maskIdx = 0;
				for ( int v = vMin; v <= vMax; v++ )
				{
					for ( int u = uMin; u <= uMax; u++ )
					{
						// Reconstruct 3D coords from axis/u/v
						var coords = new int[3];
						coords[axis] = layer;
						coords[uAxis] = u;
						coords[vAxis] = v;

						var voxel = (coords[0], coords[1], coords[2]);
						bool isSolid = solidVoxels.Contains( voxel );

						// Check neighbor in the normal direction
						var neighborCoords = new int[3];
						neighborCoords[axis] = layer + sign;
						neighborCoords[uAxis] = u;
						neighborCoords[vAxis] = v;
						var neighbor = (neighborCoords[0], neighborCoords[1], neighborCoords[2]);

						mask[maskIdx++] = isSolid && !solidVoxels.Contains( neighbor );
					}
				}

				// Greedy merge: scan the 2D mask and find maximal rectangles
				for ( int v = 0; v < vSize; v++ )
				{
					for ( int u = 0; u < uSize; )
					{
						int idx = v * uSize + u;
						if ( !mask[idx] )
						{
							u++;
							continue;
						}

						// Expand width (along u)
						int w = 1;
						while ( u + w < uSize && mask[v * uSize + u + w] )
							w++;

						// Expand height (along v) as far as the full width remains exposed
						int h = 1;
						bool canExpand = true;
						while ( v + h < vSize && canExpand )
						{
							for ( int du = 0; du < w; du++ )
							{
								if ( !mask[(v + h) * uSize + u + du] )
								{
									canExpand = false;
									break;
								}
							}
							if ( canExpand ) h++;
						}

						// Clear the merged region from the mask
						for ( int dv = 0; dv < h; dv++ )
							for ( int du = 0; du < w; du++ )
								mask[(v + dv) * uSize + u + du] = false;

						// Emit a quad for this merged rectangle
						// The quad corners in 3D, placed on the correct face of the voxel
						float faceOffset = sign > 0 ? 1f : 0f;
						float layerF = (layer + faceOffset) * voxelSize;
						float u0 = (uMin + u) * voxelSize;
						float u1 = (uMin + u + w) * voxelSize;
						float v0 = (vMin + v) * voxelSize;
						float v1 = (vMin + v + h) * voxelSize;

						var c0 = new float[3];
						var c1 = new float[3];
						var c2 = new float[3];
						var c3 = new float[3];

						// All four corners share the same position along the normal axis
						c0[axis] = layerF; c1[axis] = layerF; c2[axis] = layerF; c3[axis] = layerF;

						// Winding order depends on face sign to keep outward-facing normals
						if ( sign > 0 )
						{
							c0[uAxis] = u0; c0[vAxis] = v0;
							c1[uAxis] = u1; c1[vAxis] = v0;
							c2[uAxis] = u1; c2[vAxis] = v1;
							c3[uAxis] = u0; c3[vAxis] = v1;
						}
						else
						{
							c0[uAxis] = u0; c0[vAxis] = v0;
							c1[uAxis] = u0; c1[vAxis] = v1;
							c2[uAxis] = u1; c2[vAxis] = v1;
							c3[uAxis] = u1; c3[vAxis] = v0;
						}

						int baseVertex = vertices.Count;
						vertices.Add( new Vector3( c0[0], c0[1], c0[2] ) );
						vertices.Add( new Vector3( c1[0], c1[1], c1[2] ) );
						vertices.Add( new Vector3( c2[0], c2[1], c2[2] ) );
						vertices.Add( new Vector3( c3[0], c3[1], c3[2] ) );

						indices.Add( baseVertex + 0 );
						indices.Add( baseVertex + 1 );
						indices.Add( baseVertex + 2 );
						indices.Add( baseVertex + 0 );
						indices.Add( baseVertex + 2 );
						indices.Add( baseVertex + 3 );

						u += w;
					}
				}
			}
		}

		outVertices = vertices.ToArray();
		outIndices = indices.ToArray();
	}

	/// <summary>
	/// Extract a surface mesh with shared vertices at voxel corners, then apply
	/// Laplacian smoothing to produce natural contours instead of staircase edges.
	/// Suitable for character collision on slopes and curved surfaces.
	/// </summary>
	static void ExtractSmoothedMesh(
		HashSet<(int x, int y, int z)> solidVoxels,
		float voxelSize,
		int iterations,
		out Vector3[] outVertices,
		out int[] outIndices )
	{
		// Face definitions: neighbor offset + 4 corner offsets (CCW winding for outward normal)
		ReadOnlySpan<int> faceNeighbor =
		[
			1, 0, 0,   // +X
			-1, 0, 0,  // -X
			0, 1, 0,   // +Y
			0, -1, 0,  // -Y
			0, 0, 1,   // +Z
			0, 0, -1,  // -Z
		];

		// 4 corner offsets per face (integer coords relative to voxel origin)
		ReadOnlySpan<int> faceCorners =
		[
			// +X face
			1,0,0,  1,1,0,  1,1,1,  1,0,1,
			// -X face
			0,1,0,  0,0,0,  0,0,1,  0,1,1,
			// +Y face
			1,1,0,  0,1,0,  0,1,1,  1,1,1,
			// -Y face
			0,0,0,  1,0,0,  1,0,1,  0,0,1,
			// +Z face
			0,0,1,  1,0,1,  1,1,1,  0,1,1,
			// -Z face
			1,0,0,  0,0,0,  0,1,0,  1,1,0,
		];

		// Step 1: Extract faces with welded vertices at voxel corners.
		// Vertex key is the integer corner position, shared between adjacent faces.
		var vertexMap = new Dictionary<(int x, int y, int z), int>();
		var vertexList = new List<Vector3>();
		var indexList = new List<int>();

		foreach ( var (vx, vy, vz) in solidVoxels )
		{
			for ( int face = 0; face < 6; face++ )
			{
				int ni = face * 3;
				var neighbor = (vx + faceNeighbor[ni], vy + faceNeighbor[ni + 1], vz + faceNeighbor[ni + 2]);
				if ( solidVoxels.Contains( neighbor ) )
					continue;

				// Get or create the 4 corner vertices for this face
				int ci = face * 12;
				var vi = new int[4];

				for ( int c = 0; c < 4; c++ )
				{
					int ci2 = ci + c * 3;
					var cornerKey = (vx + faceCorners[ci2], vy + faceCorners[ci2 + 1], vz + faceCorners[ci2 + 2]);

					if ( !vertexMap.TryGetValue( cornerKey, out int idx ) )
					{
						idx = vertexList.Count;
						vertexMap[cornerKey] = idx;
						vertexList.Add( new Vector3( cornerKey.Item1 * voxelSize, cornerKey.Item2 * voxelSize, cornerKey.Item3 * voxelSize ) );
					}

					vi[c] = idx;
				}

				// Two triangles per quad
				indexList.Add( vi[0] );
				indexList.Add( vi[1] );
				indexList.Add( vi[2] );
				indexList.Add( vi[0] );
				indexList.Add( vi[2] );
				indexList.Add( vi[3] );
			}
		}

		if ( vertexList.Count == 0 )
		{
			outVertices = Array.Empty<Vector3>();
			outIndices = Array.Empty<int>();
			return;
		}

		// Step 2: Build adjacency, for each vertex, which other vertices are connected by an edge
		var adjacency = new HashSet<int>[vertexList.Count];
		for ( int i = 0; i < adjacency.Length; i++ )
			adjacency[i] = new HashSet<int>();

		for ( int i = 0; i < indexList.Count; i += 3 )
		{
			int a = indexList[i], b = indexList[i + 1], c = indexList[i + 2];
			adjacency[a].Add( b ); adjacency[a].Add( c );
			adjacency[b].Add( a ); adjacency[b].Add( c );
			adjacency[c].Add( a ); adjacency[c].Add( b );
		}

		// Step 3: Laplacian smoothing, iteratively move each vertex toward
		// the average of its edge-connected neighbors. Uses Taubin smoothing
		// (alternating positive/negative lambda) to smooth without shrinkage.
		var positions = vertexList.ToArray();
		var temp = new Vector3[positions.Length];
		float lambda = 0.5f;
		float mu = -0.53f; // slightly larger magnitude than lambda to prevent shrinkage

		for ( int iter = 0; iter < iterations; iter++ )
		{
			float factor = (iter % 2 == 0) ? lambda : mu;

			for ( int v = 0; v < positions.Length; v++ )
			{
				var neighbors = adjacency[v];
				if ( neighbors.Count == 0 )
				{
					temp[v] = positions[v];
					continue;
				}

				var avg = Vector3.Zero;
				foreach ( int n in neighbors )
					avg += positions[n];
				avg /= neighbors.Count;

				temp[v] = positions[v] + factor * (avg - positions[v]);
			}

			// Swap
			(positions, temp) = (temp, positions);
		}

		outVertices = positions;
		outIndices = indexList.ToArray();
	}

	protected override void DrawGizmos()
	{
		if ( !Gizmo.IsSelected && !Gizmo.IsHovered )
			return;

		if ( _meshVertices is null || _meshIndices is null || _meshIndices.Length < 3 )
			return;

		Gizmo.Draw.Color = Gizmo.Colors.Green.WithAlpha( Gizmo.IsSelected ? 0.5f : 0.1f );

		int totalTriangles = _meshIndices.Length / 3;

		// Batch triangle drawing to stay within the gizmo system's internal
		// vertex buffer limits. A single huge LineTriangles call gets silently
		// dropped when the mesh is large (common with dense SOG/PLY files).
		const int batchSize = 16384;

		for ( int start = 0; start < totalTriangles; start += batchSize )
		{
			int count = Math.Min( batchSize, totalTriangles - start );
			var triangles = new Triangle[count];

			for ( int i = 0; i < count; i++ )
			{
				int bi = (start + i) * 3;
				triangles[i] = new Triangle(
					_meshVertices[_meshIndices[bi]],
					_meshVertices[_meshIndices[bi + 1]],
					_meshVertices[_meshIndices[bi + 2]] );
			}

			Gizmo.Draw.LineTriangles( triangles );
		}
	}
}