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.
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;
}
}