Fixed bug when snapping polygon with hole to land boundary (#394 | GR…

…IDEDIT-1558)
Deltares · Jan 9, 2025 · 905ee14 · 905ee14
1 parent 2299246
commit 905ee14
Show file tree

Hide file tree

Showing 2 changed files with 105 additions and 6 deletions.
diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp
@@ -3373,18 +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<meshkernel::Point>& snappedPolygonPoints = polygons.Enclosure(enclosureIndex).Outer().Nodes();
 
-            for (int i = 0; i < polygon.num_coordinates; ++i)
+            // Copy the snapped polygon points
+            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
@@ -38,3 +38,84 @@ 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<int>(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<int>(xPolygon.size());
+
+    errorCode = mkernel_polygon_snap_to_landboundary(meshKernelId, land, polygon, 4, 1);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    const std::vector<double> expectedXAfterFirstSnapping{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};
+
+    const std::vector<double> expectedYAfterFirstSnapping{-0.25, -0.25, 10.0, 10.0, -0.25,
+                                                          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};
+
+    ASSERT_EQ(polygon.num_coordinates, static_cast<int>(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<double> expectedXAfterSecondSnapping{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};
+
+    const std::vector<double> expectedYAfterSecondSnapping{-0.25, -0.25, 10.0, 10.0, -0.25,
+                                                           meshkernel::constants::missing::innerOuterSeparator,
+                                                           2.0, 2.0, 7.0, 7.0, 2.0,
+                                                           meshkernel::constants::missing::doubleValue,
+                                                           -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);
+    }
+}