Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs =================================================================== diff -u -r4000 -r4052 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 4000) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 4052) @@ -24,742 +24,741 @@ using System.Linq; using Deltares.DamEngine.Data.Standard; -namespace Deltares.DamEngine.Data.Geometry +namespace Deltares.DamEngine.Data.Geometry; + +/// +/// Type of intersection +/// +public enum LineIntersection { /// - /// Type of intersection + /// The no intersection /// - public enum LineIntersection - { - /// - /// The no intersection - /// - NoIntersection = 0, + NoIntersection = 0, - /// - /// The intersects - /// - Intersects = 1, + /// + /// The intersects + /// + Intersects = 1, - /// - /// The parallel - /// - Parallel = 2 - } + /// + /// The parallel + /// + Parallel = 2 +} +/// +/// Position of point towards polygon +/// +public enum PointInPolygon +{ /// - /// Position of point towards polygon + /// The inside polygon /// - public enum PointInPolygon - { - /// - /// The inside polygon - /// - InsidePolygon = 0, + InsidePolygon = 0, - /// - /// The on polygon edge - /// - OnPolygonEdge = 1, + /// + /// The on polygon edge + /// + OnPolygonEdge = 1, - /// - /// The outside polygon - /// - OutsidePolygon = 2 - } + /// + /// The outside polygon + /// + OutsidePolygon = 2 +} +/// +/// Enum for the clockwise options +/// +public enum Clockwise +{ + IsClockwise = 0, + AntiClockwise = 1, + NotEnoughUniquePoints = 2, + PointsOnLine = 3 +} + +/// +/// Static class Routines2D +/// +public static class Routines2D +{ + private const double CEpsilon = 1.0e-4; + /// - /// Enum for the clockwise options + /// 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. /// - public enum Clockwise + /// The point1. + /// The point2. + /// The point3. + /// The point4. + /// The intersection point. + /// + /// Intersection status as well as the intersection point (when found, else null) + /// + /// + /// For connected parallel lines, the connection point will be returned as valid intersection point. + /// + public static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, + Point2D point4, out Point2D intersectionPoint) { - IsClockwise = 0, - AntiClockwise = 1, - NotEnoughUniquePoints = 2, - PointsOnLine = 3 + return DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4, + out intersectionPoint, CEpsilon); } /// - /// Static class Routines2D + /// Does the point exist in line. /// - public static class Routines2D + /// The line point1. + /// The line point2. + /// The point. + /// The tolerance. + /// + public static bool DoesPointExistInLine(Point2D linePoint1, Point2D linePoint2, Point2D point, double tolerance) { - private const double CEpsilon = 1.0e-4; + // 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 - /// - /// 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. - /// - /// The point1. - /// The point2. - /// The point3. - /// The point4. - /// The intersection point. - /// - /// Intersection status as well as the intersection point (when found, else null) - /// - /// - /// For connected parallel lines, the connection point will be returned as valid intersection point. - /// - public static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, - Point2D point4, out Point2D intersectionPoint) + double linePoint1X = linePoint1.X; + double linePoint1Z = linePoint1.Z; + double linePoint2X = linePoint2.X; + double linePoint2Z = linePoint2.Z; + double pointX = point.X; + double pointZ = point.Z; + + // first check whether the points are identical + if (Math.Abs(linePoint1X - linePoint2X) < tolerance && Math.Abs(linePoint1Z - linePoint2Z) < tolerance) { - return DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4, - out intersectionPoint, CEpsilon); + return Compute2DDistance(pointX, pointZ, linePoint1X, linePoint1Z) < tolerance; } - /// - /// Does the point exist in line. - /// - /// The line point1. - /// The line point2. - /// The point. - /// The tolerance. - /// - 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 + // 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); - double linePoint1X = linePoint1.X; - double linePoint1Z = linePoint1.Z; - double linePoint2X = linePoint2.X; - double linePoint2Z = linePoint2.Z; - double pointX = point.X; - double pointZ = point.Z; + minX -= tolerance; + minY -= tolerance; + maxX += tolerance; + maxY += tolerance; - // first check whether the points are identical - if (Math.Abs(linePoint1X - linePoint2X) < tolerance && Math.Abs(linePoint1Z - linePoint2Z) < tolerance) - { - return Compute2DDistance(pointX, pointZ, linePoint1X, linePoint1Z) < tolerance; - } + double x3 = pointX; + double y3 = pointZ; - // 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); + if (!(x3 > minX && x3 < maxX && y3 > minY && y3 < maxY)) + { + return false; + } - minX -= tolerance; - minY -= tolerance; - maxX += tolerance; - maxY += tolerance; + // 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); - double x3 = pointX; - double y3 = pointZ; + return Math.Abs(distance) < tolerance; + } - 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; + /// + /// To check if the points of the polygon are clock-wise + /// + /// + /// + public static Clockwise IsClockWise(IEnumerable aPolygon) + { + Point2D[] distinctPoints = aPolygon.Distinct().ToArray(); + if (distinctPoints.Length < 3) + { + return Clockwise.NotEnoughUniquePoints; } - /// - /// To check if the points of the polygon are clock-wise - /// - /// - /// - public static Clockwise IsClockWise(IEnumerable aPolygon) + Point2D p0 = distinctPoints.First(); + var realDistinctPoint = new Point2D[distinctPoints.Length]; + realDistinctPoint[0] = p0; + var index = 1; + foreach (Point2D distinctPoint in distinctPoints) { - Point2D[] distinctPoints = aPolygon.Distinct().ToArray(); - if (distinctPoints.Length < 3) + if (distinctPoint != p0) { - 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)) { - if (!distinctPoint.LocationEquals(p0)) - { - realDistinctPoint[index] = distinctPoint; - index++; - } + realDistinctPoint[index] = distinctPoint; + index++; } } + } - if (index < 3) - { - return Clockwise.NotEnoughUniquePoints; - } + 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; - } + 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); + clockwise = (distinctPoints[0].X - distinctPoints[distinctPoints.Length - 1].X) * + (distinctPoints[0].Z + distinctPoints[distinctPoints.Length - 1].Z); - sumClockwise += clockwise; - int result = Math.Sign(sumClockwise); + sumClockwise += clockwise; + int result = Math.Sign(sumClockwise); - if (result == 0) + if (result == 0) + { + return Clockwise.PointsOnLine; + } + + return result > 0.0 ? Clockwise.IsClockwise : Clockwise.AntiClockwise; + } + + /// + /// Find if points Ax,ay are between a boundary polygon. + /// + /// The polygon. + /// The x. + /// The z. + /// + 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 + // aproximate zero the point lies outside the polygon. + // If the sum in aproximate + 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.CalcPoints.Count > 0) + { + if (!polygon.CalcPoints[0].LocationEquals(polygon.CalcPoints[polygon.CalcPoints.Count - 1])) { - return Clockwise.PointsOnLine; + polygon.CalcPoints.Add(polygon.CalcPoints[0]); + pointAdded = true; } - - return result > 0.0 ? Clockwise.IsClockwise : Clockwise.AntiClockwise; } - /// - /// Find if points Ax,ay are between a boundary polygon. - /// - /// The polygon. - /// The x. - /// The z. - /// - public static PointInPolygon CheckIfPointIsInPolygon(GeometryLoop polygon, double x, double z) + int polygonPoints = polygon.CalcPoints.Count; + if (polygonPoints > 2) { - // 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 - // aproximate zero the point lies outside the polygon. - // If the sum in aproximate + 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.CalcPoints.Count > 0) + 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)) { - if (!polygon.CalcPoints[0].LocationEquals(polygon.CalcPoints[polygon.CalcPoints.Count - 1])) - { - polygon.CalcPoints.Add(polygon.CalcPoints[0]); - pointAdded = true; - } + UndoAddIfNeeded(polygon, pointAdded); + return PointInPolygon.OnPolygonEdge; } - int polygonPoints = polygon.CalcPoints.Count; - if (polygonPoints > 2) + double al1 = Math.Atan2(deltaZ, deltaX); + double som = 0; + var index = 1; + + while (index < polygonPoints) { - double deltaX = polygon[0].X - x; - double deltaZ = polygon[0].Z - z; - double absX = Math.Abs(deltaX); - double absZ = Math.Abs(deltaZ); + 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 al1 = Math.Atan2(deltaZ, deltaX); - double som = 0; - var index = 1; + double al2 = Math.Atan2(deltaZ, deltaX); + double al3 = al2 - al1; - while (index < polygonPoints) + if (al3 < -Math.PI) { - 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++; + al3 = al3 + (2.0 * Math.PI); } - if ((som > (1.9 * Math.PI)) || (som < (-1.9 * Math.PI))) + if (al3 > Math.PI) { - result = PointInPolygon.InsidePolygon; + al3 = al3 - (2.0 * Math.PI); } - else + + if (((Math.PI - al3) < epsilon) || ((Math.PI + al3) < epsilon)) { - result = PointInPolygon.OutsidePolygon; + UndoAddIfNeeded(polygon, pointAdded); + return PointInPolygon.OnPolygonEdge; } + + som = som + al3; + al1 = al2; + index++; } - UndoAddIfNeeded(polygon, pointAdded); - return result; + if ((som > (1.9 * Math.PI)) || (som < (-1.9 * Math.PI))) + { + result = PointInPolygon.InsidePolygon; + } + else + { + result = PointInPolygon.OutsidePolygon; + } } - /// - /// Intersection of Circle and line - /// - /// - /// - /// - /// - /// - /// - /// - public static List IntersectCircleline(double aX, double aY, double aR, double aX1, double aX2, double aY1, double aY2) + UndoAddIfNeeded(polygon, pointAdded); + return result; + } + + /// + /// Intersection of Circle and line + /// + /// + /// + /// + /// + /// + /// + /// + public static List 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(); + if ((Math.Abs(lX1) > lAbcEps) || (Math.Abs(lY1) > lAbcEps)) { - // 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). + 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; - 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(); - if ((Math.Abs(lX1) > lAbcEps) || (Math.Abs(lY1) > lAbcEps)) + if (lD > 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))) { - double lU = (-lB + Math.Sqrt(lD)) / (2 * lA); - if ((lU >= -lAbcEps) && (lU <= (1.0 + lAbcEps))) + result.Add(new Point2D { - result.Add(new Point2D - { - X = aX1 + lU * lX1, - Z = aY1 + lU * lY1 - }); - } + X = aX1 + lU * lX1, + Z = aY1 + lU * lY1 + }); + } - lU = (-lB - Math.Sqrt(lD)) / (2 * lA); + 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) + if ((lU >= -lAbcEps) && (lU <= (1.0 + lAbcEps))) { - double lU = (-lB) / (2 * lA); - if ((lU >= -lAbcEps) && (lU <= 1.0 + lAbcEps)) + result.Add(new Point2D { - result.Add(new Point2D - { - X = aX1 + lU * lX1, - Z = aY1 + lU * lY1 - }); - } + 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) + else if (Math.Abs(lD) <= lAbcEps) { - result.Add(new Point2D + double lU = (-lB) / (2 * lA); + if ((lU >= -lAbcEps) && (lU <= 1.0 + lAbcEps)) { - X = aX1, - Z = aY1 - }); + result.Add(new Point2D + { + X = aX1 + lU * lX1, + Z = aY1 + lU * lY1 + }); + } } - - return result; } - - /// - /// Determines angle between 2 lines. - /// Clockwise Input: point1, commonpoint, point2, commonpoint, normalvect(ex [0,0,-1]), validate - /// - /// - /// - /// - /// - /// - public static double FindAngle(Point2D line1Point1, Point2D line1Point2, Point2D line2Point1, - Point2D line2Point2) + else if (Math.Abs(Math.Pow(lX2, 2) + Math.Pow(lY2, 2) - Math.Pow(aR, 2)) < lAbcEps) { - double x1 = line1Point2.X - line1Point1.X; - double z1 = line1Point2.Z - line1Point1.Z; + result.Add(new Point2D + { + X = aX1, + Z = aY1 + }); + } - double x2 = line2Point2.X - line2Point1.X; - double z2 = line2Point2.Z - line2Point1.Z; + return result; + } - double angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1)); - double angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2)); + /// + /// Determines angle between 2 lines. + /// Clockwise Input: point1, commonpoint, point2, commonpoint, normalvect(ex [0,0,-1]), validate + /// + /// + /// + /// + /// + /// + 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 angle = angle2 - angle1; + double x2 = line2Point2.X - line2Point1.X; + double z2 = line2Point2.Z - line2Point1.Z; - if (angle < 0) - { - angle += 360; - } + double angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1)); + double angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2)); - return angle; - } + double angle = angle2 - angle1; - /// - /// Checks if values are equal - /// - /// - /// - /// - /// - public static bool AreEqual(double x1, double x2, double tolerance) + if (angle < 0) { - return (Math.Abs(x1 - x2) < tolerance); + angle += 360; } - /// - /// Determines the if 2D lines intersect using extrapolation when needed. - /// - /// a point1. - /// a point2. - /// a point3. - /// a point4. - /// a intersection point. - /// - 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); + return angle; + } - LineIntersection lResult = DetermineIf2DLinesIntersectWithExtrapolation(lLineConstant1, lLineConstant2, out aIntersectionPoint); + /// + /// Checks if values are equal + /// + /// + /// + /// + /// + public static bool AreEqual(double x1, double x2, double tolerance) + { + return (Math.Abs(x1 - x2) < tolerance); + } - return lResult; - } + /// + /// Determines the if 2D lines intersect using extrapolation when needed. + /// + /// a point1. + /// a point2. + /// a point3. + /// a point4. + /// a intersection point. + /// + 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); - /// - /// 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. - /// - /// - /// - /// - /// - public static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(LineConstant aLine1Constant, LineConstant aLine2Constant, out Point2D aIntersectionPoint) - { - aIntersectionPoint = new Point2D(0.0, 0.0); + LineIntersection lResult = DetermineIf2DLinesIntersectWithExtrapolation(lLineConstant1, lLineConstant2, out aIntersectionPoint); - double lA1 = aLine1Constant.X; - double lB1 = aLine1Constant.Y; - double lC1 = aLine1Constant.Z; - double lA2 = aLine2Constant.X; - double lB2 = aLine2Constant.Y; - double lC2 = aLine2Constant.Z; + return lResult; + } - if (AreLinesParallel(lA1, lA2, lB1, lB2, CEpsilon)) - { - aIntersectionPoint = null; - return LineIntersection.Parallel; - } + /// + /// 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. + /// + /// + /// + /// + /// + public static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(LineConstant aLine1Constant, LineConstant aLine2Constant, out Point2D aIntersectionPoint) + { + aIntersectionPoint = new Point2D(0.0, 0.0); - double lX = (lB2 * lC1 - lB1 * lC2) / (lA1 * lB2 - lA2 * lB1); - double lY = (lC1 * lA2 - lC2 * lA1) / (lA2 * lB1 - lA1 * lB2); + double lA1 = aLine1Constant.X; + double lB1 = aLine1Constant.Y; + double lC1 = aLine1Constant.Z; + double lA2 = aLine2Constant.X; + double lB2 = aLine2Constant.Y; + double lC2 = aLine2Constant.Z; - aIntersectionPoint.X = lX; - aIntersectionPoint.Z = lY; - - return LineIntersection.Intersects; + if (AreLinesParallel(lA1, lA2, lB1, lB2, CEpsilon)) + { + aIntersectionPoint = null; + return LineIntersection.Parallel; } - /// - /// Checks if the points coincide - /// - /// The point1 x. - /// The point1 z. - /// The point2 x. - /// The point2 z. - /// The tolerance. - /// - 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; - } + double lX = (lB2 * lC1 - lB1 * lC2) / (lA1 * lB2 - lA2 * lB1); + double lY = (lC1 * lA2 - lC2 * lA1) / (lA2 * lB1 - lA1 * lB2); - return false; - } + aIntersectionPoint.X = lX; + aIntersectionPoint.Z = lY; - private static void UndoAddIfNeeded(GeometryLoop polygon, bool needed) - { - if (needed) - { - polygon.CalcPoints.Remove(polygon.CalcPoints[polygon.CalcPoints.Count - 1]); - } - } + return LineIntersection.Intersects; + } - /// - /// Takes the cross product of two vectors - /// - /// X-coordinate of the first point - /// Y-coordinate of the first point - /// X-coordinate of the second point - /// Y-coordinate of the second point - /// - private static double CrossProduct(double pointAx, double pointAy, double pointBx, double pointBy) + /// + /// Checks if the points coincide + /// + /// The point1 x. + /// The point1 z. + /// The point2 x. + /// The point2 z. + /// The tolerance. + /// + 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 pointAx * pointBy - pointBx * pointAy; + return true; } - /// - /// Description: Finds the Normal. FPoint1 and FPoint2 are the line coordinates - /// - /// The point1. - /// The point2. - /// - private static LineConstant CalculateNormalLineConstants(Point2D point1, Point2D point2) - { - double point1X = point1.X; - double point1Z = point1.Z; - double point2X = point2.X; - double point2Z = point2.Z; + return false; + } - var lLineConstant = new LineConstant - { - X = point2Z - point1Z, - Y = -(point2X - point1X), - Z = (point2Z - point1Z) * point1X - (point2X - point1X) * point1Z - }; - - return lLineConstant; - } - - /// - /// Checks if lTwo lines are parallel. The Normal Line Constants are required. - /// - /// a line1 constant x. - /// a line1 constant y. - /// a line2 constant x. - /// a line2 constant y. - /// The tolerance. - /// - private static bool AreLinesParallel(double aLine1ConstantX, double aLine1ConstantY, double aLine2ConstantX, double aLine2ConstantY, double tolerance) + private static void UndoAddIfNeeded(GeometryLoop polygon, bool needed) + { + if (needed) { - return Math.Abs(CrossProduct(aLine1ConstantX, aLine2ConstantX, aLine1ConstantY, aLine2ConstantY)) < tolerance; + polygon.CalcPoints.Remove(polygon.CalcPoints[polygon.CalcPoints.Count - 1]); } + } - /// - /// Calculates the Euclidean distance between two points on a 2-dimensional plane - /// - /// X-coordinate of the first point - /// Y-coordinate of the first point - /// X-coordinate of the second point - /// Y-coordinate of the second point - /// The distance between the given points. - private static double Compute2DDistance(double aX1, double aY1, double aX2, double aY2) - { - double lX = aX1 - aX2; - double lY = aY1 - aY2; + /// + /// Takes the cross product of two vectors + /// + /// X-coordinate of the first point + /// Y-coordinate of the first point + /// X-coordinate of the second point + /// Y-coordinate of the second point + /// + private static double CrossProduct(double pointAx, double pointAy, double pointBx, double pointBy) + { + return pointAx * pointBy - pointBx * pointAy; + } - return Math.Sqrt(lX * lX + lY * lY); - } + /// + /// Description: Finds the Normal. FPoint1 and FPoint2 are the line coordinates + /// + /// The point1. + /// The point2. + /// + private static LineConstant CalculateNormalLineConstants(Point2D point1, Point2D point2) + { + double point1X = point1.X; + double point1Z = point1.Z; + double point2X = point2.X; + double point2Z = point2.Z; - /// - /// 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. - /// - /// The point1. - /// The point2. - /// The point3. - /// The point4. - /// The intersection point. - /// Non-inclusive maximum allowed distance between two possibly intersecting endpoints of the lines - /// - /// Intersection status as well as the intersection point (when found, else null) - /// - /// - /// For connected parallel lines, the connection point will be returned as valid intersection point. - /// - private static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, - Point2D point4, out Point2D intersectionPoint, double tolerance) + var lLineConstant = new LineConstant { - 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); + X = point2Z - point1Z, + Y = -(point2X - point1X), + Z = (point2Z - point1Z) * point1X - (point2X - point1X) * point1Z + }; - LineIntersection result = DetermineIf2DLinesIntersectWithExtrapolation(lineConstant1, lineConstant2, out intersectionPoint); + return lLineConstant; + } - 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; - } + /// + /// Checks if lTwo lines are parallel. The Normal Line Constants are required. + /// + /// a line1 constant x. + /// a line1 constant y. + /// a line2 constant x. + /// a line2 constant y. + /// The tolerance. + /// + private static bool AreLinesParallel(double aLine1ConstantX, double aLine1ConstantY, double aLine2ConstantX, double aLine2ConstantY, double tolerance) + { + return Math.Abs(CrossProduct(aLine1ConstantX, aLine2ConstantX, aLine1ConstantY, aLine2ConstantY)) < tolerance; + } - if (DetermineIfPointsCoincide(point2X, point2Z, point3X, point3Z, tolerance) || - DetermineIfPointsCoincide(point2X, point2Z, point4X, point4Z, tolerance)) - { - intersectionPoint = new Point2D(point2X, point2Z); - result = LineIntersection.Intersects; - } - } - } + /// + /// Calculates the Euclidean distance between two points on a 2-dimensional plane + /// + /// X-coordinate of the first point + /// Y-coordinate of the first point + /// X-coordinate of the second point + /// Y-coordinate of the second point + /// The distance between the given points. + private static double Compute2DDistance(double aX1, double aY1, double aX2, double aY2) + { + double lX = aX1 - aX2; + double lY = aY1 - aY2; - return result; - } + return Math.Sqrt(lX * lX + lY * lY); + } - /// - /// Determines whether the given lines at least partialy overlap. - /// Note that purely connected parallel lines are NOT considered to be overlapping. - /// - /// The point1 x. - /// The point1 z. - /// The point2 x. - /// The point2 z. - /// The point3 x. - /// The point3 z. - /// The point4 x. - /// The point4 z. - /// The tolerance. - /// - 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); + /// + /// 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. + /// + /// The point1. + /// The point2. + /// The point3. + /// The point4. + /// The intersection point. + /// Non-inclusive maximum allowed distance between two possibly intersecting endpoints of the lines + /// + /// Intersection status as well as the intersection point (when found, else null) + /// + /// + /// For connected parallel lines, the connection point will be returned as valid intersection point. + /// + 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); - // lines not parallel so can not overlap - if (!result) + 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))) { - return false; + intersectionPoint = null; + result = LineIntersection.NoIntersection; } - - if (Math.Abs(point1X - point2X) < double.Epsilon) + } + else + { + if (result == LineIntersection.Parallel && !DoLinesAtLeastPartialyOverlap(point1X, point1Z, point2X, point2Z, point3X, point3Z, + point4X, point4Z, tolerance)) { - // 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) + if (DetermineIfPointsCoincide(point1X, point1Z, point3X, point3Z, tolerance) || + DetermineIfPointsCoincide(point1X, point1Z, point4X, point4Z, tolerance)) { - result = false; + intersectionPoint = new Point2D(point1X, point1Z); + result = LineIntersection.Intersects; } - } - 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) + + if (DetermineIfPointsCoincide(point2X, point2Z, point3X, point3Z, tolerance) || + DetermineIfPointsCoincide(point2X, point2Z, point4X, point4Z, tolerance)) { - result = false; + intersectionPoint = new Point2D(point2X, point2Z); + result = LineIntersection.Intersects; } } - - return result; } - /// - /// Description: Checks if two lines are parallel within the given tolerance - /// - /// The point1 x. - /// The point1 z. - /// The point2 x. - /// The point2 z. - /// The point3 x. - /// The point3 z. - /// The point4 x. - /// The point4 z. - /// The tolerance. - /// - /// True when parallel - /// - 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; + return result; + } - double bX = point4X - point3X; - double bY = point4Z - point3Z; + /// + /// Determines whether the given lines at least partialy overlap. + /// Note that purely connected parallel lines are NOT considered to be overlapping. + /// + /// The point1 x. + /// The point1 z. + /// The point2 x. + /// The point2 z. + /// The point3 x. + /// The point3 z. + /// The point4 x. + /// The point4 z. + /// The tolerance. + /// + 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); - Normalize(ref aX, ref aY); - Normalize(ref bX, ref bY); - return AreLinesParallel(aX, aY, bX, bY, tolerance); + // lines not parallel so can not overlap + if (!result) + { + return false; } - /// - /// Normalizes this instance. - /// - private static void Normalize(ref double pointX, ref double pointY) + if (Math.Abs(point1X - point2X) < double.Epsilon) { - double q = Math.Sqrt(pointX * pointX + pointY * pointY); - - if (!q.AlmostEquals(0, 1e-8)) + // 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) { - pointX = pointX / q; - pointY = pointY / q; + 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; + } + } - public struct LineConstant + return result; + } + + /// + /// Description: Checks if two lines are parallel within the given tolerance + /// + /// The point1 x. + /// The point1 z. + /// The point2 x. + /// The point2 z. + /// The point3 x. + /// The point3 z. + /// The point4 x. + /// The point4 z. + /// The tolerance. + /// + /// True when parallel + /// + 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); + } + + /// + /// Normalizes this instance. + /// + private static void Normalize(ref double pointX, ref double pointY) + { + double q = Math.Sqrt(pointX * pointX + pointY * pointY); + + if (!q.AlmostEquals(0, 1e-8)) { - public double X { get; set; } - public double Y { get; set; } - public double Z { get; set; } + pointX = pointX / q; + pointY = pointY / q; } } + + public struct LineConstant + { + public double X { get; set; } + public double Y { get; set; } + public double Z { get; set; } + } } \ No newline at end of file