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

1 files changed, 338 insertions, 327 deletions

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()
+{
 }