Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs =================================================================== diff -u -r3893 -r4000 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 3893) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 4000) @@ -35,10 +35,12 @@ /// The no intersection /// NoIntersection = 0, + /// /// The intersects /// Intersects = 1, + /// /// The parallel /// @@ -54,10 +56,12 @@ /// The inside polygon /// InsidePolygon = 0, + /// /// The on polygon edge /// OnPolygonEdge = 1, + /// /// The outside polygon /// @@ -80,13 +84,6 @@ /// public static class Routines2D { - public struct LineConstant - { - public double X { get; set; } - public double Y { get; set; } - public double Z { get; set; } - } - private const double CEpsilon = 1.0e-4; /// @@ -105,12 +102,12 @@ /// 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) + Point2D point4, out Point2D intersectionPoint) { - return DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4, - out intersectionPoint, CEpsilon); + return DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4, + out intersectionPoint, CEpsilon); } - + /// /// Does the point exist in line. /// @@ -124,7 +121,7 @@ // 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 - + double linePoint1X = linePoint1.X; double linePoint1Z = linePoint1.Z; double linePoint2X = linePoint2.X; @@ -139,58 +136,51 @@ } // then check if point is within the bounding rectangle formed by 2 line points - var maxX = Math.Max(linePoint1X, linePoint2X); - var minX = Math.Min(linePoint1X, linePoint2X); - var maxY = Math.Max(linePoint1Z, linePoint2Z); - var minY = Math.Min(linePoint1Z, linePoint2Z); + 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; - var x3 = pointX; - var y3 = pointZ; + 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 - var a = linePoint1Z - linePoint2Z; - var b = linePoint2X - linePoint1X; - var c = -((a * linePoint1X) + (b * linePoint1Z)); - var distance = ((a * x3) + (b * y3) + c) / Compute2DDistance(linePoint1X, linePoint1Z, linePoint2X, linePoint2Z); + 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; } - private static void UndoAddIfNeeded(GeometryLoop polygon, bool needed) - { - if (needed) - { - polygon.CalcPoints.Remove(polygon.CalcPoints[polygon.CalcPoints.Count - 1]); - } - } - /// /// To check if the points of the polygon are clock-wise /// /// /// public static Clockwise IsClockWise(IEnumerable aPolygon) { - var distinctPoints = aPolygon.Distinct().ToArray(); + Point2D[] distinctPoints = aPolygon.Distinct().ToArray(); if (distinctPoints.Length < 3) { return Clockwise.NotEnoughUniquePoints; } - var p0 = distinctPoints.First(); - Point2D[] realDistinctPoint = new Point2D[distinctPoints.Length]; + + Point2D p0 = distinctPoints.First(); + var realDistinctPoint = new Point2D[distinctPoints.Length]; realDistinctPoint[0] = p0; var index = 1; - foreach (var distinctPoint in distinctPoints) + foreach (Point2D distinctPoint in distinctPoints) { if (distinctPoint != p0) { @@ -201,24 +191,26 @@ } } } + if (index < 3) { return Clockwise.NotEnoughUniquePoints; } double sumClockwise = 0; double clockwise; - for (int ii = 0; ii < distinctPoints.Length - 1; ii++) + 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; - var result = Math.Sign(sumClockwise); + int result = Math.Sign(sumClockwise); if (result == 0) { @@ -259,20 +251,20 @@ } } - var polygonPoints = polygon.CalcPoints.Count; + int polygonPoints = polygon.CalcPoints.Count; if (polygonPoints > 2) { - var deltaX = polygon[0].X - x; - var deltaZ = polygon[0].Z - z; - var absX = Math.Abs(deltaX); - var absZ = Math.Abs(deltaZ); + 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; } - var al1 = Math.Atan2(deltaZ, deltaX); + double al1 = Math.Atan2(deltaZ, deltaX); double som = 0; var index = 1; @@ -287,36 +279,41 @@ UndoAddIfNeeded(polygon, pointAdded); return PointInPolygon.OnPolygonEdge; } - var al2 = Math.Atan2(deltaZ, deltaX); - var al3 = al2 - al1; + double al2 = Math.Atan2(deltaZ, deltaX); + double al3 = al2 - al1; + if (al3 < -Math.PI) { - al3 = al3 + (2.0*Math.PI); + al3 = al3 + (2.0 * Math.PI); } + if (al3 > Math.PI) { - al3 = al3 - (2.0*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))) + if ((som > (1.9 * Math.PI)) || (som < (-1.9 * Math.PI))) { result = PointInPolygon.InsidePolygon; } else { result = PointInPolygon.OutsidePolygon; - } + } } + UndoAddIfNeeded(polygon, pointAdded); return result; } @@ -344,45 +341,63 @@ // If (x1,y1) close to circle edge, return (x1,y1). const double lAbcEps = 1E-8; - var lX1 = aX2 - aX1; - var lX2 = aX1 - aX; - var lY1 = aY2 - aY1; - var lY2 = aY1 - aY; + 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)) { - var lA = lX1*lX1 + lY1*lY1; - var lB = 2*(lX1*lX2 + lY1*lY2); - var lC = lX2*lX2 + lY2*lY2 - aR*aR; - var lD = lB*lB - 4*lA*lC; + 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) { - var lU = (-lB + Math.Sqrt(lD))/(2*lA); + 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}); + result.Add(new Point2D + { + 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}); + result.Add(new Point2D + { + X = aX1 + lU * lX1, + Z = aY1 + lU * lY1 + }); } } else if (Math.Abs(lD) <= lAbcEps) { - var lU = (-lB)/(2*lA); + double lU = (-lB) / (2 * lA); if ((lU >= -lAbcEps) && (lU <= 1.0 + lAbcEps)) { - result.Add(new Point2D { X = aX1 + lU*lX1, Z = aY1 + lU*lY1}); + 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}); + result.Add(new Point2D + { + X = aX1, + Z = aY1 + }); } + return result; } @@ -396,18 +411,18 @@ /// /// public static double FindAngle(Point2D line1Point1, Point2D line1Point2, Point2D line2Point1, - Point2D line2Point2) + Point2D line2Point2) { - var x1 = line1Point2.X - line1Point1.X; - var z1 = line1Point2.Z - line1Point1.Z; + double x1 = line1Point2.X - line1Point1.X; + double z1 = line1Point2.Z - line1Point1.Z; - var x2 = line2Point2.X - line2Point1.X; - var z2 = line2Point2.Z - line2Point1.Z; + double x2 = line2Point2.X - line2Point1.X; + double z2 = line2Point2.Z - line2Point1.Z; - var angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1)); - var angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2)); + double angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1)); + double angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2)); - var angle = angle2 - angle1; + double angle = angle2 - angle1; if (angle < 0) { @@ -430,6 +445,87 @@ } /// + /// 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); + + LineIntersection lResult = DetermineIf2DLinesIntersectWithExtrapolation(lLineConstant1, lLineConstant2, out aIntersectionPoint); + + return lResult; + } + + /// + /// 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 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, CEpsilon)) + { + 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; + } + + /// + /// 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; + } + + return false; + } + + private static void UndoAddIfNeeded(GeometryLoop polygon, bool needed) + { + if (needed) + { + polygon.CalcPoints.Remove(polygon.CalcPoints[polygon.CalcPoints.Count - 1]); + } + } + + /// /// Takes the cross product of two vectors /// /// X-coordinate of the first point @@ -439,7 +535,7 @@ /// private static double CrossProduct(double pointAx, double pointAy, double pointBx, double pointBy) { - return pointAx*pointBy - pointBx*pointAy; + return pointAx * pointBy - pointBx * pointAy; } /// @@ -479,62 +575,7 @@ return Math.Abs(CrossProduct(aLine1ConstantX, aLine2ConstantX, aLine1ConstantY, aLine2ConstantY)) < tolerance; } - /// - /// 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) - { - var lLineConstant1 = CalculateNormalLineConstants(aPoint1, aPoint2); - var lLineConstant2 = CalculateNormalLineConstants(aPoint3, aPoint4); - - var lResult = DetermineIf2DLinesIntersectWithExtrapolation(lLineConstant1, lLineConstant2, out aIntersectionPoint); - - return lResult; - } - - /// - /// 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); - - var lA1 = aLine1Constant.X; - var lB1 = aLine1Constant.Y; - var lC1 = aLine1Constant.Z; - var lA2 = aLine2Constant.X; - var lB2 = aLine2Constant.Y; - var lC2 = aLine2Constant.Z; - - if (AreLinesParallel(lA1, lA2, lB1, lB2, CEpsilon)) - { - aIntersectionPoint = null; - return LineIntersection.Parallel; - } - - var lX = (lB2*lC1 - lB1*lC2)/(lA1*lB2 - lA2*lB1); - var lY = (lC1*lA2 - lC2*lA1)/(lA2*lB1 - lA1*lB2); - - aIntersectionPoint.X = lX; - aIntersectionPoint.Z = lY; - - return LineIntersection.Intersects; - } - - /// /// Calculates the Euclidean distance between two points on a 2-dimensional plane /// /// X-coordinate of the first point @@ -544,8 +585,8 @@ /// The distance between the given points. private static double Compute2DDistance(double aX1, double aY1, double aX2, double aY2) { - var lX = aX1 - aX2; - var lY = aY1 - aY2; + double lX = aX1 - aX2; + double lY = aY1 - aY2; return Math.Sqrt(lX * lX + lY * lY); } @@ -566,51 +607,53 @@ /// /// 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) + private static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, + Point2D point4, out Point2D intersectionPoint, double tolerance) { - var point1X = point1.X; - var point1Z = point1.Z; - var point2X = point2.X; - var point2Z = point2.Z; - var point3X = point3.X; - var point3Z = point3.Z; - var point4X = point4.X; - var point4Z = point4.Z; - var lineConstant1 = CalculateNormalLineConstants(point1, point2); - var lineConstant2 = CalculateNormalLineConstants(point3, point4); + 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); - var result = DetermineIf2DLinesIntersectWithExtrapolation(lineConstant1, lineConstant2, out intersectionPoint); + 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))) + 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 (result == LineIntersection.Parallel && !DoLinesAtLeastPartialyOverlap(point1X, point1Z, point2X, point2Z, point3X, point3Z, + point4X, point4Z, tolerance)) { - if (DetermineIfPointsCoincide(point1X, point1Z, point3X, point3Z, 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) || + + if (DetermineIfPointsCoincide(point2X, point2Z, point3X, point3Z, tolerance) || DetermineIfPointsCoincide(point2X, point2Z, point4X, point4Z, tolerance)) { intersectionPoint = new Point2D(point2X, point2Z); result = LineIntersection.Intersects; } } } + return result; } @@ -629,9 +672,9 @@ /// The tolerance. /// private static bool DoLinesAtLeastPartialyOverlap(double point1X, double point1Z, double point2X, double point2Z, - double point3X, double point3Z, double point4X, double point4Z, double tolerance) + double point3X, double point3Z, double point4X, double point4Z, double tolerance) { - var result = AreLinesParallel(point1X, point1Z, point2X, point2Z, point3X, point3Z, point4X, point4Z, tolerance); + bool result = AreLinesParallel(point1X, point1Z, point2X, point2Z, point3X, point3Z, point4X, point4Z, tolerance); // lines not parallel so can not overlap if (!result) @@ -642,10 +685,10 @@ if (Math.Abs(point1X - point2X) < double.Epsilon) { // lines are vertical so check Y values - var max12 = Math.Max(point1Z, point2Z); - var min12 = Math.Min(point1Z, point2Z); - var max34 = Math.Max(point3Z, point4Z); - var min34 = Math.Min(point3Z, point4Z); + 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) @@ -656,15 +699,16 @@ else { // lines are not vertical so check X values - var max12 = Math.Max(point1X, point2X); - var min12 = Math.Min(point1X, point2X); - var max34 = Math.Max(point3X, point4X); - var min34 = Math.Min(point3X, point4X); + 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; } @@ -684,15 +728,14 @@ /// 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 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); @@ -708,27 +751,15 @@ if (!q.AlmostEquals(0, 1e-8)) { pointX = pointX / q; - pointY = pointY / q; + pointY = pointY / q; } } - /// - /// 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) + public struct LineConstant { - if ((Math.Abs(point1X - point2X)) < tolerance && Math.Abs(point1Z - point2Z) < tolerance) - { - return true; - } - - return false; + public double X { get; set; } + public double Y { get; set; } + public double Z { get; set; } } } } \ No newline at end of file