GRIDEDIT-247 First time for snapping spline to land boundary code

libs/MeshKernel/include/MeshKernel/Entities.hpp

@@ -124,35 +124,11 @@ namespace meshkernel
         /// @brief Inplace multiply point by scalar.
         Point& operator*=(const double p);
 
-        /// @brief Overloads addition with another Point
-        [[nodiscard]] Point operator+(Point const& rhs) const { return Point(x + rhs.x, y + rhs.y); }
-
-        /// @brief Overloads addition with a double
-        [[nodiscard]] Point operator+(double const& rhs) const { return Point(x + rhs, y + rhs); }
-
-        /// @brief Overloads subtraction with another Point
-        [[nodiscard]] Point operator-(Point const& rhs) const { return Point(x - rhs.x, y - rhs.y); }
-
-        /// @brief Overloads subtraction with a double
-        [[nodiscard]] Point operator-(double const& rhs) const { return Point(x - rhs, y - rhs); }
-
-        /// @brief Overloads multiplication with another Point
-        [[nodiscard]] Point operator*(Point const& rhs) const { return Point(x * rhs.x, y * rhs.y); }
-
-        /// @brief Overloads multiplication with a double
-        [[nodiscard]] Point operator*(double const& rhs) const { return Point(x * rhs, y * rhs); }
-
-        /// @brief Overloads multiplication with a double
+        /// @brief Overloads multiplication with a integer
         [[nodiscard]] Point operator*(int const& rhs) const { return Point(x * rhs, y * rhs); }
 
-        /// @brief Overloads division with another Point
-        [[nodiscard]] Point operator/(Point const& rhs) const { return Point(x / rhs.x, y / rhs.y); }
-
-        /// @brief Overloads division with a double
-        [[nodiscard]] Point operator/(double const& rhs) const { return Point(x / rhs, y / rhs); }
-
-        /// @brief Overloads equality with another Point
-        bool operator==(const Point& rhs) const { return IsEqual(x, rhs.x) && IsEqual(y, rhs.y); }
+        // /// @brief Overloads equality with another Point
+        // bool operator==(const Point& rhs) const { return IsEqual(x, rhs.x) && IsEqual(y, rhs.y); }
 
         /// @brief Transforms spherical coordinates to cartesian
         void TransformSphericalToCartesian(double referenceLatitude)
@@ -178,6 +154,48 @@ namespace meshkernel
             constants::geometric::earth_radius; ///< Factor used in the transformation from spherical to Cartesian coordinates
     };
 
+    /// @brief Add two points
+    Point operator+(const Point& p1, const Point& p2);
+
+    /// @brief Add points and scalar
+    Point operator+(const Point& p, const double x);
+
+    /// @brief Add points and scalar
+    Point operator+(const double x, const Point& p);
+
+    /// @brief Subtract two points
+    Point operator-(const Point& p1, const Point& p2);
+
+    /// @brief Subtract scalar from point
+    Point operator-(const Point& p, const double x);
+
+    /// @brief Subtract point from scalar
+    Point operator-(const double x, const Point& p);
+
+    /// @brief Multiply point by a point
+    ///
+    /// Computes a point-wise product.
+    Point operator*(const Point& p1, const Point& p2);
+
+    /// @brief Multiply point by a scalar
+    Point operator*(const double x, const Point& p);
+
+    /// @brief Multiply point by a scalar
+    Point operator*(const Point& p, const double x);
+
+    /// @brief Divide point by a point
+    ///
+    /// Computes a point-wise division.
+    /// \note No check is made to determine is divisors are zero.
+    Point operator/(const Point& p1, const Point& p2);
+
+    /// @brief Divide point by a scalar
+    /// \note No check is made to determine is divisor is zero.
+    Point operator/(const Point& p, const double x);
+
+    /// @brief Test points for equality
+    bool operator==(const Point& p1, const Point& p2);
+
     /// @brief Describes an edge with two indices
     using Edge = std::pair<UInt, UInt>;
 
@@ -368,3 +386,63 @@ inline meshkernel::Point& meshkernel::Point::operator*=(const double p)
     y *= p;
     return *this;
 }
+
+inline meshkernel::Point meshkernel::operator+(const Point& p1, const Point& p2)
+{
+    return Point(p1.x + p2.x, p1.y + p2.y);
+}
+
+inline meshkernel::Point meshkernel::operator+(const Point& p, const double x)
+{
+    return Point(p.x + x, p.y + x);
+}
+
+inline meshkernel::Point meshkernel::operator+(const double x, const Point& p)
+{
+    return Point(p.x + x, p.y + x);
+}
+
+inline meshkernel::Point meshkernel::operator-(const Point& p1, const Point& p2)
+{
+    return Point(p1.x - p2.x, p1.y - p2.y);
+}
+
+inline meshkernel::Point meshkernel::operator-(const Point& p, const double x)
+{
+    return Point(p.x - x, p.y - x);
+}
+
+inline meshkernel::Point meshkernel::operator-(const double x, const Point& p)
+{
+    return Point(x - p.x, x - p.y);
+}
+
+inline meshkernel::Point meshkernel::operator*(const Point& p1, const Point& p2)
+{
+    return Point(p1.x * p2.x, p1.y * p2.y);
+}
+
+inline meshkernel::Point meshkernel::operator*(const double x, const Point& p)
+{
+    return Point(x * p.x, x * p.y);
+}
+
+inline meshkernel::Point meshkernel::operator*(const Point& p, const double x)
+{
+    return Point(x * p.x, x * p.y);
+}
+
+inline meshkernel::Point meshkernel::operator/(const Point& p1, const Point& p2)
+{
+    return Point(p1.x / p2.x, p1.y / p2.y);
+}
+
+inline meshkernel::Point meshkernel::operator/(const Point& p, const double x)
+{
+    return Point(p.x / x, p.y / x);
+}
+
+inline bool meshkernel::operator==(const Point& p1, const Point& p2)
+{
+    return IsEqual(p1.x, p2.x) && IsEqual(p1.y, p2.y);
+}

libs/MeshKernel/include/MeshKernel/LandBoundaries.hpp

@@ -86,6 +86,14 @@ namespace meshkernel
         /// @return The number of land boundary nodes.
         auto GetNumNodes() const { return m_nodes.size(); }
 
+        // toland
+        void nearestPointOnLandBoundary(const Point& pnt,
+                                        const int jainview,
+                                        Point& nearestPoint,
+                                        double& shortestDistance,
+                                        UInt& segmentStartIndex,
+                                        double& scaledDistanceToStart) const;
+
         std::vector<UInt> m_meshNodesLandBoundarySegments; ///< Mesh nodes to land boundary mapping (lanseg_map)
 
     private:

libs/MeshKernel/include/MeshKernel/Splines.hpp

@@ -27,7 +27,9 @@
 
 #pragma once
 
+#include "MeshKernel/Utilities/LinearAlgebra.hpp"
 #include <MeshKernel/Entities.hpp>
+#include <MeshKernel/LandBoundaries.hpp>
 #include <MeshKernel/Operations.hpp>
 
 namespace meshkernel
@@ -48,11 +50,11 @@ namespace meshkernel
         Splines() = default;
 
         /// @brief Ctor, set projection
-        /// @brief[in] projection The map projection
+        /// @brief[in] projection   The map projection
         explicit Splines(Projection projection);
 
         /// @brief Ctor from grids, each gridline is converted to spline, first the first m_n horizontal lines then the m_m vertical lines
-        /// @brief[in] The curvilinear grid
+        /// @brief[in] grid         The curvilinear grid
         explicit Splines(CurvilinearGrid const& grid);
 
         /// @brief Adds a new spline to m_splineCornerPoints
@@ -75,6 +77,28 @@ namespace meshkernel
         /// @returns coordinatesDerivatives The second order derivative at corner points (x derivative or y derivative)
         [[nodiscard]] static std::vector<double> SecondOrderDerivative(const std::vector<double>& coordinates, UInt startIndex, UInt endIndex);
 
+        static lin_alg::ColumnVector<double> ComputeSplineWeights(const std::vector<Point>& splinePoints,
+                                                                  const UInt totalNumberOfPoints,
+                                                                  const Projection projection);
+
+        /// @brief Evaluate a spline function (splint)
+        static Point evaluate(const std::vector<Point>& coordinates, const std::vector<Point>& secondDerivative, const double evaluationPoint);
+
+        // sample_spline
+        static void sampleSpline(const std::vector<Point>& controlPoints,
+                                 const size_t intermediatePointCount,
+                                 std::vector<Point>& samplePoints);
+
+        // comp_afinespline
+        static void compAfinespline(const UInt n, const UInt numRef, UInt& count, lin_alg::MatrixColMajor<double>& mat);
+
+        static lin_alg::MatrixColMajor<double> ComputeInterpolationMatrix(const lin_alg::MatrixColMajor<double>& splineCoefficients,
+                                                                          const lin_alg::ColumnVector<double>& weights);
+
+        // snap_spline
+        void snapSpline(const LandBoundaries& landBoundary,
+                        const size_t splineIndex);
+
         /// @brief Computes the intersection of two splines (sect3r)
         /// @param[in] first The index of the first spline
         /// @param[in] second The index of the second spline

libs/MeshKernel/include/MeshKernel/Utilities/LinearAlgebra.hpp

@@ -13,4 +13,10 @@ namespace lin_alg
     /// @tparam T Data type
     template <class T>
     using MatrixColMajor = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
+
+    /// @brief A dynamically sized column vector
+    /// @tparam T Data type
+    template <class T>
+    using ColumnVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;
+
 } // namespace lin_alg

libs/MeshKernel/src/LandBoundaries.cpp

@@ -1030,3 +1030,50 @@ void LandBoundaries::SnapMeshToLandBoundaries()
         }
     }
 }
+
+void LandBoundaries::nearestPointOnLandBoundary(const Point& samplePoint,
+                                                const int jainview,
+                                                Point& nearestPoint,
+                                                double& minimumDistance,
+                                                UInt& segmentStartIndex,
+                                                double& scaledDistanceToStart) const
+{
+
+    nearestPoint = samplePoint;
+    segmentStartIndex = constants::missing::intValue;
+    scaledDistanceToStart = -1.0;
+
+    minimumDistance = 9.0e33;
+
+    for (size_t i = 0; i < m_nodes.size() - 1; ++i)
+    {
+        Point firstPoint = m_nodes[i];
+        Point nextPoint = m_nodes[i + 1];
+
+        bool doIt = true;
+
+        if (jainview == 2)
+        {
+            doIt = IsPointInPolygonNodes(firstPoint, m_nodes, m_mesh->m_projection) && IsPointInPolygonNodes(nextPoint, m_nodes, m_mesh->m_projection);
+        }
+
+        if (doIt)
+        {
+            auto [distance, linePoint, distanceFromFirstNode] = DistanceFromLine(samplePoint, firstPoint, nextPoint, m_mesh->m_projection);
+            // TODO What is this and from where do I get it?
+            int ja = 1;
+
+            if (ja == 1)
+            {
+
+                if (distance < minimumDistance)
+                {
+                    nearestPoint = linePoint;
+                    minimumDistance = distance;
+                    segmentStartIndex = i;
+                    scaledDistanceToStart = distanceFromFirstNode;
+                }
+            }
+        }
+    }
+}