A View Inside My Head

Jason's Random Thoughts of Interest

NAVIGATION - SEARCH

SqlGeography: Ring Orientation of Polygon Interior Rings (Holes)

I have mentioned before how the Ring Orientation for the exterior ring of a Polygon is significant when instantiating a SqlGeography object.  In this case, a Counter-Clockwise orientation is required so that as an observer walks along the path, the interior of the Polygon is always to their left.

Ring Orientation for SqlGeographyBut, what I have never really seen documented (or paid attention to, at least) is the fact that the interior rings, or holes, of a Polygon also have specific Ring Orientation requirements. 

In keeping with the "Left-handed" rule, interior rings must be defined in a Clockwise manner - the opposite orientation of the shape's exterior ring.  This is because holes within a Polygon are considered to be part of the exterior of the shape, so the observer walking in a Clockwise direction is still keeping the Polygon's interior to their left.

(I should note here that the Ring Orientation for SqlGeography is the exact opposite of ESRI's ShapeFile format, which is why Ring Orientation has been on my mind for the past few days).

 

Spatial: Determining Ring Orientation

A Ring is a list of points such that the starting point and ending point are the same (forming a closed shape).  The order the you define the points that make up a Ring  - known as Ring Orientation - is significant, for various data formats (including SQL Server's Geography type) imply special meaning for rings that are defined in a clockwise manner as opposed to a counter-clockwise manner. 

Given a list of points with no additional context, it can be difficult to determine the Ring Orientation being used. 

For example, suppose that you have a generic list of points that represent the boundary of a postal code, and that you wish to use these points in order to construct a Polygon instance using the SqlGeography type.  SqlGeography happens to use "Left-handed" ordering, so that as an observer walks along the set of points in the order defined, the "inside" of the polygon is always to their left.  This also implies that the exterior ring of a Polygon is defined in a counter-clockwise manner.

If you try to define a polygon with an area greater than a single hemisphere (this is a nice way to say "if you screw up and use the wrong orientation"), then the SqlGeography type will throw an exception.  So, aside from using Try-Catch, what can you do?

While researching solutions to this problem, I stumbled upon a paper entitled "A Winding Number and Point-in-Polygon Algorithm" from the Colorado State University.  It turns out that a simple algorithm with O(n) complexity can be used to determine if a point is within a Polygon, and a side effect also provides the Ring Orientation.  The key to this algorithm is determining the trend of the ring at each crossing of an axis.

Since I was only interested in Ring Orientation (and not point enclosure detection), I didn't need to use this particular algorithm.  Instead, I took inspiration from the winding concept, and created a simpler derivative algorithm:

Visual example of how to determine ring orientation at the extreme left point

  1. Iterate the point collection and determine the extreme "left" and "right" points
  2. Normalize the line segments connected to these points so that they each have the same "X" dimension length
  3. Compare the "Y" values of the normalized segments to establish the trend through that extreme point (i.e., is the "previous" segment above or below the "next" segment)
  4. In the spirit of the Winding algorithm, use opposite orientations for the left and right points so that the results coincide with one another
  5. A negative result (negative indicates Clockwise orientation, positive result indicates Counter-Clockwise orientation, and a result of zero would be undefined

I've actually written (and posted) several versions of this algorithm, each time discovering some edge case exception that would cause me to take down the post and rewrite the algorithm.  I believe the code below works for all simple polygons on a Cartesian coordinate system (read: I have more testing to see if this will work with an ellipsoidal model, like SqlGeography). 

Note: The following code is generic in nature, and as such, I've defined my own Point structure instead of using a SqlGeometry or SqlGeography, etc.

struct Point
{
    public double X { get; set; }
    public double Y { get; set; }
}

enum RingOrientation : int
{
    Unknown = 0,
    Clockwise = -1,
    CounterClockwise = 1
};

RingOrientation Orientation(Point[] points)
{
    // Inspired by http://www.engr.colostate.edu/~dga/dga/papers/point_in_polygon.pdf

    // This algorithm is to simply determine the Ring Orientation, so to do so, find the
    // extreme left and right points, and then check orientation

    if (points.Length < 4)
    {
        throw new ArgumentException("A polygon requires at least 4 points.");
    }

    if (points[0].X != points[points.Length - 1].X || points[0].Y != points[points.Length - 1].Y)
    {
        throw new ArgumentException("The array of points is not a polygon.  The first and last point must be identical.");
    }

    int rightmostIndex = 0;
    int leftmostIndex = 0;

    for (int i = 1; i < points.Length; i++)
    {
        if (points[i].X < points[leftmostIndex].X)
        {
            leftmostIndex = i;
        }
        if (points[i].X > points[rightmostIndex].X)
        {
            rightmostIndex = i;
        }
    }


    Point p0; // Point before the extreme
    Point p1; // The extreme point
    Point p2; // Point after the extreme

    double m; // Holds line slope

    double lenP2x;  // Length of the P1-P2 line segment's delta X
    double newP0y;  // The Y value of the P1-P0 line segment adjusted for X=lenP2x

    RingOrientation left_orientation;
    RingOrientation right_orientation;

    // Determine the orientation at the Left Point
    if (leftmostIndex == 0)
        p0 = points[points.Length - 2];
    else
        p0 = points[leftmostIndex - 1];

    p1 = points[leftmostIndex];

    if (leftmostIndex == points.Length - 1)
        p2 = points[1];
    else
        p2 = points[leftmostIndex + 1];

    m = (p1.Y - p0.Y) / (p1.X - p0.X);

    if (double.IsInfinity(m))
    {
        // This is a vertical line segment, so just calculate the dY to
        // determine orientation

        left_orientation = (RingOrientation)Math.Sign(p0.Y - p1.Y);
    }
    else
    {
        lenP2x = p2.X - p1.X;
        newP0y = p1.Y + (m * lenP2x);

        left_orientation = (RingOrientation)Math.Sign(newP0y - p2.Y);
    }



    // Determine the orientation at the Right Point
    if (rightmostIndex == 0)
        p0 = points[points.Length - 2];
    else
        p0 = points[rightmostIndex - 1];

    p1 = points[rightmostIndex];

    if (rightmostIndex == points.Length - 1)
        p2 = points[1];
    else
        p2 = points[rightmostIndex + 1];

    m = (p1.Y - p0.Y) / (p1.X - p0.X);

    if (double.IsInfinity(m))
    {
        // This is a vertical line segment, so just calculate the dY to
        // determine orientation

        right_orientation = (RingOrientation)Math.Sign(p1.Y - p0.Y);
    }
    else
    {
        lenP2x = p2.X - p1.X;
        newP0y = p1.Y + (m * lenP2x);

        right_orientation = (RingOrientation)Math.Sign(p2.Y - newP0y);
    }


    if (left_orientation == RingOrientation.Unknown)
    {
        return right_orientation;
    }
    else
    {
        return left_orientation;
    }
}


void Test()
{
    // Simple triangle - left extreme point is vertically "in between" line segments
    Point[] points = new Point[]
    {
        new Point(5,-1),
        new Point(0,0),
        new Point(5,1),
        new Point(5,-1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where both line segments are above the left extreme point
    points = new Point[]
    {
        new Point(2,1),
        new Point(0,0),
        new Point(1,1),
        new Point(2,1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where both line segments are below the left extreme point
    points = new Point[]
    {
        new Point(2,-1),
        new Point(0,0),
        new Point(1,-1),
        new Point(2,-1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);

    // Case where line segment is vertical (slope cannot be determined)
    points = new Point[]
    {
        new Point(0,0),
        new Point(0,1),
        new Point(1,1),
        new Point(1,0),
        new Point(0,0)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where angle thru left extreme point is a right angle
    points = new Point[]
    {
        new Point(0,0),
        new Point(1,1),
        new Point(1,-1),
        new Point(0,0)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Real-world case from a SHP file
    points = new Point[]
    {
        new Point(-156.92467299999998,20.738695999999997),
        new Point(-156.924636,20.738822),
        new Point(-156.924608,20.73894),
        new Point(-156.92458,20.739082),
        new Point(-156.92460599999998,20.739234),
        new Point(-156.924551,20.739326),
        new Point(-156.924507,20.739241999999997),
        new Point(-156.924482,20.739082),
        new Point(-156.924466,20.738854999999997),
        new Point(-156.924387,20.738602999999998),
        new Point(-156.924308,20.738325),
        new Point(-156.924239,20.738063999999998),
        new Point(-156.92424,20.737887999999998),
        new Point(-156.924285,20.737811999999998),
        new Point(-156.924475,20.73762),
        new Point(-156.92458299999998,20.737603999999997),
        new Point(-156.924754,20.737579),
        new Point(-156.924851,20.737731),
        new Point(-156.924956,20.738101),
        new Point(-156.924909,20.738343999999998),
        new Point(-156.924818,20.738487),
        new Point(-156.92467299999998,20.738695999999997)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

}