Merge remote-tracking branch 'origin/deform-terrain'

Two related issues remain:
* Terrain self shadowing is common and handled poorly
* Odd, but mathematically correct patterns/stripes in feature boundaries

Neither of these relate directly to deformation.

21 files changed, 954 insertions, 96 deletions

diff --git a/Jamroot.jam b/Jamroot.jam
index de27053..7a2fe16 100644
--- a/Jamroot.jam
+++ b/Jamroot.jam
@@ -17,6 +17,7 @@ lib OpenMeshCore ;
 variant coverage : debug ;
 project : requirements
 	<cxxstd>20
+	<linkflags>-Wl,-z,defs
 	<variant>release:<lto>on
 	<variant>profile:<lto>on
 	<variant>coverage:<coverage>on
diff --git a/game/geoData.cpp b/game/geoData.cpp
index 359e8c0..ed4303b 100644
--- a/game/geoData.cpp
+++ b/game/geoData.cpp
@@ -1,7 +1,10 @@
 #include "geoData.h"
+#include "collections.h"
+#include "geometricPlane.h"
 #include <fstream>
 #include <glm/gtx/intersect.hpp>
 #include <maths.h>
+#include <set>
 
 GeoData
 GeoData::loadFromAsciiGrid(const std::filesystem::path & input)
@@ -62,6 +65,8 @@ GeoData::loadFromAsciiGrid(const std::filesystem::path & input)
 	return mesh;
 };
 
+template<typename T> constexpr static T GRID_SIZE = 10'000;
+
 GeoData
 GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistance h)
 {
@@ -70,11 +75,29 @@ GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistan
 	mesh.lowerExtent = {lower, h};
 	mesh.upperExtent = {upper, h};
 
-	const auto ll = mesh.add_vertex({lower.x, lower.y, h}), lu = mesh.add_vertex({lower.x, upper.y, h}),
-			   ul = mesh.add_vertex({upper.x, lower.y, h}), uu = mesh.add_vertex({upper.x, upper.y, h});
+	std::vector<VertexHandle> vertices;
+	for (GlobalDistance row = lower.x; row < upper.x; row += GRID_SIZE<GlobalDistance>) {
+		for (GlobalDistance col = lower.y; col < upper.y; col += GRID_SIZE<GlobalDistance>) {
+			vertices.push_back(mesh.add_vertex({col, row, h}));
+		}
+	}
 
-	mesh.add_face(ll, uu, lu);
-	mesh.add_face(ll, ul, uu);
+	const auto nrows = static_cast<size_t>(std::ceil(float(upper.x - lower.x) / GRID_SIZE<RelativeDistance>));
+	const auto ncols = static_cast<size_t>(std::ceil(float(upper.y - lower.y) / GRID_SIZE<RelativeDistance>));
+	for (size_t row = 1; row < nrows; ++row) {
+		for (size_t col = 1; col < ncols; ++col) {
+			mesh.add_face({
+					vertices[ncols * (row - 1) + (col - 1)],
+					vertices[ncols * (row - 0) + (col - 0)],
+					vertices[ncols * (row - 0) + (col - 1)],
+			});
+			mesh.add_face({
+					vertices[ncols * (row - 1) + (col - 1)],
+					vertices[ncols * (row - 1) + (col - 0)],
+					vertices[ncols * (row - 0) + (col - 0)],
+			});
+		}
+	}
 
 	mesh.update_vertex_normals_only();
 
@@ -84,11 +107,11 @@ GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistan
 OpenMesh::FaceHandle
 GeoData::findPoint(GlobalPosition2D p) const
 {
-	return findPoint(p, *faces_begin());
+	return findPoint(p, *faces_sbegin());
 }
 
 GeoData::PointFace::PointFace(const GlobalPosition2D p, const GeoData * mesh) :
-	PointFace {p, mesh, *mesh->faces_begin()}
+	PointFace {p, mesh, *mesh->faces_sbegin()}
 {
 }
 
@@ -112,7 +135,7 @@ GeoData::PointFace::face(const GeoData * mesh, FaceHandle start) const
 GeoData::FaceHandle
 GeoData::PointFace::face(const GeoData * mesh) const
 {
-	return face(mesh, *mesh->faces_begin());
+	return face(mesh, *mesh->faces_sbegin());
 }
 
 namespace {
@@ -193,16 +216,16 @@ GeoData::positionAt(const PointFace & p) const
 	return positionOnTriangle(p.point, triangle<3>(p.face(this)));
 }
 
-[[nodiscard]] std::optional<GlobalPosition3D>
+[[nodiscard]] GeoData::IntersectionResult
 GeoData::intersectRay(const Ray<GlobalPosition3D> & ray) const
 {
 	return intersectRay(ray, findPoint(ray.start));
 }
 
-[[nodiscard]] std::optional<GlobalPosition3D>
+[[nodiscard]] GeoData::IntersectionResult
 GeoData::intersectRay(const Ray<GlobalPosition3D> & ray, FaceHandle face) const
 {
-	std::optional<GlobalPosition3D> out;
+	GeoData::IntersectionResult out;
 	walkUntil(PointFace {ray.start, face},
 			ray.start.xy() + (ray.direction.xy() * RelativePosition2D(upperExtent.xy() - lowerExtent.xy())),
 			[&out, &ray, this](FaceHandle face) {
@@ -210,7 +233,7 @@ GeoData::intersectRay(const Ray<GlobalPosition3D> & ray, FaceHandle face) const
 				RelativeDistance dist {};
 				const auto t = triangle<3>(face);
 				if (ray.intersectTriangle(t.x, t.y, t.z, bari, dist)) {
-					out = t * bari;
+					out.emplace(t * bari, face);
 					return true;
 				}
 				return false;
@@ -324,17 +347,321 @@ GeoData::triangleContainsPoint(const GlobalPosition2D p, FaceHandle face) const
 GeoData::HalfedgeHandle
 GeoData::findBoundaryStart() const
 {
-	return *std::find_if(halfedges_begin(), halfedges_end(), [this](const auto heh) {
+	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 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::update_vertex_normals_only()
 {
-	for (auto vh : all_vertices()) {
-		Normal3D n;
-		calc_vertex_normal_correct(vh, n);
-		this->set_normal(vh, glm::normalize(n));
+	update_vertex_normals_only(vertices_sbegin());
+}
+
+void
+GeoData::update_vertex_normals_only(VertexIter start)
+{
+	std::for_each(start, vertices_end(), [this](const auto vh) {
+		if (normal(vh) == Normal3D {}) {
+			Normal3D n;
+			calc_vertex_normal_correct(vh, n);
+			this->set_normal(vh, glm::normalize(n));
+		}
+	});
+}
+
+bool
+GeoData::triangleOverlapsTriangle(const Triangle<2> & a, const Triangle<2> & b)
+{
+	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);
+}
+
+std::array<GeoData::FaceHandle, 4>
+GeoData::split(FaceHandle _fh)
+{
+	// Collect halfedges of face
+	const HalfedgeHandle he0 = halfedge_handle(_fh);
+	const HalfedgeHandle he1 = next_halfedge_handle(he0);
+	const HalfedgeHandle he2 = next_halfedge_handle(he1);
+
+	const EdgeHandle eh0 = edge_handle(he0);
+	const EdgeHandle eh1 = edge_handle(he1);
+	const EdgeHandle eh2 = edge_handle(he2);
+
+	// Collect points of face
+	const VertexHandle p0 = to_vertex_handle(he0);
+	const VertexHandle p1 = to_vertex_handle(he1);
+	const VertexHandle p2 = to_vertex_handle(he2);
+
+	// Calculate midpoint coordinates
+	const Point new0 = centre(he0);
+	const Point new1 = centre(he1);
+	const Point new2 = centre(he2);
+
+	// Add vertices at midpoint coordinates
+	const VertexHandle v0 = add_vertex(new0);
+	const VertexHandle v1 = add_vertex(new1);
+	const VertexHandle v2 = add_vertex(new2);
+
+	const bool split0 = !is_boundary(eh0);
+	const bool split1 = !is_boundary(eh1);
+	const bool split2 = !is_boundary(eh2);
+
+	// delete original face
+	delete_face(_fh, false);
+
+	// split boundary edges of deleted face ( if not boundary )
+	if (split0) {
+		split(eh0, v0);
+	}
+
+	if (split1) {
+		split(eh1, v1);
+	}
+
+	if (split2) {
+		split(eh2, v2);
+	}
+
+	// Retriangulate
+	return {
+			add_face(v0, p0, v1),
+			add_face(p2, v0, v2),
+			add_face(v2, v1, p1),
+			add_face(v2, v0, v1),
+	};
+}
+
+void
+GeoData::setHeights(const std::span<const GlobalPosition3D> triangleStrip)
+{
+	static const RelativeDistance MAX_SLOPE = 1.5F;
+	static const RelativeDistance MIN_ARC = 0.01F;
+
+	if (triangleStrip.size() < 3) {
+		return;
+	}
+
+	const auto initialVertexCount = static_cast<unsigned int>(n_vertices());
+
+	// Create new vertices
+	std::vector<VertexHandle> newVerts;
+	newVerts.reserve(newVerts.size());
+	std::transform(triangleStrip.begin(), triangleStrip.end(), std::back_inserter(newVerts), [this](const auto tsVert) {
+		return add_vertex(tsVert);
+	});
+	//  Create new faces
+	std::for_each(strip_begin(newVerts), strip_end(newVerts), [this](const auto & newVert) {
+		const auto [a, b, c] = newVert;
+		add_face(a, b, c);
+	});
+	for (auto start = faces_sbegin(); std::any_of(start, faces_end(), [this, &start](const auto fh) {
+			 static constexpr auto MAX_FACE_AREA = 100'000'000.F;
+			 if (triangle<3>(fh).area() > MAX_FACE_AREA) {
+				 split(fh);
+				 start = FaceIter {*this, FaceHandle(fh), true};
+				 return true;
+			 }
+			 return false;
+		 });) {
+		;
 	}
+
+	// Extrude corners
+	struct Extrusion {
+		VertexHandle boundaryVertex, extrusionVertex;
+		Direction3D lowerLimit, upperLimit;
+	};
+
+	std::vector<Extrusion> extrusionExtents;
+	boundaryWalk(
+			[this, &extrusionExtents](const auto boundaryHeh) {
+				const auto prevBoundaryHeh = prev_halfedge_handle(boundaryHeh);
+				const auto prevBoundaryVertex = from_vertex_handle(prevBoundaryHeh);
+				const auto boundaryVertex = from_vertex_handle(boundaryHeh);
+				const auto nextBoundaryVertex = to_vertex_handle(boundaryHeh);
+				const auto p0 = point(prevBoundaryVertex);
+				const auto p1 = point(boundaryVertex);
+				const auto p2 = point(nextBoundaryVertex);
+				const auto e0 = glm::normalize(vector_normal(RelativePosition2D(p1 - p0)));
+				const auto e1 = glm::normalize(vector_normal(RelativePosition2D(p2 - p1)));
+
+				const auto addExtrusionFor = [this, &extrusionExtents, boundaryVertex, p1](Direction2D direction) {
+					const auto doExtrusion = [this](VertexHandle & extrusionVertex, Direction2D direction,
+													 GlobalPosition3D boundaryVertex, RelativeDistance vert) {
+						const auto extrusionDir = glm::normalize(direction || vert);
+
+						if (!extrusionVertex.is_valid()) {
+							if (const auto intersect = intersectRay({boundaryVertex, extrusionDir})) {
+								auto splitVertex = split(intersect->second, intersect->first);
+								extrusionVertex = splitVertex;
+							}
+							else if (const auto intersect
+									= intersectRay({boundaryVertex + GlobalPosition3D {1, 1, 0}, extrusionDir})) {
+								auto splitVertex = split(intersect->second, intersect->first);
+								extrusionVertex = splitVertex;
+							}
+							else if (const auto intersect
+									= intersectRay({boundaryVertex + GlobalPosition3D {1, 0, 0}, extrusionDir})) {
+								auto splitVertex = split(intersect->second, intersect->first);
+								extrusionVertex = splitVertex;
+							}
+						}
+
+						return extrusionDir;
+					};
+
+					VertexHandle extrusionVertex;
+					extrusionExtents.emplace_back(boundaryVertex, extrusionVertex,
+							doExtrusion(extrusionVertex, direction, p1, -MAX_SLOPE),
+							doExtrusion(extrusionVertex, direction, p1, MAX_SLOPE));
+					assert(extrusionVertex.is_valid());
+				};
+				if (const Arc arc {e0, e1}; arc.length() < MIN_ARC) {
+					addExtrusionFor(normalize(e0 + e1) / cosf(arc.length() / 2.F));
+				}
+				else if (arc.length() < pi) {
+					// Previous half edge end to current half end start arc tangents
+					const auto limit = std::ceil(arc.length() * 5.F / pi);
+					const auto inc = arc.length() / limit;
+					for (float step = 0; step <= limit; step += 1.F) {
+						addExtrusionFor(sincosf(arc.first + (step * inc)));
+					}
+				}
+				else {
+					// Single tangent bisecting the difference
+					addExtrusionFor(normalize(e0 + e1) / sinf((arc.length() - pi) / 2.F));
+				}
+			},
+			*voh_begin(newVerts.front()));
+
+	//  Cut existing terrain
+	extrusionExtents.emplace_back(extrusionExtents.front()); // Circular next
+	std::vector<std::vector<VertexHandle>> boundaryFaces;
+	std::adjacent_find(extrusionExtents.begin(), extrusionExtents.end(),
+			[this, &boundaryFaces](const auto & first, const auto & second) {
+				const auto p0 = point(first.boundaryVertex);
+				const auto p1 = point(second.boundaryVertex);
+				const auto bdir = RelativePosition3D(p1 - p0);
+				const auto make_plane = [p0](auto y, auto z) {
+					return GeometricPlaneT<GlobalPosition3D> {p0, crossProduct(y, z)};
+				};
+				const auto planes = ((first.boundaryVertex == second.boundaryVertex)
+						? std::array {make_plane(second.lowerLimit, first.lowerLimit),
+									make_plane(second.upperLimit, first.upperLimit),
+						}
+						: std::array {
+								  make_plane(bdir, second.lowerLimit),
+									make_plane(bdir, second.upperLimit),
+						});
+				assert(planes.front().normal.z > 0.F);
+				assert(planes.back().normal.z > 0.F);
+
+				auto & out = boundaryFaces.emplace_back();
+				out.emplace_back(first.boundaryVertex);
+				out.emplace_back(first.extrusionVertex);
+				for (auto currentVertex = first.extrusionVertex;
+						!find_halfedge(currentVertex, second.extrusionVertex).is_valid();) {
+					[[maybe_unused]] const auto n = std::any_of(
+							voh_begin(currentVertex), voh_end(currentVertex), [&](const auto currentVertexOut) {
+								const auto next = next_halfedge_handle(currentVertexOut);
+								const auto nextVertex = to_vertex_handle(next);
+								const auto startVertex = from_vertex_handle(next);
+								if (nextVertex == *++out.rbegin()) {
+									// This half edge goes back to the previous vertex
+									return false;
+								}
+								const auto edge = edge_handle(next);
+								const auto ep0 = point(startVertex);
+								const auto ep1 = point(nextVertex);
+								if (planes.front().getRelation(ep1) == GeometricPlane::PlaneRelation::Below
+										|| planes.back().getRelation(ep1) == GeometricPlane::PlaneRelation::Above) {
+									return false;
+								}
+								const auto diff = RelativePosition3D(ep1 - ep0);
+								const auto length = glm::length(diff);
+								const auto dir = diff / length;
+								const Ray r {ep1, -dir};
+								const auto dists = planes * [r](const auto & plane) {
+									RelativeDistance dist {};
+									if (r.intersectPlane(plane.origin, plane.normal, dist)) {
+										return dist;
+									}
+									return INFINITY;
+								};
+								const auto dist = *std::min_element(dists.begin(), dists.end());
+								const auto splitPos = ep1 - (dir * dist);
+								if (dist <= length) {
+									currentVertex = split(edge, splitPos);
+									out.emplace_back(currentVertex);
+									return true;
+								}
+								return false;
+							});
+					assert(n);
+				}
+				out.emplace_back(second.extrusionVertex);
+				if (first.boundaryVertex != second.boundaryVertex) {
+					out.emplace_back(second.boundaryVertex);
+				}
+				return false;
+			});
+
+	//  Remove old faces
+	std::set<FaceHandle> visited;
+	auto removeOld = [&](auto & self, const auto face) -> void {
+		if (visited.insert(face).second) {
+			std::for_each(fh_begin(face), fh_end(face), [&](const auto fh) {
+				const auto b1 = to_vertex_handle(fh);
+				const auto b2 = from_vertex_handle(fh);
+				if (opposite_face_handle(fh).is_valid()
+						&& std::none_of(boundaryFaces.begin(), boundaryFaces.end(), [b2, b1](const auto & bf) {
+							   return std::adjacent_find(bf.begin(), bf.end(), [b2, b1](const auto v1, const auto v2) {
+								   return b1 == v1 && b2 == v2;
+							   }) != bf.end();
+						   })) {
+					self(self, opposite_face_handle(fh));
+				}
+			});
+
+			delete_face(face, false);
+		}
+	};
+	removeOld(removeOld, findPoint(triangleStrip.front()));
+
+	std::for_each(boundaryFaces.begin(), boundaryFaces.end(), [&](auto & boundaryFace) {
+		std::reverse(boundaryFace.begin(), boundaryFace.end());
+		add_face(boundaryFace);
+	});
+
+	// Tidy up
+	update_vertex_normals_only(VertexIter {*this, vertex_handle(initialVertexCount), true});
 }
diff --git a/game/geoData.h b/game/geoData.h
index 15143e8..021b4c7 100644
--- a/game/geoData.h
+++ b/game/geoData.h
@@ -61,11 +61,69 @@ public:
 			});
 		}
 
+		[[nodiscard]]
 		glm::vec<Dim, GlobalDistance>
 		operator*(BaryPosition bari) const
 		{
-			const auto & t {*this};
-			return t[0] + (RelativePosition<Dim>(t[1] - t[0]) * bari.x) + (RelativePosition<Dim>(t[2] - t[1]) * bari.y);
+			return p(0) + (difference(p(0), p(1)) * bari.x) + (difference(p(0), p(2)) * bari.y);
+		}
+
+		[[nodiscard]]
+		auto
+		area() const
+			requires(Dim == 3)
+		{
+			return glm::length(crossProduct(difference(p(0), p(1)), difference(p(0), p(2)))) / 2.F;
+		}
+
+		[[nodiscard]]
+		Normal3D
+		normal() const
+			requires(Dim == 3)
+		{
+			return crossProduct(difference(p(0), p(1)), difference(p(0), p(2)));
+		}
+
+		[[nodiscard]]
+		Normal3D
+		nnormal() const
+			requires(Dim == 3)
+		{
+			return glm::normalize(normal());
+		}
+
+		[[nodiscard]]
+		auto
+		angle(glm::length_t c) const
+		{
+			return Arc {P(c), P(c + 2), P(c + 1)}.length();
+		}
+
+		template<glm::length_t D = Dim>
+		[[nodiscard]]
+		auto
+		angleAt(const GlobalPosition<D> pos) const
+			requires(D <= Dim)
+		{
+			for (glm::length_t i {}; i < 3; ++i) {
+				if (GlobalPosition<D> {p(i)} == pos) {
+					return angle(i);
+				}
+			}
+			return 0.F;
+		}
+
+		[[nodiscard]]
+		inline auto
+		p(const glm::length_t i) const
+		{
+			return base::operator[](i);
+		}
+
+		[[nodiscard]] inline auto
+		P(const glm::length_t i) const
+		{
+			return base::operator[](i % 3);
 		}
 	};
 
@@ -73,8 +131,10 @@ public:
 	[[nodiscard]] FaceHandle findPoint(GlobalPosition2D, FaceHandle start) const;
 
 	[[nodiscard]] GlobalPosition3D positionAt(const PointFace &) const;
-	[[nodiscard]] std::optional<GlobalPosition3D> intersectRay(const Ray<GlobalPosition3D> &) const;
-	[[nodiscard]] std::optional<GlobalPosition3D> intersectRay(const Ray<GlobalPosition3D> &, FaceHandle start) 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;
@@ -86,6 +146,8 @@ public:
 
 	[[nodiscard]] HalfedgeHandle findEntry(const GlobalPosition2D from, const GlobalPosition2D to) const;
 
+	void setHeights(const std::span<const GlobalPosition3D> triangleStrip);
+
 	[[nodiscard]] auto
 	getExtents() const
 	{
@@ -102,9 +164,26 @@ protected:
 
 	[[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;
+
+	template<glm::length_t D>
+	[[nodiscard]] static RelativePosition<D>
+	difference(const GlobalPosition<D> a, const GlobalPosition<D> b)
+	{
+		return b - a;
+	}
+
+	[[nodiscard]] RelativeDistance length(const HalfedgeHandle) const;
+	[[nodiscard]] GlobalPosition3D centre(const HalfedgeHandle) const;
 
 	void update_vertex_normals_only();
+	void update_vertex_normals_only(VertexIter start);
+
+	using OpenMesh::TriMesh_ArrayKernelT<GeoDataTraits>::split;
+	std::array<FaceHandle, 4> split(FaceHandle);
 
 private:
 	GlobalPosition3D lowerExtent {}, upperExtent {};
diff --git a/game/network/link.cpp b/game/network/link.cpp
index 122eaf4..248fe7d 100644
--- a/game/network/link.cpp
+++ b/game/network/link.cpp
@@ -62,7 +62,7 @@ LinkCurve::intersectRay(const Ray<GlobalPosition3D> & ray) const
 	const auto & e1p {ends[1].node->pos};
 	const auto slength = round_frac(length / 2.F, 5.F);
 	const auto segs = std::round(15.F * slength / std::pow(radius, 0.7F));
-	const auto step {glm::vec<2, RelativeDistance> {arc_length(arc), e1p.z - e0p.z} / segs};
+	const auto step {glm::vec<2, RelativeDistance> {arc.length(), e1p.z - e0p.z} / segs};
 
 	auto segCount = static_cast<std::size_t>(std::lround(segs)) + 1;
 	std::vector<GlobalPosition3D> points;
diff --git a/game/network/rail.cpp b/game/network/rail.cpp
index 34cbceb..fd07ace 100644
--- a/game/network/rail.cpp
+++ b/game/network/rail.cpp
@@ -95,9 +95,8 @@ round_sleepers(const float v)
 	return round_frac(v, sleepers);
 }
 
-RailLinkStraight::RailLinkStraight(
-		NetworkLinkHolder<RailLinkStraight> & instances, const Node::Ptr & a, const Node::Ptr & b) :
-	RailLinkStraight(instances, a, b, b->pos - a->pos)
+RailLinkStraight::RailLinkStraight(NetworkLinkHolder<RailLinkStraight> & instances, const Node::Ptr & a,
+		const Node::Ptr & b) : RailLinkStraight(instances, a, b, b->pos - a->pos)
 {
 }
 
@@ -109,16 +108,15 @@ RailLinkStraight::RailLinkStraight(
 {
 }
 
-RailLinkCurve::RailLinkCurve(
-		NetworkLinkHolder<RailLinkCurve> & instances, const Node::Ptr & a, const Node::Ptr & b, GlobalPosition2D c) :
-	RailLinkCurve(instances, a, b, c || a->pos.z, {c || 0, 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, {c, a->pos, b->pos})
 {
 }
 
 RailLinkCurve::RailLinkCurve(NetworkLinkHolder<RailLinkCurve> & instances, const Node::Ptr & a, const Node::Ptr & b,
 		GlobalPosition3D c, const Arc arc) :
 	Link({a, normalize(arc.first + half_pi)}, {b, normalize(arc.second - half_pi)},
-			glm::length(RelativePosition3D(a->pos - c)) * arc_length(arc)),
+			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)}
diff --git a/game/terrain.cpp b/game/terrain.cpp
index d2c8593..91a228f 100644
--- a/game/terrain.cpp
+++ b/game/terrain.cpp
@@ -30,14 +30,14 @@ Terrain::generateMeshes()
 	std::vector<Vertex> vertices;
 	vertices.reserve(geoData->n_vertices());
 	std::map<GeoData::VertexHandle, size_t> vertexIndex;
-	std::transform(geoData->vertices_begin(), geoData->vertices_end(), std::back_inserter(vertices),
+	std::transform(geoData->vertices_sbegin(), geoData->vertices_end(), std::back_inserter(vertices),
 			[this, &vertexIndex](const GeoData::VertexHandle v) {
 				vertexIndex.emplace(v, vertexIndex.size());
 				const auto p = geoData->point(v);
-				return Vertex {p, p / 10000, geoData->normal(v)};
+				return Vertex {p, RelativePosition2D(p) / 10000.F, geoData->normal(v)};
 			});
 	std::for_each(
-			geoData->faces_begin(), geoData->faces_end(), [this, &vertexIndex, &indices](const GeoData::FaceHandle f) {
+			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](const GeoData::VertexHandle v) {
 							return vertexIndex[v];
diff --git a/lib/collections.h b/lib/collections.h
index 943b986..dd603be 100644
--- a/lib/collections.h
+++ b/lib/collections.h
@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <array>
+#include <cstdint>
 #include <span>
 #include <utility>
 #include <vector>
@@ -129,18 +130,25 @@ vectorOfN(std::integral auto N, T start = {}, T step = 1)
 	return v;
 }
 
+template<template<typename...> typename Rtn = std::vector, typename In>
+[[nodiscard]] auto
+materializeRange(const In begin, const In end)
+{
+	return Rtn(begin, end);
+}
+
 template<template<typename...> typename Rtn = std::vector, IterableCollection In>
 [[nodiscard]] auto
-materializeRange(In && in)
+materializeRange(const In & in)
 {
-	return Rtn(in.begin(), in.end());
+	return materializeRange<Rtn>(in.begin(), in.end());
 }
 
 template<template<typename...> typename Rtn = std::vector, typename In>
 [[nodiscard]] auto
 materializeRange(const std::pair<In, In> & in)
 {
-	return Rtn(in.first, in.second);
+	return materializeRange<Rtn>(in.first, in.second);
 }
 
 template<typename T> struct pair_range {
@@ -160,3 +168,62 @@ template<typename T> struct pair_range {
 };
 
 template<typename T> pair_range(std::pair<T, T>) -> pair_range<T>;
+
+template<typename iter> struct stripiter {
+	[[nodiscard]] constexpr bool
+	operator!=(const stripiter & other) const
+	{
+		return current != other.current;
+	}
+
+	[[nodiscard]] constexpr bool
+	operator==(const stripiter & other) const
+	{
+		return current == other.current;
+	}
+
+	constexpr stripiter &
+	operator++()
+	{
+		++current;
+		off = 1 - off;
+		return *this;
+	}
+
+	constexpr stripiter &
+	operator--()
+	{
+		--current;
+		off = 1 - off;
+		return *this;
+	}
+
+	constexpr auto
+	operator-(const stripiter & other) const
+	{
+		return current - other.current;
+	}
+
+	constexpr auto
+	operator*() const
+	{
+		return std::tie(*(current - (2 - off)), *(current - off - 1), *current);
+	}
+
+	iter current;
+	uint8_t off {};
+};
+
+template<typename T> struct std::iterator_traits<stripiter<T>> : std::iterator_traits<T> { };
+
+constexpr auto
+strip_begin(IterableCollection auto & cont)
+{
+	return stripiter {cont.begin() + 2};
+}
+
+constexpr auto
+strip_end(IterableCollection auto & cont)
+{
+	return stripiter {cont.end()};
+}
diff --git a/lib/geometricPlane.h b/lib/geometricPlane.h
index 3f95d3c..d4b803d 100644
--- a/lib/geometricPlane.h
+++ b/lib/geometricPlane.h
@@ -16,18 +16,22 @@ public:
 
 template<typename PositionType> class GeometricPlaneT : public GeometricPlane {
 public:
+	GeometricPlaneT() = default;
+
+	GeometricPlaneT(PositionType origin, Normal3D normal) : origin(std::move(origin)), normal(normal) { }
+
 	struct DistAndPosition {
 		PositionType::value_type dist;
 		PositionType position;
 	};
 
-	PositionType origin;
-	Normal3D normal;
+	PositionType origin {};
+	Normal3D normal {};
 
 	[[nodiscard]] inline PlaneRelation
 	getRelation(PositionType point) const
 	{
-		const auto d = glm::dot(normal, point - origin);
+		const auto d = glm::dot(normal, RelativePosition3D(point - origin));
 		return d < 0.F ? PlaneRelation::Below : d > 0.F ? PlaneRelation::Above : PlaneRelation::On;
 	}
 
diff --git a/lib/maths.cpp b/lib/maths.cpp
index 68662fc..51e27fe 100644
--- a/lib/maths.cpp
+++ b/lib/maths.cpp
@@ -4,6 +4,11 @@
 #include <glm/gtx/rotate_vector.hpp>
 #include <glm/gtx/transform.hpp>
 
+Arc::Arc(const RelativePosition2D & dir0, const RelativePosition2D & dir1) :
+	Arc {vector_yaw(dir0), vector_yaw(dir1)} { }
+
+Arc::Arc(const Angle anga, const Angle angb) : pair {anga, (angb < anga) ? angb + two_pi : angb} { }
+
 glm::mat4
 flat_orientation(const Direction3D & diff)
 {
@@ -76,7 +81,7 @@ rotate_yp(Rotation2D a)
 }
 
 float
-vector_yaw(const Direction3D & diff)
+vector_yaw(const Direction2D & diff)
 {
 	return std::atan2(diff.x, diff.y);
 }
diff --git a/lib/maths.h b/lib/maths.h
index 5886326..63b752a 100644
--- a/lib/maths.h
+++ b/lib/maths.h
@@ -8,17 +8,28 @@
 #include <stdexcept>
 #include <utility>
 
-struct Arc : public std::pair<float, float> {
-	using std::pair<float, float>::pair;
+struct Arc : public std::pair<Angle, Angle> {
+	template<glm::length_t Lc, glm::length_t Le, typename T, glm::qualifier Q>
+		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()}}
+	{
+	}
 
-	template<typename T, glm::qualifier Q>
-	Arc(const glm::vec<3, T, Q> & centre3, const glm::vec<3, T, Q> & e0p, const glm::vec<3, T, Q> & e1p);
+	Arc(const RelativePosition2D & dir0, const RelativePosition2D & dir1);
+	Arc(const Angle angb, const Angle anga);
 
 	auto
 	operator[](bool i) const
 	{
 		return i ? second : first;
 	}
+
+	[[nodiscard]] constexpr inline float
+	length() const
+	{
+		return second - first;
+	}
 };
 
 constexpr const RelativePosition3D up {0, 0, 1};
@@ -84,9 +95,16 @@ glm::mat4 rotate_pitch(float);
 glm::mat4 rotate_yp(Rotation2D);
 glm::mat4 rotate_ypr(Rotation3D);
 
-float vector_yaw(const Direction3D & diff);
+float vector_yaw(const Direction2D & diff);
 float vector_pitch(const Direction3D & diff);
 
+template<typename T, glm::qualifier Q>
+glm::vec<2, T, Q>
+vector_normal(const glm::vec<2, T, Q> & v)
+{
+	return {-v.y, v.x};
+};
+
 float round_frac(const float & v, const float & frac);
 
 template<typename T>
@@ -111,7 +129,7 @@ template<std::integral T, glm::qualifier Q>
 inline constexpr glm::vec<3, T, Q>
 crossProduct(const glm::vec<3, T, Q> a, const glm::vec<3, T, Q> b)
 {
-	return crossProduct<int64_t, Q>(a, b);
+	return crossProduct<Q>(a, b);
 }
 
 template<std::floating_point T, glm::qualifier Q>
@@ -171,12 +189,6 @@ operator%=(glm::vec<L, T, Q> & p, const glm::mat<L + 1, L + 1, T, Q> & mutation)
 	return p = p % mutation;
 }
 
-constexpr inline float
-arc_length(const Arc & arc)
-{
-	return arc.second - arc.first;
-}
-
 float normalize(float ang);
 
 template<typename T, glm::qualifier Q>
@@ -244,21 +256,6 @@ midpoint(const std::pair<T, T> & v)
 	return std::midpoint(v.first, v.second);
 }
 
-template<typename T, glm::qualifier Q>
-Arc::Arc(const glm::vec<3, T, Q> & centre3, const glm::vec<3, T, Q> & e0p, const glm::vec<3, T, Q> & e1p) :
-	Arc([&]() -> Arc {
-		const auto diffa = e0p - centre3;
-		const auto diffb = e1p - centre3;
-		const auto anga = vector_yaw(diffa);
-		const auto angb = [&diffb, &anga]() {
-			const auto angb = vector_yaw(diffb);
-			return (angb < anga) ? angb + two_pi : angb;
-		}();
-		return {anga, angb};
-	}())
-{
-}
-
 // Conversions
 template<typename T>
 inline constexpr auto
diff --git a/lib/persistence.h b/lib/persistence.h
index c53ff99..3c95a00 100644
--- a/lib/persistence.h
+++ b/lib/persistence.h
@@ -81,11 +81,11 @@ namespace Persistence {
 			return make_s<SelectionT<T>>(value);
 		}
 
-		template<typename S>
+		template<typename S, typename... Extra>
 		[[nodiscard]] static SelectionPtr
-		make_s(T & value)
+		make_s(T & value, Extra &&... extra)
 		{
-			return std::make_unique<S>(value);
+			return std::make_unique<S>(value, std::forward<Extra>(extra)...);
 		}
 
 		T & v;
@@ -331,6 +331,65 @@ namespace Persistence {
 		}
 	};
 
+	template<typename... T> struct SelectionT<std::tuple<T...>> : public SelectionV<std::tuple<T...>> {
+		using V = std::tuple<T...>;
+		using SelectionV<V>::SelectionV;
+
+		struct Members : public SelectionV<V> {
+			template<size_t... Idx>
+			explicit Members(V & v, std::integer_sequence<size_t, Idx...>) :
+				SelectionV<V> {v}, members {SelectionV<std::tuple_element_t<Idx, V>>::make(std::get<Idx>(v))...}
+			{
+			}
+
+			void
+			beforeValue(Stack & stk) override
+			{
+				stk.push(std::move(members[idx++]));
+			}
+
+			std::size_t idx {0};
+			std::array<SelectionPtr, std::tuple_size_v<V>> members;
+		};
+
+		void
+		beginArray(Stack & stk) override
+		{
+			stk.push(this->template make_s<Members>(
+					this->v, std::make_integer_sequence<size_t, std::tuple_size_v<V>>()));
+		}
+	};
+
+	template<typename T, typename U> struct SelectionT<std::pair<T, U>> : public SelectionV<std::pair<T, U>> {
+		using V = std::pair<T, U>;
+		using SelectionV<V>::SelectionV;
+
+		struct Members : public SelectionV<V> {
+			explicit Members(V & v) :
+				SelectionV<V> {v}, members {
+										   SelectionV<T>::make(v.first),
+										   SelectionV<U>::make(v.second),
+								   }
+			{
+			}
+
+			void
+			beforeValue(Stack & stk) override
+			{
+				stk.push(std::move(members[idx++]));
+			}
+
+			std::size_t idx {0};
+			std::array<SelectionPtr, 2> members;
+		};
+
+		void
+		beginArray(Stack & stk) override
+		{
+			stk.push(this->template make_s<Members>(this->v));
+		}
+	};
+
 	template<typename Map, typename Type = typename Map::mapped_type, auto Key = &Type::element_type::id>
 	struct MapByMember : public Persistence::SelectionT<Type> {
 		MapByMember(Map & m) : Persistence::SelectionT<Type> {s}, map {m} { }
diff --git a/lib/ray.h b/lib/ray.h
index e1f43c3..a831270 100644
--- a/lib/ray.h
+++ b/lib/ray.h
@@ -43,15 +43,27 @@ public:
 	}
 
 	bool
+	intersectPlane(const PositionType orig, const Direction3D norm, RelativeDistance & distance) const
+	{
+		if constexpr (std::is_floating_point_v<typename PositionType::value_type>) {
+			return glm::intersectRayPlane(start, direction, orig, norm, distance) && distance >= 0.F;
+		}
+		else {
+			const RelativePosition3D origr = orig - start;
+			return glm::intersectRayPlane({}, direction, origr, norm, distance) && distance >= 0.F;
+		}
+	}
+
+	bool
 	intersectTriangle(const PositionType t0, const PositionType t1, const PositionType t2, BaryPosition & bary,
 			RelativeDistance & distance) const
 	{
 		if constexpr (std::is_floating_point_v<typename PositionType::value_type>) {
-			return glm::intersectRayTriangle(start, direction, t0, t1, t2, bary, distance);
+			return glm::intersectRayTriangle(start, direction, t0, t1, t2, bary, distance) && distance >= 0.F;
 		}
 		else {
 			const RelativePosition3D t0r = t0 - start, t1r = t1 - start, t2r = t2 - start;
-			return glm::intersectRayTriangle({}, direction, t0r, t1r, t2r, bary, distance);
+			return glm::intersectRayTriangle({}, direction, t0r, t1r, t2r, bary, distance) && distance >= 0.F;
 		}
 	}
 
diff --git a/lib/stream_support.h b/lib/stream_support.h
index fa536f1..57d82a1 100644
--- a/lib/stream_support.h
+++ b/lib/stream_support.h
@@ -11,17 +11,16 @@
 template<typename S>
 concept stringlike = requires(const S & s) { s.substr(0); };
 template<typename T>
-concept spanable = std::is_constructible_v<std::span<const typename T::value_type>, T> && !stringlike<T>
-		&& !std::is_same_v<std::span<typename T::value_type>, T>;
+concept NonStringIterableCollection
+		= std::is_same_v<decltype(std::declval<T>().begin()), decltype(std::declval<T>().end())> && !stringlike<T>;
 
 namespace std {
-	template<typename T, std::size_t L>
 	std::ostream &
-	operator<<(std::ostream & s, const span<T, L> v)
+	operator<<(std::ostream & s, const NonStringIterableCollection auto & v)
 	{
 		s << '(';
 		for (const auto & i : v) {
-			if (&i != &v.front()) {
+			if (&i != &*v.begin()) {
 				s << ", ";
 			}
 			s << i;
@@ -43,13 +42,6 @@ namespace std {
 		return (s << std::span {&v[0], L});
 	}
 
-	template<spanable T>
-	std::ostream &
-	operator<<(std::ostream & s, const T & v)
-	{
-		return (s << std::span {v});
-	}
-
 	template<typename First, typename Second>
 	std::ostream &
 	operator<<(std::ostream & s, const std::pair<First, Second> & v)
diff --git a/test/Jamfile.jam b/test/Jamfile.jam
index 733ef05..cce4513 100644
--- a/test/Jamfile.jam
+++ b/test/Jamfile.jam
@@ -12,6 +12,7 @@ project : requirements
 	<define>BOOST_TEST_DYN_LINK
 	<define>RESDIR=\\\"$(res)/\\\"
 	<define>FIXTURESDIR=\\\"$(fixtures)/\\\"
+	<define>GLM_FORCE_SWIZZLE
 	<variant>debug:<warnings>pedantic
 	<variant>debug:<warnings-as-errors>on
 	<variant>debug:<cflags>-Wnon-virtual-dtor
@@ -46,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 : -- : fixtures/height/SD19.asc : <library>test ;
+run 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/1.json b/test/fixtures/geoData/deform/1.json
new file mode 100644
index 0000000..6930238
--- /dev/null
+++ b/test/fixtures/geoData/deform/1.json
@@ -0,0 +1,201 @@
+[
+	[
+		[
+			[
+				70100,
+				123000,
+				6000
+			],
+			[
+				50100,
+				52300,
+				6000
+			],
+			[
+				191000,
+				283000,
+				8000
+			],
+			[
+				241000,
+				123330,
+				2000
+			]
+		],
+		[
+			[
+				[
+					[
+						20000,
+						20000,
+						90000
+					],
+					[
+						1,
+						1,
+						-1.5
+					]
+				],
+				"/tmp/geoData0.tga"
+			],
+			[
+				[
+					[
+						30000,
+						164000,
+						90000
+					],
+					[
+						1,
+						-1,
+						-1.5
+					]
+				],
+				"/tmp/geoData1.tga"
+			],
+			[
+				[
+					[
+						288000,
+						162000,
+						90000
+					],
+					[
+						-1,
+						-1,
+						-1.5
+					]
+				],
+				"/tmp/geoData2.tga"
+			]
+		]
+	],
+	[
+		[
+			[
+				3000,
+				1000,
+				500
+			],
+			[
+				3000,
+				2000,
+				500
+			],
+			[
+				2000,
+				1000,
+				500
+			]
+		],
+		[
+			[
+				[
+					[
+						-1000,
+						-3000,
+						7000
+					],
+					[
+						1,
+						1,
+						-1.5
+					]
+				],
+				"/tmp/geoData3.tga"
+			],
+			[
+				[
+					[
+						1800,
+						2500,
+						800
+					],
+					[
+						1,
+						-0.2,
+						-0.5
+					]
+				],
+				"/tmp/geoData4.tga"
+			]
+		]
+	],
+	[
+		[
+			[
+				3000,
+				1000,
+				10
+			],
+			[
+				3000,
+				2000,
+				10
+			],
+			[
+				2000,
+				1000,
+				10
+			]
+		],
+		[
+			[
+				[
+					[
+						-500,
+						-1500,
+						3000
+					],
+					[
+						1,
+						1,
+						-1.5
+					]
+				],
+				"/tmp/geoData5.tga"
+			]
+		]
+	],
+	[
+		[
+			[
+				1500,
+				2000,
+				800
+			],
+			[
+				3000,
+				2000,
+				800
+			],
+			[
+				5000,
+				4000,
+				800
+			],
+			[
+				3500,
+				700,
+				800
+			]
+		],
+		[
+			[
+				[
+					[
+						-1000,
+						-3000,
+						7000
+					],
+					[
+						1,
+						1,
+						-1.5
+					]
+				],
+				"/tmp/geoData6.tga"
+			]
+		]
+	]
+]
diff --git a/test/test-geoData.cpp b/test/test-geoData.cpp
index 302cab7..fb9aba0 100644
--- a/test/test-geoData.cpp
+++ b/test/test-geoData.cpp
@@ -1,7 +1,12 @@
 #define BOOST_TEST_MODULE terrain
+#include "game/terrain.h"
+#include "test/testMainWindow.h"
+#include "test/testRenderOutput.h"
 #include "testHelpers.h"
+#include "ui/applicationBase.h"
 #include <boost/test/data/test_case.hpp>
 #include <boost/test/unit_test.hpp>
+#include <gfx/gl/sceneRenderer.h>
 #include <stream_support.h>
 
 #include <game/geoData.h>
@@ -114,7 +119,7 @@ BOOST_DATA_TEST_CASE(findRayIntersect,
 		}),
 		p, d, i)
 {
-	BOOST_CHECK_EQUAL(fixedTerrtain.intersectRay({p, d}).value(), i);
+	BOOST_CHECK_EQUAL(fixedTerrtain.intersectRay({p, d})->first, i);
 }
 
 BOOST_AUTO_TEST_CASE(boundaryWalk)
@@ -187,6 +192,23 @@ BOOST_DATA_TEST_CASE(walkTerrainUntil,
 	BOOST_CHECK_EQUAL_COLLECTIONS(visited.begin(), visited.end(), visits.begin(), visits.end());
 }
 
+BOOST_AUTO_TEST_CASE(triangle_helpers)
+{
+	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_CHECK_CLOSE(t.area(), 12.5F, 0.01F);
+}
+
 using FindEntiesData = std::tuple<GlobalPosition2D, GlobalPosition2D, int>;
 
 BOOST_DATA_TEST_CASE(findEntries,
@@ -199,3 +221,56 @@ BOOST_DATA_TEST_CASE(findEntries,
 {
 	BOOST_CHECK_EQUAL(fixedTerrtain.findEntry(from, to).idx(), heh);
 }
+
+using DeformTerrainData = std::tuple<std::vector<GlobalPosition3D>,
+		std::vector<std::pair<std::pair<GlobalPosition3D, Direction3D>, std::string>>>;
+
+BOOST_TEST_DECORATOR(*boost::unit_test::timeout(2));
+
+BOOST_DATA_TEST_CASE(deform, loadFixtureJson<DeformTerrainData>("geoData/deform/1.json"), points, cams)
+{
+	auto gd = std::make_shared<GeoData>(GeoData::createFlat({0, 0}, {1000000, 1000000}, 100));
+	BOOST_CHECK_NO_THROW(gd->setHeights(points));
+
+	ApplicationBase ab;
+	TestMainWindow tmw;
+	TestRenderOutput tro {{640, 480}};
+
+	struct TestTerrain : public SceneProvider {
+		explicit TestTerrain(std::shared_ptr<GeoData> gd) : terrain(std::move(gd)) { }
+
+		const Terrain terrain;
+
+		void
+		content(const SceneShader & shader) const override
+		{
+			terrain.render(shader);
+		}
+
+		void
+		environment(const SceneShader &, const SceneRenderer & sr) const override
+		{
+			sr.setAmbientLight({0.1, 0.1, 0.1});
+			sr.setDirectionalLight({1, 1, 1}, south + down, *this);
+		}
+
+		void
+		lights(const SceneShader &) const override
+		{
+		}
+
+		void
+		shadows(const ShadowMapper & shadowMapper) const override
+		{
+			terrain.shadows(shadowMapper);
+		}
+	};
+
+	TestTerrain t {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));
+		BOOST_CHECK_NO_THROW(ss.render(t));
+		Texture::save(tro.outImage, cam.second.c_str());
+	});
+}
diff --git a/test/test-lib.cpp b/test/test-lib.cpp
index 58b769a..5f0b5e5 100644
--- a/test/test-lib.cpp
+++ b/test/test-lib.cpp
@@ -4,6 +4,7 @@
 #include <boost/test/data/test_case.hpp>
 #include <boost/test/unit_test.hpp>
 
+#include <collections.h>
 #include <glArrays.h>
 #include <glad/gl.h>
 #include <set>
@@ -49,3 +50,22 @@ BOOST_AUTO_TEST_CASE(generate_move_and_delete)
 	}
 	BOOST_CHECK(active.empty());
 }
+
+constexpr std::array TRIANGLE_STRIP_IN {0, 1, 2, 3, 4, 5};
+static_assert(std::distance(strip_begin(TRIANGLE_STRIP_IN), strip_end(TRIANGLE_STRIP_IN)) == 4);
+
+BOOST_AUTO_TEST_CASE(triangle_strip_iter)
+{
+	constexpr std::array TRIANGLE_STRIP_EXPECTED {0, 1, 2, 2, 1, 3, 2, 3, 4, 4, 3, 5};
+
+	std::vector<int> out;
+	std::for_each(strip_begin(TRIANGLE_STRIP_IN), strip_end(TRIANGLE_STRIP_IN), [&out](const auto & t) {
+		const auto [a, b, c] = t;
+		out.push_back(a);
+		out.push_back(b);
+		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());
+}
diff --git a/test/test-static-stream_support.cpp b/test/test-static-stream_support.cpp
index 3002ccc..6bf9ea4 100644
--- a/test/test-static-stream_support.cpp
+++ b/test/test-static-stream_support.cpp
@@ -1,11 +1,20 @@
 #include "stream_support.h"
 
 #include <array>
+#include <map>
+#include <set>
 #include <vector>
 
-static_assert(spanable<std::vector<int>>);
-static_assert(spanable<std::vector<char>>);
-static_assert(spanable<std::array<int, 1>>);
-static_assert(spanable<std::array<char, 1>>);
-static_assert(!spanable<std::string>);
-static_assert(!spanable<std::string_view>);
+static_assert(NonStringIterableCollection<std::vector<int>>);
+static_assert(NonStringIterableCollection<std::set<int>>);
+static_assert(NonStringIterableCollection<std::map<int, int>>);
+static_assert(NonStringIterableCollection<std::vector<char>>);
+static_assert(NonStringIterableCollection<std::array<int, 1>>);
+static_assert(NonStringIterableCollection<std::array<char, 1>>);
+static_assert(!NonStringIterableCollection<std::string>);
+static_assert(!NonStringIterableCollection<std::string_view>);
+
+static_assert(requires(std::vector<int> i, std::ostream & o) { o << i; });
+static_assert(requires(std::array<int, 10> i, std::ostream & o) { o << i; });
+static_assert(requires(std::set<int> i, std::ostream & o) { o << i; });
+static_assert(requires(std::map<int, int> i, std::ostream & o) { o << i; });
diff --git a/test/testHelpers.h b/test/testHelpers.h
index f2b0901..58e4372 100644
--- a/test/testHelpers.h
+++ b/test/testHelpers.h
@@ -2,11 +2,22 @@
 
 #include <boost/test/tools/context.hpp>
 #include <boost/test/tools/interface.hpp>
+#include <filesystem>
+#include <fstream>
 #include <iomanip> // IWYU pragma: keep std::setprecision
+#include <jsonParse-persistence.h>
 #include <memory>
 
 std::unique_ptr<char, decltype(&free)> uasprintf(const char * fmt, ...) __attribute__((format(printf, 1, 2)));
 
+template<typename T>
+decltype(auto)
+loadFixtureJson(const std::filesystem::path & path)
+{
+	std::ifstream in {FIXTURESDIR / path};
+	return Persistence::JsonParsePersistence {}.loadState<std::vector<T>>(in);
+}
+
 #define BOOST_CHECK_CLOSE_VEC(a_, b_) \
 	{ \
 		const auto a {a_}, b {b_}; \
diff --git a/ui/builders/freeExtend.cpp b/ui/builders/freeExtend.cpp
index 47356c3..db127e6 100644
--- a/ui/builders/freeExtend.cpp
+++ b/ui/builders/freeExtend.cpp
@@ -19,7 +19,7 @@ BuilderFreeExtend::move(
 			candidateLinks.objects = network->candidateJoins(*p1, p->pos);
 		}
 		else if (const auto p = geoData->intersectRay(ray)) {
-			candidateLinks.objects = network->candidateExtend(*p1, *p);
+			candidateLinks.objects = network->candidateExtend(*p1, p->first);
 		}
 		else {
 			candidateLinks.removeAll();
@@ -42,8 +42,8 @@ BuilderFreeExtend::click(
 					p1 = p->pos;
 				}
 				else if (const auto p = geoData->intersectRay(ray)) {
-					network->addExtend(*p1, *p);
-					p1 = *p;
+					network->addExtend(*p1, p->first);
+					p1 = p->first;
 				}
 			}
 			else {
diff --git a/ui/builders/straight.cpp b/ui/builders/straight.cpp
index 0c4a3e2..43f5ec8 100644
--- a/ui/builders/straight.cpp
+++ b/ui/builders/straight.cpp
@@ -16,7 +16,7 @@ BuilderStraight::move(
 {
 	if (p1) {
 		if (const auto p = geoData->intersectRay(ray)) {
-			candidateLinks.objects = network->candidateStraight(*p1, *p);
+			candidateLinks.objects = network->candidateStraight(*p1, p->first);
 		}
 		else {
 			candidateLinks.removeAll();
@@ -32,12 +32,12 @@ BuilderStraight::click(
 		case SDL_BUTTON_LEFT:
 			if (const auto p = geoData->intersectRay(ray)) {
 				if (p1) {
-					create(network, *p1, *p);
+					create(network, *p1, p->first);
 					candidateLinks.removeAll();
 					p1.reset();
 				}
 				else {
-					p1 = p;
+					p1 = p->first;
 				}
 			}
 			return;