Merge remote-tracking branch 'origin/terrain-for-networks'

14 files changed, 852 insertions, 504 deletions

diff --git a/game/gamestate.h b/game/gamestate.h
index 892aa69..189417d 100644
--- a/game/gamestate.h
+++ b/game/gamestate.h
@@ -6,7 +6,7 @@
 #include <special_members.h>
 
 class WorldObject;
-class GeoData;
+class Terrain;
 class Environment;
 
 class GameState {
@@ -17,7 +17,7 @@ public:
 	NO_COPY(GameState);
 
 	Collection<WorldObject> world;
-	std::shared_ptr<GeoData> geoData;
+	std::shared_ptr<Terrain> terrain;
 	std::shared_ptr<Environment> environment;
 	AssetFactory::Assets assets;
 };
diff --git a/game/geoData.cpp b/game/geoData.cpp
index a5fc4ef..4291a64 100644
--- a/game/geoData.cpp
+++ b/game/geoData.cpp
@@ -1,16 +1,15 @@
 #include "geoData.h"
 #include "collections.h"
 #include "geometricPlane.h"
+#include "util.h"
 #include <fstream>
 #include <glm/gtx/intersect.hpp>
 #include <maths.h>
 #include <ranges>
 #include <set>
-
-GeoData::GeoData()
-{
-	add_property(surface);
-}
+#ifndef NDEBUG
+#	include <stream_support.h>
+#endif
 
 GeoData
 GeoData::loadFromAsciiGrid(const std::filesystem::path & input)
@@ -110,118 +109,6 @@ GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistan
 	return mesh;
 }
 
-OpenMesh::FaceHandle
-GeoData::findPoint(GlobalPosition2D p) const
-{
-	return findPoint(p, *faces_sbegin());
-}
-
-GeoData::PointFace::PointFace(const GlobalPosition2D p, const GeoData * mesh) :
-	PointFace {p, mesh, *mesh->faces_sbegin()}
-{
-}
-
-GeoData::PointFace::PointFace(const GlobalPosition2D p, const GeoData * mesh, FaceHandle start) :
-	PointFace {p, mesh->findPoint(p, start)}
-{
-}
-
-GeoData::FaceHandle
-GeoData::PointFace::face(const GeoData * mesh, FaceHandle start) const
-{
-	if (_face.is_valid()) {
-		assert(mesh->triangleContainsPoint(point, _face));
-		return _face;
-	}
-	else {
-		return (_face = mesh->findPoint(point, start));
-	}
-}
-
-GeoData::FaceHandle
-GeoData::PointFace::face(const GeoData * mesh) const
-{
-	return face(mesh, *mesh->faces_sbegin());
-}
-
-namespace {
-	template<template<typename> typename Op>
-	[[nodiscard]] constexpr inline auto
-	pointLineOp(const GlobalPosition2D p, const GlobalPosition2D e1, const GlobalPosition2D e2)
-	{
-		return Op {}(CalcDistance(e2.x - e1.x) * CalcDistance(p.y - e1.y),
-				CalcDistance(e2.y - e1.y) * CalcDistance(p.x - e1.x));
-	}
-
-	constexpr auto pointLeftOfLine = pointLineOp<std::greater>;
-	constexpr auto pointLeftOfOrOnLine = pointLineOp<std::greater_equal>;
-
-	static_assert(pointLeftOfLine({1, 2}, {1, 1}, {2, 2}));
-	static_assert(pointLeftOfLine({2, 1}, {2, 2}, {1, 1}));
-	static_assert(pointLeftOfLine({2, 2}, {1, 2}, {2, 1}));
-	static_assert(pointLeftOfLine({1, 1}, {2, 1}, {1, 2}));
-	static_assert(pointLeftOfOrOnLine({310000000, 490000000}, {310000000, 490000000}, {310050000, 490050000}));
-	static_assert(pointLeftOfOrOnLine({310000000, 490000000}, {310050000, 490050000}, {310000000, 490050000}));
-	static_assert(pointLeftOfOrOnLine({310000000, 490000000}, {310000000, 490050000}, {310000000, 490000000}));
-
-	[[nodiscard]] constexpr inline bool
-	linesCross(
-			const GlobalPosition2D a1, const GlobalPosition2D a2, const GlobalPosition2D b1, const GlobalPosition2D b2)
-	{
-		return (pointLeftOfLine(a2, b1, b2) == pointLeftOfLine(a1, b2, b1))
-				&& (pointLeftOfLine(b1, a1, a2) == pointLeftOfLine(b2, a2, a1));
-	}
-
-	static_assert(linesCross({1, 1}, {2, 2}, {1, 2}, {2, 1}));
-	static_assert(linesCross({2, 2}, {1, 1}, {1, 2}, {2, 1}));
-
-	[[nodiscard]] constexpr inline bool
-	linesCrossLtR(
-			const GlobalPosition2D a1, const GlobalPosition2D a2, const GlobalPosition2D b1, const GlobalPosition2D b2)
-	{
-		return pointLeftOfLine(a2, b1, b2) && pointLeftOfLine(a1, b2, b1) && pointLeftOfLine(b1, a1, a2)
-				&& pointLeftOfLine(b2, a2, a1);
-	}
-
-	static_assert(linesCrossLtR({1, 1}, {2, 2}, {1, 2}, {2, 1}));
-	static_assert(!linesCrossLtR({2, 2}, {1, 1}, {1, 2}, {2, 1}));
-
-	constexpr GlobalPosition3D
-	positionOnTriangle(const GlobalPosition2D point, const GeoData::Triangle<3> & t)
-	{
-		const CalcPosition3D a = t[1] - t[0], b = t[2] - t[0];
-		const auto n = crossProduct(a, b);
-		return {point, ((n.x * t[0].x) + (n.y * t[0].y) + (n.z * t[0].z) - (n.x * point.x) - (n.y * point.y)) / n.z};
-	}
-
-	static_assert(positionOnTriangle({7, -2}, {{1, 2, 3}, {1, 0, 1}, {-2, 1, 0}}) == GlobalPosition3D {7, -2, 3});
-}
-
-OpenMesh::FaceHandle
-GeoData::findPoint(GlobalPosition2D p, OpenMesh::FaceHandle f) const
-{
-	while (f.is_valid() && !triangleContainsPoint(p, triangle<2>(f))) {
-		for (auto next = cfh_iter(f); next.is_valid(); ++next) {
-			f = opposite_face_handle(*next);
-			if (f.is_valid()) {
-				const auto e1 = point(to_vertex_handle(*next));
-				const auto e2 = point(to_vertex_handle(opposite_halfedge_handle(*next)));
-				if (pointLeftOfLine(p, e1, e2)) {
-					break;
-				}
-			}
-			f.reset();
-		}
-	}
-	return f;
-}
-
-GlobalPosition3D
-GeoData::positionAt(const PointFace & p) const
-{
-	return positionOnTriangle(p.point, triangle<3>(p.face(this)));
-}
-
 [[nodiscard]] GeoData::IntersectionResult
 GeoData::intersectRay(const Ray<GlobalPosition3D> & ray) const
 {
@@ -234,12 +121,12 @@ GeoData::intersectRay(const Ray<GlobalPosition3D> & ray, FaceHandle face) const
 	GeoData::IntersectionResult out;
 	walkUntil(PointFace {ray.start, face},
 			ray.start.xy() + (ray.direction.xy() * RelativePosition2D(upperExtent.xy() - lowerExtent.xy())),
-			[&out, &ray, this](FaceHandle face) {
+			[&out, &ray, this](const auto & step) {
 				BaryPosition bari {};
 				RelativeDistance dist {};
-				const auto t = triangle<3>(face);
+				const auto t = triangle<3>(step.current);
 				if (ray.intersectTriangle(t.x, t.y, t.z, bari, dist)) {
-					out.emplace(t * bari, face);
+					out.emplace(t * bari, step.current);
 					return true;
 				}
 				return false;
@@ -248,7 +135,7 @@ GeoData::intersectRay(const Ray<GlobalPosition3D> & ray, FaceHandle face) const
 }
 
 void
-GeoData::walk(const PointFace & from, const GlobalPosition2D to, const std::function<void(FaceHandle)> & op) const
+GeoData::walk(const PointFace & from, const GlobalPosition2D to, Consumer<WalkStep> op) const
 {
 	walkUntil(from, to, [&op](const auto & fh) {
 		op(fh);
@@ -257,41 +144,86 @@ GeoData::walk(const PointFace & from, const GlobalPosition2D to, const std::func
 }
 
 void
-GeoData::walkUntil(const PointFace & from, const GlobalPosition2D to, const std::function<bool(FaceHandle)> & op) const
+GeoData::walkUntil(const PointFace & from, const GlobalPosition2D to, Tester<WalkStep> op) const
+{
+	WalkStep step {
+			.current = from.face(this),
+	};
+	if (!step.current.is_valid()) {
+		const auto entryEdge = findEntry(from.point, to);
+		if (!entryEdge.is_valid()) {
+			return;
+		}
+		step.current = opposite_face_handle(entryEdge);
+	}
+	while (step.current.is_valid() && !op(step)) {
+		step.previous = step.current;
+		for (const auto next : fh_range(step.current)) {
+			step.current = opposite_face_handle(next);
+			if (step.current.is_valid() && step.current != step.previous) {
+				const auto nextPoints = points(toVertexHandles(next));
+				if (linesCrossLtR(from.point, to, nextPoints.second, nextPoints.first)) {
+					step.exitHalfedge = next;
+					step.exitPosition
+							= linesIntersectAt(from.point.xy(), to.xy(), nextPoints.second.xy(), nextPoints.first.xy())
+									  .value();
+					break;
+				}
+			}
+			step.current.reset();
+		}
+	}
+}
+
+void
+GeoData::walk(const PointFace & from, GlobalPosition2D to, GlobalPosition2D centre, Consumer<WalkStepCurve> op) const
+{
+	walkUntil(from, to, centre, [&op](const auto & fh) {
+		op(fh);
+		return false;
+	});
+}
+
+void
+GeoData::walkUntil(const PointFace & from, GlobalPosition2D to, GlobalPosition2D centre, Tester<WalkStepCurve> op) const
 {
-	auto f = from.face(this);
-	if (!f.is_valid()) {
+	WalkStepCurve step {WalkStep {.current = from.face(this)}};
+	if (!step.current.is_valid()) {
 		const auto entryEdge = findEntry(from.point, to);
 		if (!entryEdge.is_valid()) {
 			return;
 		}
-		f = opposite_face_handle(entryEdge);
+		step.current = opposite_face_handle(entryEdge);
 	}
-	FaceHandle previousFace;
-	while (f.is_valid() && !op(f)) {
-		for (auto next = cfh_iter(f); next.is_valid(); ++next) {
-			f = opposite_face_handle(*next);
-			if (f.is_valid() && f != previousFace) {
-				const auto e1 = point(to_vertex_handle(*next));
-				const auto e2 = point(to_vertex_handle(opposite_halfedge_handle(*next)));
-				if (linesCrossLtR(from.point, to, e1, e2)) {
-					previousFace = f;
+	ArcSegment arc {centre, from.point, to};
+	step.angle = arc.first;
+	while (step.current.is_valid() && !op(step)) {
+		step.previous = step.current;
+		for (const auto next : fh_range(step.current)) {
+			step.current = opposite_face_handle(next);
+			if (step.current.is_valid()) {
+				const auto e1 = point(to_vertex_handle(next));
+				const auto e2 = point(to_vertex_handle(opposite_halfedge_handle(next)));
+				if (const auto intersect = arc.crossesLineAt(e1, e2)) {
+					step.exitHalfedge = next;
+					arc.ep0 = step.exitPosition = intersect.value().first;
+					arc.first = std::nextafter(step.angle = intersect.value().second, INFINITY);
 					break;
 				}
 			}
-			f.reset();
+			step.current.reset();
 		}
 	}
 }
 
 void
-GeoData::boundaryWalk(const std::function<void(HalfedgeHandle)> & op) const
+GeoData::boundaryWalk(Consumer<HalfedgeHandle> op) const
 {
 	boundaryWalk(op, findBoundaryStart());
 }
 
 void
-GeoData::boundaryWalk(const std::function<void(HalfedgeHandle)> & op, HalfedgeHandle start) const
+GeoData::boundaryWalk(Consumer<HalfedgeHandle> op, HalfedgeHandle start) const
 {
 	assert(is_boundary(start));
 	boundaryWalkUntil(
@@ -303,13 +235,13 @@ GeoData::boundaryWalk(const std::function<void(HalfedgeHandle)> & op, HalfedgeHa
 }
 
 void
-GeoData::boundaryWalkUntil(const std::function<bool(HalfedgeHandle)> & op) const
+GeoData::boundaryWalkUntil(Tester<HalfedgeHandle> op) const
 {
 	boundaryWalkUntil(op, findBoundaryStart());
 }
 
 void
-GeoData::boundaryWalkUntil(const std::function<bool(HalfedgeHandle)> & op, HalfedgeHandle start) const
+GeoData::boundaryWalkUntil(Tester<HalfedgeHandle> op, HalfedgeHandle start) const
 {
 	assert(is_boundary(start));
 	if (!op(start)) {
@@ -337,45 +269,6 @@ GeoData::findEntry(const GlobalPosition2D from, const GlobalPosition2D to) const
 	return entry;
 }
 
-bool
-GeoData::triangleContainsPoint(const GlobalPosition2D p, const Triangle<2> & t)
-{
-	return pointLeftOfOrOnLine(p, t[0], t[1]) && pointLeftOfOrOnLine(p, t[1], t[2])
-			&& pointLeftOfOrOnLine(p, t[2], t[0]);
-}
-
-bool
-GeoData::triangleContainsPoint(const GlobalPosition2D p, FaceHandle face) const
-{
-	return triangleContainsPoint(p, triangle<2>(face));
-}
-
-GeoData::HalfedgeHandle
-GeoData::findBoundaryStart() const
-{
-	return *std::find_if(halfedges_sbegin(), halfedges_end(), [this](const auto heh) {
-		return is_boundary(heh);
-	});
-}
-
-[[nodiscard]] RelativePosition3D
-GeoData::difference(const HalfedgeHandle heh) const
-{
-	return ::difference(point(to_vertex_handle(heh)), point(from_vertex_handle(heh)));
-}
-
-[[nodiscard]] RelativeDistance
-GeoData::length(const HalfedgeHandle heh) const
-{
-	return glm::length(difference(heh));
-}
-
-[[nodiscard]] GlobalPosition3D
-GeoData::centre(const HalfedgeHandle heh) const
-{
-	return point(from_vertex_handle(heh)) + (difference(heh) / 2.F);
-}
-
 void
 GeoData::updateAllVertexNormals()
 {
@@ -399,178 +292,296 @@ GeoData::updateVertexNormal(VertexHandle vertex)
 	set_normal(vertex, glm::normalize(n));
 }
 
-bool
-GeoData::triangleOverlapsTriangle(const Triangle<2> & a, const Triangle<2> & b)
+OpenMesh::VertexHandle
+GeoData::setPoint(GlobalPosition3D tsPoint, const RelativeDistance nearNodeTolerance)
 {
-	return triangleContainsPoint(a.x, b) || triangleContainsPoint(a.y, b) || triangleContainsPoint(a.z, b)
-			|| triangleContainsPoint(b.x, a) || triangleContainsPoint(b.y, a) || triangleContainsPoint(b.z, a)
-			|| linesCross(a.x, a.y, b.x, b.y) || linesCross(a.x, a.y, b.y, b.z) || linesCross(a.x, a.y, b.z, b.x)
-			|| linesCross(a.y, a.z, b.x, b.y) || linesCross(a.y, a.z, b.y, b.z) || linesCross(a.y, a.z, b.z, b.x)
-			|| linesCross(a.z, a.x, b.x, b.y) || linesCross(a.z, a.x, b.y, b.z) || linesCross(a.z, a.x, b.z, b.x);
-}
-
-bool
-GeoData::triangleContainsTriangle(const Triangle<2> & a, const Triangle<2> & b)
-{
-	return triangleContainsPoint(a.x, b) && triangleContainsPoint(a.y, b) && triangleContainsPoint(a.z, b);
-}
+	const auto face = findPoint(tsPoint);
+	const auto distFromTsPoint = vertexDistanceFunction<2>(tsPoint);
+	// Check vertices
+	if (const auto nearest = std::ranges::min(fv_range(face) | std::views::transform(distFromTsPoint), {}, GetSecond);
+			nearest.second < nearNodeTolerance) {
+		point(nearest.first).z = tsPoint.z;
+		return nearest.first;
+	}
+	// Check edges
+	if (const auto nearest = std::ranges::min(fh_range(face) | std::views::transform(distFromTsPoint), {}, GetSecond);
+			nearest.second < nearNodeTolerance) {
+		const auto from = point(from_vertex_handle(nearest.first)).xy();
+		const auto to = point(to_vertex_handle(nearest.first)).xy();
+		const auto v = vector_normal(from - to);
+		const auto inter = linesIntersectAt(from, to, tsPoint.xy(), tsPoint.xy() + v);
+		if (!inter) [[unlikely]] {
+			throw std::runtime_error("Perpendicular lines do not cross");
+		}
+		return split_copy(edge_handle(nearest.first), *inter || tsPoint.z);
+	}
+	// Nothing close, split face
+	return split_copy(face, tsPoint);
+};
 
-void
+std::set<GeoData::FaceHandle>
 GeoData::setHeights(const std::span<const GlobalPosition3D> triangleStrip, const SetHeightsOpts & opts)
 {
 	if (triangleStrip.size() < 3) {
-		return;
+		return {};
 	}
 	const auto stripMinMax = std::ranges::minmax(triangleStrip, {}, &GlobalPosition3D::z);
 	lowerExtent.z = std::min(upperExtent.z, stripMinMax.min.z);
 	upperExtent.z = std::max(upperExtent.z, stripMinMax.max.z);
 
-	const auto vertexDistFrom = [this](GlobalPosition2D p) {
-		return [p, this](const VertexHandle v) {
-			return std::make_pair(v, glm::length(::difference(p, this->point(v).xy())));
-		};
-	};
+	class SetHeights {
+	public:
+		SetHeights(GeoData * geoData, const std::span<const GlobalPosition3D> triangleStrip) :
+			geoData(geoData), triangleStrip {triangleStrip},
+			strip {materializeRange(triangleStrip | triangleTriples | std::views::transform([](const auto & newVert) {
+				return std::make_from_tuple<Triangle<3>>(newVert);
+			}))}
+		{
+		}
 
-	std::set<VertexHandle> newOrChangedVerts;
-	auto addVertexForNormalUpdate = [this, &newOrChangedVerts](const VertexHandle vertex) {
-		newOrChangedVerts.emplace(vertex);
-		std::ranges::copy(vv_range(vertex), std::inserter(newOrChangedVerts, newOrChangedVerts.end()));
-	};
+		std::vector<VertexHandle>
+		createVerticesForStrip(RelativeDistance nearNodeTolerance)
+		{
+			// New vertices for each vertex in triangleStrip
+			const auto newVerts
+					= materializeRange(triangleStrip | std::views::transform([this, nearNodeTolerance](auto v) {
+						  return geoData->setPoint(v, nearNodeTolerance);
+					  }));
+			std::ranges::copy(newVerts, std::inserter(newOrChangedVerts, newOrChangedVerts.end()));
+#ifndef NDEBUG
+			geoData->sanityCheck();
+#endif
+			return newVerts;
+		}
+
+		const Triangle<3> *
+		getTriangle(const GlobalPosition2D point) const
+		{
+			if (const auto t = std::ranges::find_if(strip,
+						[point](const auto & triangle) {
+							return triangle.containsPoint(point);
+						});
+					t != strip.end()) {
+				return &*t;
+			}
+			return nullptr;
+		}
 
-	// New vertices for each vertex in triangleStrip
-	std::vector<VertexHandle> newVerts;
-	newVerts.reserve(newVerts.size());
-	std::transform(triangleStrip.begin(), triangleStrip.end(), std::back_inserter(newVerts),
-			[this, &newVerts, &vertexDistFrom, &opts](const auto tsPoint) {
-				const auto face = findPoint(tsPoint);
-				if (const auto nearest = std::ranges::min(std::views::iota(fv_begin(face), fv_end(face))
-									| std::views::transform(vertexDistFrom(tsPoint)),
-							{}, &std::pair<VertexHandle, float>::second);
-						nearest.second < opts.nearNodeTolerance && !std::ranges::contains(newVerts, nearest.first)) {
-					point(nearest.first) = tsPoint;
-					return nearest.first;
+		void
+		doBoundaryPart(VertexHandle start, VertexHandle end, const Triangle<3> & triangle,
+				const RelativeDistance nearNodeTolerance)
+		{
+			boundaryTriangles.emplace(start, &triangle);
+			const auto endPoint = geoData->point(end);
+			while (!std::ranges::contains(geoData->vv_range(start), end)) {
+				const auto startPoint = geoData->point(start);
+				const auto distanceToEndPoint = distance(startPoint.xy(), endPoint.xy());
+				if (std::ranges::any_of(geoData->vv_range(start), [&](const auto & adjVertex) {
+						const auto adjPoint = geoData->point(adjVertex);
+						if (distance(adjPoint.xy(), endPoint.xy()) < distanceToEndPoint
+								&& (Triangle<2> {startPoint, endPoint, adjPoint}.area()
+										   / distance(startPoint.xy(), endPoint.xy()))
+										< nearNodeTolerance) {
+							start = adjVertex;
+							newOrChangedVerts.emplace(start);
+							boundaryTriangles.emplace(start, &triangle);
+							geoData->point(start).z = triangle.positionOnPlane(adjPoint).z;
+							return true;
+						}
+						return false;
+					})) {
+					continue;
 				}
-				return split(face, tsPoint);
-			});
-	std::ranges::for_each(newVerts, addVertexForNormalUpdate);
-
-	// Create temporary triangles from triangleStrip
-	std::vector<Triangle<3>> strip;
-	std::transform(
-			strip_begin(triangleStrip), strip_end(triangleStrip), std::back_inserter(strip), [](const auto & newVert) {
-				const auto [a, b, c] = newVert;
-				return Triangle<3> {a, b, c};
-			});
-	auto getTriangle = [&strip](const auto point) -> const Triangle<3> * {
-		if (const auto t = std::ranges::find_if(strip,
-					[point](const auto & triangle) {
-						return triangleContainsPoint(point, triangle);
+				if (std::ranges::any_of(geoData->voh_range(start), [&](const auto & outHalf) {
+						const auto next = geoData->next_halfedge_handle(outHalf);
+						const auto nexts
+								= std::array {geoData->from_vertex_handle(next), geoData->to_vertex_handle(next)};
+						const auto nextPoints = nexts | std::views::transform([this](const auto v) {
+							return std::make_pair(v, geoData->point(v));
+						});
+						if (linesCross(startPoint, endPoint, nextPoints.front().second, nextPoints.back().second)) {
+							if (const auto intersection = linesIntersectAt(startPoint.xy(), endPoint.xy(),
+										nextPoints.front().second.xy(), nextPoints.back().second.xy())) {
+								if (const auto nextEdge = geoData->shouldFlip(next, startPoint)) {
+									geoData->flip(*nextEdge);
+									return true;
+								}
+								start = geoData->split_copy(
+										geoData->edge_handle(next), triangle.positionOnPlane(*intersection));
+								newOrChangedVerts.emplace(start);
+								boundaryTriangles.emplace(start, &triangle);
+								return true;
+							}
+							throw std::runtime_error("Crossing lines don't intersect");
+						}
+						return false;
+					})) {
+					continue;
+				}
+#ifndef NDEBUG
+				CLOG(start);
+				CLOG(startPoint);
+				CLOG(end);
+				CLOG(endPoint);
+				for (const auto v : geoData->vv_range(start)) {
+					CLOG(geoData->point(v));
+				}
+				geoData->sanityCheck();
+#endif
+				throw std::runtime_error(
+						std::format("Could not navigate to ({}, {}, {})", endPoint.x, endPoint.y, endPoint.z));
+			}
+		}
+
+		void
+		cutBoundary(const std::vector<VertexHandle> & newVerts, RelativeDistance nearNodeTolerance)
+		{
+			// Cut along each edge of triangleStrip AB, AC, BC, BD, CD, CE etc
+			std::ranges::for_each(newVerts | std::views::adjacent<3>,
+					[this, nearNodeTolerance, triangle = strip.begin()](const auto & verts) mutable {
+						const auto & [a, _, c] = verts;
+						doBoundaryPart(a, c, *triangle, nearNodeTolerance);
+						triangle++;
 					});
-				t != strip.end()) {
-			return &*t;
+			doBoundaryPart(*++newVerts.begin(), newVerts.front(), strip.front(), nearNodeTolerance);
+			doBoundaryPart(*++newVerts.rbegin(), newVerts.back(), strip.back(), nearNodeTolerance);
 		}
-		return nullptr;
-	};
 
-	// Cut along each edge of triangleStrip AB, AC, BC, BD, CD, CE etc
-	std::map<VertexHandle, const Triangle<3> *> boundaryTriangles;
-	auto doBoundaryPart = [this, &boundaryTriangles, &newVerts, &vertexDistFrom, &opts, &addVertexForNormalUpdate](
-								  VertexHandle start, VertexHandle end, const Triangle<3> & triangle) {
-		boundaryTriangles.emplace(start, &triangle);
-		const auto endPoint = point(end);
-		while (!std::ranges::contains(vv_range(start), end)
-				&& std::ranges::any_of(voh_range(start), [&](const auto & outHalf) {
-					   const auto next = next_halfedge_handle(outHalf);
-					   const auto startPoint = point(start);
-					   const auto nexts = std::array {from_vertex_handle(next), to_vertex_handle(next)};
-					   const auto nextPoints = nexts | std::views::transform([this](const auto v) {
-						   return std::make_pair(v, this->point(v));
-					   });
-					   if (linesCross(startPoint, endPoint, nextPoints.front().second, nextPoints.back().second)) {
-						   if (const auto intersection = linesIntersectAt(startPoint.xy(), endPoint.xy(),
-									   nextPoints.front().second.xy(), nextPoints.back().second.xy())) {
-							   if (const auto nextDist
-									   = std::ranges::min(nexts | std::views::transform(vertexDistFrom(*intersection)),
-											   {}, &std::pair<VertexHandle, float>::second);
-									   nextDist.second < opts.nearNodeTolerance
-									   && !boundaryTriangles.contains(nextDist.first)
-									   && !std::ranges::contains(newVerts, nextDist.first)) {
-								   start = nextDist.first;
-								   point(start) = positionOnTriangle(*intersection, triangle);
-							   }
-							   else {
-								   start = split(edge_handle(next), positionOnTriangle(*intersection, triangle));
-							   }
-							   addVertexForNormalUpdate(start);
-							   boundaryTriangles.emplace(start, &triangle);
-							   return true;
-						   }
-					   }
-					   return false;
-				   })) { }
-	};
-	auto doBoundary = [&doBoundaryPart, triangle = strip.begin()](const auto & verts) mutable {
-		const auto & [a, b, c] = verts;
-		doBoundaryPart(a, b, *triangle);
-		doBoundaryPart(a, c, *triangle);
-		triangle++;
-	};
-	std::ranges::for_each(newVerts | std::views::adjacent<3>, doBoundary);
-	doBoundaryPart(*++newVerts.rbegin(), newVerts.back(), *strip.rbegin());
-
-	std::set<HalfedgeHandle> done;
-	std::set<HalfedgeHandle> todo;
-	auto todoOutHalfEdges = [&todo, &done, this](const VertexHandle v) {
-		std::copy_if(voh_begin(v), voh_end(v), std::inserter(todo, todo.end()), [&done](const auto & h) {
-			return !done.contains(h);
-		});
-	};
-	std::ranges::for_each(newVerts, todoOutHalfEdges);
-	while (!todo.empty()) {
-		const auto heh = todo.extract(todo.begin()).value();
-		const auto fromVertex = from_vertex_handle(heh);
-		const auto toVertex = to_vertex_handle(heh);
-		const auto & fromPoint = point(fromVertex);
-		auto & toPoint = point(toVertex);
-		auto toTriangle = getTriangle(toPoint);
-		if (!toTriangle) {
-			if (const auto boundaryVertex = boundaryTriangles.find(toVertex);
-					boundaryVertex != boundaryTriangles.end()) {
-				toTriangle = boundaryVertex->second;
+		using HeightSetTodo = std::multimap<VertexHandle, VertexHandle>;
+
+		HeightSetTodo
+		setSurfaceHeights(const std::span<const VertexHandle> newVerts)
+		{
+			HeightSetTodo out;
+			auto setSurfaceVertexHeight = [this, &out](const auto & setSurfaceVertexHeight, const VertexHandle vertex,
+												  const VertexHandle previousVertex) -> void {
+				if (surfaceVerts.contains(vertex)) {
+					return;
+				}
+				auto & point = geoData->point(vertex);
+				auto triangle = getTriangle(point);
+				if (!triangle) {
+					if (const auto boundaryVertex = boundaryTriangles.find(vertex);
+							boundaryVertex != boundaryTriangles.end()) {
+						triangle = boundaryVertex->second;
+					}
+				}
+				if (triangle) { // point within the new strip, adjust vertically by triangle
+					point.z = triangle->positionOnPlane(point).z;
+					newOrChangedVerts.emplace(vertex);
+					surfaceVerts.emplace(vertex);
+					for (const auto nextVertex : geoData->vv_range(vertex)) {
+						setSurfaceVertexHeight(setSurfaceVertexHeight, nextVertex, vertex);
+					}
+				}
+				else if (previousVertex.is_valid()) {
+					out.emplace(vertex, previousVertex);
+				}
+			};
+			for (const auto vertex : newVerts) {
+				setSurfaceVertexHeight(setSurfaceVertexHeight, vertex, VertexHandle {});
+			}
+			return out;
+		}
+
+		void
+		setHalfedgeToHeight(HeightSetTodo & nexts, HeightSetTodo::const_iterator verticesBegin,
+				HeightSetTodo::const_iterator verticesEnd, const RelativeDistance maxSlope)
+		{
+			const auto vertex = verticesBegin->first;
+			auto & point = geoData->point(vertex);
+			const auto minMaxHeight = std::accumulate(verticesBegin, verticesEnd,
+					std::pair<GlobalDistance, GlobalDistance> {
+							std::numeric_limits<GlobalDistance>::min(),
+							std::numeric_limits<GlobalDistance>::max(),
+					},
+					[this, maxSlope, point](auto limit, auto previousVertexItr) {
+						const auto & fromPoint = geoData->point(previousVertexItr.second);
+						const auto maxOffset = static_cast<GlobalDistance>(
+								std::round(maxSlope * ::distance<2>(fromPoint.xy(), point.xy())));
+						limit.first = std::max(limit.first, fromPoint.z - maxOffset);
+						limit.second = std::min(limit.second, fromPoint.z + maxOffset);
+						return limit;
+					});
+
+			const auto newHeight = std::clamp(point.z, minMaxHeight.first, minMaxHeight.second);
+			if (newHeight != point.z) {
+				point.z = newHeight;
+				newOrChangedVerts.emplace(vertex);
+				for (const auto nextVertex : geoData->vv_range(vertex)) {
+					if (!std::ranges::contains(verticesBegin, verticesEnd, nextVertex, GetSecond)
+							&& !boundaryTriangles.contains(nextVertex)) {
+						nexts.emplace(nextVertex, vertex);
+					}
+				}
 			}
 		}
-		if (toTriangle) { // point within the new strip, adjust vertically by triangle
-			toPoint.z = positionOnTriangle(toPoint, *toTriangle).z;
-			addVertexForNormalUpdate(toVertex);
-			todoOutHalfEdges(toVertex);
+
+		void
+		setHeightsAsRequired(HeightSetTodo starts, RelativeDistance maxSlope)
+		{
+			while (!starts.empty()) {
+				HeightSetTodo nexts;
+
+				for (const auto chunk : starts | std::views::chunk_by([](const auto a, const auto b) {
+						 return a.first == b.first;
+					 })) {
+					setHalfedgeToHeight(nexts, chunk.begin(), chunk.end(), maxSlope);
+				}
+				starts = std::move(nexts);
+			}
 		}
-		else if (!toTriangle) { // point without the new strip, adjust vertically by limit
-			const auto maxOffset = static_cast<GlobalDistance>(opts.maxSlope * glm::length(difference(heh).xy()));
-			const auto newHeight = std::clamp(toPoint.z, fromPoint.z - maxOffset, fromPoint.z + maxOffset);
-			if (newHeight != toPoint.z) {
-				toPoint.z = newHeight;
-				addVertexForNormalUpdate(toVertex);
-				std::copy_if(voh_begin(toVertex), voh_end(toVertex), std::inserter(todo, todo.end()),
-						[this, &boundaryTriangles](const auto & heh) {
-							return !boundaryTriangles.contains(to_vertex_handle(heh));
-						});
+
+		std::set<FaceHandle>
+		setSurface(const Surface * surface)
+		{
+			std::set<FaceHandle> out;
+			auto surfaceStripWalk = [this, surface, &out](const auto & surfaceStripWalk, const auto & face) -> void {
+				if (!out.contains(face)) {
+					geoData->property(geoData->surface, face) = surface;
+					out.emplace(face);
+					std::ranges::for_each(
+							geoData->ff_range(face), [this, &surfaceStripWalk](const auto & adjacentFaceHandle) {
+								if (getTriangle(geoData->triangle<2>(adjacentFaceHandle).centroid())) {
+									surfaceStripWalk(surfaceStripWalk, adjacentFaceHandle);
+								}
+							});
+				}
+			};
+			for (const auto & triangle : strip) {
+				surfaceStripWalk(surfaceStripWalk, geoData->findPoint(triangle.centroid()));
 			}
+			return out;
 		}
-		done.insert(heh);
-	}
 
-	auto surfaceStripWalk = [this, &getTriangle, &opts](const auto & surfaceStripWalk, const auto & face) -> void {
-		if (!property(surface, face)) {
-			property(surface, face) = &opts.surface;
-			std::ranges::for_each(
-					ff_range(face), [this, &getTriangle, &surfaceStripWalk](const auto & adjacentFaceHandle) {
-						if (getTriangle(this->triangle<2>(adjacentFaceHandle).centroid())) {
-							surfaceStripWalk(surfaceStripWalk, adjacentFaceHandle);
-						}
-					});
+		std::set<GeoData::FaceHandle>
+		run(const SetHeightsOpts & opts)
+		{
+			const std::vector<VertexHandle> newVerts = createVerticesForStrip(opts.nearNodeTolerance);
+			cutBoundary(newVerts, opts.nearNodeTolerance);
+			setHeightsAsRequired(setSurfaceHeights(newVerts), opts.maxSlope);
+			const auto out = setSurface(opts.surface);
+
+			geoData->expandVerts(newOrChangedVerts);
+			geoData->updateAllVertexNormals(newOrChangedVerts);
+			geoData->afterChange();
+
+			return out;
 		}
+
+	private:
+		GeoData * geoData;
+		const std::span<const GlobalPosition3D> triangleStrip;
+		const std::vector<Triangle<3>> strip;
+		std::set<VertexHandle> newOrChangedVerts;
+		std::set<VertexHandle> surfaceVerts;
+		std::map<VertexHandle, const Triangle<3> *> boundaryTriangles;
 	};
-	surfaceStripWalk(surfaceStripWalk, findPoint(strip.front().centroid()));
 
-	updateAllVertexNormals(newOrChangedVerts);
+	return SetHeights {this, triangleStrip}.run(opts);
+}
+
+void
+GeoData::afterChange()
+{
 }
diff --git a/game/geoData.h b/game/geoData.h
index 79924d3..b2a75bd 100644
--- a/game/geoData.h
+++ b/game/geoData.h
@@ -1,89 +1,61 @@
 #pragma once
 
 #include "collections.h" // IWYU pragma: keep IterableCollection
-#include "config/types.h"
-#include "ray.h"
+#include "geoDataMesh.h"
 #include "surface.h"
-#include "triangle.h"
-#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
 #include <filesystem>
 #include <glm/vec2.hpp>
-#include <optional>
-#include <thirdparty/openmesh/glmcompat.h>
-
-struct GeoDataTraits : public OpenMesh::DefaultTraits {
-	FaceAttributes(OpenMesh::Attributes::Status);
-	EdgeAttributes(OpenMesh::Attributes::Status);
-	VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Status);
-	HalfedgeAttributes(OpenMesh::Attributes::Status);
-	using Point = GlobalPosition3D;
-	using Normal = Normal3D;
-};
 
-class GeoData : public OpenMesh::TriMesh_ArrayKernelT<GeoDataTraits> {
+class GeoData : public GeoDataMesh {
 private:
-	GeoData();
-
-	OpenMesh::FPropHandleT<const Surface *> surface;
+	const OpenMesh::Helpers::Property<const Surface *, OpenMesh::FPropHandleT> surface {this};
 
 public:
 	static GeoData loadFromAsciiGrid(const std::filesystem::path &);
 	static GeoData createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistance h);
 
-	struct PointFace {
-		// NOLINTNEXTLINE(hicpp-explicit-conversions)
-		PointFace(const GlobalPosition2D p) : point {p} { }
-
-		PointFace(const GlobalPosition2D p, FaceHandle face) : point {p}, _face {face} { }
-
-		PointFace(const GlobalPosition2D p, const GeoData *);
-		PointFace(const GlobalPosition2D p, const GeoData *, FaceHandle start);
-
-		const GlobalPosition2D point;
-		[[nodiscard]] FaceHandle face(const GeoData *) const;
-		[[nodiscard]] FaceHandle face(const GeoData *, FaceHandle start) const;
-
-		[[nodiscard]] bool
-		isLocated() const
-		{
-			return _face.is_valid();
-		}
-
-	private:
-		mutable FaceHandle _face {};
-	};
-
-	template<glm::length_t Dim> using Triangle = ::Triangle<Dim, GlobalDistance>;
-
-	[[nodiscard]] FaceHandle findPoint(GlobalPosition2D) const;
-	[[nodiscard]] FaceHandle findPoint(GlobalPosition2D, FaceHandle start) const;
-
-	[[nodiscard]] GlobalPosition3D positionAt(const PointFace &) const;
 	using IntersectionLocation = std::pair<GlobalPosition3D, FaceHandle>;
 	using IntersectionResult = std::optional<IntersectionLocation>;
 	[[nodiscard]] IntersectionResult intersectRay(const Ray<GlobalPosition3D> &) const;
 	[[nodiscard]] IntersectionResult intersectRay(const Ray<GlobalPosition3D> &, FaceHandle start) const;
 
-	void walk(const PointFace & from, const GlobalPosition2D to, const std::function<void(FaceHandle)> & op) const;
-	void walkUntil(const PointFace & from, const GlobalPosition2D to, const std::function<bool(FaceHandle)> & op) const;
+	struct WalkStep {
+		FaceHandle current;
+		FaceHandle previous {};
+		HalfedgeHandle exitHalfedge {};
+		GlobalPosition2D exitPosition {};
+	};
+
+	struct WalkStepCurve : public WalkStep {
+		Angle angle {};
+	};
+
+	template<typename T> using Consumer = const std::function<void(const T &)> &;
+	template<typename T> using Tester = const std::function<bool(const T &)> &;
+
+	void walk(const PointFace & from, GlobalPosition2D to, Consumer<WalkStep> op) const;
+	void walkUntil(const PointFace & from, GlobalPosition2D to, Tester<WalkStep> op) const;
+	void walk(const PointFace & from, GlobalPosition2D to, GlobalPosition2D centre, Consumer<WalkStepCurve> op) const;
+	void walkUntil(
+			const PointFace & from, GlobalPosition2D to, GlobalPosition2D centre, Tester<WalkStepCurve> op) const;
 
-	void boundaryWalk(const std::function<void(HalfedgeHandle)> &) const;
-	void boundaryWalk(const std::function<void(HalfedgeHandle)> &, HalfedgeHandle start) const;
-	void boundaryWalkUntil(const std::function<bool(HalfedgeHandle)> &) const;
-	void boundaryWalkUntil(const std::function<bool(HalfedgeHandle)> &, HalfedgeHandle start) const;
+	void boundaryWalk(Consumer<HalfedgeHandle>) const;
+	void boundaryWalk(Consumer<HalfedgeHandle>, HalfedgeHandle start) const;
+	void boundaryWalkUntil(Tester<HalfedgeHandle>) const;
+	void boundaryWalkUntil(Tester<HalfedgeHandle>, HalfedgeHandle start) const;
 
-	[[nodiscard]] HalfedgeHandle findEntry(const GlobalPosition2D from, const GlobalPosition2D to) const;
+	[[nodiscard]] HalfedgeHandle findEntry(GlobalPosition2D from, GlobalPosition2D to) const;
 
 	struct SetHeightsOpts {
 		static constexpr auto DEFAULT_NEAR_NODE_TOLERANACE = 500.F;
 		static constexpr auto DEFAULT_MAX_SLOPE = 0.5F;
 
-		const Surface & surface;
+		const Surface * surface = nullptr;
 		RelativeDistance nearNodeTolerance = DEFAULT_NEAR_NODE_TOLERANACE;
 		RelativeDistance maxSlope = DEFAULT_MAX_SLOPE;
 	};
 
-	void setHeights(std::span<const GlobalPosition3D> triangleStrip, const SetHeightsOpts &);
+	std::set<FaceHandle> setHeights(std::span<const GlobalPosition3D> triangleStrip, const SetHeightsOpts &);
 
 	[[nodiscard]] auto
 	getExtents() const
@@ -92,37 +64,20 @@ public:
 	}
 
 	template<typename HandleT>
+		requires(std::derived_from<HandleT, OpenMesh::BaseHandle>)
 	[[nodiscard]] auto
-	get_surface(const HandleT h)
+	getSurface(const HandleT handle) const
 	{
-		return property(surface, h);
+		assert(handle.is_valid());
+		return property(surface, handle);
 	}
 
 protected:
-	template<glm::length_t Dim>
-	[[nodiscard]] Triangle<Dim>
-	triangle(FaceHandle face) const
-	{
-		Triangle<Dim> triangle;
-		std::ranges::transform(fv_range(face), triangle.begin(), [this](auto vertex) {
-			return point(vertex);
-		});
-		return triangle;
-	}
-
-	[[nodiscard]] static bool triangleContainsPoint(const GlobalPosition2D, const Triangle<2> &);
-	[[nodiscard]] bool triangleContainsPoint(const GlobalPosition2D, FaceHandle) const;
-	[[nodiscard]] static bool triangleOverlapsTriangle(const Triangle<2> &, const Triangle<2> &);
-	[[nodiscard]] static bool triangleContainsTriangle(const Triangle<2> &, const Triangle<2> &);
-	[[nodiscard]] HalfedgeHandle findBoundaryStart() const;
-	[[nodiscard]] RelativePosition3D difference(const HalfedgeHandle) const;
-
-	[[nodiscard]] RelativeDistance length(const HalfedgeHandle) const;
-	[[nodiscard]] GlobalPosition3D centre(const HalfedgeHandle) const;
-
+	[[nodiscard]] VertexHandle setPoint(GlobalPosition3D point, RelativeDistance nearNodeTolerance);
 	void updateAllVertexNormals();
 	template<std::ranges::range R> void updateAllVertexNormals(const R &);
 	void updateVertexNormal(VertexHandle);
+	virtual void afterChange();
 
 private:
 	GlobalPosition3D lowerExtent {}, upperExtent {};
diff --git a/game/geoDataMesh.cpp b/game/geoDataMesh.cpp
new file mode 100644
index 0000000..8107a5e
--- /dev/null
+++ b/game/geoDataMesh.cpp
@@ -0,0 +1,150 @@
+#include "geoDataMesh.h"
+#include <format>
+#ifndef NDEBUG
+#	include <stream_support.h>
+#endif
+
+OpenMesh::FaceHandle
+GeoDataMesh::findPoint(GlobalPosition2D coord) const
+{
+	return findPoint(coord, *faces_sbegin());
+}
+
+GeoDataMesh::PointFace::PointFace(const GlobalPosition2D coord, const GeoDataMesh * mesh) :
+	PointFace {coord, mesh, *mesh->faces_sbegin()}
+{
+}
+
+GeoDataMesh::PointFace::PointFace(const GlobalPosition2D coord, const GeoDataMesh * mesh, FaceHandle start) :
+	PointFace {coord, mesh->findPoint(coord, start)}
+{
+}
+
+OpenMesh::FaceHandle
+GeoDataMesh::PointFace::face(const GeoDataMesh * mesh, FaceHandle start) const
+{
+	if (faceCache.is_valid() && mesh->faceContainsPoint(point, faceCache)) {
+		return faceCache;
+	}
+	return (faceCache = mesh->findPoint(point, start));
+}
+
+OpenMesh::FaceHandle
+GeoDataMesh::PointFace::face(const GeoDataMesh * mesh) const
+{
+	return face(mesh, *mesh->faces_sbegin());
+}
+
+GeoDataMesh::HalfEdgeVertices
+GeoDataMesh::toVertexHandles(HalfedgeHandle halfEdge) const
+{
+	return {from_vertex_handle(halfEdge), to_vertex_handle(halfEdge)};
+}
+
+GeoDataMesh::HalfEdgePoints
+GeoDataMesh::points(HalfEdgeVertices vertices) const
+{
+	return {point(vertices.first), point(vertices.second)};
+}
+
+OpenMesh::FaceHandle
+GeoDataMesh::findPoint(const GlobalPosition2D coord, OpenMesh::FaceHandle face) const
+{
+	while (face.is_valid() && !triangle<2>(face).containsPoint(coord)) {
+		for (auto next = cfh_iter(face); next.is_valid(); ++next) {
+			face = opposite_face_handle(*next);
+			if (face.is_valid()) {
+				const auto nextPoints = points(toVertexHandles(*next));
+				if (pointLeftOfLine(coord, nextPoints.second, nextPoints.first)) {
+					break;
+				}
+			}
+			face.reset();
+		}
+	}
+	return face;
+}
+
+GlobalPosition3D
+GeoDataMesh::positionAt(const PointFace & coord) const
+{
+	return triangle<3>(coord.face(this)).positionOnPlane(coord.point);
+}
+
+bool
+GeoDataMesh::faceContainsPoint(const GlobalPosition2D coord, FaceHandle face) const
+{
+	return triangle<2>(face).containsPoint(coord);
+}
+
+OpenMesh::HalfedgeHandle
+GeoDataMesh::findBoundaryStart() const
+{
+	return *std::find_if(halfedges_sbegin(), halfedges_end(), [this](const auto heh) {
+		return is_boundary(heh);
+	});
+}
+
+[[nodiscard]] RelativePosition3D
+GeoDataMesh::difference(const HalfedgeHandle heh) const
+{
+	return ::difference(point(to_vertex_handle(heh)), point(from_vertex_handle(heh)));
+}
+
+[[nodiscard]] GlobalPosition3D
+GeoDataMesh::centre(const HalfedgeHandle heh) const
+{
+	const auto hehPoints = points(toVertexHandles(heh));
+	return midpoint(hehPoints.first, hehPoints.second);
+}
+
+#ifndef NDEBUG
+void
+GeoDataMesh::sanityCheck(const std::source_location & loc) const
+{
+	if (const auto upSideDown = std::ranges::count_if(faces(), [this](const auto face) {
+			if (!triangle<2>(face).isUp()) {
+				for (const auto vertex : fv_range(face)) {
+					CLOG(point(vertex));
+				}
+				return true;
+			}
+			return false;
+		}) > 0) {
+		throw std::logic_error(std::format(
+				"{} upside down faces detected - checked from {}:{}", upSideDown, loc.function_name(), loc.line()));
+	}
+}
+#endif
+
+bool
+GeoDataMesh::canFlip(const HalfedgeHandle edge) const
+{
+	const auto opposite = opposite_halfedge_handle(edge);
+	const auto pointA = point(to_vertex_handle(edge));
+	const auto pointB = point(to_vertex_handle(opposite));
+	const auto pointC = point(to_vertex_handle(next_halfedge_handle(edge)));
+	const auto pointD = point(to_vertex_handle(next_halfedge_handle(opposite)));
+
+	return Triangle<2> {pointC, pointB, pointD}.isUp() && Triangle<2> {pointA, pointC, pointD}.isUp();
+};
+
+std::optional<OpenMesh::EdgeHandle>
+GeoDataMesh::shouldFlip(const HalfedgeHandle next, const GlobalPosition2D startPoint) const
+{
+	if (const auto nextEdge = edge_handle(next); is_flip_ok(nextEdge) && canFlip(next)) {
+		const auto oppositePoint = point(to_vertex_handle(next_halfedge_handle(opposite_halfedge_handle(next)))).xy();
+		if (distance<2>(startPoint, oppositePoint) < length<2>(next)) {
+			return nextEdge;
+		}
+	}
+	return std::nullopt;
+};
+
+void
+GeoDataMesh::expandVerts(std::set<VertexHandle> & verts) const
+{
+	std::ranges::for_each(std::vector<VertexHandle>(verts.begin(), verts.end()), [&verts, this](auto vertex) {
+		std::ranges::copy(vv_range(vertex), std::inserter(verts, verts.end()));
+	});
+}
diff --git a/game/geoDataMesh.h b/game/geoDataMesh.h
new file mode 100644
index 0000000..72b069e
--- /dev/null
+++ b/game/geoDataMesh.h
@@ -0,0 +1,116 @@
+#pragma once
+
+#include "config/types.h"
+#include "triangle.h"
+#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
+#include <set>
+#include <source_location>
+#include <thirdparty/openmesh/glmcompat.h>
+#include <thirdparty/openmesh/helpers.h>
+
+struct GeoDataTraits : public OpenMesh::DefaultTraits {
+	FaceAttributes(OpenMesh::Attributes::Status);
+	EdgeAttributes(OpenMesh::Attributes::Status);
+	VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Status);
+	HalfedgeAttributes(OpenMesh::Attributes::Status);
+	using Point = GlobalPosition3D;
+	using Normal = Normal3D;
+};
+
+class GeoDataMesh : public OpenMesh::TriMesh_ArrayKernelT<GeoDataTraits> {
+public:
+	struct PointFace {
+		// NOLINTNEXTLINE(hicpp-explicit-conversions)
+		PointFace(GlobalPosition2D coord) : point {coord} { }
+
+		PointFace(GlobalPosition2D coord, FaceHandle face) : point {coord}, faceCache {face} { }
+
+		PointFace(GlobalPosition2D coord, const GeoDataMesh *);
+		PointFace(GlobalPosition2D coord, GeoDataMesh const *, FaceHandle start);
+
+		const GlobalPosition2D point;
+		[[nodiscard]] FaceHandle face(const GeoDataMesh *) const;
+		[[nodiscard]] FaceHandle face(const GeoDataMesh *, FaceHandle start) const;
+
+		[[nodiscard]] bool
+		isLocated() const
+		{
+			return faceCache.is_valid();
+		}
+
+	private:
+		mutable FaceHandle faceCache;
+	};
+
+	template<glm::length_t Dim> using Triangle = ::Triangle<Dim, GlobalDistance>;
+
+	[[nodiscard]] FaceHandle findPoint(GlobalPosition2D) const;
+	[[nodiscard]] FaceHandle findPoint(GlobalPosition2D, FaceHandle) const;
+
+	[[nodiscard]] GlobalPosition3D positionAt(const PointFace &) const;
+
+protected:
+#ifndef NDEBUG
+	void sanityCheck(const std::source_location & = std::source_location::current()) const;
+#endif
+
+	[[nodiscard]] bool faceContainsPoint(GlobalPosition2D, FaceHandle) const;
+	[[nodiscard]] HalfedgeHandle findBoundaryStart() const;
+	[[nodiscard]] RelativePosition3D difference(HalfedgeHandle) const;
+	using HalfEdgeVertices = std::pair<VertexHandle, VertexHandle>;
+	[[nodiscard]] HalfEdgeVertices toVertexHandles(HalfedgeHandle) const;
+	using HalfEdgePoints = std::pair<GlobalPosition3D, GlobalPosition3D>;
+	[[nodiscard]] HalfEdgePoints points(HalfEdgeVertices) const;
+
+	template<glm::length_t D>
+	[[nodiscard]] auto
+	vertexDistanceFunction(GlobalPosition<D> point) const
+	{
+		struct DistanceCalculator {
+			[[nodiscard]] std::pair<VertexHandle, float>
+			operator()(VertexHandle compVertex) const
+			{
+				return std::make_pair(
+						compVertex, ::distance<D, GlobalDistance, glm::defaultp>(point, mesh->point(compVertex)));
+			}
+
+			[[nodiscard]]
+			std::pair<HalfedgeHandle, float>
+			operator()(const HalfedgeHandle compHalfedge) const
+			{
+				const auto edgePoints = mesh->points(mesh->toVertexHandles(compHalfedge));
+				return std::make_pair(compHalfedge, Triangle<2> {edgePoints.second, edgePoints.first, point}.height());
+			};
+
+			const GeoDataMesh * mesh;
+			GlobalPosition<D> point;
+		};
+
+		return DistanceCalculator {this, point};
+	}
+
+	[[nodiscard]] bool canFlip(HalfedgeHandle edge) const;
+	[[nodiscard]] std::optional<EdgeHandle> shouldFlip(HalfedgeHandle next, GlobalPosition2D startPoint) const;
+	void expandVerts(std::set<VertexHandle> & verts) const;
+
+	template<glm::length_t D>
+	[[nodiscard]] RelativeDistance
+	length(HalfedgeHandle heh) const
+	{
+		return ::distance<D, GlobalDistance, glm::defaultp>(
+				point(to_vertex_handle(heh)), point(from_vertex_handle(heh)));
+	}
+
+	[[nodiscard]] GlobalPosition3D centre(HalfedgeHandle) const;
+
+	template<glm::length_t Dim>
+	[[nodiscard]] Triangle<Dim>
+	triangle(FaceHandle face) const
+	{
+		Triangle<Dim> triangle;
+		std::ranges::transform(fv_range(face), triangle.begin(), [this](auto vertex) {
+			return point(vertex);
+		});
+		return triangle;
+	}
+};
diff --git a/game/network/link.cpp b/game/network/link.cpp
index 79af92a..c84524c 100644
--- a/game/network/link.cpp
+++ b/game/network/link.cpp
@@ -38,6 +38,20 @@ LinkStraight::intersectRay(const Ray<GlobalPosition3D> & ray) const
 			std::array {GlobalPosition3D {ends.front().node->pos}, GlobalPosition3D {ends.back().node->pos}}, 1000);
 }
 
+std::vector<GlobalPosition3D>
+LinkStraight::getBase(RelativeDistance width) const
+{
+	const auto start = ends.front().node->pos;
+	const auto end = ends.back().node->pos;
+	const auto direction = (vector_normal(normalize(::difference(start, end).xy())) * width / 2.F) || 0.F;
+	return {
+			start - direction,
+			start + direction,
+			end - direction,
+			end + direction,
+	};
+}
+
 Location
 LinkCurve::positionAt(float dist, unsigned char start) const
 {
@@ -73,3 +87,24 @@ LinkCurve::intersectRay(const Ray<GlobalPosition3D> & ray) const
 	}
 	return ray.passesCloseToEdges(points, 1.F);
 }
+
+std::vector<GlobalPosition3D>
+LinkCurve::getBase(RelativeDistance width) const
+{
+	const auto start = ends.front().node->pos;
+	const auto end = ends.back().node->pos;
+	const auto segs = std::ceil(std::sqrt(radius) * 0.02F * arc.length());
+	const auto step {glm::vec<2, RelativeDistance> {arc.length(), end.z - start.z} / segs};
+
+	auto segCount = static_cast<size_t>(segs) + 1;
+	std::vector<GlobalPosition3D> out;
+	out.reserve(segCount);
+	for (RelativePosition2D swing = {arc.first, centreBase.z - start.z}; segCount != 0U; swing += step, --segCount) {
+		const auto direction = sincos(swing.x);
+		const auto linkCentre = centreBase + ((direction * radius) || swing.y);
+		const auto toEdge = (direction * width / 2.F) || 0.F;
+		out.emplace_back(linkCentre + toEdge);
+		out.emplace_back(linkCentre - toEdge);
+	}
+	return out;
+}
diff --git a/game/network/link.h b/game/network/link.h
index 725e023..59bbb65 100644
--- a/game/network/link.h
+++ b/game/network/link.h
@@ -16,7 +16,7 @@ template<typename> class Ray;
 // it has location
 class Node : public StdTypeDefs<Node> {
 public:
-	explicit Node(GlobalPosition3D p) noexcept : pos(p) {};
+	explicit Node(GlobalPosition3D p) noexcept : pos(p) { };
 	virtual ~Node() noexcept = default;
 	NO_COPY(Node);
 	NO_MOVE(Node);
@@ -38,13 +38,14 @@ public:
 		Nexts nexts {};
 	};
 
-	Link(End, End, float);
+	Link(End, End, RelativeDistance length);
 	virtual ~Link() = default;
 	NO_COPY(Link);
 	NO_MOVE(Link);
 
 	[[nodiscard]] virtual Location positionAt(RelativeDistance dist, unsigned char start) const = 0;
 	[[nodiscard]] virtual bool intersectRay(const Ray<GlobalPosition3D> &) const = 0;
+	[[nodiscard]] virtual std::vector<GlobalPosition3D> getBase(RelativeDistance width) const = 0;
 
 	std::array<End, 2> ends;
 	float length;
@@ -69,6 +70,7 @@ public:
 
 	[[nodiscard]] Location positionAt(RelativeDistance dist, unsigned char start) const override;
 	[[nodiscard]] bool intersectRay(const Ray<GlobalPosition3D> &) const override;
+	[[nodiscard]] std::vector<GlobalPosition3D> getBase(RelativeDistance width) const override;
 };
 
 LinkStraight::~LinkStraight() = default;
@@ -76,12 +78,13 @@ LinkStraight::~LinkStraight() = default;
 class LinkCurve : public virtual Link {
 public:
 	inline ~LinkCurve() override = 0;
-	LinkCurve(GlobalPosition3D, RelativeDistance, Arc);
+	LinkCurve(GlobalPosition3D centreBase, RelativeDistance radius, Arc);
 	NO_COPY(LinkCurve);
 	NO_MOVE(LinkCurve);
 
 	[[nodiscard]] Location positionAt(RelativeDistance dist, unsigned char start) const override;
 	[[nodiscard]] bool intersectRay(const Ray<GlobalPosition3D> &) const override;
+	[[nodiscard]] std::vector<GlobalPosition3D> getBase(RelativeDistance width) const override;
 
 	GlobalPosition3D centreBase;
 	RelativeDistance radius;
diff --git a/game/network/network.cpp b/game/network/network.cpp
index 6ba3ed6..e67942f 100644
--- a/game/network/network.cpp
+++ b/game/network/network.cpp
@@ -8,7 +8,13 @@
 #include <stdexcept>
 #include <utility>
 
-Network::Network(const std::string & tn) : texture {std::make_shared<Texture>(tn)} { }
+Network::Network(const std::string & tn) :
+	texture {std::make_shared<Texture>(tn,
+			TextureOptions {
+					.minFilter = GL_NEAREST_MIPMAP_LINEAR,
+			})}
+{
+}
 
 Node::Ptr
 Network::nodeAt(GlobalPosition3D pos)
@@ -115,8 +121,8 @@ Network::genCurveDef(const GlobalPosition3D & start, const GlobalPosition3D & en
 	endDir += pi;
 	const auto flatStart {start.xy()}, flatEnd {end.xy()};
 	auto midheight = [&](auto mid) {
-		const auto sm = glm::length(RelativePosition2D(flatStart - mid)),
-				   em = glm::length(RelativePosition2D(flatEnd - mid));
+		const auto sm = ::distance<2>(flatStart, mid);
+		const auto em = ::distance<2>(flatEnd, mid);
 		return start.z + GlobalDistance(RelativeDistance(end.z - start.z) * (sm / (sm + em)));
 	};
 	if (const auto radii = find_arcs_radius(flatStart, startDir, flatEnd, endDir); radii.first < radii.second) {
diff --git a/game/network/network.h b/game/network/network.h
index be0900b..291c4ec 100644
--- a/game/network/network.h
+++ b/game/network/network.h
@@ -14,6 +14,8 @@
 #include <utility>
 
 class SceneShader;
+struct Surface;
+class GeoData;
 template<typename> class Ray;
 
 template<size_t... n> using GenDef = std::tuple<glm::vec<n, GlobalDistance>...>;
@@ -41,12 +43,15 @@ public:
 	virtual Link::CCollection candidateStraight(GlobalPosition3D, GlobalPosition3D) = 0;
 	virtual Link::CCollection candidateJoins(GlobalPosition3D, GlobalPosition3D) = 0;
 	virtual Link::CCollection candidateExtend(GlobalPosition3D, GlobalPosition3D) = 0;
-	virtual Link::CCollection addStraight(GlobalPosition3D, GlobalPosition3D) = 0;
-	virtual Link::CCollection addJoins(GlobalPosition3D, GlobalPosition3D) = 0;
-	virtual Link::CCollection addExtend(GlobalPosition3D, GlobalPosition3D) = 0;
+	virtual Link::CCollection addStraight(const GeoData *, GlobalPosition3D, GlobalPosition3D) = 0;
+	virtual Link::CCollection addJoins(const GeoData *, GlobalPosition3D, GlobalPosition3D) = 0;
+	virtual Link::CCollection addExtend(const GeoData *, GlobalPosition3D, GlobalPosition3D) = 0;
 
 	[[nodiscard]] virtual float findNodeDirection(Node::AnyCPtr) const = 0;
 
+	[[nodiscard]] virtual const Surface * getBaseSurface() const = 0;
+	[[nodiscard]] virtual RelativeDistance getBaseWidth() const = 0;
+
 protected:
 	static void joinLinks(const Link::Ptr & l, const Link::Ptr & ol);
 	static GenCurveDef genCurveDef(const GlobalPosition3D & start, const GlobalPosition3D & end, float startDir);
@@ -102,12 +107,12 @@ public:
 	Link::CCollection candidateStraight(GlobalPosition3D n1, GlobalPosition3D n2) override;
 	Link::CCollection candidateJoins(GlobalPosition3D, GlobalPosition3D) override;
 	Link::CCollection candidateExtend(GlobalPosition3D, GlobalPosition3D) override;
-	Link::CCollection addStraight(GlobalPosition3D n1, GlobalPosition3D n2) override;
-	Link::CCollection addJoins(GlobalPosition3D, GlobalPosition3D) override;
-	Link::CCollection addExtend(GlobalPosition3D, GlobalPosition3D) override;
+	Link::CCollection addStraight(const GeoData *, GlobalPosition3D n1, GlobalPosition3D n2) override;
+	Link::CCollection addJoins(const GeoData *, GlobalPosition3D, GlobalPosition3D) override;
+	Link::CCollection addExtend(const GeoData *, GlobalPosition3D, GlobalPosition3D) override;
 
 	[[nodiscard]] float findNodeDirection(Node::AnyCPtr) const override;
 
 protected:
-	Link::CCollection addJoins();
+	Link::CCollection addCurve(const GeoData *, const GenCurveDef &);
 };
diff --git a/game/network/network.impl.h b/game/network/network.impl.h
index ff29088..c683378 100644
--- a/game/network/network.impl.h
+++ b/game/network/network.impl.h
@@ -1,4 +1,6 @@
+#include "collections.h"
 #include "network.h"
+#include <game/geoData.h>
 #include <gfx/gl/sceneShader.h>
 #include <gfx/models/texture.h>
 
@@ -52,7 +54,7 @@ template<typename T, typename... Links>
 Link::CCollection
 NetworkOf<T, Links...>::candidateJoins(GlobalPosition3D start, GlobalPosition3D end)
 {
-	if (glm::length(RelativePosition3D(start - end)) < 2000.F) {
+	if (::distance(start, end) < 2000.F) {
 		return {};
 	}
 	const auto defs = genCurveDef(
@@ -72,28 +74,69 @@ NetworkOf<T, Links...>::candidateExtend(GlobalPosition3D start, GlobalPosition3D
 
 template<typename T, typename... Links>
 Link::CCollection
-NetworkOf<T, Links...>::addStraight(GlobalPosition3D n1, GlobalPosition3D n2)
+NetworkOf<T, Links...>::addStraight(const GeoData * geoData, GlobalPosition3D n1, GlobalPosition3D n2)
 {
-	return {addLink<typename T::StraightLink>(n1, n2)};
+	Link::CCollection out;
+	geoData->walk(n1.xy(), n2, [geoData, &out, this, &n1](const GeoData::WalkStep & step) {
+		if (step.previous.is_valid() && geoData->getSurface(step.current) != geoData->getSurface(step.previous)) {
+			const auto surfaceEdgePosition = geoData->positionAt(GeoData::PointFace(step.exitPosition, step.current));
+			out.emplace_back(addLink<typename T::StraightLink>(n1, surfaceEdgePosition));
+			n1 = surfaceEdgePosition;
+		}
+	});
+	out.emplace_back(addLink<typename T::StraightLink>(n1, n2));
+	return out;
+}
+
+template<typename T, typename... Links>
+Link::CCollection
+NetworkOf<T, Links...>::addCurve(const GeoData * geoData, const GenCurveDef & curve)
+{
+	auto [cstart, cend, centre] = curve;
+	Link::CCollection out;
+	std::set<GeoData::WalkStepCurve, SortedBy<&GeoData::WalkStepCurve::angle>> breaks;
+	const auto radiusMid = ::distance(cstart.xy(), centre);
+	for (const auto radiusOffset : {-getBaseWidth() / 2.F, 0.F, getBaseWidth() / 2.F}) {
+		const auto radius = radiusOffset + radiusMid;
+		const auto start = centre + (difference(cstart.xy(), centre) * radius) / radiusMid;
+		const auto end = centre + (difference(cend.xy(), centre) * radius) / radiusMid;
+		geoData->walk(start, end, centre, [geoData, &breaks](const GeoData::WalkStepCurve & step) {
+			if (step.previous.is_valid() && geoData->getSurface(step.current) != geoData->getSurface(step.previous)) {
+				breaks.insert(step);
+			}
+		});
+	}
+	std::vector<GlobalPosition3D> points;
+	points.reserve(breaks.size() + 2);
+	points.push_back(cstart);
+	std::ranges::transform(
+			breaks, std::back_inserter(points), [geoData, centre, radiusMid](const GeoData::WalkStepCurve & step) {
+				return (centre + (sincos(step.angle) * radiusMid))
+						|| geoData->positionAt(GeoData::PointFace(step.exitPosition, step.current)).z;
+			});
+	points.push_back(cend);
+	mergeClose(points, ::distance<3, GlobalDistance>, ::midpoint<3, GlobalDistance>, 2'000.F);
+	std::ranges::transform(points | std::views::pairwise, std::back_inserter(out), [this, centre](const auto pair) {
+		const auto [a, b] = pair;
+		return addLink<typename T::CurveLink>(a, b, centre);
+	});
+	return out;
 }
 
 template<typename T, typename... Links>
 Link::CCollection
-NetworkOf<T, Links...>::addJoins(GlobalPosition3D start, GlobalPosition3D end)
+NetworkOf<T, Links...>::addJoins(const GeoData * geoData, GlobalPosition3D start, GlobalPosition3D end)
 {
-	if (glm::length(RelativePosition3D(start - end)) < 2000.F) {
+	if (::distance(start, end) < 2000.F) {
 		return {};
 	}
 	const auto defs = genCurveDef(start, end, findNodeDirection(nodeAt(start)), findNodeDirection(nodeAt(end)));
-	const auto & [c1s, c1e, c1c] = defs.first;
-	const auto & [c2s, c2e, c2c] = defs.second;
-	return {addLink<typename T::CurveLink>(c1s, c1e, c1c), addLink<typename T::CurveLink>(c2s, c2e, c2c)};
+	return addCurve(geoData, defs.first) + addCurve(geoData, defs.second);
 }
 
 template<typename T, typename... Links>
 Link::CCollection
-NetworkOf<T, Links...>::addExtend(GlobalPosition3D start, GlobalPosition3D end)
+NetworkOf<T, Links...>::addExtend(const GeoData * geoData, GlobalPosition3D start, GlobalPosition3D end)
 {
-	const auto [cstart, cend, centre] = genCurveDef(start, end, findNodeDirection(nodeAt(start)));
-	return {addLink<typename T::CurveLink>(cstart, cend, centre)};
+	return addCurve(geoData, genCurveDef(start, end, findNodeDirection(nodeAt(start))));
 }
diff --git a/game/network/rail.cpp b/game/network/rail.cpp
index 6f04070..d7de231 100644
--- a/game/network/rail.cpp
+++ b/game/network/rail.cpp
@@ -1,4 +1,6 @@
 #include "rail.h"
+#include "game/gamestate.h"
+#include "game/geoData.h"
 #include "network.h"
 #include <game/network/network.impl.h> // IWYU pragma: keep
 #include <gfx/gl/sceneShader.h>
@@ -8,7 +10,7 @@
 template class NetworkOf<RailLink, RailLinkStraight, RailLinkCurve>;
 
 constexpr auto RAIL_CROSSSECTION_VERTICES {5U};
-constexpr Size3D RAIL_HEIGHT {0, 0, 250.F};
+constexpr Size3D RAIL_HEIGHT {0, 0, 50.F};
 
 RailLinks::RailLinks() : NetworkOf<RailLink, RailLinkStraight, RailLinkCurve> {"rails.jpg"} { }
 
@@ -38,8 +40,8 @@ RailLinks::addLinksBetween(GlobalPosition3D start, GlobalPosition3D end)
 	const auto flatStart {start.xy()}, flatEnd {end.xy()};
 	if (node2ins.second == NodeIs::InNetwork) {
 		auto midheight = [&](auto mid) {
-			const auto sm = glm::length(RelativePosition2D(flatStart - mid)),
-					   em = glm::length(RelativePosition2D(flatEnd - mid));
+			const auto sm = ::distance<2>(flatStart, mid);
+			const auto em = ::distance<2>(flatEnd, mid);
 			return start.z + GlobalDistance(RelativeDistance(end.z - start.z) * (sm / (sm + em)));
 		};
 		const float dir2 = pi + findNodeDirection(node2ins.first);
@@ -73,26 +75,31 @@ RailLinks::addLinksBetween(GlobalPosition3D start, GlobalPosition3D end)
 	return addLink<RailLinkCurve>(start, end, centre.first);
 }
 
-constexpr const std::array<RelativePosition3D, RAIL_CROSSSECTION_VERTICES> railCrossSection {{
-		{-1900.F, 0.F, 0.F},
-		{-608.F, 0.F, RAIL_HEIGHT.z},
-		{0, 0.F, RAIL_HEIGHT.z * .7F},
-		{608.F, 0.F, RAIL_HEIGHT.z},
-		{1900.F, 0.F, 0.F},
-}};
-constexpr const std::array<float, RAIL_CROSSSECTION_VERTICES> railTexturePos {
-		0.F,
-		.34F,
-		.5F,
-		.66F,
-		1.F,
-};
-constexpr auto sleepers {5.F}; // There are 5 repetitions of sleepers in the texture
-
-inline auto
-round_sleepers(const float v)
-{
-	return round_frac(v, sleepers);
+namespace {
+	constexpr const std::array<RelativePosition3D, RAIL_CROSSSECTION_VERTICES> RAIL_CROSS_SECTION {{
+			{-1330.F, 0.F, 0},
+			{-608.F, 0.F, RAIL_HEIGHT.z},
+			{0, 0.F, RAIL_HEIGHT.z / 2},
+			{608.F, 0.F, RAIL_HEIGHT.z},
+			{1330.F, 0.F, 0},
+	}};
+	constexpr const std::array<float, RAIL_CROSSSECTION_VERTICES> RAIL_TEXTURE_POS {
+			0.15F,
+			.34F,
+			.5F,
+			.66F,
+			0.85F,
+	};
+	template<std::floating_point T> constexpr T SLEEPERS_PER_TEXTURE {5};
+	template<std::floating_point T> constexpr T TEXTURE_LENGTH {2'000};
+	template<std::floating_point T> constexpr T SLEEPER_LENGTH {T {1} / SLEEPERS_PER_TEXTURE<T>};
+
+	template<std::floating_point T>
+	constexpr auto
+	roundSleepers(const T length)
+	{
+		return round_frac(length / TEXTURE_LENGTH<T>, SLEEPER_LENGTH<T>);
+	}
 }
 
 RailLinkStraight::RailLinkStraight(NetworkLinkHolder<RailLinkStraight> & instances, const Node::Ptr & a,
@@ -104,22 +111,22 @@ RailLinkStraight::RailLinkStraight(
 		NetworkLinkHolder<RailLinkStraight> & instances, Node::Ptr a, Node::Ptr b, const RelativePosition3D & diff) :
 	Link({std::move(a), vector_yaw(diff)}, {std::move(b), vector_yaw(-diff)}, glm::length(diff)),
 	instance {instances.vertices.acquire(
-			ends[0].node->pos, ends[1].node->pos, flat_orientation(diff), round_sleepers(length / 2000.F))}
+			ends[0].node->pos, ends[1].node->pos, flat_orientation(diff), roundSleepers(length))}
 {
 }
 
-RailLinkCurve::RailLinkCurve(NetworkLinkHolder<RailLinkCurve> & instances, const Node::Ptr & a, const Node::Ptr & b,
-		GlobalPosition2D c) : RailLinkCurve(instances, a, b, c || a->pos.z, {c, a->pos, b->pos})
+RailLinkCurve::RailLinkCurve(
+		NetworkLinkHolder<RailLinkCurve> & instances, const Node::Ptr & a, const Node::Ptr & b, GlobalPosition2D c) :
+	RailLinkCurve(instances, a, b, c || a->pos.z, ::distance<2>(a->pos.xy(), c), {c, a->pos, b->pos})
 {
 }
 
 RailLinkCurve::RailLinkCurve(NetworkLinkHolder<RailLinkCurve> & instances, const Node::Ptr & a, const Node::Ptr & b,
-		GlobalPosition3D c, const Arc arc) :
+		GlobalPosition3D c, RelativeDistance radius, const Arc arc) :
 	Link({a, normalize(arc.first + half_pi)}, {b, normalize(arc.second - half_pi)},
-			glm::length(RelativePosition3D(a->pos - c)) * arc.length()),
-	LinkCurve {c, glm::length(RelativePosition3D(ends[0].node->pos - c)), arc},
-	instance {instances.vertices.acquire(ends[0].node->pos, ends[1].node->pos, c, round_sleepers(length / 2000.F),
-			half_pi - arc.first, half_pi - arc.second, radius)}
+			glm::length(RelativePosition2D {radius * arc.length(), difference(a->pos, b->pos).z})),
+	LinkCurve {c, radius, arc}, instance {instances.vertices.acquire(ends[0].node->pos, ends[1].node->pos, c,
+										roundSleepers(length), half_pi - arc.first, half_pi - arc.second, radius)}
 {
 }
 
@@ -151,7 +158,7 @@ namespace {
 	renderType(const NetworkLinkHolder<LinkType> & n, auto & s)
 	{
 		if (auto count = n.vertices.size()) {
-			s.use(railCrossSection, railTexturePos);
+			s.use(RAIL_CROSS_SECTION, RAIL_TEXTURE_POS);
 			glBindVertexArray(n.vao);
 			glDrawArrays(GL_POINTS, 0, static_cast<GLsizei>(count));
 		}
@@ -163,8 +170,23 @@ RailLinks::render(const SceneShader & shader) const
 {
 	if (!links.objects.empty()) {
 		texture->bind();
+		glEnable(GL_POLYGON_OFFSET_FILL);
+		glPolygonOffset(-1, 0);
 		renderType<RailLinkStraight>(*this, shader.networkStraight);
 		renderType<RailLinkCurve>(*this, shader.networkCurve);
+		glDisable(GL_POLYGON_OFFSET_FILL);
 		glBindVertexArray(0);
 	}
 }
+
+const Surface *
+RailLinks::getBaseSurface() const
+{
+	return std::dynamic_pointer_cast<Surface>(gameState->assets.at("terrain.surface.gravel")).get();
+}
+
+RelativeDistance
+RailLinks::getBaseWidth() const
+{
+	return 5'700;
+}
diff --git a/game/network/rail.h b/game/network/rail.h
index c8effef..fa64eda 100644
--- a/game/network/rail.h
+++ b/game/network/rail.h
@@ -62,8 +62,8 @@ public:
 	};
 
 private:
-	RailLinkCurve(
-			NetworkLinkHolder<RailLinkCurve> &, const Node::Ptr &, const Node::Ptr &, GlobalPosition3D, const Arc);
+	RailLinkCurve(NetworkLinkHolder<RailLinkCurve> &, const Node::Ptr &, const Node::Ptr &, GlobalPosition3D centreBase,
+			RelativeDistance radius, Arc);
 	InstanceVertices<Vertex>::InstanceProxy instance;
 };
 
@@ -77,6 +77,9 @@ public:
 	std::shared_ptr<RailLink> addLinksBetween(GlobalPosition3D start, GlobalPosition3D end);
 	void render(const SceneShader &) const override;
 
+	[[nodiscard]] const Surface * getBaseSurface() const override;
+	[[nodiscard]] RelativeDistance getBaseWidth() const override;
+
 private:
 	void tick(TickDuration elapsed) override;
 };
diff --git a/game/terrain.cpp b/game/terrain.cpp
index bb8e3ce..01af163 100644
--- a/game/terrain.cpp
+++ b/game/terrain.cpp
@@ -1,6 +1,5 @@
 #include "terrain.h"
 #include "game/geoData.h"
-#include "gfx/models/texture.h"
 #include <algorithm>
 #include <cstddef>
 #include <gfx/gl/sceneShader.h>
@@ -15,12 +14,7 @@
 #include <utility>
 #include <vector>
 
-static constexpr RGB openSurface {-1};
-
-Terrain::Terrain(std::shared_ptr<GeoData> tm) : geoData {std::move(tm)}, grass {std::make_shared<Texture>("grass.png")}
-{
-	generateMeshes();
-}
+static constexpr RGB OPEN_SURFACE {-1};
 
 template<>
 VertexArrayObject &
@@ -35,31 +29,27 @@ Terrain::generateMeshes()
 {
 	meshes.removeAll();
 	std::vector<unsigned int> indices;
-	indices.reserve(geoData->n_faces() * 3);
+	indices.reserve(n_faces() * 3);
 	std::vector<Vertex> vertices;
-	vertices.reserve(geoData->n_vertices());
-	std::map<std::pair<GeoData::VertexHandle, const Surface *>, size_t> vertexIndex;
-	std::for_each(geoData->vertices_sbegin(), geoData->vertices_end(),
-			[this, &vertexIndex, &vertices](const GeoData::VertexHandle v) {
-				std::for_each(geoData->vf_begin(v), geoData->vf_end(v),
-						[&vertexIndex, v, this, &vertices](const GeoData::FaceHandle f) {
-							const auto surface = geoData->get_surface(f);
-							if (const auto vertexIndexRef = vertexIndex.emplace(std::make_pair(v, surface), 0);
-									vertexIndexRef.second) {
-								vertexIndexRef.first->second = vertices.size();
+	vertices.reserve(n_vertices());
+	std::map<std::pair<VertexHandle, const Surface *>, size_t> vertexIndex;
+	std::ranges::for_each(this->vertices(), [this, &vertexIndex, &vertices](const auto vertex) {
+		std::ranges::for_each(vf_range(vertex), [&vertexIndex, vertex, this, &vertices](const auto face) {
+			const auto * const surface = getSurface(face);
+			if (const auto vertexIndexRef = vertexIndex.emplace(std::make_pair(vertex, surface), 0);
+					vertexIndexRef.second) {
+				vertexIndexRef.first->second = vertices.size();
 
-								vertices.emplace_back(geoData->point(v), geoData->normal(v),
-										surface ? surface->colorBias : openSurface);
-							}
-						});
-			});
-	std::for_each(
-			geoData->faces_sbegin(), geoData->faces_end(), [this, &vertexIndex, &indices](const GeoData::FaceHandle f) {
-				std::transform(geoData->fv_begin(f), geoData->fv_end(f), std::back_inserter(indices),
-						[&vertexIndex, f, this](const GeoData::VertexHandle v) {
-							return vertexIndex[std::make_pair(v, geoData->get_surface(f))];
-						});
-			});
+				vertices.emplace_back(point(vertex), normal(vertex), surface ? surface->colorBias : OPEN_SURFACE);
+			}
+		});
+	});
+	std::ranges::for_each(faces(), [this, &vertexIndex, &indices](const auto face) {
+		std::ranges::transform(
+				fv_range(face), std::back_inserter(indices), [&vertexIndex, face, this](const auto vertex) {
+					return vertexIndex[std::make_pair(vertex, getSurface(face))];
+				});
+	});
 	meshes.create<MeshT<Vertex>>(vertices, indices);
 }
 
@@ -69,6 +59,12 @@ Terrain::tick(TickDuration)
 }
 
 void
+Terrain::afterChange()
+{
+		generateMeshes();
+}
+
+void
 Terrain::render(const SceneShader & shader) const
 {
 	shader.landmass.use();
diff --git a/game/terrain.h b/game/terrain.h
index 7d074cf..f0f9621 100644
--- a/game/terrain.h
+++ b/game/terrain.h
@@ -4,17 +4,19 @@
 #include "collection.h"
 #include "config/types.h"
 #include "game/worldobject.h"
+#include "geoData.h"
 #include "gfx/models/mesh.h"
 #include "gfx/models/texture.h"
 #include "gfx/renderable.h"
-#include <memory>
 
 class SceneShader;
-class GeoData;
 
-class Terrain : public WorldObject, public Renderable {
+class Terrain : public GeoData, public WorldObject, public Renderable {
 public:
-	explicit Terrain(std::shared_ptr<GeoData>);
+	template<typename... P> explicit Terrain(P &&... params) : GeoData {std::forward<P>(params)...}
+	{
+		generateMeshes();
+	}
 
 	void render(const SceneShader & shader) const override;
 	void shadows(const ShadowMapper &) const override;
@@ -27,10 +29,11 @@ public:
 		RGB colourBias;
 	};
 
+private:
+	void afterChange() override;
 	void generateMeshes();
 
-private:
-	std::shared_ptr<GeoData> geoData;
 	Collection<MeshT<Vertex>, false> meshes;
-	Texture::Ptr grass;
+	Texture::Ptr grass = std::make_shared<Texture>("grass.png");
+	size_t geoGeneration {};
 };