Point in Polygon using Winding Number

Question

The question is: how do you determine if a point is within a polygon?

This question has been asked and answered many times. There are multiple methods for determining whether a point is within a polygon.

I've grokked the Winding Number algorithm, ported a solid answer from another SO thread into C# and written xUnit tests around it to ensure that I could refactor ruthlessly. The goal was to take an answer, all of which seem to use a procedural programming approach and variable names that are similar to those you'd find in a mathematical formula, and refactor it into a reasonably sound set of OOP classes and methods.

So, to rephrase this question specifically to the answer that I will go on to provide:

How do I determine if a location / point (latitude and longitude) is within a polygon in OOP C#?

Answer 1

The answer that I used as my starting point was provided by Manuel Castro in the following SO thread Geo Fencing - point inside/outside polygon :

public static bool PointInPolygon(LatLong p, List<LatLong> poly)
{
    int n = poly.Count();

    poly.Add(new LatLong { Lat = poly[0].Lat, Lon = poly[0].Lon });
    LatLong[] v = poly.ToArray();

    int wn = 0;    // the winding number counter

    // loop through all edges of the polygon
    for (int i = 0; i < n; i++)
    {   // edge from V[i] to V[i+1]
        if (v[i].Lat <= p.Lat)
        {         // start y <= P.y
            if (v[i + 1].Lat > p.Lat)      // an upward crossing
                if (isLeft(v[i], v[i + 1], p) > 0)  // P left of edge
                    ++wn;            // have a valid up intersect
        }
        else
        {                       // start y > P.y (no test needed)
            if (v[i + 1].Lat <= p.Lat)     // a downward crossing
                if (isLeft(v[i], v[i + 1], p) < 0)  // P right of edge
                    --wn;            // have a valid down intersect
        }
    }
    if (wn != 0)
        return true;
    else
        return false;

}

I proceeded to write xUnit tests around an implementation that started out using the exact logic expressed in the above code. The following are a bit overkill, but I preferred more tests to ensure that refactoring would not create issues. Using inline data in xUnit theories is so easy, eh, why not. With all the tests green, I was then able to refactor to my heart's content:

public class PolygonTests
{

    public class GivenLine : PolygonTests
    {
        private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
        {
            new GeographicalPoint(1, 1),
            new GeographicalPoint(10, 1)
        });
        public class AndPointIsAnywhere : GivenLine
        {
            [Theory]
            [InlineData(5, 1)]
            [InlineData(-1, -1)]
            [InlineData(11, 11)]
            public void WhenAskingContainsLocation_ThenItShouldReturnFalse(double latitude, double longitude)
            {
                GeographicalPoint point = new GeographicalPoint(latitude, longitude);
                bool actual = _polygon.Contains(point);

                actual.Should().BeFalse();
            }
        }
    }

    public class GivenTriangle : PolygonTests
    {
        private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
        {
            new GeographicalPoint(1, 1),
            new GeographicalPoint(10, 1),
            new GeographicalPoint(10, 10)
        });

        public class AndPointWithinTriangle : GivenTriangle
        {
            private readonly GeographicalPoint _point = new GeographicalPoint(6, 4);

            [Fact]
            public void WhenAskingContainsLocation_ThenItShouldReturnTrue()
            {
                bool actual = _polygon.Contains(_point);

                actual.Should().BeTrue();
            }
        }

        public class AndPointOutsideOfTriangle : GivenTriangle
        {
            private readonly GeographicalPoint _point = new GeographicalPoint(5, 5.0001d);

            [Fact]
            public void WhenAskingContainsLocation_ThenItShouldReturnFalse()
            {
                bool actual = _polygon.Contains(_point);

                actual.Should().BeFalse();
            }
        }
    }

    public class GivenComplexPolygon : PolygonTests
    {
        private readonly Polygon _polygon = new Polygon(new List<GeographicalPoint>
        {
            new GeographicalPoint(1, 1),
            new GeographicalPoint(5, 1),
            new GeographicalPoint(8, 4),
            new GeographicalPoint(3, 4),
            new GeographicalPoint(8, 9),
            new GeographicalPoint(1, 9)
        });

        [Theory]
        [InlineData(5, 0, false)]
        [InlineData(5, 0.999d, false)]
        [InlineData(5, 1, true)]
        [InlineData(5, 2, true)]
        [InlineData(5, 3, true)]
        [InlineData(5, 4, false)]
        [InlineData(5, 5, false)]
        [InlineData(5, 5.999d, false)]
        [InlineData(5, 6, true)]
        [InlineData(5, 7, true)]
        [InlineData(5, 8, true)]
        [InlineData(5, 9, false)]
        [InlineData(5, 10, false)]
        [InlineData(0, 5, false)]
        [InlineData(0.999d, 5, false)]
        [InlineData(1, 5, true)]
        [InlineData(2, 5, true)]
        [InlineData(3, 5, true)]
        [InlineData(4.001d, 5, false)]
        //[InlineData(5, 5, false)] -- duplicate
        [InlineData(6, 5, false)]
        [InlineData(7, 5, false)]
        [InlineData(8, 5, false)]
        [InlineData(9, 5, false)]
        [InlineData(10, 5, false)]
        public void WhenAskingContainsLocation_ThenItShouldReturnCorrectAnswer(double latitude, double longitude, bool expected)
        {
            GeographicalPoint point = new GeographicalPoint(latitude, longitude);
            bool actual = _polygon.Contains(point);

            actual.Should().Be(expected);
        }
    }
}

This allowed me to refactor the original code to the following:

public interface IPolygon
{
    bool Contains(GeographicalPoint location);
}

public class Polygon : IPolygon
{
    private readonly List<GeographicalPoint> _points;

    public Polygon(List<GeographicalPoint> points)
    {
        _points = points;
    }

    public bool Contains(GeographicalPoint location)
    {
        GeographicalPoint[] polygonPointsWithClosure = PolygonPointsWithClosure();

        int windingNumber = 0;

        for (int pointIndex = 0; pointIndex < polygonPointsWithClosure.Length - 1; pointIndex++)
        {
            Edge edge = new Edge(polygonPointsWithClosure[pointIndex], polygonPointsWithClosure[pointIndex + 1]);
            windingNumber += AscendingIntersection(location, edge);
            windingNumber -= DescendingIntersection(location, edge);
        }

        return windingNumber != 0;
    }

    private GeographicalPoint[] PolygonPointsWithClosure()
    {
        // _points should remain immutable, thus creation of a closed point set (starting point repeated)
        return new List<GeographicalPoint>(_points)
        {
            new GeographicalPoint(_points[0].Latitude, _points[0].Longitude)
        }.ToArray();
    }

    private static int AscendingIntersection(GeographicalPoint location, Edge edge)
    {
        if (!edge.AscendingRelativeTo(location)) { return 0; }

        if (!edge.LocationInRange(location, Orientation.Ascending)) {  return 0; }

        return Wind(location, edge, Position.Left);
    }

    private static int DescendingIntersection(GeographicalPoint location, Edge edge)
    {
        if (edge.AscendingRelativeTo(location)) { return 0; }

        if (!edge.LocationInRange(location, Orientation.Descending)) { return 0; }

        return Wind(location, edge, Position.Right);
    }

    private static int Wind(GeographicalPoint location, Edge edge, Position position)
    {
        if (edge.RelativePositionOf(location) != position) { return 0; }

        return 1;
    }

    private class Edge
    {
        private readonly GeographicalPoint _startPoint;
        private readonly GeographicalPoint _endPoint;

        public Edge(GeographicalPoint startPoint, GeographicalPoint endPoint)
        {
            _startPoint = startPoint;
            _endPoint = endPoint;
        }

        public Position RelativePositionOf(GeographicalPoint location)
        {
            double positionCalculation =
                (_endPoint.Longitude - _startPoint.Longitude) * (location.Latitude - _startPoint.Latitude) -
                (location.Longitude - _startPoint.Longitude) * (_endPoint.Latitude - _startPoint.Latitude);

            if (positionCalculation > 0) { return Position.Left; }

            if (positionCalculation < 0) { return Position.Right; }

            return Position.Center;
        }

        public bool AscendingRelativeTo(GeographicalPoint location)
        {
            return _startPoint.Latitude <= location.Latitude;
        }

        public bool LocationInRange(GeographicalPoint location, Orientation orientation)
        {
            if (orientation == Orientation.Ascending) return _endPoint.Latitude > location.Latitude;

            return _endPoint.Latitude <= location.Latitude;
        }
    }

    private enum Position
    {
        Left,
        Right,
        Center
    }

    private enum Orientation
    {
        Ascending,
        Descending
    }
}

public class GeographicalPoint
{
    public double Longitude { get; set; }

    public double Latitude { get; set; }

    public GeographicalPoint(double latitude, double longitude)
    {
        Latitude = latitude;
        Longitude = longitude;
    }
}

The original code, of course, is far less verbose. However, it:

1. is procedural;
2. uses variable names that do not reveal intent;
3. is mutable;
4. has a cyclomatic complexity of 12.

The refactored code:

1. passes the tests;
2. reveals intention;
3. contains no duplication;
4. contains the fewest elements (given 1, 2 and 3, above)

And:

5. is object oriented;
6. does not use if except for guard clauses;
7. is immutable;
8. hides its private data;
9. has complete test coverage;
10. has one method with a cyclomatic complexity of 3, while most of the methods are at 1, and a few of them are at 2.

Now, all that said, I'm not saying that there are no additional refactors that might be suggested, or that the above refactor approaches perfection. However, I think that in terms of implementing the Winding Number algorithm for determining whether a point is in a polygon and really understanding the algorithm that this is helpful.

I hope that this helps some who, like me, had some difficulty wrapping their heads around it.

Cheers