diff options
Diffstat (limited to 'game/geoData.cpp')
| -rw-r--r-- | game/geoData.cpp | 359 |
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}); } |