// Copyright (C) Stichting Deltares 2024. All rights reserved.
// 
// This file is part of the Dam Engine.
// 
// The Dam Engine is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// 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 Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// 
// All names, logos, and references to "Deltares" are registered trademarks of
// Stichting Deltares and remain full property of Stichting Deltares at all times.
// All rights reserved.

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.Data.Geometry;

// ----------------------------------------------------------
// The main geometry regeneration manager
public class GeometryGenerator
{
    private readonly Dictionary<GeometryCurve, bool> geometryCurveForwardsIsUsed = new ();
    private readonly Dictionary<GeometryCurve, bool> geometryCurveReversedIsUsed = new ();
    private readonly GeometryData geometryData;
    
    private readonly Dictionary<GeometryLoop, List<CurveDirection>> geometryLoopDirections = new ();
    
    private readonly List<GeometrySurface> newlyDetectedSurfaceList = new ();

    /// <summary>
    /// Regenerates the geometry.
    /// </summary>
    public GeometryGenerator(GeometryData aGeometryData)
    {
        geometryData = aGeometryData;
    }

    /// <summary>
    /// Regenerates the geometry.
    /// </summary>
    public void GenerateGeometry()
    {
        FilterOutDoublePoints();

        SetupCurveSurfaceAssociations();

        var firstRegeneration = true;
        
        while (true)
        {
            // break up all curves at intersections
            RegenerateAllCurvesIntersection();

            newlyDetectedSurfaceList.Clear();
            geometryLoopDirections.Clear();

            int curvesCount = geometryData.Curves.Count;

            // initialise IsUsed of all curves to false
            for (var index1 = 0; index1 < curvesCount; index1++)
            {
                SetIsUsed(geometryData.Curves[index1], CurveDirection.Forward, false);
                SetIsUsed(geometryData.Curves[index1], CurveDirection.Reverse, false);
            }

            // detect surfaces... the plaxis algorithm
            int result = DetectSurfaces(firstRegeneration);

            if (result < 0)
            {
                if (!firstRegeneration)
                {
                    break;
                }
            }
            else
            {
                break;
            }

            firstRegeneration = false;
        }
    }

    /// <summary>
    /// Removes the curve from both the dictionaries geometryCurveForwardsIsUsed and geometryCurveReversedIsUsed.
    /// </summary>
    /// <param name="aCurve"></param>
    public void RemoveIsUsedCurve(GeometryCurve aCurve)
    {
        if ((geometryCurveForwardsIsUsed.ContainsKey(aCurve)))
        {
            geometryCurveForwardsIsUsed.Remove(aCurve);
        }

        if ((geometryCurveReversedIsUsed.ContainsKey(aCurve)))
        {
            geometryCurveReversedIsUsed.Remove(aCurve);
        }
    }

    /// <summary>
    /// Setups the curve surface associations.
    /// </summary>
    public void SetupCurveSurfaceAssociations()
    {
        SetUpGeometryLoopDirections();
        // only try to connect curves to surfaces when there are loops (i.e. surfaces)
        if (geometryLoopDirections.Count > 0)
        {
            SetupCurveSurfaceAssociation();
        }
    }

    /// <summary>
    /// Adds the curve to the list in given direction.
    /// </summary>
    /// <param name="curve"></param>
    /// <param name="direction"></param>
    /// <param name="editedLoop"></param>
    /// <returns></returns>
    private bool AddCurve(GeometryCurve curve, CurveDirection direction, GeometryLoop editedLoop)
    {
        if (FindCurveIndex(editedLoop, curve, direction) > -1)
        {
            return false;
        }

        editedLoop.CurveList.Add(curve);
        if (!geometryLoopDirections.ContainsKey(editedLoop))
        {
            geometryLoopDirections.Add(editedLoop, new List<CurveDirection>());
        }

        geometryLoopDirections[editedLoop].Add(direction);
        return true;
    }

    // -------------------------------------------------------------------------------------
    // Loop Detection (Surface Creation) Algorithm from Plaxis

    /// <summary>
    /// Creates the surface.
    /// </summary>
    /// <param name="loop">a loop.</param>
    /// <returns></returns>
    private GeometrySurface CreateSurface(GeometryLoop loop)
    {
        if (loop.IsContinuous() && loop.HasArea())
        {
            var newSurface = new GeometrySurface();
            newSurface.SetOuterLoop(loop);
            if (!geometryData.HasSurfaceWithSameOuterLoop(newSurface.OuterLoop))
            {
                geometryData.Surfaces.Add(newSurface);
            }

            return newSurface;
        }

        return null;
    }

    /// <summary>
    /// Finds the index of the curve.
    /// </summary>
    /// <param name="geometryLoop">a geometry loop.</param>
    /// <param name="curve">a curve.</param>
    /// <param name="direction">a direction.</param>
    /// <returns></returns>
    private int FindCurveIndex(GeometryLoop geometryLoop, GeometryCurve curve, CurveDirection direction)
    {
        if (!geometryLoop.CurveList.Contains(curve))
        {
            return -1;
        }

        int curvesCount = geometryLoop.CurveList.Count;

        for (var index = 0; index < curvesCount; index++)
        {
            // Note Bka: Checking the direction does allow one curve to be added to ONE loop twice!
            // This produces some strange surfaces (see LoopDetectionCase5) but that seems to be required
            if (geometryLoop.CurveList[index] == curve && geometryLoopDirections[geometryLoop][index] == direction)
            {
                return index;
            }
        }

        return -1;
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="curve"></param>
    /// <param name="direction"></param>
    /// <param name="isUsed"></param>
    public void SetIsUsed(GeometryCurve curve, CurveDirection direction, bool isUsed)
    {
        if (direction == CurveDirection.Forward)
        {
            if (!(geometryCurveForwardsIsUsed.ContainsKey(curve)))
            {
                geometryCurveForwardsIsUsed.Add(curve, isUsed);
            }
            else
            {
                geometryCurveForwardsIsUsed[curve] = isUsed;
            }
        }
        else
        {
            if (!(geometryCurveReversedIsUsed.ContainsKey(curve)))
            {
                geometryCurveReversedIsUsed.Add(curve, isUsed);
            }
            else
            {
                geometryCurveReversedIsUsed[curve] = isUsed;
            }
        }
    }

    /// <summary>
    /// 
    /// </summary>
    /// <param name="curve"></param>
    /// <param name="direction"></param>
    /// <returns></returns>
    private bool GetIsUsed(GeometryCurve curve, CurveDirection direction)
    {
        if (direction == CurveDirection.Forward)
        {
            return geometryCurveForwardsIsUsed[curve];
        }

        return geometryCurveReversedIsUsed[curve];
    }

    /// <summary>
    /// Removes all double points, replaces them in curves and newlyeffectedpoints. 
    /// Also deletes "zero" curves (beginpoint = endpoint) also from surface loops and newlyeffectedcurrves
    /// </summary>
    private void FilterOutDoublePoints()
    {
        var doublePoints = new Dictionary<Point2D, Point2D>();

        // Make sure all points (as pointers) in curves are in the point list
        foreach (GeometryCurve curve in geometryData.Curves)
        {
            if (!geometryData.Points.Contains(curve.HeadPoint))
            {
                var point = geometryData.GetPointAtLocation(curve.HeadPoint);
                if (point != null)
                {
                    curve.HeadPoint = point;
                }
                else
                {
                    geometryData.Points.Add(curve.HeadPoint);
                }
            }

            if (!geometryData.Points.Contains(curve.EndPoint))
            {
                var point = geometryData.GetPointAtLocation(curve.EndPoint);
                if (point != null)
                {
                    curve.EndPoint = point;
                }
                else
                {
                    geometryData.Points.Add(curve.EndPoint);
                }
            }
        }

        // find double points (by location!) in point list and register them with original
        for (var i = 0; i < geometryData.Points.Count; i++)
        {
            for (var j = 0; j < i; j++)
            {
                if (i != j && geometryData.Points[i].LocationEquals(geometryData.Points[j]))
                {
                    // register the double point and the original point
                    doublePoints[geometryData.Points[i]] = geometryData.Points[j];
                    break;
                }
            }
        }

        // replace double points in curves with originals
        foreach (Point2D doublePoint in doublePoints.Keys)
        {
            foreach (GeometryCurve curve in geometryData.Curves)
            {
                if (curve.HeadPoint == doublePoint)
                {
                    curve.HeadPoint = doublePoints[doublePoint];
                }

                if (curve.EndPoint == doublePoint)
                {
                    curve.EndPoint = doublePoints[doublePoint];
                }
            }
            foreach (GeometryCurve curve in geometryData.NewlyEffectedCurves)
            {
                if (curve.HeadPoint == doublePoint)
                {
                    curve.HeadPoint = doublePoints[doublePoint];
                }

                if (curve.EndPoint == doublePoint)
                {
                    curve.EndPoint = doublePoints[doublePoint];
                }
            }
        }

        // remove curves which have the same head as end point
        foreach (GeometryCurve curve in geometryData.Curves.ToArray())
        {
            if (curve.HeadPoint == curve.EndPoint)
            {
                geometryData.Curves.Remove(curve);
                if (geometryData.NewlyEffectedCurves.Contains(curve))
                {
                    geometryData.NewlyEffectedCurves.Remove(curve);
                }

                foreach (GeometrySurface surface in geometryData.Surfaces)
                {
                    surface.OuterLoop.CurveList.Remove(curve);
                    foreach (GeometryLoop loop in surface.InnerLoops)
                    {
                        loop.CurveList.Remove(curve);
                    }
                }
            }
        }

        // removing curves from loops in surfaces may have created invalid surfaces, remove those here
        foreach (GeometrySurface surface in geometryData.Surfaces.ToArray())
        {
            if (!surface.OuterLoop.HasArea())
            {
                geometryData.Surfaces.Remove(surface);
            }
        }

        // remove double points from point list
        foreach (Point2D point in doublePoints.Keys)
        {
            geometryData.Points.Remove(point);

            if (geometryData.NewlyEffectedPoints.Contains(point))
            {
                geometryData.NewlyEffectedPoints.Remove(point);
            }
        }
    }

    private int DetectSurfaces(bool firstRegeneration)
    {
        try
        {
            // declare some variables
            int curvesCount = geometryData.Curves.Count;

            var newLoopList = new List<GeometryLoop>();
            var attachedCurveList = new List<DirectionCurve>();

            // start the first iteration
            for (var index = 0; index < curvesCount * 2; index++)
            {
                int curveIndex = index / 2;
                // look for current curve
                GeometryCurve currentCurve = geometryData.Curves[curveIndex];

                if (currentCurve == null)
                {
                    continue;
                }

                // check the direction
                CurveDirection currentCurveDirection;
                if (index % 2 == 0)
                {
                    if (GetIsUsed(currentCurve, CurveDirection.Forward))
                    {
                        continue;
                    }

                    currentCurveDirection = CurveDirection.Forward; // get the direction of the current curve
                }
                else
                {
                    if (GetIsUsed(currentCurve, CurveDirection.Reverse))
                    {
                        continue;
                    }

                    currentCurveDirection = CurveDirection.Reverse; // get the direction of the current curve
                }

                // create aValue1 new loop
                var newLoop = new GeometryLoop();

                //this.geometryData.Loops.Add(newLoop);
                newLoopList.Add(newLoop);

                // initialise LoopBeginCurve
                GeometryCurve loopBeginCurve = geometryData.Curves[curveIndex];
                CurveDirection loopBeginDirection = currentCurveDirection;

                while (true)
                {
                    // set the IsUsed status
                    SetIsUsed(currentCurve, currentCurveDirection, true);

                    // add the current curve to new loop
                    if (!AddCurve(currentCurve, currentCurveDirection, newLoop))
                    {
                        // the curve wasn't added bcos the curve-direction pair was already present in loop.
                        // problem case - break here, else we'd get aValue1 hang!
                        // Todo: Solve this problem
                        if (!firstRegeneration)
                        {
                            //TODO:Show error message box
                            break;
                        }

                        return -1;
                    }

                    Point2D curveEndPoint = currentCurve.GetEndPoint(currentCurveDirection);
                    attachedCurveList.Clear();
                    var minAngle = 365.0;
                    var minIndex = 0;

                    // find all the curves that are connected to the Current Curve at curveEndPoint
                    for (var index2 = 0; index2 < curvesCount; index2++)
                    {
                        GeometryCurve curve = geometryData.Curves[index2];
                        if (curve.LocationEquals(currentCurve)) // lets not get the reverse direction of the current curve here
                        {
                            continue;
                        }

                        if (curve.HeadPoint == curveEndPoint)
                        {
                            attachedCurveList.Add(new DirectionCurve(curve, CurveDirection.Forward));
                        }
                        else if (curve.EndPoint == curveEndPoint)
                        {
                            attachedCurveList.Add(new DirectionCurve(curve, CurveDirection.Reverse));
                        }
                    }

                    if (attachedCurveList.Count == 0) // no curves found
                    {
                        CurveDirection oppCurrentDirection = currentCurveDirection == CurveDirection.Forward ? CurveDirection.Reverse : CurveDirection.Forward;
                        // if the current curve is not used in the opposite direction, it is considered in the opposite direction
                        if (!GetIsUsed(currentCurve, oppCurrentDirection))
                        {
                            currentCurveDirection = oppCurrentDirection;
                            continue;
                        }

                        break;
                    }

                    // we have aValue1 set of curves, find the one that turns right the most
                    if (attachedCurveList.Count > 1)
                    {
                        minIndex = -1;

                        Point2D point1 = currentCurve.GetEndPoint(currentCurveDirection);
                        Point2D point2 = currentCurve.GetHeadPoint(currentCurveDirection);

                        for (var index2 = 0; index2 < attachedCurveList.Count; index2++)
                        {
                            var point3 = new Point2D();
                            var point4 = new Point2D();

                            point3.X = attachedCurveList[index2].GetHeadPoint().X;
                            point3.Z = attachedCurveList[index2].GetHeadPoint().Z;
                            point4.X = attachedCurveList[index2].GetEndPoint().X;
                            point4.Z = attachedCurveList[index2].GetEndPoint().Z;

                            double angle = Routines2D.FindAngle(point1, point2, point3, point4);

                            if (angle < minAngle)
                            {
                                minAngle = angle;
                                minIndex = index2;
                            }
                        }
                    }

                    DirectionCurve pickedDirectionCurve = attachedCurveList[minIndex];

                    if (pickedDirectionCurve.Curve == loopBeginCurve && pickedDirectionCurve.Direction == loopBeginDirection)
                    {
                        break;
                    }

                    // assign the CurrentCurve from the picked one
                    currentCurve = pickedDirectionCurve.Curve;
                    currentCurveDirection = pickedDirectionCurve.Direction;
                }
            }

            // create surfaces!
            return CreateSurfaces(newLoopList);
        }
#if DEBUG
        catch (Exception ex)
        {
            Debug.WriteLine(ex);
            return 0;
        }
#else
            catch
            {
                return 0;
            }
#endif
    }

    private int CreateSurfaces(List<GeometryLoop> newLoopList)
    {
        var newSurfacesGeoDtaObjectList = new List<IGeometryObject>();

        int loopsCount = newLoopList.Count;
        int curvesCount;
        GeometrySurface newSurface;
        var newSurfaceList = new List<GeometrySurface>();

        for (var index = 0; index < loopsCount; index++)
        {
            GeometryLoop loop = newLoopList[index];
            curvesCount = loop.CurveList.Count;

            if (curvesCount < 2) // dont create aValue1 surface for loops that have less than 2 curves
            {
                continue;
            }

            if (curvesCount == 2) // if only 2 curves in loop, make sure they are distinct (non-repeated)
            {
                if (loop.CurveList[0] == loop.CurveList[1])
                {
                    continue;
                }
            }

            if (!loop.IsContinuous() || !loop.HasArea())
            {
                continue;
            }

            // if the loop is clockwise, create aValue1 surface
            if (loop.IsClockWise() && !CheckIfLoopEnclosesOpenPolyline(loop))
            {
                // TEMP: aVector1 mechanism to remember surfaces after Geometry Regeneration
                // find the surface that is most repeated in the loop's in curves
                newSurface = GetReassignmentSurfaceFromCurves(loop);

                // an existing surface (in geometryLoopDirections) has been found
                if (newSurface != null)
                { 
                    newSurface = CreateSurface(loop);
                }
                else // no existing surface found from its comprising curves... create aValue1 brand new surface!
                {
                    newSurface = CreateSurface(loop);
                    newlyDetectedSurfaceList.Add(newSurface); // populate the newly detected surface list
                    newSurfacesGeoDtaObjectList.Add(newSurface);
                }

                if (newSurface != null)
                {
                    newSurfaceList.Add(newSurface);
                    AssignSurfaceAtLeftOrRightToCurves(newSurface);
                }
            }
        }

        // clear the left and right surfaces for all curves (some curves will have redundant data)
        curvesCount = geometryData.Curves.Count;
        for (var index = 0; index < curvesCount; index++)
        {
            geometryData.Curves[index].SurfaceAtRight = null;
            geometryData.Curves[index].SurfaceAtLeft = null;
        }

        // for the new surfaces -- assign the left/right surfaces for comprising curves, and find inner loops 
        int surfacesCount = newSurfaceList.Count;
        for (var index = 0; index < surfacesCount; index++)
        {
            newSurface = newSurfaceList[index];
            AssignSurfaceAtLeftOrRightToCurves(newSurface);
            object newSurfaceObject = newSurface /*new object()*/;
            CheckAndAddInnerLoops(ref newSurfaceObject);
            //newSurface = (GeometrySurface)newSurfaceObject;
        }

        if (newSurfacesGeoDtaObjectList.Count > 0)
        {
            var lNewSurfaces = new List<GeometrySurface>();
            GetNewlyDetectedSurfaces(ref lNewSurfaces);
        }

        return surfacesCount;
    }

    /// <summary>
    /// Checks if loop encloses open polyline.
    /// </summary>
    /// <param name="loop">a loop.</param>
    /// <returns></returns>
    private bool CheckIfLoopEnclosesOpenPolyline(GeometryLoop loop)
    {
        int curvesCount = loop.CurveList.Count;
        if (curvesCount < 3)
        {
            return true;
        }

        for (var index = 0; index < curvesCount; index++)
        {
            GeometryCurve curve = loop.CurveList[index];
            CurveDirection direction = geometryLoopDirections[loop][index];
            var foundOppDirection = false;

            for (var index1 = 0; index1 < curvesCount; index1++)
            {
                if (index == index1)
                {
                    continue;
                }

                if (loop.CurveList[index1] == curve && geometryLoopDirections[loop][index] != direction)
                {
                    foundOppDirection = true;
                    break;
                }
            }

            if (!foundOppDirection)
            {
                return false;
            }
        }

        return true;
    }

    private GeometrySurface GetReassignmentSurfaceFromCurves(GeometryLoop loop)
    {
        GeometrySurface surface = null;
        GeometrySurface reassignmentSurface = null;
        int curvesCount = loop.CurveList.Count;

        if (!geometryLoopDirections.ContainsKey(loop))
        {
            return null;
        }

        for (var index = 0; index < curvesCount; index++)
        {
            GeometryCurve curve = loop.CurveList[index];
            if (geometryLoopDirections[loop][index] == CurveDirection.Forward)
            {
                if (curve.SurfaceAtRight != null)
                {
                    surface = curve.SurfaceAtRight;
                }
            }
            else
            {
                if (curve.SurfaceAtLeft != null)
                {
                    surface = curve.SurfaceAtLeft;
                }
            }

            if (surface == null)
            {
                continue;
            }

            var maxTimesSurfacesFound = 0;
            var noTimesSurfaceFound = 0;
            for (var index1 = 0; index1 < curvesCount; index1++)
            {
                if (geometryLoopDirections[loop][index] == CurveDirection.Forward)
                {
                    if (curve.SurfaceAtRight == surface)
                    {
                        noTimesSurfaceFound++;
                    }
                }
                else
                {
                    if (curve.SurfaceAtLeft == surface)
                    {
                        noTimesSurfaceFound++;
                    }
                }
            }

            if (noTimesSurfaceFound > maxTimesSurfacesFound)
            {
                maxTimesSurfacesFound = noTimesSurfaceFound;
                reassignmentSurface = surface;
            }
        }

        return reassignmentSurface;
    }

    /// <summary>
    /// Assigns the surface at left or right to its curves on the outerloop. This tells at which side of the curve the surface is present.
    /// So after this, for each curve in the outerloop of this surface, it is known at which side the surface is located.
    /// </summary>
    /// <param name="surface">The surface.</param>
    private void AssignSurfaceAtLeftOrRightToCurves(GeometrySurface surface)
    {
        GeometryLoop loop = surface.OuterLoop;
        int curvesCount = loop.CurveList.Count;
        var isClockwise = true;

        try
        {
            isClockwise = loop.IsClockWise();
        }
        catch (GeometryLoop.NotEnoughUniquePointsException e)
        {
            Debug.WriteLine(e.Message);
        }
        catch (InvalidOperationException e)
        {
            Debug.WriteLine(e.Message);
        }

        for (var index = 0; index < curvesCount; index++)
        {
            if (isClockwise)
            {
                if (geometryLoopDirections[loop][index] == CurveDirection.Forward)
                {
                    loop.CurveList[index].SurfaceAtRight = surface;
                }
                else
                {
                    loop.CurveList[index].SurfaceAtLeft = surface;
                }
            }
            else
            {
                if (geometryLoopDirections[loop][index] == CurveDirection.Forward)
                {
                    loop.CurveList[index].SurfaceAtLeft = surface;
                }
                else
                {
                    loop.CurveList[index].SurfaceAtRight = surface;
                }
            }
        }
    }

    private void SetupCurveSurfaceAssociation()
    {
        // clear the data
        int count = geometryData.Curves.Count;
        for (var i = 0; i < count; i++)
        {
            geometryData.Curves[i].SurfaceAtLeft = null;
            geometryData.Curves[i].SurfaceAtRight = null;
        }

        // reset
        count = geometryData.Surfaces.Count;
        for (var i = 0; i < count; i++)
        {
            AssignSurfaceAtLeftOrRightToCurves(geometryData.Surfaces[i]);
        }
    }

    private void SetUpGeometryLoopDirections()
    {
        geometryLoopDirections.Clear();
        foreach (GeometryLoop loop in geometryData.Loops)
        {
            if (!geometryLoopDirections.ContainsKey(loop))
            {
                SetUpGeometryLoopDirections(loop);
            }
        }
    }

    private void SetUpGeometryLoopDirections(GeometryLoop loop)
    {
        if (loop.CurveList.Count > 0)
        {
            Point2D loopPoint;

            geometryLoopDirections.Add(loop, new List<CurveDirection>());
            // get the first curve
            if (loop.CurveList[0].EndPoint == loop.CurveList[1].HeadPoint || loop.CurveList[0].EndPoint == loop.CurveList[1].EndPoint)
            {
                geometryLoopDirections[loop].Add(CurveDirection.Forward);
                loopPoint = loop.CurveList[0].EndPoint;
            }
            else
            {
                geometryLoopDirections[loop].Add(CurveDirection.Reverse);
                loopPoint = loop.CurveList[0].HeadPoint;
            }

            // the rest of the curves
            for (var index1 = 1; index1 < loop.CurveList.Count; index1++)
            {
                if (loopPoint == loop.CurveList[index1].HeadPoint)
                {
                    geometryLoopDirections[loop].Add(CurveDirection.Forward);
                    loopPoint = loop.CurveList[index1].EndPoint;
                }
                else
                {
                    geometryLoopDirections[loop].Add(CurveDirection.Reverse);
                    loopPoint = loop.CurveList[index1].HeadPoint;
                }
            }
        }
    }

    /// <summary>
    /// Checks and adds inner loop to the new surface.
    /// </summary>
    /// <param name="newSurfaceObject"></param>
    private void CheckAndAddInnerLoops(ref object newSurfaceObject)
    {
        var newSurface = (GeometrySurface) newSurfaceObject;
        int surfaceCount = geometryData.Surfaces.Count;
        GeometryLoop newLoop = newSurface.OuterLoop;
        int newPointCount = newLoop.CurveList.Count;
        List<Point2D> newPolygon = newLoop.Points;

        newSurface.RemoveAllInnerLoops();
        for (var index = 0; index < surfaceCount; index++)
        {
            if (newSurface.Equals(geometryData.Surfaces[index]))
            {
                continue;
            }

            var innerPointCount = 0;
            var outerPointCount = 0;
            var isOnPointCount = 0;
            var hasOnPointCount = 0;
            GeometryLoop loop = geometryData.Surfaces[index].OuterLoop;
            List<Point2D> polygon = loop.Points;
            int existingLoopPointCount = polygon.Count;

            // check if it is an inner loop
            for (var innerIndex = 0; innerIndex < newPointCount; innerIndex++)
            {
                PointInPolygon location = Routines2D.CheckIfPointIsInPolygon(loop, newPolygon[innerIndex].X, newPolygon[innerIndex].Z);

                if (location == PointInPolygon.InsidePolygon)
                {
                    innerPointCount++;
                }
                else if (location == PointInPolygon.OnPolygonEdge)
                {
                    isOnPointCount++;
                }
            }

            // check if it has an inner loop
            for (var innerIndex1 = 0; innerIndex1 < existingLoopPointCount; innerIndex1++)
            {
                PointInPolygon location = Routines2D.CheckIfPointIsInPolygon(newLoop, polygon[innerIndex1].X, polygon[innerIndex1].Z);

                if (location == PointInPolygon.InsidePolygon)
                {
                    outerPointCount++;
                }
                else if (location == PointInPolygon.OnPolygonEdge)
                {
                    hasOnPointCount++;
                }
            }

            //Add New Loop as inner loop to the existing Surface
            if (innerPointCount == newPointCount)    
            {
                geometryData.Surfaces[index].AddInnerLoop(newLoop);
            }

            //Add Inner Loop to the New Surface 
            if ((outerPointCount == existingLoopPointCount) || ((outerPointCount > 0) && (existingLoopPointCount == (outerPointCount + hasOnPointCount))))
            {
                newSurface.AddInnerLoop(loop);
            }
        }
    }

    /// <summary>
    /// Gets the newly detected surface from the list.
    /// </summary>
    /// <param name="newSurfaceList"></param>
    private void GetNewlyDetectedSurfaces(ref List<GeometrySurface> newSurfaceList)
    {
        newSurfaceList.AddRange(newlyDetectedSurfaceList);
    }

    #region Nested type: DirectionCurve

    internal struct DirectionCurve
    {
        internal DirectionCurve(GeometryCurve aCurve, CurveDirection aDirection)
        {
            Curve = aCurve;
            Direction = aDirection;
        }

        internal GeometryCurve Curve { get; }

        internal CurveDirection Direction { get; }

        internal Point2D GetHeadPoint()
        {
            return Curve.GetHeadPoint(Direction);
        }

        internal Point2D GetEndPoint()
        {
            return Curve.GetEndPoint(Direction);
        }
    }

    #endregion

    /// <summary>
    /// Regenerates all the curves that have intersections.
    /// Find all intersections between curves and split them at the intersection points, adding curves where needed.
    /// Find all parallel curves and split them at the intersection pointss, adding curves where needed.
    /// At the end, clean up the points and curves.
    /// </summary>
    private void RegenerateAllCurvesIntersection()
    {
        MergePoints();
        bool isCurveInserted = false;
        bool hasCurveBeenSplit = true;
        while (hasCurveBeenSplit)
        {
            hasCurveBeenSplit = false;
            List<GeometryCurve> geometryCurveList = new List<GeometryCurve>((IEnumerable<GeometryCurve>) geometryData.NewlyEffectedCurves);
            geometryData.NewlyEffectedCurves.Clear();
            foreach (GeometryCurve geometryCurve1 in geometryCurveList)
            {
                foreach (GeometryCurve geometryCurve2 in geometryData.Curves.ToArray())
                {
                    isCurveInserted = FindAndProcessCurveIntersections(geometryCurve1, geometryCurve2, isCurveInserted, ref hasCurveBeenSplit);
                }
            }

            FilterOutDoublePoints();
        }

        geometryData.NewlyEffectedCurves.Clear();
        MergePoints();
        DeleteDuplicateCurves();
    }

    private bool FindAndProcessCurveIntersections(GeometryCurve geometryCurve1, GeometryCurve geometryCurve2, bool isCurveInserted, 
                                                 ref bool hasCurveBeenSplit)
    {
        Point2D geometryCurve1HeadPoint = geometryCurve1.HeadPoint;
        Point2D geometryCurve1EndPoint = geometryCurve1.EndPoint;
        Point2D geometryCurve2HeadPoint = geometryCurve2.HeadPoint;
        Point2D geometryCurve2EndPoint = geometryCurve2.EndPoint;
        if (!geometryCurve1.LocationEquals(geometryCurve2) && 
            !RegenerateParallelCurves(geometryCurve2, geometryCurve1, ref isCurveInserted) &&
            Routines2D.DetermineIf2DLinesIntersectStrickly(geometryCurve1HeadPoint, geometryCurve1EndPoint, geometryCurve2HeadPoint, geometryCurve2EndPoint, out Point2D intersectionPoint) ==
            LineIntersection.Intersects)
        {
            if (!Routines2D.DetermineIfPointsCoincide(geometryCurve1HeadPoint.X, geometryCurve1HeadPoint.Z, intersectionPoint.X, intersectionPoint.Z, 0.001)
                && !Routines2D.DetermineIfPointsCoincide(geometryCurve1EndPoint.X, geometryCurve1EndPoint.Z, intersectionPoint.X, intersectionPoint.Z, 0.001))
            {
                Point2D point = geometryData.CreatePointAndAddItToTheProperListsWhenNeeded(new Point2D(intersectionPoint.X, intersectionPoint.Z));
                SplitCurve(geometryCurve1, point);
                hasCurveBeenSplit = true;
            }

            if (!Routines2D.DetermineIfPointsCoincide(geometryCurve2HeadPoint.X, geometryCurve2HeadPoint.Z, intersectionPoint.X, intersectionPoint.Z, 0.001) &&
                !Routines2D.DetermineIfPointsCoincide(geometryCurve2EndPoint.X, geometryCurve2EndPoint.Z, intersectionPoint.X, intersectionPoint.Z, 0.001))
            {
                Point2D point = this.geometryData.CreatePointAndAddItToTheProperListsWhenNeeded(new Point2D(intersectionPoint.X, intersectionPoint.Z));
                SplitCurve(geometryCurve2, point);
                hasCurveBeenSplit = true;
            }
        }

        return isCurveInserted;
    }

    private void MergePoints()
    {
        int count = geometryData.Points.Count;
        List<GeometryCurve> aConnectedAtHeadCurveList = new List<GeometryCurve>();
        List<GeometryCurve> aConnectedAtEndCurveList = new List<GeometryCurve>();
        List<int> list = geometryData.NewlyEffectedPoints.Select<Point2D, int>((Func<Point2D, int>) 
                                                    (newPoint => geometryData.Points.IndexOf(newPoint))).ToList<int>();
        list.Reverse();
        foreach (int index1 in list)
        {
            if (index1 >= 0 && index1 < geometryData.Points.Count)
            {
                for (int index2 = count - 1; index2 >= 0; --index2)
                {
                    if ((index2 != index1) && (Routines2D.DetermineIfPointsCoincide(
                                                      geometryData.Points[index1].X, geometryData.Points[index1].Z, 
                                                      geometryData.Points[index2].X, geometryData.Points[index2].Z, 0.001)))
                    {
                        aConnectedAtHeadCurveList.Clear();
                        aConnectedAtEndCurveList.Clear();
                        DetermineConnectedCurves(this.geometryData.Points[index1], ref aConnectedAtHeadCurveList, ref aConnectedAtEndCurveList);
                        AssignPointToCurves(aConnectedAtHeadCurveList, geometryData.Points[index2], true);
                        AssignPointToCurves(aConnectedAtEndCurveList, geometryData.Points[index2], false);
                        --count;
                        geometryData.DeletePointAndPossibleObsoleteCurveItsPartOf(geometryData.Points[index1]);
                        break;
                    }
                    
                }
            }
        }
        geometryData.NewlyEffectedPoints.Clear();
    }

    private static double DetermineSlope(GeometryCurve line)
    {
        if (line.HeadPoint.X.IsNearEqual(line.EndPoint.X))
        {
            return double.MaxValue;
        }
        return (line.HeadPoint.Z - line.EndPoint.Z)/(line.HeadPoint.X - line.EndPoint.X);
    }

    protected internal static void DeterminePointClosestToLine(Point2D point1, Point2D point2, GeometryCurve line,
                                                              ref bool isPoint1ClosestToLine,  ref bool isPoint2ClosestToLine)
    {
        double distance1 = Routines2D.CalculateDistanceToLine(point1.X, point1.Z, 
                                                                         line.HeadPoint.X, line.HeadPoint.Z, 
                                                                         line.EndPoint.X, line.EndPoint.Z);
        double distance2 = Routines2D.CalculateDistanceToLine(point2.X, point2.Z, 
                                                                         line.HeadPoint.X, line.HeadPoint.Z, 
                                                                         line.EndPoint.X, line.EndPoint.Z);
        if (distance1 > distance2 && distance1 > GeometryConstants.Accuracy)
        {
            isPoint1ClosestToLine = false;
        }

        if (distance1 < distance2 && distance2 > GeometryConstants.Accuracy)
        {
            isPoint2ClosestToLine = false;
        }
    }

    /// <summary>
    ///  Regenerates all the Parallel curves.
    /// </summary>
    /// <param name="line1"></param>
    /// <param name="line2"></param>
    /// <param name="isCurveInserted"></param>
    /// <returns></returns>
    protected internal bool RegenerateParallelCurves(GeometryCurve line1, GeometryCurve line2, ref bool isCurveInserted)
    {
        //  check if the lines are parallel and on
        //  |----o------|----o  or  o----|------o----| or
        //  |----o------o----|  or  o----|------|----o 
        //  |o----o----------|  or  |----------o----o|
        Point2D localPoint1 = line1.HeadPoint;
        Point2D localPoint2 = line1.EndPoint;
        Point2D localPoint3 = line2.HeadPoint;
        Point2D localPoint4 = line2.EndPoint;

        bool isPoint3OnLine1 = Routines2D.DoesPointExistInLine(localPoint1, localPoint2, localPoint3, GeometryConstants.Accuracy);
        bool isPoint4OnLine1 = Routines2D.DoesPointExistInLine(localPoint1, localPoint2, localPoint4, GeometryConstants.Accuracy);
        bool isPoint1OnLine2 = Routines2D.DoesPointExistInLine(localPoint3, localPoint4, localPoint1, GeometryConstants.Accuracy);
        bool isPoint2OnLine2 = Routines2D.DoesPointExistInLine(localPoint3, localPoint4, localPoint2, GeometryConstants.Accuracy);
        bool isPoint1SameAsPoint3 = localPoint1.LocationEquals(localPoint3);
        bool isPoint1SameAsPoint4 = localPoint1.LocationEquals(localPoint4);
        bool isPoint2SameAsPoint3 = localPoint2.LocationEquals(localPoint3);
        bool isPoint2SameAsPoint4 = localPoint2.LocationEquals(localPoint4);

        if ((isPoint2SameAsPoint3 && isPoint1SameAsPoint4) || (isPoint2SameAsPoint4 && isPoint1SameAsPoint3))
        {
            geometryData.Curves.Remove(line1);
            return true;
        }
        
        CheckForNearParallelLines(line1, line2, ref isPoint1OnLine2, ref isPoint2OnLine2, ref isPoint3OnLine1, 
                                                    ref isPoint4OnLine1);

        // check for parallel
        if ((isPoint1OnLine2 && isPoint2OnLine2) || (isPoint3OnLine1 && isPoint4OnLine1)
            || ((isPoint3OnLine1 || isPoint4OnLine1) && (isPoint1OnLine2 || isPoint2OnLine2)))
        {
            if (CheckCasesWhereLine1AndLine2Overlap(line1, line2, ref isCurveInserted, isPoint3OnLine1, isPoint4OnLine1, 
                                                    isPoint1SameAsPoint3, isPoint2SameAsPoint3, isPoint2OnLine2, 
                                                    isPoint1SameAsPoint4, isPoint2SameAsPoint4))
            {
                return true;
            }

            if (CheckCasesWhereLine2IsInLine1(line1, line2, ref isCurveInserted, isPoint3OnLine1, isPoint4OnLine1, 
                                              isPoint1SameAsPoint3, isPoint2SameAsPoint4, isPoint2SameAsPoint3, 
                                              isPoint1SameAsPoint4))
            {
                return true;
            }

            if (CheckCasesWhereLine1IsInLine2(line1, line2, ref isCurveInserted, isPoint1OnLine2, isPoint2OnLine2, 
                                              isPoint1SameAsPoint3, isPoint2SameAsPoint4, isPoint2SameAsPoint3, 
                                              isPoint1SameAsPoint4))
            {
                return true;
            }

            if (CheckCasesWhereLine2IsInLine1WithOnePointEqual(line1, line2, isPoint4OnLine1, isPoint3OnLine1, 
                                                               isPoint1SameAsPoint3, isPoint2SameAsPoint3, isPoint2SameAsPoint4))
            {
                return true;
            }

            if (CheckCasesWhereLine1IsInLine2WithOnePointEqual(line1, line2, isPoint1OnLine2, isPoint2OnLine2, 
                                                               isPoint1SameAsPoint3, isPoint1SameAsPoint4, isPoint2SameAsPoint4))
            {
                return true;
            }
        }

        return false;
    }

    private bool CheckCasesWhereLine1IsInLine2(GeometryCurve line1, GeometryCurve line2, ref bool isCurveInserted, 
                                               bool isPoint1OnLine2, bool isPoint2OnLine2, bool isPoint1SameAsPoint3, 
                                               bool isPoint2SameAsPoint4, bool isPoint2SameAsPoint3, bool isPoint1SameAsPoint4)
    {
        GeometryCurve newLine;
        // p3----p1------p2----p4 
        if (isPoint1OnLine2 && isPoint2OnLine2 && !isPoint1SameAsPoint3 && !isPoint2SameAsPoint4 &&
            !isPoint2SameAsPoint3 && !isPoint1SameAsPoint4)
        {
            double distance1 = Routines2D.Compute2DDistance(line2.HeadPoint.X, line2.HeadPoint.Z, line1.HeadPoint.X, line1.HeadPoint.Z);
            double distance2 = Routines2D.Compute2DDistance(line2.HeadPoint.X, line2.HeadPoint.Z, line1.EndPoint.X, line1.EndPoint.Z);

            if (distance1 < distance2)
            {
                newLine = SplitCurve(line2, line1.HeadPoint);
                newLine.HeadPoint = line1.EndPoint;
            }
            else
            {
                newLine = SplitCurve(line2, line1.EndPoint);
                newLine.HeadPoint = line1.HeadPoint;
            }
                
            isCurveInserted = true;
            return true;
        }

        return false;
    }

    private bool CheckCasesWhereLine2IsInLine1(GeometryCurve line1, GeometryCurve line2, ref bool isCurveInserted, 
                                               bool isPoint3OnLine1, bool isPoint4OnLine1, bool isPoint1SameAsPoint3, 
                                               bool isPoint2SameAsPoint4, bool isPoint2SameAsPoint3, bool isPoint1SameAsPoint4)
    {
        GeometryCurve newLine;
        //  p1----p3------p4----p2
        if (isPoint3OnLine1 && isPoint4OnLine1 && !isPoint1SameAsPoint3 && !isPoint2SameAsPoint4 &&
            !isPoint2SameAsPoint3 && !isPoint1SameAsPoint4)
        {
            double distance1 = Routines2D.Compute2DDistance(line1.HeadPoint.X, line1.HeadPoint.Z, line2.HeadPoint.X, 
                                                            line2.HeadPoint.Z);
            double distance2 = Routines2D.Compute2DDistance(line1.HeadPoint.X, line1.HeadPoint.Z, line2.EndPoint.X, 
                                                            line2.EndPoint.Z);

            if (distance1 < distance2)
            {
                newLine = SplitCurve(line1, line2.HeadPoint);
                newLine.HeadPoint = line2.EndPoint;
            }
            else
            {
                newLine = SplitCurve(line1, line2.EndPoint);
                newLine.HeadPoint = line2.HeadPoint;
            }
                
            isCurveInserted = true;
            return true;
        }

        return false;
    }

    private bool CheckCasesWhereLine1AndLine2Overlap(GeometryCurve line1, GeometryCurve line2, ref bool isCurveInserted, 
                                                     bool isPoint3OnLine1, bool isPoint4OnLine1, bool isPoint1SameAsPoint3, 
                                                     bool isPoint2SameAsPoint3, bool isPoint2OnLine2, bool isPoint1SameAsPoint4, 
                                                     bool isPoint2SameAsPoint4)
    {
        GeometryCurve newLine;
        //  p1----p3------p2----p4  
        if (isPoint3OnLine1 && !isPoint4OnLine1 && !isPoint1SameAsPoint3 && !isPoint2SameAsPoint3)
        {
            newLine = SplitCurve(line1, line2.HeadPoint);

            line2.HeadPoint = isPoint2OnLine2 ? newLine.EndPoint : line1.HeadPoint;
                
            isCurveInserted = true;
            return true;
        }

        //  p1----p4------p2----p3  
        if (isPoint4OnLine1 && !isPoint3OnLine1 && !isPoint1SameAsPoint4 && !isPoint2SameAsPoint4)
        {
            newLine = SplitCurve(line1, line2.EndPoint);
            line2.EndPoint = isPoint2OnLine2 ? newLine.EndPoint : line1.HeadPoint;
                
            isCurveInserted = true;
            return true;
        }

        return false;
    }

    private static bool CheckCasesWhereLine2IsInLine1WithOnePointEqual(GeometryCurve line1, GeometryCurve line2, 
                                                                       bool isPoint4OnLine1, bool isPoint3OnLine1, 
                                                                       bool isPoint1SameAsPoint3, bool isPoint2SameAsPoint3,
                                                                       bool isPoint2SameAsPoint4)
    {
        if (isPoint4OnLine1 && isPoint3OnLine1)
        {
            // Line2 on Line1  like p1--p3--p4--p2, p2--p3--p4--p1 etc.
            if (isPoint1SameAsPoint3)
            {
                // p1/p3----p4----p2  --> c1(p4, p2) & c2(p3, p4)
                line1.HeadPoint = line2.EndPoint;
            }
            else if (isPoint2SameAsPoint3)
            {
                // p2/p3----p4----p1  --> c1(p1, p4) & c2(p3, p4)  
                line1.EndPoint = line2.EndPoint;
            }
            else if (isPoint2SameAsPoint4)
            {
                // p2/p4----p3----p1  --> c1(p1, p3) & c2(p3, p4)
                line1.EndPoint = line2.HeadPoint;
            }
            else
            {
                // p1/p4----p3----p2  --> c1(p3, p2) & c2(p3, p4)
                line1.HeadPoint = line2.HeadPoint;
            }

            return true;
        }

        return false;
    }

    private static bool CheckCasesWhereLine1IsInLine2WithOnePointEqual(GeometryCurve line1, GeometryCurve line2, 
                                                                       bool isPoint1OnLine2, bool isPoint2OnLine2, 
                                                                       bool isPoint1SameAsPoint3, bool isPoint1SameAsPoint4, 
                                                                       bool isPoint2SameAsPoint4)
    {
        // Line1 on Line2  like p3--p1--p2--p4, p4--p1--p2--p3 etc.   
        if (isPoint1OnLine2 && isPoint2OnLine2)
        {
            if (isPoint1SameAsPoint3)
            {
                //p3/p1--p2--p4
                line2.HeadPoint = line1.EndPoint;
            }
            else if (isPoint1SameAsPoint4)
            {
                //p4/p1--p2--p3
                line2.EndPoint = line1.EndPoint;
            }
            else if (isPoint2SameAsPoint4)
            {
                //p4/p2--p1--p3
                line2.EndPoint = line1.HeadPoint;
            }
            else
            {
                //p3/p2--p1--p4
                line2.HeadPoint = line1.HeadPoint;
            }

            return true;
        }

        return false;
    }

    private static void CheckForNearParallelLines(GeometryCurve line1, GeometryCurve line2, ref bool isPoint1OnLine2,  
                                                  ref bool isPoint2OnLine2, ref bool isPoint3OnLine1, ref bool isPoint4OnLine1)
    {
        double slope1 = DetermineSlope(line1);
        double slope2 = DetermineSlope(line2);
        const double accuracyAngle = 5E-12;
        if (Math.Abs(Math.Abs(slope1) - Math.Abs(slope2)) > accuracyAngle)
        {
            if (isPoint1OnLine2 && isPoint2OnLine2)
            {
                // lines are not (exactly) parallel but within the geometry accuracy off each other
                DeterminePointClosestToLine(line1.HeadPoint, line1.EndPoint, line2, ref isPoint1OnLine2, ref isPoint2OnLine2);
            }
            if (isPoint3OnLine1 && isPoint4OnLine1)
            {
                // lines are not (exactly) parallel but within the geometry accuracy off each other
                DeterminePointClosestToLine(line2.HeadPoint, line2.EndPoint, line1, ref isPoint3OnLine1, ref isPoint4OnLine1);
            }
        }
    }

    private GeometryCurve SplitCurve(GeometryCurve srcCurve, Point2D pointOnCurve)
    {
        if (pointOnCurve.LocationEquals(srcCurve.HeadPoint) || pointOnCurve.LocationEquals(srcCurve.EndPoint))
            return null;
        GeometryCurve newCurve = geometryData.CreateCurve(pointOnCurve, srcCurve.EndPoint);
        srcCurve.EndPoint = newCurve.HeadPoint;
        newCurve.AssignSurfacesFromCurve(srcCurve);
        if (!geometryData.DoesCurveExist(newCurve))
        {
            geometryData.Curves.Add(newCurve);
        }
        return newCurve;
    }

    private void DetermineConnectedCurves(
        Point2D point,
        ref List<GeometryCurve> connectedAtHeadCurveList,
        ref List<GeometryCurve> connectedAtEndCurveList)
    {
        connectedAtHeadCurveList.Clear();
        connectedAtEndCurveList.Clear();
        int count = this.geometryData.Curves.Count;
        for (int index = 0; index < count; ++index)
        {
            GeometryCurve curve = this.geometryData.Curves[index];
            if (curve.HeadPoint == point)
                connectedAtHeadCurveList.Add(curve);
            if (curve.EndPoint == point)
                connectedAtEndCurveList.Add(curve);
        }
    }

    private static void AssignPointToCurves(List<GeometryCurve> aCurveList, Point2D aPoint, bool aHead)
    {
        int count = aCurveList.Count;
        if (count < 1)
            return;
        for (int i = 0; i < count; ++i)
        {
            if (aHead)
                aCurveList[i].HeadPoint = aPoint;
            else
                aCurveList[i].EndPoint = aPoint;
        }
    }
    
    private void DeleteDuplicateCurves()
    {
        List<GeometryCurve> curvesToDelete = geometryData.Curves.Where((Func<GeometryCurve, bool>) (curve => curve.HeadPoint == curve.EndPoint)).ToList();
        foreach (GeometryCurve aCurve in curvesToDelete)
            geometryData.DeleteCurve(aCurve, true);
        curvesToDelete.Clear();
        for (int i = 0; i < geometryData.Curves.Count; ++i)
        {
            GeometryCurve curve1 = geometryData.Curves[i];
            for (int j = i + 1; j < geometryData.Curves.Count; ++j)
            {
                GeometryCurve curve2 = geometryData.Curves[j];
                if ((curve1.HeadPoint == curve2.HeadPoint && curve1.EndPoint == curve2.EndPoint) || 
                  (curve1.HeadPoint == curve2.EndPoint && curve1.EndPoint == curve2.HeadPoint))
                {
                    curvesToDelete.Add(curve2);
                }
            }
        }
        foreach (GeometryCurve aCurve in curvesToDelete)
            geometryData.DeleteCurve(aCurve, true);
    }
    
}