I'm trying to create a fast 2D point inside polygon algorithm, for use in hit-testing (e.g. Polygon.contains(p:Point)
). Suggestions for effective techniques would be appreciated.
This question is related to
performance
graphics
collision-detection
polygon
point-in-polygon
I've made a Python implementation of nirg's c++ code:
Inputs
bounding_box_positions: candidate points to filter. (In my implementation created from the bounding box.
(The inputs are lists of tuples in the format: [(xcord, ycord), ...]
)
Returns
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
Again, the idea is taken from here
I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, etc.
We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...
If you want a great algorithm, look to the open source project PostgreSQL source code for the geo work...
I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...
UPDATE
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
In addition, my work used separate functions for "on the line" tests.
...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...
Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!
The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here
If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.
This only works for convex shapes, but Minkowski Portal Refinement, and GJK are also great options for testing if a point is in a polygon. You use minkowski subtraction to subtract the point from the polygon, then run those algorithms to see if the polygon contains the origin.
Also, interestingly, you can describe your shapes a bit more implicitly using support functions which take a direction vector as input and spit out the farthest point along that vector. This allows you to describe any convex shape.. curved, made out of polygons, or mixed. You can also do operations to combine the results of simple support functions to make more complex shapes.
More info: http://xenocollide.snethen.com/mpr2d.html
Also, game programming gems 7 talks about how to do this in 3d (:
Scala version of solution by nirg (assumes bounding rectangle pre-check is done separately):
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.
Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.
Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.
This seems to work in R (apologies for ugliness, would like to see better version!).
pnpoly <- function(nvert,vertx,verty,testx,testy){
c <- FALSE
j <- nvert
for (i in 1:nvert){
if( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))
{c <- !c}
j <- i}
return(c)}
Really like the solution posted by Nirg and edited by bobobobo. I just made it javascript friendly and a little more legible for my use:
function insidePoly(poly, pointx, pointy) {
var i, j;
var inside = false;
for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
}
return inside;
}
I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, etc.
We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...
If you want a great algorithm, look to the open source project PostgreSQL source code for the geo work...
I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...
UPDATE
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
In addition, my work used separate functions for "on the line" tests.
...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...
Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!
For the completeness, here's the lua implementation of the algorithm provided by nirg and discussed by Mecki:
function pnpoly(area, test)
local inside = false
local tx, ty = table.unpack(test)
local j = #area
for i=1, #area do
local vxi, vyi = table.unpack(area[i])
local vxj, vyj = table.unpack(area[j])
if (vyi > ty) ~= (vyj > ty)
and tx < (vxj - vxi)*(ty - vyi)/(vyj - vyi) + vxi
then
inside = not inside
end
j = i
end
return inside
end
The variable area
is a table of points which are in turn stored as 2D tables. Example:
> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true
The link to GitHub Gist.
This question is so interesting. I have another workable idea different from other answers to this post. The idea is to use the sum of angles to decide whether the target is inside or outside. Better known as winding number.
Let x be the target point. Let array [0, 1, .... n] be the all the points of the area. Connect the target point with every border point with a line. If the target point is inside of this area. The sum of all angles will be 360 degrees. If not the angles will be less than 360.
Refer to this image to get a basic understanding of the idea:
My algorithm assumes the clockwise is the positive direction. Here is a potential input:
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
The following is the python code that implements the idea:
def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
a = border[i]
b = border[i + 1]
# calculate distance of vector
A = getDistance(a[0], a[1], b[0], b[1]);
B = getDistance(target[0], target[1], a[0], a[1])
C = getDistance(target[0], target[1], b[0], b[1])
# calculate direction of vector
ta_x = a[0] - target[0]
ta_y = a[1] - target[1]
tb_x = b[0] - target[0]
tb_y = b[1] - target[1]
cross = tb_y * ta_x - tb_x * ta_y
clockwise = cross < 0
# calculate sum of angles
if(clockwise):
degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
else:
degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
if(abs(round(degree) - 360) <= 3):
return True
return False
If you are looking for a java-script library there's a javascript google maps v3 extension for the Polygon class to detect whether or not a point resides within it.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
I realize this is old, but here is a ray casting algorithm implemented in Cocoa, in case anyone is interested. Not sure it is the most efficient way to do things, but it may help someone out.
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
BOOL result;
float aggregateX = 0; //I use these to calculate the centroid of the shape
float aggregateY = 0;
NSPoint firstPoint[1];
[currentPath elementAtIndex:0 associatedPoints:firstPoint];
float olderX = firstPoint[0].x;
float olderY = firstPoint[0].y;
NSPoint interPoint;
int noOfIntersections = 0;
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
aggregateX += points[0].x;
aggregateY += points[0].y;
}
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
//line equations in Ax + By = C form
float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;
float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);
float _A_BAR = olderY - points[0].y;
float _B_BAR = points[0].x - olderX;
float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);
float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
if (det != 0) {
//intersection points with the edges
float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
if (olderX <= points[0].x) {
//doesn't matter in which direction the ray goes, so I send it right-ward.
if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {
noOfIntersections++;
}
} else {
if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
noOfIntersections++;
}
}
}
olderX = points[0].x;
olderY = points[0].y;
}
if (noOfIntersections % 2 == 0) {
result = FALSE;
} else {
result = TRUE;
}
return result;
}
I think the following piece of code is the best solution (taken from here):
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
It's both short and efficient and works both for convex and concave polygons. As suggested before, you should check the bounding rectangle first and treat polygon holes separately.
The idea behind this is pretty simple. The author describes it as follows:
I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This is called the Jordan curve theorem.
The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. So basically it's keeping track of whether the number of edges crossed are even or odd. 0 means even and 1 means odd.
.Net port:
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
int i, j, c = 0;
vertx.Add(1);
vertx.Add(2);
vertx.Add(1);
vertx.Add(4);
vertx.Add(4);
vertx.Add(1);
verty.Add(1);
verty.Add(2);
verty.Add(4);
verty.Add(4);
verty.Add(1);
verty.Add(1);
int nvert = 6; //VĂ©rtices del poligono
double testx = 2;
double testy = 5;
for (i = 0, j = nvert - 1; i < nvert; j = i++)
{
if (((verty[i] > testy) != (verty[j] > testy)) &&
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
c = 1;
}
}
David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.
Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.
Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.
Obj-C version of nirg's answer with sample method for testing points. Nirg's answer worked well for me.
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
NSUInteger nvert = [vertices count];
NSInteger i, j, c = 0;
CGPoint verti, vertj;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
c = !c;
}
return (c ? YES : NO);
}
- (void)testPoint {
NSArray *polygonVertices = [NSArray arrayWithObjects:
[NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
[NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
[NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
[NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
nil
];
CGPoint tappedPoint = CGPointMake(23.0, 70.0);
if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
NSLog(@"YES");
} else {
NSLog(@"NO");
}
}
You can do this by checking if the area formed by connecting the desired point to the vertices of your polygon matches the area of the polygon itself.
Or you could check if the sum of the inner angles from your point to each pair of two consecutive polygon vertices to your check point sums to 360, but I have the feeling that the first option is quicker because it doesn't involve divisions nor calculations of inverse of trigonometric functions.
I don't know what happens if your polygon has a hole inside it but it seems to me that the main idea can be adapted to this situation
You can as well post the question in a math community. I bet they have one million ways of doing that
Here is golang version of @nirg answer (inspired by C# code by @@m-katz)
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
When using qt (Qt 4.3+), one can use QPolygon's function containsPoint
If you're using Google Map SDK and want to check if a point is inside a polygon, you can try to use GMSGeometryContainsLocation
. It works great!! Here is how that works,
if GMSGeometryContainsLocation(point, polygon, true) {
print("Inside this polygon.")
} else {
print("outside this polygon")
}
Here is the reference: https://developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils#gaba958d3776d49213404af249419d0ffd
Surprised nobody brought this up earlier, but for the pragmatists requiring a database: MongoDB has excellent support for Geo queries including this one.
What you are looking for is:
db.neighborhoods.findOne({ geometry: { $geoIntersects: { $geometry: { type: "Point", coordinates: [ "longitude", "latitude" ] } } } })
Neighborhoods
is the collection that stores one or more polygons in standard GeoJson format. If the query returns null it is not intersected otherwise it is.
Very well documented here: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
The performance for more than 6,000 points classified in a 330 irregular polygon grid was less than one minute with no optimization at all and including the time to update documents with their respective polygon.
Here is a JavaScript variant of the answer by M. Katz based on Nirg's approach:
function pointIsInPoly(p, polygon) {
var isInside = false;
var minX = polygon[0].x, maxX = polygon[0].x;
var minY = polygon[0].y, maxY = polygon[0].y;
for (var n = 1; n < polygon.length; n++) {
var q = polygon[n];
minX = Math.min(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
VBA VERSION:
Note: Remember that if your polygon is an area within a map that Latitude/Longitude are Y/X values as opposed to X/Y (Latitude = Y, Longitude = X) due to from what I understand are historical implications from way back when Longitude was not a measurement.
CLASS MODULE: CPoint
Private pXValue As Double
Private pYValue As Double
'''''X Value Property'''''
Public Property Get X() As Double
X = pXValue
End Property
Public Property Let X(Value As Double)
pXValue = Value
End Property
'''''Y Value Property'''''
Public Property Get Y() As Double
Y = pYValue
End Property
Public Property Let Y(Value As Double)
pYValue = Value
End Property
MODULE:
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean
Dim i As Integer
Dim j As Integer
Dim q As Object
Dim minX As Double
Dim maxX As Double
Dim minY As Double
Dim maxY As Double
minX = polygon(0).X
maxX = polygon(0).X
minY = polygon(0).Y
maxY = polygon(0).Y
For i = 1 To UBound(polygon)
Set q = polygon(i)
minX = vbMin(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
isPointInPolygon = False
Exit Function
End If
' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
isPointInPolygon = False
i = 0
j = UBound(polygon)
Do While i < UBound(polygon) + 1
If (polygon(i).Y > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
ElseIf (polygon(i).Y < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
End If
j = i
i = i + 1
Loop
End Function
Function vbMax(n1, n2) As Double
vbMax = IIf(n1 > n2, n1, n2)
End Function
Function vbMin(n1, n2) As Double
vbMin = IIf(n1 > n2, n2, n1)
End Function
Sub TestPointInPolygon()
Dim i As Integer
Dim InPolygon As Boolean
' MARKER Object
Dim p As CPoint
Set p = New CPoint
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
The answer depends on if you have the simple or complex polygons. Simple polygons must not have any line segment intersections. So they can have the holes but lines can't cross each other. Complex regions can have the line intersections - so they can have the overlapping regions, or regions that touch each other just by a single point.
For simple polygons the best algorithm is Ray casting (Crossing number) algorithm. For complex polygons, this algorithm doesn't detect points that are inside the overlapping regions. So for complex polygons you have to use Winding number algorithm.
Here is an excellent article with C implementation of both algorithms. I tried them and they work well.
Java Version:
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.
I like this method better because it is more robust and less dependent on numerical precision.
Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.
EDIT: By The Way, this method works with concave and convex polygons.
EDIT: I recently found a whole Wikipedia article on the topic.
This is a presumably slightly less optimized version of the C code from here which was sourced from this page.
My C++ version uses a std::vector<std::pair<double, double>>
and two doubles as an x and y. The logic should be exactly the same as the original C code, but I find mine easier to read. I can't speak for the performance.
bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
bool in_poly = false;
auto num_verts = verts.size();
for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
double x1 = verts[i].first;
double y1 = verts[i].second;
double x2 = verts[j].first;
double y2 = verts[j].second;
if (((y1 > point_y) != (y2 > point_y)) &&
(point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
in_poly = !in_poly;
}
return in_poly;
}
The original C code is
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
Most of the answers in this question are not handling all corner cases well. Some subtle corner cases like below: This is a javascript version with all corner cases well handled.
/** Get relationship between a point and a polygon using ray-casting algorithm
* @param {{x:number, y:number}} P: point to check
* @param {{x:number, y:number}[]} polygon: the polygon
* @returns -1: outside, 0: on edge, 1: inside
*/
function relationPP(P, polygon) {
const between = (p, a, b) => p >= a && p <= b || p <= a && p >= b
let inside = false
for (let i = polygon.length-1, j = 0; j < polygon.length; i = j, j++) {
const A = polygon[i]
const B = polygon[j]
// corner cases
if (P.x == A.x && P.y == A.y || P.x == B.x && P.y == B.y) return 0
if (A.y == B.y && P.y == A.y && between(P.x, A.x, B.x)) return 0
if (between(P.y, A.y, B.y)) { // if P inside the vertical range
// filter out "ray pass vertex" problem by treating the line a little lower
if (P.y == A.y && B.y >= A.y || P.y == B.y && A.y >= B.y) continue
// calc cross product `PA X PB`, P lays on left side of AB if c > 0
const c = (A.x - P.x) * (B.y - P.y) - (B.x - P.x) * (A.y - P.y)
if (c == 0) return 0
if ((A.y < B.y) == (c > 0)) inside = !inside
}
}
return inside? 1 : -1
}
.Net port:
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
int i, j, c = 0;
vertx.Add(1);
vertx.Add(2);
vertx.Add(1);
vertx.Add(4);
vertx.Add(4);
vertx.Add(1);
verty.Add(1);
verty.Add(2);
verty.Add(4);
verty.Add(4);
verty.Add(1);
verty.Add(1);
int nvert = 6; //VĂ©rtices del poligono
double testx = 2;
double testy = 5;
for (i = 0, j = nvert - 1; i < nvert; j = i++)
{
if (((verty[i] > testy) != (verty[j] > testy)) &&
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
c = 1;
}
}
For Detecting hit on Polygon we need to test two things:
If you are looking for a java-script library there's a javascript google maps v3 extension for the Polygon class to detect whether or not a point resides within it.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.
Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.
Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.
Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.
I like this method better because it is more robust and less dependent on numerical precision.
Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.
EDIT: By The Way, this method works with concave and convex polygons.
EDIT: I recently found a whole Wikipedia article on the topic.
Scala version of solution by nirg (assumes bounding rectangle pre-check is done separately):
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
Obj-C version of nirg's answer with sample method for testing points. Nirg's answer worked well for me.
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
NSUInteger nvert = [vertices count];
NSInteger i, j, c = 0;
CGPoint verti, vertj;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
c = !c;
}
return (c ? YES : NO);
}
- (void)testPoint {
NSArray *polygonVertices = [NSArray arrayWithObjects:
[NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
[NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
[NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
[NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
nil
];
CGPoint tappedPoint = CGPointMake(23.0, 70.0);
if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
NSLog(@"YES");
} else {
NSLog(@"NO");
}
}
There is nothing more beutiful than an inductive definition of a problem. For the sake of completeness here you have a version in prolog which might also clarify the thoughs behind ray casting:
Based on the simulation of simplicity algorithm in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Some helper predicates:
exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).
inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).
get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
The equation of a line given 2 points A and B (Line(A,B)) is:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), i.e the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.
is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA));
is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).
in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines).
traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).
% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
Surprised nobody brought this up earlier, but for the pragmatists requiring a database: MongoDB has excellent support for Geo queries including this one.
What you are looking for is:
db.neighborhoods.findOne({ geometry: { $geoIntersects: { $geometry: { type: "Point", coordinates: [ "longitude", "latitude" ] } } } })
Neighborhoods
is the collection that stores one or more polygons in standard GeoJson format. If the query returns null it is not intersected otherwise it is.
Very well documented here: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
The performance for more than 6,000 points classified in a 330 irregular polygon grid was less than one minute with no optimization at all and including the time to update documents with their respective polygon.
Here is a JavaScript variant of the answer by M. Katz based on Nirg's approach:
function pointIsInPoly(p, polygon) {
var isInside = false;
var minX = polygon[0].x, maxX = polygon[0].x;
var minY = polygon[0].y, maxY = polygon[0].y;
for (var n = 1; n < polygon.length; n++) {
var q = polygon[n];
minX = Math.min(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
I think the following piece of code is the best solution (taken from here):
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
It's both short and efficient and works both for convex and concave polygons. As suggested before, you should check the bounding rectangle first and treat polygon holes separately.
The idea behind this is pretty simple. The author describes it as follows:
I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This is called the Jordan curve theorem.
The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. So basically it's keeping track of whether the number of edges crossed are even or odd. 0 means even and 1 means odd.
This is a presumably slightly less optimized version of the C code from here which was sourced from this page.
My C++ version uses a std::vector<std::pair<double, double>>
and two doubles as an x and y. The logic should be exactly the same as the original C code, but I find mine easier to read. I can't speak for the performance.
bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
bool in_poly = false;
auto num_verts = verts.size();
for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
double x1 = verts[i].first;
double y1 = verts[i].second;
double x2 = verts[j].first;
double y2 = verts[j].second;
if (((y1 > point_y) != (y2 > point_y)) &&
(point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
in_poly = !in_poly;
}
return in_poly;
}
The original C code is
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
To deal with the following special cases in Ray casting algorithm:
Check Determining Whether A Point Is Inside A Complex Polygon. The article provides an easy way to resolve them so there will be no special treatment required for the above cases.
Here is a C# version of the answer given by nirg, which comes from this RPI professor. Note that use of the code from that RPI source requires attribution.
A bounding box check has been added at the top. However, as James Brown points out, the main code is almost as fast as the bounding box check itself, so the bounding box check can actually slow the overall operation, in the case that most of the points you are checking are inside the bounding box. So you could leave the bounding box check out, or an alternative would be to precompute the bounding boxes of your polygons if they don't change shape too often.
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
minX = Math.Min( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
{
return false;
}
// https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
{
if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here
If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.
You can do this by checking if the area formed by connecting the desired point to the vertices of your polygon matches the area of the polygon itself.
Or you could check if the sum of the inner angles from your point to each pair of two consecutive polygon vertices to your check point sums to 360, but I have the feeling that the first option is quicker because it doesn't involve divisions nor calculations of inverse of trigonometric functions.
I don't know what happens if your polygon has a hole inside it but it seems to me that the main idea can be adapted to this situation
You can as well post the question in a math community. I bet they have one million ways of doing that
Most of the answers in this question are not handling all corner cases well. Some subtle corner cases like below: This is a javascript version with all corner cases well handled.
/** Get relationship between a point and a polygon using ray-casting algorithm
* @param {{x:number, y:number}} P: point to check
* @param {{x:number, y:number}[]} polygon: the polygon
* @returns -1: outside, 0: on edge, 1: inside
*/
function relationPP(P, polygon) {
const between = (p, a, b) => p >= a && p <= b || p <= a && p >= b
let inside = false
for (let i = polygon.length-1, j = 0; j < polygon.length; i = j, j++) {
const A = polygon[i]
const B = polygon[j]
// corner cases
if (P.x == A.x && P.y == A.y || P.x == B.x && P.y == B.y) return 0
if (A.y == B.y && P.y == A.y && between(P.x, A.x, B.x)) return 0
if (between(P.y, A.y, B.y)) { // if P inside the vertical range
// filter out "ray pass vertex" problem by treating the line a little lower
if (P.y == A.y && B.y >= A.y || P.y == B.y && A.y >= B.y) continue
// calc cross product `PA X PB`, P lays on left side of AB if c > 0
const c = (A.x - P.x) * (B.y - P.y) - (B.x - P.x) * (A.y - P.y)
if (c == 0) return 0
if ((A.y < B.y) == (c > 0)) inside = !inside
}
}
return inside? 1 : -1
}
I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, etc.
We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...
If you want a great algorithm, look to the open source project PostgreSQL source code for the geo work...
I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...
UPDATE
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
In addition, my work used separate functions for "on the line" tests.
...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...
Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!
Java Version:
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
C# version of nirg's answer is here: I'll just share the code. It may save someone some time.
public static bool IsPointInPolygon(IList<Point> polygon, Point 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;
}
Here is golang version of @nirg answer (inspired by C# code by @@m-katz)
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
This question is so interesting. I have another workable idea different from other answers to this post. The idea is to use the sum of angles to decide whether the target is inside or outside. Better known as winding number.
Let x be the target point. Let array [0, 1, .... n] be the all the points of the area. Connect the target point with every border point with a line. If the target point is inside of this area. The sum of all angles will be 360 degrees. If not the angles will be less than 360.
Refer to this image to get a basic understanding of the idea:
My algorithm assumes the clockwise is the positive direction. Here is a potential input:
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
The following is the python code that implements the idea:
def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
a = border[i]
b = border[i + 1]
# calculate distance of vector
A = getDistance(a[0], a[1], b[0], b[1]);
B = getDistance(target[0], target[1], a[0], a[1])
C = getDistance(target[0], target[1], b[0], b[1])
# calculate direction of vector
ta_x = a[0] - target[0]
ta_y = a[1] - target[1]
tb_x = b[0] - target[0]
tb_y = b[1] - target[1]
cross = tb_y * ta_x - tb_x * ta_y
clockwise = cross < 0
# calculate sum of angles
if(clockwise):
degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
else:
degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
if(abs(round(degree) - 360) <= 3):
return True
return False
VBA VERSION:
Note: Remember that if your polygon is an area within a map that Latitude/Longitude are Y/X values as opposed to X/Y (Latitude = Y, Longitude = X) due to from what I understand are historical implications from way back when Longitude was not a measurement.
CLASS MODULE: CPoint
Private pXValue As Double
Private pYValue As Double
'''''X Value Property'''''
Public Property Get X() As Double
X = pXValue
End Property
Public Property Let X(Value As Double)
pXValue = Value
End Property
'''''Y Value Property'''''
Public Property Get Y() As Double
Y = pYValue
End Property
Public Property Let Y(Value As Double)
pYValue = Value
End Property
MODULE:
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean
Dim i As Integer
Dim j As Integer
Dim q As Object
Dim minX As Double
Dim maxX As Double
Dim minY As Double
Dim maxY As Double
minX = polygon(0).X
maxX = polygon(0).X
minY = polygon(0).Y
maxY = polygon(0).Y
For i = 1 To UBound(polygon)
Set q = polygon(i)
minX = vbMin(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
isPointInPolygon = False
Exit Function
End If
' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
isPointInPolygon = False
i = 0
j = UBound(polygon)
Do While i < UBound(polygon) + 1
If (polygon(i).Y > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
ElseIf (polygon(i).Y < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
End If
j = i
i = i + 1
Loop
End Function
Function vbMax(n1, n2) As Double
vbMax = IIf(n1 > n2, n1, n2)
End Function
Function vbMin(n1, n2) As Double
vbMin = IIf(n1 > n2, n2, n1)
End Function
Sub TestPointInPolygon()
Dim i As Integer
Dim InPolygon As Boolean
' MARKER Object
Dim p As CPoint
Set p = New CPoint
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
This seems to work in R (apologies for ugliness, would like to see better version!).
pnpoly <- function(nvert,vertx,verty,testx,testy){
c <- FALSE
j <- nvert
for (i in 1:nvert){
if( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))
{c <- !c}
j <- i}
return(c)}
Really like the solution posted by Nirg and edited by bobobobo. I just made it javascript friendly and a little more legible for my use:
function insidePoly(poly, pointx, pointy) {
var i, j;
var inside = false;
for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
}
return inside;
}
Heres a point in polygon test in C that isn't using ray-casting. And it can work for overlapping areas (self intersections), see the use_holes
argument.
/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);
/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
const bool use_holes)
{
/* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
float angletot = 0.0;
float fp1[2], fp2[2];
unsigned int i;
const float *p1, *p2;
p1 = verts[nr - 1];
/* first vector */
fp1[0] = p1[0] - pt[0];
fp1[1] = p1[1] - pt[1];
for (i = 0; i < nr; i++) {
p2 = verts[i];
/* second vector */
fp2[0] = p2[0] - pt[0];
fp2[1] = p2[1] - pt[1];
/* dot and angle and cross */
angletot += angle_signed_v2v2(fp1, fp2);
/* circulate */
copy_v2_v2(fp1, fp2);
p1 = p2;
}
angletot = fabsf(angletot);
if (use_holes) {
const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
angletot -= nested * (float)(M_PI * 2.0);
return (angletot > 4.0f) != ((int)nested % 2);
}
else {
return (angletot > 4.0f);
}
}
/* math lib */
static float dot_v2v2(const float a[2], const float b[2])
{
return a[0] * b[0] + a[1] * b[1];
}
static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
return atan2f(perp_dot, dot_v2v2(v1, v2));
}
static void copy_v2_v2(float r[2], const float a[2])
{
r[0] = a[0];
r[1] = a[1];
}
Note: this is one of the less optimal methods since it includes a lot of calls to atan2f
, but it may be of interest to developers reading this thread (in my tests its ~23x slower then using the line intersection method).
I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, etc.
We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...
If you want a great algorithm, look to the open source project PostgreSQL source code for the geo work...
I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...
UPDATE
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
In addition, my work used separate functions for "on the line" tests.
...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...
Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!
If you're using Google Map SDK and want to check if a point is inside a polygon, you can try to use GMSGeometryContainsLocation
. It works great!! Here is how that works,
if GMSGeometryContainsLocation(point, polygon, true) {
print("Inside this polygon.")
} else {
print("outside this polygon")
}
Here is the reference: https://developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils#gaba958d3776d49213404af249419d0ffd
The Eric Haines article cited by bobobobo is really excellent. Particularly interesting are the tables comparing performance of the algorithms; the angle summation method is really bad compared to the others. Also interesting is that optimisations like using a lookup grid to further subdivide the polygon into "in" and "out" sectors can make the test incredibly fast even on polygons with > 1000 sides.
Anyway, it's early days but my vote goes to the "crossings" method, which is pretty much what Mecki describes I think. However I found it most succintly described and codified by David Bourke. I love that there is no real trigonometry required, and it works for convex and concave, and it performs reasonably well as the number of sides increases.
By the way, here's one of the performance tables from the Eric Haines' article for interest, testing on random polygons.
number of edges per polygon
3 4 10 100 1000
MacMartin 2.9 3.2 5.9 50.6 485
Crossings 3.1 3.4 6.8 60.0 624
Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787
Triangle Fan 1.2 2.1 7.3 85.4 865
Barycentric 2.1 3.8 13.8 160.7 1665
Angle Summation 56.2 70.4 153.6 1403.8 14693
Grid (100x100) 1.5 1.5 1.6 2.1 9.8
Grid (20x20) 1.7 1.7 1.9 5.7 42.2
Bins (100) 1.8 1.9 2.7 15.1 117
Bins (20) 2.1 2.2 3.7 26.3 278
For the completeness, here's the lua implementation of the algorithm provided by nirg and discussed by Mecki:
function pnpoly(area, test)
local inside = false
local tx, ty = table.unpack(test)
local j = #area
for i=1, #area do
local vxi, vyi = table.unpack(area[i])
local vxj, vyj = table.unpack(area[j])
if (vyi > ty) ~= (vyj > ty)
and tx < (vxj - vxi)*(ty - vyi)/(vyj - vyi) + vxi
then
inside = not inside
end
j = i
end
return inside
end
The variable area
is a table of points which are in turn stored as 2D tables. Example:
> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true
The link to GitHub Gist.
Heres a point in polygon test in C that isn't using ray-casting. And it can work for overlapping areas (self intersections), see the use_holes
argument.
/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);
/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
const bool use_holes)
{
/* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
float angletot = 0.0;
float fp1[2], fp2[2];
unsigned int i;
const float *p1, *p2;
p1 = verts[nr - 1];
/* first vector */
fp1[0] = p1[0] - pt[0];
fp1[1] = p1[1] - pt[1];
for (i = 0; i < nr; i++) {
p2 = verts[i];
/* second vector */
fp2[0] = p2[0] - pt[0];
fp2[1] = p2[1] - pt[1];
/* dot and angle and cross */
angletot += angle_signed_v2v2(fp1, fp2);
/* circulate */
copy_v2_v2(fp1, fp2);
p1 = p2;
}
angletot = fabsf(angletot);
if (use_holes) {
const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
angletot -= nested * (float)(M_PI * 2.0);
return (angletot > 4.0f) != ((int)nested % 2);
}
else {
return (angletot > 4.0f);
}
}
/* math lib */
static float dot_v2v2(const float a[2], const float b[2])
{
return a[0] * b[0] + a[1] * b[1];
}
static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
return atan2f(perp_dot, dot_v2v2(v1, v2));
}
static void copy_v2_v2(float r[2], const float a[2])
{
r[0] = a[0];
r[1] = a[1];
}
Note: this is one of the less optimal methods since it includes a lot of calls to atan2f
, but it may be of interest to developers reading this thread (in my tests its ~23x slower then using the line intersection method).
The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here
If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.
I've made a Python implementation of nirg's c++ code:
Inputs
bounding_box_positions: candidate points to filter. (In my implementation created from the bounding box.
(The inputs are lists of tuples in the format: [(xcord, ycord), ...]
)
Returns
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
Again, the idea is taken from here
David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.
Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.
Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.
The answer depends on if you have the simple or complex polygons. Simple polygons must not have any line segment intersections. So they can have the holes but lines can't cross each other. Complex regions can have the line intersections - so they can have the overlapping regions, or regions that touch each other just by a single point.
For simple polygons the best algorithm is Ray casting (Crossing number) algorithm. For complex polygons, this algorithm doesn't detect points that are inside the overlapping regions. So for complex polygons you have to use Winding number algorithm.
Here is an excellent article with C implementation of both algorithms. I tried them and they work well.
Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.
I like this method better because it is more robust and less dependent on numerical precision.
Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.
EDIT: By The Way, this method works with concave and convex polygons.
EDIT: I recently found a whole Wikipedia article on the topic.
Swift version of the answer by nirg:
extension CGPoint {
func isInsidePolygon(vertices: [CGPoint]) -> Bool {
guard !vertices.isEmpty else { return false }
var j = vertices.last!, c = false
for i in vertices {
let a = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
Here is a C# version of the answer given by nirg, which comes from this RPI professor. Note that use of the code from that RPI source requires attribution.
A bounding box check has been added at the top. However, as James Brown points out, the main code is almost as fast as the bounding box check itself, so the bounding box check can actually slow the overall operation, in the case that most of the points you are checking are inside the bounding box. So you could leave the bounding box check out, or an alternative would be to precompute the bounding boxes of your polygons if they don't change shape too often.
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
minX = Math.Min( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
{
return false;
}
// https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
{
if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.
I like this method better because it is more robust and less dependent on numerical precision.
Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.
EDIT: By The Way, this method works with concave and convex polygons.
EDIT: I recently found a whole Wikipedia article on the topic.
The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here
If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.
C# version of nirg's answer is here: I'll just share the code. It may save someone some time.
public static bool IsPointInPolygon(IList<Point> polygon, Point 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;
}
I realize this is old, but here is a ray casting algorithm implemented in Cocoa, in case anyone is interested. Not sure it is the most efficient way to do things, but it may help someone out.
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
BOOL result;
float aggregateX = 0; //I use these to calculate the centroid of the shape
float aggregateY = 0;
NSPoint firstPoint[1];
[currentPath elementAtIndex:0 associatedPoints:firstPoint];
float olderX = firstPoint[0].x;
float olderY = firstPoint[0].y;
NSPoint interPoint;
int noOfIntersections = 0;
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
aggregateX += points[0].x;
aggregateY += points[0].y;
}
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
//line equations in Ax + By = C form
float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;
float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);
float _A_BAR = olderY - points[0].y;
float _B_BAR = points[0].x - olderX;
float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);
float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
if (det != 0) {
//intersection points with the edges
float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
if (olderX <= points[0].x) {
//doesn't matter in which direction the ray goes, so I send it right-ward.
if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {
noOfIntersections++;
}
} else {
if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
noOfIntersections++;
}
}
}
olderX = points[0].x;
olderY = points[0].y;
}
if (noOfIntersections % 2 == 0) {
result = FALSE;
} else {
result = TRUE;
}
return result;
}
There is nothing more beutiful than an inductive definition of a problem. For the sake of completeness here you have a version in prolog which might also clarify the thoughs behind ray casting:
Based on the simulation of simplicity algorithm in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Some helper predicates:
exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).
inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).
get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
The equation of a line given 2 points A and B (Line(A,B)) is:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), i.e the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.
is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA));
is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).
in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines).
traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).
% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
The Eric Haines article cited by bobobobo is really excellent. Particularly interesting are the tables comparing performance of the algorithms; the angle summation method is really bad compared to the others. Also interesting is that optimisations like using a lookup grid to further subdivide the polygon into "in" and "out" sectors can make the test incredibly fast even on polygons with > 1000 sides.
Anyway, it's early days but my vote goes to the "crossings" method, which is pretty much what Mecki describes I think. However I found it most succintly described and codified by David Bourke. I love that there is no real trigonometry required, and it works for convex and concave, and it performs reasonably well as the number of sides increases.
By the way, here's one of the performance tables from the Eric Haines' article for interest, testing on random polygons.
number of edges per polygon
3 4 10 100 1000
MacMartin 2.9 3.2 5.9 50.6 485
Crossings 3.1 3.4 6.8 60.0 624
Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787
Triangle Fan 1.2 2.1 7.3 85.4 865
Barycentric 2.1 3.8 13.8 160.7 1665
Angle Summation 56.2 70.4 153.6 1403.8 14693
Grid (100x100) 1.5 1.5 1.6 2.1 9.8
Grid (20x20) 1.7 1.7 1.9 5.7 42.2
Bins (100) 1.8 1.9 2.7 15.1 117
Bins (20) 2.1 2.2 3.7 26.3 278
To deal with the following special cases in Ray casting algorithm:
Check Determining Whether A Point Is Inside A Complex Polygon. The article provides an easy way to resolve them so there will be no special treatment required for the above cases.
Swift version of the answer by nirg:
extension CGPoint {
func isInsidePolygon(vertices: [CGPoint]) -> Bool {
guard !vertices.isEmpty else { return false }
var j = vertices.last!, c = false
for i in vertices {
let a = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
When using qt (Qt 4.3+), one can use QPolygon's function containsPoint
This only works for convex shapes, but Minkowski Portal Refinement, and GJK are also great options for testing if a point is in a polygon. You use minkowski subtraction to subtract the point from the polygon, then run those algorithms to see if the polygon contains the origin.
Also, interestingly, you can describe your shapes a bit more implicitly using support functions which take a direction vector as input and spit out the farthest point along that vector. This allows you to describe any convex shape.. curved, made out of polygons, or mixed. You can also do operations to combine the results of simple support functions to make more complex shapes.
More info: http://xenocollide.snethen.com/mpr2d.html
Also, game programming gems 7 talks about how to do this in 3d (:
For Detecting hit on Polygon we need to test two things:
Source: Stackoverflow.com