// Copyright (C) Stichting Deltares 2023. 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.Linq;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.Data.Geometry
{
    /// <summary>
    /// Classifications of an objected with regards to a geometric line in the XZ-plane.
    /// </summary>
    public enum RelativeXzPosition
    {
        /// <summary>
        /// Indicates that the object is considered 'on' the geometric line.
        /// </summary>
        OnGeometricLine,
        /// <summary>
        /// Indicates that the object is considered 'above' the geometric line (object Z
        /// coordinate higher than that of the line).
        /// </summary>
        AboveGeometricLine,
        /// <summary>
        /// Indicates that the object is considered 'below' the geometric line (object Z
        /// coordinates lower than that of the line).
        /// </summary>
        BelowGeometricLine,
        /// <summary>
        /// Indicates that the object is considered 'outside the scope' of the geometric line.
        /// </summary>
        BeyondGeometricLine
    }

    /// <summary>
    /// Type of extrapolation
    /// </summary>
    public enum ExtraPolationMode
    {
        /// <summary>
        /// No extrapolation should be used.
        /// </summary>
        Beyond,
        /// <summary>
        /// Used constant (0th order) extrapolation at the extremes.
        /// </summary>
        Horizontal
    }

    /// <summary>
    /// Collection of points (X,Z) in world coordinates used for describing Surface lines
    /// </summary>
    public class GeometryPointString : GeometryObject
    {
        /// <summary>
        /// Matching distance where a point within this range is considered on the same point.
        /// </summary>
        private const double epsilon = GeometryConstants.Accuracy;
        private readonly List<GeometryPoint> points = new List<GeometryPoint>();
        /// <summary>
        /// The calculate points as protected field (to be able to prevent recursive calls to CalcPoints)
        /// </summary>
        protected readonly List<Point2D> calcPoints = new List<Point2D>();
        private bool isFrozen;
        private bool hasNaNx;
        private double frozenMaxZ = double.NaN;

        // sortedPoints must never be used outside this class. Either the GPS concerned must have sorted points but then they already are 
        // (eg. surfaceline, headline) or they may be unsorted in which case using the sorted list in other classes leads to errors (eg.
        // geometrysurfaces, waternetlines) 
        private List<Point2D> sortedPoints;

        /// <summary>
        /// The calculate points (to be used in calcualtion instead of Points for better performance)
        /// </summary>
        public virtual List<Point2D> CalcPoints
        {
            get
            {
                return calcPoints;
            }
        }
        
        /// <summary>
        /// Freezes this instance.
        /// </summary>
        public void Freeze()
        {
            if (!isFrozen)
            {
                sortedPoints = new List<Point2D>(calcPoints.Count);
                foreach (Point2D point in calcPoints)
                {
                    sortedPoints.Add(point);
                    hasNaNx = hasNaNx || double.IsNaN(point.X);
                    frozenMaxZ = Math.Max(frozenMaxZ, point.Z);
                }
                sortedPoints.Sort();
            }
            isFrozen = true;
        }

        /// <summary>
        /// List of points that describe the physical surface line or surface.
        /// </summary>
        /// <value>
        /// The points.
        /// </value>
        public virtual IList<GeometryPoint> Points
        {
            get
            {
                return points;
            }
        }

        /// <summary>
        /// Gets the <see cref="GeometryPoint"/> at the specified index.
        /// </summary>
        /// <value>
        /// The <see cref="GeometryPoint"/>.
        /// </value>
        /// <param name="index">The index.</param>
        /// <exception cref="ArgumentOutOfRangeException">When <paramref name="index"/> less
        /// than zero or is greater or equals to <see cref="Count"/>.</exception>
        public Point2D this[int index]
        {
            get
            {
                return calcPoints[index];
            }
        }

        /// <summary>
        /// Gets the count of the points.
        /// </summary>
        /// <value>
        /// The count.
        /// </value>
        public int Count
        {
            get
            {
                return calcPoints.Count;
            }
        }


        /// <summary>
        /// Clones this object except the points
        /// </summary>
        /// <returns></returns>
        public virtual GeometryPointString Clone()
        {
            var clone = new GeometryPointString();
            this.CloneProperties(clone); // exludes the points !
            clone.Points.Clear();
            foreach (var point in Points)
            {
                clone.Points.Add(new GeometryPoint(point));
            }
            clone.SyncCalcPoints();
            return clone;
        }

        /// <summary>
        /// Gets the minimum z.
        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
        /// </summary>
        /// <returns></returns>
        public double GetMinZ()
        {
            return CalcPoints.Any()
                       ? CalcPoints.Min(p => p.Z)
                       : double.NaN;
        }

        /// <summary>
        /// Gets the maximum z.
        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
        /// </summary>
        /// <returns></returns>
        public double GetMaxZ()
        {
            if (isFrozen) return frozenMaxZ;
            return CalcPoints.Any()
                       ? CalcPoints.Max(p => p.Z)
                       : double.NaN;
        }

        /// <summary>
        /// The minimal X value among <see cref="CalcPoints"/>.
        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
        /// </summary>
        /// <returns>The minimal X value or <see cref="double.NaN"/> in case there are no points.</returns>
        public double GetMinX()
        {
            if (isFrozen && !hasNaNx) return sortedPoints[0].X;
            return CalcPoints.Any(p => !double.IsNaN(p.X))
                       ? CalcPoints.Where(p => !double.IsNaN(p.X)).Min(p => p.X)
                       : double.NaN;
        }

        /// <summary>
        /// Gets the maximum x.
        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
        /// </summary>
        /// <returns></returns>
        public double GetMaxX()
        {
            if (isFrozen && !hasNaNx) return sortedPoints[sortedPoints.Count - 1].X;
            return CalcPoints.Any(p => !double.IsNaN(p.X))
                       ? CalcPoints.Where(p => !double.IsNaN(p.X)).Max(p => p.X)
                       : double.NaN;
        }

        /// <summary>
        /// Finds all vertical XZ intersections with this geometric point string, returning
        /// the highest value among the intersections.
        /// </summary>
        /// <param name="x">X coordinate</param>
        /// <returns>
        /// The z value determined by <see cref="CalcPoints" /> or <see cref="double.NaN" />
        /// when <see cref="CalcPoints" /> is empty.
        /// </returns>
        /// <remarks>
        /// Uses constant extrapolation and linear interpolation.
        /// </remarks>        
        public double GetZAtUnsortedX(double x)
        {
            if (calcPoints.Count > 0)
            {
                var verticalLineAtX = new Line(new Point2D {X = x, Z= GetMaxZ()},
                                               new Point2D { X = x, Z = GetMinZ() });
                if (Math.Abs(verticalLineAtX.BeginPoint.Z - verticalLineAtX.EndPoint.Z) < GeometryConstants.Accuracy)
                {
                    verticalLineAtX.BeginPoint.Z += 1.0;
                    verticalLineAtX.EndPoint.Z -= 1.0;
                }
                var intersectionPoints = IntersectionPointsXzWithLineXz(verticalLineAtX);
                if (intersectionPoints.Count != 0)
                {
                    return intersectionPoints.Max(gp => gp.Z);
                }

                // Remain consistent with GetZAtX, do constant extrapolation:
                double zAtX;
                if (DoConstantExtrapolationXz(x, out zAtX))
                {
                    return zAtX;
                }
            }
            return double.NaN;
        }

        /// <summary>
        /// Retrieves the Z value for the given x.
        /// </summary>
        /// <param name="x">X coordinate</param>
        /// <returns>
        /// The z value determined by <see cref="CalcPoints"/> or <see cref="double.NaN"/>
        /// when <see cref="CalcPoints"/> is empty.
        /// </returns>
        /// <remarks>
        /// <para>Uses constant extrapolation and linear interpolation.</para>
        /// </remarks>
        public virtual double GetZatX(double x)
        {
            if (calcPoints.Any())
            {
                double zAtX;
                if (DoConstantExtrapolationXz(x, out zAtX))
                {
                    return zAtX;
                }
            }
            for (int i = 0; i < calcPoints.Count - 1; i++)
            {
                var current = calcPoints[i];
                var next = calcPoints[i + 1];

                var leftOffset = x - current.X;
                var rightOffset = next.X - x;

                if (Math.Abs(leftOffset) < epsilon)
                {
                    return current.Z;
                }
                if (Math.Abs(rightOffset) < epsilon)
                {
                    return next.Z;
                }
                if (leftOffset >= 0 && rightOffset >= 0)
                {
                    var fraction = leftOffset / (leftOffset + rightOffset);

                    return (1.0 - fraction) * current.Z + fraction * next.Z;
                }
            }

            return Double.NaN;
        }

        /// <summary>
        /// Gets the z at x starting from index i in the point list.
        /// This is only meant as a fast(er) version for GetZAtX for lists that are sorted by X, so use with care.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="i">The i.</param>
        /// <returns></returns>
        public virtual double GetZatX(double x, ref int i)
        {
            for (; i < calcPoints.Count - 1; i++)
            {
                var current = calcPoints[i];
                var leftOffset = x - current.X;
                var next = calcPoints[i + 1];
                var rightOffset = next.X - x;

                if (leftOffset < epsilon)
                {
                    return current.Z;
                }
                if (rightOffset >= epsilon)
                {
                    var fraction = leftOffset / (leftOffset + rightOffset);
                    return (1.0 - fraction) * current.Z + fraction * next.Z;
                }
                if (i + 1 == calcPoints.Count - 1)
                {
                    return next.Z;
                }
            }

            return Double.NaN;
        }

        /// <summary>
        /// Gets all z-values at x for a line.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <returns></returns>
        public List<double> GetAllZatXForLine(double x)
        {
            return GetAllZatX(x, false);
        }

        /// <summary>
        /// Gets all z values at x for a surface.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <returns></returns>
        public List<double> GetAllZatXForSurface(double x)
        {
            return GetAllZatX(x, true);
        }

        /// <summary>
        /// Finds the first intersection of the line with a given horizontal line.
        /// </summary>
        /// <param name="z">The height level of the horizontal line.</param>
        /// <param name="xmin">Optional: Discard intersections whose X value is less then 
        /// or equal to this value.</param>
        /// <returns>The first intersection point matching the criteria.</returns>
        /// <seealso cref="GetXatZStartingAt"/>
        public double GetXatZ(double z, double xmin = Double.MinValue)
        {
            return GetXatZStartingAt(0, z, xmin);
        }

        /// <summary>
        /// Finds the first intersection of the line with a given horizontal line.
        /// </summary>
        /// <param name="start">Start considering <see cref="GeometryPoint"/> instances from
        /// <see cref="CalcPoints"/> starting from this index.</param>
        /// <param name="z">The height level of the horizontal line.</param>
        /// <param name="xmin">Optional: Discard intersections whose X value is less then 
        /// or equal to this value.</param>
        /// <returns>The first intersection point matching the criteria.</returns>
        /// <seealso cref="GetXatZ"/>
        public double GetXatZStartingAt(int start, double z, double xmin = Double.MinValue)
        {
            for (int i = start; i < calcPoints.Count - 1; i++)
            {
                var current = calcPoints[i];
                var next = calcPoints[i + 1];

                if (IsSegmentNotCrossingZ(z, current, next)) // Performance micro-optimization
                {
                    continue;
                }

                var x = GetXIntersectingZ(z, current, next);
                if (x > xmin)
                {
                    return x;
                }
            }

            return Double.NaN;
        }

        /// <summary>
        /// Gets the points at x.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <returns></returns>
        public IEnumerable<Point2D> GetPointsAtX(double x)
        {
            return (from Point2D point in calcPoints
                    where point.X.AlmostEquals(x, GeometryPoint.Precision)
                    select point);
        }

        /// <summary>
        /// Gets the Point2D at the specified coordinates.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="z">The z.</param>
        /// <returns></returns>
        public Point2D GetPointAt(double x, double z)
        {
            return calcPoints.FirstOrDefault(
                point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
                         point.Z.AlmostEquals(z, GeometryPoint.Precision));
        }

        /// <summary>
        /// Gets the GeometryPoint at the specified coordinates.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="z">The z.</param>
        /// <returns></returns>
        public GeometryPoint GetGeometryPointAt(double x, double z)
        {
            return Points.FirstOrDefault(
                point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
                         point.Z.AlmostEquals(z, GeometryPoint.Precision));
        }

        //#Bka: this should never be 1 method! Split it into 2 methods
        // It is used for sliplanes only and therefor uses Points instead of calcPoints
        /// <summary>
        /// Adds the or remove list object.
        /// </summary>
        /// <param name="aPoint">a point.</param>
        /// <param name="aToAdd">if set to <c>true</c> [a to add].</param>
        /// <param name="aIndex">a index.</param>
        public virtual void AddOrRemoveListObject(object aPoint, bool aToAdd, int aIndex)
        {
            var point = aPoint as GeometryPoint;
            if (point == null)
            {
                return;
            }

            if (aToAdd)
            {
                Points.Insert(aIndex, point); //#bka what happens if index already exists? what if point exists? What if both???                
            }
            else
            {
                if (Points.Contains(point))
                {
                    Points.Remove(point);                    
                }
            }
        }

        /// <summary>
        /// Returns ALL intersection points that are found for a line costructed at level z, 
        /// including the double ones as the number of points can be relevant!
        /// </summary>
        /// <param name="z">the level at which the line for intersections is constructed.</param>
        /// <returns></returns>
        public IList<double> IntersectionsXAtZ(double z)
        {
            var intersectionsX = new List<double>();

            var lineAtZ = new Line
            {
                BeginPoint = new Point2D(GetMinX(), z),
                EndPoint = new Point2D(GetMaxX(), z)
            };

            intersectionsX.AddRange(IntersectionPointsXzWithLineXzWithAllPoints(lineAtZ).Select(gp => gp.X));
            return intersectionsX;
        }

        /// <summary>
        /// Returns ALL intersection points that are found, including the double ones as the number of points can be relevant!
        /// </summary>
        /// <param name="line"></param>
        /// <returns></returns>
        public IList<Point2D> IntersectionPointsXzWithLineXzWithAllPoints(Line line)
        {
            return IntersectionPointsWithLineCore(line, true);
        }

        /// <summary>
        /// Returns the UNIQUE intersectionpoints (so can not be used where number of interscections is relevant)
        /// </summary>
        /// <param name="line"></param>
        /// <returns></returns>
        public IList<Point2D> IntersectionPointsXzWithLineXz(Line line)
        {
            return IntersectionPointsWithLineCore(line, false);
        }

        private IList<Point2D> IntersectionPointsXzClosedStringWithLineXz(Line line)
        {
            IList<Point2D> intersectionPointsWithLine = IntersectionPointsXzWithLineXz(line);

            // close the poly line 
            if (calcPoints.Count > 0)
            {
                DoIntersectAndAddToCollection(line, calcPoints[calcPoints.Count - 1], calcPoints[0], 
                    intersectionPointsWithLine, false);
            }

            return intersectionPointsWithLine;
        }

        /// <summary>
        /// Find intersection (xz-plane) from this surface with another surface
        /// or phratic line 
        /// </summary>
        /// <param name="list"></param>
        /// <returns></returns>
        /// <remarks>Considers all <see cref="GeometryPoint"/> in <see cref="CalcPoints"/> to
        /// such that they describe a closed loop.</remarks>
        /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
        public IList<Point2D> ClosedGeometryIntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
        {
            return IntersectWithPointsListCore(list, true);
        }

        /// <summary>
        /// Finds all intersections in the XZ-plane the given <see cref="GeometryPoint"/> list.
        /// </summary>
        /// <param name="list">The list.</param>
        /// <returns></returns>
        /// <seealso cref="ClosedGeometryIntersectionXzPointsWithGeometryPointList"/>
        /// <seealso cref="IntersectionXzPointsWithGeometryString"/>
        private List<Point2D> IntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
        {
            return IntersectWithPointsListCore(list, false);
        }

        /// <summary>
        /// Finds all intersections in the XZ-plane the given <see cref="GeometryPointString"/>.
        /// </summary>
        /// <param name="externalSurface">The geometry point string.</param>
        /// <returns></returns>
        /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
        public List<Point2D> IntersectionXzPointsWithGeometryString(GeometryPointString externalSurface)
        {
            return IntersectionXzPointsWithGeometryPointList(externalSurface.CalcPoints);
        }

        /// <summary>
        /// Sorts the points by x ascending (only to be used for Surface lines).
        /// </summary>
        public virtual void SortPointsByXAscending()
        {
            calcPoints.Sort();
            points.Sort();
        }

        /// <summary>
        /// Removes the non ascending points on surface line.
        /// </summary>
        public void RemoveNonAscendingPointsOnSurfaceLine()
        {
            for (int i = points.Count - 1; i > 0; i--)
            {
                if (Math.Abs(points[i].X - points[i - 1].X) < GeometryConstants.Accuracy)
                {
                    if (Math.Abs(points[i].Z - points[i - 1].Z) > GeometryConstants.Accuracy)
                    {
                        points.Remove(points[i]);
                    }
                }
            }
        }

        /// <summary>
        /// Gets a part of the geometry point string defined by a begin and end point
        /// </summary>
        /// <param name="begin"></param>
        /// <param name="end"></param>
        /// <returns></returns>
        public GeometryPointString GetPart(double begin, double end)
        {
            var part = new GeometryPointString();
            bool filling = false;
            bool filled = false;

            for (int i = 0; i < calcPoints.Count; i++)
            {
                if (!filling && !filled)
                {
                    filling = calcPoints[i].X >= begin - epsilon;
                }
                else
                {
                    filled = calcPoints[i].X >= end - epsilon;
                    if (filled)
                    {
                        filling = false;
                    }
                }

                if (filling)
                {
                    part.calcPoints.Add(calcPoints[i]);
                }
            }

            var beginPoint = new Point2D(begin, GetZatX(begin));
            var endPoint = new Point2D(end, GetZatX(end));

            if (part.calcPoints.Count == 0)
            {
                part.calcPoints.Add(beginPoint);
                part.calcPoints.Add(endPoint);
            }

            if (!part.calcPoints[0].LocationEquals(beginPoint))
            {
                part.calcPoints.Insert(0, beginPoint);
            }

            if (!part.calcPoints[part.calcPoints.Count - 1].LocationEquals(endPoint))
            {
                part.calcPoints.Add(endPoint);
            }
            SyncPoints();
            return part;
        }

        /// <summary>
        /// Removes all double points at a location, if consecutive
        /// </summary>
        public void CondensePoints()
        {
            for (int i = calcPoints.Count - 1; i > 0; i--)
            {
                if (calcPoints[i].LocationEquals(calcPoints[i - 1]))
                {
                    calcPoints.RemoveAt(i);
                }
            }
        }

        /// <summary>
        /// Removes all points which don't influence the shape of the line
        /// </summary>
        public void RemoveUnnecessaryPoints()
        {
            const double slopeTolerance = 1E-10;

            CondensePoints();

            for (int i = calcPoints.Count - 2; i > 1; i--)
            {
                // if the slope of the line before the point is equal to the slope after the point, the point can be removed
                double slopeBefore = (calcPoints[i].Z - calcPoints[i - 1].Z) / (calcPoints[i].X - calcPoints[i - 1].X);
                double slopeAfter = (calcPoints[i + 1].Z - calcPoints[i].Z) / (calcPoints[i + 1].X - calcPoints[i].X);

                if (Routines2D.AreEqual(slopeBefore, slopeAfter, slopeTolerance))
                {
                    calcPoints.RemoveAt(i);
                }
            }
        }

        /// <summary>
        /// Gets the surrounding rectangle around the geometry point string
        /// </summary>
        /// <returns></returns>
        public override GeometryBounds GetGeometryBounds()
        {
            if (!Points.Any())
            {
                // Sync with calcPoints
                SyncPoints();
                // if still no points, then return null
                if (!Points.Any()) return null;
            }

            GeometryBounds bounds = Points[0].GetGeometryBounds();
            for (var i = 1; i < Points.Count; i++)
            {
                GeometryPoint point = Points[i];

                bounds.Left = Math.Min(bounds.Left, point.X);
                bounds.Right = Math.Max(bounds.Right, point.X);
                bounds.Top = Math.Max(bounds.Top, point.Z);
                bounds.Bottom = Math.Min(bounds.Bottom, point.Z);
            }

            return bounds;
        }

        /// <summary>
        /// Checks if constant extrapolation can be applied, and if so set the Z value.
        /// </summary>
        /// <param name="x">The evaluated X coordinate.</param>
        /// <param name="z">Output param: Extrapolated Z value, or <see cref="double.NaN"/>
        /// when no extrapolation is possible.</param>
        /// <returns>True if extrapolation possible; false otherwise.</returns>
        private bool DoConstantExtrapolationXz(double x, out double z)
        {
            if (calcPoints.Count > 0)
            {
                var first = calcPoints[0];
                if (x < first.X || Math.Abs(x - first.X) < epsilon)
                {
                    z = first.Z;
                    return true;
                }

                var last = calcPoints[calcPoints.Count - 1];
                if (x > last.X)
                {
                    z = last.Z;
                    return true;
                }
            }

            z = double.NaN;
            return false;
        }

        /// <summary>
        /// Can be used for a Line or for a Surface where a surface is supposed to closed.
        /// In case of a waternet line it is possible
        /// to have more z values at a give x coor
        /// Furthermore a z is not needed at al x values
        /// </summary>
        /// <param name="x"></param>
        /// <param name="asSurface"></param>
        /// <returns></returns>
        private List<double> GetAllZatX(double x, bool asSurface)
        {
            var result = new List<double>();

            if (calcPoints.Count == 1 && Math.Abs(calcPoints[0].X - x) < epsilon)
            {
                result.Add(calcPoints[0].Z);
            }

            var pointsCount = calcPoints.Count - 1;
            if (asSurface)
            {
                pointsCount = calcPoints.Count;
            }
            for (int i = 0; i < pointsCount; i++)
            {
                Point2D current;
                Point2D next;
                if (i == calcPoints.Count - 1)
                {
                    current = calcPoints[i];
                    next = calcPoints[0];
                }
                else
                {
                    current = calcPoints[i];
                    next = calcPoints[i + 1];
                }

                var leftOffset = x - current.X;
                var rightOffset = next.X - x;

                var matchedWithAPoint = false;
                if (Math.Abs(leftOffset) < epsilon)
                {
                    result.Add(current.Z);
                    matchedWithAPoint = true;
                }
                if (Math.Abs(rightOffset) < epsilon)
                {
                    result.Add(next.Z);
                    matchedWithAPoint = true;
                }
                if (!matchedWithAPoint)
                {
                    if (leftOffset > 0 && rightOffset > 0)
                    {
                        var fraction = leftOffset / (leftOffset + rightOffset);

                        result.Add((1.0 - fraction) * current.Z + fraction * next.Z);
                    }

                    // if both ofsets are negative the waterline goes back
                    if ((leftOffset < 0) && (rightOffset < 0))
                    {
                        var fraction = rightOffset / (rightOffset + leftOffset);

                        result.Add((1.0 - fraction) * next.Z + fraction * current.Z);
                    }
                }
            }

            return result.Distinct().ToList();
        }

        private static bool IsSegmentNotCrossingZ(double z, Point2D current, Point2D next)
        {
            if (double.IsNaN(z))
            {
                return true;
            }

            var leftOffset = Math.Abs(current.Z - z);
            var rightOffset = Math.Abs(next.Z - z);

            int currentSign = leftOffset < epsilon ? 0 : Math.Sign(current.Z - z);
            int nextSign = rightOffset < epsilon ? 0 : Math.Sign(next.Z - z);
            return currentSign == nextSign && currentSign != 0;
        }

        private static double GetXIntersectingZ(double z, Point2D current, Point2D next)
        {
            var leftOffset = Math.Abs(current.Z - z);
            if (leftOffset < epsilon)
            {
                return current.X;
            }
            var rightOffset = Math.Abs(next.Z - z);
            if (rightOffset < epsilon)
            {
                return next.X;
            }

            var fraction = leftOffset / (leftOffset + rightOffset);
            return GeneralMathRoutines.LinearInterpolate(current.X, next.X, fraction);
        }

        private IList<Point2D> IntersectionPointsWithLineCore(Line line, bool allowDuplicates)
        {
            var intersectionPointsWithLine = new List<Point2D>();

            for (int pointIndex = 0; pointIndex < calcPoints.Count - 1; pointIndex++)
            {
                DoIntersectAndAddToCollection(line, calcPoints[pointIndex], calcPoints[pointIndex + 1], 
                    intersectionPointsWithLine, allowDuplicates);
            }
            return intersectionPointsWithLine;
        }

        private static void DoIntersectAndAddToCollection(Line line, Point2D begin, Point2D end, ICollection<Point2D> intersectionPointsWithLine, bool allowDuplicates)
        {
            var lineInPoly = new Line
                {
                    BeginPoint = begin,
                    EndPoint = end
                };
            Point2D intersectionPoint = lineInPoly.GetIntersectPointXz(line);
            if (intersectionPoint != null && (allowDuplicates || NoPointSameXzLocation(intersectionPointsWithLine, intersectionPoint)))
            {
                intersectionPointsWithLine.Add(intersectionPoint);
            }
        }

        private static bool NoPointSameXzLocation(IEnumerable<Point2D> collection, Point2D point)
        {
            return !collection.Any(
                p => Math.Abs(p.X - point.X) < GeometryConstants.Accuracy &&
                     Math.Abs(p.Z - point.Z) < GeometryConstants.Accuracy);
        }

        private List<Point2D> IntersectWithPointsListCore(IList<Point2D> list, bool closePointString)
        {
            var result = new List<Point2D>();
            var line = new Line();
            for (int externalPointIndex = 0; externalPointIndex < list.Count - 1; externalPointIndex++)
            {
                line.BeginPoint = list[externalPointIndex];
                line.EndPoint = list[externalPointIndex + 1];

                result.AddRange(!closePointString
                                    ? IntersectionPointsXzWithLineXzWithAllPoints(line)
                                    : IntersectionPointsXzClosedStringWithLineXz(line));
            }

            return result;
        }

        /// <summary>
        /// Synchronizes the calculation points.
        /// </summary>
        public void SyncCalcPoints()
        {
            calcPoints.Clear();
            foreach (var geometryPoint in Points)
            {
                var p2D = new Point2D
                {
                    X = geometryPoint.X,
                    Z = geometryPoint.Z
                        
                };
                calcPoints.Add(p2D);
            }
        }

        /// <summary>
        /// Synchronizes the points.
        /// </summary>
        public void SyncPoints()
        {
            points.Clear();
            foreach (var p2D in calcPoints)
            {
                var geometryPoint = new GeometryPoint
                {
                    X = p2D.X,
                    Z = p2D.Z

                };
                points.Add(geometryPoint);
            }
        }
    }
}