// 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.Linq;
using Deltares.DamEngine.Data.Standard;

namespace Deltares.DamEngine.Data.Geometry;

/// <summary>
/// Type of intersection
/// </summary>
public enum LineIntersection
{
    /// <summary>
    /// The no intersection
    /// </summary>
    NoIntersection = 0,

    /// <summary>
    /// The intersects
    /// </summary>
    Intersects = 1,

    /// <summary>
    /// The parallel
    /// </summary>
    Parallel = 2
}

/// <summary>
/// Position of point towards polygon
/// </summary>
public enum PointInPolygon
{
    /// <summary>
    /// The inside polygon
    /// </summary>
    InsidePolygon = 0,

    /// <summary>
    /// The on polygon edge
    /// </summary>
    OnPolygonEdge = 1,

    /// <summary>
    /// The outside polygon
    /// </summary>
    OutsidePolygon = 2
}

/// <summary>
/// Enum for the clockwise options
/// </summary>
public enum Clockwise
{
    IsClockwise = 0,
    AntiClockwise = 1,
    NotEnoughUniquePoints = 2,
    PointsOnLine = 3
}

/// <summary>
/// Static class Routines2D
/// </summary>
public static class Routines2D
{
    private const double tolerance = 1.0e-4;

    /// <summary>
    /// Checks if the 2D Lines strictly intersect. Strictly means that the intersection point must be part of both lines so
    /// extrapolated points do not count.
    /// </summary>
    /// <param name="point1">The point1.</param>
    /// <param name="point2">The point2.</param>
    /// <param name="point3">The point3.</param>
    /// <param name="point4">The point4.</param>
    /// <param name="intersectionPoint">The intersection point.</param>
    /// <returns>
    /// Intersection status as well as the intersection point (when found, else null)
    /// </returns>
    /// <remarks>
    /// For connected parallel lines, the connection point will be returned as valid intersection point.
    /// </remarks>
    public static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3,
                                                                       Point2D point4, out Point2D intersectionPoint)
    {
        return DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4,
                                                   out intersectionPoint, tolerance);
    }

    /// <summary>
    /// Does the point exist in line (within the given tolerance).
    /// </summary>
    /// <param name="linePoint1">The line point1.</param>
    /// <param name="linePoint2">The line point2.</param>
    /// <param name="point">The point.</param>
    /// <param name="tolerance">The tolerance.</param>
    /// <returns>True when point is on line or equals a line point</returns>
    public static bool DoesPointExistInLine(Point2D linePoint1, Point2D linePoint2, Point2D point, double tolerance)
    {
        // Function to find if lies on or close to a line (within a given tolerance)
        // Input - a line defined by lTwo points, a third point to test and tolerance aValue
        // Output - TRUE or FALSE

        // first check whether the point is equal to a line point
        if (linePoint1.LocationEquals(point, tolerance) || linePoint2.LocationEquals(point, tolerance))
        {
            return true;
        }
        
        double linePoint1X = linePoint1.X;
        double linePoint1Z = linePoint1.Z;
        double linePoint2X = linePoint2.X;
        double linePoint2Z = linePoint2.Z;
        double pointX = point.X;
        double pointZ = point.Z;
        
        // then check if point is within the bounding rectangle formed by 2 line points
        double maxX = Math.Max(linePoint1X, linePoint2X);
        double minX = Math.Min(linePoint1X, linePoint2X);
        double maxY = Math.Max(linePoint1Z, linePoint2Z);
        double minY = Math.Min(linePoint1Z, linePoint2Z);

        minX -= tolerance;
        minY -= tolerance;
        maxX += tolerance;
        maxY += tolerance;

        double x3 = pointX;
        double y3 = pointZ;

        if (!(x3 > minX && x3 < maxX && y3 > minY && y3 < maxY))
        {
            return false;
        }

        // check if perpendicular distance between point and line is within tolerance
        double a = linePoint1Z - linePoint2Z;
        double b = linePoint2X - linePoint1X;
        double c = -((a * linePoint1X) + (b * linePoint1Z));
        double distance = ((a * x3) + (b * y3) + c) / Compute2DDistance(linePoint1X, linePoint1Z, linePoint2X, linePoint2Z);

        return Math.Abs(distance) < tolerance;
    }

    /// <summary>
    /// To check if the points of the polygon are clock-wise
    /// </summary>
    /// <param name="aPolygon"></param>
    /// <returns></returns>
    public static Clockwise IsClockWise(IEnumerable<Point2D> aPolygon)
    {
        Point2D[] distinctPoints = aPolygon.Distinct().ToArray();
        if (distinctPoints.Length < 3)
        {
            return Clockwise.NotEnoughUniquePoints;
        }

        Point2D p0 = distinctPoints.First();
        var realDistinctPoint = new Point2D[distinctPoints.Length];
        realDistinctPoint[0] = p0;
        var index = 1;
        foreach (Point2D distinctPoint in distinctPoints)
        {
            if (distinctPoint != p0)
            {
                if (!distinctPoint.LocationEquals(p0))
                {
                    realDistinctPoint[index] = distinctPoint;
                    index++;
                }
            }
        }

        if (index < 3)
        {
            return Clockwise.NotEnoughUniquePoints;
        }

        double sumClockwise = 0;
        double clockwise;
        for (var ii = 0; ii < distinctPoints.Length - 1; ii++)
        {
            clockwise = (distinctPoints[ii + 1].X - distinctPoints[ii].X) *
                        (distinctPoints[ii + 1].Z + distinctPoints[ii].Z);
            sumClockwise += clockwise;
        }

        clockwise = (distinctPoints[0].X - distinctPoints[distinctPoints.Length - 1].X) *
                    (distinctPoints[0].Z + distinctPoints[distinctPoints.Length - 1].Z);

        sumClockwise += clockwise;
        int result = Math.Sign(sumClockwise);

        if (result == 0)
        {
            return Clockwise.PointsOnLine;
        }

        return result > 0.0 ? Clockwise.IsClockwise : Clockwise.AntiClockwise;
    }

    /// <summary>
    /// Find if points Ax,ay are between a boundary polygon.
    /// </summary>
    /// <param name="polygon">The polygon.</param>
    /// <param name="x">The x.</param>
    /// <param name="z">The z.</param>
    /// <returns></returns>
    public static PointInPolygon CheckIfPointIsInPolygon(GeometryLoop polygon, double x, double z)
    {
        // This is done by calculating and adding the angles from the point to all points
        // of the polygon (using the ArcTan2 function). If the sum of these angles is
        // approximate zero the point lies outside the polygon.
        // If the sum in approximate + or - 2Pi the point lies in the polygon;
        // Performance tests have proven that this algorithm performs better than the 
        // polynomial winding number algorithm described at:
        // http://geomalgorithms.com/a03-_inclusion.html

        const double epsilon = 1e-10;
        var result = PointInPolygon.OutsidePolygon;

        // add the last point at the end, cos the polygon must be closed, so the last node equals the first
        var pointAdded = false;
        if (polygon.Points.Count > 0)
        {
            if (!polygon.Points[0].LocationEquals(polygon.Points[polygon.Points.Count - 1]))
            {
                polygon.Points.Add(polygon.Points[0]);
                pointAdded = true;
            }
        }

        if (polygon.Points.Count > 2)
        {
            double deltaX = polygon[0].X - x;
            double deltaZ = polygon[0].Z - z;
            double absX = Math.Abs(deltaX);
            double absZ = Math.Abs(deltaZ);
            if ((absX < epsilon) && (absZ < epsilon))
            {
                UndoAddIfNeeded(polygon, pointAdded);
                return PointInPolygon.OnPolygonEdge;
            }

            double al1 = Math.Atan2(deltaZ, deltaX);
            double som = 0;
            var index = 1;

            while (index < polygon.Points.Count)
            {
                if (polygon[index] != null)
                {
                    deltaX = polygon[index].X - x;
                    deltaZ = polygon[index].Z - z;
                    absX = Math.Abs(deltaX);
                    absZ = Math.Abs(deltaZ);
                    if ((absX < epsilon) && (absZ < epsilon))
                    {
                        UndoAddIfNeeded(polygon, pointAdded);
                        return PointInPolygon.OnPolygonEdge;
                    }

                    double al2 = Math.Atan2(deltaZ, deltaX);
                    double al3 = al2 - al1;

                    if (al3 < -Math.PI)
                    {
                        al3 = al3 + (2.0 * Math.PI);
                    }

                    if (al3 > Math.PI)
                    {
                        al3 = al3 - (2.0 * Math.PI);
                    }

                    if (((Math.PI - al3) < epsilon) || ((Math.PI + al3) < epsilon))
                    {
                        UndoAddIfNeeded(polygon, pointAdded);
                        return PointInPolygon.OnPolygonEdge;
                    }

                    som = som + al3;
                    al1 = al2;
                }

                index++;
            }

            if ((som > (1.9 * Math.PI)) || (som < (-1.9 * Math.PI)))
            {
                result = PointInPolygon.InsidePolygon;
            }
            else
            {
                result = PointInPolygon.OutsidePolygon;
            }
        }

        UndoAddIfNeeded(polygon, pointAdded);
        return result;
    }

    /// <summary>
    /// Intersection of Circle and line
    /// </summary>
    /// <param name="aX"></param>
    /// <param name="aY"></param>
    /// <param name="aR"></param>
    /// <param name="aX1"></param>
    /// <param name="aX2"></param>
    /// <param name="aY1"></param>
    /// <param name="aY2"></param>
    public static List<Point2D> IntersectCircleline(double aX, double aY, double aR, double aX1, double aX2, double aY1, double aY2)
    {
        // Solve by filling in parametric equation of line :
        // X = x1 + u*(x2 - x1)
        // Y = y1 + u*(y2 - y1)
        // Compare it to that of a cricle:
        // (X - xm)^2 + (Y - ym)^2 = r^2
        // And check solve for parameter u in a quadratic equation.
        // - no solutions => no intersections
        // - one or two solutions => intersection when -eps <= u <= 1+eps.
        // If (x1,y1) close to circle edge, return (x1,y1).

        const double lAbcEps = 1E-8;
        double lX1 = aX2 - aX1;
        double lX2 = aX1 - aX;
        double lY1 = aY2 - aY1;
        double lY2 = aY1 - aY;
        var result = new List<Point2D>();
        if ((Math.Abs(lX1) > lAbcEps) || (Math.Abs(lY1) > lAbcEps))
        {
            double lA = lX1 * lX1 + lY1 * lY1;
            double lB = 2 * (lX1 * lX2 + lY1 * lY2);
            double lC = lX2 * lX2 + lY2 * lY2 - aR * aR;
            double lD = lB * lB - 4 * lA * lC;

            if (lD > lAbcEps)
            {
                double lU = (-lB + Math.Sqrt(lD)) / (2 * lA);
                if ((lU >= -lAbcEps) && (lU <= (1.0 + lAbcEps)))
                {
                    result.Add(new Point2D
                    {
                        X = aX1 + lU * lX1,
                        Z = aY1 + lU * lY1
                    });
                }

                lU = (-lB - Math.Sqrt(lD)) / (2 * lA);

                if ((lU >= -lAbcEps) && (lU <= (1.0 + lAbcEps)))
                {
                    result.Add(new Point2D
                    {
                        X = aX1 + lU * lX1,
                        Z = aY1 + lU * lY1
                    });
                }
            }
            else if (Math.Abs(lD) <= lAbcEps)
            {
                double lU = (-lB) / (2 * lA);
                if ((lU >= -lAbcEps) && (lU <= 1.0 + lAbcEps))
                {
                    result.Add(new Point2D
                    {
                        X = aX1 + lU * lX1,
                        Z = aY1 + lU * lY1
                    });
                }
            }
        }
        else if (Math.Abs(Math.Pow(lX2, 2) + Math.Pow(lY2, 2) - Math.Pow(aR, 2)) < lAbcEps)
        {
            result.Add(new Point2D
            {
                X = aX1,
                Z = aY1
            });
        }

        return result;
    }

    /// <summary>
    ///  Determines angle between 2 lines.
    ///  Clockwise Input: point1, commonpoint, point2, commonpoint, normalvect(ex [0,0,-1]), validate  
    /// </summary>
    /// <param name="line1Point1"></param>
    /// <param name="line1Point2"></param>
    /// <param name="line2Point1"></param>
    /// <param name="line2Point2"></param>
    /// <returns></returns>
    public static double FindAngle(Point2D line1Point1, Point2D line1Point2, Point2D line2Point1,
                                   Point2D line2Point2)
    {
        double x1 = line1Point2.X - line1Point1.X;
        double z1 = line1Point2.Z - line1Point1.Z;

        double x2 = line2Point2.X - line2Point1.X;
        double z2 = line2Point2.Z - line2Point1.Z;

        double angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1));
        double angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2));

        double angle = angle2 - angle1;

        if (angle < 0)
        {
            angle += 360;
        }

        return angle;
    }

    /// <summary>
    /// Checks if values are equal
    /// </summary>
    /// <param name="x1"></param>
    /// <param name="x2"></param>
    /// <param name="tolerance"></param>
    /// <returns></returns>
    public static bool AreEqual(double x1, double x2, double tolerance)
    {
        return (Math.Abs(x1 - x2) < tolerance);
    }

    /// <summary>
    /// Determines the if 2D lines intersect using extrapolation when needed.
    /// </summary>
    /// <param name="aPoint1">a point1.</param>
    /// <param name="aPoint2">a point2.</param>
    /// <param name="aPoint3">a point3.</param>
    /// <param name="aPoint4">a point4.</param>
    /// <param name="aIntersectionPoint">a intersection point.</param>
    /// <returns></returns>
    public static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(Point2D aPoint1, Point2D aPoint2, Point2D aPoint3, Point2D aPoint4, out Point2D aIntersectionPoint)
    {
        LineConstant lLineConstant1 = CalculateNormalLineConstants(aPoint1, aPoint2);
        LineConstant lLineConstant2 = CalculateNormalLineConstants(aPoint3, aPoint4);

        LineIntersection lResult = DetermineIf2DLinesIntersectWithExtrapolation(lLineConstant1, lLineConstant2, out aIntersectionPoint);

        return lResult;
    }

    /// <summary>
    /// Finds the Intersection Points for two given 2D Lines.
    /// The intersection point does not need to be on either given line, extrapolation is used to find it beyond the length of the given lines.
    /// Note that parallel lines never return an intersection point, not even when there two parallel lines are neatly connected.
    /// </summary>
    /// <param name="aLine1Constant"></param>
    /// <param name="aLine2Constant"></param>
    /// <param name="aIntersectionPoint"></param>
    /// <returns></returns>
    public static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(LineConstant aLine1Constant, LineConstant aLine2Constant, out Point2D aIntersectionPoint)
    {
        aIntersectionPoint = new Point2D(0.0, 0.0);

        double lA1 = aLine1Constant.X;
        double lB1 = aLine1Constant.Y;
        double lC1 = aLine1Constant.Z;
        double lA2 = aLine2Constant.X;
        double lB2 = aLine2Constant.Y;
        double lC2 = aLine2Constant.Z;

        if (AreLinesParallel(lA1, lA2, lB1, lB2, tolerance))
        {
            aIntersectionPoint = null;
            return LineIntersection.Parallel;
        }

        double lX = (lB2 * lC1 - lB1 * lC2) / (lA1 * lB2 - lA2 * lB1);
        double lY = (lC1 * lA2 - lC2 * lA1) / (lA2 * lB1 - lA1 * lB2);

        aIntersectionPoint.X = lX;
        aIntersectionPoint.Z = lY;

        return LineIntersection.Intersects;
    }

    /// <summary>
    /// Checks if the points coincide
    /// </summary>
    /// <param name="point1X">The point1 x.</param>
    /// <param name="point1Z">The point1 z.</param>
    /// <param name="point2X">The point2 x.</param>
    /// <param name="point2Z">The point2 z.</param>
    /// <param name="tolerance">The tolerance.</param>
    /// <returns>true when points coincide</returns>
    public static bool DetermineIfPointsCoincide(double point1X, double point1Z, double point2X, double point2Z, double tolerance)
    {
        if ((Math.Abs(point1X - point2X)) < tolerance && Math.Abs(point1Z - point2Z) < tolerance)
        {
            return true;
        }

        return false;
    }

    private static void UndoAddIfNeeded(GeometryLoop polygon, bool needed)
    {
        if (needed)
        {
            polygon.Points.Remove(polygon.Points[polygon.Points.Count - 1]);
        }
    }

    /// <summary>
    /// Takes the cross product of two vectors
    /// </summary>
    /// <param name="pointAx">X-coordinate of the first point</param>
    /// <param name="pointAy">Y-coordinate of the first point</param>
    /// <param name="pointBx">X-coordinate of the second point</param>
    /// <param name="pointBy">Y-coordinate of the second point</param>
    /// <returns></returns>
    private static double CrossProduct(double pointAx, double pointAy, double pointBx, double pointBy)
    {
        return pointAx * pointBy - pointBx * pointAy;
    }

    /// <summary>
    /// Description: Finds the Normal. FPoint1 and FPoint2 are the line coordinates
    /// </summary>
    /// <param name="point1">The point1.</param>
    /// <param name="point2">The point2.</param>
    /// <returns></returns>
    private static LineConstant CalculateNormalLineConstants(Point2D point1, Point2D point2)
    {
        double point1X = point1.X;
        double point1Z = point1.Z;
        double point2X = point2.X;
        double point2Z = point2.Z;

        var lLineConstant = new LineConstant
        {
            X = point2Z - point1Z,
            Y = -(point2X - point1X),
            Z = (point2Z - point1Z) * point1X - (point2X - point1X) * point1Z
        };

        return lLineConstant;
    }

    /// <summary>
    /// Checks if lTwo lines are parallel. The Normal Line Constants are required.
    /// </summary>
    /// <param name="aLine1ConstantX">a line1 constant x.</param>
    /// <param name="aLine1ConstantY">a line1 constant y.</param>
    /// <param name="aLine2ConstantX">a line2 constant x.</param>
    /// <param name="aLine2ConstantY">a line2 constant y.</param>
    /// <param name="tolerance">The tolerance.</param>
    /// <returns></returns>
    private static bool AreLinesParallel(double aLine1ConstantX, double aLine1ConstantY, double aLine2ConstantX, double aLine2ConstantY, double tolerance)
    {
        return Math.Abs(CrossProduct(aLine1ConstantX, aLine2ConstantX, aLine1ConstantY, aLine2ConstantY)) < tolerance;
    }

    /// <summary>
    /// Calculates the Euclidean distance between two points on a 2-dimensional plane
    /// </summary>
    /// <param name="aX1">X-coordinate of the first point</param>
    /// <param name="aY1">Y-coordinate of the first point</param>
    /// <param name="aX2">X-coordinate of the second point</param>
    /// <param name="aY2">Y-coordinate of the second point</param>
    /// <returns>The distance between the given points.</returns>
    public static double Compute2DDistance(double aX1, double aY1, double aX2, double aY2)
    {
        double lX = aX1 - aX2;
        double lY = aY1 - aY2;

        return Math.Sqrt(lX * lX + lY * lY);
    }

    /// <summary>
    /// Checks if the 2D Lines strictly intersect. Strictly means that the intersection point must be part of both lines so
    /// extrapolated points do not count.
    /// </summary>
    /// <param name="point1">The point1.</param>
    /// <param name="point2">The point2.</param>
    /// <param name="point3">The point3.</param>
    /// <param name="point4">The point4.</param>
    /// <param name="intersectionPoint">The intersection point.</param>
    /// <param name="tolerance">Non-inclusive maximum allowed distance between two possibly intersecting endpoints of the lines</param>
    /// <returns>
    /// Intersection status as well as the intersection point (when found, else null)
    /// </returns>
    /// <remarks>
    /// For connected parallel lines, the connection point will be returned as valid intersection point.
    /// </remarks>
    private static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3,
                                                                        Point2D point4, out Point2D intersectionPoint, double tolerance)
    {
        double point1X = point1.X;
        double point1Z = point1.Z;
        double point2X = point2.X;
        double point2Z = point2.Z;
        double point3X = point3.X;
        double point3Z = point3.Z;
        double point4X = point4.X;
        double point4Z = point4.Z;
        LineConstant lineConstant1 = CalculateNormalLineConstants(point1, point2);
        LineConstant lineConstant2 = CalculateNormalLineConstants(point3, point4);

        LineIntersection result = DetermineIf2DLinesIntersectWithExtrapolation(lineConstant1, lineConstant2, out intersectionPoint);

        if (result == LineIntersection.Intersects)
        {
            //check if lies on the line and is not somewhere outside !
            if (!(DoesPointExistInLine(point1, point2, intersectionPoint, tolerance) &&
                  DoesPointExistInLine(point3, point4, intersectionPoint, tolerance)))
            {
                intersectionPoint = null;
                result = LineIntersection.NoIntersection;
            }
        }
        else
        {
            if (result == LineIntersection.Parallel && !DoLinesAtLeastPartialyOverlap(point1X, point1Z, point2X, point2Z, point3X, point3Z,
                                                                                      point4X, point4Z, tolerance))
            {
                if (DetermineIfPointsCoincide(point1X, point1Z, point3X, point3Z, tolerance) ||
                    DetermineIfPointsCoincide(point1X, point1Z, point4X, point4Z, tolerance))
                {
                    intersectionPoint = new Point2D(point1X, point1Z);
                    result = LineIntersection.Intersects;
                }

                if (DetermineIfPointsCoincide(point2X, point2Z, point3X, point3Z, tolerance) ||
                    DetermineIfPointsCoincide(point2X, point2Z, point4X, point4Z, tolerance))
                {
                    intersectionPoint = new Point2D(point2X, point2Z);
                    result = LineIntersection.Intersects;
                }
            }
        }

        return result;
    }

    /// <summary>
    /// Determines whether the given lines at least partialy overlap.
    /// Note that purely connected parallel lines are NOT considered to be overlapping.
    /// </summary>
    /// <param name="point1X">The point1 x.</param>
    /// <param name="point1Z">The point1 z.</param>
    /// <param name="point2X">The point2 x.</param>
    /// <param name="point2Z">The point2 z.</param>
    /// <param name="point3X">The point3 x.</param>
    /// <param name="point3Z">The point3 z.</param>
    /// <param name="point4X">The point4 x.</param>
    /// <param name="point4Z">The point4 z.</param>
    /// <param name="tolerance">The tolerance.</param>
    /// <returns></returns>
    private static bool DoLinesAtLeastPartialyOverlap(double point1X, double point1Z, double point2X, double point2Z,
                                                      double point3X, double point3Z, double point4X, double point4Z, double tolerance)
    {
        bool result = AreLinesParallel(point1X, point1Z, point2X, point2Z, point3X, point3Z, point4X, point4Z, tolerance);

        // lines not parallel so can not overlap
        if (!result)
        {
            return false;
        }

        if (Math.Abs(point1X - point2X) < double.Epsilon)
        {
            // lines are vertical so check Y values
            double max12 = Math.Max(point1Z, point2Z);
            double min12 = Math.Min(point1Z, point2Z);
            double max34 = Math.Max(point3Z, point4Z);
            double min34 = Math.Min(point3Z, point4Z);
            // When max of point1/point2 does not exceed min of point3/point4 lines do not overlap
            // When max of point3/point4 does not exceed min of point1/point2 lines do not overlap
            if (max12 <= min34 || min12 >= max34)
            {
                result = false;
            }
        }
        else
        {
            // lines are not vertical so check X values
            double max12 = Math.Max(point1X, point2X);
            double min12 = Math.Min(point1X, point2X);
            double max34 = Math.Max(point3X, point4X);
            double min34 = Math.Min(point3X, point4X);
            if (max12 <= min34 || min12 >= max34)
            {
                result = false;
            }
        }

        return result;
    }

    /// <summary>
    /// Description: Checks if two lines are parallel within the given tolerance
    /// </summary>
    /// <param name="point1X">The point1 x.</param>
    /// <param name="point1Z">The point1 z.</param>
    /// <param name="point2X">The point2 x.</param>
    /// <param name="point2Z">The point2 z.</param>
    /// <param name="point3X">The point3 x.</param>
    /// <param name="point3Z">The point3 z.</param>
    /// <param name="point4X">The point4 x.</param>
    /// <param name="point4Z">The point4 z.</param>
    /// <param name="tolerance">The tolerance.</param>
    /// <returns>
    /// True when parallel
    /// </returns>
    private static bool AreLinesParallel(double point1X, double point1Z, double point2X, double point2Z,
                                         double point3X, double point3Z, double point4X, double point4Z, double tolerance)
    {
        double aX = point2X - point1X;
        double aY = point2Z - point1Z;

        double bX = point4X - point3X;
        double bY = point4Z - point3Z;

        Normalize(ref aX, ref aY);
        Normalize(ref bX, ref bY);
        return AreLinesParallel(aX, aY, bX, bY, tolerance);
    }

    /// <summary>
    /// Normalizes this instance.
    /// </summary>
    private static void Normalize(ref double pointX, ref double pointY)
    {
        double q = Math.Sqrt(pointX * pointX + pointY * pointY);

        if (!q.AlmostEquals(0, 1e-8))
        {
            pointX = pointX / q;
            pointY = pointY / q;
        }
    }

    /// <summary>
    /// Structure for line constant in 3D
    /// </summary>
    public struct LineConstant
    {
        public double X { get; set; }
        public double Y { get; set; }
        public double Z { get; set; }
    }
    
    /// <summary>
    /// Calculates the distance between a point and a line segment.
    /// </summary>
    /// <param name="pointX"></param>
    /// <param name="pointY"></param>
    /// <param name="lineStartX"></param>
    /// <param name="lineStartY"></param>
    /// <param name="lineEndX"></param>
    /// <param name="lineEndY"></param>
    /// <returns></returns>
    public static double CalculateDistanceToLine(
        double pointX,
        double pointY,
        double lineStartX,
        double lineStartY,
        double lineEndX,
        double lineEndY)
    {
        return Math.Sqrt(CalculateSquaredDistance(new Point2D(pointX, pointY), new Point2D(lineStartX, lineStartY), 
                                                  new Point2D(lineEndX, lineEndY)));
    }
    
    private static double CalculateSquaredDistance(Point2D point, Point2D lineStart, Point2D lineEnd)
    {
        Point2D distanceVectorLineStartToEnd  = lineEnd - lineStart;
        Point2D distanceVectorLineStartToPoint = point - lineStart;
        double DotProduct(Point2D first, Point2D second) => first.X * second.X + first.Z * second.Z;
        double dotProductLineStartAndEnd = DotProduct(distanceVectorLineStartToPoint, distanceVectorLineStartToEnd );
        double dotProductPointAndLineStart = DotProduct(distanceVectorLineStartToEnd , distanceVectorLineStartToEnd );
        Point2D point2D3;
        if (dotProductLineStartAndEnd <= 0.0)
            point2D3 = point - lineStart;
        else if (dotProductPointAndLineStart <= dotProductLineStartAndEnd)
        {
            point2D3 = point - lineEnd;
        }
        else
        {
            double vectorToClosestPointOnLine = dotProductLineStartAndEnd / dotProductPointAndLineStart;
            Point2D point2D4 = lineStart + vectorToClosestPointOnLine * distanceVectorLineStartToEnd ;
            point2D3 = point - point2D4;
        }
        return DotProduct(point2D3, point2D3);
    }
}