diff options
Diffstat (limited to 'game/geoData.cpp')
| -rw-r--r-- | game/geoData.cpp | 821 |
1 files changed, 374 insertions, 447 deletions
diff --git a/game/geoData.cpp b/game/geoData.cpp index 72aa056..b886efd 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) @@ -36,16 +35,16 @@ GeoData::loadFromAsciiGrid(const std::filesystem::path & input) std::vector<VertexHandle> vertices; vertices.reserve(ncols * nrows); GeoData mesh; - mesh.lowerExtent = {xllcorner, yllcorner, std::numeric_limits<GlobalDistance>::max()}; - mesh.upperExtent = {xllcorner + (cellsize * (ncols - 1)), yllcorner + (cellsize * (nrows - 1)), - std::numeric_limits<GlobalDistance>::min()}; + mesh.extents = {{xllcorner, yllcorner, std::numeric_limits<GlobalDistance>::max()}, + {xllcorner + (cellsize * (ncols - 1)), yllcorner + (cellsize * (nrows - 1)), + std::numeric_limits<GlobalDistance>::min()}}; for (size_t row = 0; row < nrows; ++row) { for (size_t col = 0; col < ncols; ++col) { float heightf = 0; f >> heightf; const auto height = static_cast<GlobalDistance>(std::round(heightf * 1000.F)); - mesh.upperExtent.z = std::max(mesh.upperExtent.z, height); - mesh.lowerExtent.z = std::min(mesh.lowerExtent.z, height); + mesh.extents.max.z = std::max(mesh.extents.max.z, height); + mesh.extents.min.z = std::min(mesh.extents.min.z, height); vertices.push_back(mesh.add_vertex({xllcorner + (col * cellsize), yllcorner + (row * cellsize), height})); } } @@ -66,7 +65,7 @@ GeoData::loadFromAsciiGrid(const std::filesystem::path & input) }); } } - mesh.update_vertex_normals_only(); + mesh.updateAllVertexNormals(); return mesh; }; @@ -79,8 +78,7 @@ GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistan assert((upper - lower) % GRID_SIZE == GlobalPosition2D {}); GeoData mesh; - mesh.lowerExtent = {lower, h}; - mesh.upperExtent = {upper, h}; + mesh.extents = {{lower, h}, {upper, h}}; std::vector<VertexHandle> vertices; for (GlobalDistance row = lower.x; row <= upper.x; row += GRID_SIZE) { @@ -105,123 +103,11 @@ GeoData::createFlat(GlobalPosition2D lower, GlobalPosition2D upper, GlobalDistan } } - mesh.update_vertex_normals_only(); + mesh.updateAllVertexNormals(); 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 { @@ -233,13 +119,11 @@ 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) { - BaryPosition bari {}; - RelativeDistance dist {}; - const auto t = triangle<3>(face); - if (ray.intersectTriangle(t.x, t.y, t.z, bari, dist)) { - out.emplace(t * bari, face); + ray.start.xy() + (ray.direction.xy() * ::difference(extents.max.xy(), extents.min.xy())), + [&out, &ray, this](const auto & step) { + const auto t = triangle<3>(step.current); + if (const auto inter = ray.intersectTriangle(t.x, t.y, t.z)) { + out.emplace(t * inter->bary, step.current); return true; } return false; @@ -248,7 +132,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 +141,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 +232,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,339 +266,337 @@ 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 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() +GeoData::updateAllVertexNormals() { - update_vertex_normals_only(vertices_sbegin()); + updateAllVertexNormals(vertices()); } +template<std::ranges::range R> void -GeoData::update_vertex_normals_only(VertexIter start) +GeoData::updateAllVertexNormals(const R & range) { - 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)); - } + std::ranges::for_each(range, [this](const auto vertex) { + updateVertexNormal(vertex); }); } -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) +void +GeoData::updateVertexNormal(VertexHandle vertex) { - return triangleContainsPoint(a.x, b) && triangleContainsPoint(a.y, b) && triangleContainsPoint(a.z, b); + Normal3D normal {}; + { // Lifted from calc_vertex_normal_correct in PolyMeshT_impl.hh but doesn't use Scalar to normalise + ConstVertexIHalfedgeIter cvih_it = this->cvih_iter(vertex); + if (!cvih_it.is_valid()) { // don't crash on isolated vertices + return; + } + Normal in_he_vec; + calc_edge_vector(*cvih_it, in_he_vec); + for (; cvih_it.is_valid(); ++cvih_it) { // calculates the sector normal defined by cvih_it and adds it to normal + if (this->is_boundary(*cvih_it)) { + continue; + } + HalfedgeHandle out_heh(this->next_halfedge_handle(*cvih_it)); + Normal out_he_vec; + calc_edge_vector(out_heh, out_he_vec); + normal += cross(in_he_vec, out_he_vec); // sector area is taken into account + in_he_vec = out_he_vec; + in_he_vec *= -1; // change the orientation + } + } // End lift + set_normal(vertex, glm::normalize(normal)); } -void -GeoData::split(FaceHandle _fh) +OpenMesh::VertexHandle +GeoData::setPoint(GlobalPosition3D tsPoint, const RelativeDistance nearNodeTolerance) { - // 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); + 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; } - - if (split2) { - split(eh2, v2); + // 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); +}; - // Retriangulate - 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, const Surface & newFaceSurface) +std::set<GeoData::FaceHandle> +GeoData::setHeights(const std::span<const GlobalPosition3D> triangleStrip, const SetHeightsOpts & opts) { - static const RelativeDistance MAX_SLOPE = 1.5F; - static const RelativeDistance MIN_ARC = 0.01F; - if (triangleStrip.size() < 3) { - return; + return {}; + } + for (const auto & vertex : triangleStrip) { + extents += vertex; } - const auto initialVertexCount = static_cast<unsigned int>(n_vertices()); + 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); + }))} + { + } - // 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 - const auto initialFaceCount = static_cast<int>(n_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 fhi = FaceIter {*this, FaceHandle {initialFaceCount}, true}; fhi != faces_end(); fhi++) { - static constexpr auto MAX_FACE_AREA = 100'000'000.F; - const auto fh = *fhi; - if (triangle<3>(fh).area() > MAX_FACE_AREA) { - split(fh); + 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; } - } - std::vector<FaceHandle> newFaces; - std::copy_if(FaceIter {*this, FaceHandle {initialFaceCount}, true}, faces_end(), std::back_inserter(newFaces), - [this](FaceHandle fh) { - return !this->status(fh).deleted(); - }); - // Extrude corners - struct Extrusion { - VertexHandle boundaryVertex, extrusionVertex; - Direction3D lowerLimit, upperLimit; - }; + 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; + } - 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; + 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; + } + 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++; + }); + doBoundaryPart(*++newVerts.begin(), newVerts.front(), strip.front(), nearNodeTolerance); + doBoundaryPart(*++newVerts.rbegin(), newVerts.back(), strip.back(), nearNodeTolerance); + } + + using HeightSetTodo = std::multimap<VertexHandle, VertexHandle>; - 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)); + 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; } - 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))); + auto & point = geoData->point(vertex); + auto triangle = getTriangle(point); + if (!triangle) { + if (const auto boundaryVertex = boundaryTriangles.find(vertex); + boundaryVertex != boundaryTriangles.end()) { + triangle = boundaryVertex->second; } } - else { - // Single tangent bisecting the difference - addExtrusionFor(normalize(e0 + e1) / sinf((arc.length() - pi) / 2.F)); + 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); + } } - }, - *voh_begin(newVerts.front())); - - // Cut existing terrain - extrusionExtents.emplace_back(extrusionExtents.front()); // Circular next - std::vector<std::vector<VertexHandle>> boundaryFaces; - for (const auto & [first, second] : extrusionExtents | std::views::adjacent<2>) { - 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); + 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); + } + } + } } - out.emplace_back(second.extrusionVertex); - if (first.boundaryVertex != second.boundaryVertex) { - out.emplace_back(second.boundaryVertex); + + 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); + } } - } - // 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)); + 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; + } - delete_face(face, false); + 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; } - }; - removeOld(removeOld, findPoint(triangleStrip.front())); - std::for_each(boundaryFaces.begin(), boundaryFaces.end(), [&](auto & boundaryFace) { - std::reverse(boundaryFace.begin(), boundaryFace.end()); - add_face(boundaryFace); - }); + 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; + }; - // Tidy up - update_vertex_normals_only(VertexIter {*this, vertex_handle(initialVertexCount), true}); + return SetHeights {this, triangleStrip}.run(opts); +} - std::for_each(newFaces.begin(), newFaces.end(), [&newFaceSurface, this](const auto fh) { - property(surface, fh) = &newFaceSurface; - }); +void +GeoData::afterChange() +{ } |