Point in Polygon algorithm giving wrong results sometimes

Question

I saw on StackOverflow a "point in polygon" raytracing algorithm that I implemented in my PHP Code. Most of the time, it works well, but in some complicated cases, with complex polygons and vicious points, it fails and it says that point in not in polygon when it is.

For example:
You will find here my Polygon and Point classes: pointInPolygon method is in Polygon class. At the end of the file, there are two points that are supposed to lie inside the given polygon (True on Google Earth). The second one works well, but the first one is buggy :( .

You can easily check the polygon on Google Earth using this KML file.

Answer 1

Have been there :-) I also travelled through Stackoverflow's PiP-suggestions, including your reference and this thread. Unfortunately, none of the suggestions (at least those I tried) were flawless and sufficient for a real-life scenario: like users plotting complex polygons on a Google map in freehand, "vicious" right vs left issues, negative numbers and so on.

The PiP-algorithm must work in all cases, even if the polygon consists of hundreds of thousands of points (like a county-border, nature park and so on) - no matter how "crazy" the polygon is.

So I ended up building a new algorithm, based on some source from an astronomy-app:

//Point class, storage of lat/long-pairs
class Point {
    public $lat;
    public $long;
    function Point($lat, $long) {
        $this->lat = $lat;
        $this->long = $long;
    }
}

//the Point in Polygon function
function pointInPolygon($p, $polygon) {
    //if you operates with (hundred)thousands of points
    set_time_limit(60);
    $c = 0;
    $p1 = $polygon[0];
    $n = count($polygon);

    for ($i=1; $i<=$n; $i++) {
        $p2 = $polygon[$i % $n];
        if ($p->long > min($p1->long, $p2->long)
            && $p->long <= max($p1->long, $p2->long)
            && $p->lat <= max($p1->lat, $p2->lat)
            && $p1->long != $p2->long) {
                $xinters = ($p->long - $p1->long) * ($p2->lat - $p1->lat) / ($p2->long - $p1->long) + $p1->lat;
                if ($p1->lat == $p2->lat || $p->lat <= $xinters) {
                    $c++;
                }
        }
        $p1 = $p2;
    }
    // if the number of edges we passed through is even, then it's not in the poly.
    return $c%2!=0;
}

Illustrative test :

$polygon = array(
    new Point(1,1), 
    new Point(1,4),
    new Point(4,4),
    new Point(4,1)
);

function test($lat, $long) {
    global $polygon;
    $ll=$lat.','.$long;
    echo (pointInPolygon(new Point($lat,$long), $polygon)) ? $ll .' is inside polygon<br>' : $ll.' is outside<br>';
}

test(2, 2);
test(1, 1);
test(1.5333, 2.3434);
test(400, -100);
test(1.01, 1.01);

Outputs :

2,2 is inside polygon 
1,1 is outside
1.5333,2.3434 is inside polygon 
400,-100 is outside
1.01,1.01 is inside polygon

It is now more than a year since I switched to the above algorithm on several sites. Unlike the "SO-algorithms" there have not been any complaints so far. See it in action here (national mycological database, sorry for the Danish). You can plot a polygon, or select a "kommune" (a county) - ultimately compare a polygon with thousands of points to thousands of records).

Update Note, this algorithm is targeting geodata / lat,lngs which can be very precise (n'th decimal), therefore considering "in polygon" as inside polygon - not on border of polygon. 1,1 is considered outside, since it is on the border. 1.0000000001,1.01 is not.