Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Dev/nich/point.in polygon bug #419

Merged
merged 4 commits into from
Aug 7, 2023
Merged

Dev/nich/point.in polygon bug #419

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
28c3e36
WIP
NishantChandrashekhar Jul 6, 2023
2f35635
Implemented
NishantChandrashekhar Jul 7, 2023
3e78484
commented code
NishantChandrashekhar Jul 8, 2023
88b5b6c
reverting back an unnecessary change
NishantChandrashekhar Jul 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
241 changes: 237 additions & 4 deletions src/GShark.Test.XUnit/Geometry/Point3Tests.cs
View file
Open in desktop
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
Open in desktop
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;
}

}
}

Loading