C# Point in polygon
See this it's in c++ and can be done in c# in a same way.
for convex polygon is too easy:
If the polygon is convex then one can
consider the polygon as a "path" from
the first vertex. A point is on the
interior of this polygons if it is
always on the same side of all the
line segments making up the path.Given a line segment between P0
(x0,y0) and P1 (x1,y1), another point
P (x,y) has the following relationship
to the line segment. Compute (y - y0)
(x1 - x0) - (x - x0) (y1 - y0)if it is less than 0 then P is to the
right of the line segment, if greater
than 0 it is to the left, if equal to
0 then it lies on the line segment.
Here is its code in c#, I didn't check edge cases.
public static bool IsInPolygon(Point[] poly, Point point)
{
var coef = poly.Skip(1).Select((p, i) =>
(point.Y - poly[i].Y)*(p.X - poly[i].X)
- (point.X - poly[i].X) * (p.Y - poly[i].Y))
.ToList();
if (coef.Any(p => p == 0))
return true;
for (int i = 1; i < coef.Count(); i++)
{
if (coef[i] * coef[i - 1] < 0)
return false;
}
return true;
}
I test it with simple rectangle works fine:
Point[] pts = new Point[] { new Point { X = 1, Y = 1 },
new Point { X = 1, Y = 3 },
new Point { X = 3, Y = 3 },
new Point { X = 3, Y = 1 } };
IsInPolygon(pts, new Point { X = 2, Y = 2 }); ==> true
IsInPolygon(pts, new Point { X = 1, Y = 2 }); ==> true
IsInPolygon(pts, new Point { X = 0, Y = 2 }); ==> false
Explanation on the linq query:
poly.Skip(1) ==> creates a new list started from position 1 of the poly list and then by (point.Y - poly[i].Y)*(p.X - poly[i].X) - (point.X - poly[i].X) * (p.Y - poly[i].Y) we'll going to calculate the direction (which mentioned in referenced paragraph).
similar example (with another operation):
lst = 2,4,8,12,7,19
lst.Skip(1) ==> 4,8,12,7,19
lst.Skip(1).Select((p,i)=>p-lst[i]) ==> 2,4,4,-5,12
Is point inside polygon?
Solved it this way:
public static bool IsPointInPolygon4(PointF[] polygon, PointF testPoint)
{
bool result = false;
int j = polygon.Count() - 1;
for (int i = 0; i < polygon.Count(); i++)
{
if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y)
{
if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X)
{
result = !result;
}
}
j = i;
}
return result;
}
Find point in polygon
Draw an horizontal line at some ordinate between the extreme values (to be sure to cross the polygon). Compute all intersections of this line with the polygon outline. This is done in a simple loop on the polygon edges, where you detect the changes of side.
You will find an even number of intersections. Sort them by increasing abscissa and pick any point between an even and an odd intersection.
This takes time O(N), which is optimal.
As you can check, the point-in-polygon test by ray casting reports true.
Erratum:
In the worst case, there can be O(N) intersections, and the sorting will take O(N Log(N)) operations. In practice, this is rarely met.
Side note:
If you have the extra requirement that the interior point should be at a sufficient distance from the edges, you can perform the offsetting of the polygon (inner side), and pick a point inside the offset using the same procedure. Offset computation is by no means trivial, though.
The possible offset distances are bounded, but I have no idea how this bound can be found. It will correspond to "deepest" point(s).
How to add a point in somewhere between two points of a complete polygon
public static double GetAngleBetweenPoints(Point p1, Point p2)
{
double xDiff = p2.X - p1.X;
double yDiff = p1.Y - p2.Y;
double Angle = Math.Atan2(yDiff, xDiff) * (180 / Math.PI);
if (Angle < 0.0)
Angle = Angle + 360.0;
return Angle;
}
public static GenericList<Point> PointAddingByAngle(GenericList<Point> LinePoints2)
{
List<int> AnglefromCenterPointsofShapePoint = new List<int>();
List<Point> Linepoints3 = new List<Point>(LinePoints2);
List<Point> LinePoints = Linepoints3.Distinct().ToList();
Point CenterPoint = new Point();
CenterPoint.X = PolygonCenterX;
CenterPoint.Y = PolygonCenterY;
for (int index = 0; index < LinePoints.Count; index++)
{
AnglefromCenterPointsofShapePoint.Add((int) Math.Round(GetAngleBetweenPoints(CenterPoint, LinePoints[index])));
}
AnglefromCenterPointsofShapePoint.Sort();
int currentPointAngle = 0;
for (int AngleListIndex = 0; AngleListIndex < AnglefromCenterPointsofShapePoint.Count; AngleListIndex++)
{
currentPointAngle = AnglefromCenterPointsofShapePoint[AngleListIndex];
for (int index = 0; index < LinePoints.Count; index++)
{
if ((Math.Round(GetAngleBetweenPoints(CenterPoint, LinePoints[index]))) == currentPointAngle)
{
SortedPoints.Add(LinePoints[index]);
break;
}
}
}
SortedPoints.Add(SortedPoints[0]);
return SortedPoints;
}
private void pictureBox_MouseDoubleClick(object sender, MouseEventArgs e)
{
if (PathPoints.Count > 2)
{
PathPoints.Add(PathPoints[0]);
IsPathComplete = true;
IsPointSelected = false;
StartingPoint = PathPoints[0];
EndingPoint = PathPoints[PathPoints.Count - 1];
}
if (IsPathComplete && NewPointAdd)
{
System.Drawing.Point Centerpoint = new System.Drawing.Point();
PathPoints.Add(e.Location);
MinMaxFinder(PathPoints);
Centerpoint.X = PolygonCenterX;
Centerpoint.Y = PolygonCenterY;
lblAngle.Text = Convert.ToString(Math.Round(ShapeDirectory.GetAngleBetweenPoints(Centerpoint, e.Location)));
PathPoints3.Clear();
PathPoints3 = ShapeDirectory.PointAddingByAngle(PathPoints);
PathPoints.Clear();
foreach (System.Drawing.Point p in PathPoints3)
{
PathPoints.Add(p);
}
PathPoints.Add(PathPoints3[0]);
NewPointAddDraw = true;
}
}
Point in Polygon using Winding Number
The answer that I used as my starting point was provided by Manuel Castro in the following SO thread Geo Fencing - point inside/outside polygon
:
public static bool PointInPolygon(LatLong p, List<LatLong> poly)
{
int n = poly.Count();
poly.Add(new LatLong { Lat = poly[0].Lat, Lon = poly[0].Lon });
LatLong[] v = poly.ToArray();
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (int i = 0; i < n; i++)
{ // edge from V[i] to V[i+1]
if (v[i].Lat <= p.Lat)
{ // start y <= P.y
if (v[i + 1].Lat > p.Lat) // an upward crossing
if (isLeft(v[i], v[i + 1], p) > 0) // P left of edge
++wn; // have a valid up intersect
}
else
{ // start y > P.y (no test needed)
if (v[i + 1].Lat <= p.Lat) // a downward crossing
if (isLeft(v[i], v[i + 1], p) < 0) // P right of edge
--wn; // have a valid down intersect
}
}
if (wn != 0)
return true;
else
return false;
}
I proceeded to write xUnit tests around an implementation that started out using the exact logic expressed in the above code. The following are a bit overkill, but I preferred more tests to ensure that refactoring would not create issues. Using inline data in xUnit theories is so easy, eh, why not. With all the tests green, I was then able to refactor to my heart's content:
public class PolygonTests
{
public class GivenLine : PolygonTests
{
private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
{
new GeographicalPoint(1, 1),
new GeographicalPoint(10, 1)
});
public class AndPointIsAnywhere : GivenLine
{
[Theory]
[InlineData(5, 1)]
[InlineData(-1, -1)]
[InlineData(11, 11)]
public void WhenAskingContainsLocation_ThenItShouldReturnFalse(double latitude, double longitude)
{
GeographicalPoint point = new GeographicalPoint(latitude, longitude);
bool actual = _polygon.Contains(point);
actual.Should().BeFalse();
}
}
}
public class GivenTriangle : PolygonTests
{
private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
{
new GeographicalPoint(1, 1),
new GeographicalPoint(10, 1),
new GeographicalPoint(10, 10)
});
public class AndPointWithinTriangle : GivenTriangle
{
private readonly GeographicalPoint _point = new GeographicalPoint(6, 4);
[Fact]
public void WhenAskingContainsLocation_ThenItShouldReturnTrue()
{
bool actual = _polygon.Contains(_point);
actual.Should().BeTrue();
}
}
public class AndPointOutsideOfTriangle : GivenTriangle
{
private readonly GeographicalPoint _point = new GeographicalPoint(5, 5.0001d);
[Fact]
public void WhenAskingContainsLocation_ThenItShouldReturnFalse()
{
bool actual = _polygon.Contains(_point);
actual.Should().BeFalse();
}
}
}
public class GivenComplexPolygon : PolygonTests
{
private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
{
new GeographicalPoint(1, 1),
new GeographicalPoint(5, 1),
new GeographicalPoint(8, 4),
new GeographicalPoint(3, 4),
new GeographicalPoint(8, 9),
new GeographicalPoint(1, 9)
});
[Theory]
[InlineData(5, 0, false)]
[InlineData(5, 0.999d, false)]
[InlineData(5, 1, true)]
[InlineData(5, 2, true)]
[InlineData(5, 3, true)]
[InlineData(5, 4, false)]
[InlineData(5, 5, false)]
[InlineData(5, 5.999d, false)]
[InlineData(5, 6, true)]
[InlineData(5, 7, true)]
[InlineData(5, 8, true)]
[InlineData(5, 9, false)]
[InlineData(5, 10, false)]
[InlineData(0, 5, false)]
[InlineData(0.999d, 5, false)]
[InlineData(1, 5, true)]
[InlineData(2, 5, true)]
[InlineData(3, 5, true)]
[InlineData(4.001d, 5, false)]
//[InlineData(5, 5, false)] -- duplicate
[InlineData(6, 5, false)]
[InlineData(7, 5, false)]
[InlineData(8, 5, false)]
[InlineData(9, 5, false)]
[InlineData(10, 5, false)]
public void WhenAskingContainsLocation_ThenItShouldReturnCorrectAnswer(double latitude, double longitude, bool expected)
{
GeographicalPoint point = new GeographicalPoint(latitude, longitude);
bool actual = _polygon.Contains(point);
actual.Should().Be(expected);
}
}
}
This allowed me to refactor the original code to the following:
public interface IPolygon
{
bool Contains(GeographicalPoint location);
}
public class Polygon : IPolygon
{
private readonly List<GeographicalPoint> _points;
public Polygon(List<GeographicalPoint> points)
{
_points = points;
}
public bool Contains(GeographicalPoint location)
{
GeographicalPoint[] polygonPointsWithClosure = PolygonPointsWithClosure();
int windingNumber = 0;
for (int pointIndex = 0; pointIndex < polygonPointsWithClosure.Length - 1; pointIndex++)
{
Edge edge = new Edge(polygonPointsWithClosure[pointIndex], polygonPointsWithClosure[pointIndex + 1]);
windingNumber += AscendingIntersection(location, edge);
windingNumber -= DescendingIntersection(location, edge);
}
return windingNumber != 0;
}
private GeographicalPoint[] PolygonPointsWithClosure()
{
// _points should remain immutable, thus creation of a closed point set (starting point repeated)
return new List<GeographicalPoint>(_points)
{
new GeographicalPoint(_points[0].Latitude, _points[0].Longitude)
}.ToArray();
}
private static int AscendingIntersection(GeographicalPoint location, Edge edge)
{
if (!edge.AscendingRelativeTo(location)) { return 0; }
if (!edge.LocationInRange(location, Orientation.Ascending)) { return 0; }
return Wind(location, edge, Position.Left);
}
private static int DescendingIntersection(GeographicalPoint location, Edge edge)
{
if (edge.AscendingRelativeTo(location)) { return 0; }
if (!edge.LocationInRange(location, Orientation.Descending)) { return 0; }
return Wind(location, edge, Position.Right);
}
private static int Wind(GeographicalPoint location, Edge edge, Position position)
{
if (edge.RelativePositionOf(location) != position) { return 0; }
return 1;
}
private class Edge
{
private readonly GeographicalPoint _startPoint;
private readonly GeographicalPoint _endPoint;
public Edge(GeographicalPoint startPoint, GeographicalPoint endPoint)
{
_startPoint = startPoint;
_endPoint = endPoint;
}
public Position RelativePositionOf(GeographicalPoint location)
{
double positionCalculation =
(_endPoint.Longitude - _startPoint.Longitude) * (location.Latitude - _startPoint.Latitude) -
(location.Longitude - _startPoint.Longitude) * (_endPoint.Latitude - _startPoint.Latitude);
if (positionCalculation > 0) { return Position.Left; }
if (positionCalculation < 0) { return Position.Right; }
return Position.Center;
}
public bool AscendingRelativeTo(GeographicalPoint location)
{
return _startPoint.Latitude <= location.Latitude;
}
public bool LocationInRange(GeographicalPoint location, Orientation orientation)
{
if (orientation == Orientation.Ascending) return _endPoint.Latitude > location.Latitude;
return _endPoint.Latitude <= location.Latitude;
}
}
private enum Position
{
Left,
Right,
Center
}
private enum Orientation
{
Ascending,
Descending
}
}
public class GeographicalPoint
{
public double Longitude { get; set; }
public double Latitude { get; set; }
public GeographicalPoint(double latitude, double longitude)
{
Latitude = latitude;
Longitude = longitude;
}
}
The original code, of course, is far less verbose. However, it:
- is procedural;
- uses variable names that do not reveal intent;
- is mutable;
- has a cyclomatic complexity of 12.
The refactored code:
- passes the tests;
- reveals intention;
- contains no duplication;
- contains the fewest elements (given 1, 2 and 3, above)
And:
- is object oriented;
- does not use if except for guard clauses;
- is immutable;
- hides its private data;
- has complete test coverage;
- has one method with a cyclomatic complexity of 3, while most of the methods are at 1, and a few of them are at 2.
Now, all that said, I'm not saying that there are no additional refactors that might be suggested, or that the above refactor approaches perfection. However, I think that in terms of implementing the Winding Number algorithm for determining whether a point is in a polygon and really understanding the algorithm that this is helpful.
I hope that this helps some who, like me, had some difficulty wrapping their heads around it.
Cheers
Point inside compound polygon
With this line
p1.Y < testY != p2.Y < testY
you only count segments that approach the query line from the bottom (or cross it). And this is exactly what you need to do for the arcs.
For this, it may be desirable to split the arc into monotone parts (with respect to the y-axis). In your example, the lower arc already is monotone. The upper arc should be split into two (along a vertical line through the its center). Then you have minY and maxY for each segment and you can apply the above formula:
minY < testY != maxY < testY
Alternatively, you could check if the intersection is at the arc end (equal y-coordinate) and evaluate if the arc continues downwards based on the angle and arc direction. But depending on the implementation, this might have some numerical stability issues. The first option (with splitting into monotone parts) is probably easier to implement. And it generalizes to other primitives.
How to determine if a Point is inside of a 2D Polygon that has Arcs
I asked this question in the AutoDesk .NET Forum and recieved this solution:
The Answer
How can I determine whether a 2D Point is within a Polygon?
For graphics, I'd rather not prefer integers. Many systems use integers for UI painting (pixels are ints after all), but macOS, for example, uses float for everything. macOS only knows points and a point can translate to one pixel, but depending on monitor resolution, it might translate to something else. On retina screens half a point (0.5/0.5) is pixel. Still, I never noticed that macOS UIs are significantly slower than other UIs. After all, 3D APIs (OpenGL or Direct3D) also work with floats and modern graphics libraries very often take advantage of GPU acceleration.
Now you said speed is your main concern, okay, let's go for speed. Before you run any sophisticated algorithm, first do a simple test. Create an axis aligned bounding box around your polygon. This is very easy, fast and can already save you a lot of calculations. How does that work? Iterate over all points of the polygon and find the min/max values of X and Y.
E.g. you have the points (9/1), (4/3), (2/7), (8/2), (3/6). This means Xmin is 2, Xmax is 9, Ymin is 1 and Ymax is 7. A point outside of the rectangle with the two edges (2/1) and (9/7) cannot be within the polygon.
// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
This is the first test to run for any point. As you can see, this test is ultra fast but it's also very coarse. To handle points that are within the bounding rectangle, we need a more sophisticated algorithm. There are a couple of ways how this can be calculated. Which method works also depends on whether the polygon can have holes or will always be solid. Here are examples of solid ones (one convex, one concave):
And here's one with a hole:
The green one has a hole in the middle!
The easiest algorithm, that can handle all three cases above and is still pretty fast is named ray casting. The idea of the algorithm is pretty simple: Draw a virtual ray from anywhere outside the polygon to your point and count how often it hits a side of the polygon. If the number of hits is even, it's outside of the polygon, if it's odd, it's inside.
The winding number algorithm would be an alternative, it is more accurate for points being very close to a polygon line but it's also much slower. Ray casting may fail for points too close to a polygon side because of limited floating point precision and rounding issues, but in reality that is hardly a problem, as if a point lies that close to a side, it's often visually not even possible for a viewer to recognize if it is already inside or still outside.
You still have the bounding box of above, remember? Just pick a point outside the bounding box and use it as starting point for your ray. E.g. the point (Xmin - e/p.y) is outside the polygon for sure.
But what is e? Well, e (actually epsilon) gives the bounding box some padding. As I said, ray tracing fails if we start too close to a polygon line. Since the bounding box might equal the polygon (if the polygon is an axis aligned rectangle, the bounding box is equal to the polygon itself!), we need some padding to make this safe, that's all. How big should you choose e? Not too big. It depends on the coordinate system scale you use for drawing. If your pixel step width is 1.0, then just choose 1.0 (yet 0.1 would have worked as well)
Now that we have the ray with its start and end coordinates, the problem shifts from "is the point within the polygon" to "how often does the ray intersects a polygon side". Therefore we can't just work with the polygon points as before, now we need the actual sides. A side is always defined by two points.
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
You need to test the ray against all sides. Consider the ray to be a vector and every side to be a vector. The ray has to hit each side exactly once or never at all. It can't hit the same side twice. Two lines in 2D space will always intersect exactly once, unless they are parallel, in which case they never intersect. However since vectors have a limited length, two vectors might not be parallel and still never intersect because they are too short to ever meet each other.
// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
// Test if current side intersects with ray.
// If yes, intersections++;
}
if ((intersections & 1) == 1) {
// Inside of polygon
} else {
// Outside of polygon
}
So far so well, but how do you test if two vectors intersect? Here's some C code (not tested), that should do the trick:
#define NO 0
#define YES 1
#define COLLINEAR 2
int areIntersecting(
float v1x1, float v1y1, float v1x2, float v1y2,
float v2x1, float v2y1, float v2x2, float v2y2
) {
float d1, d2;
float a1, a2, b1, b2, c1, c2;
// Convert vector 1 to a line (line 1) of infinite length.
// We want the line in linear equation standard form: A*x + B*y + C = 0
// See: http://en.wikipedia.org/wiki/Linear_equation
a1 = v1y2 - v1y1;
b1 = v1x1 - v1x2;
c1 = (v1x2 * v1y1) - (v1x1 * v1y2);
// Every point (x,y), that solves the equation above, is on the line,
// every point that does not solve it, is not. The equation will have a
// positive result if it is on one side of the line and a negative one
// if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
// 2 into the equation above.
d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
d2 = (a1 * v2x2) + (b1 * v2y2) + c1;
// If d1 and d2 both have the same sign, they are both on the same side
// of our line 1 and in that case no intersection is possible. Careful,
// 0 is a special case, that's why we don't test ">=" and "<=",
// but "<" and ">".
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// The fact that vector 2 intersected the infinite line 1 above doesn't
// mean it also intersects the vector 1. Vector 1 is only a subset of that
// infinite line 1, so it may have intersected that line before the vector
// started or after it ended. To know for sure, we have to repeat the
// the same test the other way round. We start by calculating the
// infinite line 2 in linear equation standard form.
a2 = v2y2 - v2y1;
b2 = v2x1 - v2x2;
c2 = (v2x2 * v2y1) - (v2x1 * v2y2);
// Calculate d1 and d2 again, this time using points of vector 1.
d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
d2 = (a2 * v1x2) + (b2 * v1y2) + c2;
// Again, if both have the same sign (and neither one is 0),
// no intersection is possible.
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// If we get here, only two possibilities are left. Either the two
// vectors intersect in exactly one point or they are collinear, which
// means they intersect in any number of points from zero to infinite.
if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;
// If they are not collinear, they must intersect in exactly one point.
return YES;
}
The input values are the two endpoints of vector 1 (v1x1/v1y1 and v1x2/v1y2) and vector 2 (v2x1/v2y1 and v2x2/v2y2). So you have 2 vectors, 4 points, 8 coordinates. YES and NO are clear. YES increases intersections, NO does nothing.
What about COLLINEAR? It means both vectors lie on the same infinite line, depending on position and length, they don't intersect at all or they intersect in an endless number of points. I'm not absolutely sure how to handle this case, I would not count it as intersection either way. Well, this case is rather rare in practice anyway because of floating point rounding errors; better code would probably not test for == 0.0f but instead for something like < epsilon, where epsilon is a rather small number.
If you need to test a larger number of points, you can certainly speed up the whole thing a bit by keeping the linear equation standard forms of the polygon sides in memory, so you don't have to recalculate these every time. This will save you two floating point multiplications and three floating point subtractions on every test in exchange for storing three floating point values per polygon side in memory. It's a typical memory vs computation time trade off.
Last but not least: If you may use 3D hardware to solve the problem, there is an interesting alternative. Just let the GPU do all the work for you. Create a painting surface that is off screen. Fill it completely with the color black. Now let OpenGL or Direct3D paint your polygon (or even all of your polygons if you just want to test if the point is within any of them, but you don't care for which one) and fill the polygon(s) with a different color, e.g. white. To check if a point is within the polygon, get the color of this point from the drawing surface. This is just a O(1) memory fetch.
Of course this method is only usable if your drawing surface doesn't have to be huge. If it cannot fit into the GPU memory, this method is slower than doing it on the CPU. If it would have to be huge and your GPU supports modern shaders, you can still use the GPU by implementing the ray casting shown above as a GPU shader, which absolutely is possible. For a larger number of polygons or a large number of points to test, this will pay off, consider some GPUs will be able to test 64 to 256 points in parallel. Note however that transferring data from CPU to GPU and back is always expensive, so for just testing a couple of points against a couple of simple polygons, where either the points or the polygons are dynamic and will change frequently, a GPU approach will rarely pay off.
Related Topics
How to Read Data of an Excel File Using C#
Gradient Circles for Map Generator
Oledbcommand Parameters Order and Priority
Adding Multiple Parameterized Variables to a Database in C#
What Are the True Benefits of Expandoobject
Get Mime Type from Filename Extension
How to Know If a Datetime Is Between a Daterange in C#
Asynchronous Locking Based on a Key
How to Compile and Execute New Code at Runtime in .Net
Deserialize Collection of Interface-Instances
Write a Well Designed Async/Non-Async API
When Is a Static Constructor Called in C#
Two Different Dll with Same Namespace
Regular Expression to Remove HTML Tags