using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.Geometry;
namespace Deltares.Geotechnics.WaternetCreator
{
public static class LineHelper
{
///
/// Calculate intersection between two lines
///
///
///
///
///
public static GeometryPoint GetIntersectPoint(Line line1, Line line2, PointType pointType)
{
if (pointType == PointType.XZ)
{
double dx1 = line1.EndPoint.X - line1.BeginPoint.X;
double dz1 = line1.EndPoint.Z - line1.BeginPoint.Z;
double dx2 = line2.EndPoint.X - line2.BeginPoint.X;
double dz2 = line2.EndPoint.Z - line2.BeginPoint.Z;
double noem = dx1*dz2 - dz1*dx2;
double u, v;
if (Math.Abs(noem) > 0.0)
{
u = (dx2*(line1.BeginPoint.Z - line2.BeginPoint.Z) - dz2*(line1.BeginPoint.X - line2.BeginPoint.X))/
noem;
if ((u >= 0.0) && (u <= 1.0))
{
v = (dx1*(line1.BeginPoint.Z - line2.BeginPoint.Z) -
dz1*(line1.BeginPoint.X - line2.BeginPoint.X))/
noem;
if ((v >= 0.0) && (v <= 1.0))
{
return new GeometryPoint(line1.BeginPoint.X + u*dx1, 0, line1.BeginPoint.Z + u*dz1);
}
}
}
return new GeometryPoint(double.NaN, 0, double.NaN);
}
else if (pointType == PointType.XY)
{
double dx1 = line1.EndPoint.X - line1.BeginPoint.X;
double dy1 = line1.EndPoint.Y - line1.BeginPoint.Y;
double dx2 = line2.EndPoint.X - line2.BeginPoint.X;
double dy2 = line2.EndPoint.Y - line2.BeginPoint.Y;
double noem = dx1*dy2 - dy1*dx2;
double u, v;
if (Math.Abs(noem) > 0.0)
{
u = (dx2*(line1.BeginPoint.Y - line2.BeginPoint.Y) - dy2*(line1.BeginPoint.X - line2.BeginPoint.X))/
noem;
if ((u >= 0.0) && (u <= 1.0))
{
v = (dx1*(line1.BeginPoint.Y - line2.BeginPoint.Y) -
dy1*(line1.BeginPoint.X - line2.BeginPoint.X))/
noem;
if ((v >= 0.0) && (v <= 1.0))
{
return new GeometryPoint(line1.BeginPoint.X + u*dx1, 0, line1.BeginPoint.Y + u*dy1);
}
}
}
return new GeometryPoint(double.NaN, double.NaN, 0);
}
else
{
return new GeometryPoint(double.NaN, double.NaN, double.NaN);
}
}
public static GeometryPoint GetIntersectPoint(GeometryPoint p1, GeometryPoint p2, GeometryPoint p3, GeometryPoint p4)
{
return IntersectionPoint(p1, p2, p3, p4);
}
public static double GetSlope(GeometryPoint point1, GeometryPoint point2)
{
return (Math.Abs(point1.X - point2.X) > 0) ? (point1.Z - point2.Z)/(point1.X - point2.X) : 0;
}
public static GeometryPoint IntersectionPoint(this Line thisLine, Line otherLine)
{
return IntersectionPoint(thisLine.BeginPoint, thisLine.EndPoint, otherLine.BeginPoint, otherLine.EndPoint);
}
///
/// Calculates the intersection point between two line (segments)
///
/// Starting point line segment 1
/// End point line segment 1
/// Starting point line segment 2
/// End point line segment 2
/// The intersection point
///
/// The point can contain infinite values if the lines are parallel
///
public static GeometryPoint IntersectionPoint(GeometryPoint p1, GeometryPoint p2, GeometryPoint p3, GeometryPoint p4)
{
if (p1 == null)
{
throw new ArgumentNullException("p1", "The point is null or empty which is not allowed");
}
if (p2 == null)
{
throw new ArgumentNullException("p2", "The point is null or empty which is not allowed");
}
if (p3 == null)
{
throw new ArgumentNullException("p3", "The point is null or empty which is not allowed");
}
if (p4 == null)
{
throw new ArgumentNullException("p4", "The point is null or empty which is not allowed");
}
if (p1.LocationEquals(p3))
{
return new GeometryPoint
{
X = p1.X, Z = p1.Z
};
}
if (p1.LocationEquals(p4))
{
return new GeometryPoint
{
X = p1.X, Z = p1.Z
};
}
if (p2.LocationEquals(p3))
{
return new GeometryPoint
{
X = p2.X, Z = p2.Z
};
}
if (p2.LocationEquals(p4))
{
return new GeometryPoint
{
X = p4.X, Z = p4.Z
};
}
var t = (p1.X*(p4.Z - p3.Z) + p3.X*(p1.Z - p4.Z) + p4.X*(p3.Z - p1.Z))/((p2.Z - p1.Z)*(p4.X - p3.X) -
(p2.X - p1.X)*(p4.Z - p3.Z));
var s = (p1.X*(p3.Z - p2.Z) + p2.X*(p2.Z - p3.Z) + p3.X*(p2.Z - p1.Z))/((p4.Z - p3.Z)*(p2.X - p1.X) -
(p4.X - p3.X)*(p2.Z - p1.Z));
var xi = p1.X + t*(p2.X - p1.X);
var zi = p3.Z + s*(p4.Z - p3.Z);
var pi = new GeometryPoint
{
X = xi, Z = zi
};
return pi;
}
///
/// Determines whether the given point is above, beneath or on the surfaceline.
///
///
///
public static PLLinePointPositionXzType PositionXzOfPointRelatedToPLLine(GeometryPoint point, IList points)
{
// if point is out of scope of the surface line, return beyond
if ((point.X < points[0].X) || (point.X > points[points.Count - 1].X))
{
return PLLinePointPositionXzType.BeyondPLLine;
}
double z = ZFromX(point.X, points);
if (Math.Abs(point.Z - z) < GeometryPoint.Precision)
{
return PLLinePointPositionXzType.OnPLLine;
}
else
{
if (point.Z > z)
{
return PLLinePointPositionXzType.AbovePLLine;
}
else
{
return PLLinePointPositionXzType.BelowPLLine;
}
}
}
public static double ZFromX(double X, IList points)
{
return YZFromX(X, points, PointType.XZ);
}
///
/// Check if all the points are in ascending X order
///
///
///
public static bool IsXAscending(IList points)
{
bool isAscending = true;
for (int pointIndex = 0; pointIndex < points.Count - 1; pointIndex++)
{
if (points[pointIndex + 1].X < points[pointIndex].X)
{
isAscending = false;
break;
}
}
return isAscending;
}
///
/// Deletes the coinsiding points.
///
///
/// The tolerance.
public static void DeleteCoinsidingPoints(IList points, double tolerance)
{
// First build a list of indices of points that have to be removed (from end of list to start)
var indicesToDelete = new List();
for (int pointIndex = points.Count - 1; pointIndex > 0; pointIndex--)
{
for (int i = pointIndex - 1; i >= 0; i--)
{
if (points[pointIndex].LocationEquals(points[i], tolerance))
{
indicesToDelete.Add(i);
}
}
}
// Remove duplicate points beginning from the end
for (int index = 0; index < indicesToDelete.Count; index++)
{
points.RemoveAt(indicesToDelete[index]);
}
}
///
/// Gets point at X if present within tolerance; creates one there if not present.
///
///
public static GeometryPoint EnsurePointAtX(double X, IList points, double tolerance)
{
var point = GetPointAtX(X, points, tolerance);
if (point == null)
{
point = InsertPointAtX(X, points);
}
return point;
}
///
/// Gets point at X if present within tolerance or null if not present.
///
///
public static GeometryPoint GetPointAtX(double X, IList points, double tolerance)
{
return (points.Cast().Where(point => Math.Abs(point.X - X) < tolerance)).FirstOrDefault();
}
///
/// Gets the points in the segment between starting x and ending x
///
///
///
///
///
///
///
public static IEnumerable GetPointSegmentBetween(double startX, double endX, IList points)
{
if (endX < startX)
{
throw new ArgumentException("End value is smaller then the start value");
}
return PointsOrderdByX(points).Where(
point => point != null && (point.X > startX && point.X < endX)).OrderBy(point => point.X);
}
///
/// Removes the index of all points of GeometryPointString after the start index.
///
///
///
public static void RemoveAllPointsOfGeometryPointStringAfterIndex(GeometryPointString pointString, int startIndex)
{
for (int pointIndex = pointString.Points.Count - 1; pointIndex > startIndex; pointIndex--)
{
pointString.Points.RemoveAt(pointIndex);
}
}
private static double YZFromX(double X, IList points, PointType pointType)
{
if (!IsXAscending(points))
{
throw new System.Exception("Interpolation only possible with ascending polyline");
}
double valueYZ = 0.0;
// Less then first X
if (X <= points[0].X)
{
valueYZ = pointType == PointType.XZ ? points[0].Z : points[0].Y;
}
// More then last X
else if (X >= points[points.Count - 1].X)
{
valueYZ = pointType == PointType.XZ ? points[points.Count - 1].Z : points[points.Count - 1].Y;
}
else
{
// X is inside boundaries
for (int pointIndex = 0; pointIndex < points.Count - 1; pointIndex++)
{
if ((X > points[pointIndex].X) && (X <= points[pointIndex + 1].X))
{
// interpolate in this section
double fractionX = (X - points[pointIndex].X)/(points[pointIndex + 1].X - points[pointIndex].X);
if (pointType == PointType.XZ)
{
valueYZ = points[pointIndex].Z + fractionX*(points[pointIndex + 1].Z - points[pointIndex].Z);
}
else
{
valueYZ = points[pointIndex].Y + fractionX*(points[pointIndex + 1].Y - points[pointIndex].Y);
}
break;
}
}
}
return valueYZ;
}
private static GeometryPoint InsertPointAtX(double X, IList points)
{
var newPoint = new GeometryPoint();
newPoint.X = X;
try
{
GeometryPoint pointAfter = (from GeometryPoint point in points
where point.X > X
select point).First();
points.Insert(points.IndexOf(pointAfter), newPoint);
}
catch
{
points.Add(newPoint);
}
return newPoint;
}
private static IEnumerable PointsOrderdByX(IList points)
{
return points.OrderBy(p => p.X);
}
}
}