Splines.hpp

//---- GPL ---------------------------------------------------------------------
//
// Copyright (C)  Stichting Deltares, 2011-2021.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation version 3.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
// contact: delft3d.support@deltares.nl
// Stichting Deltares
// P.O. Box 177
// 2600 MH Delft, The Netherlands
//
// All indications and logos of, and references to, "Delft3D" and "Deltares"
// are registered trademarks of Stichting Deltares, and remain the property of
// Stichting Deltares. All rights reserved.
//
//------------------------------------------------------------------------------
 
#pragma once
 
#include "MeshKernel/BoundingBox.hpp"
#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/LandBoundary.hpp"
#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Utilities/LinearAlgebra.hpp"
 
namespace meshkernel
{
    class CurvilinearGrid;
 
    class Splines
    {
 
    public:
        Splines() = default;
 
        explicit Splines(Projection projection);
 
        explicit Splines(CurvilinearGrid const& grid);
 
        // void AddSpline(const std::vector<Point>& splines, UInt start, UInt size);
        void AddSpline(const std::vector<Point>& splines, UInt start, UInt size);
 
        void AddSpline(const std::vector<Point>& splines);
 
        void Replace(const UInt splineIndex, const std::vector<Point>& splinePoints);
 
        void SwapSplines(const UInt firstSpline, const UInt secondSpline);
 
        void Reverse(const UInt splineId);
 
        void SnapSpline(const size_t splineIndex,
                        const LandBoundary& landBoundary,
                        const int numberOfIterations = constants::numeric::defaultSnappingIterations);
 
        bool GetSplinesIntersection(UInt first,
                                    UInt second,
                                    double& crossProductIntersection,
                                    Point& intersectionPoint,
                                    double& firstSplineRatio,
                                    double& secondSplineRatio) const;
 
        [[nodiscard]] double ComputeSplineLength(UInt index,
                                                 double startAdimensionalCoordinate,
                                                 double endAdimensionalCoordinate,
                                                 UInt numSamples = 100,
                                                 bool accountForCurvature = false,
                                                 double height = 1.0,
                                                 double assignedDelta = -1.0) const;
 
        std::tuple<std::vector<Point>, std::vector<double>> ComputePointOnSplineFromAdimensionalDistance(UInt index,
                                                                                                         double maximumGridHeight,
                                                                                                         bool isSpacingCurvatureAdapted,
                                                                                                         const std::vector<double>& distances) const;
 
        Point ComputeClosestPointOnSplineSegment(UInt index, double startSplineSegment, double endSplineSegment, Point point) const;
 
        Point ComputeClosestPoint(UInt index, Point point) const;
 
        Point Evaluate(const UInt whichSpline, const double lambda) const;
 
        BoundingBox GetBoundingBox(const UInt splineIndex) const;
 
        auto GetNumSplines() const { return static_cast<UInt>(m_splineNodes.size()); }
 
        UInt Size(const UInt whichSpline) const;
 
        UInt MaxSizeIndex() const;
 
        std::vector<std::vector<Point>> m_splineNodes;       
        std::vector<std::vector<Point>> m_splineDerivatives; 
        std::vector<double> m_splinesLength;                 
        Projection m_projection = Projection::cartesian;     
 
    private:
        void AddPointInExistingSpline(UInt splineIndex, const Point& point);
 
        std::tuple<Point, Point, double> ComputeCurvatureOnSplinePoint(UInt splineIndex,
                                                                       double adimensionalPointCoordinate) const;
 
        void DeleteSpline(UInt splineIndex);
 
        void AllocateSplinesProperties();
 
        static std::vector<Point> ComputeSplineDerivative(const std::vector<Point>& splinesNodes);
    };
 
    struct FuncAdimensionalToDimensionalDistanceOnSpline
    {
        FuncAdimensionalToDimensionalDistanceOnSpline(const Splines& splines,
                                                      UInt splineIndex,
                                                      bool isSpacingCurvatureAdapted,
                                                      double h) : m_spline(splines),
                                                                  m_splineIndex(splineIndex),
                                                                  m_isSpacingCurvatureAdapted(isSpacingCurvatureAdapted),
                                                                  m_h(h) {}
 
        FuncAdimensionalToDimensionalDistanceOnSpline(const Splines& splines,
                                                      UInt splineIndex,
                                                      bool isSpacingCurvatureAdapted,
                                                      double h,
                                                      double distance) : m_spline(splines),
                                                                         m_splineIndex(splineIndex),
                                                                         m_isSpacingCurvatureAdapted(isSpacingCurvatureAdapted),
                                                                         m_h(h),
                                                                         m_DimensionalDistance(distance) {}
 
        void SetDimensionalDistance(double distance)
        {
            m_DimensionalDistance = distance;
        }
 
        double operator()(double adimensionalDistanceReferencePoint) const
        {
            auto distanceFromReferencePoint = m_spline.ComputeSplineLength(m_splineIndex, 0.0, adimensionalDistanceReferencePoint, m_numSamples, m_isSpacingCurvatureAdapted, m_h, 0.1);
            distanceFromReferencePoint = std::abs(distanceFromReferencePoint - m_DimensionalDistance);
            return distanceFromReferencePoint;
        }
 
        const Splines& m_spline;            
        UInt m_splineIndex;                 
        bool m_isSpacingCurvatureAdapted;   
        double m_h;                         
        UInt m_numSamples = 10;             
        double m_DimensionalDistance = 0.0; 
    };
 
    struct FuncDistanceFromAPoint
    {
        FuncDistanceFromAPoint(const Splines& splines,
                               UInt splineIndex,
                               Point point) : m_spline(splines),
                                              m_splineIndex(splineIndex),
                                              m_point(point) {}
 
        double operator()(double adimensionalDistanceReferencePoint) const
        {
            const auto pointOnSpline = ComputePointOnSplineAtAdimensionalDistance(m_spline.m_splineNodes[m_splineIndex], m_spline.m_splineDerivatives[m_splineIndex], adimensionalDistanceReferencePoint);
            return ComputeDistance(m_point, pointOnSpline, Projection::cartesian);
        }
 
        const Splines& m_spline;            
        UInt m_splineIndex;                 
        Point m_point;                      
        double m_DimensionalDistance = 0.0; 
    };
 
} // namespace meshkernel