Operations.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 <iostream>
#include <numeric>
#include <span>
#include <vector>
 
#include "MeshKernel/BoundingBox.hpp"
#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/Utilities/LinearAlgebra.hpp"
#include "MeshKernel/Utilities/RTreeBase.hpp"
#include "MeshKernel/Vector.hpp"
 
namespace meshkernel
{
    template <typename T>
    void ResizeAndFill2DVector(std::vector<std::vector<T>>& v, UInt const& firstDimension, UInt const& secondDimension, bool fill = false, const T& fillValue = {})
    {
        v.resize(firstDimension);
        for (auto& e : v)
        {
            e.resize(secondDimension);
            if (fill)
            {
                std::fill(e.begin(), e.end(), fillValue);
            }
        }
    }
    void ResizeAndFill2DVector(std::vector<std::vector<T>>& v, UInt const& firstDimension, UInt const& secondDimension, bool fill = false, const T& fillValue = {}) {…}
 
    template <typename T>
    void ResizeAndFill3DVector(std::vector<std::vector<std::vector<T>>>& v, UInt const& firstDimension, UInt const& secondDimension, UInt const& thirdDim, bool fill = false, const T& fillValue = {})
    {
        v.resize(firstDimension);
        for (auto& e : v)
        {
            e.resize(secondDimension);
            for (auto& ee : e)
            {
                ee.resize(thirdDim);
                if (fill)
                {
                    std::fill(ee.begin(), ee.end(), fillValue);
                }
            }
        }
    }
    void ResizeAndFill3DVector(std::vector<std::vector<std::vector<T>>>& v, UInt const& firstDimension, UInt const& secondDimension, UInt const& thirdDim, bool fill = false, const T& fillValue = {}) {…}
 
    template <typename T>
    [[nodiscard]] UInt FindIndex(const std::vector<T>& vec, T el)
    {
        for (UInt n = 0; n < vec.size(); n++)
        {
            if (vec[n] == el)
            {
                return n;
            }
        }
 
        return constants::missing::uintValue;
    }
    [[nodiscard]] UInt FindIndex(const std::vector<T>& vec, T el) {…}
 
    template <typename T>
    UInt FindNextIndex(const std::vector<T>& vec, UInt current)
    {
        if (vec.size() == 0) [[unlikely]]
        {
            return constants::missing::uintValue;
        }
 
        return current == vec.size() - 1 ? 0 : current + 1;
    }
    UInt FindNextIndex(const std::vector<T>& vec, UInt current) {…}
 
    template <typename T>
    UInt FindPreviousIndex(const std::vector<T>& vec, UInt current)
    {
        if (vec.size() == 0) [[unlikely]]
        {
            return constants::missing::uintValue;
        }
 
        return current == 0 ? static_cast<UInt>(vec.size()) - 1 : current - 1;
    }
    UInt FindPreviousIndex(const std::vector<T>& vec, UInt current) {…}
 
    std::vector<std::pair<UInt, UInt>> FindIndices(const std::vector<Point>& vec, size_t start, size_t end, double separator);
 
    UInt InvalidPointCount(const std::vector<Point>& points);
 
    UInt InvalidPointCount(const std::vector<Point>& points, size_t start, size_t end);
 
    template <typename T>
    [[nodiscard]] std::vector<UInt> SortedIndices(const std::vector<T>& v)
    {
        std::vector<UInt> indices(v.size());
        iota(indices.begin(), indices.end(), 0);
        std::ranges::stable_sort(indices.begin(), indices.end(), [&v](UInt i1, UInt i2)
                                 { return v[i1] < v[i2]; });
        return indices;
    }
    [[nodiscard]] std::vector<UInt> SortedIndices(const std::vector<T>& v) {…}
 
    template <typename T>
    auto ReorderVector(const std::vector<T>& v, const std::vector<UInt>& order)
    {
        std::vector<T> ordered;
        ordered.reserve(v.size());
        for (const auto& value : order)
        {
            ordered.emplace_back(v[value]);
        }
        return ordered;
    }
    auto ReorderVector(const std::vector<T>& v, const std::vector<UInt>& order) {…}
 
    template <typename F>
    [[nodiscard]] double FindFunctionRootWithGoldenSectionSearch(F func, double min, double max)
    {
        // golden distance factors
        const double c = 0.38196602;
        const double r = 0.61803399;
        const double tolerance = 1e-5;
 
        double left = min;
        double middle = (min + max) * 0.5;
        double right = max;
 
        double x0 = left;
        double x1 = middle - c * (middle - left);
        double x2 = middle;
        double x3 = right;
        if (std::abs(right - middle) > std::abs(middle - left))
        {
            x1 = middle;
            x2 = middle + c * (right - left);
        }
 
        double f1 = func(x1);
        double f2 = func(x2);
 
        while (std::abs(x3 - x0) > tolerance * std::max(std::abs(x1) + std::abs(x2), 1e-8)) // tolerance * std::max(std::abs(x1) + std::abs(x2), 1e-10)
        {
            if (f2 < f1)
            {
                x0 = x1;
                x1 = x2;
                x2 = r * x1 + c * x3;
 
                f1 = f2;
                f2 = func(x2);
            }
            else
            {
                x3 = x2;
                x2 = x1;
                x1 = r * x2 + c * x0;
 
                f2 = f1;
                f1 = func(x1);
            }
        }
 
        return f1 < f2 ? x1 : x2;
    }
    [[nodiscard]] double FindFunctionRootWithGoldenSectionSearch(F func, double min, double max) {…}
 
    [[nodiscard]] UInt NextCircularForwardIndex(UInt currentIndex, UInt size);
 
    [[nodiscard]] UInt NextCircularBackwardIndex(UInt currentIndex, UInt size);
 
    [[nodiscard]] bool IsPointOnPole(const Point& point);
 
    [[nodiscard]] double crossProduct(const Point& firstSegmentFirstPoint, const Point& firstSegmentSecondPoint, const Point& secondSegmentFirstPoint, const Point& secondSegmentSecondPoint, const Projection& projection);
 
    [[nodiscard]] bool IsPointInPolygonNodes(const Point& point,
                                             const std::vector<Point>& polygonNodes,
                                             const Projection& projection,
                                             Point polygonCenter = {constants::missing::doubleValue, constants::missing::doubleValue},
                                             UInt startNode = constants::missing::uintValue,
                                             UInt endNode = constants::missing::uintValue);
 
    [[nodiscard]] bool IsPointInPolygonNodes(const Point& point,
                                             const std::vector<Point>& polygonNodes,
                                             const Projection& projection,
                                             const BoundingBox& boundingBox,
                                             Point polygonCenter = {constants::missing::doubleValue, constants::missing::doubleValue},
                                             UInt startNode = constants::missing::uintValue,
                                             UInt endNode = constants::missing::uintValue);
 
    [[nodiscard]] bool IsPointInTriangle(const Point& point,
                                         const std::span<const Point> triangleNodes,
                                         const Projection& projection);
 
    void ComputeThreeBaseComponents(const Point& point, std::array<double, 3>& exxp, std::array<double, 3>& eyyp, std::array<double, 3>& ezzp);
 
    void ComputeTwoBaseComponents(const Point& point, double (&elambda)[3], double (&ephi)[3]);
 
    [[nodiscard]] double GetDx(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    Vector GetDelta(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    Vector ComputeNormalToline(const Point& start, const Point& end, const Projection projection);
 
    [[nodiscard]] double GetDy(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    [[nodiscard]] double OuterProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment,
                                                 const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection);
 
    [[nodiscard]] Point ComputeMiddlePoint(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    [[nodiscard]] Point ComputeMiddlePointAccountingForPoles(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    [[nodiscard]] Point NormalVector(const Point& firstPoint, const Point& secondPoint, const Point& insidePoint, const Projection& projection);
 
    void TransformGlobalVectorToLocal(const Point& reference, const Point& globalCoordinates, const Point& globalComponents, const Projection& projection, Point& localComponents);
 
    Point NormalVectorOutside(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    void NormalVectorInside(const Point& firstPoint, const Point& secondPoint, const Point& insidePoint, Point& normal, bool& flippedNormal, const Projection& projection);
 
    void AddIncrementToPoint(const Point& normal, double increment, const Point& referencePoint, const Projection& projection, Point& point);
 
    void TranslateSphericalCoordinates(std::vector<Point>& polygon);
 
    [[nodiscard]] Point ReferencePoint(const std::vector<Point>& polygon, const Projection& projection);
 
    [[nodiscard]] Point ReferencePoint(const std::vector<Point>& nodes,
                                       const std::vector<UInt>& polygonIndices,
                                       const Projection& projection);
 
    [[nodiscard]] double ComputeSquaredDistance(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    [[nodiscard]] double ComputeDistance(const Point& firstPoint, const Point& secondPoint, const Projection& projection);
 
    [[maybe_unused]] std::tuple<double, Point, double> DistanceFromLine(const Point& point, const Point& firstNode, const Point& secondNode, const Projection& projection);
 
    [[nodiscard]] double InnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection);
 
    [[nodiscard]] double NormalizedInnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection);
 
    [[nodiscard]] Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection& projection);
 
    UInt CountNumberOfValidEdges(const std::vector<UInt>& edgesNumFaces, UInt numNodes);
 
    void ComputeMidPointsAndNormals(const std::vector<Point>& polygon,
                                    const std::vector<UInt>& edgesNumFaces,
                                    const UInt numNodes,
                                    std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
                                    std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals,
                                    UInt& pointCount,
                                    const Projection& projection);
 
    Point ComputeCircumCenter(const Point& centerOfMass,
                              const UInt pointCount,
                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& middlePoints,
                              const std::array<Point, constants::geometric::maximumNumberOfNodesPerFace>& normals,
                              const Projection& projection);
 
    Point ComputeFaceCircumenter(std::vector<Point>& polygon,
                                 const std::vector<UInt>& edgesNumFaces,
                                 const Projection& projection);
 
    [[nodiscard]] std::tuple<bool, Point, double, double, double> AreSegmentsCrossing(const Point& firstSegmentFirstPoint,
                                                                                      const Point& firstSegmentSecondPoint,
                                                                                      const Point& secondSegmentFirstPoint,
                                                                                      const Point& secondSegmentSecondPoint,
                                                                                      bool adimensionalCrossProduct,
                                                                                      const Projection& projection);
 
    template <typename T>
    [[nodiscard]] T ComputePointOnSplineAtAdimensionalDistance(const std::vector<T>& coordinates,
                                                               const std::vector<T>& coordinatesDerivatives,
                                                               double pointAdimensionalCoordinate)
    {
        T pointCoordinate{};
        if (pointAdimensionalCoordinate < 0)
        {
            return pointCoordinate;
        }
 
        const double eps = 1e-5;
        const double splFac = 1.0;
        const auto coordinate = std::floor(pointAdimensionalCoordinate);
 
        if (pointAdimensionalCoordinate - coordinate < eps)
        {
            return pointCoordinate = coordinates[static_cast<UInt>(coordinate)];
        }
 
        const UInt low = static_cast<UInt>(coordinate);
        const UInt high = low + 1;
        const double a = high - pointAdimensionalCoordinate;
        const double b = pointAdimensionalCoordinate - low;
 
        pointCoordinate = coordinates[low] * a + coordinates[high] * b +
                          (coordinatesDerivatives[low] * (pow(a, 3) - a) + coordinatesDerivatives[high] * (pow(b, 3) - b)) / 6.0 * splFac;
 
        return pointCoordinate;
    }
    [[nodiscard]] T ComputePointOnSplineAtAdimensionalDistance(const std::vector<T>& coordinates, {…}
 
    std::tuple<std::vector<double>, double> ComputeAdimensionalDistancesFromPointSerie(const std::vector<Point>& v, const Projection& projection);
 
    template <typename T>
    [[nodiscard]] int sgn(T val)
    {
        return (T(0) < val ? 1 : 0) - (val < T(0) ? 1 : 0);
    }
    [[nodiscard]] int sgn(T val) {…}
 
    [[nodiscard]] lin_alg::Matrix<Point> DiscretizeTransfinite(const std::vector<Point>& leftDiscretization,
                                                               const std::vector<Point>& rightDiscretization,
                                                               const std::vector<Point>& bottomDiscretization,
                                                               const std::vector<Point>& upperDiscretization,
                                                               const Projection& projection,
                                                               UInt numM,
                                                               UInt numN);
 
    [[nodiscard]] std::vector<Point> ComputeEdgeCenters(const std::vector<Point>& nodes, const std::vector<Edge>& edges);
 
    [[nodiscard]] double LinearInterpolationInTriangle(const Point& interpolationPoint, const std::vector<Point>& polygon, const std::vector<double>& values, const Projection& projection);
 
    [[nodiscard]] Point ComputeAverageCoordinate(const std::span<const Point> points, const Projection& projection);
 
    std::tuple<Point, double, bool> OrthogonalProjectionOnSegment(Point const& firstNode,
                                                                  Point const& secondNode,
                                                                  Point const& pointToProject);
 
    UInt AbsoluteDifference(UInt number_1, UInt number_2);
 
    [[nodiscard]] std::vector<Point> ComputePolyLineDiscretization(std::vector<Point> const& polyline, std::vector<double>& chainages, Projection projection);
 
    [[nodiscard]] std::vector<double> ComputePolyLineNodalChainages(std::vector<Point> const& polyline, Projection projection);
 
    [[nodiscard]] std::vector<double> ComputePolyLineEdgesLengths(std::vector<Point> const& polyline, Projection projection);
 
    [[nodiscard]] double MatrixNorm(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& matCoefficients);
 
    inline void IncrementValidValue(UInt& value, const UInt increment)
    {
        if (value != constants::missing::uintValue) [[likely]]
        {
            value += increment;
        }
    }
    inline void IncrementValidValue(UInt& value, const UInt increment) {…}
 
} // namespace meshkernel