From 2ea20df4b47b3de4b6ca137134db317f78e2e83a Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 7 Jan 2025 17:05:10 +0100 Subject: [PATCH] GRIDEDIT-1558 Fixed bug when snapping polygon with hole to land boundary --- libs/MeshKernel/tests/src/PolygonsTests.cpp | 53 ++++++++++++++++++ libs/MeshKernelApi/src/MeshKernel.cpp | 29 +++++++--- libs/MeshKernelApi/tests/src/PolygonTests.cpp | 54 +++++++++++++++++++ 3 files changed, 130 insertions(+), 6 deletions(-) diff --git a/libs/MeshKernel/tests/src/PolygonsTests.cpp b/libs/MeshKernel/tests/src/PolygonsTests.cpp index c50e1bab1..1a057d3f9 100644 --- a/libs/MeshKernel/tests/src/PolygonsTests.cpp +++ b/libs/MeshKernel/tests/src/PolygonsTests.cpp @@ -1114,3 +1114,56 @@ TEST(Polygons, LinearRefinePolygonOutOfBoundsCheck) // polygon end node index incorrect EXPECT_THROW(auto refinedPolygon = polygons.LinearRefinePolygon(0, 1, 9), meshkernel::ConstraintError); } + +TEST(Polygons, SnapPolygonWithHoleToLandBoundary) +{ + // Test the algorithm for snapping multi-polygons to land boundaries. + + //constexpr double tolerance = 1.0e-8; + + // The land boundary to which the polygon is to be snapped. + // This land boundary is made up of two sections + std::vector landBoundaryPoints{{-1.0, -0.2}, + {11.0, -0.2}}; + + meshkernel::LandBoundary landBoundary(landBoundaryPoints); + + // The original polygon points. + // The points make up two distinct polygons + std::vector polygonPoints{{0.0, 0.0}, + {10.0, 0.0}, + {10.0, 10.0}, + {0.0, 10.0}, + {0.0, 0.0}, + {meshkernel::constants::missing::innerOuterSeparator, meshkernel::constants::missing::innerOuterSeparator}, + {2.0, 2.0}, + {7.0, 2.0}, + {7.0, 7.0}, + {2.0, 7.0}, + {2.0, 2.0}}; + + meshkernel::Polygons polygon(polygonPoints, meshkernel::Projection::cartesian); + + // Snap the polygon to the land boundary + // Only snaps the first polygon + polygon.SnapToLandBoundary(landBoundary, 0, 0); + + const meshkernel::Polygon& firstPolygon = polygon.Enclosure(0).Outer(); + const meshkernel::Polygon& innerPolygon = polygon.Enclosure(0).Inner(0); + + std::cout << "first poly "; + + for (size_t i = 0; i < firstPolygon.Size (); ++i) { + std::cout << "{" << firstPolygon.Node (i).x << " " << firstPolygon.Node (i).y << "} "; + } + + std::cout << std::endl; + + std::cout << "inner poly "; + + for (size_t i = 0; i < innerPolygon.Size (); ++i) { + std::cout << "{" << innerPolygon.Node (i).x << " " << innerPolygon.Node (i).y << "} "; + } + + std::cout << std::endl; +} diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index 00e4833aa..9ca4c9f8c 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -3373,19 +3373,36 @@ namespace meshkernelapi meshkernel::LandBoundary landBoundary(landBoundaryPoints); meshkernel::Polygons polygons(polygonPoints, meshKernelState[meshKernelId].m_mesh2d->m_projection); - polygons.SnapToLandBoundary(landBoundary, startIndex, endIndex); const auto [enclosureIndex, polygonStartNode, polygonEndNode] = polygons.PolygonIndex(startIndex, endIndex); + polygons.SnapToLandBoundary(landBoundary, startIndex, endIndex); - //-------------------------------- - // Now copy back the polygon values + meshkernel::UInt enclosureStartPoint = 0; + meshkernel::UInt enclosurePreviousStartPoint = 0; + meshkernel::UInt enclosureCount = 0; + + // Find the index of the first point in the polygnal enclosure with index enclosureIndex + for (meshkernel::UInt i = 0; i < polygonPoints.size (); ++i) + { + if ((polygonPoints [i].x == meshkernel::constants::missing::doubleValue && polygonPoints [i].y == meshkernel::constants::missing::doubleValue) || i == polygonPoints.size () - 1) + { + if (enclosureCount == enclosureIndex) { + enclosureStartPoint = enclosurePreviousStartPoint; + break; + } + + enclosurePreviousStartPoint = i + 1; + ++enclosureCount; + } + } const std::vector& snappedPolygonPoints = polygons.Enclosure(enclosureIndex).Outer().Nodes(); - for (int i = 0; i < polygon.num_coordinates; ++i) + for (size_t i = 0; i < snappedPolygonPoints.size (); ++i) { - polygon.coordinates_x[i] = snappedPolygonPoints[i].x; - polygon.coordinates_y[i] = snappedPolygonPoints[i].y; + polygon.coordinates_x[i + enclosureStartPoint] = snappedPolygonPoints[i].x; + polygon.coordinates_y[i + enclosureStartPoint] = snappedPolygonPoints[i].y; } + } catch (...) { diff --git a/libs/MeshKernelApi/tests/src/PolygonTests.cpp b/libs/MeshKernelApi/tests/src/PolygonTests.cpp index db6056eb9..8d81511e9 100644 --- a/libs/MeshKernelApi/tests/src/PolygonTests.cpp +++ b/libs/MeshKernelApi/tests/src/PolygonTests.cpp @@ -38,3 +38,57 @@ TEST(PolygonTests, PolygonRefinementTests) ASSERT_NEAR(10.0, yPolygon[3], tolerance); ASSERT_NEAR(10.0, yPolygon[4], tolerance); } + +TEST(PolygonTests, PolygonWithHoleSnappingTest) +{ + constexpr double tolerance = 1.0e-8; + + int meshKernelId; + auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::GeometryList land; + land.geometry_separator = meshkernel::constants::missing::doubleValue; + std::vector xLand{-1.0, 31.0}; + std::vector yLand{-0.25, -0.25}; + land.coordinates_x = xLand.data(); + land.coordinates_y = yLand.data(); + land.num_coordinates = static_cast(xLand.size()); + + meshkernelapi::GeometryList polygon; + polygon.geometry_separator = meshkernel::constants::missing::doubleValue; + polygon.inner_outer_separator = meshkernel::constants::missing::innerOuterSeparator; + + std::vector xPolygon{0.0, 10.0, 10.0, 0.0, 0.0, meshkernel::constants::missing::innerOuterSeparator, 2.0, 7.0, 7.0, 2.0, 2.0, meshkernel::constants::missing::doubleValue, 20.0, 30.0, 30.0, 20.0, 20.0}; + std::vector yPolygon{0.0, 0.0, 10.0, 10.0, 0.0, meshkernel::constants::missing::innerOuterSeparator, 2.0, 2.0, 7.0, 7.0, 2.0, meshkernel::constants::missing::doubleValue, 0.0, 0.0, 10.0, 10.0, 0.0}; + + polygon.coordinates_x = xPolygon.data(); + polygon.coordinates_y = yPolygon.data(); + polygon.num_coordinates = static_cast(xPolygon.size()); + + errorCode = mkernel_polygon_snap_to_landboundary(meshKernelId, land, polygon, 4, 1); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector expectedXAfterFirstSnapping{0.0, 10.0, 10.0, 0.0, 0.0, -998.0, 2.0, 7.0, 7.0, 2.0, 2.0, -999.0, 20.0, 30.0, 30.0, 20.0, 20.0}; + const std::vector expectedYAfterFirstSnapping{-0.25, -0.25, 10.0, 10.0, -0.25, -998.0, 2.0, 2.0, 7.0, 7.0, 2.0, -999.0, 0.0, 0.0, 10.0, 10.0, 0.0}; + + ASSERT_EQ(polygon.num_coordinates, static_cast(expectedXAfterFirstSnapping.size())); + + for (size_t i = 0; i < expectedXAfterFirstSnapping.size(); ++i) + { + EXPECT_NEAR(polygon.coordinates_x[i], expectedXAfterFirstSnapping[i], tolerance); + EXPECT_NEAR(polygon.coordinates_y[i], expectedYAfterFirstSnapping[i], tolerance); + } + + errorCode = mkernel_polygon_snap_to_landboundary(meshKernelId, land, polygon, 16, 13); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector expectedXAfterSecondSnapping{0.0, 10.0, 10.0, 0.0, 0.0, -998.0, 2.0, 7.0, 7.0, 2.0, 2.0, -999.0, 20.0, 30.0, 30.0, 20.0, 20.0}; + const std::vector expectedYAfterSecondSnapping{-0.25, -0.25, 10.0, 10.0, -0.25, -998.0, 2.0, 2.0, 7.0, 7.0, 2.0, -999.0, -0.25, -0.25, 10.0, 10.0, -0.25}; + + for (size_t i = 0; i < expectedXAfterSecondSnapping.size(); ++i) + { + EXPECT_NEAR(polygon.coordinates_x[i], expectedXAfterSecondSnapping[i], tolerance); + EXPECT_NEAR(polygon.coordinates_y[i], expectedYAfterSecondSnapping[i], tolerance); + } +}