Code/Building/CoastScanner.cs

CoastScanner is a static utility that takes a height field scan and extracts coastlines as closed, simplified 2D loops in world space. It marks the scan edges as sea, floods reachable sea cells, builds a contour via marching squares, follows edges into closed loops, simplifies them with Douglas-Peucker, filters by orientation and minimum area, and returns arrays of Vector2 points.

File Access
using System;
using System.Collections.Generic;
using Sandbox;

namespace RedSnail.WaterTool2;

/// <summary>
/// Turns a scan of the ground into coastlines: closed loops round every piece of land the open sea reaches. The scan is
/// a grid of heights relative to the coast's height (land where positive); the open sea is whatever below that the
/// scan's edge connects to, so a basin inland, however low, is land and the sea never floods it.
/// </summary>
internal static class CoastScanner
{
	/// <summary>
	/// The coastlines in a scan: <paramref name="_Field"/> holds, at every grid point (row after row, <paramref
	/// name="_Width"/> to a row), how far the ground is above the coast's height. The scan's edge is taken as open sea.
	/// Each loop is in the world (in plan), simplified to within a quarter of a grid step; loops round less land than
	/// <paramref name="_MinArea"/> are left out. The field is changed in place.
	/// </summary>
	public static List<Vector2[]> Extract(float[] _Field, int _Width, int _Height, Vector2 _Origin, float _Step, float _MinArea)
	{
		var loops = new List<Vector2[]>();

		if (_Width < 3 || _Height < 3)
			return loops;

		// The scan's edge is the open sea, whatever is there.
		for (int i = 0; i < _Width; i++)
		{
			_Field[i] = MathF.Min(_Field[i], -1.0f);
			_Field[(_Height - 1) * _Width + i] = MathF.Min(_Field[(_Height - 1) * _Width + i], -1.0f);
		}

		for (int j = 0; j < _Height; j++)
		{
			_Field[j * _Width] = MathF.Min(_Field[j * _Width], -1.0f);
			_Field[j * _Width + _Width - 1] = MathF.Min(_Field[j * _Width + _Width - 1], -1.0f);
		}

		FillBasins(_Field, _Width, _Height);

		var next = Contour(_Field, _Width, _Height);
		var visited = new HashSet<int>();

		foreach (int start in next.Keys)
		{
			if (visited.Contains(start))
				continue;

			var loop = new List<Vector2>();
			int edge = start;
			bool closed = false;

			while (visited.Add(edge))
			{
				loop.Add(Crossing(_Field, _Width, _Origin, _Step, edge));

				if (!next.TryGetValue(edge, out int following))
					break;

				edge = following;
				closed = edge == start;
			}

			if (!closed || loop.Count < 3)
				continue;

			var simplified = SimplifyClosed(loop, _Step * 0.25f);

			// Land is on the right of the way round, so a loop round land runs clockwise; one the other way round would be
			// sea shut in by land, which the fill has already made land.
			float area = SignedArea(simplified);

			if (simplified.Count < 3 || area > 0.0f || -area < _MinArea)
				continue;

			loops.Add(simplified.ToArray());
		}

		return loops;
	}



	/// <summary>
	/// Everything below the coast's height that the open sea doesn't reach (through the grid, edge to edge) becomes land.
	/// </summary>
	private static void FillBasins(float[] _Field, int _Width, int _Height)
	{
		var sea = new bool[_Field.Length];
		var queue = new Queue<int>();

		for (int i = 0; i < _Field.Length; i++)
		{
			int x = i % _Width;
			int y = i / _Width;

			if ((x == 0 || y == 0 || x == _Width - 1 || y == _Height - 1) && _Field[i] <= 0.0f)
			{
				sea[i] = true;
				queue.Enqueue(i);
			}
		}

		while (queue.Count > 0)
		{
			int i = queue.Dequeue();
			int x = i % _Width;
			int y = i / _Width;

			Visit(x - 1, y);
			Visit(x + 1, y);
			Visit(x, y - 1);
			Visit(x, y + 1);
		}

		for (int i = 0; i < _Field.Length; i++)
		{
			if (!sea[i] && _Field[i] <= 0.0f)
				_Field[i] = 0.01f;
		}

		void Visit(int _X, int _Y)
		{
			if (_X < 0 || _Y < 0 || _X >= _Width || _Y >= _Height)
				return;

			int index = _Y * _Width + _X;

			if (sea[index] || _Field[index] > 0.0f)
				return;

			sea[index] = true;
			queue.Enqueue(index);
		}
	}



	/// <summary>
	/// Marching squares: every grid cell the coast crosses gives a piece of it, from edge to edge, with the land on its
	/// right. Edges are numbered so a piece leaving one cell through an edge is the piece that enters the next cell
	/// through it: the result maps every crossed edge to the next one along the coast.
	/// </summary>
	private static Dictionary<int, int> Contour(float[] _Field, int _Width, int _Height)
	{
		var next = new Dictionary<int, int>();
		var edges = new int[4];
		var corners = new bool[4];
		var crossEdges = new int[4];
		var crossIntoLand = new bool[4];

		for (int y = 0; y < _Height - 1; y++)
		{
			for (int x = 0; x < _Width - 1; x++)
			{
				float f00 = _Field[y * _Width + x];
				float f10 = _Field[y * _Width + x + 1];
				float f11 = _Field[(y + 1) * _Width + x + 1];
				float f01 = _Field[(y + 1) * _Width + x];

				// Round the cell counter-clockwise: bottom, right, top, left.
				corners[0] = f00 > 0.0f;
				corners[1] = f10 > 0.0f;
				corners[2] = f11 > 0.0f;
				corners[3] = f01 > 0.0f;

				if (corners[0] == corners[1] && corners[1] == corners[2] && corners[2] == corners[3])
					continue;

				edges[0] = HorizontalEdge(x, y, _Width);
				edges[1] = VerticalEdge(x + 1, y, _Width);
				edges[2] = HorizontalEdge(x, y + 1, _Width);
				edges[3] = VerticalEdge(x, y, _Width);

				int crossings = 0;

				for (int k = 0; k < 4; k++)
				{
					bool from = corners[k];
					bool to = corners[(k + 1) % 4];

					if (from == to)
						continue;

					crossEdges[crossings] = edges[k];
					crossIntoLand[crossings] = to;
					crossings++;
				}

				// A piece runs from where the way round goes into the land to where it comes out, which keeps the land on
				// its right. A saddle has two: they cut off the land corners when the middle is sea, the sea corners when
				// it is land.
				bool middleLand = f00 + f10 + f11 + f01 > 0.0f;

				for (int k = 0; k < crossings; k++)
				{
					if (!crossIntoLand[k])
						continue;

					int partner = crossings == 4 && middleLand ? (k + crossings - 1) % crossings : (k + 1) % crossings;

					next[crossEdges[k]] = crossEdges[partner];
				}
			}
		}

		return next;
	}



	private static int HorizontalEdge(int _X, int _Y, int _Width)
	{
		return (_Y * _Width + _X) * 2;
	}



	private static int VerticalEdge(int _X, int _Y, int _Width)
	{
		return (_Y * _Width + _X) * 2 + 1;
	}



	/// <summary>
	/// Where the coast crosses a grid edge, between its two points' heights.
	/// </summary>
	private static Vector2 Crossing(float[] _Field, int _Width, Vector2 _Origin, float _Step, int _Edge)
	{
		int cell = _Edge / 2;
		int x = cell % _Width;
		int y = cell / _Width;
		bool vertical = (_Edge & 1) == 1;

		float a = _Field[y * _Width + x];
		float b = vertical ? _Field[(y + 1) * _Width + x] : _Field[y * _Width + x + 1];
		float t = Math.Clamp(a / (a - b), 0.01f, 0.99f);

		return vertical
			? new Vector2(_Origin.x + x * _Step, _Origin.y + (y + t) * _Step)
			: new Vector2(_Origin.x + (x + t) * _Step, _Origin.y + y * _Step);
	}



	/// <summary>
	/// Douglas-Peucker on a closed loop: split at the point furthest from the first, each half simplified on its own.
	/// </summary>
	private static List<Vector2> SimplifyClosed(List<Vector2> _Loop, float _Tolerance)
	{
		if (_Loop.Count < 8)
			return _Loop;

		int far = 0;
		float farthest = -1.0f;

		for (int i = 1; i < _Loop.Count; i++)
		{
			float distance = _Loop[i].Distance(_Loop[0]);

			if (distance <= farthest)
				continue;

			farthest = distance;
			far = i;
		}

		var keep = new bool[_Loop.Count + 1];
		keep[0] = true;
		keep[far] = true;
		keep[_Loop.Count] = true;

		Simplify(_Loop, 0, far, _Tolerance, keep);
		Simplify(_Loop, far, _Loop.Count, _Tolerance, keep);

		var result = new List<Vector2>();

		for (int i = 0; i < _Loop.Count; i++)
		{
			if (keep[i])
				result.Add(_Loop[i]);
		}

		return result;
	}



	/// <summary>
	/// Marks which points between two kept ones stay (index <c>Count</c> is the first point again), without recursion.
	/// </summary>
	private static void Simplify(List<Vector2> _Loop, int _From, int _To, float _Tolerance, bool[] _Keep)
	{
		var stack = new Stack<(int From, int To)>();
		stack.Push((_From, _To));

		while (stack.Count > 0)
		{
			var (from, to) = stack.Pop();

			if (to - from < 2)
				continue;

			Vector2 a = _Loop[from];
			Vector2 b = _Loop[to % _Loop.Count];
			int worst = -1;
			float worstDistance = _Tolerance;

			for (int i = from + 1; i < to; i++)
			{
				float distance = DistanceToSegment(_Loop[i], a, b);

				if (distance <= worstDistance)
					continue;

				worst = i;
				worstDistance = distance;
			}

			if (worst < 0)
				continue;

			_Keep[worst] = true;
			stack.Push((from, worst));
			stack.Push((worst, to));
		}
	}



	private static float DistanceToSegment(Vector2 _Point, Vector2 _A, Vector2 _B)
	{
		Vector2 span = _B - _A;
		float lengthSquared = span.LengthSquared;
		float t = lengthSquared > 0.0f ? Math.Clamp(Vector2.Dot(_Point - _A, span) / lengthSquared, 0.0f, 1.0f) : 0.0f;

		return (_A + span * t).Distance(_Point);
	}



	private static float SignedArea(List<Vector2> _Loop)
	{
		float area = 0.0f;

		for (int i = 0; i < _Loop.Count; i++)
		{
			Vector2 a = _Loop[i];
			Vector2 b = _Loop[(i + 1) % _Loop.Count];

			area += a.x * b.y - b.x * a.y;
		}

		return area * 0.5f;
	}
}