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.

1 files changed, 343 insertions, 16 deletions

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});
 }