Skip to content

Commit

Permalink
Merge pull request #419 from NishantChandrashekhar/dev/nich/Point.InP…
Browse files Browse the repository at this point in the history 
…olygonBug

Dev/nich/point.in polygon bug
  • Loading branch information
@d3ssy
d3ssy authored Aug 7, 2023
2 parents 76819e0 + 88b5b6c commit 8eb033d
Show file tree
Hide file tree
Showing 2 changed files with 300 additions and 51 deletions.
241 changes: 237 additions & 4 deletions src/GShark.Test.XUnit/Geometry/Point3Tests.cs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
using System.Collections.Generic;
using System.Linq;
using FluentAssertions;
using GShark.Geometry;
using Xunit;
Expand Down Expand Up @@ -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
Expand All @@ -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);
}

/// <summary>
/// x is the point of interest
///
/// (0,2)___(1,2)
/// | |
/// | |____(2,1) x(3,1)
/// | (1,1) |
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}

/// <summary>
/// x is the point of interest
///
/// (0,2)___(1,2)
/// | |
/// | |____x(2,1)
/// | (1,1) |
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}

/// <summary>
/// x is the point of interest (0.5, 1)
///
/// (0,2)___(1,2)
/// | |
/// | x |____(2,1)
/// | (1,1)|
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}
/// <summary>
/// x is the point of interest (-1, 1)
///
/// (0,2)___(1,2)
/// | |
/// x(-1,1) | |____(2,1)
/// | (1,1)|
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}
/// <summary>
/// x is the point of interest (0.5, 1.5)
///
/// (0,2)___(1,2)
/// | x |
/// | |____(2,1)
/// | (1,1)|
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}
/// <summary>
/// x is the point of interest (-1, 1)
///
/// x(-1, 1) (0,2)___(1,2)
/// | |
/// | |____(2,1)
/// | (1,1)|
/// |________|
/// (0,0) (2,0)
///
/// </summary>
[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);
}
}
}
110 changes: 63 additions & 47 deletions src/GShark/Geometry/Point3.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -657,67 +657,83 @@ public bool IsOnLine(Line line, double tolerance = GSharkMath.MaxTolerance)
/// <para>See https://stackoverflow.com/a/63436180</para>
/// </summary>
/// <param name="polygon">The polygon to test against.</param>
/// <returns>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>
/// <returns>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>
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<double, double, double, bool> 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<double, double, double, bool> 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;
}

}
}

0 comments on commit 8eb033d

Please sign in to comment.