GeometryUtilities.hpp

Go to the documentation of this file.

// _LICENSE_HEADER_
//
// Copyright (C) 2019 - 2025.
// Terms register on the GPL-3.0 license.
//
// This file can be redistributed and/or modified under the license terms.
//
// See top level LICENSE file for more details.
//
// This file can be used citing references in CITATION.cff file.
 
#ifndef __GEOMETRYUTILITIES_H
#define __GEOMETRYUTILITIES_H
 
#include "Eigen/Eigen"
#include "IOUtilities.hpp"
#include <iostream>
#include <numbers>
 
namespace Gedim
{
struct GeometryUtilitiesConfig final
{
    static constexpr double DefaultMinTolerance()
    {
        return 10.0 * std::numeric_limits<double>::epsilon();
    }
 
    double MinTolerance = Gedim::GeometryUtilitiesConfig::DefaultMinTolerance();
    double Tolerance1D = Gedim::GeometryUtilitiesConfig::DefaultMinTolerance();
    double Tolerance2D = Gedim::GeometryUtilitiesConfig::DefaultMinTolerance();
    double Tolerance3D = Gedim::GeometryUtilitiesConfig::DefaultMinTolerance();
};
 
class GeometryUtilities final
{
  private:
    const GeometryUtilitiesConfig &_configuration;
 
  public:
    enum struct CompareTypes
    {
        Unknown = 0,
        FirstBeforeSecond = 1,
        Coincident = 2,
        SecondBeforeFirst = 3
    };
 
    enum struct PolygonTypes
    {
        Unknown = 0,
        Triangle = 1,
        Quadrilateral_Convex = 2,
        Quadrilateral_Concave = 3,
        Generic_Convex = 4,
        Generic_Concave = 5
    };
 
    enum struct PolygonOrientations
    {
        Unknown = 0,
        Clockwise = 1,
        CounterClockwise = 2
    };
 
    enum struct PointSegmentPositionTypes
    {
        Unknown = 0,
        OnSegmentLineBeforeOrigin = 1,
        OnSegmentOrigin = 2,
        InsideSegment = 3,
        OnSegmentEnd = 4,
        OnSegmentLineAfterEnd = 5,
        LeftTheSegment = 6,
        RightTheSegment = 7
    };
 
    enum struct PointPlanePositionTypes
    {
        Unknown = 0,
        Negative = 1,
        OnPlane = 2,
        Positive = 3
    };
 
    enum struct PolygonCirclePositionTypes
    {
        Unknown = 0,
        PolygonOutsideCircleNoIntersection = 1,
        PolygonOutsideCircleOneIntersectionOnVertex = 2,
        PolygonOutsideCircleOneIntersectionTangentOnEdge = 3,
        CircleInsidePolygonNoIntersection = 4,
        CircleInsidePolygonOneIntersectionTangentOnEdge = 5,
        PolygonInsideCircleNoIntersection = 6,
        PolygonInsideCircleOneVertexIntersection = 7,
        PolygonInsideCircleIntersectionOnlyOnVertices = 8,
        CirclePolygonMultipleIntersections = 9
    };
 
    struct IntersectionPolygonCircleResult final
    {
        struct Intersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Secant = 1,
                Tangent = 2
            };
 
            enum struct IndexTypes
            {
                Unknown = 0,
                Vertex = 1,
                Edge = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolygonCircleResult::Intersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolygonCircleResult::Intersection::Types::Unknown;
            Gedim::GeometryUtilities::IntersectionPolygonCircleResult::Intersection::IndexTypes IndexType =
                Gedim::GeometryUtilities::IntersectionPolygonCircleResult::Intersection::IndexTypes::Unknown;
            unsigned int Index;
            double CurvilinearCoordinate; 
        };
 
        std::vector<Gedim::GeometryUtilities::IntersectionPolygonCircleResult::Intersection> Intersections = {}; 
    };
 
    struct PolygonDivisionByAngleQuadrantResult final
    {
        enum struct Types
        {
            Unknown = 0,
            ExternalOrigin = 1,
            Internal = 2,
            ExternalEnd = 3
        };
 
        Eigen::MatrixXd Points;                             
        std::vector<std::vector<unsigned int>> SubPolygons; 
        std::vector<Gedim::GeometryUtilities::PolygonDivisionByAngleQuadrantResult::Types> SubPolygonTypes; 
    };
 
    struct PolygonDivisionByCircleResult final
    {
        Eigen::MatrixXd Points;                              
        std::vector<std::vector<unsigned int>> SubTriangles; 
        std::vector<std::vector<unsigned int>> InternalTriangles; 
        std::vector<std::vector<unsigned int>> SubPolygons;       
    };
 
    struct CircleDivisionByPolygonResult final
    {
        Eigen::MatrixXd Points;                              
        std::vector<std::vector<unsigned int>> SubTriangles; 
        std::vector<std::vector<unsigned int>> InternalTriangles; 
        std::vector<std::vector<unsigned int>> SubPolygons;       
    };
 
    struct SplitPolygonInput final
    {
        struct AlignedEdge
        {
            unsigned int OriginVertexIndex;
            unsigned int EndVertexIndex;
        };
 
        struct SplitSegment
        {
            struct Vertex final
            {
                enum struct Types
                {
                    Unknown = 0,
                    Vertex = 1,
                    Edge = 2
                };
 
                Gedim::GeometryUtilities::SplitPolygonInput::SplitSegment::Vertex::Types Type;
                unsigned int Index;
            };
 
            Gedim::GeometryUtilities::SplitPolygonInput::SplitSegment::Vertex Origin;
            Gedim::GeometryUtilities::SplitPolygonInput::SplitSegment::Vertex End;
        };
 
        unsigned int NumberPolygonVertices;
        std::vector<Gedim::GeometryUtilities::SplitPolygonInput::AlignedEdge> AlignedEdges;
        Gedim::GeometryUtilities::SplitPolygonInput::SplitSegment Segment;
    };
 
    struct SplitPolygonWithSegmentResult final
    {
        enum struct Types
        {
            Unknown = 0,
            NoAction = 1,
            PolygonUpdate = 2,
            PolygonCreation = 3
        };
 
        struct NewVertex final
        {
            enum struct Types
            {
                Unknown = 0,
                SegmentOrigin = 1,
                SegmentEnd = 2
            };
 
            Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewVertex::Types Type =
                Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewVertex::Types::Unknown;
        };
 
        struct NewEdge final
        {
            enum struct Types
            {
                Unknown = 0,
                EdgeNew = 1,
                EdgeUpdate = 2
            };
 
            Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewEdge::Types Type =
                Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewEdge::Types::Unknown;
            unsigned int OldEdgeId = 0;
            unsigned int OriginId = 0;
            unsigned int EndId = 0;
            std::vector<unsigned int> Cell2DNeighbours = {};
        };
 
        struct NewPolygon final
        {
            std::list<unsigned int> Vertices = {};
            std::list<unsigned int> Edges = {};
        };
 
        Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::Types Type =
            Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::Types::Unknown;
        std::list<Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewVertex> NewVertices = {};
        std::list<Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewEdge> NewEdges = {};
        std::vector<Gedim::GeometryUtilities::SplitPolygonWithSegmentResult::NewPolygon> NewPolygons = {};
    };
 
    struct SplitPolygonWithCircleResult final
    {
        enum struct Types
        {
            Unknown = 0,
            NoAction = 1,
            PolygonUpdate = 2,
            PolygonCreation = 3
        };
 
        struct NewVertex final
        {
            enum struct Types
            {
                Unknown = 0,
                PolygonVertex = 1,
                CircleIntersection = 2,
                Both = 3
            };
 
            Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewVertex::Types Type =
                Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewVertex::Types::Unknown;
            unsigned int PolygonIndex;      
            unsigned int IntersectionIndex; 
        };
 
        struct NewEdge final
        {
            enum struct Types
            {
                Unknown = 0,
                Segment = 1,
                Arc = 2
            };
 
            enum struct ArcTypes
            {
                Unknown = 0,
                InsidePolygon = 1,
                OutsidePolygon = 2
            };
 
            Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewEdge::Types Type =
                Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewEdge::Types::Unknown;
            Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewEdge::ArcTypes ArcType =
                Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewEdge::ArcTypes::Unknown; 
            std::vector<unsigned int> VertexIndices; 
            unsigned int PolygonIndex;               
        };
 
        struct NewPolygon final
        {
            enum struct Types
            {
                Unknown = 0,
                InsideOnlyCircle = 1,
                InsideOnlyPolygon = 2,
                InsideCircleAndPolygon = 3
            };
 
            Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewPolygon::Types Type =
                Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewPolygon::Types::Unknown;
            std::vector<unsigned int> Vertices = {};
            std::vector<unsigned int> Edges = {};
        };
 
        Gedim::GeometryUtilities::SplitPolygonWithCircleResult::Types Type =
            Gedim::GeometryUtilities::SplitPolygonWithCircleResult::Types::Unknown;
        std::vector<unsigned int> PolygonVerticesNewVerticesPosition = {};
        std::vector<unsigned int> CircleIntersectionsNewVerticesPosition = {};
        std::vector<Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewVertex> NewVertices = {};
        std::vector<Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewEdge> NewEdges = {};
        std::vector<Gedim::GeometryUtilities::SplitPolygonWithCircleResult::NewPolygon> NewPolygons = {};
    };
 
    struct IntersectionSegmentSegmentResult final
    {
        enum struct IntersectionLineTypes
        {
            Unknown = 0,
            OnDifferentPlanes = 1,
            CoPlanarParallel = 2,
            CoPlanarIntersecting = 3
        };
 
        enum struct IntersectionSegmentTypes
        {
            Unknown = 0,
            NoIntersection = 1,
            SingleIntersection = 2,
            MultipleIntersections = 3
        };
 
        struct IntersectionPosition
        {
            Gedim::GeometryUtilities::PointSegmentPositionTypes Type = Gedim::GeometryUtilities::PointSegmentPositionTypes::Unknown;
            double CurvilinearCoordinate = 0.0;
        };
 
        Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionLineTypes IntersectionLinesType =
            Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionLineTypes::Unknown;
        Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionSegmentTypes IntersectionSegmentsType =
            Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionSegmentTypes::Unknown;
        std::vector<unsigned int> SecondIntersectionRelation;
        std::vector<Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionPosition> FirstSegmentIntersections;
        std::vector<Gedim::GeometryUtilities::IntersectionSegmentSegmentResult::IntersectionPosition> SecondSegmentIntersections; 
    };
 
    struct IntersectionSegmentCircleResult final
    {
        enum struct Types
        {
            Unknown = 0,
            NoIntersection = 1,
            TangentIntersection = 2,
            TwoIntersections = 3
        };
 
        struct IntersectionPosition
        {
            Gedim::GeometryUtilities::PointSegmentPositionTypes Type = Gedim::GeometryUtilities::PointSegmentPositionTypes::Unknown;
            double CurvilinearCoordinate = 0.0;
        };
 
        Gedim::GeometryUtilities::IntersectionSegmentCircleResult::Types Type =
            Gedim::GeometryUtilities::IntersectionSegmentCircleResult::Types::Unknown;
        std::vector<Gedim::GeometryUtilities::IntersectionSegmentCircleResult::IntersectionPosition> SegmentIntersections = {};
    };
 
    struct IntersectionSegmentPlaneResult final
    {
        enum struct Types
        {
            Unknown = 0,
            SingleIntersection = 1,
            NoIntersection = 2,
            MultipleIntersections = 3
        };
 
        struct Intersection
        {
            Gedim::GeometryUtilities::PointSegmentPositionTypes Type = Gedim::GeometryUtilities::PointSegmentPositionTypes::Unknown;
            double CurvilinearCoordinate = 0.0;
        };
 
        Gedim::GeometryUtilities::IntersectionSegmentPlaneResult::Types Type =
            Gedim::GeometryUtilities::IntersectionSegmentPlaneResult::Types::Unknown; 
        Gedim::GeometryUtilities::IntersectionSegmentPlaneResult::Intersection SingleIntersection; 
    };
 
    struct IntersectionPolyhedronLineResult final
    {
        enum struct Types
        {
            Unknown = 0,
            None = 1,                 
            OneIntersection = 2,      
            TwoIntersections = 3,     
            MultipleIntersections = 4 
        };
 
        struct PolyhedronFaceIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Intersection = 1,
                NoIntersection = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronFaceIntersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronFaceIntersection::Types::Unknown;
            unsigned int LineIntersectionIndex = 0; 
        };
 
        struct PolyhedronEdgeIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Intersection = 1,
                NoIntersection = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronEdgeIntersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronEdgeIntersection::Types::Unknown;
            unsigned int LineIntersectionIndex = 0; 
        };
 
        struct PolyhedronVertexIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Intersection = 1,
                NoIntersection = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronVertexIntersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronVertexIntersection::Types::Unknown;
            unsigned int LineIntersectionIndex = 0; 
        };
 
        struct LineIntersection
        {
            enum struct Types
            {
                Unknown = 0,
                OnVertex = 1, 
                OnEdge = 2,   
                OnFace = 3,   
                Inside = 4,   
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::LineIntersection::Types PolyhedronType =
                Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::LineIntersection::Types::Unknown; 
            unsigned int PolyhedronIndex = 0; 
            double CurvilinearCoordinate = 0.0; 
        };
 
        Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::Types Type =
            Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::Types::Unknown;
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::LineIntersection> LineIntersections; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronVertexIntersection> PolyhedronVertexIntersections = {}; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronEdgeIntersection> PolyhedronEdgeIntersections = {}; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronLineResult::PolyhedronFaceIntersection> PolyhedronFaceIntersections = {}; 
    };
 
    struct IntersectionPolyhedronsSegmentResult final
    {
        struct IntersectionPoint final
        {
            std::vector<unsigned int> Cell3DIndices = {};
        };
 
        struct IntersectionSegment final
        {
            std::vector<double> Points = {};
            std::vector<unsigned int> Cell3DIndices = {};
        };
 
        std::map<double, Gedim::GeometryUtilities::IntersectionPolyhedronsSegmentResult::IntersectionPoint> Points;
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronsSegmentResult::IntersectionSegment> Segments;
    };
 
    struct IntersectionPolyhedronPlaneResult final
    {
        enum struct Types
        {
            Unknown = 0,
            None = 1,      
            OnVertex = 2,  
            OnEdge = 3,    
            OnFace = 4,    
            NewPolygon = 5 
        };
 
        struct FaceIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Intersection = 1,
                NoIntersection = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::FaceIntersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::FaceIntersection::Types::Unknown;
        };
 
        struct EdgeIntersection final
        {
            Gedim::GeometryUtilities::IntersectionSegmentPlaneResult Intersection; 
        };
 
        struct VertexIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Intersection = 1,
                NoIntersection = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::VertexIntersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::VertexIntersection::Types::Unknown;
        };
 
        struct Intersection final
        {
            enum struct Types
            {
                Unknown = 0,
                Vertex = 1,
                Edge = 2
            };
 
            Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::Intersection::Types Type =
                Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::Intersection::Types::Unknown;
            unsigned int EdgeId = 0;   
            unsigned int VertexId = 0; 
        };
 
        Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::Types Type =
            Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::Types::Unknown; 
        unsigned int IntersectionId = 0; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::VertexIntersection> VertexIntersections = {}; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::EdgeIntersection> EdgeIntersections = {}; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::FaceIntersection> FaceIntersections = {}; 
        std::vector<Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult::Intersection> Intersections = {}; 
        Eigen::MatrixXd IntersectionCoordinates; 
    };
 
    enum struct PointCirclePositionResult
    {
        Unknown = 0,
        Outside = 1,
        OnBorder = 2,
        Inside = 3
    };
 
    struct PointPolygonPositionResult final
    {
        enum struct Types
        {
            Unknown = 0,
            Outside = 1,
            BorderEdge = 2,
            BorderVertex = 3,
            Inside = 4
        };
 
        unsigned int BorderIndex = 0; 
        Gedim::GeometryUtilities::PointPolygonPositionResult::Types Type =
            Gedim::GeometryUtilities::PointPolygonPositionResult::Types::Unknown;
    };
 
    struct LinePolygonPositionResult final
    {
        enum struct Types
        {
            Unknown = 0,
            Outside = 1,
            Intersecting = 2
        };
 
        struct EdgeIntersection final
        {
            enum struct Types
            {
                Unknown = 0,
                OnEdgeOrigin = 1,
                InsideEdge = 2,
                OnEdgeEnd = 3,
                Parallel = 4
            };
 
            Gedim::GeometryUtilities::LinePolygonPositionResult::EdgeIntersection::Types Type =
                Gedim::GeometryUtilities::LinePolygonPositionResult::EdgeIntersection::Types::Unknown;
            unsigned int Index = 0;
            double CurvilinearCoordinate = 0.0;
        };
 
        std::vector<Gedim::GeometryUtilities::LinePolygonPositionResult::EdgeIntersection> EdgeIntersections = {};
        Gedim::GeometryUtilities::LinePolygonPositionResult::Types Type =
            Gedim::GeometryUtilities::LinePolygonPositionResult::Types::Unknown;
    };
 
    struct PointPolyhedronPositionResult final
    {
        enum struct Types
        {
            Unknown = 0,
            Outside = 1,
            BorderFace = 2,
            BorderEdge = 3,
            BorderVertex = 4,
            Inside = 5
        };
 
        unsigned int BorderIndex = 0;               
        std::vector<unsigned int> Internal_indices; 
        Gedim::GeometryUtilities::PointPolyhedronPositionResult::Types Type =
            Gedim::GeometryUtilities::PointPolyhedronPositionResult::Types::Unknown;
    };
 
    struct SegmentPolyhedronPositionResult final
    {
        enum struct Types
        {
            Unknown = 0,
            Outside = 1,
            BorderFace = 2,
            BorderEdge = 3,
            Inside = 4
        };
 
        unsigned int BorderIndex = 0; 
        Gedim::GeometryUtilities::SegmentPolyhedronPositionResult::Types Type =
            Gedim::GeometryUtilities::SegmentPolyhedronPositionResult::Types::Unknown;
    };
 
    struct Polyhedron final
    {
        Eigen::MatrixXd Vertices;           
        Eigen::MatrixXi Edges;              
        std::vector<Eigen::MatrixXi> Faces; 
    };
 
    struct SplitPolygonWithPlaneResult final
    {
        enum struct Types
        {
            Unknown = 0,
            Split = 1,
            Positive = 2,
            Negative = 3,
            OnPlane = 4
        };
 
        std::vector<unsigned int> PositiveVertices;     
        std::vector<unsigned int> NegativeVertices;     
        std::vector<unsigned int> PointsOnPlane;        
        std::vector<Eigen::Vector3d> NewVertices;       
        std::vector<unsigned int> NewVerticesEdgeIndex; 
        Gedim::GeometryUtilities::SplitPolygonWithPlaneResult::Types Type =
            Gedim::GeometryUtilities::SplitPolygonWithPlaneResult::Types::Unknown; 
    };
 
    struct SplitPolyhedronWithPlaneResult
    {
        enum struct Types
        {
            Unknown = 0,
            Split = 1,
            None = 2
        };
 
        struct NewVertices
        {
            Eigen::MatrixXd Vertices;                          
            std::vector<unsigned int> NewVerticesOriginalEdge; 
        };
 
        struct NewEdges
        {
            Eigen::MatrixXi Edges;                  
            std::vector<int> NewEdgesOriginalEdges; 
            std::vector<int> NewEdgesOriginalFace; 
        };
 
        struct NewFaces
        {
            std::vector<Eigen::MatrixXi> Faces;
            std::vector<int> NewFacesOriginalFaces; 
        };
 
        struct NewPolyhedron
        {
            std::vector<unsigned int> Vertices;
            std::vector<unsigned int> Edges;
            std::vector<unsigned int> Faces;
        };
 
        std::vector<std::vector<unsigned int>> OriginalEdgesNewEdges;
        std::vector<std::vector<unsigned int>> OriginalFacesNewFaces;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::NewVertices Vertices;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::NewEdges Edges;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::NewFaces Faces;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::NewPolyhedron PositivePolyhedron;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::NewPolyhedron NegativePolyhedron;
        Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::Types Type =
            Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult::Types::Unknown; 
    };
 
    struct AlignedPolyhedronEdgesResult
    {
        std::vector<std::vector<unsigned int>> AlignedEdgesVertices;
        std::vector<std::vector<unsigned int>> AlignedEdgesEdges;
    };
 
    struct MergePolyhedronsInput
    {
        enum struct MergeTypes
        {
            None = 0,
            Common = 1,
            Remove = 2
        };
 
        constexpr static unsigned int none = UINT_MAX;
 
        using MergeTypesList = std::vector<std::pair<Gedim::GeometryUtilities::MergePolyhedronsInput::MergeTypes, unsigned int>>;
 
        std::array<Gedim::GeometryUtilities::MergePolyhedronsInput::MergeTypesList, 2> Vertices_Type;
        std::array<Gedim::GeometryUtilities::MergePolyhedronsInput::MergeTypesList, 2> Edges_Type;
        std::array<Gedim::GeometryUtilities::MergePolyhedronsInput::MergeTypesList, 2> Faces_Type;
        std::vector<std::array<unsigned int, 2>> Common_vertices;
        std::vector<std::array<unsigned int, 2>> Common_edges;
        std::vector<std::array<unsigned int, 2>> Common_faces;
    };
 
    struct MergePolyhedronsResult
    {
        constexpr static unsigned int none = UINT_MAX;
        std::array<std::vector<unsigned int>, 2> OriginalToMergedVertices;
        std::array<std::vector<unsigned int>, 2> OriginalToMergedEdges;
        std::array<std::vector<unsigned int>, 2> OriginalToMergedFaces;
        std::vector<std::array<unsigned int, 2>> MergedToOriginalVertices;
        std::vector<std::array<unsigned int, 2>> MergedToOriginalEdges;
        std::vector<std::array<unsigned int, 2>> MergedToOriginalFaces;
        Gedim::GeometryUtilities::Polyhedron MergedPolyhedron;
    };
 
  public:
    GeometryUtilities(const GeometryUtilitiesConfig &configuration);
    ~GeometryUtilities();
 
    inline double Tolerance1D() const
    {
        return std::max(_configuration.MinTolerance, _configuration.Tolerance1D);
    }
    inline double Tolerance1DSquared() const
    {
        return std::max(_configuration.MinTolerance * _configuration.MinTolerance,
                        _configuration.Tolerance1D * _configuration.Tolerance1D);
    }
    inline double Tolerance2D() const
    {
        return std::max(_configuration.MinTolerance,
                        std::max(_configuration.Tolerance2D, 0.25 * std::sqrt(3.0) * Tolerance1D() * Tolerance1D()));
    }
    inline double Tolerance3D() const
    {
        return std::max(_configuration.MinTolerance,
                        std::max(_configuration.Tolerance3D,
                                 std::max(std::sqrt(2.0) / 12.0 * Tolerance1D() * Tolerance1D() * Tolerance1D(),
                                          std::sqrt(6.0) / 9.0 * Tolerance2D() * Tolerance1D())));
    }
 
    inline double RelativeDifference(const double &first, const double &second, const double &tolerance) const
    {
        const double max_tolerance = std::max(std::abs(tolerance), _configuration.MinTolerance);
        return std::abs(second - first) / ((std::abs(first) <= max_tolerance) ? 1.0 : std::abs(first));
    }
 
    CompareTypes CompareValues(const double &first, const double &second, const double &tolerance) const;
 
    inline bool AreValuesEqual(const double &first, const double &second, const double &tolerance) const
    {
        return CompareValues(first, second, tolerance) == CompareTypes::Coincident;
    }
 
    inline bool IsValueGreater(const double &first, const double &second, const double &tolerance) const
    {
        return CompareValues(first, second, tolerance) == CompareTypes::SecondBeforeFirst;
    }
 
    inline bool IsValueGreaterOrEqual(const double &first, const double &second, const double &tolerance) const
    {
        const CompareTypes result = CompareValues(first, second, tolerance);
 
        return result == CompareTypes::SecondBeforeFirst || result == CompareTypes::Coincident;
    }
 
    inline bool IsValueLower(const double &first, const double &second, const double &tolerance) const
    {
        return !IsValueGreaterOrEqual(first, second, tolerance);
    }
 
    inline bool IsValueLowerOrEqual(const double &first, const double &second, const double &tolerance) const
    {
        return !IsValueGreater(first, second, tolerance);
    }
 
    inline bool IsValuePositive(const double &value, const double &tolerance) const
    {
        return CompareValues(0.0, value, tolerance) == CompareTypes::FirstBeforeSecond;
    }
 
    inline bool IsValueNegative(const double &value, const double &tolerance) const
    {
        return CompareValues(0.0, value, tolerance) == CompareTypes::SecondBeforeFirst;
    }
 
    inline bool IsValueZero(const double &value, const double &tolerance) const
    {
        return CompareValues(0.0, value, tolerance) == CompareTypes::Coincident;
    }
 
    Eigen::MatrixXd fibonacci_sphere(const unsigned int num_points) const;
    Eigen::MatrixXd generate_uniform_random_points_in_sphere(const unsigned int num_points, const double radius = 1.0) const;
 
    std::vector<double> EquispaceCoordinates(const double &step, const bool &insertExtremes) const;
 
    std::vector<double> EquispaceCoordinates(const unsigned int &size, const double &origin, const double &end, const bool &insertExtremes) const;
    std::vector<double> RandomCoordinates(const unsigned int size,
                                          const bool insertExtremes,
                                          const double &minDistance,
                                          const unsigned int seed = time(nullptr)) const;
 
    inline double PolarAngle(const Eigen::Vector3d &v_prev,
                             const Eigen::Vector3d &v,
                             const Eigen::Vector3d &v_next,
                             const double &norm_v_prev_v,
                             const double &norm_v_next_v) const
    {
        return IsValueZero(norm_v_prev_v, Tolerance1D()) || IsValueZero(norm_v_next_v, Tolerance1D())
                   ? 0.0
                   : (v.x() - v_prev.x()) * (v_next.y() - v_prev.y()) / (norm_v_prev_v * norm_v_next_v) -
                         (v_next.x() - v_prev.x()) * (v.y() - v_prev.y()) / (norm_v_prev_v * norm_v_next_v);
    }
 
    inline double Angle(const Eigen::Vector3d &v_prev,
                        const Eigen::Vector3d &v,
                        const Eigen::Vector3d &v_next,
                        const double &norm_v_prev_v,
                        const double &norm_v_next_v) const
    {
        return IsValueZero(norm_v_prev_v, Tolerance1D()) || IsValueZero(norm_v_next_v, Tolerance1D())
                   ? 0.0
                   : std::atan2(((v_next.x() - v.x()) * (v_prev.y() - v.y()) - (v_next.y() - v.y()) * (v_prev.x() - v.x())),
                                (v_prev.x() - v.x()) * (v_next.x() - v.x()) + (v_prev.y() - v.y()) * (v_next.y() - v.y()));
    }
 
    inline double AngleRadiansToDegree(const double &rad_angle) const
    {
        return std::atan2(std::sin(rad_angle), std::cos(rad_angle)) * 180.0 / std::numbers::pi;
    }
 
    inline double AngleDegreeToRadians(const double &deg_angle) const
    {
        return IsValueZero(deg_angle, Tolerance1D()) ? 0.0
               : IsValueZero(std::fmod(deg_angle, 180.0), Tolerance1D())
                   ? IsValuePositive(deg_angle, Tolerance1D()) ? std::numbers::pi : -std::numbers::pi
                   : std::atan2(std::sin(deg_angle * std::numbers::pi / 180.0), std::cos(deg_angle * std::numbers::pi / 180.0));
    }
 
    inline double PointDistance(const Eigen::Vector3d &firstPoint, const Eigen::Vector3d &secondPoint) const
    {
        return (secondPoint - firstPoint).norm();
    }
 
    Eigen::VectorXd PointDistances(const Eigen::MatrixXd &points, const Eigen::Vector3d &point) const;
 
    Eigen::MatrixXd PointsDistance(const Eigen::MatrixXd &points) const;
 
    Eigen::MatrixXd PointsBoundingBox(const Eigen::MatrixXd &points) const;
 
    inline bool IsPointInBoundingBox(const Eigen::Vector3d &point, const Eigen::MatrixXd &boudingBox) const
    {
        return (IsValueGreaterOrEqual(point.x(), boudingBox(0, 0), Tolerance1D()) &&
                IsValueGreaterOrEqual(boudingBox(0, 1), point.x(), Tolerance1D())) &&
               (IsValueGreaterOrEqual(point.y(), boudingBox(1, 0), Tolerance1D()) &&
                IsValueGreaterOrEqual(boudingBox(1, 1), point.y(), Tolerance1D())) &&
               (IsValueGreaterOrEqual(point.z(), boudingBox(2, 0), Tolerance1D()) &&
                IsValueGreaterOrEqual(boudingBox(2, 1), point.z(), Tolerance1D()));
    }
 
    inline bool BoundingBoxesIntersects(const Eigen::MatrixXd &boudingBox_1, const Eigen::MatrixXd &boudingBox_2) const
    {
        return IsValueLowerOrEqual(boudingBox_1(0, 0), boudingBox_2(0, 1), Tolerance1D()) &&
               IsValueGreaterOrEqual(boudingBox_1(0, 1), boudingBox_2(0, 0), Tolerance1D()) &&
               IsValueLowerOrEqual(boudingBox_1(1, 0), boudingBox_2(1, 1), Tolerance1D()) &&
               IsValueGreaterOrEqual(boudingBox_1(1, 1), boudingBox_2(1, 0), Tolerance1D()) &&
               IsValueLowerOrEqual(boudingBox_1(2, 0), boudingBox_2(2, 1), Tolerance1D()) &&
               IsValueGreaterOrEqual(boudingBox_1(2, 1), boudingBox_2(2, 0), Tolerance1D());
    }
 
    double PointsMaxDistance(const Eigen::MatrixXd &points) const;
 
    inline bool IsPointZero(const Eigen::Vector3d &point) const
    {
        return IsValueZero(PointDistance(Eigen::Vector3d::Zero(), point), Tolerance1D());
    }
 
    inline bool PointsAreCoincident(const Eigen::Vector3d &firstPoint, const Eigen::Vector3d &secondPoint) const
    {
        return IsValueZero(PointDistance(firstPoint, secondPoint), Tolerance1D());
    }
 
    std::vector<unsigned int> FindPointInPoints(const Eigen::MatrixXd &points, const Eigen::Vector3d &point) const;
 
    inline double PointLineDistance(const Eigen::Vector3d &point, const Eigen::Vector3d &lineOrigin, const Eigen::Vector3d &normalToLine) const
    {
        return std::abs(normalToLine.dot(point - lineOrigin)) / normalToLine.norm();
    }
 
    inline bool PointIs2D(const Eigen::Vector3d &point) const
    {
        return PointsAre2D(point);
    }
 
    inline bool PointsAre2D(const Eigen::MatrixXd &points) const
    {
        Gedim::Output::Assert(points.rows() == 3 && points.cols() > 0);
        return points.row(2).isZero(_configuration.Tolerance1D);
    }
 
    inline double PointCurvilinearCoordinate(const Eigen::Vector3d &point,
                                             const Eigen::Vector3d &segmentOrigin,
                                             const Eigen::Vector3d &segmentEnd) const
    {
        const Eigen::Vector3d segmentTangent = (segmentEnd - segmentOrigin);
        return PointLineCurvilinearCoordinate(point, segmentOrigin, segmentTangent, segmentTangent.squaredNorm());
    }
 
    inline double PointLineCurvilinearCoordinate(const Eigen::Vector3d &point,
                                                 const Eigen::Vector3d &lineOrigin,
                                                 const Eigen::Vector3d &lineTangent,
                                                 const double &lineTangentSquaredLength) const
    {
        return (point - lineOrigin).dot(lineTangent) / lineTangentSquaredLength;
    }
 
    bool IsPointOnLine(const Eigen::Vector3d &point,
                       const Eigen::Vector3d &lineOrigin,
                       const Eigen::Vector3d &lineTangent,
                       const double &lineTangentSquaredLength) const;
 
    PointSegmentPositionTypes PointSegmentPosition(const Eigen::Vector3d &point,
                                                   const Eigen::Vector3d &segmentOrigin,
                                                   const Eigen::Vector3d &segmentEnd) const;
 
    PointSegmentPositionTypes PointSegmentPosition(const double &curvilinearCoordinate) const;
 
    inline double PointSegmentProjection(const Eigen::Vector3d &point,
                                         const Eigen::Vector3d &segmentOrigin,
                                         const Eigen::Vector3d &segmentEnd) const
    {
        return PointCurvilinearCoordinate(point, segmentOrigin, segmentEnd);
    }
 
    double PointPlaneDistance(const Eigen::Vector3d &point, const std::array<Eigen::Vector3d, 3> &planePoints) const;
    double PointPlaneDistance(const Eigen::Vector3d &point, const Eigen::Vector3d &planeNormal, const Eigen::Vector3d &planeOrigin) const;
 
    PointPlanePositionTypes PointPlanePosition(const double &pointPlaneDistance) const;
 
    inline bool IsPointOnPlane(const double &pointPlaneDistance) const
    {
        return PointPlanePosition(pointPlaneDistance) == PointPlanePositionTypes::OnPlane;
    }
 
    inline bool IsPointOnPlane(const Eigen::Vector3d &point, const Eigen::Vector3d &planeNormal, const Eigen::Vector3d &planeOrigin) const
    {
        return PointPlanePosition(PointPlaneDistance(point, planeNormal, planeOrigin)) == PointPlanePositionTypes::OnPlane;
    }
 
    inline double SegmentLength(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        return (segmentEnd - segmentOrigin).norm();
    }
 
    inline Eigen::Vector3d SegmentTangent(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        return segmentEnd - segmentOrigin;
    }
 
    inline Eigen::Vector3d SegmentNormal(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        Gedim::Output::Assert(PointsAre2D(segmentOrigin) && PointsAre2D(segmentEnd));
        Eigen::Vector3d tangent = SegmentTangent(segmentOrigin, segmentEnd).normalized();
        return Eigen::Vector3d(tangent.y(), -tangent.x(), 0.0);
    }
 
    inline double SegmentSlope(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        Gedim::Output::Assert(!AreValuesEqual(segmentEnd.x(), segmentOrigin.x(), Tolerance1D()));
        return (segmentEnd.y() - segmentOrigin.y()) / (segmentEnd.x() - segmentOrigin.x());
    }
 
    inline double SegmentIntercept(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        Gedim::Output::Assert(!AreValuesEqual(segmentEnd.x(), segmentOrigin.x(), Tolerance1D()));
        return segmentOrigin.y() -
               segmentOrigin.x() * (segmentEnd.y() - segmentOrigin.y()) / (segmentEnd.x() - segmentOrigin.x());
    }
 
    Eigen::MatrixXi MakeConcatenation(const Eigen::MatrixXi &segments, const unsigned int starting_vertex) const;
 
    inline bool CheckNoSpheresIntersection(const Eigen::Vector3d &firstSphereCenter,
                                           const Eigen::Vector3d &secondSphereCenter,
                                           const double &firstSphereDiameter,
                                           const double &secondSphereDiameter) const
    {
        return IsValueGreater(2.0 * PointDistance(firstSphereCenter, secondSphereCenter),
                              firstSphereDiameter + secondSphereDiameter,
                              Tolerance1D());
    }
 
    bool CheckTrianglesIntersection(const Eigen::MatrixX3d &triangle_one,
                                    const Eigen::MatrixX3d &triangle_two,
                                    const bool admit_boundary = true) const;
 
    bool AreLineCoplanar(const Eigen::Vector3d &firstLineOrigin,
                         const Eigen::Vector3d &firstLineTangent,
                         const Eigen::Vector3d &secondLineOrigin,
                         const Eigen::Vector3d &secondLineTangent) const;
 
    bool IsPolygonCoplanar(const Eigen::Vector3d &planeNormal,
                           const Eigen::Vector3d &planeOrigin,
                           const Eigen::MatrixXd &polygonVertices,
                           const std::vector<unsigned int> &polygonUnalignedVertices) const;
 
    IntersectionSegmentSegmentResult IntersectionSegmentSegment(const Eigen::Vector3d &firstSegmentOrigin,
                                                                const Eigen::Vector3d &firstSegmentEnd,
                                                                const Eigen::Vector3d &secondSegmentOrigin,
                                                                const Eigen::Vector3d &secondSegmentEnd) const;
 
    std::vector<std::list<double>> IntersectionsBetweenSegments(const std::vector<Eigen::MatrixXd> &segmentsVertices,
                                                                const std::vector<Eigen::Vector3d> &segmentsTangent,
                                                                const std::vector<Eigen::Vector3d> &segmentsBarycenter,
                                                                const std::vector<double> &segmentsLength) const;
 
    Gedim::GeometryUtilities::IntersectionSegmentCircleResult IntersectionSegmentCircle(const Eigen::Vector3d &segmentOrigin,
                                                                                        const Eigen::Vector3d &segmentEnd,
                                                                                        const Eigen::Vector3d &circleCenter,
                                                                                        const double &circleRadius) const;
 
    Gedim::GeometryUtilities::IntersectionSegmentPlaneResult IntersectionSegmentPlane(const Eigen::Vector3d &segmentOrigin,
                                                                                      const Eigen::Vector3d &segmentEnd,
                                                                                      const Eigen::Vector3d &planeNormal,
                                                                                      const Eigen::Vector3d &planeOrigin) const;
 
    Gedim::GeometryUtilities::IntersectionPolyhedronPlaneResult IntersectionPolyhedronPlane(
        const Eigen::MatrixXd &polyhedronVertices,
        const Eigen::MatrixXi &polyhedronEdges,
        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
        const Eigen::Vector3d &planeNormal,
        const Eigen::Vector3d &planeOrigin,
        const Eigen::Matrix3d &planeRotationMatrix,
        const Eigen::Vector3d &planeTranslation) const;
 
    Gedim::GeometryUtilities::SplitPolyhedronWithPlaneResult SplitPolyhedronWithPlane(
        const Eigen::MatrixXd &polyhedronVertices,
        const Eigen::MatrixXi &polyhedronEdges,
        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
        const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
        const std::vector<Eigen::MatrixXd> &polyhedronFaceEdgeTangents,
        const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
        const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices,
        const Eigen::Vector3d &planeNormal,
        const Eigen::Vector3d &planeOrigin,
        const Eigen::Matrix3d &planeRotationMatrix,
        const Eigen::Vector3d &planeTranslation) const;
 
    std::vector<GeometryUtilities::Polyhedron> SplitPolyhedronWithPlaneResultToPolyhedra(const SplitPolyhedronWithPlaneResult &result) const;
 
    IntersectionPolyhedronLineResult IntersectionPolyhedronLine(const Eigen::MatrixXd &polyhedronVertices,
                                                                const Eigen::MatrixXi &polyhedronEdges,
                                                                const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                                const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
                                                                const std::vector<bool> &polyhedronFaceNormalDirections,
                                                                const Eigen::Vector3d &lineTangent,
                                                                const Eigen::Vector3d &lineOrigin) const;
 
    Gedim::GeometryUtilities::IntersectionPolyhedronLineResult IntersectionPolyhedronSegment(
        const Eigen::MatrixXd &polyhedronVertices,
        const Eigen::MatrixXi &polyhedronEdges,
        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
        const Eigen::Vector3d &segmentOrigin,
        const Eigen::Vector3d &segmentEnd,
        const Eigen::Vector3d &segmentTangent,
        const IntersectionPolyhedronLineResult &polyhedronLineIntersections) const;
 
    Gedim::GeometryUtilities::IntersectionPolyhedronsSegmentResult IntersectionPolyhedronsSegment(
        const std::vector<Polyhedron> &polyhedrons,
        const std::vector<std::vector<Eigen::Vector3d>> &polyhedronFaceNormals,
        const std::vector<std::vector<bool>> &polyhedronFaceNormalDirections,
        const Eigen::Vector3d &segmentOrigin,
        const Eigen::Vector3d &segmentEnd,
        const Eigen::Vector3d &segmentTangent) const;
 
    Gedim::GeometryUtilities::PointPolygonPositionResult PointPolygonPosition(const Eigen::Vector3d &point,
                                                                              const Eigen::MatrixXd &polygonVertices) const;
 
    Gedim::GeometryUtilities::PointPolygonPositionResult PointPolygonPosition_RayCasting(const Eigen::Vector3d &point,
                                                                                         const Eigen::MatrixXd &polygonVertices) const;
 
    inline bool IsPointInsidePolygon(const Eigen::Vector3d &point, const Eigen::MatrixXd &polygonVertices) const
    {
        return !(PointPolygonPosition(point, polygonVertices).Type ==
                 Gedim::GeometryUtilities::PointPolygonPositionResult::Types::Outside);
    }
 
    inline bool IsPointInsidePolygon_RayCasting(const Eigen::Vector3d &point, const Eigen::MatrixXd &polygonVertices) const
    {
        return !(PointPolygonPosition_RayCasting(point, polygonVertices).Type ==
                 Gedim::GeometryUtilities::PointPolygonPositionResult::Types::Outside);
    }
 
    LinePolygonPositionResult LinePolygonPosition(const Eigen::Vector3d &lineTangent,
                                                  const Eigen::Vector3d &lineOrigin,
                                                  const Eigen::MatrixXd &polygonVertices) const;
 
 
    Gedim::GeometryUtilities::PointPolyhedronPositionResult PointPolyhedronPosition(
        const Eigen::Vector3d &point,
        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
        const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
        const std::vector<Eigen::MatrixXd> &polyhedronFaceRotatedVertices,
        const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
        const std::vector<bool> &polyhedronFaceNormalDirections,
        const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
        const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices) const;
    Gedim::GeometryUtilities::PointPolyhedronPositionResult PointPolyhedronPosition(
        const Eigen::Vector3d &point,
        const std::vector<Eigen::MatrixXi> &polyhedron_faces,
        const std::vector<Eigen::MatrixXd> &polyhedron_faces_3D_vertices,
        const std::vector<Eigen::MatrixXd> &polyhedron_faces_2D_vertices,
        const std::vector<Eigen::Vector3d> &polyhedron_faces_normals,
        const std::vector<bool> &polyhedron_faces_normal_direction,
        const std::vector<Eigen::Vector3d> &polyhedron_faces_translation,
        const std::vector<Eigen::Matrix3d> &polyhedron_faces_rotation_matrix,
        const std::vector<Eigen::MatrixXd> &polyhedron_tetrahedrons) const;
 
    bool IsPointInsideTetrahedron(const Eigen::MatrixXd &tetrahedron, const Eigen::Vector3d &point) const;
 
    PointCirclePositionResult PointCirclePosition(const Eigen::Vector3d &point,
                                                  const Eigen::Vector3d &circleCenter,
                                                  const double &circleRadius) const;
 
    std::vector<PointCirclePositionResult> PointCirclePositions(const Eigen::MatrixXd &points,
                                                                const Eigen::Vector3d &circleCenter,
                                                                const double &circleRadius) const;
 
    Gedim::GeometryUtilities::PolygonCirclePositionTypes PolygonCirclePosition(
        const Eigen::MatrixXd &polygonVertices,
        const Eigen::Vector3d &circleCenter,
        const double &circleRadius,
        const std::vector<PointCirclePositionResult> &vertexPositions,
        const IntersectionPolygonCircleResult &polygonCircleIntersections) const;
 
    Gedim::GeometryUtilities::IntersectionPolygonCircleResult IntersectionPolygonCircle(const Eigen::MatrixXd &polygonVertices,
                                                                                        const Eigen::Vector3d &circleCenter,
                                                                                        const double &circleRadius) const;
 
    std::vector<unsigned int> PolygonTriangulationByFirstVertex(const Eigen::MatrixXd &polygonVertices) const;
 
    std::vector<unsigned int> PolygonTriangulationByEarClipping(const Eigen::MatrixXd &polygonVertices) const;
 
    std::vector<unsigned int> PolygonTriangulationByInternalPoint(const Eigen::MatrixXd &polygonVertices,
                                                                  const Eigen::Vector3d &internalPoint) const;
 
    PolygonDivisionByAngleQuadrantResult PolygonOutsideCircleDivisionByAngleQuadrant(const Eigen::MatrixXd &polygonVertices,
                                                                                     const Eigen::MatrixXd &polygonEdgeTangents,
                                                                                     const Eigen::Vector3d &circleCenter,
                                                                                     const double &circleRadius,
                                                                                     const unsigned int &curvedEdgeIndex) const;
 
    PolygonDivisionByAngleQuadrantResult PolygonInsideCircleDivisionByAngleQuadrant(const Eigen::MatrixXd &polygonVertices,
                                                                                    const Eigen::MatrixXd &polygonEdgeTangents,
                                                                                    const Eigen::Vector3d &circleCenter,
                                                                                    const double &circleRadius,
                                                                                    const unsigned int &curvedEdgeIndex) const;
 
    PolygonDivisionByCircleResult PolygonDivisionByCircle(const Eigen::MatrixXd &polygonVertices,
                                                          const Eigen::MatrixXd &polygonEdgeTangents,
                                                          const Eigen::Vector3d &circleCenter,
                                                          const double &circleRadius,
                                                          const unsigned int &curvedEdgeIndex) const;
 
    CircleDivisionByPolygonResult CircleDivisionByPolygon(const Eigen::MatrixXd &polygonVertices,
                                                          const Eigen::MatrixXd &polygonEdgeTangents,
                                                          const Eigen::Vector3d &circleCenter,
                                                          const double &circleRadius,
                                                          const unsigned int &curvedEdgeIndex) const;
 
    double PolygonArea(const Eigen::MatrixXd &polygonVertices) const;
 
    double PolygonArea3D(const Eigen::MatrixXd &polygonVertices) const;
 
    Eigen::Matrix2d PolygonMass(const Eigen::Vector3d &polygonCentroid,
                                const std::vector<Eigen::Matrix3d> &polygonTriangulationPoints) const;
 
    Eigen::Matrix3d PolygonInertia(const Eigen::Vector3d &polygonCentroid,
                                   const std::vector<Eigen::Matrix3d> &polygonTriangulationPoints) const;
 
    Gedim::GeometryUtilities::SplitPolygonWithSegmentResult SplitPolygonWithSegment(const SplitPolygonInput &input) const;
 
    Gedim::GeometryUtilities::SplitPolygonWithCircleResult SplitPolygonWithCircle(
        const Eigen::MatrixXd &polygonVertices,
        const Eigen::Vector3d &circleCenter,
        const double &circleRadius,
        const std::vector<PointCirclePositionResult> &vertexPositions,
        const IntersectionPolygonCircleResult &polygonCircleIntersections,
        const PolygonCirclePositionTypes &polygonCirclePosition) const;
 
    Eigen::MatrixXd SplitPolygonWithCircleBuildSubPolygon(const SplitPolygonWithCircleResult &splitResult,
                                                          const unsigned int &subPolygonIndex,
                                                          const Eigen::MatrixXd &polygonVertices,
                                                          const Gedim::GeometryUtilities::IntersectionPolygonCircleResult &polygonCircleIntersections) const;
 
    Gedim::GeometryUtilities::SplitPolygonWithPlaneResult SplitPolygonWithPlane(const Eigen::MatrixXd &polygonVertices,
                                                                                const Eigen::MatrixXd &polygonEdgeTangents,
                                                                                const Eigen::Vector3d &planeNormal,
                                                                                const Eigen::Vector3d &planeOrigin,
                                                                                const Eigen::Vector3d &polygonTranslation,
                                                                                const Eigen::Matrix3d &polygonRotationMatrix) const;
 
    Eigen::Vector3d PolygonNormal(const Eigen::MatrixXd &polygonVertices) const;
 
    std::array<Eigen::Vector3d, 2> PolygonTangents(const Eigen::MatrixXd &polygonVertices, const Eigen::Vector3d &polygonNormal) const;
 
    Eigen::MatrixXd PolygonEdgesCentroid(const Eigen::MatrixXd &polygonVertices) const;
 
    Eigen::VectorXd PolygonEdgeLengths(const Eigen::MatrixXd &polygonVertices) const;
 
    Eigen::MatrixXd PolygonEdgeTangents(const Eigen::MatrixXd &polygonVertices) const;
 
    Eigen::MatrixXd PolygonEdgeNormals(const Eigen::MatrixXd &polygonVertices) const;
 
    inline Eigen::Vector3d SimplexBarycenter(const Eigen::MatrixXd &vertices) const
    {
        Gedim::Output::Assert(vertices.rows() == 3);
        return vertices.rowwise().mean();
    }
 
    inline double SimplexMeasure(const Eigen::MatrixXd &vertices) const
    {
        Gedim::Output::Assert(vertices.rows() == 3 && vertices.cols() < 5);
 
        Eigen::MatrixXd measure_matrix(vertices.cols(), 4);
        measure_matrix.block(0, 0, vertices.cols(), 3) = vertices.transpose();
        measure_matrix.col(3).setOnes();
        const double div_dim = vertices.cols() == 3 ? 0.5 : 1.0 / 6.0;
 
        return div_dim * std::abs(measure_matrix.determinant());
    }
 
    inline Eigen::Vector3d SegmentBarycenter(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd) const
    {
        return SimplexBarycenter((Eigen::MatrixXd(3, 2) << segmentOrigin, segmentEnd).finished());
    }
 
    inline Eigen::Vector3d PolygonBarycenter(const Eigen::MatrixXd &polygonVertices) const
    {
        Gedim::Output::Assert(polygonVertices.rows() == 3 && polygonVertices.cols() > 2);
        return SimplexBarycenter(polygonVertices);
    }
 
    inline Eigen::Vector3d PolyhedronBarycenter(const Eigen::MatrixXd &polyhedronVertices) const
    {
        Gedim::Output::Assert(polyhedronVertices.rows() == 3 && polyhedronVertices.cols() > 2);
        return SimplexBarycenter(polyhedronVertices);
    }
 
    Eigen::Vector3d PolygonCentroid(const Eigen::MatrixXd &polygonVertices, const double &polygonArea) const;
 
    inline Eigen::Vector3d PolygonCentroid(const Eigen::MatrixXd &subPolygonCentroids,
                                           const Eigen::VectorXd &subPolygonAreas,
                                           const double &polygonArea) const
    {
        Gedim::Output::Assert(subPolygonCentroids.rows() == 3 && subPolygonCentroids.cols() > 0 &&
                              subPolygonCentroids.cols() == subPolygonAreas.size());
 
        return subPolygonCentroids * subPolygonAreas / polygonArea;
    }
 
    double PolygonAreaByBoundaryIntegral(
        const Eigen::MatrixXd &polygonVertices,
        const Eigen::VectorXd &edgeLengths,
        const Eigen::MatrixXd &edgeTangents,
        const Eigen::MatrixXd &edgeNormals,
        const Eigen::MatrixXd &referenceQuadraturePoints = (Eigen::MatrixXd(3, 1) << 0.5, 0.0, 0.0).finished(),
        const Eigen::VectorXd &referenceQuadratureWeights = Eigen::VectorXd::Ones(1)) const;
 
    double PolygonAreaByInternalIntegral(const std::vector<Eigen::Matrix3d> &polygonTriangulationPoints,
                                         const Eigen::VectorXd &referenceQuadratureWeights = Eigen::VectorXd::Constant(1, 0.5)) const;
 
    Eigen::Vector3d PolygonCentroidByIntegral(const Eigen::MatrixXd &polygonVertices,
                                              const Eigen::VectorXd &edgeLengths,
                                              const Eigen::MatrixXd &edgeTangents,
                                              const Eigen::MatrixXd &edgeNormals,
                                              const double &polygonArea) const;
 
    Eigen::VectorXd PolygonCentroidEdgesDistance(const Eigen::MatrixXd &polygonVertices,
                                                 const Eigen::Vector3d &polygonCentroid,
                                                 const Eigen::MatrixXd &polygonEdgeNormals) const;
 
    Eigen::VectorXd PolygonCentroidVerticesDistance(const Eigen::MatrixXd &polygonVertices, const Eigen::Vector3d &polygonCentroid) const;
 
    inline double PolygonInRadius(const Eigen::VectorXd &polygonCentroidEdgesDistance) const
    {
        return polygonCentroidEdgesDistance.minCoeff();
    }
 
    inline double PolygonAspectRatio(const double &polygonDiameter, const double &polygonInRadius) const
    {
        return polygonDiameter / (2.0 * polygonInRadius);
    };
 
    inline double PolygonDiameter(const Eigen::MatrixXd &polygonVertices) const
    {
        return PointsMaxDistance(polygonVertices);
    }
 
    inline Eigen::Vector3d PolygonTranslation(const Eigen::MatrixXd &polygonVertices) const
    {
        return polygonVertices.col(0);
    }
 
    Eigen::Matrix3d PolygonRotationMatrix(const Eigen::MatrixXd &polygonVertices,
                                          const Eigen::Vector3d &polygonNormal,
                                          const Eigen::Vector3d &polygonTranslation) const;
 
    bool PolygonIsConvex(const Eigen::MatrixXd &polygonVertices, const Eigen::MatrixXd &convexHull) const;
 
    PolygonTypes PolygonType(const unsigned int &numPolygonVertices, const bool &isPolygonConvex) const;
 
    PolygonOrientations PolygonOrientation(const std::vector<unsigned int> &convexHull) const;
 
    inline std::vector<unsigned int> ChangePolygonOrientation(const unsigned int numPolygonVertices) const
    {
        std::vector<unsigned int> newVertices(numPolygonVertices);
        newVertices[0] = 0;
        for (unsigned int i = 1; i < numPolygonVertices; i++)
            newVertices[i] = numPolygonVertices - i;
        return newVertices;
    }
 
    Eigen::Matrix3d PlaneRotationMatrix(const Eigen::Vector3d &planeNormal) const;
 
    inline Eigen::Vector3d PlaneTranslation(const Eigen::Vector3d &planeOrigin) const
    {
        return planeOrigin;
    }
 
    Eigen::Matrix3d PlaneReflectionMatrix(const Eigen::Vector3d &plane_normal) const;
    Eigen::Vector3d PlaneReflectionTranslation(const Eigen::Vector3d &plane_normal, const Eigen::Vector3d &planeOrigin) const;
 
    Eigen::MatrixXd PlaneToPolygon(const Eigen::Vector3d &plane_normal, const Eigen::Vector3d &plane_origin) const;
 
    inline Eigen::MatrixXd RotatePoints(const Eigen::MatrixXd &points,
                                        const Eigen::Matrix3d &rotationMatrix,
                                        const Eigen::Vector3d &translation = Eigen::Vector3d::Zero()) const
    {
        Gedim::Output::Assert(points.rows() == 3 && points.cols() > 0);
        return (rotationMatrix * points).colwise() + translation;
    }
 
    inline Eigen::MatrixXd RotatePointsFrom2DTo3D(const Eigen::MatrixXd &points,
                                                  const Eigen::Matrix3d &rotationMatrix,
                                                  const Eigen::Vector3d &translation = Eigen::Vector3d::Zero()) const
    {
        Gedim::Output::Assert(points.rows() == 3 && points.cols() > 0 && PointsAre2D(points));
        return (rotationMatrix * points).colwise() + translation;
    }
    inline Eigen::MatrixXd RotatePointsFrom3DTo2D(const Eigen::MatrixXd &points,
                                                  const Eigen::Matrix3d &rotationMatrix,
                                                  const Eigen::Vector3d &translation = Eigen::Vector3d::Zero()) const
    {
        Gedim::Output::Assert(points.rows() == 3 && points.cols() > 0);
        Eigen::MatrixXd rotatedPoints = rotationMatrix * (points.colwise() - translation);
        rotatedPoints.row(2).setZero();
        return rotatedPoints;
    }
 
    std::vector<unsigned int> ConvexHull(const Eigen::MatrixXd &points, const bool &includeCollinear = true) const;
 
    inline std::vector<bool> PointsAreAligned(const Eigen::Vector3d &segmentOrigin,
                                              const Eigen::Vector3d &segmentEnd,
                                              const Eigen::MatrixXd &points) const
    {
        return PointsAreOnLine(points, segmentOrigin, SegmentTangent(segmentOrigin, segmentEnd));
    }
 
    std::vector<bool> PointsAreOnLine(const Eigen::MatrixXd &points,
                                      const Eigen::Vector3d &lineOrigin,
                                      const Eigen::Vector3d &lineTangent) const;
 
    inline bool PointIsAligned(const Eigen::Vector3d &segmentOrigin, const Eigen::Vector3d &segmentEnd, const Eigen::Vector3d &point) const
    {
        return PointsAreAligned(segmentOrigin, segmentEnd, point)[0];
    }
 
    inline bool PointIsOnLine(const Eigen::Vector3d &point, const Eigen::Vector3d &lineOrigin, const Eigen::Vector3d &lineTangent) const
    {
        return PointsAreOnLine(point, lineOrigin, lineTangent)[0];
    }
 
    std::array<double, 3> PointToBarycentricCoordinates2D(const Eigen::Matrix3d &triangle, const Eigen::Vector3d &point) const;
    Eigen::Vector3d BarycentricCoordinatesToPoint2D(const Eigen::Matrix3d &triangle,
                                                    const std::array<double, 3> &barycentric_coordinates) const;
 
    std::vector<unsigned int> UnalignedPoints(const Eigen::MatrixXd &points, const unsigned int numDesiredUnalignedPoints = 0) const;
 
    std::vector<unsigned int> UnalignedPolyhedronPoints(const Eigen::MatrixXd &polyhedronVertices,
                                                        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                        const std::vector<Eigen::Vector3d> &polyhedronFacesTranslation,
                                                        const std::vector<Eigen::Matrix3d> &polyhedronFacesRotationMatrix,
                                                        const std::vector<std::vector<unsigned int>> &polyhedronUnaligedFaces,
                                                        const std::vector<std::vector<unsigned int>> &polyhedronFacesUnalignedVertices) const;
 
    AlignedPolyhedronEdgesResult AlignedPolyhedronEdges(const Eigen::MatrixXd &polyhedronVertices,
                                                        const std::vector<std::vector<unsigned int>> &verticesAdjacency,
                                                        const std::vector<std::vector<unsigned int>> &edgesAdjacency,
                                                        const std::vector<std::unordered_map<unsigned int, unsigned int>> &adjacencyVerticesMap,
                                                        const Eigen::MatrixXd &polyhedronEdgeTangents,
                                                        const Eigen::VectorXd &polyhedronEdgeSquaredLenghts) const;
 
    Eigen::MatrixXd ExtractPoints(const Eigen::MatrixXd &points, const std::vector<unsigned int> &filter) const;
 
    std::vector<Eigen::Matrix3d> ExtractTriangulationPoints(const Eigen::MatrixXd &points,
                                                            const std::vector<unsigned int> &pointsTriangulation) const;
 
    std::vector<Eigen::Matrix3d> ExtractTriangulationPointsByInternalPoint(const Eigen::MatrixXd &points,
                                                                           const Eigen::Vector3d &internalPoint,
                                                                           const std::vector<unsigned int> &pointsTriangulation) const;
 
    Eigen::MatrixXd CreateEllipse(const double &axisMajorLength, const double &axisMinorLength, const unsigned int &resolution) const;
 
    Eigen::MatrixXd CreateTriangle(const Eigen::Vector3d &p1, const Eigen::Vector3d &p2, const Eigen::Vector3d &p3) const;
    Eigen::MatrixXd CreateParallelogram(const Eigen::Vector3d &origin,
                                        const Eigen::Vector3d &lengthVector,
                                        const Eigen::Vector3d &widthVector) const;
    inline Eigen::MatrixXd CreateRectangle(const Eigen::Vector3d &origin, const double &base, const double &height) const
    {
        return CreateParallelogram(origin, Eigen::Vector3d(base, 0.0, 0.0), Eigen::Vector3d(0.0, height, 0.0));
    }
 
    inline Eigen::MatrixXd CreateSquare(const Eigen::Vector3d &origin, const double &edgeLength) const
    {
        return CreateRectangle(origin, edgeLength, edgeLength);
    }
 
    Polyhedron CreateTetrahedronWithOrigin(const Eigen::Vector3d &origin,
                                           const Eigen::Vector3d &lengthVector,
                                           const Eigen::Vector3d &heightVector,
                                           const Eigen::Vector3d &widthVector) const;
 
    Polyhedron CreateTetrahedronWithVertices(const Eigen::Vector3d &v1,
                                             const Eigen::Vector3d &v2,
                                             const Eigen::Vector3d &v3,
                                             const Eigen::Vector3d &v4) const;
 
    Polyhedron CreateParallelepipedWithOrigin(const Eigen::Vector3d &origin,
                                              const Eigen::Vector3d &lengthVector,
                                              const Eigen::Vector3d &heightVector,
                                              const Eigen::Vector3d &widthVector) const;
 
    inline Polyhedron CreatePolyhedronWithExtrusion(const Eigen::MatrixXd &polygonVertices, const Eigen::Vector3d &heightVector) const
    {
        return CreatePolyhedronWithExtrusion(polygonVertices, std::vector<Eigen::Vector3d>(polygonVertices.cols(), heightVector));
    }
 
    Polyhedron CreatePolyhedronWithExtrusion(const Eigen::MatrixXd &polygonVertices,
                                             const std::vector<Eigen::Vector3d> &heightVectors) const;
 
    inline Polyhedron CreateCubeWithOrigin(const Eigen::Vector3d &origin, const double &edgeLength) const
    {
        return CreateParallelepipedWithOrigin(origin,
                                              Eigen::Vector3d(edgeLength, 0.0, 0.0),
                                              Eigen::Vector3d(0.0, 0.0, edgeLength),
                                              Eigen::Vector3d(0.0, edgeLength, 0.0));
    }
 
    double PolyhedronVolumeByBoundaryIntegral(
        const std::vector<std::vector<Eigen::Matrix3d>> &polyhedronFaceRotatedTriangulationPoints,
        const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
        const std::vector<bool> &polyhedronFaceNormalDirections,
        const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
        const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices,
        const Eigen::MatrixXd &referenceQuadraturePoints = (Eigen::MatrixXd(3, 1) << 1.0 / 3.0, 1.0 / 3.0, 0.0).finished(),
        const Eigen::VectorXd &referenceQuadratureWeights = Eigen::VectorXd::Constant(1, 0.5)) const;
 
    double PolyhedronVolumeByInternalIntegral(const std::vector<Eigen::MatrixXd> &polyhedronTetrahedronVertices,
                                              const Eigen::VectorXd &referenceQuadratureWeights = Eigen::VectorXd::Constant(1, 1.0 / 6.0)) const;
 
    Eigen::Vector3d PolyhedronCentroid(const std::vector<std::vector<Eigen::Matrix3d>> &polyhedronFaceRotatedTriangulationPoints,
                                       const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
                                       const std::vector<bool> &polyhedronFaceNormalDirections,
                                       const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
                                       const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices,
                                       const double &polyhedronVolume) const;
 
    inline double PolyhedronDiameter(const Eigen::MatrixXd &polyhedronVertices) const
    {
        return PointsMaxDistance(polyhedronVertices);
    }
 
    Eigen::MatrixXd PolyhedronEdgesCentroid(const Eigen::MatrixXd &polyhedronVertices, const Eigen::MatrixXi &polyhedronEdges) const;
 
    Eigen::VectorXd PolyhedronEdgesLength(const Eigen::MatrixXd &polyhedronVertices, const Eigen::MatrixXi &polyhedronEdges) const;
 
    Eigen::MatrixXd PolyhedronEdgeTangents(const Eigen::MatrixXd &polyhedronVertices, const Eigen::MatrixXi &polyhedronEdges) const;
 
    std::vector<Eigen::MatrixXd> PolyhedronFaceVertices(const Eigen::MatrixXd &polyhedronVertices,
                                                        const std::vector<Eigen::MatrixXi> &polyhedronFaces) const;
 
    std::vector<std::vector<bool>> PolyhedronFaceEdgeDirections(const Eigen::MatrixXd &polyhedronVertices,
                                                                const Eigen::MatrixXi &polyhedronEdges,
                                                                const std::vector<Eigen::MatrixXi> &polyhedronFaces) const;
 
    std::vector<Eigen::MatrixXd> PolyhedronFaceEdgeTangents(const Eigen::MatrixXd &polyhedronVertices,
                                                            const Eigen::MatrixXi &polyhedronEdges,
                                                            const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                            const std::vector<std::vector<bool>> &polyhedronFaceEdgeDirections,
                                                            const Eigen::MatrixXd &polyhedronEdgeTangents) const;
 
    std::vector<Eigen::Matrix3d> PolyhedronFaceRotationMatrices(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
                                                                const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
                                                                const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations) const;
 
    std::vector<Eigen::Vector3d> PolyhedronFaceTranslations(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices) const;
 
    std::vector<Eigen::MatrixXd> PolyhedronFaceRotatedVertices(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
                                                               const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
                                                               const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices) const;
 
    std::vector<std::vector<unsigned int>> PolyhedronFacesUnalignedVertices(const std::vector<Eigen::MatrixXd> &polyhedronFacesRotatedVertices) const;
 
    std::vector<Eigen::Vector3d> PolyhedronFaceNormals(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices) const;
 
    std::vector<std::array<Eigen::Vector3d, 2>> PolyhedronFaceTangents(const std::vector<Eigen::MatrixXd> &polyhedronFacesVertices,
                                                                       const std::vector<Eigen::Vector3d> &polyhedronFacesNormal,
                                                                       const std::vector<bool> &polyhedronFacesNormalDirection) const;
 
    std::vector<Eigen::Vector3d> PolyhedronFaceBarycenter(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices) const;
 
    Eigen::Matrix3d PolyhedronMass(const Eigen::Vector3d &polyhedronCentroid,
                                   const std::vector<Eigen::MatrixXd> &polyhedronTetrahedronPoints) const;
 
    Eigen::Matrix3d PolyhedronInertia(const Eigen::Vector3d &polyhedronCentroid,
                                      const std::vector<Eigen::MatrixXd> &polyhedronTetrahedraPoints) const;
 
    Eigen::VectorXd PolyhedronCentroidFacesDistance(const Eigen::Vector3d &polyhedronCentroid,
                                                    const std::vector<Eigen::Vector3d> &polyhedronFacesNormal,
                                                    const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices) const;
 
    inline double PolyhedronInRadius(const Eigen::VectorXd &polyhedronCentroidFacesDistance) const
    {
        return polyhedronCentroidFacesDistance.minCoeff();
    }
 
    bool PolyhedronIsConvex(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
                            const std::vector<Eigen::MatrixXd> &polyhedronFaceRotatedVertices,
                            const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints,
                            const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
                            const std::vector<bool> &polyhedronFaceNormalDirections,
                            const Eigen::Vector3d &polyhedronInternalPoint) const;
 
    std::vector<bool> PolyhedronFaceNormalDirections(const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
                                                     const Eigen::Vector3d &pointInsidePolyhedron,
                                                     const std::vector<Eigen::Vector3d> &polyhedronFaceNormals) const;
    std::vector<bool> PolyhedronFaceNormalDirections(const Eigen::MatrixXd &polyhedronVertices,
                                                     const Eigen::MatrixXi &polyhedronEdges,
                                                     const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                     const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
                                                     const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints,
                                                     const std::vector<Eigen::MatrixXd> &polyhedronFaceRotatedVertices,
                                                     const std::vector<Eigen::Vector3d> &polyhedronFaceNormals,
                                                     const std::vector<Eigen::Vector3d> &polyhedronFaceTranslations,
                                                     const std::vector<Eigen::Matrix3d> &polyhedronFaceRotationMatrices) const;
 
    std::vector<std::vector<unsigned int>> PolyhedronFaceTriangulations(const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                                        const std::vector<std::vector<unsigned int>> &localFaceTriangulations) const;
    std::vector<std::vector<unsigned int>> PolyhedronFaceTriangulationsByFirstVertex(const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                                                     const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices) const;
 
    std::vector<std::vector<unsigned int>> PolyhedronFaceTriangulationsByEarClipping(const unsigned int numPolyhedronFaces,
                                                                                     const std::vector<Eigen::MatrixXd> &polyhedronFaces2DVertices) const;
 
    std::vector<std::vector<Eigen::Matrix3d>> PolyhedronFaceExtractTriangulationPoints(
        const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
        const std::vector<std::vector<unsigned int>> &polyhedronFaceTriangulations) const;
 
    std::vector<std::vector<Eigen::Matrix3d>> PolyhedronFaceTriangulationPointsByInternalPoint(
        const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
        const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints,
        const std::vector<std::vector<unsigned int>> &polyhedronFaceTriangulations) const;
 
    std::vector<std::vector<unsigned int>> PolyhedronFaceTriangulationsByInternalPoint(
        const Eigen::MatrixXd &polyhedronVertices,
        const std::vector<Eigen::MatrixXi> &polyhedronFaces,
        const std::vector<Eigen::MatrixXd> &polyhedronFaceVertices,
        const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints) const;
 
    std::vector<unsigned int> PolyhedronTetrahedronsByFaceTriangulations(const Eigen::MatrixXd &polyhedronVertices,
                                                                         const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                                         const std::vector<std::vector<unsigned int>> &polyhedronFaceTriangulations,
                                                                         const Eigen::Vector3d &polyhedronInternalPoint) const;
 
    std::vector<unsigned int> PolyhedronTetrahedronsByFaceTriangulations(const Eigen::MatrixXd &polyhedronVertices,
                                                                         const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                                                                         const std::vector<std::vector<unsigned int>> &polyhedronFaceTriangulations,
                                                                         const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints,
                                                                         const Eigen::Vector3d &polyhedronInternalPoint) const;
 
    std::vector<Eigen::MatrixXd> ExtractTetrahedronPoints(const Eigen::MatrixXd &polyhedronVertices,
                                                          const Eigen::Vector3d &polyhedronInternalPoint,
                                                          const std::vector<unsigned int> &pointTetrahedrons) const;
    std::vector<Eigen::MatrixXd> ExtractTetrahedronPoints(const Eigen::MatrixXd &polyhedronVertices,
                                                          const Eigen::Vector3d &polyhedronInternalPoint,
                                                          const std::vector<Eigen::Vector3d> &polyhedronFaceInternalPoints,
                                                          const std::vector<unsigned int> &pointTetrahedrons) const;
 
    std::vector<unsigned int> PolyhedronCoordinateSystem(const Eigen::MatrixXd &polyhedronVertices,
                                                         const Eigen::MatrixXi &polyhedronEdges);
 
    void ExportPolyhedronToVTU(const Eigen::MatrixXd &polyhedronVertices,
                               const Eigen::MatrixXi &polyhedronEdges,
                               const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                               const std::string &exportFolder) const;
 
    inline void ExportPolyhedronToVTU(const Polyhedron &polyhedron, const std::string &exportFolder) const
    {
        ExportPolyhedronToVTU(polyhedron.Vertices, polyhedron.Edges, polyhedron.Faces, exportFolder);
    }
 
    void ExportPolyhedronToVTU(const unsigned int &index,
                               const Eigen::MatrixXd &polyhedronVertices,
                               const Eigen::MatrixXi &polyhedronEdges,
                               const std::vector<Eigen::MatrixXi> &polyhedronFaces,
                               const std::vector<Eigen::MatrixXd> &polyhedronTetra,
                               const double &polyhedronVolume,
                               const Eigen::Vector3d &polyhedronCentroid,
                               const std::vector<Eigen::MatrixXd> &polyhedronFaces3DVertices,
                               const std::vector<double> &polyhedronFacesArea,
                               const std::vector<Eigen::Vector3d> &polyhedronFaces2DCentroid,
                               const std::vector<Eigen::Vector3d> &polyhedronFacesTranslation,
                               const std::vector<Eigen::Matrix3d> &polyhedronFacesRotationMatrix,
                               const std::vector<std::vector<Eigen::Matrix3d>> &polyhedronFaces3DTriangles,
                               const std::vector<Eigen::Vector3d> &polyhedronFaces3DInternalPoint,
                               const std::vector<Eigen::Vector3d> &polyhedronFaces3DNormal,
                               const std::vector<bool> &polyhedronFaces3DNormalDirection,
                               const std::string &exportFolder) const;
 
    void ExportPolygonToVTU(const unsigned int &index,
                            const Eigen::MatrixXd &polygon,
                            const std::vector<Eigen::Matrix3d> &polygon_triangles,
                            const double &polygon_volume,
                            const Eigen::Vector3d &polygon_centroid,
                            const Eigen::MatrixXd &polygon_edges_centroid,
                            const Eigen::MatrixXd &polygon_edges_normal,
                            const std::vector<bool> &polygon_edges_normal_direction,
                            const std::string &exportFolder) const;
 
    void ExportPlanteToVTU(const Eigen::Vector3d &plane_origin,
                           const Eigen::Vector3d &plane_normal,
                           const Eigen::Vector3d &plane_translation,
                           const Eigen::Matrix3d &plane_rotation_matrix,
                           const std::string &exportFolder) const;
 
    MergePolyhedronsInput MergePolyhedronByFace(const std::array<Polyhedron, 2> &polyhedrons,
                                                const std::array<unsigned int, 2> polyhedrons_common_face_index,
                                                const bool remove_common_face) const;
 
    MergePolyhedronsResult MergePolyhedrons(const std::array<Polyhedron, 2> &polyhedrons,
                                            const MergePolyhedronsInput &merge_information = {}) const;
 
    std::map<unsigned int, unsigned int> FacetsToVertices(const std::vector<std::vector<unsigned int>> &facets) const;
    std::map<std::pair<unsigned int, unsigned int>, unsigned int> FacetsToEdges(const std::vector<std::vector<unsigned int>> &facets) const;
 
    Polyhedron FacetsToPolyhedron(const Eigen::MatrixXd &points, const std::vector<std::vector<unsigned int>> &facets) const;
    std::vector<std::vector<unsigned int>> PolyhedronToFacets(const Polyhedron &polyhedron) const;
 
    void PolygonChebyshevCenter(const Eigen::MatrixXd &polygonVertices,
                                const Eigen::MatrixXd &edgeNormals,
                                Eigen::Vector3d &in_center,
                                double &in_radius,
                                const double r0,
                                const bool use_r0) const;
 
    Eigen::MatrixXi ChangePolyhedronFaceOrientation(const Gedim::GeometryUtilities::Polyhedron &polyhedron,
                                                    const unsigned int face_index) const;
    std::vector<Eigen::MatrixXi> ChangePolyhedronFacesOrientation(const Gedim::GeometryUtilities::Polyhedron &polyhedron,
                                                                  const std::vector<bool> &faces_normal_direction) const;
};
 
} // namespace Gedim
 
#endif // __GEOMETRYUTILITIES_H