58

I would like to determine a polygon and implement an algorithm which would check if a point is inside or outside the polygon.

Does anyone know if there is any example available of any similar algorithm?

tskuzzy
  • 35,812
  • 14
  • 73
  • 140
Niko Gamulin
  • 66,025
  • 95
  • 221
  • 286

16 Answers16

68

If i remember correctly, the algorithm is to draw a horizontal line through your test point. Count how many lines of of the polygon you intersect to reach your point.

If the answer is odd, you're inside. If the answer is even, you're outside.

Edit: Yeah, what he said (Wikipedia):

alt text

Community
  • 1
  • 1
Ian Boyd
  • 246,734
  • 253
  • 869
  • 1,219
  • 7
    I appreciate your inclusion of the picture for quick reference. It really says it all. – BigBeagle Sep 16 '09 at 20:05
  • 2
    Don't forget the logic to process if the ray passes through a vertice. And to determine if a point is on the edge, is it in or out of the poly. The links above provide more detail on this. – htm11h Sep 18 '13 at 12:54
  • "Unfortunately, this method won't work if the point is on the edge of the polygon." - Does this mean that sometimes this algorithm won't work? – Ciaran Gallagher Aug 31 '17 at 18:38
  • @CiaranGallagher If a point is neither inside nor outside the polygon, then the algorithm that answers if the point is inside or outside the polygon will not be able to answer you. – Ian Boyd Aug 31 '17 at 23:15
42

C# code

bool IsPointInPolygon(List<Loc> poly, Loc point)
{
    int i, j;
    bool c = false;
    for (i = 0, j = poly.Count - 1; i < poly.Count; j = i++)
    {
        if ((((poly[i].Lt <= point.Lt) && (point.Lt < poly[j].Lt)) 
                || ((poly[j].Lt <= point.Lt) && (point.Lt < poly[i].Lt))) 
                && (point.Lg < (poly[j].Lg - poly[i].Lg) * (point.Lt - poly[i].Lt) 
                    / (poly[j].Lt - poly[i].Lt) + poly[i].Lg))
        {

            c = !c;
        }
    }

    return c;
}

Location class

public class Loc
{
    private double lt;
    private double lg;

    public double Lg
    {
        get { return lg; }
        set { lg = value; }
    }

    public double Lt
    {
        get { return lt; }
        set { lt = value; }
    }

    public Loc(double lt, double lg)
    {
        this.lt = lt;
        this.lg = lg;
    }
}
live2
  • 3,771
  • 2
  • 37
  • 46
kober
  • 832
  • 7
  • 13
  • 1
    Just mention that I've implemented this and tested it and it seems to work. Just wish the author had provided a few comments and an indication of which algorithm this implements. – RenniePet Jun 11 '13 at 06:20
  • 7
    For anyone still wondering, this is obviously the even-odd rule algorithm. – Alexios Nov 04 '13 at 18:32
  • I tested the above code with WGS84 coordinates in C# and it seemed to work for me too! thanks – programmer Apr 26 '14 at 05:22
  • @Alexios, Is this algorithm acceptable for non-convex polygons? – CAMOBAP Jun 02 '14 at 13:15
  • 1
    I am converting this into vbscript. I get a divide by zero error here - " / (poly[j].Lt - poly[i].Lt)" when poly[j].Lt and poly[i].Lt are the same. – MrVimes Jul 27 '17 at 13:53
  • ... I've crudely solved my divide by zero problem by detecting equality and adding 0.0000001 to one of the values, but I am wondering if I'm doing something wrong to be in a divide by zero situation. It just means two points have the same latitude, which is entirely feasible. – MrVimes Jul 27 '17 at 14:30
  • Working perfectly. Thanks a bunch! – Sohail Ahmad Jun 21 '18 at 09:09
  • @MrVimes The first two && conditions check that poly[i].Lt and poly[j].Lt aren't equal indirectly by comparing them to point.Lt. Not sure about VBScript, but it means that in C# it never gets to the divide by zero scenario. – Bloopy Apr 30 '19 at 05:13
  • 1
    this algorithm failes if you work with a polygon where any line is crossing the date line on 180 degree longitude. So we add to every longitude that is < 0 +360 degrees then it works well! I added our version also as an answer. – Nelly Oct 14 '19 at 13:05
16

After searching the web and trying various implementations and porting them from C++ to C# I finally got my code straight:

        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;

    }

    private static int isLeft(LatLong P0, LatLong P1, LatLong P2)
    {
        double calc = ((P1.Lon - P0.Lon) * (P2.Lat - P0.Lat)
                - (P2.Lon - P0.Lon) * (P1.Lat - P0.Lat));
        if (calc > 0)
            return 1;
        else if (calc < 0)
            return -1;
        else
            return 0;
    }

The isLeft function was giving me rounding problems and I spent hours without realizing that I was doing the conversion wrong, so forgive me for the lame if block at the end of that function.

BTW, this is the original code and article: http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm

Manuel Castro
  • 2,460
  • 1
  • 19
  • 14
  • 6
    Not that I know of, but any seasoned PHP programmer should be able to port the original code to PHP – Manuel Castro Apr 23 '11 at 20:27
  • This is on the money. I wrote xUnit tests around the above and it checks out perfectly. I have also provided a refactor of the above code as OOP C# classes. I hope you find this helpful: https://stackoverflow.com/questions/46144205/point-in-polygon-using-winding-number – Jim Speaker Sep 10 '17 at 19:02
5

I think there is a simpler and more efficient solution.

Here is the code in C++. I should be simple to convert it to C#.

int pnpoly(int npol, float *xp, float *yp, float x, float y)
{
  int i, j, c = 0;
  for (i = 0, j = npol-1; i < npol; j = i++) {
    if ((((yp[i] <= y) && (y < yp[j])) ||
         ((yp[j] <= y) && (y < yp[i]))) &&
        (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
      c = !c;
  }
  return c;
}
Shawn Chin
  • 84,080
  • 19
  • 162
  • 191
wael
  • 51
  • 1
  • 1
5

By far the best explanation and implementation can be found at Point In Polygon Winding Number Inclusion

There is even a C++ implementation at the end of the well explained article. This site also contains some great algorithms/solutions for other geometry based problems.

I have modified and used the C++ implementation and also created a C# implementation. You definitely want to use the Winding Number algorithm as it is more accurate than the edge crossing algorithm and it is very fast.

eesh
  • 1,384
  • 2
  • 13
  • 25
2

The complete solution in asp.Net C#, you can see the complete detail here, you can see how to find point(lat,lon) whether its inside or Outside the Polygon using the latitude and longitudes ? Article Reference Link

private static bool checkPointExistsInGeofencePolygon(string latlnglist, string lat, string lng) {

    List<Loc> objList = new List<Loc>();
    // sample string should be like this strlatlng = "39.11495,-76.873259|39.114588,-76.872808|39.112921,-76.870373|";
    string[] arr = latlnglist.Split('|');
    for (int i = 0; i <= arr.Length - 1; i++)
    {
        string latlng = arr[i];
        string[] arrlatlng = latlng.Split(',');

        Loc er = new Loc(Convert.ToDouble(arrlatlng[0]), Convert.ToDouble(arrlatlng[1]));
        objList.Add(er);
    }
    Loc pt = new Loc(Convert.ToDouble(lat), Convert.ToDouble(lng));

    if (IsPointInPolygon(objList, pt) == true)
    {
          return true;
    }
    else
    {
           return false;
    }
}
private static bool IsPointInPolygon(List<Loc> poly, Loc point)
{
    int i, j;
    bool c = false;
    for (i = 0, j = poly.Count - 1; i < poly.Count; j = i++)
    {
        if ((((poly[i].Lt <= point.Lt) && (point.Lt < poly[j].Lt)) |
            ((poly[j].Lt <= point.Lt) && (point.Lt < poly[i].Lt))) &&
            (point.Lg < (poly[j].Lg - poly[i].Lg) * (point.Lt - poly[i].Lt) / (poly[j].Lt - poly[i].Lt) + poly[i].Lg))
            c = !c;
    }
    return c;
}
2

Just a heads up (using answer as I can't comment), if you want to use point-in-polygon for geo fencing, then you need to change your algorithm to work with spherical coordinates. -180 longitude is the same as 180 longitude and point-in-polygon will break in such situation.

Justin Zhang
  • 4,433
  • 1
  • 12
  • 9
2

Relating to kobers answer I worked it out with more readable clean code and changed the longitudes that crosses the date border:

public bool IsPointInPolygon(List<PointPosition> polygon, double latitude, double longitude)
{
  bool isInIntersection = false;
  int actualPointIndex = 0;
  int pointIndexBeforeActual = polygon.Count - 1;

  var offset = calculateLonOffsetFromDateLine(polygon);
  longitude = longitude < 0.0 ? longitude + offset : longitude;

  foreach (var actualPointPosition in polygon)
  {
    var p1Lat = actualPointPosition.Latitude;
    var p1Lon = actualPointPosition.Longitude;

    var p0Lat = polygon[pointIndexBeforeActual].Latitude;
    var p0Lon = polygon[pointIndexBeforeActual].Longitude;

    if (p1Lon < 0.0) p1Lon += offset;
    if (p0Lon < 0.0) p0Lon += offset;

    // Jordan curve theorem - odd even rule algorithm
    if (isPointLatitudeBetweenPolyLine(p0Lat, p1Lat, latitude)
    && isPointRightFromPolyLine(p0Lat, p0Lon, p1Lat, p1Lon, latitude, longitude))
    {
      isInIntersection = !isInIntersection;
    }

    pointIndexBeforeActual = actualPointIndex;
    actualPointIndex++;
  }

  return isInIntersection;
}

private double calculateLonOffsetFromDateLine(List<PointPosition> polygon)
{
  double offset = 0.0;
  var maxLonPoly = polygon.Max(x => x.Longitude);
  var minLonPoly = polygon.Min(x => x.Longitude);
  if (Math.Abs(minLonPoly - maxLonPoly) > 180)
  {
    offset = 360.0;
  }

  return offset;
}

private bool isPointLatitudeBetweenPolyLine(double polyLinePoint1Lat, double polyLinePoint2Lat, double poiLat)
{
  return polyLinePoint2Lat <= poiLat && poiLat < polyLinePoint1Lat || polyLinePoint1Lat <= poiLat && poiLat < polyLinePoint2Lat;
}

private bool isPointRightFromPolyLine(double polyLinePoint1Lat, double polyLinePoint1Lon, double polyLinePoint2Lat, double polyLinePoint2Lon, double poiLat, double poiLon)
{
  // lon <(lon1-lon2)*(latp-lat2)/(lat1-lat2)+lon2
  return poiLon < (polyLinePoint1Lon - polyLinePoint2Lon) * (poiLat - polyLinePoint2Lat) / (polyLinePoint1Lat - polyLinePoint2Lat) + polyLinePoint2Lon;
}
Nelly
  • 522
  • 6
  • 15
1

I add one detail to help people who live in the... south of earth!! If you're in Brazil (that's my case), our GPS coord are all negatives. And all these algo give wrong results.

The easiest way if to use the absolute values of the Lat and Long of all point. And in that case Jan Kobersky's algo is perfect.

Peter
  • 1,247
  • 19
  • 33
0

the polygon is defined as a sequential list of point pairs A, B, C .... A. no side A-B, B-C ... crosses any other side

Determine box Xmin, Xmax, Ymin, Ymax

case 1 the test point P lies outside the box

case 2 the test point P lies inside the box:

Determine the 'diameter' D of the box {[Xmin,Ymin] - [Xmax, Ymax]} ( and add a little extra to avoid possible confusion with D being on a side)

Determine the gradients M of all sides

Find a gradient Mt most different from all gradients M

The test line runs from P at gradient Mt a distance D.

Set the count of intersections to zero

For each of the sides A-B, B-C test for the intersection of P-D with a side from its start up to but NOT INCLUDING its end. Increment the count of intersections if required. Note that a zero distance from P to the intersection indicates that P is ON a side

An odd count indicates P is inside the polygon

0

I translated c# method in Php and I added many comments to understand code.

Description of PolygonHelps:
Check if a point is inside or outside of a polygon. This procedure uses gps coordinates and it works when polygon has a little geographic area.


INPUT:
$poly: array of Point: polygon vertices list; [{Point}, {Point}, ...];
$point: point to check; Point: {"lat" => "x.xxx", "lng" => "y.yyy"}


When $c is false, the number of intersections with polygon is even, so the point is outside of polygon;
When $c is true, the number of intersections with polygon is odd, so the point is inside of polygon;
$n is the number of vertices in polygon;
For each vertex in polygon, method calculates line through current vertex and previous vertex and check if the two lines have an intersection point.
$c changes when intersection point exists.
So, method can return true if point is inside the polygon, else return false.

class PolygonHelps {

    public static function isPointInPolygon(&$poly, $point){

        $c = false; 
        $n = $j = count($poly);


        for ($i = 0, $j = $n - 1; $i < $n; $j = $i++){

            if ( ( ( ( $poly[$i]->lat <= $point->lat ) && ( $point->lat < $poly[$j]->lat ) ) 
                || ( ( $poly[$j]->lat <= $point->lat ) && ( $point->lat < $poly[$i]->lat ) ) ) 

            && ( $point->lng <   ( $poly[$j]->lng - $poly[$i]->lng ) 
                               * ( $point->lat    - $poly[$i]->lat ) 
                               / ( $poly[$j]->lat - $poly[$i]->lat ) 
                               +   $poly[$i]->lng ) ){

                $c = !$c;
            }
        }

        return $c;
    }
}
0

Jan's answer is great.

Here is the same code using the GeoCoordinate class instead.

using System.Device.Location;

...

public static bool IsPointInPolygon(List<GeoCoordinate> poly, GeoCoordinate point)
{
    int i, j;
    bool c = false;
    for (i = 0, j = poly.Count - 1; i < poly.Count; j = i++)
    {
        if ((((poly[i].Latitude <= point.Latitude) && (point.Latitude < poly[j].Latitude))
                || ((poly[j].Latitude <= point.Latitude) && (point.Latitude < poly[i].Latitude)))
                && (point.Longitude < (poly[j].Longitude - poly[i].Longitude) * (point.Latitude - poly[i].Latitude)
                    / (poly[j].Latitude - poly[i].Latitude) + poly[i].Longitude))
            c = !c;
    }

    return c;
}
Fidel
  • 7,027
  • 11
  • 57
  • 81
0

Check if a point is inside a polygon or not -

Consider the polygon which has vertices a1,a2,a3,a4,a5. The following set of steps should help in ascertaining whether point P lies inside the polygon or outside.

Compute the vector area of the triangle formed by edge a1->a2 and the vectors connecting a2 to P and P to a1. Similarly, compute the vector area of the each of the possible triangles having one side as the side of the polygon and the other two connecting P to that side.

For a point to be inside a polygon, each of the triangles need to have positive area. Even if one of the triangles have a negative area then the point P stands out of the polygon.

In order to compute the area of a triangle given vectors representing its 3 edges, refer to http://www.jtaylor1142001.net/calcjat/Solutions/VCrossProduct/VCPATriangle.htm

Arnkrishn
  • 29,828
  • 40
  • 114
  • 128
0

The problem is easier if your polygon is convex. If so, you can do a simple test for each line to see if the point is on the inside or outside of that line (extending to infinity in both directions). Otherwise, for concave polygons, draw an imaginary ray from your point out to infinity (in any direction). Count how many times it crosses a boundary line. Odd means the point is inside, even means the point is outside.

This last algorithm is trickier than it looks. You will have to be very careful about what happens when your imaginary ray exactly hits one of the polygon's vertices.

If your imaginary ray goes in the -x direction, you can choose only to count lines that include at least one point whose y coordinate is strictly less than the y coordinate of your point. This is how you get most of the weird edge cases to work correctly.

Dietrich Epp
  • 205,541
  • 37
  • 345
  • 415
0

If you have a simple polygon (none of the lines cross) and you don't have holes you can also triangulate the polygon, which you are probably going to do anyway in a GIS app to draw a TIN, then test for points in each triangle. If you have a small number of edges to the polygon but a large number of points this is fast.

For an interesting point in triangle see link text

Otherwise definately use the winding rule rather than edge crossing, edge crossing has real problems with points on edges, which if your data is generated form a GPS with limited precision is very likely.

Martin Beckett
  • 94,801
  • 28
  • 188
  • 263
-1

You can try this simple class https://github.com/xopbatgh/sb-polygon-pointer

It is easy to deal with it

  1. You just insert polygon coordinates into array
  2. Ask library is desired point with lat/lng inside the polygon
$polygonBox = [
    [55.761515, 37.600375],
    [55.759428, 37.651156],
    [55.737112, 37.649566],
    [55.737649, 37.597301],
];

$sbPolygonEngine = new sbPolygonEngine($polygonBox);

$isCrosses = $sbPolygonEngine->isCrossesWith(55.746768, 37.625605);

// $isCrosses is boolean

(answer was returned from deleted by myself because initially it was formatted wrong)