diff --git a/libs/MeshKernel/include/MeshKernel/BoundingBox.hpp b/libs/MeshKernel/include/MeshKernel/BoundingBox.hpp index d09c00863..184100502 100644 --- a/libs/MeshKernel/include/MeshKernel/BoundingBox.hpp +++ b/libs/MeshKernel/include/MeshKernel/BoundingBox.hpp @@ -106,6 +106,11 @@ namespace meshkernel point.y >= m_lowerLeft.y && point.y <= m_upperRight.y; } + /// @brief Checks if two bounding boxes overlaps + /// @param[in] boundingBox The input bounding box + /// @return True if the point if the this bounding box overlaps with another, false otherwise + bool Overlaps(const BoundingBox& boundingBox) const; + /// @brief Returns the lower left corner of the bounding box /// @return The lower left corner of the bounding box [[nodiscard]] auto& lowerLeft() const { return m_lowerLeft; } @@ -141,6 +146,18 @@ namespace meshkernel /// @brief Return the delta of the bounding box. Vector Delta() const; + /// @brief Create a bounding box from two points + template T> + static BoundingBox CreateBoundingBox(const T& first, const T& second) + { + const auto lowerLeftX = std::min(first.x, second.x); + const auto lowerLeftY = std::min(first.y, second.y); + const auto upperRightX = std::max(first.x, second.x); + const auto upperRightY = std::max(first.y, second.y); + + return BoundingBox({lowerLeftX, lowerLeftY}, {upperRightX, upperRightY}); + } + private: Point m_lowerLeft; ///< The lower left corner of the bounding box Point m_upperRight; ///< The upper right corner of the bounding box @@ -203,3 +220,19 @@ inline meshkernel::Vector meshkernel::BoundingBox::Delta() const { return Vector(m_upperRight.x - m_lowerLeft.x, m_upperRight.y - m_lowerLeft.y); } + +inline bool meshkernel::BoundingBox::Overlaps(const BoundingBox& other) const +{ + + const auto& otherLowerleft = other.lowerLeft(); + const auto& otherUpperRight = other.upperRight(); + if (m_upperRight.x < otherLowerleft.x || + otherUpperRight.x < m_lowerLeft.x || + m_upperRight.y < otherLowerleft.y || + otherUpperRight.y < m_lowerLeft.y) + { + return false; + } + + return true; +} diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index dcd3f6ed3..8080f9064 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -308,6 +308,16 @@ namespace meshkernel /// it may be required to call Administrate after merging static Mesh2D Merge(const Mesh2D& mesh1, const Mesh2D& mesh2); + /// @brief Get the mesh bounding box + /// + /// @return The mesh bounding box + [[nodiscard]] BoundingBox GetBoundingBox() const; + + /// @brief Get the bounding boxes of the mesh edges + /// + /// @return The mesh edges bounding boxes + [[nodiscard]] std::vector GetEdgesBoundingBoxes() const; + private: // orthogonalization static constexpr double m_minimumEdgeLength = 1e-4; ///< Minimum edge length diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp index 6e853b972..1de9de516 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2DIntersections.hpp @@ -107,11 +107,13 @@ namespace meshkernel /// @returns The intersection seed std::tuple GetIntersectionSeed(const Mesh2D& mesh, const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, const std::vector& vistedEdges) const; /// @brief Gets the next edge intersection /// @returns The intersection seed std::tuple GetNextEdgeIntersection(const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, UInt edgeIndex, UInt firstIndex, UInt secondIndex, @@ -120,6 +122,7 @@ namespace meshkernel /// @brief Gets the next edge intersection /// @returns The intersection seed void IntersectFaceEdges(const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, const std::vector& cumulativeLength, UInt currentCrossingEdge, UInt currentFaceIndex, @@ -164,6 +167,8 @@ namespace meshkernel std::vector m_facesIntersectionsCache; ///< A cache for saving the local face intersections of one inner or outer std::vector m_edgesIntersections; ///< A vector collecting all edge intersection results std::vector m_faceIntersections; ///< A vector collecting all face intersection results + BoundingBox m_meshBoundingBox; ///< The mesh bounding box + std::vector m_meshEdgesBoundingBoxes; ///< The mesh edges bounding boxes static constexpr UInt maxSearchSegments = 1000; ///< max number of steps in polyline intersection algorithm }; diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 931255baa..5c43cb11a 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2019,3 +2019,33 @@ meshkernel::Mesh2D Mesh2D::Merge(const Mesh2D& mesh1, const Mesh2D& mesh2) return mergedMesh; } + +meshkernel::BoundingBox Mesh2D::GetBoundingBox() const +{ + + Point lowerLeft(std::numeric_limits::max(), std::numeric_limits::max()); + Point upperRight = -lowerLeft; + + const auto numNodes = GetNumNodes(); + for (UInt e = 0; e < numNodes; ++e) + { + lowerLeft.x = std::min(lowerLeft.x, m_nodes[e].x); + lowerLeft.y = std::min(lowerLeft.y, m_nodes[e].y); + upperRight.x = std::max(upperRight.x, m_nodes[e].x); + upperRight.y = std::max(upperRight.y, m_nodes[e].y); + } + + return BoundingBox(lowerLeft, upperRight); +} + +std::vector Mesh2D::GetEdgesBoundingBoxes() const +{ + std::vector result; + result.reserve(GetNumEdges()); + for (const auto& e : m_edges) + { + result.emplace_back(BoundingBox::CreateBoundingBox(m_nodes[e.first], m_nodes[e.second])); + } + + return result; +} diff --git a/libs/MeshKernel/src/Mesh2DIntersections.cpp b/libs/MeshKernel/src/Mesh2DIntersections.cpp index 06fe15ea6..9d5682ca1 100644 --- a/libs/MeshKernel/src/Mesh2DIntersections.cpp +++ b/libs/MeshKernel/src/Mesh2DIntersections.cpp @@ -27,10 +27,14 @@ Mesh2DIntersections::Mesh2DIntersections(Mesh2D& mesh) : m_mesh(mesh) // Declare local caches m_edgesIntersectionsCache.resize(mesh.GetNumEdges()); m_facesIntersectionsCache.resize(mesh.GetNumFaces()); + + m_meshBoundingBox = m_mesh.GetBoundingBox(); + m_meshEdgesBoundingBoxes = m_mesh.GetEdgesBoundingBoxes(); } std::tuple Mesh2DIntersections::GetIntersectionSeed(const Mesh2D& mesh, const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, const std::vector& vistedEdges) const { UInt crossedSegmentIndex = constants::missing::uintValue; @@ -48,6 +52,16 @@ std::tuple Mesh2DIntersections::GetIntersectionSeed(const Mesh2D& me continue; } + if (!m_meshBoundingBox.Overlaps(polyLineBoundingBoxes[segmentIndex])) + { + continue; + } + + if (!m_meshEdgesBoundingBoxes[edgeIndex].Overlaps(polyLineBoundingBoxes[segmentIndex])) + { + continue; + } + const auto [isEdgeCrossed, intersectionPoint, crossProductValue, @@ -79,6 +93,7 @@ std::tuple Mesh2DIntersections::GetIntersectionSeed(const Mesh2D& me } std::tuple Mesh2DIntersections::GetNextEdgeIntersection(const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, UInt edgeIndex, UInt firstIndex, UInt secondIndex, @@ -122,6 +137,18 @@ std::tuple Mesh2DIntersections::GetNex firstIndex = firstIndex - 1; } + if (!m_meshBoundingBox.Overlaps(polyLineBoundingBoxes[firstIndex])) + { + numSteps++; + continue; + } + + if (!m_meshEdgesBoundingBoxes[edgeIndex].Overlaps(polyLineBoundingBoxes[firstIndex])) + { + numSteps++; + continue; + } + auto intersection = AreSegmentsCrossing(polyLine[firstIndex], polyLine[secondIndex], m_mesh.m_nodes[m_mesh.m_edges[edgeIndex].first], @@ -141,6 +168,7 @@ std::tuple Mesh2DIntersections::GetNex } void Mesh2DIntersections::IntersectFaceEdges(const std::vector& polyLine, + const std::vector& polyLineBoundingBoxes, const std::vector& cumulativeLength, UInt currentCrossingEdge, UInt currentFaceIndex, @@ -179,7 +207,7 @@ void Mesh2DIntersections::IntersectFaceEdges(const std::vector& polyLine, secondIndex, crossProductValue, adimensionalPolylineSegmentDistance, - adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Forward); + adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, polyLineBoundingBoxes, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Forward); } if (!intersectionFound) @@ -189,7 +217,7 @@ void Mesh2DIntersections::IntersectFaceEdges(const std::vector& polyLine, secondIndex, crossProductValue, adimensionalPolylineSegmentDistance, - adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Backward); + adimensionalEdgeDistance) = GetNextEdgeIntersection(polyLine, polyLineBoundingBoxes, edgeIndex, segmentIndex, segmentIndex + 1, Direction::Backward); } // none of the polyline intersect the current edge @@ -228,9 +256,12 @@ void Mesh2DIntersections::Compute(const std::vector& polyLine) const auto polyLineSize = static_cast(polyLine.size()); std::vector cumulativeLength(polyLine.size(), 0.0); + std::vector polyLineBoundingBoxes(polyLine.size() - 1); for (UInt i = 1; i < polyLineSize; ++i) { cumulativeLength[i] = cumulativeLength[i - 1] + ComputeDistance(polyLine[i], polyLine[i - 1], m_mesh.m_projection); + + polyLineBoundingBoxes[i - 1] = BoundingBox::CreateBoundingBox(polyLine[i - 1], polyLine[i]); } std::queue> crossingEdges; @@ -244,6 +275,7 @@ void Mesh2DIntersections::Compute(const std::vector& polyLine) const auto [crossedEdgeIndex, crossedSegmentIndex] = GetIntersectionSeed( m_mesh, polyLine, + polyLineBoundingBoxes, vistedEdges); // no valid seed found in the entire mesh, we are done @@ -267,6 +299,7 @@ void Mesh2DIntersections::Compute(const std::vector& polyLine) continue; } IntersectFaceEdges(polyLine, + polyLineBoundingBoxes, cumulativeLength, currentCrossingEdge, currentFaceIndex,