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

44 files changed, 1690 insertions, 667 deletions

diff --git a/application/main.cpp b/application/main.cpp
index db42a63..723f3d2 100644
--- a/application/main.cpp
+++ b/application/main.cpp
@@ -11,6 +11,7 @@
 #include <game/network/link.h>
 #include <game/network/rail.h>
 #include <game/objective.h>
+#include <game/objectives/freeroam.h>
 #include <game/objectives/goto.h>
 #include <game/orders.h>
 #include <game/scenary/foliage.h>
@@ -25,8 +26,13 @@
 #include <glm/gtx/transform.hpp> // IWYU pragma: keep
 #include <memory>
 #include <random>
+#include <ranges>
 #include <special_members.h>
+#include <stream_support.h>
 #include <ui/applicationBase.h>
+#include <ui/builders/freeExtend.h>
+#include <ui/builders/join.h>
+#include <ui/builders/straight.h>
 #include <ui/gameMainWindow.h>
 #include <ui/window.h>
 
@@ -38,46 +44,54 @@ public:
 	int
 	run()
 	{
-		geoData = std::make_shared<GeoData>(GeoData::loadFromAsciiGrid("test/fixtures/height/SD19.asc"));
-
 		windows.create<MainWindow>(DISPLAY_WIDTH, DISPLAY_HEIGHT)->setContent<GameMainWindow>();
 
-		world.create<Terrain>(geoData);
-		world.create<Water>(geoData);
+		terrain = world.create<Terrain>(GeoData::loadFromAsciiGrid("test/fixtures/height/SD19.asc"));
+		world.create<Water>(terrain);
 		assets = AssetFactory::loadAll("res");
 
 		{
 			auto rl = world.create<RailLinks>();
-			const GlobalPosition3D j {-1120000, -1100000, 3000}, k {-1100000, -1000000, 15000},
-					l {-1000000, -800000, 20000}, m {-900000, -600000, 30000}, n {-600000, -500000, 32000},
-					o {-500000, -800000, 30000}, p {-600000, -900000, 25000}, q {-1025000, -1175000, 10000},
-					r {-925000, -1075000, 10000}, s {-1100000, -500000, 15000}, t {-1100000, -450000, 15000},
-					u {-1000000, -400000, 15000};
-			auto l3 = rl->addLinksBetween(j, k);
-			rl->addLinksBetween(k, l);
-			rl->addLinksBetween(l, m);
-			rl->addLinksBetween(m, n);
-			rl->addLinksBetween(n, o);
-			rl->addLinksBetween(o, p);
-			// branch 1
-			rl->addLinksBetween(p, q);
-			rl->addLinksBetween(q, j);
-			// branch 2
-			rl->addLinksBetween(p, r);
-			rl->addLinksBetween(r, j);
-			// early loop
-			rl->addLinksBetween(s, t);
-			rl->addLinksBetween(l, s);
-			rl->addLinksBetween(t, u);
-			rl->addLinksBetween(u, m);
-			const std::shared_ptr<Train> train = world.create<Train>(l3);
+			const auto nodes = materializeRange(std::vector<GlobalPosition2D> {
+														{315103000, 491067000},
+														{315977000, 490777000},
+														{316312000, 490557000},
+														{316885000, 491330000},
+														{316510934, 491255979},
+														{316129566, 490893054},
+														{315825622, 490833929},
+														{315106182, 491073714},
+												}
+					| std::views::transform([this](const auto n) {
+						  return terrain->positionAt(n);
+					  }));
+			auto l3 = BuilderStraight {}.create(rl.get(), terrain.get(), *nodes.begin(), *++nodes.begin()).front();
+			for (const auto & [from, to] : nodes | std::views::drop(1) | std::views::pairwise) {
+				const auto links = BuilderFreeExtend {}.createExtend(rl.get(), terrain.get(), from, to);
+			}
+			for (const auto & [from, to] : std::initializer_list<std::pair<GlobalPosition2D, GlobalPosition2D>> {
+						 {{315103000, 491067000}, {315003434, 491076253}},
+						 {{315103000, 491067000}, {315016495, 491019224}},
+						 {{315016495, 491019224}, {314955393, 490999023}},
+				 }) {
+				const auto links = BuilderFreeExtend {}.createExtend(
+						rl.get(), terrain.get(), terrain->positionAt(from), terrain->positionAt(to));
+			}
+			for (const auto & [from, to] : std::initializer_list<std::pair<GlobalPosition2D, GlobalPosition2D>> {
+						 {{315106182, 491073714}, {314955393, 490999023}},
+				 }) {
+				auto p1 = rl->intersectRayNodes({from || 0, up})->pos;
+				auto p2 = rl->intersectRayNodes({to || 0, up})->pos;
+				const auto links = BuilderFreeExtend {}.createJoin(rl.get(), terrain.get(), p1, p2);
+			}
+
+			const std::shared_ptr<Train> train = world.create<Train>(l3, 800000);
 			auto b47 = std::dynamic_pointer_cast<RailVehicleClass>(assets.at("brush-47"));
 			for (int N = 0; N < 6; N++) {
 				train->create<RailVehicle>(b47);
 			}
 			train->orders.removeAll();
-			train->orders.create<GoTo>(
-					&train->orders, l3->ends[1], l3->length, rl->findNodeAt({-1100000, -450000, 15000}));
+			train->orders.create<FreeRoam>(&train->orders);
 			train->currentActivity = train->orders.current()->createActivity();
 
 			std::random_device randomdev {};
@@ -89,7 +103,7 @@ public:
 				for (auto y = 491100000; y < 491130000; y += 5000) {
 					world.create<Plant>(std::dynamic_pointer_cast<Foliage>(assets.at(std::format("Tree-{:#02}-{}",
 												treeDistribution(randomdev), treeVariantDistribution(randomdev)))),
-							Location {geoData->positionAt({{x + positionOffsetDistribution(randomdev),
+							Location {terrain->positionAt({{x + positionOffsetDistribution(randomdev),
 											  y + positionOffsetDistribution(randomdev)}}),
 									{0, rotationDistribution(randomdev), 0}});
 				}
diff --git a/assetFactory/modelFactoryMesh.cpp b/assetFactory/modelFactoryMesh.cpp
index 3d4b5f3..3660fb7 100644
--- a/assetFactory/modelFactoryMesh.cpp
+++ b/assetFactory/modelFactoryMesh.cpp
@@ -1,13 +1,5 @@
 #include "modelFactoryMesh.h"
 
-ModelFactoryMesh::ModelFactoryMesh()
-{
-	add_property(smoothFaceProperty);
-	add_property(materialFaceProperty);
-	add_property(nameFaceProperty);
-	add_property(nameAdjFaceProperty);
-}
-
 void
 ModelFactoryMesh::configNamedFace(const std::string & name, OpenMesh::FaceHandle handle)
 {
diff --git a/assetFactory/modelFactoryMesh.h b/assetFactory/modelFactoryMesh.h
index 299986e..6a18155 100644
--- a/assetFactory/modelFactoryMesh.h
+++ b/assetFactory/modelFactoryMesh.h
@@ -8,6 +8,7 @@
 #include <glm/vec3.hpp>
 #include <glm/vec4.hpp>
 #include <thirdparty/openmesh/glmcompat.h>
+#include <thirdparty/openmesh/helpers.h>
 
 struct ModelFactoryTraits : public OpenMesh::DefaultTraits {
 	FaceAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Status | OpenMesh::Attributes::Color);
@@ -21,13 +22,11 @@ struct ModelFactoryTraits : public OpenMesh::DefaultTraits {
 };
 
 struct ModelFactoryMesh : public OpenMesh::PolyMesh_ArrayKernelT<ModelFactoryTraits> {
-	ModelFactoryMesh();
-
 	bool normalsProvidedProperty {};
-	OpenMesh::FPropHandleT<bool> smoothFaceProperty;
-	OpenMesh::FPropHandleT<GLuint> materialFaceProperty;
-	OpenMesh::FPropHandleT<std::string> nameFaceProperty;
-	OpenMesh::HPropHandleT<std::string> nameAdjFaceProperty;
+	const OpenMesh::Helpers::Property<bool, OpenMesh::FPropHandleT> smoothFaceProperty {this};
+	const OpenMesh::Helpers::Property<GLuint, OpenMesh::FPropHandleT> materialFaceProperty {this};
+	const OpenMesh::Helpers::Property<std::string, OpenMesh::FPropHandleT> nameFaceProperty {this};
+	const OpenMesh::Helpers::Property<std::string, OpenMesh::HPropHandleT> nameAdjFaceProperty {this};
 
 	template<typename... Vs>
 	std::pair<std::string, OpenMesh::FaceHandle>
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 {};
 };
diff --git a/gfx/models/texture.cpp b/gfx/models/texture.cpp
index a508421..3457fb5 100644
--- a/gfx/models/texture.cpp
+++ b/gfx/models/texture.cpp
@@ -50,6 +50,16 @@ Texture::Texture(GLsizei width, GLsizei height, const void * data, TextureOption
 	glTexParameter(type, GL_TEXTURE_MIN_FILTER, to.minFilter);
 	glTexParameter(type, GL_TEXTURE_MAG_FILTER, to.magFilter);
 	glTexImage2D(type, 0, GL_RGBA, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, data);
+	auto isMimmap = [](auto value) {
+		auto eqAnyOf = [value](auto... test) {
+			return (... || (value == test));
+		};
+		return eqAnyOf(
+				GL_NEAREST_MIPMAP_NEAREST, GL_LINEAR_MIPMAP_NEAREST, GL_NEAREST_MIPMAP_LINEAR, GL_LINEAR_MIPMAP_LINEAR);
+	};
+	if (isMimmap(to.minFilter) || isMimmap(to.magFilter)) {
+		glGenerateMipmap(type);
+	}
 }
 
 void
diff --git a/lib/collections.h b/lib/collections.h
index ea5d5dc..27eff0a 100644
--- a/lib/collections.h
+++ b/lib/collections.h
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <array>
 #include <cstdint>
+#include <ranges>
 #include <span>
 #include <tuple>
 #include <utility>
@@ -106,6 +107,15 @@ operator+=(std::vector<T...> & in, std::vector<T...> && src)
 	return in;
 }
 
+template<typename... T>
+constexpr auto
+operator+(std::vector<T...> in1, std::vector<T...> in2)
+{
+	in1.reserve(in1.size() + in2.size());
+	std::move(in2.begin(), in2.end(), std::back_inserter(in1));
+	return in1;
+}
+
 template<typename... T, typename Vn>
 [[nodiscard]] constexpr auto
 operator+(const std::vector<T...> & in, Vn && vn)
@@ -187,6 +197,14 @@ template<typename iter> struct stripiter {
 		return *this;
 	}
 
+	constexpr stripiter
+	operator++(int)
+	{
+		auto out {*this};
+		++*this;
+		return out;
+	}
+
 	constexpr stripiter &
 	operator--()
 	{
@@ -195,6 +213,14 @@ template<typename iter> struct stripiter {
 		return *this;
 	}
 
+	constexpr stripiter
+	operator--(int)
+	{
+		auto out {*this};
+		--*this;
+		return out;
+	}
+
 	constexpr auto
 	operator-(const stripiter & other) const
 	{
@@ -224,3 +250,54 @@ strip_end(IterableCollection auto & cont)
 {
 	return stripiter {cont.end()};
 }
+
+inline constexpr auto dereference = std::views::transform([](const auto & iter) -> decltype(auto) {
+	return *iter;
+});
+
+struct TriangleTriples : public std::ranges::range_adaptor_closure<TriangleTriples> {
+	decltype(auto)
+	operator()(const auto & triangleStrip) const
+	{
+		return std::views::iota(strip_begin(triangleStrip), strip_end(triangleStrip)) | dereference;
+	}
+};
+
+inline constexpr TriangleTriples triangleTriples;
+
+template<typename T, typename Dist, typename Merger>
+void
+mergeClose(std::vector<T> & range, const Dist & dist, const Merger & merger,
+		decltype(dist(range.front(), range.front())) tolerance)
+{
+	using DistanceType = decltype(tolerance);
+	std::vector<DistanceType> distances;
+	distances.reserve(range.size() - 1);
+	std::ranges::transform(range | std::views::pairwise, std::back_inserter(distances), [&dist](const auto & pair) {
+		return (std::apply(dist, pair));
+	});
+	while (distances.size() > 1) {
+		const auto closestPair = std::ranges::min_element(distances);
+		if (*closestPair > tolerance) {
+			return;
+		}
+		const auto offset = std::distance(distances.begin(), closestPair);
+		const auto idx = static_cast<std::size_t>(offset);
+		if (closestPair == distances.begin()) {
+			// Remove second element
+			range.erase(range.begin() + 1);
+			distances.erase(distances.begin());
+		}
+		else if (closestPair == --distances.end()) {
+			// Remove second from last element
+			range.erase(range.end() - 2);
+			distances.erase(distances.end() - 1);
+		}
+		else {
+			range[idx] = merger(range[idx], range[idx + 1]);
+			range.erase(range.begin() + offset + 1);
+			distances.erase(distances.begin() + offset);
+		}
+		distances[idx] = dist(range[idx], range[idx + 1]);
+	}
+}
diff --git a/lib/maths.h b/lib/maths.h
index 3959896..2049c78 100644
--- a/lib/maths.h
+++ b/lib/maths.h
@@ -1,6 +1,8 @@
 #pragma once
 
 #include "config/types.h"
+#include <algorithm>
+#include <array>
 #include <cmath>
 #include <glm/glm.hpp>
 #include <glm/gtc/constants.hpp>
@@ -12,8 +14,11 @@
 template<typename T>
 concept Arithmetic = std::is_arithmetic_v<T>;
 
+template<Arithmetic T> using CalcType = std::conditional_t<std::is_floating_point_v<T>, T, int64_t>;
+template<Arithmetic T> using DifferenceType = std::conditional_t<std::is_floating_point_v<T>, T, float>;
+
 struct Arc : public std::pair<Angle, Angle> {
-	template<glm::length_t Lc, glm::length_t Le, Arithmetic T, glm::qualifier Q>
+	template<glm::length_t Lc, glm::length_t Le, Arithmetic T, glm::qualifier Q = glm::defaultp>
 		requires(Lc >= 2, Le >= 2)
 	Arc(const glm::vec<Lc, T, Q> & centre, const glm::vec<Le, T, Q> & e0p, const glm::vec<Le, T, Q> & e1p) :
 		Arc {RelativePosition2D {e0p.xy() - centre.xy()}, RelativePosition2D {e1p.xy() - centre.xy()}}
@@ -36,6 +41,26 @@ struct Arc : public std::pair<Angle, Angle> {
 	}
 };
 
+template<typename T, glm::qualifier Q = glm::defaultp> struct ArcSegment : public Arc {
+	using PointType = glm::vec<2, T, Q>;
+
+	constexpr ArcSegment(PointType centre, PointType ep0, PointType ep1);
+
+	PointType centre;
+	PointType ep0;
+	PointType ep1;
+	RelativeDistance radius;
+
+	[[nodiscard]] constexpr std::optional<std::pair<glm::vec<2, T, Q>, Angle>> crossesLineAt(
+			const glm::vec<2, T, Q> & lineStart, const glm::vec<2, T, Q> & lineEnd) const;
+
+	[[nodiscard]] constexpr bool
+	angleWithinArc(Angle angle) const
+	{
+		return first <= angle && angle <= second;
+	}
+};
+
 constexpr const RelativePosition3D up {0, 0, 1}; // NOLINT(readability-identifier-length)
 constexpr const RelativePosition3D down {0, 0, -1};
 constexpr const RelativePosition3D north {0, 1, 0};
@@ -79,13 +104,33 @@ operator-(const GlobalPosition<D> & global, const CalcPosition<D> & relative)
 	return global - GlobalPosition<D>(relative);
 }
 
-template<glm::length_t D, std::integral T, glm::qualifier Q>
-constexpr RelativePosition<D>
+template<glm::length_t D, Arithmetic T, glm::qualifier Q = glm::defaultp>
+using DifferenceVector = glm::vec<D, DifferenceType<T>, Q>;
+
+template<glm::length_t D, Arithmetic T, glm::qualifier Q = glm::defaultp>
+constexpr DifferenceVector<D, T, Q>
 difference(const glm::vec<D, T, Q> & globalA, const glm::vec<D, T, Q> & globalB)
 {
 	return globalA - globalB;
 }
 
+template<glm::length_t D, Arithmetic T, glm::qualifier Q = glm::defaultp>
+using CalcVector = glm::vec<D, CalcType<T>, Q>;
+
+template<glm::length_t D, Arithmetic T, glm::qualifier Q = glm::defaultp>
+constexpr CalcVector<D, T, Q>
+calcDifference(const glm::vec<D, T, Q> & globalA, const glm::vec<D, T, Q> & globalB)
+{
+	return globalA - globalB;
+}
+
+template<glm::length_t D, Arithmetic T, glm::qualifier Q = glm::defaultp>
+constexpr auto
+distance(const glm::vec<D, T, Q> & pointA, const glm::vec<D, T, Q> & pointB)
+{
+	return glm::length(difference(pointA, pointB));
+}
+
 glm::mat4 flat_orientation(const Rotation3D & diff);
 
 namespace {
@@ -122,7 +167,8 @@ namespace {
 	}
 
 	// Helper to lookup into a matrix given an xy vector coordinate
-	template<glm::length_t C, glm::length_t R, Arithmetic T, glm::qualifier Q, std::integral I = glm::length_t>
+	template<glm::length_t C, glm::length_t R, Arithmetic T, glm::qualifier Q = glm::defaultp,
+			std::integral I = glm::length_t>
 	constexpr auto &
 	operator^(glm::mat<C, R, T, Q> & matrix, const glm::vec<2, I> rowCol)
 	{
@@ -130,7 +176,7 @@ namespace {
 	}
 
 	// Create a matrix for the angle, given the targets into the matrix
-	template<glm::length_t D, std::floating_point T, glm::qualifier Q, std::integral I = glm::length_t>
+	template<glm::length_t D, std::floating_point T, glm::qualifier Q = glm::defaultp, std::integral I = glm::length_t>
 	constexpr auto
 	rotation(const T angle, const glm::vec<2, I> cos1, const glm::vec<2, I> sin1, const glm::vec<2, I> cos2,
 			const glm::vec<2, I> negSin1)
@@ -220,7 +266,7 @@ vector_pitch(const glm::vec<D, T, Q> & diff)
 	return std::atan(diff.z);
 }
 
-template<std::floating_point T, glm::qualifier Q>
+template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<2, T, Q>
 vector_normal(const glm::vec<2, T, Q> & vector)
 {
@@ -242,7 +288,7 @@ sq(T value)
 	return value * value;
 }
 
-template<glm::qualifier Q>
+template<glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<3, int64_t, Q>
 crossProduct(const glm::vec<3, int64_t, Q> & valueA, const glm::vec<3, int64_t, Q> & valueB)
 {
@@ -253,14 +299,14 @@ crossProduct(const glm::vec<3, int64_t, Q> & valueA, const glm::vec<3, int64_t,
 	};
 }
 
-template<std::integral T, glm::qualifier Q>
+template<std::integral T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<3, T, Q>
 crossProduct(const glm::vec<3, T, Q> & valueA, const glm::vec<3, T, Q> & valueB)
 {
 	return crossProduct<Q>(valueA, valueB);
 }
 
-template<std::floating_point T, glm::qualifier Q>
+template<std::floating_point T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<3, T, Q>
 crossProduct(const glm::vec<3, T, Q> & valueA, const glm::vec<3, T, Q> & valueB)
 {
@@ -275,35 +321,35 @@ ratio(const Ta valueA, const Tb valueB)
 	return static_cast<R>((static_cast<Common>(valueA) / static_cast<Common>(valueB)));
 }
 
-template<Arithmetic R = float, Arithmetic T, glm::qualifier Q>
+template<Arithmetic R = float, Arithmetic T, glm::qualifier Q = glm::defaultp>
 constexpr auto
 ratio(const glm::vec<2, T, Q> & value)
 {
 	return ratio<R>(value.x, value.y);
 }
 
-template<glm::length_t L = 3, std::floating_point T, glm::qualifier Q>
+template<glm::length_t L = 3, std::floating_point T, glm::qualifier Q = glm::defaultp>
 constexpr auto
 perspective_divide(const glm::vec<4, T, Q> & value)
 {
 	return value / value.w;
 }
 
-template<glm::length_t L1, glm::length_t L2, Arithmetic T, glm::qualifier Q>
+template<glm::length_t L1, glm::length_t L2, Arithmetic T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<L1 + L2, T, Q>
 operator||(const glm::vec<L1, T, Q> valueA, const glm::vec<L2, T, Q> valueB)
 {
 	return {valueA, valueB};
 }
 
-template<glm::length_t L, Arithmetic T, glm::qualifier Q>
+template<glm::length_t L, Arithmetic T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<L + 1, T, Q>
 operator||(const glm::vec<L, T, Q> valueA, const T valueB)
 {
 	return {valueA, valueB};
 }
 
-template<glm::length_t L, std::floating_point T, glm::qualifier Q>
+template<glm::length_t L, std::floating_point T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<L, T, Q>
 perspectiveMultiply(const glm::vec<L, T, Q> & base, const glm::mat<L + 1, L + 1, T, Q> & mutation)
 {
@@ -311,7 +357,7 @@ perspectiveMultiply(const glm::vec<L, T, Q> & base, const glm::mat<L + 1, L + 1,
 	return mutated / mutated.w;
 }
 
-template<glm::length_t L, std::floating_point T, glm::qualifier Q>
+template<glm::length_t L, std::floating_point T, glm::qualifier Q = glm::defaultp>
 constexpr glm::vec<L, T, Q>
 perspectiveApply(glm::vec<L, T, Q> & base, const glm::mat<L + 1, L + 1, T, Q> & mutation)
 {
@@ -331,8 +377,6 @@ normalize(T ang)
 	return ang;
 }
 
-template<Arithmetic T> using CalcType = std::conditional_t<std::is_floating_point_v<T>, T, int64_t>;
-
 template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 [[nodiscard]] constexpr std::optional<glm::vec<2, T, Q>>
 linesIntersectAt(const glm::vec<2, T, Q> Aabs, const glm::vec<2, T, Q> Babs, const glm::vec<2, T, Q> Cabs,
@@ -359,7 +403,7 @@ linesIntersectAt(const glm::vec<2, T, Q> Aabs, const glm::vec<2, T, Q> Babs, con
 	return Aabs + CVec {(b1 * c2) / -determinant, (a1 * c2) / determinant};
 }
 
-template<Arithmetic T, glm::qualifier Q>
+template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 std::pair<glm::vec<2, T, Q>, bool>
 find_arc_centre(glm::vec<2, T, Q> start, Rotation2D startDir, glm::vec<2, T, Q> end, Rotation2D endDir)
 {
@@ -372,7 +416,7 @@ find_arc_centre(glm::vec<2, T, Q> start, Rotation2D startDir, glm::vec<2, T, Q>
 	throw std::runtime_error("no intersection");
 }
 
-template<Arithmetic T, glm::qualifier Q>
+template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 std::pair<glm::vec<2, T, Q>, bool>
 find_arc_centre(glm::vec<2, T, Q> start, Angle entrys, glm::vec<2, T, Q> end, Angle entrye)
 {
@@ -382,7 +426,7 @@ find_arc_centre(glm::vec<2, T, Q> start, Angle entrys, glm::vec<2, T, Q> end, An
 	return find_arc_centre(start, sincos(entrys + half_pi), end, sincos(entrye - half_pi));
 }
 
-template<Arithmetic T, glm::qualifier Q>
+template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 Angle
 find_arcs_radius(glm::vec<2, T, Q> start, Rotation2D ad, glm::vec<2, T, Q> end, Rotation2D bd)
 {
@@ -407,7 +451,7 @@ find_arcs_radius(glm::vec<2, T, Q> start, Rotation2D ad, glm::vec<2, T, Q> end,
 			/ (2 * (sq(X) - 2 * X * Z + sq(Z) + sq(Y) - 2 * Y * W + sq(W) - 4));
 }
 
-template<Arithmetic T, glm::qualifier Q>
+template<Arithmetic T, glm::qualifier Q = glm::defaultp>
 std::pair<Angle, Angle>
 find_arcs_radius(glm::vec<2, T, Q> start, Angle entrys, glm::vec<2, T, Q> end, Angle entrye)
 {
@@ -424,6 +468,13 @@ midpoint(const std::pair<T, T> & v)
 	return std::midpoint(v.first, v.second);
 }
 
+template<glm::length_t D, std::integral T, glm::qualifier Q = glm::defaultp>
+auto
+midpoint(const glm::vec<D, T, Q> & valueA, const glm::vec<D, T, Q> & valueB)
+{
+	return valueA + (valueB - valueA) / 2;
+}
+
 // std::pow is not constexpr
 template<Arithmetic T>
 	requires requires(T n) { n *= n; }
@@ -461,3 +512,78 @@ operator"" _degrees(long double degrees)
 {
 	return static_cast<float>(degrees) * degreesToRads;
 }
+
+// Late implementations due to dependencies
+template<typename T, glm::qualifier Q>
+constexpr ArcSegment<T, Q>::ArcSegment(PointType centre, PointType ep0, PointType ep1) :
+	Arc {centre, ep0, ep1}, centre {centre}, ep0 {ep0}, ep1 {ep1}, radius {::distance(centre, ep0)}
+{
+}
+
+template<typename T, glm::qualifier Q>
+[[nodiscard]] constexpr std::optional<std::pair<glm::vec<2, T, Q>, Angle>>
+ArcSegment<T, Q>::crossesLineAt(const glm::vec<2, T, Q> & lineStart, const glm::vec<2, T, Q> & lineEnd) const
+{
+	// Based on formulas from https://mathworld.wolfram.com/Circle-LineIntersection.html
+	const auto lineDiff = difference(lineEnd, lineStart);
+	const auto lineLen = glm::length(lineDiff);
+	const auto lineRelStart = difference(lineStart, centre);
+	const auto lineRelEnd = difference(lineEnd, centre);
+	const auto determinant = (lineRelStart.x * lineRelEnd.y) - (lineRelEnd.x * lineRelStart.y);
+	const auto discriminant = (radius * radius * lineLen * lineLen) - (determinant * determinant);
+	if (discriminant < 0) {
+		return std::nullopt;
+	}
+
+	const auto rootDiscriminant = std::sqrt(discriminant);
+	const auto drdr = lineLen * lineLen;
+	const RelativeDistance sgn = (lineDiff.y < 0 ? -1 : 1);
+	std::array<std::pair<RelativePosition2D, Angle>, 2> points;
+	std::ranges::transform(std::initializer_list {1, -1}, points.begin(), [&](RelativeDistance N) {
+		const auto point = RelativePosition2D {((determinant * lineDiff.y) + sgn * lineDiff.x * rootDiscriminant * N),
+								   ((-determinant * lineDiff.x) + std::abs(lineDiff.y) * rootDiscriminant * N)}
+				/ drdr;
+		return std::make_pair(point, vector_yaw(point));
+	});
+	const auto end
+			= std::remove_if(points.begin(), points.end(), [this, lineRelStart, lineDiff, drdr](const auto point) {
+				  const auto dot = glm::dot(lineDiff, point.first - lineRelStart);
+				  return !angleWithinArc(point.second) || dot < 0 || dot > drdr;
+			  });
+	if (points.begin() == end) {
+		return std::nullopt;
+	}
+	const auto first = *std::ranges::min_element(points.begin(), end, {}, [lineRelStart](const auto point) {
+		return glm::distance(lineRelStart, point.first);
+	});
+	return std::make_pair(centre + first.first, first.second);
+}
+
+namespace {
+	template<template<typename> typename Op>
+	[[nodiscard]] constexpr auto
+	pointLineOp(const GlobalPosition2D point, const GlobalPosition2D end1, const GlobalPosition2D end2)
+	{
+		return Op {}(CalcDistance(end2.x - end1.x) * CalcDistance(point.y - end1.y),
+				CalcDistance(end2.y - end1.y) * CalcDistance(point.x - end1.x));
+	}
+}
+
+constexpr auto pointLeftOfLine = pointLineOp<std::greater>;
+constexpr auto pointLeftOfOrOnLine = pointLineOp<std::greater_equal>;
+
+[[nodiscard]] constexpr bool
+linesCross(const GlobalPosition2D lineAend1, const GlobalPosition2D lineAend2, const GlobalPosition2D lineBend1,
+		const GlobalPosition2D lineBend2)
+{
+	return (pointLeftOfLine(lineAend2, lineBend1, lineBend2) == pointLeftOfLine(lineAend1, lineBend2, lineBend1))
+			&& (pointLeftOfLine(lineBend1, lineAend1, lineAend2) == pointLeftOfLine(lineBend2, lineAend2, lineAend1));
+}
+
+[[nodiscard]] constexpr bool
+linesCrossLtR(const GlobalPosition2D lineAend1, const GlobalPosition2D lineAend2, const GlobalPosition2D lineBend1,
+		const GlobalPosition2D lineBend2)
+{
+	return pointLeftOfLine(lineAend2, lineBend1, lineBend2) && pointLeftOfLine(lineAend1, lineBend2, lineBend1)
+			&& pointLeftOfLine(lineBend1, lineAend1, lineAend2) && pointLeftOfLine(lineBend2, lineAend2, lineAend1);
+}
diff --git a/lib/ray.h b/lib/ray.h
index a831270..793e21e 100644
--- a/lib/ray.h
+++ b/lib/ray.h
@@ -27,8 +27,7 @@ public:
 		const auto n2 = crossProduct(direction, n);
 		const auto c1 = p1 + PositionType((glm::dot(RelativePosition3D(start - p1), n2) / glm::dot(d1, n2)) * d1);
 		const auto difflength = glm::length(diff);
-		if (glm::length(RelativePosition3D(c1 - p1)) > difflength
-				|| glm::length(RelativePosition3D(c1 - e1)) > difflength) {
+		if (::distance(c1, p1) > difflength || ::distance(c1, e1) > difflength) {
 			return std::numeric_limits<typename PositionType::value_type>::infinity();
 		}
 		return static_cast<PositionType::value_type>(glm::abs(glm::dot(n, RelativePosition3D(p1 - start))));
diff --git a/lib/sorting.h b/lib/sorting.h
index 777de00..be5a7e2 100644
--- a/lib/sorting.h
+++ b/lib/sorting.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <functional>
 #include <glm/fwd.hpp>
 #include <type_traits>
 
@@ -30,6 +31,14 @@ template<typename T, auto M> struct PtrMemberSorter : public PtrSorter<T> {
 	}
 };
 
+template<auto Proj> struct SortedBy {
+	auto
+	operator()(const auto & left, const auto & right) const
+	{
+		return (std::invoke(Proj, left) < std::invoke(Proj, right));
+	}
+};
+
 struct CompareBy {
 	glm::length_t index;
 
diff --git a/lib/stream_support.h b/lib/stream_support.h
index f21622a..5f276fd 100644
--- a/lib/stream_support.h
+++ b/lib/stream_support.h
@@ -4,6 +4,7 @@
 #include <glm/glm.hpp>
 #include <iostream>
 #include <maths.h>
+#include <optional>
 #include <source_location>
 #include <span>
 #include <sstream>
@@ -82,6 +83,16 @@ namespace std {
 	{
 		return s << EnumTypeDetails<E>::typeName << "::" << EnumDetails<E>::to_string(e).value();
 	}
+
+	template<typename T>
+	inline std::ostream &
+	operator<<(std::ostream & s, const std::optional<T> & v)
+	{
+		if (v) {
+			return s << *v;
+		}
+		return s << "nullopt";
+	}
 }
 
 template<typename T>
@@ -104,3 +115,4 @@ namespace {
 }
 
 #define CLOG(x) clogImpl(x, #x)
+#define CLOGf(...) clogImpl(std::format(__VA_ARGS__), "msg")
diff --git a/lib/triangle.h b/lib/triangle.h
index 812bfab..abd697c 100644
--- a/lib/triangle.h
+++ b/lib/triangle.h
@@ -31,6 +31,13 @@ struct Triangle : public glm::vec<3, glm::vec<Dim, T, Q>> {
 		return glm::length(crossProduct(sideDifference(1), sideDifference(2))) / T {2};
 	}
 
+	[[nodiscard]] constexpr auto
+	area() const
+		requires(Dim == 2)
+	{
+		return std::abs((sideDifference(1).x * sideDifference(2).y) - (sideDifference(2).x * sideDifference(1).y)) / 2;
+	}
+
 	[[nodiscard]] constexpr Normal3D
 	normal() const
 		requires(Dim == 3)
@@ -38,6 +45,12 @@ struct Triangle : public glm::vec<3, glm::vec<Dim, T, Q>> {
 		return crossProduct(sideDifference(1), sideDifference(2));
 	}
 
+	[[nodiscard]] constexpr auto
+	height()
+	{
+		return (area() * 2) / ::distance(p(0), p(1));
+	}
+
 	[[nodiscard]] constexpr Normal3D
 	nnormal() const
 		requires(Dim == 3)
@@ -52,6 +65,12 @@ struct Triangle : public glm::vec<3, glm::vec<Dim, T, Q>> {
 	}
 
 	[[nodiscard]] constexpr auto
+	calcSideDifference(glm::length_t side) const
+	{
+		return calcDifference(p(side), p(0));
+	}
+
+	[[nodiscard]] constexpr auto
 	angle(glm::length_t corner) const
 	{
 		return Arc {P(corner), P(corner + 2), P(corner + 1)}.length();
@@ -71,6 +90,14 @@ struct Triangle : public glm::vec<3, glm::vec<Dim, T, Q>> {
 	}
 
 	[[nodiscard]] constexpr auto
+	isUp() const
+	{
+		const auto edgeAB = sideDifference(1);
+		const auto edgeAC = sideDifference(2);
+		return edgeAB.x * edgeAC.y >= edgeAB.y * edgeAC.x;
+	}
+
+	[[nodiscard]] constexpr auto
 	p(const glm::length_t idx) const
 	{
 		return Base::operator[](idx);
@@ -105,4 +132,25 @@ struct Triangle : public glm::vec<3, glm::vec<Dim, T, Q>> {
 	{
 		return begin() + 3;
 	}
+
+	[[nodiscard]]
+	constexpr auto
+	positionOnPlane(const glm::vec<2, T, Q> coord2d) const
+		requires(Dim == 3)
+	{
+		const auto edgeCrossProduct = crossProduct(calcSideDifference(1), calcSideDifference(2));
+		return coord2d
+				|| static_cast<T>(
+						((edgeCrossProduct.x * p(0).x) + (edgeCrossProduct.y * p(0).y) + (edgeCrossProduct.z * p(0).z)
+								- (edgeCrossProduct.x * coord2d.x) - (edgeCrossProduct.y * coord2d.y))
+						/ edgeCrossProduct.z);
+	}
+
+	[[nodiscard]]
+	constexpr bool
+	containsPoint(const GlobalPosition2D coord) const
+	{
+		return pointLeftOfOrOnLine(coord, p(0), p(1)) && pointLeftOfOrOnLine(coord, p(1), p(2))
+				&& pointLeftOfOrOnLine(coord, p(2), p(0));
+	}
 };
diff --git a/lib/util.h b/lib/util.h
index 290492f..cd7971b 100644
--- a/lib/util.h
+++ b/lib/util.h
@@ -3,6 +3,7 @@
 #include <algorithm> // IWYU pragma: keep
 #include <array>
 #include <cstddef>
+#include <tuple>
 
 template<typename T, std::size_t N>
 constexpr auto
@@ -12,3 +13,23 @@ transform_array(const std::array<T, N> & in, auto && transform)
 	std::transform(in.begin(), in.end(), out.begin(), transform);
 	return out;
 }
+
+namespace {
+	template<size_t... N> struct GetNth {
+		decltype(auto)
+		operator()(const auto & tup) const
+		{
+			if constexpr (sizeof...(N) == 1) {
+				return std::get<N...>(tup);
+			}
+			else {
+				return std::tie(std::get<N>(tup)...);
+			}
+		}
+	};
+}
+
+template<size_t... N> inline constexpr auto Nth = GetNth<N...> {};
+inline constexpr auto GetFirst = Nth<0>;
+inline constexpr auto GetSecond = Nth<1>;
+inline constexpr auto GetSwapped = Nth<0, 1>;
diff --git a/test/Jamfile.jam b/test/Jamfile.jam
index 3ab4c4c..c5c49d2 100644
--- a/test/Jamfile.jam
+++ b/test/Jamfile.jam
@@ -47,7 +47,7 @@ lib test : [ glob *.cpp : test-*.cpp perf-*.cpp ] ;
 run test-collection.cpp ;
 run test-maths.cpp ;
 run test-lib.cpp ;
-run test-geoData.cpp : -- : [ sequence.insertion-sort [ glob-tree $(fixtures)/geoData : *.json ] fixtures/height/SD19.asc ] : <library>test ;
+run [ glob test-geoData*.cpp ] : -- : [ sequence.insertion-sort [ glob-tree $(fixtures)/geoData : *.json ] fixtures/height/SD19.asc ] : <library>test ;
 run test-network.cpp : : : <library>test ;
 run test-persistence.cpp : -- : [ sequence.insertion-sort [ glob-tree $(fixtures)/json : *.json ] ] : <library>test ;
 run test-text.cpp : -- : test-glContainer : <library>test ;
diff --git a/test/fixtures/geoData/deform/multi1.json b/test/fixtures/geoData/deform/multi1.json
new file mode 100644
index 0000000..c7456b6
--- /dev/null
+++ b/test/fixtures/geoData/deform/multi1.json
@@ -0,0 +1,21 @@
+[
+	[
+		[
+			[
+				100,
+				100,
+				100
+			],
+			[
+				150,
+				100,
+				100
+			],
+			[
+				100,
+				150,
+				100
+			]
+		]
+	]
+]
diff --git a/test/perf-geoData.cpp b/test/perf-geoData.cpp
index 4d4505e..d9ea8c6 100644
--- a/test/perf-geoData.cpp
+++ b/test/perf-geoData.cpp
@@ -33,10 +33,25 @@ namespace {
 			});
 		}
 	}
+
+	void
+	terrain_deform(benchmark::State & state)
+	{
+		std::array<GlobalPosition3D, 3> points {{
+				{315555000, 495556000, 0},
+				{315655000, 495556000, 0},
+				{315655000, 495557000, 0},
+		}};
+		for (auto _ : state) {
+			auto geoData {tm};
+			benchmark::DoNotOptimize(geoData.setHeights(points, GeoData::SetHeightsOpts {.surface = nullptr}));
+		}
+	}
 }
 
 BENCHMARK(terrain_findPoint);
 BENCHMARK(terrain_walk);
 BENCHMARK(terrain_walkBoundary);
+BENCHMARK(terrain_deform);
 
 BENCHMARK_MAIN();
diff --git a/test/test-geoData-counts.cpp b/test/test-geoData-counts.cpp
new file mode 100644
index 0000000..446a68a
--- /dev/null
+++ b/test/test-geoData-counts.cpp
@@ -0,0 +1,65 @@
+#include <OpenMesh/Core/IO/MeshIO.hh>
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/test_tools.hpp>
+#include <boost/test/unit_test_suite.hpp>
+
+#include "game/geoData.h"
+
+using GeoMutation = std::function<void(GeoData &)>;
+using Something = std::tuple<const char *, GeoMutation, size_t, size_t, size_t>;
+BOOST_TEST_DONT_PRINT_LOG_VALUE(GeoMutation);
+
+BOOST_TEST_DECORATOR(*boost::unit_test::timeout(2));
+
+BOOST_DATA_TEST_CASE(deformLogical,
+		boost::unit_test::data::make<Something>({
+				{"nochange", [](GeoData &) {}, 16, 33, 18}, // No change base case
+				{"simple",
+						[](GeoData & geoData) {
+							Surface surface;
+							// Basic triangle, no crossing, simple case
+							geoData.setHeights(std::array<GlobalPosition3D, 3> {{
+													   {2000, 8000, 1000},
+													   {2000, 4000, 1000},
+													   {6000, 8000, 1000},
+											   }},
+									{.surface = &surface, .nearNodeTolerance = 0});
+						},
+						19, 42, 24},
+				{"simple-cross",
+						[](GeoData & geoData) {
+							Surface surface;
+							// Basic triangle, with crossing, reasonably simple case
+							geoData.setHeights(std::array<GlobalPosition3D, 3> {{
+													   {2000, 8000, 1000},
+													   {3000, 4000, 1000},
+													   {4000, 9000, 1000},
+											   }},
+									{.surface = &surface, .nearNodeTolerance = 0});
+						},
+						19, 42, 24},
+				{"quad-multi-cross",
+						[](GeoData & geoData) {
+							Surface surface;
+							// Basic quad, with crossing, spans into adjacent original triangle, should remove that edge
+							geoData.setHeights(std::array<GlobalPosition3D, 4> {{
+													   {2000, 8000, 1000},
+													   {3000, 4000, 1000},
+													   {4000, 9000, 1000},
+													   {8000, 2000, 1000},
+											   }},
+									{.surface = &surface, .nearNodeTolerance = 0});
+						},
+						20, 45, 26},
+		}),
+		name, func, expVertices, expEdges, expFaces)
+{
+	auto geoData = GeoData::createFlat({0, 0}, {30'000, 30'000}, 1000);
+
+	BOOST_REQUIRE_NO_THROW(func(geoData));
+	OpenMesh::IO::write_mesh(geoData, std::format("/tmp/mesh-{}.obj", name));
+
+	BOOST_CHECK_EQUAL(geoData.n_vertices(), expVertices);
+	BOOST_CHECK_EQUAL(geoData.n_edges(), expEdges);
+	BOOST_CHECK_EQUAL(geoData.n_faces(), expFaces);
+}
diff --git a/test/test-geoData.cpp b/test/test-geoData.cpp
index bd1ff87..2332513 100644
--- a/test/test-geoData.cpp
+++ b/test/test-geoData.cpp
@@ -29,26 +29,9 @@ BOOST_AUTO_TEST_CASE(loadSuccess)
 	BOOST_CHECK_EQUAL(upper, GlobalPosition3D(319950000, 499950000, 571600));
 }
 
-BOOST_AUTO_TEST_CASE(normalsAllPointUp)
+BOOST_AUTO_TEST_CASE(sanityCheck)
 {
-	BOOST_CHECK(std::ranges::all_of(vertices(), [this](auto && vertex) {
-		return normal(vertex).z > 0;
-	}));
-}
-
-BOOST_AUTO_TEST_CASE(trianglesContainsPoints)
-{
-	const auto face = face_handle(0);
-
-	BOOST_TEST_CONTEXT(this->triangle<2>(face)) {
-		BOOST_CHECK(triangleContainsPoint(GlobalPosition2D {xllcorner, yllcorner}, face));
-		BOOST_CHECK(triangleContainsPoint(GlobalPosition2D {xllcorner + cellsize, yllcorner + cellsize}, face));
-		BOOST_CHECK(triangleContainsPoint(GlobalPosition2D {xllcorner, yllcorner + cellsize}, face));
-		BOOST_CHECK(triangleContainsPoint(GlobalPosition2D {xllcorner + 1, yllcorner + 1}, face));
-		BOOST_CHECK(triangleContainsPoint(GlobalPosition2D {xllcorner + 1, yllcorner + 2}, face));
-		BOOST_CHECK(!triangleContainsPoint(GlobalPosition2D {xllcorner + 3, yllcorner + 2}, face));
-		BOOST_CHECK(!triangleContainsPoint(GlobalPosition2D {xllcorner + cellsize, yllcorner}, face));
-	}
+	BOOST_CHECK_NO_THROW(sanityCheck());
 }
 
 BOOST_AUTO_TEST_SUITE_END();
@@ -105,7 +88,7 @@ BOOST_DATA_TEST_CASE(findPositionAt,
 		}),
 		p, h)
 {
-	BOOST_CHECK_EQUAL(fixedTerrtain.positionAt(p), GlobalPosition3D(p, h));
+	BOOST_CHECK_EQUAL(fixedTerrtain.positionAt(p), p || h);
 }
 
 using FindRayIntersectData = std::tuple<GlobalPosition3D, Direction3D, GlobalPosition3D>;
@@ -150,8 +133,11 @@ BOOST_DATA_TEST_CASE(walkTerrain,
 		from, to, visits)
 {
 	std::vector<int> visited;
-	BOOST_CHECK_NO_THROW(fixedTerrtain.walk(from, to, [&visited](auto fh) {
-		visited.emplace_back(fh.idx());
+	BOOST_CHECK_NO_THROW(fixedTerrtain.walk(from, to, [&visited](auto step) {
+		if (!visited.empty()) {
+			BOOST_CHECK_EQUAL(step.previous.idx(), visited.back());
+		}
+		visited.emplace_back(step.current.idx());
 	}));
 	BOOST_CHECK_EQUAL_COLLECTIONS(visited.begin(), visited.end(), visits.begin(), visits.end());
 }
@@ -183,28 +169,57 @@ BOOST_DATA_TEST_CASE(walkTerrainUntil,
 		from, to, visits)
 {
 	std::vector<int> visited;
-	BOOST_CHECK_NO_THROW(fixedTerrtain.walkUntil(from, to, [&visited](auto fh) {
-		visited.emplace_back(fh.idx());
+	BOOST_CHECK_NO_THROW(fixedTerrtain.walkUntil(from, to, [&visited](const auto & step) {
+		visited.emplace_back(step.current.idx());
 		return visited.size() >= 5;
 	}));
 	BOOST_CHECK_EQUAL_COLLECTIONS(visited.begin(), visited.end(), visits.begin(), visits.end());
 }
 
-BOOST_AUTO_TEST_CASE(triangle_helpers)
+using WalkTerrainCurveData = std::tuple<GlobalPosition2D, GlobalPosition2D, GlobalPosition2D, std::vector<int>,
+		std::vector<GlobalPosition2D>>;
+
+BOOST_TEST_DECORATOR(*boost::unit_test::timeout(1))
+
+BOOST_DATA_TEST_CASE(walkTerrainCurveSetsFromFace,
+		boost::unit_test::data::make<WalkTerrainCurveData>({
+				{{310002000, 490003000}, {310002000, 490003000}, {310002000, 490003000}, {0}, {}},
+				{{310003000, 490002000}, {310003000, 490002000}, {310003000, 490002000}, {1}, {}},
+				{{310202000, 490203000}, {310002000, 490003000}, {310002000, 490203000},
+						{1600, 1601, 1202, 1201, 802, 803, 404, 403, 4, 3, 2, 1, 0},
+						{
+								{310201997, 490201997},
+								{310201977, 490200000},
+								{310200000, 490174787},
+								{310194850, 490150000},
+								{310192690, 490142690},
+								{310173438, 490100000},
+								{310150000, 490068479},
+								{310130806, 490050000},
+								{310100000, 490028656},
+								{310062310, 490012310},
+								{310050000, 490008845},
+								{310003003, 490003003},
+						}},
+				{{310999999, 490205000}, {310999999, 490203000}, {310999000, 490204000}, {1631, 1632, 1631},
+						{
+								{311000000, 490204999},
+								{311000000, 490203001},
+						}},
+		}),
+		from, to, centre, visits, exits)
 {
-	constexpr static GeoData::Triangle<3> t {{0, 0, 0}, {5, 0, 0}, {5, 5, 0}};
-
-	BOOST_CHECK_EQUAL(t.nnormal(), up);
-	BOOST_CHECK_CLOSE(t.angle(0), quarter_pi, 0.01F);
-	BOOST_CHECK_CLOSE(t.angleAt({0, 0, 0}), quarter_pi, 0.01F);
-	BOOST_CHECK_CLOSE(t.angle(1), half_pi, 0.01F);
-	BOOST_CHECK_CLOSE(t.angleAt({5, 0, 0}), half_pi, 0.01F);
-	BOOST_CHECK_CLOSE(t.angle(2), quarter_pi, 0.01F);
-	BOOST_CHECK_CLOSE(t.angleAt({5, 5, 0}), quarter_pi, 0.01F);
-
-	BOOST_CHECK_CLOSE(t.angleAt({0, 1, 0}), 0.F, 0.01F);
+	BOOST_REQUIRE_EQUAL(visits.size(), exits.size() + 1);
 
-	BOOST_CHECK_CLOSE(t.area(), 12.5F, 0.01F);
+	std::vector<int> visited;
+	std::vector<GlobalPosition2D> exited;
+	BOOST_CHECK_NO_THROW(fixedTerrtain.walk(from, to, centre, [&](const auto & step) {
+		visited.emplace_back(step.current.idx());
+		BOOST_REQUIRE(!std::ranges::contains(exited, step.exitPosition));
+		exited.emplace_back(step.exitPosition);
+	}));
+	BOOST_CHECK_EQUAL_COLLECTIONS(visited.begin(), visited.end(), visits.begin(), visits.end());
+	BOOST_CHECK_EQUAL_COLLECTIONS(exited.begin() + 1, exited.end(), exits.begin(), exits.end());
 }
 
 using FindEntiesData = std::tuple<GlobalPosition2D, GlobalPosition2D, int>;
@@ -229,18 +244,13 @@ BOOST_DATA_TEST_CASE(deform, loadFixtureJson<DeformTerrainData>("geoData/deform/
 {
 	Surface surface;
 	surface.colorBias = RGB {0, 0, 1};
-	auto gd = std::make_shared<GeoData>(GeoData::createFlat({0, 0}, {1000000, 1000000}, 100));
-	BOOST_CHECK_NO_THROW(gd->setHeights(points, {.surface = surface}));
-	BOOST_CHECK(std::ranges::all_of(gd->vertices(), [&gd](auto && vertex) {
-		return gd->normal(vertex).z > 0;
-	}));
 
 	ApplicationBase ab;
 	TestMainWindow tmw;
 	TestRenderOutput tro {{640, 480}};
 
 	struct TestTerrain : public SceneProvider {
-		explicit TestTerrain(std::shared_ptr<GeoData> gd) : terrain(std::move(gd)) { }
+		explicit TestTerrain(GeoData gd) : terrain(std::move(gd)) { }
 
 		const Terrain terrain;
 
@@ -269,7 +279,11 @@ BOOST_DATA_TEST_CASE(deform, loadFixtureJson<DeformTerrainData>("geoData/deform/
 		}
 	};
 
-	TestTerrain t {gd};
+	TestTerrain t {[&points, &surface]() {
+		auto gd = GeoData::createFlat({0, 0}, {1000000, 1000000}, 100);
+		BOOST_CHECK_NO_THROW(gd.setHeights(points, {.surface = &surface}));
+		return gd;
+	}()};
 	SceneRenderer ss {tro.size, tro.output};
 	std::for_each(cams.begin(), cams.end(), [&ss, &t, &tro](const auto & cam) {
 		ss.camera.setView(cam.first.first, glm::normalize(cam.first.second));
@@ -277,3 +291,17 @@ BOOST_DATA_TEST_CASE(deform, loadFixtureJson<DeformTerrainData>("geoData/deform/
 		Texture::save(tro.outImage, cam.second.c_str());
 	});
 }
+
+BOOST_TEST_DECORATOR(*boost::unit_test::timeout(2));
+
+BOOST_DATA_TEST_CASE(
+		deformMulti, loadFixtureJson<std::vector<std::vector<GlobalPosition3D>>>("geoData/deform/multi1.json"), points)
+{
+	BOOST_REQUIRE(!points.empty());
+	Surface surface;
+	auto gd = std::make_shared<GeoData>(GeoData::createFlat({0, 0}, {1000000, 1000000}, 100));
+	for (const auto & strip : points) {
+		BOOST_REQUIRE_GE(strip.size(), 3);
+		BOOST_CHECK_NO_THROW(gd->setHeights(strip, {.surface = &surface, .nearNodeTolerance = 50}));
+	}
+}
diff --git a/test/test-lib.cpp b/test/test-lib.cpp
index 5f0b5e5..17c0f63 100644
--- a/test/test-lib.cpp
+++ b/test/test-lib.cpp
@@ -3,6 +3,7 @@
 #include "testHelpers.h"
 #include <boost/test/data/test_case.hpp>
 #include <boost/test/unit_test.hpp>
+#include <stream_support.h>
 
 #include <collections.h>
 #include <glArrays.h>
@@ -66,6 +67,46 @@ BOOST_AUTO_TEST_CASE(triangle_strip_iter)
 		out.push_back(c);
 	});
 	BOOST_REQUIRE_EQUAL(out.size(), (TRIANGLE_STRIP_IN.size() - 2) * 3);
-	BOOST_CHECK_EQUAL_COLLECTIONS(
-			out.begin(), out.end(), TRIANGLE_STRIP_EXPECTED.begin(), TRIANGLE_STRIP_EXPECTED.end());
+	BOOST_CHECK_EQUAL_COLCOL(out, TRIANGLE_STRIP_EXPECTED);
+}
+
+BOOST_AUTO_TEST_CASE(triangle_strip_range_adapter)
+{
+	using TriTuple = std::tuple<int, int, int>;
+	std::vector<TriTuple> outRange;
+	std::ranges::copy(TRIANGLE_STRIP_IN | triangleTriples, std::back_inserter(outRange));
+	constexpr std::array<TriTuple, 4> TRIANGLE_STRIP_EXPECTED_TUPLES {{{0, 1, 2}, {2, 1, 3}, {2, 3, 4}, {4, 3, 5}}};
+	BOOST_CHECK_EQUAL_COLCOL(outRange, TRIANGLE_STRIP_EXPECTED_TUPLES);
+}
+
+using MergeCloseData = std::tuple<std::vector<int>, int, std::vector<int>>;
+
+BOOST_DATA_TEST_CASE(mergeCloseInts,
+		boost::unit_test::data::make<MergeCloseData>({
+				{{0}, 0, {0}},
+				{{0, 1}, 0, {0, 1}},
+				{{0, 1}, 2, {0, 1}},
+				{{0, 1, 2}, 2, {0, 2}},
+				{{0, 1, 4}, 2, {0, 4}},
+				{{0, 1, 2}, 4, {0, 2}},
+				{{0, 4, 8}, 4, {0, 8}},
+				{{0, 4, 10, 14}, 4, {0, 14}},
+				{{0, 3, 6}, 2, {0, 3, 6}},
+				{{0, 3, 4}, 2, {0, 4}},
+				{{0, 5, 7, 12}, 4, {0, 6, 12}},
+				{{0, 3, 4, 5, 10, 17, 18, 19}, 2, {0, 4, 10, 19}},
+		}),
+		collection, tolerance, expected)
+{
+	auto mutableCollection {collection};
+	BOOST_REQUIRE_NO_THROW((mergeClose(
+			mutableCollection,
+			[](int left, int right) {
+				return std::abs(left - right);
+			},
+			[](int left, int right) {
+				return (left + right) / 2;
+			},
+			tolerance)));
+	BOOST_CHECK_EQUAL_COLCOL(mutableCollection, expected);
 }
diff --git a/test/test-maths.cpp b/test/test-maths.cpp
index b9d08bb..a6b780a 100644
--- a/test/test-maths.cpp
+++ b/test/test-maths.cpp
@@ -12,6 +12,7 @@
 #include <gfx/gl/camera.h>
 #include <glm/glm.hpp>
 #include <maths.h>
+#include <triangle.h>
 #include <tuple>
 
 using vecter_and_angle = std::tuple<glm::vec3, float>;
@@ -341,3 +342,98 @@ static_assert(linesIntersectAt(GlobalPosition2D {311000100, 491100100}, {3110500
 					  .value()
 		== GlobalPosition2D {311000100, 491100100});
 static_assert(!linesIntersectAt(glm::dvec2 {0, 1}, {0, 4}, {1, 8}, {1, 4}).has_value());
+
+BOOST_AUTO_TEST_CASE(triangle2d_helpers)
+{
+	constexpr static Triangle<2, float> t {{0, 0}, {5, 0}, {5, 5}};
+
+	BOOST_CHECK_CLOSE(t.angle(0), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({0, 0}), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angle(1), half_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({5, 0}), half_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angle(2), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({5, 5}), quarter_pi, 0.01F);
+
+	BOOST_CHECK_CLOSE(t.angleAt({0, 1}), 0.F, 0.01F);
+
+	BOOST_CHECK_CLOSE(t.area(), 12.5F, 0.01F);
+}
+
+BOOST_AUTO_TEST_CASE(triangle3d_helpers)
+{
+	constexpr static Triangle<3, float> t {{0, 0, 0}, {5, 0, 0}, {5, 5, 0}};
+
+	BOOST_CHECK_EQUAL(t.nnormal(), up);
+	BOOST_CHECK_CLOSE(t.angle(0), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({0, 0, 0}), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angle(1), half_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({5, 0, 0}), half_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angle(2), quarter_pi, 0.01F);
+	BOOST_CHECK_CLOSE(t.angleAt({5, 5, 0}), quarter_pi, 0.01F);
+
+	BOOST_CHECK_CLOSE(t.angleAt({0, 1, 0}), 0.F, 0.01F);
+
+	BOOST_CHECK_CLOSE(t.area(), 12.5F, 0.01F);
+}
+
+using ArcLineIntersectExp = std::pair<GlobalPosition2D, Angle>;
+using ArcLineIntersectData = std::tuple<GlobalPosition2D, GlobalPosition2D, GlobalPosition2D, GlobalPosition2D,
+		GlobalPosition2D, std::optional<ArcLineIntersectExp>>;
+
+BOOST_DATA_TEST_CASE(arcline_intersection,
+		boost::unit_test::data::make<ArcLineIntersectData>({
+				{{0, 0}, {0, 100}, {100, 0}, {200, 0}, {0, 200}, std::nullopt},
+				{{0, 0}, {0, 100}, {100, 0}, {0, 0}, {10, 10}, std::nullopt},
+				{{0, 0}, {0, 100}, {100, 0}, {0, 0}, {100, 100}, ArcLineIntersectExp {{71, 71}, quarter_pi}},
+				{{15, 27}, {15, 127}, {115, 27}, {15, 27}, {115, 127}, ArcLineIntersectExp {{86, 98}, quarter_pi}},
+				{{0, 0}, {0, 100}, {100, 0}, {0, 0}, {-100, -100}, std::nullopt},
+				{{0, 0}, {0, 100}, {100, 0}, {-10, 125}, {125, -10}, ArcLineIntersectExp {{16, 99}, 0.164F}},
+				{{0, 0}, {0, 100}, {100, 0}, {125, -10}, {-10, 125}, ArcLineIntersectExp {{99, 16}, 1.407F}},
+				{{0, 0}, {0, 100}, {100, 0}, {10, 125}, {125, -10}, ArcLineIntersectExp {{38, 93}, 0.385F}},
+				{{0, 0}, {0, 100}, {100, 0}, {12, 80}, {125, -10}, ArcLineIntersectExp {{99, 10}, 1.467F}},
+				{{0, 0}, {0, 100}, {100, 0}, {40, 80}, {125, -10}, ArcLineIntersectExp {{98, 18}, 1.387F}},
+				{{0, 0}, {0, 100}, {100, 0}, {40, 80}, {80, 20}, std::nullopt},
+				{{0, 0}, {0, 100}, {100, 0}, {40, 80}, {80, 80}, ArcLineIntersectExp {{60, 80}, 0.6435F}},
+				{{0, 0}, {0, 100}, {100, 0}, {80, 40}, {80, 80}, ArcLineIntersectExp {{80, 60}, 0.9273F}},
+				{{310002000, 490203000}, {310202000, 490203000}, {310002000, 490003000}, {310200000, 490150000},
+						{310150000, 490150000}, ArcLineIntersectExp {{310194850, 490150000}, 1.839F}},
+		}),
+		centre, arcStart, arcEnd, lineStart, lineEnd, expected)
+{
+	const ArcSegment arc {centre, arcStart, arcEnd};
+	BOOST_TEST_INFO(arc.first);
+	BOOST_TEST_INFO(arc.second);
+	BOOST_TEST_INFO(arc.length());
+
+	const auto intersection = arc.crossesLineAt(lineStart, lineEnd);
+	BOOST_REQUIRE_EQUAL(expected.has_value(), intersection.has_value());
+	if (expected.has_value()) {
+		BOOST_CHECK_EQUAL(expected->first, intersection->first);
+		BOOST_CHECK_CLOSE(expected->second, intersection->second, 1.F);
+	}
+}
+
+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}));
+
+static_assert(linesCross({1, 1}, {2, 2}, {1, 2}, {2, 1}));
+static_assert(linesCross({2, 2}, {1, 1}, {1, 2}, {2, 1}));
+
+static_assert(linesCrossLtR({1, 1}, {2, 2}, {1, 2}, {2, 1}));
+static_assert(!linesCrossLtR({2, 2}, {1, 1}, {1, 2}, {2, 1}));
+
+static_assert(Triangle<3, GlobalDistance> {{1, 2, 3}, {1, 0, 1}, {-2, 1, 0}}.positionOnPlane({7, -2})
+		== GlobalPosition3D {7, -2, 3});
+static_assert(Triangle<3, GlobalDistance> {
+					  {310000000, 490000000, 32800}, {310050000, 490050000, 33000}, {310000000, 490050000, 32700}}
+					  .positionOnPlane({310000000, 490000000})
+		== GlobalPosition3D {310000000, 490000000, 32800});
+static_assert(Triangle<3, GlobalDistance> {
+					  {310750000, 490150000, 58400}, {310800000, 490200000, 55500}, {310750000, 490200000, 57600}}
+					  .positionOnPlane({310751000, 490152000})
+		== GlobalPosition3D {310751000, 490152000, 58326});
diff --git a/test/test-network.cpp b/test/test-network.cpp
index 59eebae..e7419b5 100644
--- a/test/test-network.cpp
+++ b/test/test-network.cpp
@@ -47,8 +47,8 @@ struct TestLinkS : public TestLink, public LinkStraight {
 	}
 };
 
-constexpr GlobalPosition3D p000 {0, 0, 0}, p100 {10000, 0, 0}, p200 {20000, 0, 0}, p300 {30000, 0, 0};
-constexpr GlobalPosition3D p110 {10000, 10000, 0};
+constexpr GlobalPosition3D p000 {0, 0, 500}, p100 {10500, 0, 1000}, p200 {20100, 0, 2000}, p300 {30700, 0, 3000};
+constexpr GlobalPosition3D p110 {10300, 10400, 4000};
 
 template<> NetworkLinkHolder<TestLinkS>::NetworkLinkHolder() = default;
 
@@ -73,6 +73,18 @@ struct TestNetwork : public NetworkOf<TestLink, TestLinkS> {
 	render(const SceneShader &) const override
 	{
 	}
+
+	const Surface *
+	getBaseSurface() const override
+	{
+		return nullptr;
+	}
+
+	RelativeDistance
+	getBaseWidth() const override
+	{
+		return 5'700;
+	}
 };
 
 const auto VALID_NODES = boost::unit_test::data::make<GlobalPosition3D>({
@@ -208,6 +220,14 @@ BOOST_AUTO_TEST_CASE(routeTo_downStream_3to300)
 
 BOOST_AUTO_TEST_SUITE_END()
 
+namespace std {
+	std::ostream &
+	operator<<(std::ostream & s, const Link::End & e)
+	{
+		return s << std::format("End[dir: {}, loc: ({}, {}, {})]", e.dir, e.node->pos.x, e.node->pos.y, e.node->pos.z);
+	}
+}
+
 BOOST_FIXTURE_TEST_CASE(test_rail_network, RailLinks)
 {
 	//       0        1        2
@@ -221,7 +241,7 @@ BOOST_FIXTURE_TEST_CASE(test_rail_network, RailLinks)
 	//    --------
 	auto l0 = addLinksBetween(p000, p100);
 	BOOST_CHECK(dynamic_cast<RailLinkStraight *>(l0.get()));
-	BOOST_CHECK_EQUAL(l0->length, 10000);
+	BOOST_CHECK_EQUAL(l0->length, ::distance(p000, p100));
 	BOOST_CHECK_CLOSE(l0->ends[0].dir, half_pi, 0.1F);
 	BOOST_CHECK_CLOSE(l0->ends[1].dir, -half_pi, 0.1F);
 	BOOST_CHECK(l0->ends[0].nexts.empty());
@@ -229,7 +249,7 @@ BOOST_FIXTURE_TEST_CASE(test_rail_network, RailLinks)
 
 	auto l1 = addLinksBetween(p200, p100);
 	BOOST_CHECK(dynamic_cast<RailLinkStraight *>(l1.get()));
-	BOOST_CHECK_EQUAL(l1->length, 10000);
+	BOOST_CHECK_EQUAL(l1->length, ::distance(p200, p100));
 	BOOST_CHECK_CLOSE(l1->ends[0].dir, half_pi, 0.1F);
 	BOOST_CHECK_CLOSE(l1->ends[1].dir, -half_pi, 0.1F);
 	BOOST_CHECK(l0->ends[0].nexts.empty());
@@ -241,7 +261,7 @@ BOOST_FIXTURE_TEST_CASE(test_rail_network, RailLinks)
 
 	auto l2 = addLinksBetween(p200, p300);
 	BOOST_CHECK(dynamic_cast<RailLinkStraight *>(l2.get()));
-	BOOST_CHECK_EQUAL(l2->length, 10000);
+	BOOST_CHECK_EQUAL(l2->length, ::distance(p200, p300));
 	BOOST_CHECK_CLOSE(l2->ends[0].dir, half_pi, 0.1F);
 	BOOST_CHECK_CLOSE(l2->ends[1].dir, -half_pi, 0.1F);
 	BOOST_CHECK(l0->ends[0].nexts.empty());
@@ -251,19 +271,28 @@ BOOST_FIXTURE_TEST_CASE(test_rail_network, RailLinks)
 	BOOST_CHECK_EQUAL(l2->ends[0].nexts.at(0).second, 1);
 	BOOST_CHECK(l2->ends[1].nexts.empty());
 
-	auto l3 = addLinksBetween(p000, p110);
-	BOOST_CHECK(dynamic_cast<RailLinkCurve *>(l3.get()));
-	BOOST_CHECK_CLOSE(l3->length, (pi + half_pi) * 10000.F, 0.1F);
-	BOOST_CHECK_CLOSE(l3->ends[0].dir, -half_pi, 0.1F);
-	BOOST_CHECK_CLOSE(l3->ends[1].dir, 0, 0.1F);
-	BOOST_CHECK_EQUAL(l0->ends[0].nexts.at(0).first.lock(), l3);
-	BOOST_CHECK_EQUAL(l0->ends[0].nexts.at(0).second, 0);
-	BOOST_CHECK_EQUAL(l3->ends[0].nexts.at(0).first.lock(), l0);
-	BOOST_CHECK_EQUAL(l3->ends[0].nexts.at(0).second, 0);
-	BOOST_CHECK(l3->ends[1].nexts.empty());
-
-	auto l4 = addLinksBetween(p110, p300);
-	BOOST_CHECK_CLOSE(l4->length, 30400.F, 0.1F);
-	BOOST_CHECK_BETWEEN(l4->ends[0].dir, .23F, .24F);
-	BOOST_CHECK_CLOSE(l4->ends[1].dir, half_pi, 0.1F);
+	BOOST_CHECK_IF(l3, addLinksBetween(p000, p110)) {
+		BOOST_CHECK_IF(l3c, dynamic_cast<RailLinkCurve *>(l3.get())) {
+			BOOST_CHECK_CLOSE(l3c->radius, 10'300.F, 0.1F);
+			BOOST_CHECK_CLOSE(l3c->arc.length(), pi + half_pi, 0.5F);
+			BOOST_CHECK_CLOSE(l3->length, 48'563.F, 0.1F);
+			BOOST_CHECK_CLOSE(l3->ends[0].dir, -half_pi, 0.5F);
+			BOOST_CHECK_CLOSE(l3->ends[1].dir, -0.0097F, 0.5F);
+			BOOST_CHECK_EQUAL(l0->ends[0].nexts.at(0).first.lock(), l3);
+			BOOST_CHECK_EQUAL(l0->ends[0].nexts.at(0).second, 0);
+			BOOST_CHECK_EQUAL(l3->ends[0].nexts.at(0).first.lock(), l0);
+			BOOST_CHECK_EQUAL(l3->ends[0].nexts.at(0).second, 0);
+			BOOST_CHECK(l3->ends[1].nexts.empty());
+		}
+	}
+
+	BOOST_CHECK_IF(l4, addLinksBetween(p110, p300)) {
+		BOOST_CHECK_IF(l4c, dynamic_cast<RailLinkCurve *>(l4.get())) {
+			BOOST_CHECK_CLOSE(l4c->radius, 6950.F, 0.1F);
+			BOOST_CHECK_CLOSE(l4c->arc.length(), 4.456F, 0.1F);
+			BOOST_CHECK_CLOSE(l4->length, 30'981.F, 0.1F);
+			BOOST_CHECK_BETWEEN(l4->ends[0].dir, .25F, .26F);
+			BOOST_CHECK_CLOSE(l4->ends[1].dir, half_pi, 0.1F);
+		}
+	}
 }
diff --git a/test/test-render.cpp b/test/test-render.cpp
index 3966f28..3c453bd 100644
--- a/test/test-render.cpp
+++ b/test/test-render.cpp
@@ -30,11 +30,10 @@ class TestScene : public SceneProvider {
 	RailVehicleClassPtr brush47rvc;
 	std::shared_ptr<RailVehicle> train1, train2;
 	RailLinks rail;
-	std::shared_ptr<GeoData> gd = std::make_shared<GeoData>(GeoData::createFlat({0, 0}, {1000000, 1000000}, 1));
 	std::shared_ptr<Environment> env = std::make_shared<Environment>();
 
-	Terrain terrain {gd};
-	Water water {gd};
+	std::shared_ptr<Terrain> terrain = std::make_shared<Terrain>(GeoData::createFlat({0, 0}, {1000000, 1000000}, 1));
+	Water water {terrain};
 
 public:
 	TestScene()
@@ -71,7 +70,7 @@ public:
 	void
 	content(const SceneShader & shader) const override
 	{
-		terrain.render(shader);
+		terrain->render(shader);
 		water.render(shader);
 		rail.render(shader);
 		std::ranges::for_each(gameState->assets, [&shader](const auto & asset) {
@@ -95,7 +94,7 @@ public:
 	void
 	shadows(const ShadowMapper & shadowMapper) const override
 	{
-		terrain.shadows(shadowMapper);
+		terrain->shadows(shadowMapper);
 		std::ranges::for_each(gameState->assets, [&shadowMapper](const auto & asset) {
 			if (const auto renderable = std::dynamic_pointer_cast<const Renderable>(asset.second)) {
 				renderable->shadows(shadowMapper);
@@ -167,15 +166,14 @@ BOOST_AUTO_TEST_CASE(terrain)
 	ss.camera.setView({310000000, 490000000, 600000}, glm::normalize(glm::vec3 {1, 1, -0.5F}));
 
 	class TestTerrain : public SceneProvider {
-		std::shared_ptr<GeoData> gd
-				= std::make_shared<GeoData>(GeoData::loadFromAsciiGrid(FIXTURESDIR "height/SD19.asc"));
-		Terrain terrain {gd};
-		Water water {gd};
+		std::shared_ptr<Terrain> terrain
+				= std::make_shared<Terrain>(GeoData::loadFromAsciiGrid(FIXTURESDIR "height/SD19.asc"));
+		Water water {terrain};
 
 		void
 		content(const SceneShader & shader) const override
 		{
-			terrain.render(shader);
+			terrain->render(shader);
 			water.render(shader);
 		}
 
@@ -194,7 +192,7 @@ BOOST_AUTO_TEST_CASE(terrain)
 		void
 		shadows(const ShadowMapper & shadowMapper) const override
 		{
-			terrain.shadows(shadowMapper);
+			terrain->shadows(shadowMapper);
 		}
 	};
 
diff --git a/thirdparty/openmesh/helpers.h b/thirdparty/openmesh/helpers.h
new file mode 100644
index 0000000..d148c06
--- /dev/null
+++ b/thirdparty/openmesh/helpers.h
@@ -0,0 +1,28 @@
+#pragma once
+#include <OpenMesh/Core/Mesh/BaseKernel.hh>
+#include <OpenMesh/Core/Mesh/PolyConnectivity.hh>
+#include <ranges>
+
+namespace OpenMesh {
+	template<typename Iter, typename CenterEntityHandle>
+	using IteratorFunction = Iter (OpenMesh::PolyConnectivity::*)(CenterEntityHandle) const;
+
+	template<typename Iter, typename CenterEntityHandle, IteratorFunction<Iter, CenterEntityHandle> BeginFunc,
+			IteratorFunction<Iter, CenterEntityHandle> EndFunc, typename Adaptor>
+	auto
+	operator|(const OpenMesh::PolyConnectivity::CirculatorRange<OpenMesh::PolyConnectivity, Iter, CenterEntityHandle,
+					  BeginFunc, EndFunc> & range,
+			Adaptor && adaptor)
+	{
+		return std::views::iota(range.begin(), range.end()) | std::forward<Adaptor>(adaptor);
+	}
+
+	namespace Helpers {
+		template<typename Type, template<typename> typename PropertyT> struct Property : public PropertyT<Type> {
+			template<typename... Params> explicit Property(OpenMesh::BaseKernel * kernel, Params &&... params)
+			{
+				kernel->add_property(*this, std::forward<Params>(params)...);
+			}
+		};
+	}
+}
diff --git a/ui/builders/freeExtend.cpp b/ui/builders/freeExtend.cpp
index db127e6..ab5a998 100644
--- a/ui/builders/freeExtend.cpp
+++ b/ui/builders/freeExtend.cpp
@@ -38,11 +38,11 @@ BuilderFreeExtend::click(
 		case SDL_BUTTON_LEFT:
 			if (p1) {
 				if (const auto p = network->intersectRayNodes(ray)) {
-					network->addJoins(*p1, p->pos);
+					createJoin(network, geoData, *p1, p->pos);
 					p1 = p->pos;
 				}
 				else if (const auto p = geoData->intersectRay(ray)) {
-					network->addExtend(*p1, p->first);
+					createExtend(network, geoData, *p1, p->first);
 					p1 = p->first;
 				}
 			}
@@ -57,3 +57,21 @@ BuilderFreeExtend::click(
 			return;
 	}
 }
+
+Link::CCollection
+BuilderFreeExtend::createJoin(
+		Network * network, const GeoData * geoData, GlobalPosition3D p1, GlobalPosition3D p2) const
+{
+	const auto links = network->addJoins(geoData, p1, p2);
+	setHeightsFor(network, links);
+	return links;
+}
+
+Link::CCollection
+BuilderFreeExtend::createExtend(
+		Network * network, const GeoData * geoData, GlobalPosition3D p1, GlobalPosition3D p2) const
+{
+	const auto links = network->addExtend(geoData, p1, p2);
+	setHeightsFor(network, links);
+	return links;
+}
diff --git a/ui/builders/freeExtend.h b/ui/builders/freeExtend.h
index 0d5f327..6f28493 100644
--- a/ui/builders/freeExtend.h
+++ b/ui/builders/freeExtend.h
@@ -5,11 +5,17 @@ class Network;
 class GeoData;
 
 class BuilderFreeExtend : public EditNetwork::Builder {
+private:
 	std::string hint() const override;
 	void click(Network * network, const GeoData * geoData, const SDL_MouseButtonEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 	void move(Network * network, const GeoData * geoData, const SDL_MouseMotionEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 
+public:
+	Link::CCollection createJoin(Network * network, const GeoData *, GlobalPosition3D, GlobalPosition3D) const;
+	Link::CCollection createExtend(Network * network, const GeoData *, GlobalPosition3D, GlobalPosition3D) const;
+
+private:
 	std::optional<GlobalPosition3D> p1;
 };
diff --git a/ui/builders/join.cpp b/ui/builders/join.cpp
index 7474c5b..6941e23 100644
--- a/ui/builders/join.cpp
+++ b/ui/builders/join.cpp
@@ -25,13 +25,13 @@ BuilderJoin::move(Network * network, const GeoData *, const SDL_MouseMotionEvent
 
 void
 BuilderJoin::click(
-		Network * network, const GeoData *, const SDL_MouseButtonEvent & e, const Ray<GlobalPosition3D> & ray)
+		Network * network, const GeoData * geoData, const SDL_MouseButtonEvent & e, const Ray<GlobalPosition3D> & ray)
 {
 	switch (e.button) {
 		case SDL_BUTTON_LEFT:
 			if (const auto p = network->intersectRayNodes(ray)) {
 				if (p1) {
-					create(network, p1, p);
+					create(network, geoData, p1, p);
 					p1.reset();
 					candidateLinks.removeAll();
 				}
@@ -46,8 +46,10 @@ BuilderJoin::click(
 	}
 }
 
-void
-BuilderJoin::create(Network * network, const Node::Ptr & p1, const Node::Ptr & p2) const
+Link::CCollection
+BuilderJoin::create(Network * network, const GeoData * geoData, const Node::Ptr & p1, const Node::Ptr & p2) const
 {
-	network->addJoins(p1->pos, p2->pos);
+	const auto links = network->addJoins(geoData, p1->pos, p2->pos);
+	setHeightsFor(network, links);
+	return links;
 }
diff --git a/ui/builders/join.h b/ui/builders/join.h
index dd57895..326d23d 100644
--- a/ui/builders/join.h
+++ b/ui/builders/join.h
@@ -5,13 +5,14 @@ class Network;
 class GeoData;
 
 class BuilderJoin : public EditNetwork::Builder {
+private:
 	std::string hint() const override;
 	void click(Network * network, const GeoData * geoData, const SDL_MouseButtonEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 	void move(Network * network, const GeoData * geoData, const SDL_MouseMotionEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 
-	void create(Network * network, const Node::Ptr & p1, const Node::Ptr & p2) const;
+	Link::CCollection create(Network * network, const GeoData *, const Node::Ptr & p1, const Node::Ptr & p2) const;
 
 	Node::Ptr p1;
 };
diff --git a/ui/builders/straight.cpp b/ui/builders/straight.cpp
index 43f5ec8..338aa8a 100644
--- a/ui/builders/straight.cpp
+++ b/ui/builders/straight.cpp
@@ -1,4 +1,5 @@
 #include "straight.h"
+#include "stream_support.h"
 #include <game/geoData.h>
 
 std::string
@@ -32,7 +33,7 @@ BuilderStraight::click(
 		case SDL_BUTTON_LEFT:
 			if (const auto p = geoData->intersectRay(ray)) {
 				if (p1) {
-					create(network, *p1, p->first);
+					create(network, geoData, *p1, p->first);
 					candidateLinks.removeAll();
 					p1.reset();
 				}
@@ -48,8 +49,10 @@ BuilderStraight::click(
 	}
 }
 
-void
-BuilderStraight::create(Network * network, GlobalPosition3D p1, GlobalPosition3D p2) const
+Link::CCollection
+BuilderStraight::create(Network * network, const GeoData * geoData, GlobalPosition3D p1, GlobalPosition3D p2) const
 {
-	network->addStraight(p1, p2);
+	const auto links = network->addStraight(geoData, p1, p2);
+	setHeightsFor(network, links);
+	return links;
 }
diff --git a/ui/builders/straight.h b/ui/builders/straight.h
index 28eb66e..0a6f290 100644
--- a/ui/builders/straight.h
+++ b/ui/builders/straight.h
@@ -5,13 +5,16 @@ class Network;
 class GeoData;
 
 class BuilderStraight : public EditNetwork::Builder {
+private:
 	std::string hint() const override;
 	void click(Network * network, const GeoData * geoData, const SDL_MouseButtonEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 	void move(Network * network, const GeoData * geoData, const SDL_MouseMotionEvent & e,
 			const Ray<GlobalPosition3D> & ray) override;
 
-	void create(Network * network, GlobalPosition3D p1, GlobalPosition3D p2) const;
+public:
+	Link::CCollection create(Network * network, const GeoData *, GlobalPosition3D p1, GlobalPosition3D p2) const;
 
+private:
 	std::optional<GlobalPosition3D> p1;
 };
diff --git a/ui/editNetwork.cpp b/ui/editNetwork.cpp
index ac2d93d..c4c0297 100644
--- a/ui/editNetwork.cpp
+++ b/ui/editNetwork.cpp
@@ -4,7 +4,7 @@
 #include "builders/straight.h"
 #include "text.h"
 #include <game/gamestate.h>
-#include <game/geoData.h>
+#include <game/terrain.h>
 #include <gfx/gl/sceneShader.h>
 #include <gfx/models/texture.h>
 
@@ -26,7 +26,7 @@ bool
 EditNetwork::click(const SDL_MouseButtonEvent & e, const Ray<GlobalPosition3D> & ray)
 {
 	if (builder && (e.button == SDL_BUTTON_LEFT || e.button == SDL_BUTTON_MIDDLE)) {
-		builder->click(network, gameState->geoData.get(), e, ray);
+		builder->click(network, gameState->terrain.get(), e, ray);
 		return true;
 	}
 	return false;
@@ -36,7 +36,7 @@ bool
 EditNetwork::move(const SDL_MouseMotionEvent & e, const Ray<GlobalPosition3D> & ray)
 {
 	if (builder) {
-		builder->move(network, gameState->geoData.get(), e, ray);
+		builder->move(network, gameState->terrain.get(), e, ray);
 	}
 	return false;
 }
@@ -64,6 +64,17 @@ EditNetwork::Builder::render(const SceneShader & shader) const
 }
 
 void
+EditNetwork::Builder::setHeightsFor(Network * network, const Link::CCollection & links, GeoData::SetHeightsOpts opts)
+{
+	opts.surface = network->getBaseSurface();
+	const auto width = network->getBaseWidth();
+
+	for (const auto & link : links) {
+		gameState->terrain->setHeights(link->getBase(width), opts);
+	}
+}
+
+void
 EditNetwork::render(const UIShader & shader, const UIComponent::Position & parentPos) const
 {
 	if (builder) {
diff --git a/ui/editNetwork.h b/ui/editNetwork.h
index ec06fa7..2ae467d 100644
--- a/ui/editNetwork.h
+++ b/ui/editNetwork.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include "game/geoData.h"
 #include "gameMainSelector.h"
 #include "modeHelper.h"
 #include "toolbar.h"
@@ -30,6 +31,8 @@ public:
 		virtual void click(Network *, const GeoData *, const SDL_MouseButtonEvent &, const Ray<GlobalPosition3D> &) = 0;
 		virtual void move(Network *, const GeoData *, const SDL_MouseMotionEvent &, const Ray<GlobalPosition3D> &) = 0;
 
+		static void setHeightsFor(Network *, const Link::CCollection &, GeoData::SetHeightsOpts = {});
+
 		using Ptr = std::unique_ptr<Builder>;
 
 	protected:
diff --git a/ui/gameMainSelector.cpp b/ui/gameMainSelector.cpp
index 5bef48d..e9642ec 100644
--- a/ui/gameMainSelector.cpp
+++ b/ui/gameMainSelector.cpp
@@ -4,8 +4,8 @@
 #include "ui/uiComponent.h"
 #include <SDL2/SDL.h>
 #include <game/gamestate.h>
-#include <game/geoData.h>
 #include <game/selectable.h>
+#include <game/terrain.h>
 #include <game/worldobject.h> // IWYU pragma: keep
 #include <gfx/gl/camera.h>
 #include <optional>
@@ -83,7 +83,7 @@ GameMainSelector::defaultClick(const Ray<GlobalPosition3D> & ray)
 		const auto & ref = *selected.base()->get();
 		clicked = typeid(ref).name();
 	}
-	else if (const auto pos = gameState->geoData->intersectRay(ray)) {
+	else if (const auto pos = gameState->terrain->intersectRay(ray)) {
 		clicked = streamed_string(*pos);
 	}
 	else {