using System;
using System.Collections.Generic;
using Sandbox;

namespace RedSnail.WaterTool2;

/// <summary>
/// A point of an outline a surface is meshed inside. Points that belong to a river's row carry the exact state of the
/// water on that row, and the edge between two points of the same row is never split: the surface meets the river on
/// the river's own vertices, with the river's own data on them, so nothing can open up between the two.
/// </summary>
internal struct OutlineVertex
{
	public Vector3 Position;

	/// <summary>
	/// The shared row this point belongs to, or -1 for a free point (a shore, a bank, a fillet).
	/// </summary>
	public int Row;

	public Vector2 Velocity;
	public float Waves;
	public float Foam;

	/// <summary>
	/// A point where the surface meets the sea's own (the edge of a river mouth's water, against the coast zone): it is
	/// the sea there exactly, and the edge between two such points is never split either, so it stays on the grid the sea
	/// is drawn on.
	/// </summary>
	public bool SeaEdge;

	public readonly bool IsFixed => Row >= 0;



	public static OutlineVertex Free(Vector3 _Position)
	{
		return new OutlineVertex { Position = _Position, Row = -1 };
	}



	public static OutlineVertex OnSeaEdge(Vector3 _Position)
	{
		return new OutlineVertex { Position = _Position, Row = -1, SeaEdge = true };
	}



	public static OutlineVertex OnRow(RiverRow _Row, int _Index, int _RowId)
	{
		return new OutlineVertex
		{
			Position = _Row.Points[_Index],
			Row = _RowId,
			Velocity = _Row.Velocities[_Index],
			Waves = _Row.Waves[_Index],
			Foam = _Row.Foams[_Index]
		};
	}
}



/// <summary>
/// Meshes the inside of a closed outline for the water: the outline's own points plus a grid of points inside it, joined
/// into a constrained Delaunay triangulation. Waves need vertices all over a surface, not just round its edge, and the
/// triangulation keeps every edge of the outline as it is, which is what lets a lake or a confluence share its edge
/// with the rows of the rivers running into it.
/// </summary>
internal static class SurfaceMesher
{
	/// <summary>
	/// Interior points closer than this many spacings to the outline are left out: they would only make slivers.
	/// </summary>
	private const float EdgeClearance = 0.5f;

	/// <summary>
	/// The interior grid is offset from the world's axes by these fractions of a spacing, so that outlines drawn on
	/// round numbers never put a grid point exactly on an edge.
	/// </summary>
	private const double GridOffsetX = 0.31;
	private const double GridOffsetY = 0.57;



	/// <summary>
	/// Splits every edge of an outline that is not part of a river's row into pieces no longer than the spacing, so
	/// the surface has vertices all along its shore. The new points are free; their heights are interpolated.
	/// </summary>
	public static List<OutlineVertex> Subdivide(IReadOnlyList<OutlineVertex> _Outline, float _Spacing)
	{
		var result = new List<OutlineVertex>(_Outline.Count * 2);
		float spacing = MathF.Max(1.0f, _Spacing);

		for (int i = 0; i < _Outline.Count; i++)
		{
			var from = _Outline[i];
			var to = _Outline[(i + 1) % _Outline.Count];

			result.Add(from);

			// An edge between two points of the same row is the river's own, one along the sea's edge the sea's: neither
			// is ever split.
			if ((from.IsFixed && to.IsFixed && from.Row == to.Row) || (from.SeaEdge && to.SeaEdge))
				continue;

			float length = new Vector2(to.Position.x - from.Position.x, to.Position.y - from.Position.y).Length;
			int pieces = (int)MathF.Ceiling(length / spacing);

			for (int k = 1; k < pieces; k++)
				result.Add(OutlineVertex.Free(Vector3.Lerp(from.Position, to.Position, (float)k / pieces)));
		}

		return result;
	}



	/// <summary>
	/// Removes points that sit on the one before them, which would only make edges of no length. Points of a river's
	/// row (or of the sea's edge) are always kept; the free point next to one goes instead.
	/// </summary>
	public static List<OutlineVertex> RemoveDuplicates(IReadOnlyList<OutlineVertex> _Outline, float _Distance = 0.5f)
	{
		var result = new List<OutlineVertex>(_Outline.Count);

		foreach (var vertex in _Outline)
		{
			if (result.Count > 0 && Flat(result[^1].Position).Distance(Flat(vertex.Position)) < _Distance)
			{
				if (IsKept(result[^1]) || !IsKept(vertex))
					continue;

				result[^1] = vertex;
				continue;
			}

			result.Add(vertex);
		}

		// The loop closes back on its first point.
		while (result.Count > 1 && Flat(result[^1].Position).Distance(Flat(result[0].Position)) < _Distance)
		{
			if (IsKept(result[^1]) && !IsKept(result[0]))
				result.RemoveAt(0);
			else
				result.RemoveAt(result.Count - 1);
		}

		return result;
	}



	private static bool IsKept(in OutlineVertex _Vertex)
	{
		return _Vertex.IsFixed || _Vertex.SeaEdge;
	}



	/// <summary>
	/// Triangulates the inside of a closed outline (either winding) together with a grid of interior points every
	/// <paramref name="_Spacing"/>. The result lists the outline's points first, in their order, then the interior
	/// points; every triangle winds counter-clockwise seen from above. False when the outline is not a simple polygon.
	/// </summary>
	public static bool Triangulate(IReadOnlyList<Vector2> _Outline, float _Spacing, out List<Vector2> _Points, out List<int> _Triangles)
	{
		_Points = null;
		_Triangles = null;

		int count = _Outline.Count;

		if (count < 3)
			return false;

		// Everything is worked out relative to the outline's middle, in doubles: a lake can be tens of thousands of
		// units across and its predicates must still agree with each other.
		var bounds = new BBox(new Vector3(_Outline[0].x, _Outline[0].y, 0.0f), new Vector3(_Outline[0].x, _Outline[0].y, 0.0f));

		foreach (var point in _Outline)
			bounds = bounds.AddPoint(new Vector3(point.x, point.y, 0.0f));

		double originX = Math.Round(bounds.Center.x);
		double originY = Math.Round(bounds.Center.y);

		double area = 0.0;

		for (int i = 0; i < count; i++)
		{
			var a = _Outline[i];
			var b = _Outline[(i + 1) % count];

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

		if (Math.Abs(area) < 1e-6)
			return false;

		// The triangulation wants the outline counter-clockwise; a clockwise one is walked backwards and mapped back.
		bool reversed = area < 0.0;
		var mesh = new Triangulation(Math.Max(bounds.Size.x, bounds.Size.y));

		for (int i = 0; i < count; i++)
		{
			var point = _Outline[reversed ? count - 1 - i : i];
			mesh.AddPoint(point.x - originX, point.y - originY);
		}

		if (!mesh.EarClip(count))
			return false;

		mesh.BuildAdjacency(count);
		mesh.MakeDelaunay();

		InsertGrid(mesh, _Outline, reversed, originX, originY, _Spacing);

		// Back to the caller's order: outline points first, as given, then the interior points.
		_Points = new List<Vector2>(mesh.PointCount);

		for (int i = 0; i < count; i++)
			_Points.Add(_Outline[i]);

		for (int i = count; i < mesh.PointCount; i++)
			_Points.Add(new Vector2((float)(mesh.X[i] + originX), (float)(mesh.Y[i] + originY)));

		_Triangles = new List<int>(mesh.TriangleCount * 3);

		for (int t = 0; t < mesh.TriangleCount; t++)
		{
			for (int k = 0; k < 3; k++)
			{
				int index = mesh.V[t * 3 + k];

				if (reversed && index < count)
					index = count - 1 - index;

				_Triangles.Add(index);
			}
		}

		return true;
	}



	/// <summary>
	/// Adds the interior grid: every point of it inside the outline and clear of its edges, row by row and back and
	/// forth along the rows, so each one is found a short walk from the last.
	/// </summary>
	private static void InsertGrid(Triangulation _Mesh, IReadOnlyList<Vector2> _Outline, bool _Reversed, double _OriginX, double _OriginY, float _Spacing)
	{
		double spacing = Math.Max(1.0, _Spacing);
		int count = _Outline.Count;

		// The outline relative to the origin, in the triangulation's own (counter-clockwise) order.
		var xs = new double[count];
		var ys = new double[count];
		double minX = double.MaxValue, minY = double.MaxValue, maxX = double.MinValue, maxY = double.MinValue;

		for (int i = 0; i < count; i++)
		{
			xs[i] = _Mesh.X[i];
			ys[i] = _Mesh.Y[i];

			minX = Math.Min(minX, xs[i]);
			minY = Math.Min(minY, ys[i]);
			maxX = Math.Max(maxX, xs[i]);
			maxY = Math.Max(maxY, ys[i]);
		}

		// The grid is anchored in the world, not in the outline, so it does not swim as the outline is edited.
		double offsetX = Mod(GridOffsetX * spacing - _OriginX, spacing);
		double offsetY = Mod(GridOffsetY * spacing - _OriginY, spacing);
		double clearance = EdgeClearance * spacing;

		var edges = new EdgeGrid(xs, ys, spacing, clearance);
		var crossings = new List<double>();
		var row = new List<double>();
		var hints = new Dictionary<long, int>();

		int hint = 0;
		int firstRow = (int)Math.Ceiling((minY - offsetY) / spacing);
		int lastRow = (int)Math.Floor((maxY - offsetY) / spacing);

		for (int j = firstRow; j <= lastRow; j++)
		{
			double y = j * spacing + offsetY;

			// Where the row crosses the outline: inside between each pair of crossings.
			crossings.Clear();

			for (int i = 0; i < count; i++)
			{
				int next = (i + 1) % count;
				double y0 = ys[i];
				double y1 = ys[next];

				if ((y0 > y) == (y1 > y))
					continue;

				crossings.Add(xs[i] + (y - y0) * (xs[next] - xs[i]) / (y1 - y0));
			}

			crossings.Sort();
			row.Clear();

			for (int c = 0; c + 1 < crossings.Count; c += 2)
			{
				int from = (int)Math.Ceiling((crossings[c] - offsetX) / spacing);
				int to = (int)Math.Floor((crossings[c + 1] - offsetX) / spacing);

				for (int i = from; i <= to; i++)
				{
					double x = i * spacing + offsetX;

					if (!edges.IsNear(x, y))
						row.Add(x);
				}
			}

			// Back and forth, so the walk to the first point of a row starts right next to it.
			if (((j - firstRow) & 1) == 1)
				row.Reverse();

			foreach (double x in row)
			{
				long bucket = BucketKey(x, y, spacing * 2.0);
				int triangle = _Mesh.Locate(x, y, hint);

				if (triangle < 0 && hints.TryGetValue(bucket, out int bucketHint))
					triangle = _Mesh.Locate(x, y, bucketHint);

				if (triangle < 0)
					triangle = FindNearHint(_Mesh, hints, x, y, spacing * 2.0);

				if (triangle < 0)
					triangle = _Mesh.Scan(x, y);

				if (triangle < 0)
					continue;

				int inserted = _Mesh.InsertAt(triangle, x, y, spacing);

				if (inserted < 0)
					continue;

				hint = inserted;
				hints[bucket] = inserted;
			}
		}
	}



	/// <summary>
	/// A walk starting from the hint of any bucket around the point.
	/// </summary>
	private static int FindNearHint(Triangulation _Mesh, Dictionary<long, int> _Hints, double _X, double _Y, double _BucketSize)
	{
		long bx = (long)Math.Floor(_X / _BucketSize);
		long by = (long)Math.Floor(_Y / _BucketSize);

		for (int ring = 1; ring <= 2; ring++)
		{
			for (long dy = -ring; dy <= ring; dy++)
			{
				for (long dx = -ring; dx <= ring; dx++)
				{
					if (Math.Max(Math.Abs(dx), Math.Abs(dy)) != ring)
						continue;

					if (!_Hints.TryGetValue(Key(bx + dx, by + dy), out int start))
						continue;

					int triangle = _Mesh.Locate(_X, _Y, start);

					if (triangle >= 0)
						return triangle;
				}
			}
		}

		return -1;
	}



	private static long BucketKey(double _X, double _Y, double _BucketSize)
	{
		return Key((long)Math.Floor(_X / _BucketSize), (long)Math.Floor(_Y / _BucketSize));
	}



	private static long Key(long _X, long _Y)
	{
		return (_X << 32) ^ (_Y & 0xFFFFFFFFL);
	}



	private static double Mod(double _Value, double _Modulo)
	{
		double result = _Value % _Modulo;

		return result < 0.0 ? result + _Modulo : result;
	}



	private static Vector2 Flat(Vector3 _Vector)
	{
		return new Vector2(_Vector.x, _Vector.y);
	}



	/// <summary>
	/// The outline's edges bucketed on a grid, for the "too close to the outline" test of the interior points.
	/// </summary>
	private sealed class EdgeGrid
	{
		private readonly double[] m_X;
		private readonly double[] m_Y;
		private readonly double m_Cell;
		private readonly double m_Clearance;
		private readonly Dictionary<long, List<int>> m_Buckets = new();



		public EdgeGrid(double[] _X, double[] _Y, double _Cell, double _Clearance)
		{
			m_X = _X;
			m_Y = _Y;
			m_Cell = _Cell;
			m_Clearance = _Clearance;

			int count = _X.Length;

			for (int i = 0; i < count; i++)
			{
				int next = (i + 1) % count;

				long fromX = (long)Math.Floor((Math.Min(_X[i], _X[next]) - _Clearance) / _Cell);
				long toX = (long)Math.Floor((Math.Max(_X[i], _X[next]) + _Clearance) / _Cell);
				long fromY = (long)Math.Floor((Math.Min(_Y[i], _Y[next]) - _Clearance) / _Cell);
				long toY = (long)Math.Floor((Math.Max(_Y[i], _Y[next]) + _Clearance) / _Cell);

				for (long y = fromY; y <= toY; y++)
				{
					for (long x = fromX; x <= toX; x++)
					{
						long key = Key(x, y);

						if (!m_Buckets.TryGetValue(key, out var list))
						{
							list = new List<int>();
							m_Buckets[key] = list;
						}

						list.Add(i);
					}
				}
			}
		}



		public bool IsNear(double _X, double _Y)
		{
			long key = Key((long)Math.Floor(_X / m_Cell), (long)Math.Floor(_Y / m_Cell));

			if (!m_Buckets.TryGetValue(key, out var list))
				return false;

			double limit = m_Clearance * m_Clearance;
			int count = m_X.Length;

			foreach (int i in list)
			{
				int next = (i + 1) % count;
				double ex = m_X[next] - m_X[i];
				double ey = m_Y[next] - m_Y[i];
				double lengthSquared = ex * ex + ey * ey;
				double t = lengthSquared > 0.0 ? Math.Clamp(((_X - m_X[i]) * ex + (_Y - m_Y[i]) * ey) / lengthSquared, 0.0, 1.0) : 0.0;
				double dx = m_X[i] + ex * t - _X;
				double dy = m_Y[i] + ey * t - _Y;

				if (dx * dx + dy * dy < limit)
					return true;
			}

			return false;
		}
	}



	/// <summary>
	/// A constrained triangulation: triangles winding counter-clockwise, each knowing its three neighbours (across the
	/// edge from its k-th vertex to the next) and which of its edges belong to the outline and must never be flipped.
	/// </summary>
	private sealed class Triangulation
	{
		public readonly List<double> X = new();
		public readonly List<double> Y = new();
		public readonly List<int> V = new();

		private readonly List<int> m_N = new();
		private readonly List<bool> m_C = new();
		private readonly Stack<int> m_Stack = new();
		private readonly double m_AreaEpsilon;

		public int PointCount => X.Count;

		public int TriangleCount => V.Count / 3;



		public Triangulation(double _Extent)
		{
			// Twice a triangle's area, below which it counts as flat.
			m_AreaEpsilon = Math.Max(1e-9, _Extent * _Extent * 1e-12);
		}



		public int AddPoint(double _X, double _Y)
		{
			X.Add(_X);
			Y.Add(_Y);

			return X.Count - 1;
		}



		private double Orient(int _A, int _B, int _C)
		{
			return (X[_B] - X[_A]) * (Y[_C] - Y[_A]) - (Y[_B] - Y[_A]) * (X[_C] - X[_A]);
		}



		private double Orient(int _A, int _B, double _X, double _Y)
		{
			return (X[_B] - X[_A]) * (_Y - Y[_A]) - (Y[_B] - Y[_A]) * (_X - X[_A]);
		}



		private int AddTriangle(int _A, int _B, int _C)
		{
			V.Add(_A);
			V.Add(_B);
			V.Add(_C);

			for (int k = 0; k < 3; k++)
			{
				m_N.Add(-1);
				m_C.Add(false);
			}

			return TriangleCount - 1;
		}



		private void SetTriangle(int _T, int _A, int _B, int _C, int _NAB, int _NBC, int _NCA, bool _CAB, bool _CBC, bool _CCA)
		{
			int i = _T * 3;

			V[i] = _A;
			V[i + 1] = _B;
			V[i + 2] = _C;

			m_N[i] = _NAB;
			m_N[i + 1] = _NBC;
			m_N[i + 2] = _NCA;

			m_C[i] = _CAB;
			m_C[i + 1] = _CBC;
			m_C[i + 2] = _CCA;
		}



		private void ReplaceNeighbour(int _Triangle, int _Old, int _New)
		{
			if (_Triangle < 0)
				return;

			for (int k = 0; k < 3; k++)
			{
				if (m_N[_Triangle * 3 + k] != _Old)
					continue;

				m_N[_Triangle * 3 + k] = _New;
				return;
			}
		}



		/// <summary>
		/// Ear clipping of the outline (points 0 to count - 1, counter-clockwise). An ear is a convex corner with no
		/// other point inside or on its triangle; testing "on" as well is what keeps a corner from being cut along a
		/// diagonal that runs through a point of a straight stretch, like the points across a river's row.
		/// </summary>
		public bool EarClip(int _Count)
		{
			int n = _Count;
			var prev = new int[n];
			var next = new int[n];
			var removed = new bool[n];
			var reflex = new bool[n];
			var reflexList = new List<int>();

			for (int i = 0; i < n; i++)
			{
				prev[i] = (i + n - 1) % n;
				next[i] = (i + 1) % n;
			}

			for (int i = 0; i < n; i++)
			{
				reflex[i] = Orient(prev[i], i, next[i]) <= m_AreaEpsilon;

				if (reflex[i])
					reflexList.Add(i);
			}

			int remaining = n;
			int current = 0;
			int misses = 0;

			while (remaining > 3)
			{
				if (IsEar(current, prev, next, reflex, removed, reflexList))
				{
					int a = prev[current];
					int c = next[current];

					AddTriangle(a, current, c);

					removed[current] = true;
					next[a] = c;
					prev[c] = a;
					remaining--;

					// Cutting an ear can only make its neighbours more convex.
					if (reflex[a] && Orient(prev[a], a, next[a]) > m_AreaEpsilon)
						reflex[a] = false;

					if (reflex[c] && Orient(prev[c], c, next[c]) > m_AreaEpsilon)
						reflex[c] = false;

					current = c;
					misses = 0;
					continue;
				}

				current = next[current];

				// A full lap without an ear: the outline crosses itself.
				if (++misses > remaining)
					return false;
			}

			int last = current;

			if (Orient(prev[last], last, next[last]) <= m_AreaEpsilon)
				return false;

			AddTriangle(prev[last], last, next[last]);

			return true;
		}



		private bool IsEar(int _I, int[] _Prev, int[] _Next, bool[] _Reflex, bool[] _Removed, List<int> _ReflexList)
		{
			int a = _Prev[_I];
			int c = _Next[_I];

			if (_Reflex[_I] || Orient(a, _I, c) <= m_AreaEpsilon)
				return false;

			foreach (int r in _ReflexList)
			{
				if (_Removed[r] || !_Reflex[r] || r == a || r == _I || r == c)
					continue;

				if (Orient(a, _I, r) >= -m_AreaEpsilon && Orient(_I, c, r) >= -m_AreaEpsilon && Orient(c, a, r) >= -m_AreaEpsilon)
					return false;
			}

			return true;
		}



		/// <summary>
		/// Links every triangle to its neighbours, and marks the outline's edges as constraints.
		/// </summary>
		public void BuildAdjacency(int _OutlineCount)
		{
			var edges = new Dictionary<long, int>(V.Count);

			for (int t = 0; t < TriangleCount; t++)
			{
				for (int k = 0; k < 3; k++)
					edges[EdgeKey(V[t * 3 + k], V[t * 3 + (k + 1) % 3])] = t;
			}

			for (int t = 0; t < TriangleCount; t++)
			{
				for (int k = 0; k < 3; k++)
				{
					int a = V[t * 3 + k];
					int b = V[t * 3 + (k + 1) % 3];

					m_N[t * 3 + k] = edges.TryGetValue(EdgeKey(b, a), out int other) ? other : -1;
					m_C[t * 3 + k] = a < _OutlineCount && b < _OutlineCount && (b == (a + 1) % _OutlineCount || a == (b + 1) % _OutlineCount);
				}
			}
		}



		private static long EdgeKey(int _A, int _B)
		{
			return ((long)_A << 32) | (uint)_B;
		}



		/// <summary>
		/// Lawson's flips over the whole triangulation: every edge that is not the outline's and fails the empty circle
		/// test is flipped, until none is left. The result is the constrained Delaunay triangulation of the outline.
		/// </summary>
		public void MakeDelaunay()
		{
			m_Stack.Clear();

			for (int t = 0; t < TriangleCount; t++)
			{
				for (int k = 0; k < 3; k++)
				{
					if (m_N[t * 3 + k] >= 0 && !m_C[t * 3 + k])
						m_Stack.Push(t * 3 + k);
				}
			}

			Legalize();
		}



		private void Legalize()
		{
			int guard = 0;
			int limit = V.Count * 40 + 1000;

			while (m_Stack.Count > 0 && guard++ < limit)
			{
				int edge = m_Stack.Pop();
				Flip(edge / 3, edge % 3);
			}

			m_Stack.Clear();
		}



		/// <summary>
		/// Flips the edge from the k-th vertex of a triangle to the next when it is not locally Delaunay: the triangles
		/// (a, b, c) and (b, a, d) become (a, d, c) and (d, b, c). The four edges round the pair are queued again.
		/// </summary>
		private bool Flip(int _T, int _K)
		{
			int edge = _T * 3 + _K;
			int u = m_N[edge];

			if (u < 0 || m_C[edge])
				return false;

			int a = V[_T * 3 + _K];
			int b = V[_T * 3 + (_K + 1) % 3];
			int c = V[_T * 3 + (_K + 2) % 3];

			int j = -1;

			for (int m = 0; m < 3; m++)
			{
				if (V[u * 3 + m] == b && V[u * 3 + (m + 1) % 3] == a)
				{
					j = m;
					break;
				}
			}

			if (j < 0)
				return false;

			int d = V[u * 3 + (j + 2) % 3];

			if (!InCircle(a, b, c, d))
				return false;

			// Never flip into a folded pair (only possible right at the precision limit).
			if (Orient(a, d, c) <= m_AreaEpsilon || Orient(d, b, c) <= m_AreaEpsilon)
				return false;

			int nbc = m_N[_T * 3 + (_K + 1) % 3];
			bool cbc = m_C[_T * 3 + (_K + 1) % 3];
			int nca = m_N[_T * 3 + (_K + 2) % 3];
			bool cca = m_C[_T * 3 + (_K + 2) % 3];
			int nad = m_N[u * 3 + (j + 1) % 3];
			bool cad = m_C[u * 3 + (j + 1) % 3];
			int ndb = m_N[u * 3 + (j + 2) % 3];
			bool cdb = m_C[u * 3 + (j + 2) % 3];

			SetTriangle(_T, a, d, c, nad, u, nca, cad, false, cca);
			SetTriangle(u, d, b, c, ndb, nbc, _T, cdb, cbc, false);

			ReplaceNeighbour(nad, u, _T);
			ReplaceNeighbour(nbc, _T, u);

			m_Stack.Push(_T * 3);
			m_Stack.Push(_T * 3 + 2);
			m_Stack.Push(u * 3);
			m_Stack.Push(u * 3 + 1);

			return true;
		}



		/// <summary>
		/// Whether d lies inside the circle through the counter-clockwise triangle (a, b, c), by more than the
		/// arithmetic can be trusted to (cocircular points, like a square of grid points, are left alone).
		/// </summary>
		private bool InCircle(int _A, int _B, int _C, int _D)
		{
			double adx = X[_A] - X[_D];
			double ady = Y[_A] - Y[_D];
			double bdx = X[_B] - X[_D];
			double bdy = Y[_B] - Y[_D];
			double cdx = X[_C] - X[_D];
			double cdy = Y[_C] - Y[_D];

			double ad = adx * adx + ady * ady;
			double bd = bdx * bdx + bdy * bdy;
			double cd = cdx * cdx + cdy * cdy;

			double t1 = ad * (bdx * cdy - cdx * bdy);
			double t2 = bd * (cdx * ady - adx * cdy);
			double t3 = cd * (adx * bdy - bdx * ady);

			double scale = Math.Abs(t1) + Math.Abs(t2) + Math.Abs(t3);

			return t1 + t2 + t3 > scale * 1e-10;
		}



		/// <summary>
		/// Walks from a triangle towards a point, across whichever edge has the point on its far side, until the point
		/// is inside. -1 when the walk runs off the outline (the straight line to the point leaves the surface).
		/// </summary>
		public int Locate(double _X, double _Y, int _Start)
		{
			if (TriangleCount == 0)
				return -1;

			int t = _Start >= 0 && _Start < TriangleCount ? _Start : 0;
			int limit = TriangleCount + 16;
			int rotate = 0;

			for (int step = 0; step < limit; step++)
			{
				int crossed = -1;

				for (int e = 0; e < 3; e++)
				{
					int k = (e + rotate) % 3;

					if (Orient(V[t * 3 + k], V[t * 3 + (k + 1) % 3], _X, _Y) < 0.0)
					{
						crossed = k;
						break;
					}
				}

				if (crossed < 0)
					return t;

				int neighbour = m_N[t * 3 + crossed];

				if (neighbour < 0)
					return -1;

				t = neighbour;
				rotate = (rotate + 1) % 3;
			}

			return -1;
		}



		/// <summary>
		/// The triangle a point is in, by testing them all. Only for when every walk failed.
		/// </summary>
		public int Scan(double _X, double _Y)
		{
			for (int t = 0; t < TriangleCount; t++)
			{
				if (Orient(V[t * 3], V[t * 3 + 1], _X, _Y) >= 0.0 &&
					Orient(V[t * 3 + 1], V[t * 3 + 2], _X, _Y) >= 0.0 &&
					Orient(V[t * 3 + 2], V[t * 3], _X, _Y) >= 0.0)
				{
					return t;
				}
			}

			return -1;
		}



		/// <summary>
		/// Splits a triangle at a point inside it and restores the Delaunay property round it. A point sitting on (or
		/// a hair off) an edge is nudged into the triangle first, so the split never makes a flat triangle. Returns
		/// the triangle the point ended up in, or -1 when it could not be placed.
		/// </summary>
		public int InsertAt(int _Triangle, double _X, double _Y, double _Spacing)
		{
			int t = _Triangle;
			int a = V[t * 3];
			int b = V[t * 3 + 1];
			int c = V[t * 3 + 2];

			double x = _X;
			double y = _Y;
			double clearance = Math.Max(1e-4, _Spacing * 1e-3);

			for (int attempt = 0; attempt < 4; attempt++)
			{
				if (EdgeDistance(a, b, x, y) > clearance && EdgeDistance(b, c, x, y) > clearance && EdgeDistance(c, a, x, y) > clearance)
					break;

				if (attempt == 3)
					return -1;

				// Towards the middle of the triangle, a little further each time.
				double cx = (X[a] + X[b] + X[c]) / 3.0;
				double cy = (Y[a] + Y[b] + Y[c]) / 3.0;
				double fraction = 0.01 * (attempt + 1);

				x += (cx - x) * fraction;
				y += (cy - y) * fraction;
			}

			int p = AddPoint(x, y);

			int nab = m_N[t * 3];
			int nbc = m_N[t * 3 + 1];
			int nca = m_N[t * 3 + 2];
			bool cab = m_C[t * 3];
			bool cbc = m_C[t * 3 + 1];
			bool cca = m_C[t * 3 + 2];

			int t1 = AddTriangle(b, c, p);
			int t2 = AddTriangle(c, a, p);

			SetTriangle(t, a, b, p, nab, t1, t2, cab, false, false);
			SetTriangle(t1, b, c, p, nbc, t2, t, cbc, false, false);
			SetTriangle(t2, c, a, p, nca, t, t1, cca, false, false);

			ReplaceNeighbour(nbc, t, t1);
			ReplaceNeighbour(nca, t, t2);

			m_Stack.Clear();
			m_Stack.Push(t * 3);
			m_Stack.Push(t1 * 3);
			m_Stack.Push(t2 * 3);

			Legalize();

			return t;
		}



		/// <summary>
		/// How far a point is inside a triangle's edge (negative outside it).
		/// </summary>
		private double EdgeDistance(int _A, int _B, double _X, double _Y)
		{
			double length = Math.Sqrt((X[_B] - X[_A]) * (X[_B] - X[_A]) + (Y[_B] - Y[_A]) * (Y[_B] - Y[_A]));

			return length > 0.0 ? Orient(_A, _B, _X, _Y) / length : 0.0;
		}
	}
}