diff --git a/src/GShark.Test.XUnit/Geometry/Point3Tests.cs b/src/GShark.Test.XUnit/Geometry/Point3Tests.cs
index 0190c21f..1ea3622f 100644
--- a/src/GShark.Test.XUnit/Geometry/Point3Tests.cs
+++ b/src/GShark.Test.XUnit/Geometry/Point3Tests.cs
@@ -1,4 +1,5 @@
using System.Collections.Generic;
+using System.Linq;
using FluentAssertions;
using GShark.Geometry;
using Xunit;
@@ -197,10 +198,10 @@ public void It_Returns_Whether_A_Point_Is_Inside_Outside_Or_Coincident_With_A_Po
var testPolygon = new GShark.Geometry.Polygon(new Point3[] {
new Point3(-27.1829592472304,-12.3049979552339,59.1652925745504),
- new Point3(-40.9982936339814,20.739935073677,34.4162859330554),
- new Point3(19.082346606075,37.5838503530052,21.8010335204645),
- new Point3(49.9258838416119,19.3255605082373,35.4755819371661),
- new Point3(34.18282837764,-33.7777754487285,75.2473319096853)
+ new Point3(-40.9982936339814,20.739935073677 ,34.4162859330554),
+ new Point3(19.082346606075 ,37.5838503530052 ,21.8010335204645),
+ new Point3(49.9258838416119 ,19.3255605082373 ,35.4755819371661),
+ new Point3(34.18282837764 ,-33.7777754487285,75.2473319096853)
});
//Act
@@ -213,5 +214,237 @@ public void It_Returns_Whether_A_Point_Is_Inside_Outside_Or_Coincident_With_A_Po
pointInside.Should().Be(1);
pointCoincident.Should().Be(0);
}
+
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Inside_Triangle()
+ {
+ // Arrange
+ var testPolygon2 = new GShark.Geometry.Polygon(new Point3[] {
+ new Point3(10.00,0,0.00),
+ new Point3(-5.00,8.660,0.00),
+ new Point3(-5.00,-8.660,0.00),
+
+ });
+ var p3 = new Point3(5, 0, 0);
+
+ // Act
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside2.Should().Be(1);
+ }
+
+ ///
+ /// x is the point of interest
+ ///
+ /// (0,2)___(1,2)
+ /// | |
+ /// | |____(2,1) x(3,1)
+ /// | (1,1) |
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Outside_Polygon_Right_Side()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(3, 1, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(-1);
+ pointIsInside2.Should().Be(-1);
+ }
+
+ ///
+ /// x is the point of interest
+ ///
+ /// (0,2)___(1,2)
+ /// | |
+ /// | |____x(2,1)
+ /// | (1,1) |
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_On_Polygon_Right_Side()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(2, 1, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(0);
+ pointIsInside2.Should().Be(0);
+ }
+
+ ///
+ /// x is the point of interest (0.5, 1)
+ ///
+ /// (0,2)___(1,2)
+ /// | |
+ /// | x |____(2,1)
+ /// | (1,1)|
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Inside_Polygon()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(0.5, 1, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(1);
+ pointIsInside2.Should().Be(1);
+ }
+ ///
+ /// x is the point of interest (-1, 1)
+ ///
+ /// (0,2)___(1,2)
+ /// | |
+ /// x(-1,1) | |____(2,1)
+ /// | (1,1)|
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Outside_Polygon_Left_Side()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(-1, 1, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(-1);
+ pointIsInside2.Should().Be(-1);
+ }
+ ///
+ /// x is the point of interest (0.5, 1.5)
+ ///
+ /// (0,2)___(1,2)
+ /// | x |
+ /// | |____(2,1)
+ /// | (1,1)|
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Inside_Polygon_Top()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(0.5, 1.5, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(1);
+ pointIsInside2.Should().Be(1);
+ }
+ ///
+ /// x is the point of interest (-1, 1)
+ ///
+ /// x(-1, 1) (0,2)___(1,2)
+ /// | |
+ /// | |____(2,1)
+ /// | (1,1)|
+ /// |________|
+ /// (0,0) (2,0)
+ ///
+ ///
+ [Fact]
+ public void It_Returns_Whether_A_Point_Is_Outside_Polygon_Top()
+ {
+ // Arrange
+ Point3[] pts = new Point3[] {
+ new Point3(0,0,0),
+ new Point3(2,0,0),
+ new Point3(2,1,0),
+ new Point3(1,1,0),
+ new Point3(1,2,0),
+ new Point3(0,2,0),
+ };
+ var testPolygon1 = new GShark.Geometry.Polygon(pts); // CCW
+ var testPolygon2 = new GShark.Geometry.Polygon(pts.Reverse().ToList()); // CW
+ var p3 = new Point3(-1, 1, 0);
+
+ // Act
+ var pointIsInside1 = p3.InPolygon(testPolygon1);
+ var pointIsInside2 = p3.InPolygon(testPolygon2);
+
+ // Assert
+ pointIsInside1.Should().Be(-1);
+ pointIsInside2.Should().Be(-1);
+ }
}
}
diff --git a/src/GShark/Geometry/Point3.cs b/src/GShark/Geometry/Point3.cs
index 12ce29e3..04051619 100644
--- a/src/GShark/Geometry/Point3.cs
+++ b/src/GShark/Geometry/Point3.cs
@@ -657,67 +657,83 @@ public bool IsOnLine(Line line, double tolerance = GSharkMath.MaxTolerance)
/// See https://stackoverflow.com/a/63436180
///
/// The polygon to test against.
- /// Returns -1 if point is outside the polygon, 0 if it is coincident with a polygon edge, or 1 if it is inside the polygon.
+ /// Returns -1 if point is outside the polygon, 0 if it is coincident with a polygon edge, or 1 if it is inside the polygon.
public int InPolygon(Polygon polygon)
{
- //check if point lies on polygon plane, else return
- var polygonPlane = polygon.Plane;
-
- if (!this.IsOnPlane(polygonPlane)) return -1;
-
- //translate polygon and point to XY plane for 2d calculations to account for rotated polygons and 3d points
+ // Tests whether a value is isBetween two other values.
+ Func isValueBetween = (p, a, b) =>
+ (p > a) &&
+ (p < b) ||
+ (p < a) &&
+ (p > b);
+
+ Plane polygonPlane = polygon.Plane;
+
+ // Check if point lies on the polygon plane.
+ if (!this.IsOnPlane(polygonPlane))
+ return -1;
+
+ // Transform polygon and point to XY plane for 2d calculations to account for rotated polygons and 3d points
var xForm = Core.Transform.PlaneToPlane(polygonPlane, Plane.PlaneXY);
var polygonOriented = polygon.Transform(xForm);
var pointOriented = this.Transform(xForm);
- //tests whether a value is isBetween two other values
- Func isValueBetween = (p, a, b) =>
- ((p - a) >= double.Epsilon) &&
- ((p - b) <= double.Epsilon) ||
- ((p - a) <= double.Epsilon) &&
- ((p - b) >= double.Epsilon);
-
- bool inside = false;
+ // The winding number
+ int wn = 0;
- for (int i = polygonOriented.ControlPointLocations.Count - 1, j = 0; j < polygonOriented.ControlPointLocations.Count; i = j, j++)
+ // Looping through all the edges.
+ int n = polygonOriented.ControlPointLocations.Count;
+ for (int i = 0; i < n - 1; i++)
{
- Point3 A = polygonOriented.ControlPointLocations[i];
- Point3 B = polygonOriented.ControlPointLocations[j];
+ Line edge = new Line(polygonOriented.ControlPointLocations[i], polygonOriented.ControlPointLocations[i + 1]);
+
+ // If point on edge, we return.
+ if (pointOriented.IsOnLine(edge))
+ return 0;
+
+ double Py = pointOriented.Y;
+ double x1 = polygonOriented.ControlPointLocations[i].X;
+ double y1 = polygonOriented.ControlPointLocations[i].Y;
+ double x2 = polygonOriented.ControlPointLocations[i + 1].X;
+ double y2 = polygonOriented.ControlPointLocations[i + 1].Y;
+
+ Point3 start = new Point3(x1, y1, 0);
+ Point3 end = new Point3(x2, y2, 0);
- // corner cases
- if (
- (Math.Abs(pointOriented.X - A.X) <= double.Epsilon) &&
- (Math.Abs(pointOriented.Y - A.Y) <= double.Epsilon) ||
- (Math.Abs(pointOriented.X - B.X) <= double.Epsilon) &&
- (Math.Abs(pointOriented.Y - B.Y) <= double.Epsilon)) return 0;
+ // Corner case when an edge is horizontal, we just treat it like a edge that is a bit lower
+ // than the point (hence doesn't intersect).
+ if (y2 - y1 == 0)
+ {
+ continue;
+ }
- if (
- Math.Abs(A.Y - B.Y) <= double.Epsilon &&
- Math.Abs(pointOriented.Y - A.Y) <= double.Epsilon &&
- isValueBetween(pointOriented.X, A.X, B.X)) return 0;
-
- if (isValueBetween(pointOriented.Y, A.Y, B.Y))
+ // Calculating the intersection of the ray casted with the edge
+ double intX = (Py * (x2 - x1) + x1 * y2 - x2 * y1) / (y2 - y1);
+ Point3 intersectionPoint = new Point3(intX, Py, 0);
+
+ // Incase the ray intersects one of the vertices, we see if the previous vertex and the next vertex
+ // are on the same side or the opposite side of the horizontal ray. If on same side we just carry on
+ // else we increament the winding number. Since the end vertex of the edge will be the start vertex
+ // of the next edge, we don't have to do this explicitly for both start and end vertex.
+ if (intersectionPoint.Equals(end))
{
- // if P inside the vertical range
- // filter out "ray pass vertex" problem by treating the line a little lower
- if (
- Math.Abs(pointOriented.Y - A.Y) <= double.Epsilon &&
- (B.Y - A.Y) >= double.Epsilon ||
- Math.Abs(pointOriented.Y - B.Y) <= double.Epsilon &&
- (A.Y - B.Y) <= double.Epsilon) continue;
-
- // calc cross product `PA X PB`, P lays on left side of AB if c > 0
- double c = (A.X - pointOriented.X) * (B.Y - pointOriented.Y) - (B.X - pointOriented.X) * (A.Y - pointOriented.Y);
-
- if (c > 0 && c < GSharkMath.MinTolerance ) return 0;
-
- if ((A.Y < B.Y) == (c > 0))
- inside = !inside;
+ Point3 previous = start; // Nothing but the start vertex of this edge.
+ Point3 next = polygonOriented.ControlPointLocations[(i + 2) % n]; // The end vertex of the next edge.
+ if (isValueBetween(Py, previous.Y, next.Y))
+ wn++;
+ continue;
}
+
+ // The general case when the ray intersects the edge. Since the intersection point calculated before
+ // considered that the line was NOT a half-line, we need to add additional check to see if the inter-
+ // section point is to the right of this point.
+ if (intersectionPoint.IsOnLine(edge) && intersectionPoint.X > pointOriented.X)
+ wn++;
}
- return inside ? 1 : -1;
+
+ // Applying the even-odd rule.
+ return wn % 2 == 0 ? -1 : 1;
}
}
}
-